{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"import glob\nfrom skimage import measure\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"ACKNOWLEDGEMENT\n* https://www.kaggle.com/gzuidhof/full-preprocessing-tutorial\n* https://www.raddq.com/dicom-processing-segmentation-visualization-in-python/\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"%matplotlib inline\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport pydicom\nimport os\nimport scipy.ndimage\nimport matplotlib.pyplot as plt\n\nfrom skimage import measure, morphology\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In this notebook, we learn\n\n* Loading DICOM data using pydicom\n\n* Used 1D (histogram), 2D, and 3D plots to display DICOM images.\n* Pre-processed data for future machine learning projects\n* Conversion of pixel value to Hundsfeld units\n* Resampling for isotropy\n* normalization\n* zero centering","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Loading the files\n\n\nDicom is the de-facto file standard in medical imaging.  These files contain a lot of metadata (such as the pixel size, so how long one pixel is in every dimension in the real world).\n\nThis pixel size/coarseness of the scan differs from scan to scan (e.g. the distance between slices may differ), which can hurt performance of CNN approaches. We can deal with this by isomorphic resampling, which we will do later.\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\ntest = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\nsubmission = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def myfunc(e):\n    return e[-6:] ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_list = train.Patient.unique()\ndata_path = '../input/osic-pulmonary-fibrosis-progression/train/'\n\npatient = pd.DataFrame()\npid = []\ncount = []\npath = []\nfor pat in patient_list:\n    \n    data_list = glob.glob(data_path + pat + '/*.dcm')\n    data_list.sort(key=myfunc)\n    pid.append(pat)\n    path.append(data_list)\n    count.append(len(data_list))\n    \npatient['pid']=pid\npatient['path']= path\npatient['count']=count","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"patient.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true},"cell_type":"code","source":"patient.path[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"slices = [pydicom.read_file(s) for s in patient.path[0]] #lets read metadeta of pydiacom file\nprint('The total no of ct scan associated with 1st patient',len(slices))\nprint(slices[3])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"slices.sort(key = lambda x: int(x.InstanceNumber)) # order the slice serially\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Slice Thickness is a parameter that can be selected by the technologist.  This will change the thickness of our slice in millimeters.  By increasing the slice thickness, many more different types of tissues will be collected in our 2D slice.  This can cause blurring in our image also known as partial voluming.\nPicture\nCoverage\n\nBy increasing the slice thickness we will increase our coverage. \n![1.jpg](attachment:1.jpg)","attachments":{"1.jpg":{"image/jpeg":"/9j/2wBDAAIBAQIBAQICAgICAgICAwUDAwMDAwYEBAMFBwYHBwcGBwcICQsJCAgKCAcHCg0KCgsMDAwMBwkODw0MDgsMDAz/2wBDAQICAgMDAwYDAwYMCAcIDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAz/wAARCACcAYYDASIAAhEBAxEB/8QAHgABAAEEAwEBAAAAAAAAAAAAAAkBBgcIAgUKBAP/xABhEAABAgQCBAMPDA4HBQkAAAABAgMABAUGBxEICRIhEzHSFBUXGBkiQVFUVVeTlJXRChoyM1NWYXGRs7TUIyQ0Nzg5YnJzdaOxxNUoQlh2d7LBFidDR4E1NkRFZnSCg6T/xAAdAQEAAAcBAQAAAAAAAAAAAAAAAQIDBAUGBwgJ/8QASREAAQIDAwQMCQoGAwEAAAAAAQACAwQRBSExBhJBUQcTFFJTVGFxkaGx0RUWFyIycoHB8AgYIzM0NUKSotI2YnOy4eIkVcIl/9oADAMBAAIRAxEAPwCfyEIQRIQhBEhCEESEIQRIQhBEhCEESEIQRIQhBEhCEESEIQRIQhBEhFt4sKuo4aV0WQugou8yTooyq2l1VNTNbP2IzAa+yFray2tjrss8oietrXcaTktguxiBVrawcqtGVikxhmKZR5CoJqc1Mheb6muGmQ0AtGSWlKWMlnr05RRfHYxwa7T3gdpCzlk5PzVosc+XLfNIFCaGprSnQehTCwiOxGuepuJuNmCNQpFQr1k2XcMxdkrdVFrFvyszMpco0upx3bmW5g8CEBO0Cyl3bz2Tsx02kBr9n0aH1dxEw9w6uu35qnuUmdps1e9FUaTW6ZOThl1PsOyr5SVpyJ4NS0uJzB2CAYk3XDzS6t3+K16FdMyQtR0SHC2uheQL7qEuLaGumrThVSWQjX/Rn09ZDS70fLrxIsa07nnbepCp1qguTTbTLl3KlUr2lSrSVKdQhbqODTwqULJ/q7o0Y0d9eZjFilg9itc1XYwJRUrMs+p19q2pZ+pSdw0CclVdaxOSMyUqmWikp2nWFpCCRmQolImiTLIfpaq+wKjJ5MT8wYgY0Asc1pBIBDnGgFO3VhjcpZ4RGjYmupmLAuS4brxOqj6rapWEVs3m7b9Gt5pOVSqfM6eDlplUwXXCtx7ZDbiUpQDmXDsknJV+a9mwsMLAo89cWH+J9v3VWq1N0Vm1qxJydLnAZZpLzszzRMTDcmqW2Fo2XEvEKUsJGZiO6Id9Th307VGJkrabXhjYRdXCnNnHGmAvJw5VvJCNQtIzT7ua69WWjSGwAXaM7ItU1dxvtXfKzJSqnsod5oYS3LrBE0lxAQNpex1q8zxGNUrZ16WNeEtSwsq+KFl2XctqX/YlQxAmZexpR9upyNPZQpSCTOzSWgpGwtTgBVmnLYJO6IPmWNcWu0dGnuKnkslZ+bhGLBAqHFuaTR2c0VIpropaoRHW7rqqJYOMeK1zVyrV6p4bUKzrZr1BoEtbsu1PrfqyGy00iYMxtOurU5kpLiUIbyJClAZn4NIfX2tYJXphNUatZ1wWFZdcqtdpl6SV2UB9Fcpy5BhC2uY+CcLLqXVuIAcTwjZ2vZJ2VEQM1DAqTpp15vaowskbUivbDhwiS4VHL5mfQVvJpqurdVSRwjQLFnX+WhgfVESlx4TYqSrsva1PvGqBk0x40anTrzbTK3gJrevada2kI2iC4BvyOV01LXq4N0zSuThatu4DlXGLZXX/ALVEgipPNKcQxwJe5rUnNJQXUsFtK8klQzGc26IdaV+MO0U51b+LNq5m2CA4ihNReKCmrnHLet1IRotaOvsw0uzAu+cTE2XiVL2JZsqZlFY5kk32akozYlEMZNTC1Sr63Ckpbmg0SglXEI2S0Q9KVnS3wxVc0vat0WmymaXLIZrAlnEzaQAQ/LzEq68w+yoHILbcOSkqSciInZEa70TXA+w4K1m7GnZVhfMQy0A5prrpWnOKiurSsrwhCJ1jEhCEESEIQRIQhBEhCEESEIQRIQiilhPHBFWEWBfelZhfhbczlFubEaw7crDKELXIVS4JSTmUJWM0EtuOJUAobwSN/YjtaJjjZt0WZVbipV1W9VqDQi8KhUJCoNTUtJKZRtupcW2pQSpCcipJOYBGY3xLnDGqrGXihocWmh5FdUItvCTF+2ceMOaVd1nVqQuO2a20XpCpSLnCy80gKUgqQrsgKSofGDH1XxiLQMMrbfrNyVqk29SJUpD09U5xuUlmto5J2nHFJSMychmd5iaqkzHZ2ZS/Cmld1CLNZ0irAmadQpxq97SelboLqaK83WJdaKwWwS4JYhZ4cpAJIb2iMt8frgpjvZ+kbh3J3bYtxUu6raqC3W5apU54PS7ym3FNuBKhxlK0qSfhBiAINwUzoMRrc5zSBzfGo9Cu2EceGBGe/IRYmIGlRhjhNcRpF1YiWLbNWS0l8yVWr8pJTAbVnsr4N1xKtk5HI5ZHIwJClhwnvNGAk8iv2EWva+NlnXzbtTq9Dum3q3S6KtxuoTdOqLM2zJLbQHFocU2pQQpKCFFJ3gEHLfHHBvG209IPDinXfZNfptzWzVg4ZOpU93hZeY2HFNr2VdnZWhST2ikwB0KLoL2glzSALjdp1K6oR0ltYlW9edaq9No9do9VqNvvplqpKyc60+9TXVDMNvoQoqaWRvCVgHKO54UfDAEKRwINCuUI48KPhz7UC6Act/yRGqguUIoFgjOKjeIIkIQgiQhCCLitsKGXFGi9G1HdBo+CdOstOIFdWxT8Wk4riaNNYC1vpKTzGU55BvrfZ+y38Ub1RTa35RTdCa45xF/+Qe0BX8naczKAtl35oNK4aKgY856VoZZeodtG3bltGaqF712ryFs1S6p96SVIMspqCK8yWX2StJzRwaVEpUnMk8YjqRqFhUdGKpYU1nHbEeu2tsU2TokjNy0vzBQpWTmTMBKJVJDTjzhOyt9WSikZZZHKJCSrKG0M4piVhAZubdh1U7Fk/G21qh23GoIIubcQ4uqLrr3E3Y1oblrXor6tih6LOGWKlg0y5a1MYf4jVOdqElRmUiRNrtziFIeYlH2jtpRkUlB3FBQCN5MYFoPqfukurrkxdeMl/XvUnrNnLIoM5VJSV4eiSUylSFKWtAC5txLathJeVuGfYCQmQ/MdsQzHbETOl4ZpUYCnsw+NSoQspLShRHxWRSHPIJNBUkacOmmOmq0FujUHWbfFuXDS6tfFxvsVvD2g2M0tqTYackl0hTC5eeScyFLUthJU2Rs5LUM+Ij7r91KT2J1t2tP3Djdfdx4j2hXJyr065a5T5OqS0s1NMoZXJN015KpZuXSlttSUgdasFXZAG9mY7YhmO2Ibnh3mmPfWvT8UVUZVWqKfTG7C4aW5pGGBbcRgdNVgZ/QQpfSAVDANFy116n1K3JigPV2cDcxPuKfCuEmFDJKCoqWohIASNw4hGCrz1Glv3jbFjUxeIFdYRZGF0/hi0tNNYJmmZtpxszahtdatPCEhA605DfG92Y7YhmO2Ii+Cx1c4Y/57yraWt6flzWDEI84u0ekbicFH/XtQNaF1WvdlMn79uNw3FatuW9LTDcjLodpj9FS3zPOJ3kLKy31zahs5LUARuI7HEfUcUzH1uzF4o4sXxiZNW1N1ucnF19ll9E8qoy6WUtst57Mo1LlCHG228wFpz3ZgDe/MdsQzHbESuloRqCMbz017fiium5WWsHB4jGowNBUVbm3GlR5t13Pjeo4a76nqp942TV6RW8Yboqb1TsKnWAmcXRpZLrUnIz7E2w4evO0sJYSzv/q7+MRkCiak23rX0nZy/KTiDc1Ko1Wrstc9RocrT5NLkzUWGSgFM8UGYbl1uHhVMJOyVgb8gnLd7MdsQzHbERECGDUC/wDyT2kqD8qrVcww3RvNNRSjaX00U/lFDopdRaDWtqNhb9/3neS8Y7i/20umhOUBNYk7UoskXWnJnh3HJ+XQxzPUHVgJbUp1tOaB2FBKk5x1d+ryomr4sa6aXSa1MVmYvCuu16eKJBimSEs6pIQGpWSYAal2glIGynjy7AAA2H2hDaETMhMaatHJ7+1W03b09MwjBjPq00qKNGGGAFw1YKsIoDnArAGcVFiFWEUSsK4t8VzgiQhCCJCEIIkIQgiQhCCJHFSczHKEEUZOtx0dLPvbWmaIEzP2Lb9XTddfnpW4nn6M3MCqsMssBlubUUEOIQCdkOZgZnKNLNFWrXvgji5ItWtW76olGvK+sSqfV6JLTMyijvsytKbXLL5ly4JLgWvc7ltHg0AHrAI9A60FRG8iOAbIG8nd8MWD5Kri4OpWvWG/tW8SOW0SXlGyb4We0NDbzdcXkGlDv/0heeWt41YxOYC2O5J4hYu0RVF0e5i6GkUmtz8gh6qt1x1lDrqGiA4rglkEKGZCU57gI3q1ylLlMWdXXgXUr1dvNkc8qLV56uyNvIr1NpU0ZMZvVankhT0qpTiwUoSTtHLZVtbCpMS2cuM/LDg8x7I7vhyicytWlpNa06iT7wPYFJHyxz5iBMMgBphOcbiKnOAurm6MbwRfgoG9GupVS7rh0T6iMM6LY0mxindDknO29RpumU25G+dzQ55tSz+apZK1J2dkBCM0dahGRSLfwt0i8SKTo84CN37iNjjY2Hs9bl0TKavaCZrnpP19urPplmFpQBw6w2Gghh0hCkqJ7JMegTgs+yc+zvO+BaIAyUr5YlbKFv4vigHuqrl2XDDUGWFL6X1IqXm4kG/z9RF2GqE2e0idIt7WdopdQv8AxCo9QYu23Jag25OSFRTz9ojssS+pdOl08wBDic1zL7qtph0ZJIAVs7Da5zR6tG+tPbQ+nZ+x6BWHLmvg0yvvv0duZNTkkBjYYmlFB4RlO0vJLnWjaV8MSWFGXb+WCm9/GRl2on3KCwNca0IPQQfcse/KwiZhzMGEGZjHMo26tWltbgMK1ANb9KgO0Z5m8MC9KuYNpVa9LdpFz42XrQ6lRJKYmGaLMyTNLSpkqlAOB2gpZCXdnPJCAD1gy+LDTSMxKpWCeADN+4j47WNYU5YFVnKbUrNbmXKlVbjTXHkIYU2kATTgaDQSw8QjYVn2czP7wf5R+UxXg93slbvhinDki0AZ2rq+K86ykfL0R3l8SWaSaX1Fbs6+pab/ADqXg3DogeszEi/MHtabe1YrdUvSj4M1/FenyN2uU1btJnpupO0x0yS5wy6QsSoc2i80hYb2sweIR1WjdpC6R9HpFDqdtX1i5ct33bYl9zJp9Uqk5VGFTkjMLRKKYl3ypCHUIQSjZGalADfnkZ9g1s5b1D/rAo7GZO/txKJEgAB5ur1gDq0KD8umOZR0owmjQa4HMDheKYmorffmhQO4H6QeM1e0Sb/cOkDd0lTTL2s8/Vqgi5JpinVJw/bVPdqpZVMSSnxufLALEupKQFJCwTzxW0v8RZzRbwvuNN94w0GUkU11Btieu2sStRu3YmUBucpdeZk8pwIKthiXnkrUQM1cMneZ3w3u9kr5YBrLsq+WKplifxfFfi7BUXZZwnRNsMsPSriBi3NIuYMca4jQRfW2cEa7N3Tg7atTn5OtU+eqNHlJmYlqwECoyzi2UKUiYCEpSHgSQvZSkbQOQA3C6huEcUIIPwRyi8rU1WjPIJJCQhCClSEIQRM8o+OuVyTtykzVQqE1KyMhIMrmJmZmHUtMy7aAVLWtaiAlKQCSSQABvjFWmVpz4b6B2FT124jV5mmSpChJSDRS7Uau4CkFqVY2gp1Q2k5kZJQDmpSRviJnErFrH/Xr150Tb1VwV0bipPBU5pRXN3GlKlAlask80naGeS8pdvJOQcWNo0IkcNOaL3avjBbXk5kjN2qdt+rgg3vOHM0YudyD2kLO+lZr/LmxHxmcwx0RrMksSq9TXiqp3HUWyukIbQsoXwKQ42FNFRSOaVuIQd4SleYVFmHWD6yTsYUYPZ//AFfzKM4aPujXZei9YLFu2VQ5WjyLeyp5xI25iecCQC6+4eucWcs8zuHYAG6L7iUQYhve815MF1uWyfsWXhiFDlWvA/E+pcTrNCAOYCgWqnVCNZN4KMHv2X8yh1QjWTeCjB79l/Mo2rhDc5356u5XHgqyOJQuh37lqp1QjWTeCjB79l/ModUI1k3gowe/ZfzKNq4Q3Od+eruTwVZHEoXQ79y1U6oRrJvBRg9+y/mUOqEaybwUYPfsv5lG1cIbnO/PV3J4KsjiULod+5aqdUI1k3gowe/ZfzKHVCNZN4KMHv2X8yjauENznfnq7k8FWRxKF0O/ctU+qEayUf8AKnB75Wv5lDqhOsl8FOD3ytfzKNrMs+OMEaf2nPR9BPCSWrs5TXK3V61MLkqTTkPpaDzobUsuOEnaDKMkhSkAnNaRuJESvhZjS5zzQc3cqkKxbLivEOHIwy44Cjv3LGl860XWFYYWnO1648PsD6JRacjbmZ2dfZZYZHYzUqp5ZniAG8ncATGuSfVU2kogbaqLhKpCRtlPOKbG0Bvy+692caXaVGmjiBphXeuqXlWnX5RtZMlSZclqn09OZyDbIORUAci4rNashmo5RiZ45su/o1f5TGGjTz6/RuNOWnct1ldjyxdqzpqUh15Abuk3r2Z4Z3U7e+HdCrT7SGHqtTpedW0glSWlOtIWUgneQCrKMf6aOlHI6J+Cc/cD3AzFXmM5SkSalZGamlA7OfGdhO9SjlkAMuMiLswCOzgVZnw0KQ+jNxE1rENIyt6Q+P089PSdTpVDt916mUiRnGFsqaSheTrqkq3cI4pIJy4kpbHYi8tm0dySxe30jcO9eBNkXKfwHIPfBH0jyWs5OX2dtFthqttYHWcYa/O2JftXVUq++VzlHn3W0NrnEAbTsuQgAbSACtO72GY37O/eZBz7MQA2lXanbF006pUV+YlqvT5lEzJPMJ23GnkKCkqSnI5kEcWWR357om70U8ZJrH7AS27snqY/R56qypMzKuoUjYdQtTa1JCgDsKUkqSct6SIscnLUdMQzBi3ubp1jvWr7FeV8a05Z0hOEuiw7w435za6TrBNOUEaisiQhCNmXXEhCEESEIQRIQhBEhCEESGQhCCJkO1DLfCEETKBGcIQRMoZQhBFTZHaiuW+EIIkMoQgiZQhCCJCEIIkUUchFYEZwRcFuAcZ+SI/dZVr1rb0YLlfwwwkp6cUca5iaFPRTZRC36fSHyoJUiYW2dpbwzy4FveDntqRlkbJ9U26UmImj1gDh/RrCuqoWuzf9WnKTWFyRS0/NMCXSUth7IuNJJWcy2Ukg5ZxEFo83Xi/oqT1QmbHq9q0moVIkTE87TpacnFJPG2HnmVLS2SMylJAJ3kExhbRteFLvEJ7s3l7l3LYy2JJi3pXwsW7awGgYKipB/GdA5BedYUg2COrVuXHjFF/FzSprfRAvyolLjNCU8ldMpiQCA24hADSgnJJDTOTScszwhJjdhpCW20pSlKEpASlKQAEgbgABxADiHYiHTqgmldmP94Fv7t3/AGPI/V4HWCaV5/5g0DzPI/V4toWUNnQx5ru9dtibGeUT6N2mjRcABQAagAKAKY2EQ5dUE0r/AAg0DzPI/V4dUE0r/CDQPM8j9Xip4zWfvuxUvJdlDwPb3KY2EQ5dUE0r/CDQPM8j9Xh1QTSv8INA8zyP1eHjNIb7sUfJdlDwPb3KY2EQ5dUE0r/CDQPM8j9Xh1QTSv8ACDQPM8j9Xh4zSG+7E8l2UPA9vcpjYRDl1QTSv8INA8zyP1eHVBNK/wAINA8zyP1eHjNIb7sTyXZQ8D29ymNgTl/0iHLqgmlf4QaB5nkfq8dRfWsr0pLStl+bncQqUmXOTKlS1Ikg6krzSCk8AMiOPMHMbjEzco5FxDWuvVKNsaW9CYYkSFRoxN93Ut+dYdrR7W0LqPMUSmcBcWIs1LFcpTUq25anEnJLs4pJBSOMhodevZ37AO1EMWOmO106R2JNQuu76q/VaxUVlSlqOy0wjiS00jibbSAAEp3ZDfmczFr1OpzFXqMxNzb703NTbqnn33llx15xRJUtSjvUokkkneSY+fOLaZm3xjfhqWYsmxIMk2uL9J7tQTOKO+0O/o1fuMVj7bdtioXpXZOj0mSm6lVKs8mTk5SVaU6/NPOHYQ22hOZUtSiAEgZkxaEVFFl4jg1pc40AXsS0fRngXZmfeGQ+jNxpLr2fb8MRmcsqkeP/ANtG6ejbL1qR0fLHZuKk84a+zQJBupUzmhMzzvmEy7aXGeET1q9hQKdobjlujSzXs/dGGPxVL+GjL2/9gf7O0L5hbK33BN05P7wtXNX2f6a+Gv64HzLsTWNb0A9mIU9X3+Gvhp+uB8y7E1jPtYjHZJfZX+t7gtS2FfuqP/UP9rVyhCEbWuypCEIIkIQgiQhCCJCEM8oIkIZwgiQhCCJCEIIkIQgiQhCCJCEIIkIQgiQhCCJCEIgUURXqsUZ4bYEf3qnPo7URyuH7Kr84/viRr1WJuw3wI/vVOfR2ojlcH2VX5x/fHNMsvtTebuX0K+S7/C7/AFz2lUhCEaevSyQhCCJCEIIkIQgiQhCCJFkaQ/3rpv8ATs/5oveLI0h/vXTX6dn/ADRdyP2hnOFh8oPu2P6p7Frsd8IcUI6CvP6/aQlhOzzLJdZYDy0o4R1Wy23mQNpR35AcZPajJejHpdXxoVYgzNyYd1CiU64Cgy7NVmKFJ1F6XR1wKpdU0ysslQJzUgJUU5A7sxGLoK4jEWuINQrealYUzDMGYaHMOIIqDzhev3QjxQq+N+h1hXedwOsPV26rTplWqC2Ggy0uYflW3HClA3JBUo5AbhGo+vZ+6MMfiqX8NGy+rI/Fy4E/3Aon0FmNaNe190YY/FUv4aMpbp/+e/2doXzN2XWNbYk61ooAf/YWrmr7/DXw0/XA+Zdiaxn2sRCnq+x/TWw0/XA+ZdiatogNiLDJL7M/1vcFpuwr91R/6h/taucIptDKG2O3G1LsirCKbQyhtCCKsIptCK5iCJCKFYHGYQqoVCrHxSlblKhPzUtLzUs+/JKSiYbbdSpbCiMwFgHNJI378t0Yv0ztKOn6KGCc/cL/AAMxVn85WkSS1ZGcmlA7IIBz2E+yWRxJHbIiLnRK0y61o/aSbl61SbmanK3G+pNzbWRXPNuL2lPZAH7I2rr05DPIFAyCoxE9bMGVjMgv048gWjZSZeyNjz8CRj3l58419BpuBOu/RqqdSmiT2YrHw23XpK6KHK1KnTLc5IT7KJiXebOaHW1pCkqB7RBBj7HFBIzJyjL1BvC3hrgQCLwuUUKgOONTtYDrkcG9XzKLp9eqzlzXs4CJe1qEtuYqAVl1pf67ZlkEkDNzrjv2UqIyiPC/5nSp1yAS9iBVDghgpMzPDS9tSDbjVQqLO8JLqVZOPZpPsnihonrksnIEUHxwDmtFTqHv1Lc7ByJn7RaI7/ooO/dp9UYu9l3KpwNsfD8hhwg/K+QxCC1qg7/YaQhGlZjChCEhKUpmJoJSAMgAObdwjl1IfEH+1djJ5TN/XYl22LwfWFsnkzZx1v5H9ym84QflfIYcIPyvkMQh9SHxB/tXYyeUzf12HUh8Qf7V2MnlM39dhtsXg+sJ5M2cdb+R/cpvOEH5XyGHCD8r5DEIfUh8Qf7V2MnlM39dh1IfEH+1djJ5TN/XYbbF4PrCeTNnHW/kf3KbzhB+V8hhtj4fkMQQYq6uOu4HWLPXNdumPivQaDTgDMTkzNzmykqICUgJnSpaiSAEpBUewIjrvfS1xBp131Bm3MYMXZ+isTCkyM3PXHPS0zMNj2Li20zCg2Tx7IUcsxvzijGnTC9NvWFfSWxBGm67RNtNP5HAdJXr0BzEIw5q9r/rOK+gxg/c9xVB6rV64LNpU/UZ14JDk2+5KtqW4oJAGalEk5AbzGY4vgaiq5FMQTBiuhOxaSOi5IQhEVRSEI4rO6CKIv1WQoIw0wJJIAF0zhJJyA+12ojdcrcjwivt6S9kf/EI7fxxKP6rGsmRqehJh/cLod540m9m5KXIcybDczIzRc2k9k5sN5Hsb+3ELKNFiZUEnnzKZEA5czK7XxxoOVEvBfMNMZ+bdq+NS92fJxtGcg5NOZJQds881vpS896y7z7ke7pLyhHphz7ke7pHyhHpjEfSqzPfmV8mV6Yp0qsz35lPJlemNX3JKcN1L0J4Xtvif6gsu8+5Hu6R8oR6Yc+5Hu6R8oR6YxF0qsz35lPJlemHSqzPfmU8mV6YbjlOG6k8L23xP9QWXefcj3dI+UI9MOfcj3dI+UI9MYi6VWZ78ynkyvTDpVZnvzKeTK9MNxynDdSeF7b4n+oLLvPuR7ukfKEemHPuR7ukfKEemMRdKrM9+ZTyZXph0qsz35lPJlemG45ThupPC9t8T/UFl3n3I93SPlCPTDn3I93SPlCPTGIulVme/Mp5Mr0w6VWZ78ynkyvTDccpw3Unhe2+J/qCy7z7ke7pHyhHpiy9ICqSs1hjNIampV1ZfZySh5Kj7LtAxa/SqzPfmV8mV6Y6y8tHx+zrYnKmupy0wJRIUUJYUkqzUE7jn8MV5WWlWxmlsWpqNCsLVtO2IsnFhxJXNaWmpzgaCixwRkYRU8faikbeuTJA8RhHa2TZNXxIu2nUGgUyoVqtVeYRKSUjIsKfmJt1Z2UtoQkEqUSQAIcykiRGsaXvNAMV6x9WR+LlwJ/uBRPoLUa0a9r7owx+Kpfw0bS6vbDW7MGtCHCu0L4lKdI3Ta9tSVKn5eReLzTCmGw2hO3vClhCUBZSSnbC9klOUas69j2/DH82pfw0ZS3fu9/MO0L5k7L5BsaeLTUZ2IwP0gvC1c1fn4a2Gv65HzLsSk6SOlejR0vqzZCcpjMxSLhTNLn59cyWzTktBsIUE7JCtpxxCd5GW0OOIt9X3+Gthr+uR8y7EtmMujdbGPapQ3HLzUwJSVclW0tvltISt1h45jLeduXbyPYG0OzGIycEUybhBNDne4LnexSybfYEy2RcGxNsuJwwZWvIRVWDhNp7UK6bathVzyb9v1+4ph5gyUs07OsSGzOLlGy8+EBLYcdTsJKgAVZj4Y7y09OnDy9Llt6lSc/Ukzd0FSZAP0x9lDmS1NpJUpOSQpxC0JJ3KUhQHFHJzQisZ24aZUeCq6HabOOzpbRUFpZnFLnFTqUvIG5xDcwtS0JPsTu3jdH0W1obWbadfteoyIrTUxakqJOXAqLgRMtodW62H0jIO7Dji1Jz3ddvB3RsUMTlW59KVvxwuXSoDLfZRjjDIFLyDU4VwoNa6iV0tZ83ndE1M2wmXw6tScnaZP3AJ4Kfln5RoOOuLltkHgSc20lKlLKgM0gHOPsp2nTh9PyiXXJyqySsiFszdMeYeZUl0tOIWhQBSps9coZdaghXEY/e49CqxLsvCt1eflak+K8iY5qkOeDokA/MMhl6ZQyDspfW2AkrHwkAEkx+dN0H7CkpSlNzMpU6q5THKg6X5+fW+9OrnW+CmFPrO9ZKAAOLZyGUQDZ0XCmGnX7KXKfMtxriGFhFT6WqopSgFABoNTXE6vxtjTvw/vNmnmlzFdnnqjLzU02wzRphx5Dctsl0rQEkp3LQRnxhacuOPjf1g2Hoo1GmpZ2szj1cfXLS8oinqbfbWh9EuoOBZSlsBxaRmTxZkR9khoL2RI2y7S1LuKaZfo81QnHJirOuPKlphTSlDbO8FPAthOWQSE5ZHM52jXtAESV2UGYtq45+RkKY+uZe5tmFzL63XX0uTCzmNlzhW0lvZITs7RIJ3JEHGdFKBvx7fiqt4sTKJsMEBhOmgP8ALgCR/Nppgst4IaQtsaQ1BnKla849Ny0jMczu8LLrYWCUhSVBKgCUqSQpKuyDCOeCOA1CwCtl6k0JVTVJuuBSROzq5pTKEpCG2kFXsW0JASlI7A35k5wi/h52YM/HTz6VskgZnaG7rA2zTm1pjoryLQjXCYUYk1W+3rzq5pq8PqPzPTqOlmc+yNKeA21LZP8AxFOZgrG7ZSgdgxpAy2p15CE5bS1BI+MnIRLBrkU7OhjMb/8Azun/ADioijp2+oy/6ZH+YRzvKOA2HO3GucAT0ry3sp2bDlrfOY4nbAHGpreSRQXYAAABTAat3B3EbADBCctm/jTksyE6VUVqXmxMrl5dSQpaCoAJCeEKikbyNpXYyERi67bXmYr2ZpH3fgdYMwcO6Na8+1I1K4qaozFZn0qYbcXwROylhA4XcEHhCUD7IgEiJv1DJjP8n/SPMnrXsQpPDrW949TM6zMvomaoy0kMpSSDzJKnskbshG2Wi58pJhsCppdjefavfuwHklZke02WdOnOhQofml/nUoQASLgaVoAbuhdtok6aGjDou1Fq45uzMUb4xDW4uZmbnrUpKOzBfWTtuNNqmlBtRCiCslTh3kr35RsifVB+EpJJtDEok8Z5jk8z/wDoiO/pnaD3BVfFt8qHTPUHuCq+Lb5Ua+y3p5gzWQKBeyYux7k/EdnxJ6p5j33DkUiHrg7CT3oYleRyf1iHrg7CT3oYleRyf1iI7+meoPcFV8W3yodM9Qe4Kr4tvlRP4xWjwPWpPJxk5x3qPepEPXB2EnvQxK8jk/rEPXB2EnvQxK8jk/rER39M9Qe4Kr4tvlQ6Z6g9wVXxbfKh4xWjwPWnk4yc471HvUiHrg7CT3oYleRyf1iPkr/qhTDKXoc2umWTf8zUkMqMqzNNyjDDruXWpW4l5RSknLMhKiBnkDEfXTPUHuCq+Lb5UW7ipjdSr5s56nSkpPtPLdbcCnUICcknM8RJieHb9oOcGug0GtW81seZPw4LokObziBcKEV5MV+2l1ppXzpnX+azdtQPMkqtznZSZfrZKlNqOew2n+srizcVmtWW85ZAYmJKgRnvPEe1BWRO6CQMiSrZCQSTlnkBviu95dVziraDBhwYYZDFANS9b2rQmpSe1euB70hJqp8k7YtHUzLKfMwWE8xtZJLhAK8u2QM4zjGJtBPD/oT6GGE9r83tVTnBaFLkebGmVMomtiUaHCBCuuSFceR3jsxlmNubgF4In3B01Fc3AuPakIQiKtEgRnCECijF9Vbfi8bR/wAQZH6BPxFO17Uj80fuESseqtvxeFof4gyP0CoRFO17Uj80fuEc4y0+uZzL398lf+H43rlcoQhGlL1EkIQgiQhCCJCEIIkIQgiRauN33qqz+iR84iLqi1cbvvVVn9Ej5xEXEp9eznHaFjbYP/Aj+o7sWtSuMxSCuMwjogXnlI7Wz77reHlZFSoFYq1DqCW1NCap047KPhCh1yQttSVZHsjPfHVRRXEYKV7GvaWPFQdBXrq1dty1G89AjBesVeem6pVqpZFHmpycmnS6/NOrkmlLcWs71KUSSSd5JjVvXs/dGGPxVL+GjZfVkfi5cCf7gUT6CzGtGva+6MMfiqX8NGUt37vf7O0L5kbL4AsaeA1/+wtXNX5v018NP1wPmXYmsa65sZ9qIU9X3+Gvhp+uB8y7E1jPtYixyS+yv9b3BaXsK/dUf+of7WquyIZZCKwjal2VUCQIbA+SKwgiZZxTZGWUVhCiKmyIRWEQoi1R1ySs9DGY+Ct0/wCcVEUdO3VGW/TI/wAwjdrXCYoYl0i/HbQq65AYd1rgKhR+BkxtOLZA20reO/hEuZkpBy2VoPZMaQsuKaeQtPskKCk/GDmI5vlFHbEnbgRmihrznqXk/ZTtKHM2+c1pG1gNNRS8Emo5CCCCvQSogy//AMf9I8ymtdl7en9bzjyLhXJCWTVGeBL7pQnb5klQciCN+We6J5dWxi7iRjrglM3RiCuQcZqM6ecrjEomWW9LJSEqWUpOyUlwKCTuJ2VcY2THnd0yMKmcTdYrpDJem5iV5iv6rkcE2F7W1Pv8efF7GNltqYZEkRFJLWmh5Qvo98mKO60rX3ZKw87Oh1DXilQSDeK6hUa1anOTDH3SieVL5UOcmGPulE8qXyo+XpUpLvtUfJ0emHSpSXfao+To9MaNtktw7+te8dy2l/18L9Pevq5yYY+6UTypfKhzkwx90onlS+VHy9KlJd9qj5Oj0w6VKS77VHydHpiG2S3Dv603LaX/AF8L9Pevq5yYY+6UTypfKhzkwx90onlS+VHy9KlJd9qj5Oj0w6VKS77VHydHphtktw7+tNy2l/18L9Pevq5yYY+6UTypfKi3MVKZZErZzqqEulmoh1sIDD6lr2c+u3EndlHc9KlJd9qj5Oj0xauLeCkvhvbrE81PTUyp2YDBS60lAAKVKz3dndFxKvgOitDIzidRqsfasCfZKRHRZKExtLyM2o5ReseFO/447KzbIrGJN0yFv0CmTtZrdaeTJSEjJtF1+bfc61DaEjeVEn/XiBi7dGvRgvnS7xXp1lYf29PXBXqitKQhlB4GUQTkX33MtlllPGpasgOIZkgH0b6qDUs2Nq3rZTWJ4U+88U5sKM1cjsnsCnoUkpMvJJUSplvZJCl57buZ2sk7KE7hLSjo51DWvPeWuXklYMEsd58Zw81g7Xah1nQNK2vwOt+ctPBq0qXUGeZp+m0WSlZhoqCi243LtoWnMEg5KSRmCRui6o4oGUco2YLxs95e4uOlIQhBSpCEIgUUYvqrc5au+0SfCDI/QKhERLeKFuJbSDXKYCEj/jDtRLr6q4AOrvtEHi6IEj9AqEQyN6M1AU2DzTVuuAJ+zI7X5saHlW2XMdu3Ei7QvdPyaotoMyfibhY13nmucaK6Oilbff2meOEOilbXf2meOEWz0slvd0VXxqOTDpZbe7oqvjUcmNT2uR37ugL0dum3+Bh9J71c3RStrv7TPHCHRStrv7TPHCLZ6WW3u6Kr41HJh0stvd0VXxqOTENrkd+7oCbpt/gYfSe9XN0Ura7+0zxwh0Ura7+0zxwi2ellt7uiq+NRyYdLLb3dFV8ajkw2uR37ugJum3+Bh9J71c3RStrv7TPHCHRStrv7TPHCLZ6WW3u6Kr41HJh0stvd0VXxqOTDa5Hfu6Am6bf4GH0nvVzdFK2u/tM8cIdFK2u/tM8cItnpZbe7oqvjUcmHSy293RVfGo5MNrkd+7oCbpt/gYfSe9XN0Ura7+0zxwi28XMQKHWcN6rLSlWkJmYebSENtuhSl9ek7h8QMcellt/uiq+NRyY6TEfAijWlZFQqMs/UVvyiEqQHHUlJJWlO8BI7BMV5aHJba3Nc6tRo5VZWlMW4ZSKIsFgbmurQmtKX6Vh5XGYRVXsopG5rjqQUnd8cVG49nfGSNFrRLv7TMxWlLNw7t2duCszGyt7g0lMvT2StKDMTLvsWmUlQzUfiAJyERAJNBereamoMtCdHjuDWNvJJoAvUnqx8jq5sCR/6Aon0FqNaNez90YY/FUv4aNuNCrAip6MGijh/h5WK8bmqNm0SXpTtS4AMJf4JOykIQOJCE7KE59cUoSVdcTGo2vY9vwx+Kpfw0ZO3a+D311DtC+Zuy69j7FnXQzUE3HWM8X33rV3V9/hr4afrgfMuxNYz7WIhT1fY/prYa/rgfMuxNWz7WPiiwyS+zP8AW9wWm7Cv3VH/AKh/taucIQja12RIQhBEhCEESEIQRYk0ydF+naV+Ck/bsyWpeqNfbNJnVozMlNJB2Sct+wr2KwONJPZAiLvRI0Lq1j9pKvWXVpOYpspbMwo3Pteyk221bKmMwfZuKGwMjnslSxmExM8UgjeAY/CVpMrJTUw+zLS7T02oLfcQ2EqeIGQKiBmogbt/YjET1jQZqMyM/wDDjyjUVo2UeQcjbE9Ano9xZ6Qp6bReAeY6dVRqp+NBocpbNFlKdIS7UpIyLKGJdlsZIabQkJSkDtAACIxcfPUxVq47Y+XtfjmL96UeavauTlcflJWmSxbl1TD63i2FbQKkpKyATvyESkZQjIxYEOI3NeKhdPsW3Z+yIhiWbFMMkUqKYalEb604tXw54g+bZblw9acWr4c8QfNsty4lyhFv4Nlt4FsflMyo46/q7lEb604tXw54g+bZblw9acWr4c8QfNsty4lyhDwbLbwJ5TMqOOv6u5RG+tOLV8OeIPm2W5cPWnFq+HPEHzbLcuJcoQ8Gy28CeUzKjjr+ruURnrTi1fDniB5tluXHW3d6kroVRt99um453WKgciwqpUViYlkKzGZUhDqVHrcwMlDeR2N0TCwyEPB0tvAojZMyorfOOPPQ+5YG0CNXbhvq68KDbNh01YmJ/g3axWJtQcn6y8hGyHHV9hI3lLSckI2lbIzJJzwiKkZwyyi8a0NGa0XLTJqajTMZ0xMOLnuvJN5KQhCJlQSEIQRIQhECijF9Vbfi8LQ/xBkfoFQiKdr2pH5o/cIlY9Vbfi8LQ/xBkfoFQiKdr2pH5o/cI5zlp9czmXv75K/8PxvXK5QhCNKXqJIQhBEhCEESEIQRIQhBEi1cbvvVVn9Ej5xEXVFq43feqrP6JHziIuJT69nOO0LG2z9gj+o7sWtKuMxUbu1FdnM/GY201WmqKvzWZ3tMO05xNtWBQ30NVm45plS0BRKSqWlU8T0zsEqyJCEDZKyNpKVdGZDc85rRevNNq2rK2bLOm5x4axuk9gGknQFYWr71et+axXHKVtCzZUy8iyUvVquzDSjJUOWzGbjhGW04eJDQO0tXaSFLT6X9ATV44dauvB42pYUg+p2ddExVqxPlDlRrDwzCVPLSlI2UAkIQkBKATkMypSr00XNF2y9DzBij2FYdHao9vUdvJKAdt6adOXCTD6+Nx5Z3qWePiGQAAyHsgdiNhlJNsIVN7tfcvJOXOyBNW/F2qHVkBuDdfK7l1DAcpvXBwbuxlEPWsX0gr9xWxbXbN+Uij0eZsmcmmpRqRacG208UKSsrUo8IC2hshQCeNW7sCYjZB7AjT7WwaGxxrw5/24t+T4W6rUYVw7bSM3KlIglSkZAEqW2SVoG7cXB2RFjb8tFjShEIm68jWPi9ee9kux52fsZ7ZJ5q28tH4wMRzjEdCjUwjxRqGCmJlFu2lIk3Kjb8xzXLpmkFbBUEqT14BSSMlHiIianRUvi6sSsBLbr16UuVo1xVWW5omJSXbW2hCFKPBHYWSpClN7Cikk5EkRGjqvtD0aSOLYuCsS6XrNtJ5D0wFAFuoTW5bUvv40jrVrBBBTkn+tEt7KEto2UgJA3AAZRjslpWMyE6K8+acB71q2w3Y87AlIk7FcRCefNboOFX9VB7eRc4QhG2LtaQhCCJCEIIkIQgiQhCCJCEIIkIQgiQhCCJCEIIkIQgiQhCCJCEIIkIQgiRRRyEViixnlBCoxPVXTyGdXdaRWttA6IMjkVKAB+0Kh24hmRpP2qlCQXZzMAA/YkZcX58esK4LZp11SqZap0+RqMuhfCJbmpdDyEqAICgFAjPInf8MdOcGLQ2v+6tt+a5fkRgrVsKFPPD4jqUXctjbZtmckLPdIS8uH1JJJNPZSh7V5VemhtT3ad8Ujlw6aG1Pdp3xSOXHqr6DFn+9W2/NbHIivQYs/3q215rY5EYrxMlt+V0X52FpcTb+b/VeVPpobU92nfFI5cOmhtT3ad8Ujlx6rOgxZ/vVtrzWxyIdBiz/erbXmtjkQ8TJbflPnYWlxNv5v8AVeVPpobU92nfFI5cOmhtT3ad8Ujlx6rOgxZ/vVtrzWxyIdBiz/erbXmtjkQ8TJbflPnYWlxNv5v9V5U+mhtT3ad8Ujlw6aG1Pdp3xSOXHqs6DFn+9W2vNbHIh0GLP96ttea2ORDxMlt+U+dhaXE2/m/1XlT6aG1Pdp3xSOXDpobU92nfFI5ceqzoMWf71ba81sciHQYs/wB6ttea2ORDxMlt+U+dhaXE2/m/1XlT6aG1Pdp3xSOXHSYk4/25dlj1GnSj0xzTNISlG2hKU7lpO87W7cDHrH6DFn+9W2vNbHIjgvBi0Ejdatt8feuX5ETQ8kJdjg8PNR8a1RmPlU2hGhOgvk20cCD52g3b1ed/U96je6tPi4aTe97Sk5bmDCHVOqnOEDU3cvBqyLEoPZJaJzCpjLZACkoJVvT6KMK8JLZwSsKnWvaFCpVtW9SG+Ck6fTpZMvLsJJJOSUjLMkkkneokkkkkx3dNp7FKpzMtKstS0vLoDbTTSAhDaQMglKRuAA4gI/dBzzjZ5aWbBFG4rz1lbllPW/M7bMHNYPRYMB3nWeiguVQkCKwhFytTSOK2kuJKVAFKuMEZgxyhBF1Vo2NRrAo4p1CpVPo0glxTolpGXQw0FqO0pWykAZknMnsx2gTlFYRAAAUClYxrQGtFAEhCERUyQhCCJCEIIkIQgi//2Q=="}},"execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\nprint(slice_thickness)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"#lets create a function for above task\ndef load_scan(path):\n    slices = [pydicom.read_file(path +'/'+ s) for s in os.listdir(path)]\n    slices.sort(key = lambda x: int(x.InstanceNumber))\n    try:\n        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\n    except:\n        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n        \n    for s in slices:\n        s.SliceThickness = slice_thickness\n        \n    return slices \n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The unit of measurement in CT scans is the Hounsfield Unit (HU), which is a measure of radiodensity. CT scanners are carefully calibrated to accurately measure this. From Wikipedia:\n\n![4rlyReh.png](attachment:4rlyReh.png)\n\nBy default however, the returned values are not in this unit. Let's fix this.\n\nSome scanners have cylindrical scanning bounds, but the output image is square. The pixels that fall outside of these bounds get the fixed value -2000. The first step is setting these values to 0, which currently corresponds to air. Next, let's go back to HU units, by multiplying with the rescale slope and adding the intercept (which are conveniently stored in the metadata of the scans!).","attachments":{"4rlyReh.png":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAdgAAAGgCAYAAAAJqHY1AAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAC3oSURBVHhe7d3NauRMf/fxXs3x3EGLOYI5iiTgBsMcRwIxuPFRPJiQJzDgnqV33pgQyMa4vTLMPniR5fVkU4+qSiXVm1SSWypLXd+GD7dbpZeS+q/6dcme+9oJXrx48eLFi9fsLwKWFy9evHjxWuBFwPLixYsXL14LvAhYXrx48eLFa4EXAcuLFy9evHgt8CJgefHixYsXrwVevQH7f/7vEQAAjBB7DQbsf3/8D/Ap1A+AUhCwyIr6AVAKAhZZUT8ASkHAIivqB0ApCFhkRf0AKAUBi6yoHwClIGCRFfUDoBQELLKifgCUgoBFVtQPgFJcUMD+u/j73U7sjL/9i/iP6Hp9/kv809/ktv8o/jXaPt6//sNO/P2/xdtKR8Bifube9+9d/54eux4wjwsJWC9cjX/498i6fea5yWS4ymMTsHEELOZHwGKdLiNg//NfxN/kjdPOWv0bKXYDNeu021jr/Ns/1v8rf7bbG+ZYLf8YFhPw9v6UvxP/9J9mnz3Hdfpqr9fX/j/iP/7576z2dYY8AYv5EbBYp8sK2N5Qid1AfQEbEQR3rD2yvQrYoW0Sx21n4H3rdOdjZs6+v/3zfzX7WAcCFvMjYLFOlxGw0QCyb5bYDdQfsF0omRuymXEGM+VQ+hGxf5NHjmtmsuY47ReIrv/OcaL9WuegQcBifuae6mPuAQIWeV1IwGrhLM7cMFMC1r3J3MCM3chD69v8bf2+2Y+Nvb41gds7Gw0eQdvs/X49Ahbzi92XNnOvmfXce5aAxVIuKmBbwSPjuQK2Zu2704XYuPUlc5wRfTMB2vdHWwQsijY2OPvW61sOnOciArb94x4rgNygMzdaF3ztNkHADjwibvZt8wPVf2+O0+4zeNw7ImCDbbz9mvb2XNaLgMX8xgZs7B6v+b+SAWZyEQGbniV2wReIBGzvOlOPIwO/d3bp3/T24OAFbG/fppzfOhCwmN/YgA3/0t7W+ysY4JMuI2ClIPz8Wae5CbW//7e+EKtvRjsU/YBKHcdub7a1w0/exOYmd2fXQwGruSEazqqDkF1ZuEoELOY3PmCVyJdewhVLuJyAxSZQPwBKQcAiK+oHQCkIWGRF/QAoBQGLrKgfAKUgYJEV9QOgFAQssqJ+AJSCgEVW1A+AUhCwyIr6AVAKAhZZUT8ASvGpgAUAAGmx12DA/vX//hf4FOoHQCkIWGRF/QAoBQGLrKgfAKUgYJEV9QOgFAQssqJ+AJSCgEVW1A+AUhCwyIr6AVAKAhZZUT8ASkHAIivqB0ApFgvY0/2j+Pbjl7h+8tv+iOsfj+Lu3V+OEswXsG/ittqJq2N8+W6nVYe3mdsBYJyFAlaG6C/x/Wcdsj9fxCm6Dko0T8B2IegGbLN8b45xFFdOSJ7bDgDjLROwT8/17PVZ/H5/Ed+D2ao9g/0Qd3UI393L9X/VmNleunMD9vVQ6dlltRdXdRg6AXs6iGpXiduTtey4r9c9iNc52gFggkUC9vdNPXu9/6h/lgFqfjbtfsDWwcostxjnB+xBPKif9WzTCdhYGKrQ3Ottzm0HgAnmD1hv1qp/F1vPZtt1woB1AxiX7NyA7YQBq2a30YDUs9Jz29tlADDC7AGrAtWZkerfx3Z/7BQGbPiHULhUg/UjZ5Dy8W8gFnAELIB1mzlgdZjq36d62tAlYEuW+oI2Ho+IAazbvAFr/rgputyEKgFbskUDNhaGdmie2w4AE8wasPKPm77d/Im02b9rJWBLtmjANsv4ZzoA1mC+gI3+k5xO98dOBGzJlg3Ybrn5/W0Yjue2A8A48wUsMAL1A6AUBCyyon4AlIKARVbUD4BSELDIivoBUAoCFllRPwBKQcAiK+oHQCkIWGRF/QAoBQGLrKgfAKUgYJEV9QOgFJ8KWAAAkBZ7MYPFIqgfAKUgYJEV9QOgFAQssqJ+AJSCgEVW1A+AUhCwyIr6AVAKAhZZUT8ASkHAIivqB0ApCFhkRf0AKAUBi6yoHwClmDlgP8Tdz1/i+/1HpA2YIWBPB1HtdmJnqw7itV3nTdxWXVt1eHO3T7YDwDwIWGR1dsAe92K379tHE55t+1FcOSGaageA+WQO2D/i+scvcf1kLXt6Ft9+PIvfbfujuHt6Ed/r9b4p9ft3a/1mH23bvb091u7cgH3Y78TVMd6mZ7eVuD1Zy2Qgmxluqh0AZrTCgJXB2QXm75v6/c8XcbLa2/2/myAmYLfivICVM9BKVNYjXiccY2GpQnUvHsa0A8CMVhmwve3y5zZstdP9IwG7IecFbOSRrgzIJjRfD1VPgOpZa6q9XQYAM9hUwKowvfnTtXntznKs0mD9yBmmmZk6hgJQh658bEzAAlgTAhZZDdfPZ+g/XFK/l+URMYAV2VTAqp95RLxpZwWsmuH6YShnsM0MNBaWdqim2gFgRpkDVrd3s9Dm/diADfav1ydgt+OsgLVnq2aZE5D8Mx0A67FIwOp/QuNpQ9WEomT+Sc7YgO3WabeX/0zHm9Vivc4LWEmHYvv72WD22YRo0x6GZ6odAOYxc8B+gchjY6zX6uoHABayrYBVs1n7/3gi9Ugaa0PAAijFtgK2pv+oyTwiJly3hoAFUIrNBSy2jfoBUAoCFllRPwBKQcAiK+oHQCkIWGRF/QAoBQGLrKgfAKUgYJEV9QOgFJ8KWAAAkBZ7MYPFIqgfAKUgYJEV9QOgFAQssqJ+AJSCgEVW1A+AUhCwyIr6AVAKAhZZUT8ASkHAIivqB0ApCFhkRf0AKMVqA/b0zn/n9RJtMWBfT2/R5QAwZOaA/RB3P7v/GLrv+im2jU/vg/+Q+mWaL2DfxG21E1fH+PLdTqsOfjim2m163eF1znQ6iKrpS6s6iNd2nSn9BbAmiwTseeFIwF6yeQK2Cx03YJvle3OMo7hyQinV7ssQsMe91R/f1P4CWJOvCdj3F/E9Orv1ZsA3f8JtsWnnBuzroWpmeXtxVYePE7BqNliJ25O1TAaYmRGm2h1diCttyHnLo9uO97CPzcIbk/oLYG2+IGD/iOs2ULXfNzJQn8Vv9Z4Z7CU7P2AP4kH9rIPOCadY+KiQ2uttUu0Bfwbrzyh1QH4+8OT+KlH1Bfbk/gJYk0UC1p6ZGoOB+fRcr0PAluDcgO2EAatmt9FA0rPAVHu7rOUFbDTc9GPb3lnooMgjX3mMpo/T+wtgTRYJ2HHhqGeyXQgTsCUYrB85YzMzOUcsUL4gYKOPZ711bJPOx+gCm4AFtu0LArYL1nY9ZrDFGK6fKcKATT5SnfzI9cyA/RTrvCb3F8Ca5A9YGaY/X8TJX0bAFmHRgI2Fjx1SqfaAF57RcDvjEbE8dnR/zQx1cn8BrMnXBGwbprX2L4ofxd27XEbAXrJFA7ZZ1v/PWlLtPi9gg+3n+COnoVn41P4CWJP8AVvTfzVs/+7V/cvi0/2jbuOf6VycZQO2W25+3xmGUardpX4PKtdtQ87d/vzZpA7N/v1N6y+A9Zg5YIFh1A+AUhCwyIr6AVAKAhZZUT8ASkHAIivqB0ApCFhkRf0AKAUBi6yoHwClIGCRFfUDoBQELLKifgCUgoBFVtQPgFJ8KmABAEBa7MUMFougfgCUgoBFVtQPgFIQsMiK+gFQCgIWWVE/AEpBwCIr6gdAKQhYZEX9ACgFAYusqB8ApSBgkRX1A6AUmwjY0/tHdDm2h4AFUIpZA/b3zS/x7eaPu/zpWXz78Ut8v3dD8nT/KL79fBEna1noQ9z9DLfFdqUD9k3cVjux22nV4S2yDgCs36wBGwtNGbrXN3XIOsE7NjgJ2EszHLBNuO7NOkdxRcgC2KhZA/av9xfx/cez+N0ukwH5KO7e/4hrZ7l8L5ebbeqZr+X6yWxrLW8D2lvuBHpzvHs9a/5mjtG246sN1s/pIKpdJW5P1rLjXuyqg3i11wOADZg3YJvw0wFZk+GpAjCyXAWuDFpreU09Zm7DWG/XzWD99/6sWbenHz3jqwzWTyxMVejuxYO9DAA2YOaA1QHZBuDTc/uzDEL75+B3tYb6nW1PwDpthh3eYQBjXYbq5/VQ9QSsN6sFgA2YPWBVCDbhqX7/GsxmveWKnsm2j317AlYFs7NeR69jhy3WiIAFUIr5A9Z5/GvPNs17GYLmd6NdsMZnqZGAHXz8S8Cu3WD98IgYwAWZP2BNgD51M1bTpmau99ZyGaZ+YE5+RGwjYNdusH5iYcofOQHYqAUCtvk9bB2ybTA21O9h7eV+YLZ/UWxmuF7ANu+dUFb7cNcnYNdruH74ZzoALsciAatDLxJ0TYDay2UYd79LlWHr/mVx+3vXvn+m4/xTHAJ27dL1w//RBIDLsEzAAj2oHwClIGCRFfUDoBQELLKifgCUgoBFVtQPgFIQsMiK+gFQCgIWWVE/AEpBwCIr6gdAKQhYZEX9ACgFAYusqB8ApfhUwAIAgLTYixksFkH9ACgFAYusqB8ApSBgkRX1A6AUBCyyon4AlIKARVbUD4BSELDIivoBUAoCFllRPwBKQcAiK+oHQCkIWGRF/QAoxUIB+0dc//glvlm+338E6/2+cdf5dvPHav8Qdz+9dsv1k7svbMP4gH0Tt9VOXB295aeDqHY7sbNVB/FqrwMAKzB/wL6/iO9BADaBawWoCtdYoLbL9PtYMGO7xgWsDlcZnkHAHvdit2cWDGD9Zg7YgVB8eq5nno/i7l2+l4FrfvbXeRa/1XsC9hKlAvb1UDWz0r24isxgH/aR0AWAFZo3YNXsNRKcAX+2GkPAXqJ0wB7Eg/o59ohYLqtE1cxueTwMYM3mDVhnBpri/Z7254s4Oe1NCNvrNAjd7UoFbCcWsEdxVYdqdXjrlsnfyRKyAFboCwPWZoUtv4O9aOcFbIwOXR4bA1ibeQN28BHxhzilHh072xOwl2j+gB27HgDkNW/ADoWi/dfFvTNdPZPVf4FMwF6iswJW/gXxbt/8jtaQM9hK3J7sZQDw9WYO2FoTpG4w+v9Mp/n9qvd719P9oxW8BOwlOitg+0KX38ECWKH5A1bx/oApCFwt+D+acAKXgL1E5wWspH/nyl8RA1i7hQIWiKN+AJSCgEVW1A+AUhCwyIr6AVAKAhZZUT8ASkHAIivqB0ApCFhkRf0AKAUBi6yoHwClIGCRFfUDoBQELLKifgCU4lMBCwAA0mIvZrBYBPUDoBQELLKifgCUgoBFVtQPgFIQsMiK+gFQCgIWWVE/AEpBwCIr6gdAKQhYZEX9ACgFAYusqB8ApSBgkRX1A6AUMwfsH3H945e4fvKWv7+I7/Xybzd/mnUexd27t47x9Cy+/XgWv2Nt2Lx0wL6J22ondrvG3l/fba8Ob147AKxDhoDtCd0+BOxFSwXsw94OVR2mXYg24dq2H8UVIQtgpRYO2L7AtWewep1vSr383g7YZt2nZgZs1nFmvx/i7qdpq/18Eaem7fdN/V7Nmrv1T/ePwTLkM1g/p4OodpW4PXXLXg+V2FUH8drT/tdx37UDwIosGLA6+MKZqx2wep3v9x9WmwxKO2Dt901otiHqb98EqGkPZsN9fUIuw/Xj82awsTBVobsXD/YyAFiBxQJWBWEbtv46TcDGHgc7y7r9Rdujj5PtEPW2V78L9tdHTuMC1vo9qxWozmzWiM1qAWAFFglYE6w6ZP1A6wLWmW0aTggOB6zavjmez8xq7cfEPB7+euMC1iJnrc0MlYAFsCWLBGzwyNcJtZkD1t/e1+6Px8NrMDlgm9ns1bH+mUfEADZkkYANA9Fe1gVs9BGvs2w4YKPbB5pgvefx8BoM1o81W+2W678UVgEbC1P+yAnASi0fsDX9KNcOTfePnLoZrt5+dMCa7e1ZrGo3+9faR8k8Hv5yw/VjzVbNMidA+Wc6ALYjS8C2QaiC0Q5Yu02qlwf/TGcoYCV7+2YfVrgqzf/RRdgv5DZcP5L1B07eHznF2glXAGs1c8CuFH89vBqbrB8A+IQiAlY+Irb/rSy+DgELoBSXHbDNo+HkXxojGwIWQCkuO2CxOtQPgFIQsMiK+gFQCgIWWVE/AEpBwCIr6gdAKQhYZEX9ACgFAYusqB8ApfhUwAIAgLTYixksFkH9ACgFAYusqB8ApSBgkRX1A6AUBCyyon4AlIKARVbUD4BSELDIivoBUAoCFllRPwBKQcAiK+oHQCkIWGS1xvp5Pb1FlwPAOWYO2D/iWv4HzgPP4rezzqO4e7e3m9HTs3c8rMn5Afsmbqud2O0ae39/bnt1GApPve7wOjM6HUS1q8TtyV9mnY9UHcSrvR2ATVokYK+f3OW/b+qQ/fkiTu06BGypzg3Yh30dQG2o+gHZhGvbfhRXgyGbM2BN8HsBe9xHviQAuARZAtYNPT9gP8TdT2u22wbx2HZ9TN1e7/eegF2zsQHrBmkjMgN8PVTdjC82Q5QBFp0RmsBreKE9dTYZ7a9F9rPa74P+ye2uju66AC5DloCVM9humR2wTXje/HHW7UJ0XPv3+4+m3YQtAbtWZwVswJuBxsJUhe5ePNjLWv4MtglX67iqHyNCdrC/sg9yH8EXAHm8Ong/EegA1m+RgG1nmxY3BJuAfX8R34MwtEI61R57HMwj4lWbJ2CtWaYVSM5s1ojNaltewEbDWD9mTs0y+/urQ1QdP+hL5BG2CeN2ewBbtUjA+jPY0/1jHXpm1moFrAzDnkfCKpAT7Wq/fns0lLEWg/UjZ6BmJufpDTi1jQ7FswM2+jjZW8c2or/q0bAT4H19McYFOoD1yxKw7nICtmTD9dMZ94hY0gGoAuncR8RTA9YS7a8/Gx0VsNb5RNsBbMXXBiyPiItzVsBas9VuuTXji4VpNDQNLzyjYfz5R8RqRt3MaH3tF4Lo8VIhDGALsgSsO9O0AraZjZ77R05duz4+AbteZwVsbHbnBKhu77aL/I7T4c9O/e2bfvQGdCfeX08wg02dD4AtWyRgwz9ysgPPDlipCUmzbs8j4XHt9X75ZzqrNjZg+zUhaGaDQRi57alHu+0ssw3H1P7PEH1ErL8ELHI8AF9q5oAFhlE/AEpBwCIr6gdAKQhYZEX9ACgFAYusqB8ApSBgkRX1A6AUBCyyon4AlIKARVbUD4BSELDIivoBUAoCFllRPwBK8amABQAAabEXM1gsgvoBUAoCFllRPwBKQcAiK+oHQCkIWGRF/QAoBQGLrKgfAKUgYJEV9QOgFAQssqJ+AJSCgEVW1A+AUhCwyGqN9fN6eosuB4BzzBywf8T1j1/i249n8Tto+xB3P/va5iT78Cju3mNt+GrzBeybuK124uoYX77badVhKDz1usPrzOh0ENWuErcne/mU/gLYkgUC9lF8r4P0+slre3+pl9dtBGzR5gnYLpTcgG2W780xjuJqMLRyBqzpsx2wU/sLYEsWCdi7+2fx7eaP03a6fxTX93XItgGrZ7tOED/V21kBLLf5pmbEkheaMrBN288XcTLLg4A1M+fYushtbMA+7O3g6bweKj3bq/biqg4nJ2BjM8Tjvl73IF7N+1YX0kp7LG95dNtQX38N2e9qv3f7N6m/ALZmmYB9l/9rz1RlyNXvVSiODFgvbN339rY6QL/ffzT7sQPWb2tCm5D9MucH7EE8qJ91EDoBGwsnFWL7Zhuf3kc3Y/RnlE0/RgTeYMDKPsh9+IE6ub8AtmShgP1f8fvGCk8ZrHJGe07A2oba7ICNrqdD1zkusjk3YDthwKrZbTSwvFliywvYaLjpx7bh73pd/f2Vx2iO7/Vlen8BbMliAStnimbm2P48JWC9R7v2esOzULcP7aNhjz2rRT6D9SNndObRrCcMuAUCNvp41lvHNqK/6tGwE+AELFCKxQJW/1GTDEEZlNay0QHbkbNhHYxWcI4NWB4Hr8pw/XQ+M4M9+xHx1IC1RPsrj23vzw/Pyf0FsCXLBawJ1icTtPWyTwasZv0+9axHxPhKiwZsLJyioWl44RkNt88/IlYzVG9ma6j9Te4vgC1ZMGC7R7Tt41gnYJtHwO1fG5tHwk27H45qW7NvN5zVcdr9+CFf79Oexar9dn1EXosGbLOs227qP9Pxt2/6MSLw0v2t+TPYyf0FsCWLBqwbiua9PaPUQdk+/pWzXau9ezSsObNdtS/T5u/TDlAT3Abh+pXGBmxaLGC75WammAqrdpbZhpy7/ayzySBgpWn9BbAdMwcsMIz6AVAKAhZZUT8ASkHAIivqB0ApCFhkRf0AKAUBi6yoHwClIGCRFfUDoBQELLKifgCUgoBFVtQPgFIQsMiK+gFQik8FLAAASIu9mMFiEdQPgFIQsMiK+gFQCgIWWVE/AEpBwCIr6gdAKQhYZEX9ACgFAYusqB8ApSBgkRX1A6AUBCyyon4AlIKAneD19Db4HmlrrB8+RwBLWCxgH/Y7sdvZKnF7iq+7fm/ittqJ6mAGYv89xjo/YPW1b+tqb+3vuLfqzXV1tPdhZP4cTwdR+feBWub1tzqIV3s7AJu0QMAexZUcJOyBT2oGv/hAt3YE7FzODVj1xa2trfTnoNbvDaycn6M+VvBFU94X/r0C4CLMHrDuAOh6PVTWYCcHnHqwOZhZhxl4zEDUsAbH2L7VPmPHM7OFdv/NtvYsx94uMpPQXwa8/uwP3nuzj/5+959recYGbLSOIjNAt6Y86rPei4dYW/C5jvkc+w3VvST7We33Qf/ldtv80gkgZeaAlbPXseHRDGRBELkzCmcADQZMvX50gDKBaQY9E6zmvTNY61m3vR81YLbH8vuVeh/7MuGfa5nOCthAeN3HtRnxz9E+rurHiM9tsL+y1uQ+gi8I8nh18MpjTgx0AOs3b8CqAaRvxuCLDIDRGYdeT4efF4RDx1Ntdmj6IRqGqsPpS3wgbt8n+x0510LNE7D6eg4GUvQz8XmfS7SeEnXS6O+vPEYTqkHA6n07dSHXIWSBi7B8wDZB135DbwcYO4A0Netz1u2YQcgeyNT6Q7OGYDALBzd34NTLuuOOC9h0v8NzLdVg/ZinDBG9127gy40TXFHeenJfQbgN7GtEf2VttNsGNRkzLtABrN+8ARuEmMcZYMLQcR+r9mhDPBFakwK2C1ZnsJ0SsIP9TvS1IMP10xmewdoi13ZUkEne5zo1YC3R/sp+2Psb1S9qBbgUMwdsYmB0BpjIQBKdjfia7Q5yXwPrTgnY2MA6IWDT/WbQNM4K2Oh1tr8oNWKfZ5T3OaqaGbH/iFh/h55stHUXPd6YLwcA1m72gDUDUjDAqcFLDi4DAdssc7ZVg5A74LQD19AMZ2rA2gNdT197AzbZ79i5lumsgI1dx0iYDn7Jc/R8jta2al8jwnrUMYOaHHc+ALZpgYDVYt/eY2Eahk4zyLXb2QNSownAwcCaErD1ezVAtseUYeu2+6EehvxQv/vOtTxj66efd50jYSQ/yy40hyU/xznDLqhJSdfZIscD8KUWC9hFqYHKf7SGLVhF/QBABpsMWDnrGDtDwboQsABKsa2AbR4N8xhtuwhYAKXYVsBi86gfAKUgYJEV9QOgFAQssqJ+AJSCgEVW1A+AUhCwyIr6AVAKAhZZUT8ASvGpgAUAAGmxFzNYLIL6AVAKAhZZUT8ASkHAIivqB0ApCFhkRf0AKAUBi6yoHwClIGCRFfUDoBQELLKifgCUgoBFVtQPgFJsMmBfT+v5j637fVlT39ZojQHLZwZgCbMG7MN+J3Z7r/24F7vdTlQHL4gOVfMfTj+Kq7r96mht45Dtlbg9yZ/fxG0V7utr+H1ZU9/W6/yA1fUia0rbiwenXX8Opn3488jxmbn9Ce6PSf0FsCWzBmwXmt0yGbpX+zpknYHFHthSAWvLMSCO5fdlTX1br/MCVteKfY1VzbUh24RVW2vh+q7lPzP3S2e8Zsb3F8CWzBqwf50OonJmFHIAkbNPOXDYy+V7MyvVg8rVQW5rvsmbNnvdZjAy63iDVrvcC3iH6l+9r4OeVbf7aWbZ7Xtn/WZ5Q38R8Pty+ETfmmvT9sU+58s1NmD7n4b4M1arlszna19HuU20JmaoJ0u0v5H+OF9CJ/UXwNbMG7DN4NTORuUAogaLyPJ2oNQBaw+carBqBxk7jPV+/BmA/Y0/NotuqeNaA6EJVvPeGfCa4Ldm1qpfbT9TfUn1TbeXNpieFbAx9mcWCyen1nzxz8w+rluL/cb11zve5P4C2JKZA1YPNPYAYn6W4WL/3A1GYZDp4LMDuCdgnfUMvY6zP6MJ2K7NP3akLzbnePHBeXzf/O3LMG/A6mto1ot+uXK+NPm8zyAabomaaAz3t+lnvR+7f9P7C2BLZg9YFSzNQCMHnXZgkgNHM5g4y88IWDVAyUErIhpcweBl79u89wdTvazb97iATfdNr58auC/NYP2ozz1+zcLr1ISWFWpnB6w8vr998DlbJvW3YdU2AQtctvkDVg0QcgCRwWTPBsx7OWAlQm1KwAYD4oBJAat/loOlMwC3/fIH3ql90+sTsHGDM0L1OUZCLxaQbT1ay1reZxjbPvic48bNuCXrc5/cXwBbMn/AqgGkDq1jPVB4g4eauco/ZnKWfz5g3fVGmBKwscHPOZ4/8E7tmzXQRtsv09kBq66r/ZlZYuEU+xxb3mcWDbdIfUZE+xutAWt/k/sLYEsWCFg92FR1yPrf+uWsLlweGcCcgckOwXioOQNScgCeELD24Ke2rY+V6Mv4vul2AjYuGljB5+drrnm7nf48/Trs9Hxm1nFVP0YEXvwLQeQzdgJ0an8BbMkiAauDJBIeTUi5y61QM8uccHNDUD16rdfvBqVmkJLLlIEBeErA1u/VoNnuV/bHbff7Mq1vkcG3AGMDNqa9vhHddXSveSqskp/Z2bPJ1P6m9RfAdiwTsEAP6gdAKQhYZEX9ACgFAYusqB8ApSBgkRX1A6AUBCyyon4AlIKARVbUD4BSELDIivoBUAoCFllRPwBKQcAiK+oHQCk+FbAAACAt9mIGi0VQPwBKQcAiK+oHQCkIWGRF/QAoBQGLrKgfAKUgYJEV9QOgFAQssqJ+AJSCgEVW1A+AUhCwyIr6AVAKAhZZrbF+Xk9v0eUAcI6ZA/ZD3P38Jb7ff0TapD/i+sejuHuPtaEEswbs6SCqXSVuT/byN3Fb7cRup1WHofDU6w6vc76Hfdef3Vn9BbAlmQMWpZsvYE0w2YHVLNubYxzF1WBoLR+wMlyd/R/3dZ/34kG9n9pfAFuSOWC7Gezvm1/i280fp91dpvf17Ufj54s4tevKtno/989NO7PirRgbsGrW1wZP6PVQiWq/d2ewsRmtDLTqIF7N+5Y7c+yO5S2PbhuK9jc6w7ZM6i+ArfmygP3rSYbjs/jttP0S10/y53A/p/tHK2Sb8HVCF1swS8DKYJIh5AdULJzUOmbG6PNnsP6MsunHiMCL9jcVlpP7C2BLvi5gnUCt2YEbhK+k990XwNiG8wNWhmATql7AylltPLD6ZpFewEbDTT+2vTray0Kx/pr+PMj/NTNia//T+wtgS74wYN1HwvbParZqHg179L7tsMWWDNaPnNG1QeQyAaceDTuBOGPARmec/izXkuiv6k/9s7Ot2kaHLAELXLYvDdhupurOZt3HwTEE7FYN108nOoOV4WMHkh9Gkx+5nhmwlqEZbGx/6gvD5P4C2JKvDVgTrPcv4rv9SDj6iNhGwG7VOQFrZoQxKrBi4RQNTcMLz2i4ff4R8VBgf66/ALbkiwPWehzs/EWx3o8zi1Wha7YlYLfqnIANqICyH6fq8Oq2m/rPdPztm36MCLx4f60wNctkgLahOrW/ALZkkYD1f2/aBWgYsH+9y9lrLCz9fdnbEbBbNTZgRwkCVmpCq5nZpsKqnRW3Ieduf/5s0ttfMEOe1l8A2zFzwALDqB8ApSBgkRX1A6AUBCyyon4AlIKARVbUD4BSELDIivoBUAoCFllRPwBKQcAiK+oHQCkIWGRF/QAoBQGLrKgfAKX4VMACAIC02IsZLBZB/QAoBQGLrKgfAKUgYJEV9QOgFAQssqJ+AJSCgEVW1A+AUhCwyIr6AVAKAhZZUT8ASkHAIivqB0ApNhOwr6e36HJMt/S1HNr/FgOW2gPwGYsE7MN+J3Y7y37ioHrct9tWh6O4reT/xge510PlHstzdTyKq10lbk/htusj++r2v++8P+dt8FqeL73/+QJWH+vqGF/ef/1S7balr5fk9ie8V+Y8HwA5zR6wKlydQaIZACaErLuPCYPc6SCqzYSpR/XdD4wmcKd+Qek14Vp+Snr/8wRsFyru9fJrTV+/rj+pdt/S1ytV63OfD4CcZg5YeYNHAk7NSPfioV3WDAz1YKBUB/HatLmz30pU9nrtQNIjGrBun9wZr7tuf5seuJzBfMI5pfkDq0Udx+7L0HGacz3qsHbPw9tOXUu5rG4/mCcGzbpN2Lfr1uxzj1+n2P67bYyxARt+UdPaY1d7cVUfz/lMYp+/vH7mGqXaHX3nM3T9+0XPJ9IfdX5j+zvpfADkNnPANoNPz+CqheuowccaFIa/1Q+IDTh2wKqwskLRfj/UlgzYsI/OQJkS7XdM6trpftrn4bb7/Wz25/QzPFe1j1HXKf1ZnR+wB+dYwWfiX3N1ba2+D7UHeq7XQO326Tsfl3e82c8HQE4zB6xkBvnGqAHAHdTdwSg9aLfUvicErG2oLRI6zvrRbSMB0Gfw2JbktUv0sycwktc2ea5Gen/nBmwnvL7RLzVWTaTa22Ut73xG1G6f4fPRx/Hvl/nPB0BOCwSsTQ8+auAwg4scoP1BwRvIFgtYeyALBsWhtuHgUgNds51vVL8HQ8uSvHafC9h4OFifneLuwyyPXcOhcx6sH9VX+5h9x5HCvi8esMnrby+vTTqfhvV5EbDAti0csA37ph8xSC0XsB11DDXYjWkbEbDBOU2QGORfzfLktZsjYLtgbddz9tEJr1P6sxpbP24NxET6Hrs+6to2fU+1B7zzSV5/e7krfT6GdV6znw+AnOYNWHnDR29ua+CPDgBuMLiD0bgBTFH79oNK7rs/vPr3bbclgqv3vMca6Ic6p7HXLtHP4Dj6fbC+P2gPnp+9z4HzaCwasLHrY59Pqj3gnU/y+veLnk/0ulr7m/18AOQ0b8A2A5J/g6sZXjsQNOtYg40afKxtFgtYf0Cz1x9qC/rcvPfPyT5vtb++YI9Qx/PPUw+2wXF7r91MARtcB3mu46/TlwVscH309fDPt7/d55+Pv31Yu33i5xM5B3l92/3NfT4Acpo5YDU1mKhBuREMQM3A0NPuD0Y6oN1lUc5gb8hBp1vm980e3IbazOCl2+r9qX8KYwWNf05eP+S+0wOffQwt3Gbo2unthwLTvZZ6X+55+tdBbuvud+g6pT6rsQGbFu+7f31S1y/1mYTnM3T9PyO1v3nPB0A+iwQsIurwv2Lwo34AFIOAzeT1sPdm1mWifgCUgoBFVtQPgFIQsMiK+gFQCgIWWVE/AEpBwCIr6gdAKQhYZEX9ACgFAYusqB8ApSBgkRX1A6AUnwpYAACQFnsxg8UiqB8ApSBgkRX1A6AUBCyyon4AlIKARVbUD4BSELDIivoBUAoCFllRPwBKQcAiK+oHQCkIWGRF/QAoBQGLrKgfrMHr6S26HJ/HNQ3NHLB/xPWPX+Jb4Fn8jq6P0pwVsMe92O12UVdHs96buK265dXBv+lT7XEPe/sYa3QUV7tK3J78n3OSx/3K6zTmvPXnP/Zz/xy3xnZ7v+bPrdG11bB/Teeqv6H9bKHWFgrY6yd3+e+bOmR/voiTtQxlOitgI+SgsasO4lW9bwaedkDTN2F346fa414PVWSQXBv7hh93889PX88v/SJyOoiqrYcYPwzmp2qyrRf/eOfW6BpreKlrOlTHW6i1TAH719OzN4v9EHc/rRmuE75yH4/i7ulFfDft8v27aU9tjzUbG7DuINVDzWj34sG8lwXv35ByHXMTpNqj5I1sHUNS+6n7J/nb2m2NbhBoBoyjvY7Xn6F9m8E12m4PRv7AlNrOG6i866oGZ7Ot319Hs6+DfX7etRtzDkPXZ3B7TdZOfLD3tm3rK73PmGiNRmpMXT+zz1QNntseJa/rhBoevB6yrT7+QdaItY6hrodXf0PHstsawf1in2tr7bWmZQlYOYPtljXhePPHae9C0jxm7gLZbdfbf7//aLc/3T8SshsxX8Dq4neKOzbQqBu4ufFS7RHhN//mxlaDgN8Hu01T59HuX7fbA4Fqb/s0tO/wfJ2B2wwYaoCwf24GCesc+o/ZkNfJvmb29fHfO8Lz0+Fs3o/ri7292566Bo3Y59yKX9f+PvVT6w3WqOQdL1WD57ZHTKvh1PVo2p0++Puw62/oWHabpo7Vnou9H5/edt21tlDAtjNLSxuI73Jm6v9O1g7mSEjbM+BgNizp0A1mzlid2QJWFrY3qERvADX46Js01d4ua4U3Wey4g5z1wwHFaR/ad7RN90/vzx6MrJ+jg6/djzP6FIjsy162SF/sa2DI/Yz8TJN96jdco/o4/szn3BpdvIaT1yOyv2CZdf0n1U/NWX/oc4x9Rtay5HlEtp+91jLNYNUM0zzmlQEZzDbtWelwwOp9JUIcqzVYP6qom0HJ4xZ17CZfYnAKb8LoPgJ6u67/9oDRf1MP7Vu1Ofvs6Otg3+jWz3L/wT7t6zfcJ7OuOZb7OfjsPhjWsc7sS/oauPuM99U+Xv0+2Sd7eU31J96H3mvjn4N/vAk1mmpvl7XCaxrdh5G8HrFr61+vrg4Gj9XSfeyupX2/DJ2X32b1I3ke4XUJPienTx23LmLXo5MlYJ3lcwQsj4M3a7h+OoOzg74BJXZTqXWbGzbVHghv4uEBoxso2pvQumlH3dQ9+04PVHZfrZ9j5zxhoOm2aT4TdX6Ra6/YfTBmHvQGr4Fh7zPRluyTvdyVfMrS0vtT5xU73pQaTbUHws9k8Domr4d1LtF2+b475vn3y5Zr7SsCdpFHxNiKWQI2evPUYgONvW6qPTB8E7rrNm3+vpz1E/tL7buvTbEHHOvn6OBr92PiOQ4OKJF92cvO7UvyGhi6j24/3La2/8k+9YvWaLSPiWsgtxlbo6n2QOKa+pLXI3Zt/ZqQ6zf1N3SsWL+d9a392Ou0bX4/rGXJ84hsbx97qN+OoVrLFLDurHPcHzn1BqzZ3p7Fqnb/L42xRnMEbH+bLvauTd9E3Y2favf5A4fk3pjqm67Zn39Tqpu8Pp4TfAM39dC+Td/tAUltGwlV52f/nJvr1+7Hb2/e9w006pyGBz17/eFj+e2p65O6BoZ9/j69j/6a8PvUL16Hen/BOfReA33O42s01e7T64+u4WD/8c/QOb/gGPb1HzhWtLbqYznbbrnWFgrY8Pej/oyzCUnT7jzyTQWs5G1PuG7G2IAdIm+U1IAS/33JmHaXO/g02oFAsgaImrqJnTb7Rk7d1LWBfft9d292+0b3b3pvO3vQaLe19qn+6UJ3bPec/MHV1hzX/qcPwbGG+jLi+gxeg4bcJjhuR32mcls/VMw+B7YdJ7U/t316jabaXVNreLj/us2vAfeaevV31v0S+XyVbdTazAELDNte/Qzd5Fij4S9gJaKGl5KqNQIWWW2xfqIzAKyTnC0NzChKRQ0vYEStEbDIapv1E38shrVhptaPGp7XuFojYJEV9QOgFAQssqJ+AJSCgEVW1A+AUhCwyIr6AVAKAhZZUT8ASkHAIivqB0ApPhWwAAAgLfbqDVhevHjx4sWL1+dfBCwvXrx48eK1wIuA5cWLFy9evGZ/CfH/AT+Qsl3t6cjNAAAAAElFTkSuQmCC"}},"execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_pixels_hu(slices):\n    image = np.stack([s.pixel_array for s in slices])\n    # Convert to int16 (from sometimes int16), \n    # should be possible as values should always be low enough (<32k)\n    image = image.astype(np.int16)\n\n    # Set outside-of-scan pixels to 0\n    # The intercept is usually -1024, so air is approximately 0\n    image[image == -2000] = 0\n    \n    # Convert to Hounsfield units (HU)\n    for slice_number in range(len(slices)):\n        \n        intercept = slices[slice_number].RescaleIntercept\n        slope = slices[slice_number].RescaleSlope\n        \n        if slope != 1:\n            image[slice_number] = slope * image[slice_number].astype(np.float64)\n            image[slice_number] = image[slice_number].astype(np.int16)\n            \n        image[slice_number] += np.int16(intercept)\n    \n    return np.array(image, dtype=np.int16)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#lets take alook at patient\npatient = load_scan('../input/osic-pulmonary-fibrosis-progression/train/ID00007637202177411956430')\nimgs = get_pixels_hu(patient)\nplt.hist(imgs.flatten(), bins=80, color='c')\nplt.xlabel(\"Hounsfield Units (HU)\")\nplt.ylabel(\"Frequency\")\nplt.show()\n\n\nplt.imshow(imgs[29], cmap=plt.cm.gray)\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Resampling\n\nA scan may have a pixel spacing of [2.5, 0.5, 0.5], which means that the distance between slices is 2.5 millimeters. For a different scan this may be [1.5, 0.725, 0.725], this can be problematic for automatic analysis (e.g. using ConvNets)!\n\nA common method of dealing with this is resampling the full dataset to a certain isotropic resolution. If we choose to resample everything to 1mm1mm1mm pixels we can use 3D convnets without worrying about learning zoom/slice thickness invariance.\n\nWhilst this may seem like a very simple step, it has quite some edge cases due to rounding. Also, it takes quite a while.\n\nBelow code worked well for us (and deals with the edge cases):","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Slice Thickness: %f\" % patient[0].SliceThickness)\nprint(\"Pixel Spacing (row, col): (%f, %f) \"% (patient[0].PixelSpacing[0], patient[0].PixelSpacing[1]))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Although we have each individual slices, it is not immediately clear how thick each slice is.\n\nSlice Thickness: 10.000000\nPixel Spacing (row, col): (0.652344, 0.652344) \nThis means we have 10.00 mm slices, and each voxel represents 0.652344 mm.\n\nBecause a CT slice is typically reconstructed at 512 x 512 voxels, each slice represents approximately 370 mm of data in length and width.\n\nUsing the metadata from the DICOM we can figure out the size of each voxel as the slice thickness. In order to display the CT in 3D isometric form (which we will do below), and also to compare between different scans, it would be useful to ensure that each slice is resampled in 1x1x1 mm pixels and slices.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def resample(image, scan, new_spacing=[1,1,1]):\n    # Determine current pixel spacing\n    spacing = map(float, ([scan[0].SliceThickness] + list(scan[0].PixelSpacing)))\n    spacing = np.array(list(spacing))\n    \n    resize_factor = spacing / new_spacing\n    new_real_shape = image.shape * resize_factor\n    new_shape = np.round(new_real_shape)\n    real_resize_factor = new_shape / image.shape\n    new_spacing = spacing / real_resize_factor\n    \n    image = scipy.ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest')\n    \n    return image, new_spacing","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pix_resampled, spacing = resample(imgs, patient, [1,1,1])\nprint(\"Shape before resampling\\t\", imgs.shape)\nprint(\"Shape after resampling\\t\", pix_resampled.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.hist(pix_resampled.flatten(), bins=80, color='c')\nplt.xlabel(\"Hounsfield Units (HU)\")\nplt.ylabel(\"Frequency\")\nplt.show()\n\n\nplt.imshow(imgs[29], cmap=plt.cm.gray)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"   \n# 3D Plotting","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_3d(image, threshold=-300):\n    \n    # Position the scan upright, \n    # so the head of the patient would be at the top facing the camera\n    \n    p = image.transpose(2,1,0) #get image in order(h,w,c)\n    \n    verts, faces = measure.marching_cubes_classic(p, threshold)\n\n    fig = plt.figure(figsize=(10, 10))\n    ax = fig.add_subplot(111, projection='3d')\n\n    # Fancy indexing: `verts[faces]` to generate a collection of triangles\n    mesh = Poly3DCollection(verts[faces], alpha=0.70)\n    face_color = [0.45, 0.45, 0.75]\n    mesh.set_facecolor(face_color)\n    ax.add_collection3d(mesh)\n\n    ax.set_xlim(0, p.shape[0])\n    ax.set_ylim(0, p.shape[1])\n    ax.set_zlim(0, p.shape[2])\n\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_3d(pix_resampled, 400)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Normalization\n\nOur values currently range from -1024 to around 2000. Anything above 400 is not interesting to us, as these are simply bones with different radiodensity. A commonly used set of thresholds in the LUNA16 competition to normalize between are -1000 and 400. Here's some code you can use:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"MIN_BOUND = -1000.0\nMAX_BOUND = 400.0\n    \ndef normalize(image):\n    image = (image - MIN_BOUND) / (MAX_BOUND - MIN_BOUND)\n    image[image>1] = 1.\n    image[image<0] = 0.\n    return image","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Zero centering\n\nAs a final preprocessing step, it is advisory to zero center your data so that your mean value is 0. To do this you simply subtract the mean pixel value from all pixels.\n\nTo determine this mean you simply average all images in the whole dataset. If that sounds like a lot of work, we found this to be around 0.25 in the LUNA16 competition.\n\nWarning: Do not zero center with the mean per image (like is done in some kernels on here). The CT scanners are calibrated to return accurate HU measurements. There is no such thing as an image with lower contrast or brightness like in normal pictures.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"PIXEL_MEAN = 0.25\n\ndef zero_center(image):\n    image = image - PIXEL_MEAN\n    return image","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"With these steps your images are ready for consumption by your CNN or other ML method ","execution_count":null}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}