{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"%env JOBLIB_TEMP_FOLDER=/tmp","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!conda install -c conda-forge gdcm -y","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Table of Contents\n* [Coordinate System For Medical Imaging](#coordinate-system)\n* [About data](#about-data)\n* [What is DICOM?](#dicom)\n* [Conversion to Hounsfield Units](#hounsfield)\n* [Windowing](#windowing)\n* [Histogram Analysis](#histogram-analysis)\n* [Storing metadata in dataframe](#metadata)\n* [Voxel Size and Volume](#voxel-size)\n* [Resampling](#resample)\n* [3D Plotting](#3d-plot)\n* [Segmentation](#segmentation)"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport scipy.ndimage\nfrom skimage import measure, morphology\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nfrom sklearn.cluster import KMeans\nimport shutil\nimport cv2\nimport pydicom\nimport matplotlib.pyplot as plt\nfrom matplotlib.ticker import (MultipleLocator, FormatStrFormatter,\n                               AutoMinorLocator)\nimport seaborn as sns\nfrom IPython.display import HTML\nimport gdcm\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\n'''\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n'''\n# You can write up to 5GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def set_options():\n    pd.set_option('display.max_columns', 100)\n    pd.set_option('display.max_colwidth', None)\n    pd.set_option('display.max_rows', 1000)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"set_options()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"coordinate-system\"></a>\n# Coordinate Systems in medical imaging\n\nCoordinate system is used for identifying the location of a point. Three types coordinate systems commonly used in imaging applications: the world, anatomical and the medical image coordinate system.\n\nWe will talk about Anatomical coordinate system. This system has three planes and **dataset we will be using is on Axial Plane.**\n\n* Axial plane -  The axial plane is actually when you place point of view above the patient and look down. Depending on the region of the 3D medical image you will observe different anatomical structures. For a 3D total body scan, if you had a control-bar over this 2D view you would start from a 2D slice of the head, and by increasing you would end up in the legs. Let’s practically call this view the “drone plane” or “top-view”. Slices near to head is known as superior and towards feet is known as inferior. Below is axial plane of lungs CT Scans.\n\n* Sagittal plane - Basically, this is a side view. Instead of looking from above the patient, now we look from the side. The side can be either right or left. Which side and direction is the positive one, depends on the coordinate system.\n\n* Coronal plane – In this point of view is either in front of eyes(anterior plane) or back of the patient(posterior plane)"},{"metadata":{},"cell_type":"markdown","source":"<a id=\"about-data\"></a>\n# About data\nMedical dataset containing lungs CT scans of patients diagnosed with pulmonary fibrosis a disorder with no known cause and no known cure, created by scarring of the lungs. Prognosis of the troubling disease becomes frightening for the patients because outcomes can range from long-term stability to rapid deterioration, but doctors aren’t easily able to tell where an individual may fall on that spectrum. This is where data science can help in predicting the detoriating condition of the patients. Detailed description about data is found [here.](https://www.kaggle.com/c/osic-pulmonary-fibrosis-progression/data)"},{"metadata":{"trusted":true},"cell_type":"code","source":"TRAIN_PATH = '../input/osic-pulmonary-fibrosis-progression/train.csv'\nTRAIN_IMG_PATH = '../input/osic-pulmonary-fibrosis-progression/train'\nTEST_PATH = '../input/osic-pulmonary-fibrosis-progression/test.csv'\nTEST_IMG_PATH = '../input/osic-pulmonary-fibrosis-progression/test'\nSUBMISSION_PATH = '../input/osic-pulmonary-fibrosis-progression/sample_submission.csv'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"'''Load data'''\ntrain = pd.read_csv(TRAIN_PATH)\ntest = pd.read_csv(TEST_PATH)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Duplicates on basis of patient and weeks\ntrain[train.duplicated(['Patient','Weeks'], keep=False)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Remove duplicates\ntrain = train.drop_duplicates()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"dicom\"></a>\n# What is DICOM?\n\nA DICOM image file is an outcome of the Digital Imaging and Communications in Medicine standard and represented as .dcm. Because of its ease of integration and continuous evolution this communication standard has over the years achieved a nearly universal level of acceptance among vendors of radiological equipment.DICOM differs from other image formats because it groups information into datasets. DICOM file consist of header and image data collectively in one file. We will see later how these group of information looks like and interpreted.\nDetailed information about DICOM can be read [here.](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3354356/)"},{"metadata":{},"cell_type":"markdown","source":"Video for better understanding about DICOM images."},{"metadata":{"trusted":true},"cell_type":"code","source":"HTML('<iframe width=\"600\" height=\"400\" src=\"https://www.youtube.com/embed/KZld-5W99cI\" frameborder=\"0\" allow=\"accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## In this dataset ImagePosition is not available so we will sort the slices on InstanceNumber.\ndef load_slices(path):\n    filenames = os.listdir(path)\n    slices = [pydicom.dcmread(f'{path}/{file}') for file in filenames]\n    slices.sort(key = lambda x: int(x.InstanceNumber), reverse=True)\n    return slices","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scans = load_slices(f'{TRAIN_IMG_PATH}/ID00007637202177411956430')\nscans[0]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"hounsfield\"></a>\n# Conversion to Hounsfield Units\nThe Hounsfield scale is a quantitative scale for describing radiodensity in medical CT scan and provides an accurate density for the type of tissue. Plain x-rays which only displays 5 densities (i.e. air/fat/soft tissue/bone/metal), CT displays a huge range of densities ranging from air (black) to bone (white). On the Hounsfield scale, water represented by 0 value, air is represented by a value of −1000 (black on the grey scale) and bone between +700 (cancellous bone) to +3000 (dense bone) (white on the grey scale). As bones are much denser than surrounding soft tissues, they show up very clearly in CT images. Raw pixel values of images gets convereted into Hounsfield Units because the spectral composition of the x-rays depends on the measurement settings like acquisition parameters and tube voltage. By normalizing to values of water and air (water has HU 0 and air -1000) the images of different measurements are becoming comparable.\nRead more about Hounsfield units on [Hounsfield Scale](https://www.sciencedirect.com/topics/medicine-and-dentistry/hounsfield-scale)\n\n![image.png](attachment:image.png)","attachments":{"image.png":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAewAAAFRCAYAAAC2fp7LAAAgAElEQVR4Ae2dvWsb2/fu958xrcHFz5Ai7qwyhhQRN8URuPgaUgRximBO8UW4CCKNGVIYkSKIUwSR4oBcHK5SBJTigFzccOUioBQHnCKgXEjhIoWKFCpSPJc9r3tmljSSRqMZ2U/A8Wi0X9b+rJn1zH6ZbQX+IwESIAESIAESKD0BVXoLaSAJkAAJkAAJkAAo2LwISIAESIAESGALCFCwt8BJNJEESIAESIAEKNi8BkiABEiABEhgCwgkBFspBf6QAa8BXgO8BngN8Boo9hqIP0OIgh1PxM/FEdA3DP+RAAmQAAncLQJS7E+ogZTobmEqV2vpj3L5g9aQAAmQwCYISLGfgr0J8hnqkJyWoThmJQESIAES2AICUuynYJfccZLTSm4yzSMBEiABEshIQIr9FOyMUPPOLjkt7zpZPgmQAAmQQLEEpNhPwS7WJ6m1S05LzcQEJEACJEACW01Aiv0U7JK7VHJayU2meSRAAiRAAhkJSLGfgp0Rat7ZJaflXSfLJwESIAESKJaAFPsp2MX6JLV2yWmpmZiABEiABEhgqwlIsZ+CXXKXSk4ruck0jwRIgARIICMBKfavINhD2ML2pfbVqtaF5dUublYtRMh3g+6RwnrLFKrJ+ZTkNLdKt336e3U2FKwQuF7Z3razNXS/J7MMz7xt+I66WKcnkjXxDAmQAAmQwDwCUuxfTrCDgC/vr7qaOArCMq8VC32XR5kLVbz2RJLT3Eoo2GuHzQJJgARIoCQEpNi/hGAbAiH0sJ2enpJ7bvPbn4e45lHm/Fbk9a3kNLcuwx/sYeeFn+WSAAmQQCEEpNi/hGCHIpgYgjV63n4ve9bw6s1FLTYsG5ar84bf61689AAQpncfEnQ6G8GgsGFL+L1CMGT/vYta/IEjNgQcsT2eXhLHeBpl1Ge4OijXq99nZSRJHEpOcxNRsBOweIIESIAEbgkBKfavJtiikEYpBeIUE8NQkH0xlgQ4OuQeiC3mpfVEe55gC8IaiLphZ2B7XNj9z6Zoz6hPlxsKsiGufhn+b7OsKELnk+Q0N5lRplhGyCqwI7DVZx+tMGi3wSKagp9IgARIgAQ2QUCK/UsINmK931BUA0EwWjEr+KcJtizOcTE2BMcQ4dAOQawABDYZApe0x0gXEd2wzLBHb5wLRM4QUq/nL9WBFPH0UUpOc78z6wl9odPHfwIuKXUGfIK2+FbwNwmQAAmQwCYJSLF/KcHWxgZBXRAGc6g8SBcL/knxMkTPEFIHTFxggs9alIxh8ATFsMxArOakcQUufAgIbI/VEdo+5wEC3oNN0G5DWCPtC8/PthGO+CZMd06E+eMCHf8clB/wC9tqlh20O7Dd/JbHJEACJEACmyKwFsH2jQ3FK9qj88VhVvAP8/miMUdcjd6z2/OWRcqv07cNxtB59LuwrriomfPl6bbPF+zQDn00r06PXUTIo7klp7kpZBbJdhlD8xTsKFx+IgESIIGSEpBi/9I9bLlthih5vbN00VtFsP3ajfqMnr40nB4KdkzgPJFMPkAYowixnmaYtmSCLQp+yChgEAh2yqK4WLt96vxNAiRAAiSwGQKZBDsUKynYm2Loitkqgm0OqTtIAoHxxV0CZdQdCI0gVkZPNxR2c14+rCPd9rhgx5g4dvvlGfaJwiq1KTwnOc39Nq1cgYExYpGYUghYK6iAY2gHj0iABEiABDZHQIr9i/ewI8E+OgyuCw5+Ej1XQ8wiZfiCFgqLLiMUU+O8JyDhQ4Of14UXCGww5xzmDXqXhmCHDwZhupWGxM0yDZEL7XHtFO02WIRtTl4MktPcVCsItmmv6bPYccgsaQ/PkAAJkAAJ5E9Aiv2LC7a/mCoW3AOhds4bQmr22MQ8flpTNA3hN/KEgjY/bSg0hph55egyQiGV6vHtMdIZIqzdEwqv18PWJ+e1M+hRz7E7Vkf8MpCc5qYx2hjUY+YO6wy5pNjrsDLaZhbHYxIgARIggY0RkGL/UoLtWhoKQUSsJdGIi9nZ0BA9XyDD8rSwhKKoRdVPE2UkCW8o6l7aWN2uaBkip8XJEcto/Tp3UH5MTEPbYqJm9JZ9JhGR9EwKyvUfRiRm0aYutkpcLCfZrrDo8DvfXue3WE6Yi0ckQAIkQAKbIaBjcvxf4oyUKJ6JnzdHgP7YHGvWRAIkQAJlISDFfgp2Wbwzww7JaTOS8jQJkAAJkMAtISDFfgp2yZ0rOa3kJtM8EiABEiCBjASk2E/Bzgg17+yS0/Kuk+WTAAmQAAkUS0CK/RTsYn2SWrvktNRMTEACJEACJLDVBKTYT8EuuUslp5XcZJpHAiRAAiSQkYAU+ynYGaHmnV1yWt51snwSIAESIIFiCUixn4JdrE9Sa5eclpqJCUiABEiABLaagBT7Kdgld6nktJKbTPNIgARIgAQyEpBiPwU7I9S8s0tOy7tOlk8CJEACJFAsASn2U7CL9Ulq7ZLTUjMxAQmQAAmQwFYTkGI/BbvkLpWcVnKTaR4JkAAJkEBGAlLsFwVbJ+QPGfAa4DXAa4DXAK+B4q6BuOaLgh1PxM/FEdA3C/+RAAmQAAncLQJS7E+ogZTobmEqV2vpj3L5g9aQAAmQwCYISLGfgr0J8hnqkJyWoThmJQESIAES2AICUuynYJfccZLTSm4yzSMBEiABEshIQIr9FOyMUPPOLjkt7zpZPgmQAAmQQLEEpNhPwS7WJ6m1S05LzcQEJEACJEACW01Aiv0U7JK7VHJayU2meSRAAiRAAhkJSLGfgp0Rat7ZJaflXSfLJwESIAESKJaAFPsp2MX6JLV2yWmpmZiABEiABEhgqwlIsZ+CXXKXSk4ruck0jwRIgARIICMBKfavLNiTd3V3+9KnPUwEw24ualCqhu534UueWpiA5LSFM4sJx+g+sVC7uBG+neL6ooHqPcvx7c5BHfaHsZBuguF5HZVdvWWfhb3HJ+h8Eq6CH0O0nlawo7e6tfZQfdbBSEgmVMBTJEACJHCnCUixf0XBnqD3VGH/oAJLVdH5muRKwU4yWeWM5LRVynHy/Bqj90fFEWNJsEev9Hc7qJ110b/sof37ISxl4fhvU9wn6D+zoKxD1F/3MPjQhX20A6UqsD9NQ9N+9HFiKVgP6mi/G6B/YaOmBf7AxshIFmbgEQmQAAmQgE9Aiv2rCfbXDqpKofGuj6alsH8+8uvg7zUTkJy2ShWTz100Hrg9Z11mQrB/9FBXCtU3Zo96iuGLfSiriYEvsp9b2FcWGv/4J7Q1Y3SPFNSjDvzco3Odr4HBT8Pab13UdB1v/VTGdzwkARIgARIICEixfyXBHr85hFIn6E+mGDzXvS0joHvVxXvY7mcb/asWqrqntVtBMxL0Azt5YBCQnGZ8vdjhd1co1W4N7c992IJgT96fyFMYnkDbV25Vo1daiG0MYzVH84/Quq9gnSVSof9MQR11YfbZY0XxIwmQAAnceQJS7F9BsN1grJ713bnrKxuWUqi/i05OyoJtwdqtwnaGSLsY/rjzPkkFIDktNVM8wfceGqddXDsuGoqCPXqpe99JIcZ0gIYW+L90r/gGvSczBPdrB4dKoflRJ+vhWHgo0Ga5D3vNhODHTeZnEiABErjLBKTYv7xgewIdDokOYVsKKrb4TBZshZMP5lDqXXbHYm2XnLZYzlmpZMEenimo37rBkHaY202vnN7yjTv0/SLecwbg9eKdoXbv2BHvsCDnKH5dxL7mRxIgARIgAcBZaxQHsaRgy0PgznxlbPFZPDDHP8cN4WeZwEzB9oe59QrsGT/+MHa0ZEmwpxiczug5wxTsMTqPFFzxjpYaEWxvjYNUP6+DGDd+JAESIAGBgBT7lxNsb2GS9XyASD/ZmeeMLj6LB+b4Z8E+nhIISE5zkk2G6Jy30Jrz0xdW78MT4PiiM/awBfg8RQIkQAIFEZBi/1KCHbx7PaNHZy4+iwt0/HNBDLauWslp2Roh9bABzmFno8rcJEACJLBOAlLsX0KwveFQ6wSdywEGsZ/eWdUZmvUXn8UFOv55nQ27zWVJTsvWXlmwo6u8jRq80RN/Pjp1lfg3nTdllbg4V27UyUMSIAESuOMEpNi/uGB7q4BnvnPtDZf77+LGBTr++Y77YuHmS05bOLOYUBZs8D1skRZPkgAJkEARBKTYv7BgSwvLoo3wFqSpfbQ+A3GBjn+O5uWnWQQkp81Ku9j5GYKNKUYv9U5nFqrPYzudXZgbnehXu/QrYJXkTmdXxsqG7z0c67cHDpI7nQ3NzVQWM5qpSIAESOBOEZBi/2KCPR04O5rFX91K0DMWn8UFOv45kZcnRAKS08SEC5+cJdi6gMlie4n/usFgkb3EbwbcS3xhvzAhCZAACYQEpNi/mGCHZfBowwQkp23YBFZHAiRAAiSwYQJS7Kdgb9gJy1YnOW3ZMpieBEiABEhguwhIsZ+CXXIfSk4ruck0jwRIgARIICMBKfZTsDNCzTu75LS862T5JEACJEACxRKQYj8Fu1ifpNYuOS01ExOQAAmQAAlsNQEp9lOwS+5SyWklN5nmkQAJkAAJZCQgxX4KdkaoeWeXnJZ3nSyfBEiABEigWAJS7KdgF+uT1Nolp6VmYgISIAESIIGtJiDFfgp2yV0qOa3kJtM8EiABEiCBjASk2E/Bzgg17+yS0/Kuk+WTAAmQAAkUS0CK/RTsYn2SWrvktNRMTEACJEACJLDVBKTYT8EuuUslp5XcZJpHAiRAAiSQkYAU+0XB1gn5Qwa8BngN8BrgNcBroLhrIK75omDHE/FzcQT0zcJ/JEACJEACd4uAFPsTaiAluluYytVa+qNc/qA1JEACJLAJAlLsp2BvgnyGOiSnZSiOWUmABEiABLaAgBT7Kdgld5zktJKbTPNIgARIgAQyEpBiPwU7I9S8s0tOy7tOlk8CJEACJFAsASn2U7CL9Ulq7ZLTUjMxAQmQAAmQwFYTkGI/BbvkLpWcVnKTaR4JkAAJkEBGAlLsp2BnhJp3dslpedfJ8kmABEiABIolIMV+CnaxPkmtXXJaaiYmIAESIAES2GoCUuynYJfcpZLTSm4yzSMBEiABEshIQIr9Swj2DbpHCuqoi5uMhjD74gQkpy2e2085gm3J2+vVLkxvTnF90UD1nuVsTbtzUIf9YewXYvyeYHheR2VXl2lh7/EJOp8mxvc8JAESIAESyEJAiv0U7CxEN5BXctrS1X7voqYUqs9aaJ1HfzpXodCOXlWg1A5qZ130L3to/34IS1k4/tsU9Qn6zywo6xD11z0MPnRhH+1AqQrsT9OlTWMGEiABEiCBJAEp9lOwk5xKdUZy2tIGfmxCqRq63+bk/NFDXYv6G7NHPcXwxT6U1cTA1+LPLewrC41//BO6zLE7+vKoAzP3nNr4FQmQAAmQwBwCUuxfu2DfXNSglI1hzJDoeW94/WyIyccW6ge6h6Zg3auicXENUwqAKa7f2Tj2h2kfNdD9cu0KxFm8llilt+Cj5LRlmzX+S/aJWc7k/Ykr6t/NswA8gbav3POjV1rAk/6dmT9WHD+SAAmQAAmkE5Bif7GCfVBBZbcG+6JvDK1aaF6Gkj2+OHaGZQ9/b6N32Uf3rIYdq4LKfQVFwU73OoDhCwX1pIOBnp925p13UHnawsAY6R691PPWSSHGdICGUqj9pfvON+g9mbGO4WsHh0qh+XEhk5iIBEiABEhgDoHyCXZimHboLo7yhdgbpq2cDSO9blfEKdhzfG18NUbnobs4zH3oGaB/YaOmhds6DobJh2cK6reuMKQ9hK3/PrrjE29k5IUwsuHNk0cXsRlm8JAESIAESGBhAuUT7IRARFeiu8Osh+h8jbdxiGYgIvHvbtdnyWlOCz2B1N/P+nGGsX8N0TrYQfyhB5M+Gnrl+NMeJphicDqj5wxTsMfoPPLFO8aZgh0Dwo8kQAIksDoBKfYXOySeeEUsKtju3GszMR+uh2adV8z8nvjqTEqfU3KaY/RkiE5sxXd8BXg/8aATbe7wTA+DN5wFZexhR9nwEwmQAAkUSUCK/eUW7DeHUIqCnddF4z4Q1dD9DnAOOy/KLJcESIAEliewQcE+QT98vdexNCoI0Z502JTo+XlD4s5GIHe5hx1Cm390ZTtD5tHXsHQWbxjcsjECMHOVt7NKPFxMlrpKfN6rY/Mt5bckQAIkQAIegY0Ithv499H6bHD/OURTr+oOViFHhTlMGTs/Y9HZ9J8GrLs+hx1Cm380HaBpKVhPYgvKvrmbqeyfa7kGwPew53PktyRAAiSwQQLrEez7x2iKc6cdDHWv+kcfJ3ox00Ed7XcDDN61UX+wg+Mn5rvAMWEOICTPj9/WjNe6Buid69e6LAp2wCz9wF9Vv3PUQu9SrxJvour4yMbwp59/itFLvdOZherz2E5nF+Z2KPrVLj33XUnudHYVvo7nl8rfJEACJEACyxNYj2DPXJXszoVqs6Zfumg83nNENdgMxRma9d/zTQqz2xzpfGx/a2fjlL7xqtHyILYph+S0Vey/uQw3qFG7FdTPergOxNovcbLYXuK/bjDgXuI+NP4mARIggbUTkGL/EovO1m5PhgLd17rczTwyFLMFWSWnbYHZNJEESIAESCADASn2l1uwP7dQuXeI1qdoq6cfm85+1v52mdFvb9cnyWm3q4VsDQmQAAmQQJyAFPvLLdgYoXWgoPztSy8H6L2u41BaRBVv7S35LDntljSNzSABEiABEphBQIr9JRds/b7RNbqnNe9vL7t/IOTk9QA3v2a08padlpx2y5rI5pAACZAACcQISLG//IIda8Rd+yg57a4xYHtJgARI4K4RkGI/BbvkV4HktJKbTPNIgARIgAQyEpBiPwU7I9S8s0tOy7tOlk8CJEACJFAsASn2U7CL9Ulq7ZLTUjMxAQmQAAmQwFYTkGI/BbvkLpWcVnKTaR4JkAAJkEBGAlLsp2BnhJp3dslpedfJ8kmABEiABIolIMV+CnaxPkmtXXJaaiYmIAESIAES2GoCUuynYJfcpZLTSm4yzSMBEiABEshIQIr9omDrhPwhA14DvAZ4DfAa4DVQ3DUQ13xRsOOJ+Lk4Avpm4T8SIAESIIG7RUCK/Qk1kBLdLUzlai39US5/0BoSIAES2AQBKfZTsDdBPkMdktMyFMesJEACJEACW0BAiv0U7JI7TnJayU2meSRAAiRAAhkJSLGfgp0Rat7ZJaflXSfLJwESIAESKJaAFPsp2MX6JLV2yWmpmZiABEiABEhgqwlIsZ+CXXKXSk4ruck0jwRIgARIICMBKfZTsDNCzTu75LS862T5JEACJEACxRKQYj8Fu1ifpNYuOS01ExOQAAmQAAlsNQEp9lOwS+5SyWklN5nmkQAJkAAJZCQgxf6lBHt4lrZFm43hMkZOrtH9bxujZfLcsbSS07IhGKP7xELt4kYoZorriwaq9yxna9qdgzrsD2Mh3QTD8zoqu/p6sLD3+ASdTxMhXexUUf7+ZMMSt9utofvdtHHFdplF8JgESIAE1kBAiv0rCHYVJ+cttMSfPqTwPsv2m4salFpS5GcVdkvPS05buam/xuj9UXHEWBLs0Sv93Q5qZ130L3to/34IS1k4/tsU9wn6zywo6xD11z0MPnRhH+1AqQrsT9O5phXlb7de6brtYBg8Z6zerrmN5pckQAIksAIBKfavINjrE9iiAvgK7ArLIjltFWMmn7toPHB7zrrMhGD/6KGuFKpvzEeuKYYv9qGsJga+Fn9uYV9ZaPzjn9DWjNE9UlCPOnMf2Iry9/CFgvqtO9c2ZGjXKv5gHhIgARKYR0CK/bkJ9vRLD/bTQ+xZ3jD6bgW10y6uf7omJobXz5YaTJ/Xzlv1neS0pRv4vYuaHhLeraH9uQ9bEOzJ+xMoFR8iBnwhs6/cWkevtIAnH9pm5veMnefvm8s2Th7vucPW1h6qz9oY/li6lTMyjNH9TUGlXF+rtmtGpTxNAiRAApkISLE/H8H+2kFVKewc2eh+GGBw2UPnv1UnIFvP+tCjkJMvA/ReHEKpOtqXAwy+BGOTmRp52zJLTlu6jd97aOiHJQfxUBTs0Uvd+04KMaYDNLTA/6V73jfoPVFQR12Yg+SOPV87OFQKzY+ydbP8Pb44dq6LytM2epcD9C9s1PTcuHWM7je5rOXODtFUCsdvBuieVrHjPLhUUD8f4OaXX9Lq7fJL4G8SIAESWCcBKfbnItijPw+xd2BjaI6aYorB86goFDVEuk6oeZclOS1bnbJgOz1gcdjYTe/2UG/coe8XwmiI14tPDLUbxib87Q3D778YInKpTPpoWArW80H0vFHWwofeg0Qw537ZR/es5gi39cQfJs/WroVtYUISIAESWJCAFPtXEGxviDux6lboncUMcwN2I5gPTQTwWHp+hLNATOTgD3Mn/BD6xx/GjuaXBHuKwemMnjNMwR6j82jG8PIKgu0Oox+i8zVqof40Ot+HUuG1Ek+RGGKPc/BHAa5aqOxWYF9FHgkw+dBwevb1d3rYIVu74rbxMwmQAAlkJbAmwZZW2+pV48kV4tPJDa6v9DBnG81nVW8+O5wnpWCnu1RympNrMkRHXKkfruDvC0IIT4DjPeEietjxBziThvtdeK2Y3+nj8fuwneIbC2+HztRLPF/4eQhbr6841b149rBDLjwiARIoAwEp9q/Qw07vSeNb31iRvIPKw0PUT9toneo56zAIU7DTLwvJaem55qWQethA3nPY2qK4v7MI9rwWJr4L5qrNb7zFaE5PnHPYJhkekwAJFE9Aiv05CPYEvacK6n4D/e/RYcjxGwr2speB5LRly4imlwV75ipv53WncDFZ6mrqOQvF4oKdPiTeXG4jnmhDnU/u0LkwtO4tprNeutv2ZGmXUC1PkQAJkEAmAlLsz0GwvTnP+MKkX9do6/lP9rCXcqLktKUKSCSWBRtFvIedtujsD/eNgkQTljgxvWy6m79cmO+XA+O/9KY9+2h99grje9hLUGVSEiCBvAlIsT8HwZ64O2HpHbPOexh4r+oc37Ows6t3xAqHxN0e1j5O3vT5WtcM70tOm5F0wdMzBBtTjF7qnc4sVJ/HdjqLiJ0ePtar/SvJnc5iC7viBkn+zv+1LncrVmcHN+d67KP73H3FsHJmrk5fvV3xdvIzCZAACWQlIMX+HAQbwM9r551Xd9MUC3sP6mh9GGPi9WKal95Q+c8R2v/xNsx40ku+25u1xbcgv+S0bM2aJdi61Mlie4n/usFglb3EZ/g7sXHKaQejtW2cAiBir4KzR/q76+QrY5F0S+yRns0hzE0CJEACCQJS7F9KsBMl8kTuBCSn5V4pKyABEiABEiiUgBT7KdiFuiS9cslp6bmYggRIgARIYJsJSLGfgl1yj0pOK7nJNI8ESIAESCAjASn2U7AzQs07u+S0vOtk+SRAAiRAAsUSkGI/BbtYn6TWLjktNRMTkAAJkAAJbDUBKfZTsEvuUslpJTeZ5pEACZAACWQkIMV+CnZGqHlnl5yWd50snwRIgARIoFgCUuynYBfrk9TaJaelZmICEiABEiCBrSYgxX4KdsldKjmt5CbTPBIgARIggYwEpNhPwc4INe/sktPyrpPlkwAJkAAJFEtAiv0U7GJ9klq75LTUTExAAiRAAiSw1QSk2E/BLrlLJaeV3GSaRwIkQAIkkJGAFPtFwdYJ+UMGvAZ4DfAa4DXAa6C4ayCu+aJgxxPxc3EE9M3CfyRAAiRAAneLgBT7E2ogJbpbmMrVWvqjXP6gNSRAAiSwCQJS7Kdgb4J8hjokp2UojllJgARIgAS2gIAU+ynYJXec5LSSm0zzSIAESIAEMhKQYj8FOyPUvLNLTsu7TpZPAiRAAiRQLAEp9lOwi/VJau2S01IzMQEJkAAJkMBWE5BiPwW75C6VnFZyk2keCZAACZBARgJS7KdgZ4Sad3bJaXnXyfJJgARIgASKJSDFfgp2sT5JrV1yWmomJiABEiABEthqAlLsp2CX3KWS00puMs0jARIgARLISECK/QsK9g26RwrqYQfjhBHed0rh8E3yW3zt4FAp1C5uEjnlExNcXzTQ/iR/e9fOSk4zGUy/dNF4vAdLbye7W0H9rI/xLzMFj0mABEiABLaNgBT7FxRsYPzmEEodoxfX3UkfJ8qCZSmooy4SX78/gVKH6HxdENf3LmpKwb5aMP0tTyY5LWjy5xYqSmHnyEb3wwC913UcWgrWk17CD0EeHpAACZAACZSegBT7FxZsfG5hXyk0/plGGjr9pwGlTtB5q4X5BP1J5GsMXyio+y2Moqdnf6JgR9hITnMTTNB7qqAeRUc9ph+b2FcWmpdRP0UK5QcSIAESIIFSE5Bi/+KCjSGaSmH/VVR6h2eW27P2hDYq6CO07itYZ0MXzI8h2s9qqOx6f/3E2sPh0xaGPzxuV3bsr4TZ8HICNwO0nlaw4/wlsR1UntrofTFFaQhbD73/NUTvdz1EbGHvqIPrUrsk3TjJaU4uZ2RDmmqIMU+vgilIgARIgARKRkCK/UsI9gT9Z/Fhb1cc3LnrMToPDXHWjTdF/OcADT1c+6CO9rsBBpd9dM/rzpCuOvB64D+uMfi76cx51//Uaa6hO+zTTy1nqNcf+h186MI+2oFSFdhXvmi7gm1ZFip/dNC/7KHzftvlGs4DjHgdfbKdeevk1MEUg1MF9VtXWG8glsSTJEACJEACJSOQUbCBm4salGpg4Guks6BsH63PbktHr/Yjw98TZ/7aHSafvG+gsltD93uUyvhtFUoZ5xND4mN0HilYp31HvMPcWpgsY0jYFWz1sL31veqwjXME2xmNqKH7zUztHg/P9AiGMTqRTMIzJEACJEACJSaQWbDdFd+hQDsCbhnCcKV7feECs2C4fB4UR3jCPH6vPOg5fnMXodXO+xhc6l53+NN7oRfC+WLvCbY//D6vzi36TnKaNn96qdcO+G2PNoiCHeXBTyRAAiSwbQSk2L/EkLhu7gi2peeJ9UAoSIUAACAASURBVOtb3hD56QB+hxvTARrBK1zuEHniVa+fE4w/DzB410HrtB7MZwcCHe9hJ+a1vflvZy5bH1uwnVfA7pZggz3sbbv/aC8JkAAJLExgDYLtzY9qkfbE+eS9uSzcEPGbHo4DMQXwc4S2M+/siuzeg0PUnjXROas787Rpgh18P7O5d0ywOYc980rgFyRAAiSw7QTWINiAMy9t2Rg5r3klh2SdYfL7LQyd173C+e7R+b4zhNv+bAo84L4WZrx3He9hexuvJHrqCW/cMcFOWSWuXgTr6xOkeIIESIAESKDcBNYi2O4ccw3N5zV55zNHYI9x8sc+1DN/oZi3G1pi5fIEfb1wzNwoJS7YcBedqftNDH+agKcYvtiHda+JgXP+jgk2+B62eTXwmARIgARuE4H1CDbcV7l0YfF3sl1Y+nt35zNzO9LRq4oz33z4X/3KlZ7DbuPk0Q6s3Z3o60lez3H/mU7nvdZ1Zbuvf+3W0NKvhH3oovX7ofOu9fGFvx3qXRNs/bqby8V63IztdMZXum7Tjcu2kAAJ3D0CaxJswFn9reenZ2wf6uxuZqwWd1FPMNTvXXubpuwc1NB4O8TNjd7aVKH61hfeKUavj7Gntzo1tkKdfunBDjZOsbD3oI7WpbkR6t0TbM118i/3Er97tzJbTAIkcNsJrE2wbzuoMrVPclqZ7KMtJEACJEAC6ycgxf4lX+tav1EscT4ByWnzc/BbEiABEiCBbScgxX4Kdsm9Kjmt5CbTPBIgARIggYwEpNhPwc4INe/sktPyrpPlkwAJkAAJFEtAiv0U7GJ9klq75LTUTExAAiRAAiSw1QSk2E/BLrlLJaeV3GSaRwIkQAIkkJGAFPsp2Bmh5p1dclredbJ8EiABEiCBYglIsZ+CXaxPUmuXnJaaiQlIgARIgAS2moAU+ynYJXep5LSSm0zzSIAESIAEMhKQYj8FOyPUvLNLTsu7TpZPAiRAAiRQLAEp9lOwi/VJau2S01IzMQEJkAAJkMBWE5BiPwW75C6VnFZyk2keCZAACZBARgJS7BcFWyfkDxnwGuA1wGuA1wCvgeKugbjmi4IdT8TPxRHQNwv/kQAJkAAJ3C0CUuxPqIGU6G5hKldr6Y9y+YPWkAAJkMAmCEixn4K9CfIZ6pCclqE4ZiUBEiABEtgCAlLsp2CX3HGS00puMs0jARIgARLISECK/RTsjFDzzi45Le86WT4JkAAJkECxBKTYT8Eu1ieptUtOS83EBCRAAiRAAltNQIr9FOySu1RyWslNpnkkQAIkQAIZCUixn4KdEWre2SWn5V0nyycBEiABEiiWgBT7KdjF+iS1dslpqZmYgARIgARIYKsJSLGfgl1yl0pOK7nJNI8ESIAESCAjASn2LyHYN+geyVu07RzUcPJ6iIlh4M1FDUrZGBrn8j8cwlYKtYub/KvaUA2S01apevKxjZNHO+6Ws7sV1M8HuPkVL2mK64sGqvcsJ93OQR32h3E8kfB54uRrfxK+yuvUlZ2yfa557a3arryMZ7kkQAIkMJ+AFPuXF+z7x2iet9AKfpo4ebwHSylYz/qBaFOw5ztj0W8lpy2a1083vbJRUQqVp230Lgfova7j0FKwnnRhyvHoVQVK7aB21kX/sof274ewlIXjv1MegL53UVMK9pVf4wZ+f+0b16BxPZ7WsKevxSc9+Fav3K4NNINVkAAJkIBEQIr9ywv2UTcIhGElUwxf7EOpY/S8KEnBDulkOZKctlx5I7TuK6iY36b/NBwxtv1e8Y8e6kqh+saUcM+vVhOD6ZxaixBs0RzP3gMbw59egiztEuvgSRIgARLIn4AU+9ck2IAr0DV0v7kNEQV7co3eWR2VXXdo3RlyfXeNuBZM/u3BflrBjvNXw3ZQeWqj9yWeyhzmtLD3nzZGP+/ukPjwLCnKjidu+mg82MPJe3PCAsDnFvaNXvHk/QmUqqH7PXYhOums2b3nxNB0OBS9mB9j9cU+zmxXLJ3+6D6E7KP5MbxWVm6XUD5PkQAJkMCmCOQm2NPvA9iPosOQCcHWYnqgoKwqmhd9DIIhV4XK2TAQbX/41nrcRPfDAIN3bdQf6DnVCuyrMBC7w5wWDv/bQf+yj+5ZDXtPjp2h2bs4h72UsH0bwH5sQR3YGHlIRy8141Bsg4tyOkBDrwv4y+x5B98CP64x+LuJQ6VQ/3OAweW1My2yqB+NksTDhds1HcK+H52W0QWu3C7RGp4kARIggc0QWI9gz/pb2f/pYBzqqdfjDgVgdK6HzMMeuN/k6z+rUKqKzld9Jhy+jcjDr2u0HymoRx13ztUb5jSFXuce/6UXut3NRWcLCZs3dK0Zqf9poO9P8gJw8v8WndN2feSOWqizOcsHE0PiC/rRvwjm/F6oXdr3b83rKCwwU7vCYnhEAiRAAhslsB7Bnrfo7HELI2/uMNrD9gL46SDoSQct98T3UM+desO0jX8M5fcSTt7VodShI+x66NM/DsrRB35v8C6sEjfFd8ZDVGIR2LcR+pcD9C9aqDujHcfeFMYUg9MZQ+pYQbAX9GPEd/6HVdrlPeiZix7d4jK2y7eJv0mABEhgwwTWI9ixxUt+GyaXTWdOtPL62jkVFWw36Dui7GcIfhuC4MyHuqIcfO0feHOlWoTcshvCQij3weBODIlPhugEK/VbONEjELGHqb4zauEDjP3+OUDDUlBPe84QdqaeaLyHvaAfYxa5H1do1/Sy6Syga14mH/QytUs0kCdJgARIIH8CuQq2OZytR1qLEewxOg85JL7opTR8oRf/udMWmeZ61ynYMePTh8SnGDy3oKxw+sUsIlO7zIJ4TAIkQAIbJFCgYKcPiTuLmlKHUt05cA6JJ6+amcL2qY3K7j5an+N5vOFiy8YIwPzV1ArNj/H8xue4YC/oR6OEmYcz2xXkGMHW75TPmGPP1K6gDh6QAAmQwGYJ5CrY04/ukPj+Kx3+4z1sYP6iswra/+pcKYuVDtpwBtwnfZxYCvsvwtXlOvf44tjZwOVODIk7lMP/ZgqbZqU3EomvH/jWxbEWuufeuoIs7yvHBXtRP4bmzzya2S4/xzd305bEa2v+91na5ZfB3yRAAiSwYQLrEezYPKne8azhvzNtnaD/w21VdEgcwJpf63LF2cLh797uXec1773tuzkkPu9airPqX9io6XfhzQ1GMMXopd7pzEL1eWyns4vImv1kVd5Dwf4z/Yrdel/rSlYWO/OxCaWkEQQ/XYZ2+UXwNwmQAAlsmMB6BFtYkWzdO0T9rIdrY2+OhGDrxgobp7Q+jBMrx5MbbrTQ/5ZcUHRzqVc7u/tj601YWh97aN7R17rSriWTldJ7icf85eZ39wRffi/xKUavj7GnF7EZu90t6sc02+d9715nwoYvkUyrtitSCD+QAAmQwMYIZBTsjdnJigwCktOMr3lIAiRAAiRwCwlIsX+JrUlvIZEtaJLktC0wmyaSAAmQAAlkICDFfgp2BqCbyCo5bRP1sg4SIAESIIHiCEixn4JdnD8Wqlly2kIZmYgESIAESGBrCUixn4JdcndKTiu5yTSPBEiABEggIwEp9lOwM0LNO7vktLzrZPkkQAIkQALFEpBiPwW7WJ+k1i45LTUTE5AACZAACWw1ASn2U7BL7lLJaSU3meaRAAmQAAlkJCDFfgp2Rqh5Z5eclnedLJ8ESIAESKBYAlLsp2AX65PU2iWnpWZiAhIgARIgga0mIMV+CnbJXSo5reQm0zwSIAESIIGMBKTYT8HOCDXv7JLT8q6T5ZMACZAACRRLQIr9omDrhPwhA14DvAZ4DfAa4DVQ3DUQf2QQBTueiJ+LI6BvFv4jARIgARK4WwSk2J9QAynR3cJUrtbSH+XyB60hARIggU0QkGI/BXsT5DPUITktQ3HMSgIkQAIksAUEpNhPwS654ySnldxkmkcCJEACJJCRgBT7KdgZoeadXXJa3nWyfBIgARIggWIJSLGfgl2sT1Jrl5yWmokJSIAESIAEtpqAFPsp2CV3qeS0kptM80iABEiABDISkGI/BTsj1LyzS07Lu06WTwIkQAIkUCwBKfZTsIv1SWrtktNSMzEBCZAACZDAVhOQYj8Fu+QulZxWcpNpHgmQAAmQQEYCUuxfQrBv0D3ytmh72MF4pjFD2JabrnZxMzNVbl9c2VCqhu733GrYaMGS07IZMEb3iQXZN1NcXzRQvWc5W9PuHNRhf5jt6dCOiZOv/Sk8U8TR5P0JLNH3q7ariFawThIgARKAE4PjHFYTbHWIztd4Ud7nKxuWtxe5LAoz8q3rNAV7NslfY/T+qDgXguSb0Sv93Q5qZ130L3to/34IS1k4/jvlwet7FzWlYF/Nrjr3b751cew8KCYf1lZuV+5GswISIAESkAlInbXlBfu3GmqWwuEbuec1PLNgPariUKkZvTjZuLWdpWCLKCefu2g8cHvO+kJICPaPHupKoRrx6xTDF/tQVhODqVise7JwwR47oz+VA/3AERPsLO2a02R+RQIkQAJ5EliPYB910T2zoKRh8ekATcuCfeH2uHxRuLmoQSkbw1jrEud/XqN7WkNl1x1S3zmooXFxjUks3+RjGyeP99yevLWH6rM2hj+8RJJg/7rB4LxulFuH/e4a8zQoVmVhHyWnScYMzxTUURdiX9gTVLVbQ/tzH7Yg2Ho4OSF2uqLPLewra3bv2eFt/jWb0M+Tf3uwn1aw44y47KDy1Ebvy3LU57bLATHF8KzitH38MTkdsnK7JMg8RwIkQAIbIiDF/uV72EddjJ0gnhwWn142YVk2ht9WEWxvjvygjtZFH4PLPrrPq44oV9+GvfnxxbFzbufIRvfDAP0LGzUt8Ac2RloL4oI9HaGle5a7Ndh+uWc1R0QqZ8PSi7bkNOl6mSts33tonHZx7Tz5DEXBHr3Uve9QbIM6pgM0tMD/Ffog+E4f/LjG4O+mM6JS/3OAwaX7gDW9slFRCtbjpuOnwbs26k4PvwL7anHRntsuAE491jG63wTfA1i5XZFG8gMJkAAJbJaAFPtXEuwbjNC6Hx8Wn2Lw3IJ1NgS8Ht1SPWwvT/OjCWWK/h8W9p523UVuTg9ewTodRIX2awfV3Qqa/0wTgj1+W4WyGujHuunTfxqwVHX2XLxpRoHHktMkc9KELcwjC7aT/zePc5gYgJteab/O+pcYEnevD93jj8j8r2u0HymoR/MWLUYrmduuH32cWBZO3nvOjT+saev1yMOq7Yqawk8kQAIksDECUuxfUbCB0fl+dFjcEdN9tD5jNcH2xfixjf6XSVSQfUSf9IK2OcOzOl0kaI/R/U0PFbfQv9S9P+PH6xX6DxV+FWX7LTnNsdEf5vYW+Ol08R95EZgk2FMMTmcNqa8g2M4IjEJDP0DF/k3e1aHmLVpcuF3uavfIKEnE97rijO2K2c6PJEACJLApAlLsX1mw3bnNcFhc91jV/RZGujVe0PXFMDFX7bU4fn78voFD75UwZe3h8KmN7uU4mMO++ftYnmc1CUaCtic2gpj54ma9dCw2SyjVseQ0x8DJEJ3zFlrez4nuud4/RtM41xdX8kuCnbEnGu9hOz4Ir40IUOe7OSvKF2zX+E3VmQYZ/jRKj/jePc8etsGHhyRAAltDQIr9qwt2ZFjc7cnsn3vit6JgOyR/TXB92UXrWbj4zHriDq26Ah9bBRzHHwnaC/QO4/lL9llymmTi3KHjSAZZsDPN9a5TsCO2eg8SicV06Q9i/nx8pnbFbOFHEiABEtgUASn2ZxBsb1hcz0c6C5O84XDdGlGwTxLzyDODaUBkivFbvcLcK3vmkPgI9u4eam+vE0PinYcqOnQflL0dB5LTJMuzCvb81dQK0bUFMQvigp06JF5zF4nFipE+yu3SD3XG9IZ//Kc73N78O1z8lqldkkE8RwIkQAIbICDF/kyC7Q6LV9F+bQyH64bEBNsNmoag6zQ/h2je1/Ou3srkqxYO71XcOXADxvSDft3Iy+vPc8cWnU0/NsNXjyI9bMBZdKb20fwYnU918lh77kI1o76yHUpOk2yUhU1MKa4SR5b3leOC7Y2+zFx0dtDGtWSacG7xdsXXL3iFZWmXYA9PkQAJkMAmCEixP5tg+4FZKey/MuaCY4INZzWvfvWqjva7AdxXfHZw/MR4P3s6gn2gnNevGm97zgKx3tsGqpaC9aQXvF/sv9ZVedpG73KA/lt33jtIExNs/WDglKt38DrX5fbRPa87c+X+UPsm4K9ah+S0Vcty88lD4nqB1uil3njEQvV5bKezi8ha72T1kz5O9DXwrIP+ml/rSlY250zc907SDO2aUxW/IgESIIE8CUixP6Nge8Pifg/Ytz4u2Hq97pcuGt5mJ9a9qrMhin6HNuhh67yTETrG3LXaraB+PsDNL79g93dk4xQnzTBYmBZdJe7l+3mN3lm4cYp171AsN1pLOT5JTstm2SzB1qW6e4Ivv5f4FKPXx9hzFgweo+ft3pLcOKWF/rfoSEe2tsRyi4KdpV2x8vmRBEiABDZEQIr9Swj2hqxkNRECktMiCfiBBEiABEjg1hGQYj8Fu+RulpxWcpNpHgmQAAmQQEYCUuynYGeEmnd2yWl518nySYAESIAEiiUgxX4KdrE+Sa1dclpqJiYgARIgARLYagJS7Kdgl9ylktNKbjLNIwESIAESyEhAiv0U7IxQ884uOS3vOlk+CZAACZBAsQSk2E/BLtYnqbVLTkvNxAQkQAIkQAJbTUCK/RTskrtUclrJTaZ5JEACJEACGQlIsZ+CnRFq3tklp+VdJ8snARIgARIoloAU+ynYxfoktXbJaamZmIAESIAESGCrCUixn4JdcpdKTiu5yTSPBEiABEggIwEp9lOwM0LNO7vktLzrZPkkQAIkQALFEpBivyjYOiF/yIDXAK8BXgO8BngNFHcNxB8ZRMGOJ+Ln4gjom4X/SIAESIAE7hYBKfYn1EBKdLcwlau19Ee5/EFrSIAESGATBKTYT8HeBPkMdUhOy1Acs5IACZAACWwBASn2U7BL7jjJaSU3meaRAAmQAAlkJCDFfgp2Rqh5Z5eclnedLJ8ESIAESKBYAlLsp2AX65PU2iWnpWZiAhIgARIgga0mIMV+CnbJXSo5reQm0zwSIAESIIGMBKTYT8HOCDXv7JLT8q6T5ZMACZAACRRLQIr9FOxifZJau+S01ExMQAIkQAIksNUEpNhPwS65SyWnldxkmkcCJEACJJCRgBT7lxDsG3SPFNRRFzdzDLm5qEGpGrrf5yTiVwsTkJy2cOYZCSfvT2CJPpri+qKB6j3L2Zp256AO+8N4Rinm6YmTr/3JPLeh4x9DtJ9VsWfp7QMt7D1uoPvvJFb5qu2KFcOPJEACJLAhAlLsp2BvCP6q1UhOW7UsJ9+3Lo4dcUs+VI1eVaDUDmpnXfQve2j/fghLWTj+e94jGoDvXdSUgn2VybLlM3ttsR7U0X43wOBDF83H+mGjgtbnsLiV2xUWwSMSIAES2CgBKfavXbA32qI7UJnktNWbPXZGSSoHWphjgv2jh7pSqL4xe9RTDF/sQ1lNDKZzai1EsKcYPLegDmwMf5q2jdC6r6CeD9yTWdplFstjEiABEtggASn2r12wI0Pin1vYVwqNf5LRfnhmQd1vYeQDuBmg9bSCHecvhe2g8tRG74uZbwhbKdT+GqL3+57T89s76uDaz39Lf0tOk5o6PEubrphieFZxpjTGH+2EYOth8oSI64ocH1qze89Xuizzr9nYGHoGTv7twZ7rU6kl0XMz2zXp42TGtWWWsHK7zEJ4TAIkQAIbJiDF/nwFG2N0Hiqo0wFM6QWGsC2F/VeuXE8/tXBoKewc2eh+cIc27aMdZ2jTvvJzuoJtWRYqf3ScIdvO+9su13DEcJHrZKaweZmnVzYq1jG63wA4IhvtYY9e6qHkUGyDOqcDNJwHJbPnHXwL/LjG4O8mDpVC/c8BBpfX0DPITn1KwXrcdH36ro36A3e4OvSpUc6Mw5nt8tvwbYLh6xNUd/VDg37Qa2FgjOCv3K4Z9vA0CZAACWyCQAGCDYzfVqFUIzqkemXDUofofNXNHqPzSME67TuBPgQxxeDUgnrUgSsVrmCrh+1b36sOGaxJsH/0cWJZOHnvLcbyxc5YGOgI429dj7Vpgcf9zO83m995x4khcW9Y+ihW3q9rtB8pw6dCWbFTswT75u9jKLWPwwc7xkNBCzUt3NYxel7bMrUrZgs/kgAJkMCmCBQi2PjaQTU2dOkMh/tC/M1dsFQ772NwqXto4U/vxaExTLuAcGyK5AbrkZzmVO+JpP5+1o+7CGyM7hMLlbNhOMqREGz9cDRrSH0B7nHBnjMVMnlXhwoe1gSQC7bLnXpxbY70/W96OLEUrOd6VCdjuwTzeIoESIAENkFAiv05D4nrZnmvgwXD4u5wePWtF2Yd8ZgtOvpVHdt5XWgB4dgExQ3XITnNMWEyROe8hZb3c6J7rveP0TTO9b8C4zfV5MKshGADmXqiccF2yvdHUGLAPH/PXFG+YLt8wQ5GDYxqhi90L9t21kdkapdRJg9JgARIYJMEpNi/AcEGnF6Vv9LYGQ6vesPh/nzqIq8EUbDnXSzy0LHHbE4v3J+3zjTXu07BjjVSbhcw/afhjCxIwu/k8ebjM7UrZgs/kgAJkMCmCBQm2HBerbHQvNQrla3o5itfO86CpcPI60QSEgq2RMU/JwvbBNfGFEMw3fCnOyzd/DtcJDZ/NbVC86Nfk/A7LtipQ+I1d/GbUFT8lNwuQE+16IVutQtjhZmT2RsGf+iufcjUrrgx/EwCJEACGyJQnGBjgv4zBeuFjaalUH9n7kTlLjpT95ux92ndd4Cte00MnPdsKdjzrpOZwiZlEobE3Yeqdb2HnbLo7GDxhYOz2+VdN/H3sL01Efvn3guDfA9bugJ4jgRIoOQE1iPYsXlSfw61dd7BcAJE3sM2gLg9HT1XXUfvh/GF8QqQ2q2h5e1Y1fJ32brwlxRRsKPUMnySBBtTjF7qDVUsVJ/HdjoLfDCjTu+d6P1n+nW79b7WNaNG53Tw6pi/09k7b5V4RMQztGte5fyOBEiABHIksB7Bnjkf6r7XO0uw4b3Pq57FX99yWzz9Ym6yYWHvQR2tS3O4k4K9tmtDFGxdursn+PJ7iU8xen3s7ed9jJ7ntuTGKS30v/nv1a+nNZHrxtpD9VkHI3MAx6lm1Xatx0aWQgIkQALLEsgo2MtWx/TrICA5bR3lsgwSIAESIIHyEpBi/xKrxMvbsNtsmeS029xeto0ESIAESEDeNIuCXfIrg4JdcgfRPBIgARLIgYAU+ynYOYBeZ5GS09ZZPssiARIgARIoHwEp9lOwy+eniEWS0yIJ+IEESIAESODWEZBiPwW75G6WnFZyk2keCZAACZBARgJS7KdgZ4Sad3bJaXnXyfJJgARIgASKJSDFfgp2sT5JrV1yWmomJiABEiABEthqAlLsp2CX3KWS00puMs0jARIgARLISECK/RTsjFDzzi45Le86WT4JkAAJkECxBKTYT8Eu1ieptUtOS83EBCRAAiRAAltNQIr9FOySu1RyWslNpnkkQAIkQAIZCUixXxRsnZA/ZMBrgNcArwFeA7wGirsG4povCnY8ET8XR0DfLPxHAiRAAiRwtwhIsT+hBlKiu4WpXK2lP8rlD1pDAiRAApsgIMV+CvYmyGeoQ3JahuKYlQRIgARIYAsISLGfgl1yx0lOK7nJNI8ESIAESCAjASn2U7AzQs07u+S0vOtk+SRAAiRAAsUSkGI/BbtYn6TWLjktNRMTkAAJkAAJbDUBKfZTsEvuUslpJTeZ5pEACZAACWQkIMV+CnZGqHlnl5yWd50snwRIgARIoFgCUuynYBfrk9TaJaelZmICEiABEiCBrSYgxX4KdsldKjmt5CbTPBIgARIggYwEpNi/hGDfoHskb9G2c1DDyeshJhkNZPYkAclpyVRzzlzZKdvM2hgG2ae4vmiges9y8uwc1GF/GAffhgcTDM/rqOzq68HC3uMTdD4t4f1vXRzvNjD4GZa4lUcO2xq6313rby5qUMrkWXyrXJtCG4u3yLDgcwuVRx1IV5iRyj2cXKP73zZGiS/yPjHB8PUJqs61rqDvidbljVDpgvfEjyFaTyvY0ds/W3uoPutgJNw60y9dNB7vwdLpdiuon/Ux/iVUK54ao/tkB41/puK3eZy8uWzBfhdyye+6G8JWCrWLsK5Ee7yYZ18lvin+xI8+ThaMfVLsX16w7x+jed5CK/hp4sS7sKxnfYr2mi8JyWlLVfG1b/jK8NtpDXtKwXrSg3/pj15VoNQOamdd9C97aP9+CEtZOP7bT6FrnqD/zIKyDlF/3cPgQxf20Q6UqsD+tEiAuEHviTX/hluqgQUmpmBnhD/F4NRC9U26ZBfzMDTF8EzfExX3WjfviQvT5gXvCR2sLQXrQR3tdwP0L2zU9IPAgY2ReevoBxmlsHNko/thgN7rOg51PuNenQf+5u9jWEfd4L6el3Y93yVFlII9m+z4TRXW6QCmy6XUUuxfXrDFC2GK4Yt9KHWMnhnbJSt4bikCktOWKkBM7PnrwMbQ7+X+6KGuVCx4eumsJgb+1fW5hX1lxZ7ex+7oywK9pellE9Z9G0O/PNG+LTlJwc7uqK8dVFUdvR/ziypEsJ1rPd6bcx8ylGWHvf0F74nR+T6UFRtZ+tZFTd93b/0HgAl6TxVU7F6afmw6913zMuXGmQ7QtPZhX6Wkm497yW+Tgr1kAUskX6CuMvewdUt/DtCw9tH6PL/ZUuxfk2ADwRPVN8OIyTV6Z/7QqTucZL+7Np4svGH2syEmH1uoH+iemoJ1r4rGhZnOLXPyrzFM5AwntTFMudENa7byUHKa1JDhmYISH6aSqaf/NGCpfTQ/hjf15P0JlBKGTr1g5A8vjV7poJMc9p2ZP1L9GJ1HCvuvkgObk4/tYKTGHSqM+fbHmxSTxQAAD51JREFUEO1nNW8Y3h1OPHzaivg/COrxYcfTLq79BxPPnumXHmx/aFLtoPLURu9LyEMnS73eFhHsJe6BCCqE90Zw3mnXIfYsbyriwazh2SBHeF9+GqL1H3eI1b+/4iOxk3/TmPjB8joxddL6mCgtkkbXKU+bue3cP09eE34rnGvb/AuCZ+Ekzs1lynXjFzLr93dXMP3r20x2876Bw3sn6Mea5twDxtTHYvfECK37CpZhu1uX7p0b9+6kjxNxyHdWftNiYPy2CnW/FT5M+F//TIvF+npf0P9/DdH7XV9LFvb+1+94aPrGiyGBHnjTRY4ZqTZMcf3ORv2BNxWgRxkOajEt8K/BOT1DX7A/XKN7WnWnH/S0wvkAN/FphSXuz3Vq1PDMgnramzsiLcX+tQj29PsA9qPYkM3PIewDHViraF70MQiGkxQqZ0NPtL2gdFBBZbcGW6cLhlgtmE+TN+/q2FEWDn9vo3c5wOBdG/UHemj2GF3zIcG/QG/Jb8lpUtMWFuzpELYOHLHpi9FLPW+dFGJMB2joAPKX7gHo4WwjuJiGfO3gUCk0P5onY8dOmuST5fji2Jmr84cAE0OFzhNpOJQ4uOyjq+fQdaA4CIOTGyQqqBz4w/p9dM9qzg1rPQ+HoKZXtpPXH5o0r6WePx+9yPWWJtjL3AOJQB4TbN9vjxvovBu499NTPVwbvU9ixD3BtmBZFg7/20Ffs3tedXhXXo6Ch+eAyeOmMwwbMNFTHUFvzQ2W+wcV7HjDumG6Kjpf/dpv0Hu6E06bXJrDut3EnLXjN3MUxy/G+z35MkDvxSGUqqOt7/0vroL6103lqRsTgutmmZgwR7BjZgC/phhf2qhaCiG7Be+Jmx6ORSEGxm9025ruWpJPtuOb5AOE7tkrqN+S/EI7x+g8FB6Ig+vQm8YKYnE43bWM//W1VPlDX0s9dP73/8Hosu2Mzh2+6GFwOcLN1OjA+YK9gA16qDicktMxvoPGYx2XLJy895+aFhdsbefOUcvRCz2toONFZFohsKkAjbrSfjbvmdCL/pEU+5cX7MjTlH7S937+08HY6KA4wz/6aSsmptd/aqf4hnpBKZFuCFv3Ivwg5gzzKFT/vPbb4v12e2xmMI4l2PqPktOkRi0q2M4TeMA/LMnJLwYD9wZxfeH560XYwwlK8ALfvMUgkcDkZ/R8m5jT0UOluxU0/5li8r7hPND5i7v8rG5bwlEBV7D9hws/FeA8zQYPI+41Ex9ydIapdvdQe3sN/ZDStBa43lIEe6l7wL/WA7Njgu3XFbmfrtF+uIPKueAPrxyfSXye2A2OPju396ZHaPyBWSf7r2u0H5nDs961cL8ZTqXohDHfO9Meqor2l6Ax7oEz/C08YPijOJ9i6Y2PbjuMB0pvCmf/hf/w7yWe9NHQ873GA5pRTPJwQcH2Oer7ce+/faOntuA94dUjPdC6ZXu+EP3smu2ONBgM4q2Z8dAs3/P6AWAHe0cdXGNJ/z9sIxqJkyIaaRO8nn8i7kRtaD/cMzpzXuO8ezHQAiTrimOAw1Al5ondkUUL/sPQUvfnujXKG0mJrg+KtkSK/csL9rxFZ49bGDlDj94FIE2sezfaobPQxLvYE0LhnfeHeD82odQeGm91zyL60/5dPzDMuYijDLbuk+Q0pxFeANDfz/rxL8yw0d6wWqx3DXhP7z7vMAPg3SDuDeOJXUJckkE7UoTzYUYdTo8ivImS+eaccW7Mw6BnFwSJiKjFnvbn9HSCmha93vzg6vfKI6vEl7wHEkxjgu2InXJ6NsNvxpNxYLR84DIR5oi9uVMnYHhztdKq4sm7OpTyGXuCnXhg8857bRi+UFD/00Andq8OvJ5YGHw9m2cOA4dtctsR3ufuFIxvV5hOH7mBuBGuu4h+DVf4Zt830tTS+JMeJQxHdqwn/sPNgveE57/kPRm9PqeXDXlqSt+JetprTqzTgpSc1pozAuBzWdb/iWs1KaLBvejcGwvY4NuS+O3dB4GWJOtKZHHuy+RIno5lTd3LfqmnX5a8P9euUV79CZZha6TYv7xgi0EdmFzqRREKldf62cuF6opyaIB7ZN7cMWEOkkbPu86fc4Pd4sVuktMcTJMhOsFK/RZOdE8o9jDVD4YoXbBuz0fo4fjBIHFR6nyCvxIBexHB9nwae1jQK1qTQSa4EKIHPycYf3aHylqn4doIPwhGg0SYNXJ+5pBjPP0C19tcwV7yHkjcuB6v4PwUo9fu8L7zgLZbQe1ZG73Pc+by9CRG5CEibKPvV2c9QezBx0zl91Zcxua1YKYyz3t2z3mQVInVzmZ+s9zwON4O97MsyhF/h0UER+P3xtsSz4+duFV9Zpx7O/8VVbenplB/p4dpvfam3RMb6GG77Y7PuY+io5UBBeNgTf43R9eiPljAhsCcKSY31xjqh6PX4RtI4UOUe62YdQVZ/YPYfemfDnzlxCC3nOI0yrtuZuiptlmK/WsT7OCJxTEgDxj+8F2I/y4cSU6T2u08fc9xvtOLfq7n/MNeillO/nPYcQFya4/e2KZFxvHPEdrOq2NaRC3sPThE7VkTnTPd+1PBENessiLnnZs5zGPUEhxG0gdnhYNYYHDz+XyXvAcCYfbrkXlheoORntt76i8+M9eE+HnD31GbwvO+YDsBa00BOzJtMvdaNO3Qxy6rRM/bSBZvh/t5NcE2ig2G8/2Hvsh3Mz+4PTW/vQut65gzshOZKpr5QOmNUIkP1a6hcUbu2XS27kOZPFqxzAObKaKuLX7MXsAGPWyuF/k5Cyrdd88PH9bReN1C46G5bsYty6wr4abYfRl+791TTm/dLWe9gu23N6xx9lGpBDt9uMFfxORsxpK4uaONcYd6VOxVotkobtM36xNs9yk3uUrVpTVzlbc3XObPvaWuiI0NR4e+8Hwa62Fj5pD4CLY3p+zPNbU/+wtP3FL968IPttEgYdTs9DK9m2lm4NTDdjuovBhg4gwtLnC9xQJDNGAueQ/EuWCBXsmvCQbO65TxXlW87cL33pC4s5gndUjUX4syK/Ca5z1hUbKYhpaZR2Z+83x4HGULpA+Je4u4wiLkozlz2KPXFexIK669hZju0Cqw2D3hXg/J+2/irhL3hXjm9IB3PUk9ea9lLqO4r2cPR+vRrZ2DJgb/V7+qKV/v7pTIYv43RTR6Ly5gw/9zXyvdP+07i9ZCZ7kL6ZbvYc8eEnffUlny/ly7RkU1LmxveCTF/rX1sN33BMMVin6QlRedVdD+Vxs2y+jYeX/hQXxRDNz3f3ce64UTt/Of5DSppak9bDNASwXk/h72jB7CjEVn/nun9pV3LfgBLbB9gv6puyPbUoKNGYvOPD5O0Fn0epsr2P5cqh/sAsPhLrz07wEvYMeEwb+f/F7n+EJPASTnoq//1CuM40E6rMsNnPGFeP6GIH55XvCK31/+orMDf5HRLGGNnnenXuJ1AnAY76CqF/aZ/2aKVJgoLthIW3T2x4KbOM0RbPehIL7nAOCuTjemlpwHnni65N4ETkzM8T3sYA479tAsLzrz7HOEaD3+ny3YsxadGTY495LwlsmXNqp6eiUQTPdaM+sKrxLvyCsrviDR9Vso5MVqlMc8MbIWtkaK/csLdmyeVO941vDfZ7VO0Pffi15myXzgDN/YmGDr4RLv1R/rQQOdD96Sf3+HreC1Ez//7fktOW2l1jkLqcKLNVnGFKOX7mtC1eexnc4iuzrpp2UtlN7uT8FreObrP8nS9ZnxX3rrzmTPy/dt8HrOW3dozH8Fw92BzX8tyX2l7+TRDqzdncgrMNGn+tCG+PnkKyytxI5Tvk1zr7cUwcZC94DbW9RbUPrtd3a22j3Gsd4K2L+hv/dw7O2S1XJekwwXQIWvGIVt9o/ctpuvuhk72Bl+TTLxXpsUXusKbPIrSQxp660x3ddx3FfJBui9bQSMg816/PyO4M27Nv0e9T5O3vTX+1qXb4P422uHv6ufnlf1XhMMX03VGRe8JzwfqoPkTmcmk+kn77VD7xW7cKczf6GbaKz3QCT0lIXrsOXHTm93wmz+98TnkY1e6mtd4StUERu8XeDUbg0t57VFzfoYe9YOdvRucIFGLCrY+87rne5rwDP8JnBxd3c0p5mSWuTST55fKGb4rvMeUsPX1fwvwt9S7F9esIXFJNa9Q9TPeriOjlgCwkvprQ/j4N3PhXvYXhsim2vozS6OGuj+G680bPBtOJKctkq74qIllzGJbHYxcy/xXzcYrLKX+Iz3sLUtEd86mxyYC3/MfZq9zRTeDnFz424y4e8SNauN0vnIxil6Ex5hc5WITdL1libYTsOSG1ZE7wGdyNzDXe/N3kD3i+4FG4Ktk33rwzbmrmdtgKKT+v/cttsYfusHe1Nb947FPeKTG2e00I+sSI/2pP065Dlobw9ub196vR92TTMWbldnDnfG2oqgDr2Owdv4xVy0ltg45bSDkd9pCDJnOIhc6+HmT4lmRNLN2V//ZrDQXuKRTXsW3kvcew9b2oQmsmlJeI2ZZFb3vzv/7O63bsH+5C92jM3pptgQ2T/d2oPeGKn/deJNOfi7LS4q2DV0v4zQ9juTu1WcvB0lNyopSKOC0RD/PXXTEd6xFPuXEGyhRJ7KnYDktNwrza0C96l03q5WuVXNgktKwJ2i4DWxHvc4D2ix6ZX1lMxS1kdgioFeAJzvTmfrM5clLU7gdgk24MxvztnVanEyTHkrCDjD4f5c+q1oUbGNcNZfGPPrxVrD2iUCzvqL+VNAOpsU+9nDloCW6JzktBKZt4Ip7nzf3EUjK5TKLNtIQC9EtFAR9pbfxtaUxebN/7WusrR8O+xw/lpX4q2QpO1S7KdgJzmV6ozktFIZuIoxt+XvYa/SduYJCeg/I2nsBR9+waNsBDb/97Cz2XuHcjt/D9tYnD2n6VLsp2DPAVaGrySnlcEu2kACJEACJJAfASn2U7Dz472WkiWnraVgFkICJEACJFBaAlLsp2CX1l2uYZLTSm4yzSMBEiABEshIQIr9FOyMUPPOLjkt7zpZPgmQAAmQQLEEpNhPwS7WJ6m1S05LzcQEJEACJEACW01Aiv0U7JK7VHJayU2meSRAAiRAAhkJSLGfgp0Rat7ZJaflXSfLJwESIAESKJaAFPsp2MX6JLV2yWmpmZiABEiABEhgqwlIsV8UbJ2QP2TAa4DXAK8BXgO8Boq7BuJPHAnBjifgZxIgARIgARIggeIJULCL9wEtIAESIAESIIFUAhTsVERMQAIkQAIkQALFE6BgF+8DWkACJEACJEACqQQo2KmImIAESIAESIAEiidAwS7eB7SABEiABEiABFIJULBTETEBCZAACZAACRRP4P8D5YHtyQdBsbsAAAAASUVORK5CYII="}}},{"metadata":{},"cell_type":"markdown","source":"Below methods we consider eveything below -1000 as air. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# Rescale intercept, (0028|1052), and rescale slope (0028|1053) are DICOM tags that specify the linear \n# transformation from pixels in their stored on disk representation to their in memory representation.\n# Whenever the values stored in each voxel have to be scaled to different units, \n# Dicom makes use of a scale factor using two fields into the header \n# defining the slope and the intercept of the linear transformation to be used to \n# convert pixel values to real world values.\n\ndef 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 <= -1000] = 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)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"windowing\"></a>\n# Windowing\nWindowing is the process in which the grayscale of a particular image can be adjusted. Windowing uses two values window width and window level for getting a right/interested window. The window width is the range of the grayscale that can be displayed. The center of grayscale range is referred to as the window level.Window width controls contrast and window level controls brightness. That is why windowing also known as grey-level mapping, contrast stretching, histogram modification or contrast enhancement. A large window width means there is a long grayscale and the transition black to white will take longer and vice versa for smaller window width. An example to better explain:\nWW of 100 HU could mean the grayscale only ranges from 0HU to +100 HU, with a WL of +50 HU.\nFor Lungs window width is 1500 and window level is -600, so grayscale ranges from 150HU to -1350HU.\nMore about windowing can be read [here.](https://radiopaedia.org/articles/windowing-ct?lang=us)\nHow to get window width and level is explained briefly in [this](https://www.researchgate.net/publication/272179277_Determining_effective_window_width_and_center_using_different_windowing_techniques_for_radio_therapy_images) article."},{"metadata":{"trusted":true},"cell_type":"code","source":"def apply_window(hu_image, center, width):\n    hu_image = hu_image.copy()\n    min_value = center - width // 2\n    max_value = center + width // 2\n    hu_image[hu_image < min_value] = min_value\n    hu_image[hu_image > max_value] = max_value\n    return hu_image","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"histogram-analysis\"></a>\n# Histogram Analysis\nLets plot histogram for image pixels after converting to HU and raw pixel values. After converting to HU we can see there is lot of air(-1000) in the scan. Some fat and muscle is also seen."},{"metadata":{"trusted":true},"cell_type":"code","source":"train.loc[0]['Patient']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig,ax = plt.subplots(1,2,figsize=(20,5))\nexample = train.loc[0]['Patient']\nscans = load_slices(f'{TRAIN_IMG_PATH}/{example}')\nrescaled_images=get_pixels_hu(scans)\nimages = [scan.pixel_array for scan in scans]\nfor i in range(10):\n    sns.distplot(images[i].flatten(), ax=ax[0])\n    sns.distplot(rescaled_images[i].flatten(), ax=ax[1])\nax[0].set_title(\"Raw pixel array distributions for 10 examples\")\nax[1].set_title(\"HU unit distributions for 10 examples\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"metadata\"></a>\n# Storing metadata in dataframe\nEvery DICOM image has lot of metadata as we saw earlier for one scan. Lets put this metadata in dataframe for easy access. "},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_dicom_raw(dicom):\n    return ({attr:getattr(dicom, attr) for attr in dir(dicom) if attr[0].isupper() and attr not in ['PixelData']})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n# Get dicom metadata\n# Image features like lung volume are implementation from a detailed discussion \"Domain expert's insight\" by Dr. Konya.\n# https://www.kaggle.com/c/osic-pulmonary-fibrosis-progression/discussion/165727\n\ndef get_dicom_metadata(df):\n    patients = df.Patient.unique()\n    dicom_metadata = []\n    for patient in patients:\n        path = f'{TRAIN_IMG_PATH}/{patient}'\n        img_list = os.listdir(path)\n        for img in img_list:\n            image = pydicom.dcmread(f'{path}/{img}')\n            record = get_dicom_raw(image)\n            raw = image.pixel_array\n            pixelspacing_r, pixelspacing_c = image.PixelSpacing[0], image.PixelSpacing[1]\n            row_distance = pixelspacing_r * image.Rows\n            col_distance = pixelspacing_c * image.Columns\n            record.update({'raw_min':raw.min(),\n                        'raw_max':raw.max(),\n                        'raw_mean':raw.mean(),\n                        'raw_std':raw.std(),\n                        'raw_diff':raw.max()-raw.min(),\n                        'pixel_spacing_area':pixelspacing_r * pixelspacing_c,\n                        'img_area':image.Rows * image.Columns,\n                        'pixel_row_distance':row_distance,\n                        'pixel_col_distance':col_distance,\n                        'slice_area_cm2':(0.1 * row_distance) * (0.1 * col_distance),\n                        'slice_vol_cm3':(0.1 * image.SliceThickness) * (0.1 * row_distance) * (0.1 * col_distance),\n                        'patient_img_path':f'{path}/{img}'})\n\n            dicom_metadata.append(record)\n            \n    metadata_df = pd.DataFrame(dicom_metadata)\n    metadata_df.to_pickle('metadata_df.pkl')\n    return metadata_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nmetadata_df = get_dicom_metadata(train.copy())\nmetadata_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"voxel-size\"></a>\n# Voxel Size and Volume\nVolume of scans vary highly. Does it mean that lungs of images having more voxel volume are bigger than others? NO certainly not. Let's understand what is voxel? Voxel is 3D pixel having pixel spacing shows distance travelled by a pixel in x and y coordinates and slice thickness as z coordinate. Voxel we can imagine like a cuboid. So if pixel spacing(size) is more than definitely volume of image will be more. Lets plot variance of area, volume and pixel spacing to have an idea about number of images having large volumes. "},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.tight_layout()\nfig, ax = plt.subplots(2, 2, figsize=(20,10))\nsns.distplot(metadata_df.pixel_row_distance, ax=ax[0,0], color='green')\nsns.distplot(metadata_df.pixel_col_distance, ax=ax[0,1], color='blue')\nsns.distplot(metadata_df.slice_area_cm2, ax=ax[1,0], color='pink')\nsns.distplot(metadata_df.slice_vol_cm3, ax=ax[1,1], color='magenta')\nax[0,0].set_title(\"Pixel Rows Distance\")\nax[0,0].set_xlabel(\"Pixel Rows\")\nax[0,1].set_title(\"Pixel Column Distance\")\nax[0,1].set_xlabel(\"Pixel Columns\")\nax[1,0].set_title(\"CT-slice area in $cm^{2}$\")\nax[1,0].set_xlabel(\"Area in $cm^{2}$\")\nax[1,1].set_title(\"CT-slice volume in $cm^{3}$\")\nax[1,1].set_xlabel(\"Volume in $cm^{3}$\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# It is clearly visible that area and volume of lungs vary a lot. Let's show images with maximum volume and minimum volume.\nhighest_vol_patients = list(metadata_df[metadata_df.slice_vol_cm3 == max(metadata_df.slice_vol_cm3)]['PatientID'])\nlowest_vol_patients = list(metadata_df[metadata_df.slice_vol_cm3 == min(metadata_df.slice_vol_cm3)]['PatientID'])\n# Load scans for highest and lowest volume lung patients\nmax_vol_scans = load_slices(f\"{TRAIN_IMG_PATH}/{highest_vol_patients[0]}\")\nmin_vol_scans = load_slices(f\"{TRAIN_IMG_PATH}/{lowest_vol_patients[0]}\")\n# Convert to HU\nmax_vol_hu_imgs = get_pixels_hu(max_vol_scans)\nmin_vol_hu_imgs = get_pixels_hu(min_vol_scans)\n# Apply windowing]\n# We can try with different window width and levels.\nmax_vol_window_img = apply_window(max_vol_hu_imgs[20], -600, 1200)\nmin_vol_window_img = apply_window(min_vol_hu_imgs[18], -600, 1200)\nfig, ax = plt.subplots(1, 2, figsize=(20, 10))\nax[0].imshow(max_vol_window_img, cmap=\"YlGnBu\")\nax[0].set_title(\"CT with large volume\")\nax[1].imshow(min_vol_window_img, cmap=\"YlGnBu\")\nax[1].set_title(\"CT with small volume\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"resample\"></a>\n# Resampling\nVoxel size resampling is an appropriate preprocessing step for image data sets acquired with variable voxel sizes in order to obtain more reproducible CT features. We found that some of radiomics features were voxel size and gray level discretization dependent. The introduction of normalizing factors in their definitions greatly reduced or removed these dependencies. In computed tomography, voxel size in a region of interest depends on both pixel dimensions (x-y plane) and slice thickness (z-axis), assuming slice thickness equals interslice distance. Any change in these two parameters changes CT image resolution or voxel size. A minimally curation step may be to resample image sets so that all have the same voxel size. In this paper, voxel size resampling was investigated as a way to minimize the variability in feature values due to differing voxel sizes.\n\nVoxel intensities within a region of interest (ROI) are typically resampled into a limited number of discrete values or bin sizes before calculating feature values. Different studies have used different gray level resampling before extracting texture features. Later normalization is also done to improve robustness of these features.\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"metadata_df.SliceThickness.unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Lets see thickness for slices before thinking about resampling.\npatient1 = train.Patient.unique()[0]\npatient2 = train.Patient.unique()[5]\nscans1 = load_slices(f\"{TRAIN_IMG_PATH}/{patient1}\")\nscans2 = load_slices(f\"{TRAIN_IMG_PATH}/{patient2}\")\nprint(f\"{scans1[0].SliceThickness}, {scans1[0].PixelSpacing}\")\nprint(f\"{scans2[0].SliceThickness}, {scans2[0].PixelSpacing}\")\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"patient1_hu_scans = get_pixels_hu(scans1)\npatient2_hu_scans = get_pixels_hu(scans2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def resample(image, scan, new_spacing=[1,1,1]):\n    # Determine current pixel spacing\n    spacing = np.array([scan[0].SliceThickness] + list(scan[0].PixelSpacing), dtype=np.float32)\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    image = scipy.ndimage.interpolation.zoom(image, real_resize_factor)\n    return image, new_spacing","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image1, rounded_new_spacing1 = resample(patient1_hu_scans, scans1, [1,1,1])\nimage2, rounded_new_spacing2 = resample(patient2_hu_scans, scans2, [1,1,1])\nprint(f\"Original shape : {patient2_hu_scans.shape}\")\nprint(f\"Shape after resampling : {image2.shape}\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"3d-plot\"></a>\n# 3D Plotting"},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_3d(image,threshold=800):\n    \n    # Position the scan upright, \n    # so the head of the patient would be at the top facing the   \n    # camera\n    p = image.transpose(2,1,0)\n    \n    verts, faces, _, _ = measure.marching_cubes_lewiner(p, threshold)\n    fig = plt.figure(figsize=(10, 10))\n    ax = fig.add_subplot(111, projection='3d')\n    # Fancy indexing: `verts[faces]` to generate a collection of    \n    # 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    ax.set_xlim(0, p.shape[0])\n    ax.set_ylim(0, p.shape[1])\n    ax.set_zlim(0, p.shape[2])\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_3d(image1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_3d(patient1_hu_scans)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"segmentation\"></a>\n# Segmentation\n\nSegmentation is most important part of medical image processing as it extracts region of interest. Segmentation defines narrowly what algorithm want to look at, so definitely CNN will perform better on segmented images rather than on whole chest image.\nSegmentation is done by many ways and clustering is most common among all. Clustering has several techniques such as K-means clustering, hierarchical clustering, divisive clustering, and mean shift clustering. Moreover, due to the irregular and fuzzy borders in most of the medical images, fuzzy set and neutrosophic set theories become important in the segmentation process to handle uncertainty in the medical images. Read more about medical image segmentation [here.](https://www.sciencedirect.com/topics/engineering/medical-image-segmentation)\nAfter clustering images are morphed using erosion(contraction) and dialation(expansion) to remove unwanted border areas and label different reasons separately. So the steps goes as:\n1. Normalization of image.\n2. Clustering for separating lung with everything else.\n3. Threshold image.\n4. Morphology - Erosion followed by dialation.\n5. Label different regions and define regions with different colors.\n6. Create lung mask.\n7. Apply mask on original image and get final masked image."},{"metadata":{"trusted":true},"cell_type":"code","source":"#Standardize the pixel values\ndef make_lungmask(img, display=False):\n    row_size= img.shape[0]\n    col_size = img.shape[1]\n    mean = np.mean(img)\n    std = np.std(img)\n    img = img-mean\n    img = img/std\n    \n    # Find the average pixel value near the lungs\n    # to renormalize washed out images\n    middle = img[int(col_size/5):int(col_size/5*4),int(row_size/5):int(row_size/5*4)] \n    mean = np.mean(middle)  \n    max = np.max(img)\n    min = np.min(img)\n    # To improve threshold finding, I'm moving the \n    # underflow and overflow on the pixel spectrum\n    img[img==max]=mean\n    img[img==min]=mean\n    \n    # Using Kmeans to separate foreground (soft tissue / bone) and background (lung/air)\n    #\n    kmeans = KMeans(n_clusters=2).fit(np.reshape(middle,[np.prod(middle.shape),1]))\n    centers = sorted(kmeans.cluster_centers_.flatten())\n    threshold = np.mean(centers)\n    \n    # Threshold the image and the output will be a binary image. Morphology workes either on binary or gray images.\n    thresh_img = np.where(img<threshold,1.0,0.0)\n    \n    # First erode away the finer elements, then dilate to include some of the pixels surrounding the lung.  \n    # We don't want to accidentally clip the lung.\n\n    eroded = morphology.erosion(thresh_img,np.ones([3,3]))\n    dilation = morphology.dilation(eroded,np.ones([8,8]))\n\n    labels = measure.label(dilation) # Different labels are displayed in different colors\n    label_vals = np.unique(labels)\n    regions = measure.regionprops(labels)\n    good_labels = []\n    for prop in regions:\n        B = prop.bbox\n        if B[2]-B[0]<row_size/10*9 and B[3]-B[1]<col_size/10*9 and B[0]>row_size/5 and B[2]<col_size/5*4:\n            good_labels.append(prop.label)\n    mask = np.ndarray([row_size,col_size],dtype=np.int8)\n    mask[:] = 0\n\n    #  After just the lungs are left, we do another large dilation\n    #  in order to fill in and out the lung mask \n    for N in good_labels:\n        mask = mask + np.where(labels==N,1,0)\n    mask = morphology.dilation(mask,np.ones([10,10])) # one last dilation\n\n    if (display):\n        fig, ax = plt.subplots(3, 2, figsize=[12, 12])\n        ax[0, 0].set_title(\"Original\")\n        ax[0, 0].imshow(img, cmap='gray')\n        ax[0, 0].axis('off')\n        ax[0, 1].set_title(\"Threshold\")\n        ax[0, 1].imshow(thresh_img, cmap='gray')\n        ax[0, 1].axis('off')\n        ax[1, 0].set_title(\"After Erosion and Dilation\")\n        ax[1, 0].imshow(dilation, cmap='gray')\n        ax[1, 0].axis('off')\n        ax[1, 1].set_title(\"Color Labels\")\n        ax[1, 1].imshow(labels)\n        ax[1, 1].axis('off')\n        ax[2, 0].set_title(\"Final Mask\")\n        ax[2, 0].imshow(mask, cmap='gray')\n        ax[2, 0].axis('off')\n        ax[2, 1].set_title(\"Apply Mask on Original\")\n        ax[2, 1].imshow(mask*img, cmap='gray')\n        ax[2, 1].axis('off')\n        \n        plt.show()\n    return mask*img","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"make_lungmask(image1[14], True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_rows_cols(size):\n    cols = 6\n    rows = size // cols\n    if (int(size%cols) != 0):\n        rows = rows+1\n    return rows,cols","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_stack(stack, start_with=10, show_every=3):\n    size = (len(stack) - (start_with - 1))//show_every\n    rows, cols = get_rows_cols(size)\n    plt.tight_layout()\n    fig,ax = plt.subplots(rows,cols,figsize=[12,12])\n    for i in range(size-1):\n        ind = start_with + i*show_every\n        ax[int(i/cols),int(i % cols)].set_title('slice %d' % ind)\n        ax[int(i/cols),int(i % cols)].imshow(stack[ind],cmap='gray')\n        ax[int(i/cols),int(i % cols)].axis('off')\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_stack(patient1_hu_scans, start_with=0, show_every=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"masked_lung = []\n\nfor img in image1:\n    masked_lung.append(make_lungmask(img))\n    \nplot_stack(masked_lung, start_with=0, show_every=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Code for running processing on the whole data all together. Outcome of the code will \n\n# be .npz file for all patients. .npz files can be loaded using np.load() function for further use.\n# If you want to store images in .png files remove the comments from below code and comment out code mentioned below.\n'''\npath = \"./segmented-images\"\nif not shutil.os.path.isdir(path):\n    shutil.os.mkdir(path)\n\npatients = train.Patient.unique()[0:10]\nfor patient in patients:\n    #if not shutil.os.path.isdir(path + \"/\" + patient):\n    #    shutil.os.mkdir(path + \"/\" + patient)\n    scans = load_slices(f'{TRAIN_IMG_PATH}/{patient}')\n    hu_imgs = get_pixels_hu(scans)\n    rescaled_images, spacing = resample(hu_imgs, scans,[1,1,1])\n\n    masked_lung = []\n    for img_number in range(len(rescaled_images)):\n        window_img = apply_window(rescaled_images[img_number], -600, 1200)\n        masked_img = make_lungmask(window_img)\n        masked_lung.append(masked_img)\n        #cv2.imwrite(f'{path}/{patient}/{img_number + 1}.png', masked_img)\n    # Comment the below line if images required to store in .png format.\n    np.savez(f'{path}/{patient}',masked_lung)\n    #plot_stack(masked_lung, start_with=0, show_every=1)\n'''","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"References:\n1. [Medical Coordinate System](https://theaisummer.com/medical-image-coordinates/)\n2. [Managing DICOM Images](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3354356/)\n3. [All about radiology - Radiopedia.org](https://radiopaedia.org/)\n4. [Domain expert's insight](https://www.kaggle.com/c/osic-pulmonary-fibrosis-progression/discussion/165727)\n5. [Intrinsic dependencies of CT radiomic features on voxel size and number of gray levels](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5462462/)\n6. [Hounsfield Scale](https://www.sciencedirect.com/topics/medicine-and-dentistry/hounsfield-scale)\n7. [Determining effective window width and center using different windowing techniques](https://www.researchgate.net/publication/272179277_Determining_effective_window_width_and_center_using_different_windowing_techniques_for_radio_therapy_images)\n8. [Resampling](https://www.dl-c.com/Temp/downloads/Whitepapers/Resampling.pdf)\n9. [Marching cubes](http://www.cs.carleton.edu/cs_comps/0405/shape/marching_cubes.html)\n10. [Image Segmentation](https://www.sciencedirect.com/topics/engineering/medical-image-segmentation)\n11. [Morphological Filtering](https://scikit-image.org/docs/dev/auto_examples/applications/plot_morphology.html)\n12. [DICOM Processing Segmentation Visualization in Python](https://www.raddq.com/dicom-processing-segmentation-visualization-in-python/)\n13. [Data Science Bowl 2017 - Preprocessing Tutorial by Guido Zuidhof](https://www.kaggle.com/gzuidhof/full-preprocessing-tutorial)\n14. [Pulmonary DICOM Preprocessing](https://www.kaggle.com/allunia/pulmonary-dicom-preprocessing)\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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}