{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.7.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":22307,"databundleVersionId":1502524,"sourceType":"competition"}],"dockerImageVersionId":30008,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"In this notebook, I will be going through a step-by-step guide on how to apply statistical clustering methods, computer graphics algorithms, and image processing techniques to medical images to help understand and visualize the data in both 2D and 3D.\n\n\nNOTE: In order to use plotly you will need to: \n1. Install all necessary packages/extensions following these [instructions](https://plotly.com/python/getting-started/). \n2. Sign up for a free account and activate your key following [these instructions](https://plotly.com/python/getting-started/) (only read the section ‘Initialization for Online Plotting’).","metadata":{}},{"cell_type":"markdown","source":"## Step 1: Install Necessary Packages\n\nNOTE: In order to use plotly you will need to: \n1. Install all necessary packages/extensions following these [instructions](https://plotly.com/python/getting-started/). \n2. Sign up for a free account and activate your key following [these instructions](https://plotly.com/python/getting-started/) (only read the section ‘Initialization for Online Plotting’).","metadata":{}},{"cell_type":"code","source":"!pip install  chart_studio\n","metadata":{"trusted":true,"_kg_hide-output":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Now let's load the necessary packages for the whole notebook. You might need to ","metadata":{}},{"cell_type":"code","source":"# common packages \nimport numpy as np \nimport os\nimport copy\nfrom math import *\nimport matplotlib.pyplot as plt\nfrom functools import reduce\nfrom glob import glob\n\n# reading in dicom files\nimport pydicom\n\n# skimage image processing packages\nfrom skimage import measure, morphology\nfrom skimage.morphology import ball, binary_closing\nfrom skimage.measure import label, regionprops\n\n# scipy linear algebra functions \nfrom scipy.linalg import norm\nimport scipy.ndimage\n\n# ipywidgets for some interactive plots\nfrom ipywidgets.widgets import * \nimport ipywidgets as widgets\n\n# plotly 3D interactive graphs \nimport plotly\nfrom plotly.graph_objs import *\nimport chart_studio\nchart_studio.tools.set_credentials_file(username='redwankarimsony', api_key='aEbXWsleQv7PJrAtOkBk')\n# set plotly credentials here \n# this allows you to send results to your account plotly.tools.set_credentials_file(username=your_username, api_key=your_key)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Step 2: Loading DICOM Data\nLet's set the paths of the dicom files and then we will be able to have a look at them. ","metadata":{}},{"cell_type":"code","source":"patient_id = '6897fa9de148'\npatient_folder = f'../input/rsna-str-pulmonary-embolism-detection/train/{patient_id}/'\ndata_paths = glob(patient_folder + '/*/*.dcm')\n\n# Print out the first 5 file names to verify we're in the right folder.\nprint (f'Total of {len(data_paths)} DICOM images.\\nFirst 5 filenames:' )\ndata_paths[:5]","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def load_scan(paths):\n    slices = [pydicom.read_file(path ) for path in paths]\n    slices.sort(key = lambda x: int(x.InstanceNumber), reverse = True)\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\ndef get_pixels_hu(scans):\n    image = np.stack([s.pixel_array for s in scans])\n    image = image.astype(np.int16)\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    intercept = scans[0].RescaleIntercept\n    slope = scans[0].RescaleSlope\n    \n    if slope != 1:\n        image = slope * image.astype(np.float64)\n        image = image.astype(np.int16)\n        \n    image += np.int16(intercept)\n    \n    return np.array(image, dtype=np.int16)","metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Run the following code to extract DICOM pixels for each slice location and display a single slice:","metadata":{}},{"cell_type":"code","source":"# set path and load files \npatient_dicom = load_scan(data_paths)\npatient_pixels = get_pixels_hu(patient_dicom)\n#sanity check\nplt.imshow(patient_pixels[80], cmap=plt.cm.bone)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Step 3: Image Processing\nLets use some thresholding and morphological operations to segment just the lung from the chest:","metadata":{}},{"cell_type":"code","source":"def largest_label_volume(im, bg=-1):\n    vals, counts = np.unique(im, return_counts=True)\n    counts = counts[vals != bg]\n    vals = vals[vals != bg]\n    if len(counts) > 0:\n        return vals[np.argmax(counts)]\n    else:\n        return None\n    \ndef segment_lung_mask(image, fill_lung_structures=True):\n    # not actually binary, but 1 and 2. \n    # 0 is treated as background, which we do not want\n    binary_image = np.array(image >= -700, dtype=np.int8)+1\n    labels = measure.label(binary_image)\n \n    # Pick the pixel in the very corner to determine which label is air.\n    # Improvement: Pick multiple background labels from around the  patient\n    # More resistant to “trays” on which the patient lays cutting the air around the person in half\n    background_label = labels[0,0,0]\n \n    # Fill the air around the person\n    binary_image[background_label == labels] = 2\n \n    # Method of filling the lung structures (that is superior to \n    # something like morphological closing)\n    if fill_lung_structures:\n        # For every slice we determine the largest solid structure\n        for i, axial_slice in enumerate(binary_image):\n            axial_slice = axial_slice - 1\n            labeling = measure.label(axial_slice)\n            l_max = largest_label_volume(labeling, bg=0)\n \n            if l_max is not None: #This slice contains some lung\n                binary_image[i][labeling != l_max] = 1\n    binary_image -= 1 #Make the image actual binary\n    binary_image = 1-binary_image # Invert it, lungs are now 1\n \n    # Remove other air pockets inside body\n    labels = measure.label(binary_image, background=0)\n    l_max = largest_label_volume(labels, bg=0)\n    if l_max is not None: # There are air pockets\n        binary_image[labels != l_max] = 0\n \n    return binary_image","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"By running the code below, you are using skimage functions from above to create a mask that covers the lung. We will use both `fill_lung_structures=True` and `fill_lung_structures=False`, to isolate the lung and the internal structures. Let’s run the code below, and display an example of isolating the lung from the chest:","metadata":{}},{"cell_type":"code","source":"# get masks \nsegmented_lungs = segment_lung_mask(patient_pixels, fill_lung_structures=False)\nsegmented_lungs_fill = segment_lung_mask(patient_pixels, fill_lung_structures=True)\ninternal_structures = segmented_lungs_fill - segmented_lungs\n\n# isolate lung from chest\ncopied_pixels = copy.deepcopy(patient_pixels)\nfor i, mask in enumerate(segmented_lungs_fill): \n    get_high_vals = mask == 0\n    copied_pixels[i][get_high_vals] = 0\nseg_lung_pixels = copied_pixels\n# sanity check\nf, ax = plt.subplots(1,2, figsize=(10,6))\nax[0].imshow(patient_pixels[80], cmap=plt.cm.bone)\nax[0].axis(False)\nax[0].set_title('Original')\nax[1].imshow(seg_lung_pixels[80], cmap=plt.cm.bone)\nax[1].axis(False)\nax[1].set_title('Segmented')\nplt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"If it looks like this, then perfect! You just isolated the lung from the rest of the scan. Don’t worry about the masks yet, I will show you some visualizations later on.","metadata":{}},{"cell_type":"markdown","source":"## Step 3: 2D Visualizations Techniques\n**Non Interactive:**\nWhen visualizing data, I find very beneficial to visualize each process of your script. This will not only help you understand each step of your code, but it makes for a very nice, clean presentation.\nThe code below will, in essence, visually convey a story on how you segmented your data: \n\n&#9632; original image <br>\n&#9632; the binary mask that covers the lung <br>\n&#9632; highlighting internal structures of the lung using one of the masks <br>\n&#9632; extracting internal structures using GK clustering","metadata":{}},{"cell_type":"code","source":"f, ax = plt.subplots(2,2, figsize = (10,10))\n\n# pick random slice \nslice_id = 80\n\nax[0,0].imshow(patient_pixels[slice_id], cmap=plt.cm.bone)\nax[0,0].set_title('Original Dicom')\nax[0,0].axis(False)\n\n\nax[0,1].imshow(segmented_lungs_fill[slice_id], cmap=plt.cm.bone)\nax[0,1].set_title('Lung Mask')\nax[0,1].axis(False)\n\nax[1,0].imshow(seg_lung_pixels[slice_id], cmap=plt.cm.bone)\nax[1,0].set_title('Segmented Lung')\nax[1,0].axis(False)\n\nax[1,1].imshow(seg_lung_pixels[slice_id], cmap=plt.cm.bone)\nax[1,1].imshow(internal_structures[slice_id], cmap='jet', alpha=0.7)\nax[1,1].set_title('Segmentation with \\nInternal Structure')\nax[1,1].axis(False)\n\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Interactive(1):\nThere are really only a few ways to display multiple images on jupyter notebook — manually plot each image one by one, or make a plot with multiple columns and rows to display in one figure. Unfortunately, this is not helpful for those that want to scan through each image and get a better understanding of the data. With the code below you create an interactive slide bar that lets you scroll through the images <font color=red>(But the catch is that you have to run the code in **interactive window**. So if you want to use the slider window to browse through the slices, fork the code and run it manually in interactive mode)</font>:","metadata":{}},{"cell_type":"code","source":"# slide through dicom images using a slide bar \nplt.figure(1)\ndef dicom_animation(x):\n    plt.imshow(patient_pixels[x], cmap = plt.cm.gray)\n    return x\ninteract(dicom_animation, x=(0, len(patient_pixels)-1))","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Interactive(2):\nAnother way to visualize the CT-Angiograms in a bit lively fashion is to use gif images. Basically GIFs are a series of images shown at an preconfigured interval automatically. <font color=blue>However, the upper side of GIF is that, it stores images in lossless compression. </blue>. So it is highly efficent and accurate for storing CT scan slides. ","metadata":{"trusted":true}},{"cell_type":"code","source":"import imageio\nfrom IPython import display\nprint('Original Image Slices before processing')\nimageio.mimsave(f'./{patient_id}.gif', patient_pixels, duration=0.1)\ndisplay.Image(f'./{patient_id}.gif', format='png')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print('Lung Segmentation Mask')\nimageio.mimsave(f'./{patient_id}.gif', segmented_lungs_fill, duration=0.1)\ndisplay.Image(f'./{patient_id}.gif', format='png')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print('Segmented Part of Lung Tissue')\nimageio.mimsave(f'./{patient_id}.gif', seg_lung_pixels, duration=0.1)\ndisplay.Image(f'./{patient_id}.gif', format='png')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"However, among the previous 3 GIFs, one of the most important is the Lung Segmentation Mask. You can see there are several images where we can see that the mask only selects the tissue portion of the lung but it is not considering the air vessels thorugh the lungs i.e. Bronchioles.\n\n<br>\n<font color=blue>(**Bronchioles** are air passages inside the lungs that branch off like tree limbs from the bronchi—the two main air passages into which air flows from the trachea (windpipe) after being inhaled through the nose or mouth. The bronchioles deliver air to tiny sacs called alveoli where oxygen and carbon dioxide are exchanged. They are vulnerable to conditions like asthma, bronchiolitis, cystic fibrosis, and emphysema that can cause constriction and/or obstruction of the airways.) </font>\n\nWe need to close those gaps so that we can segment the whole lung portion. With this view in mind, we can run closing operation so that it will fill up those portions. ","metadata":{}},{"cell_type":"code","source":"from skimage.morphology import opening, closing\nfrom skimage.morphology import disk\n\ndef plot_comparison(original, filtered, filter_name):\n\n    fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(8, 4), sharex=True,\n                                   sharey=True)\n    ax1.imshow(original, cmap=plt.cm.gray)\n    ax1.set_title('original')\n    ax1.axis('off')\n    ax2.imshow(filtered, cmap=plt.cm.gray)\n    ax2.set_title(filter_name)\n    ax2.axis('off')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"original = segmented_lungs_fill[65]\n\nrows = 4\ncols = 4\nf, ax = plt.subplots(rows, cols, figsize = (15,12))\n\nfor i in range(rows*cols):\n    if i==0:\n        ax[0,0].imshow(original, cmap = plt.cm.gray)\n        ax[0,0].set_title('Original')\n        ax[0,0].axis(False)\n    else:\n        closed = closing(original, disk(i))\n        ax[int(i/rows),int(i % rows)].set_title(f'closed disk({i})')\n        ax[int(i/rows),int(i % rows)].imshow(closed, cmap = plt.cm.gray)\n        ax[int(i/rows),int(i % rows)].axis('off')\nplt.show()   \n    ","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Therefore we can select the desired filters from here and just multiply that with the original image to get lung segmentation from these. Before selecting any particular filter, let's see the segmentation quality of the filters and then we can select the desired filter. ","metadata":{"trusted":true}},{"cell_type":"code","source":"original_image = patient_pixels[65]\noriginal = segmented_lungs_fill[65]\nf, ax = plt.subplots(rows, cols, figsize = (15,15))\n\nfor i in range(rows*cols):\n    if i==0:\n        ax[0,0].imshow(original_image, cmap = plt.cm.gray)\n        ax[0,0].set_title('Original')\n        ax[0,0].axis(False)\n    else:\n        closed = closing(original, disk(i))\n        ax[int(i/rows),int(i % rows)].set_title(f'closed disk({i})')\n        ax[int(i/rows),int(i % rows)].imshow(original_image * closed, cmap = plt.cm.gray)\n        ax[int(i/rows),int(i % rows)].axis('off')\nplt.show()   ","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Now it is super clear that filter size 15 is a good filter for lung segmentation. <font color=red></font>","metadata":{}},{"cell_type":"markdown","source":"![](https://www.clipartmax.com/png/middle/265-2655834_work-in-progress-icon.png)\n\n\n\n\n\n### In the meantime, check out my other ongoing works in this same competition: \n💥 [RSNA-STR Pulmonary Embolism [Dummy Sub]](https://www.kaggle.com/redwankarimsony/rsna-str-pulmonary-embolism-dummy-sub)<br>\n💥 [CT-Scans, DICOM files, Windowing Explained](https://www.kaggle.com/redwankarimsony/ct-scans-dicom-files-windowing-explained)<br>\n💥 [RSNA-STR-PE [Gradient & Sigmoid Windowing]](https://www.kaggle.com/redwankarimsony/rsna-str-pe-gradient-sigmoid-windowing)<br>\n💥 [RSNA-STR [✔️3D Stacking ✔️3D Plot ✔️Segmentation]](https://www.kaggle.com/redwankarimsony/rsna-str-3d-stacking-3d-plot-segmentation/edit/run/42517982)<br>\n💥 [RSNA-STR [DICOM 👉 GIF 👉 npy]](https://www.kaggle.com/redwankarimsony/rsna-str-dicom-gif-npy)<br>\n💥 [RSNA-STR Pulmonary Embolism [EDA]](https://www.kaggle.com/redwankarimsony/rsna-str-pulmonary-embolism-eda)<br>","metadata":{}},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}