{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os\nimport pydicom\n\nimport matplotlib.pyplot as plt\nimport matplotlib.animation as animation\nfrom matplotlib.widgets import Slider\n\nfrom IPython.display import HTML","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We take one of the patients and load all its dicom images as numpy arrays. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"?np.sort","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_id = \"ID00035637202182204917484\"\n\ndicom_path = \"/kaggle/input/osic-pulmonary-fibrosis-progression/train\"\n\nfiles = np.array([f.replace(\".dcm\",\"\") for f in os.listdir(f\"{dicom_path}/{patient_id}/\")])\nfiles = -np.sort(-files.astype(\"int\"))\ndicoms = [f\"{dicom_path}/{patient_id}/{f}.dcm\" for f in files]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In these Dicoms the slope and intercept seems to be 1 and 0 respectively so there is no need for the transformation, but putting the whole code just in case and also for recyclability!","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"images = []\nfor dcm in dicoms:\n    tmp = pydicom.dcmread(dcm)\n    slope = tmp.RescaleSlope\n    intercept = tmp.RescaleIntercept\n    final = tmp.pixel_array*slope + intercept\n    images.append(final)\n    \nimages = np.array(images) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## From the top","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Here is how to create an animation out of an array of images.","execution_count":null},{"metadata":{"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure()\n\nims = []\nfor image in range(0,images.shape[0],10):\n    im = plt.imshow(images[image,:,:], \n                    animated=True, cmap=plt.cm.bone)\n    plt.axis(\"off\")\n    ims.append([im])\n\nani = animation.ArtistAnimation(fig, ims, interval=100, blit=False,\n                                repeat_delay=1000)\n\nplt.close()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Two ways of animating the images! \nThe JavaScript one is way more iteractive while the second html5 one is basically a video.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"HTML(ani.to_jshtml())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"HTML(ani.to_html5_video())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Frontwise","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure()\n\nims = []\nfor image in range(0,images.shape[1],5):\n    im = plt.imshow(images[:,image,:], animated=True, cmap=plt.cm.bone)\n    plt.axis(\"off\")\n    ims.append([im])\n\nani = animation.ArtistAnimation(fig, ims, interval=100, blit=False,\n                                repeat_delay=1000)\n\nplt.close()\n\nHTML(ani.to_jshtml())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Sidewise","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure()\n\nims = []\nfor image in range(0,images.shape[2],5):\n    im = plt.imshow(images[:,:,image], animated=True, cmap=plt.cm.bone)\n    plt.axis(\"off\")\n    ims.append([im])\n\nani = animation.ArtistAnimation(fig, ims, interval=100, blit=False,\n                                repeat_delay=1000)\n\nplt.close()\n\nHTML(ani.to_jshtml())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From the internet I found that in a CT Scan, the HU (Hounsfield's units) is calcluated based on the linear attenuation coefficient $\\mu$, with the following fomula:\n\n$$\nHU = 1000 \\cdot \\frac{\\mu_{X} - \\mu_{water}}{\\mu_{air} - \\mu_{water}}\n$$\n\nAnd the most common values are:\n\n|Substance       |  \t HU      |\n|----------------|---------------|\n|Air \t         |    -1000      |\n|Lung \t         |    -700       |\n|Fat \t         |    -84        |\n|Water \t         |     0         |\n|CSF \t         |     15        |\n|Blood \t         | +30 to +45    |\n|Muscle          |\t+40          |\n|Soft Tissue     | \t+100 to +300 |\n|Cancellous Bone | \t+700         |\n|Dense Bone \t |  +3000        |\n\nWhile -2048 would indicate a missing value (since the scan has a circular shape while our plot is rectangular.\n\nHere is the histogram of the CT scan of the Patient in consideration.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.hist(np.array(images).reshape(-1,), bins=50)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Unfortunately this last part won't work in the kaggle view mode, but on edit mode of the notebook. \nBut truth to be told the the java script version above can be seen as a interactive mode, where you can go frame by frame anche see all the slices.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from ipywidgets import interact\nimport ipywidgets as widgets","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%matplotlib notebook\n\nfig = plt.figure(figsize=(5,5))\n\nimg_plot = plt.imshow(images[0], cmap=\"Greys\")\nplt.axis(\"off\")\n\n@interact(slice = widgets.IntSlider(min=0, max=len(images), step=1, value=0))\ndef update(slice):\n    global img_plot\n    img_plot.set_data(images[int(slice)])\n    plt.draw()","execution_count":null,"outputs":[]},{"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}