{"cells":[{"metadata":{},"cell_type":"markdown","source":"This notebook was created using [this notebook.](https://www.kaggle.com/gzuidhof/full-preprocessing-tutorial/data?)"},{"metadata":{"trusted":true},"cell_type":"code","source":"import plotly.express as px\nimport plotly.graph_objects as go\nfrom plotly.subplots import make_subplots\nfrom plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\nfrom plotly import figure_factory as FF\n\nimport scipy.ndimage\nfrom skimage import measure, morphology\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","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\n\nfrom IPython.display import HTML","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Some constants \npath_input = '../input/osic-pulmonary-fibrosis-progression/train/'\npatients = os.listdir(path_input)\npatients.sort()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load files\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"},{"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, sourced from Wikipedia.\n\n|Substance|     HU    |\n|:-------:|:---------:|\n|Air      |   −1000   |\n|Lung     |−600 to -400 |\n|Fat\t  |−100 to −50|\n|Water\t  |0          |\n|Blood\t  |+30 to +70 |\n|Muscle\t  |+10 to +40 |\n|Liver\t  |+40 to +60 |\n|Bone\t  |+700 to +3000|\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"def 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","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"first_patient = load_scan(path_input + patients[1])\n\nprint(type(first_patient[8].pixel_array))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_pixels_hu(scans):\n    image = np.stack([s.pixel_array for s in scans])\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 1\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)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"imgs = get_pixels_hu(first_patient)\nplt.hist(imgs.flatten(), bins=50, color='c')\nplt.xlabel(\"Hounsfield Units (HU)\")\nplt.ylabel(\"Frequency\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Standardize the pixel values\nfrom sklearn.cluster import KMeans\nfrom scipy import ndimage\n\ndef make_lungmask(img, display=False):\n    row_size= img.shape[0]\n    col_size = img.shape[1]\n    \n    \n    mean = np.mean(img)\n    std = np.std(img)\n    img = img-mean\n    img = img/std\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    thresh_img = np.where(img<threshold,1.0,0.0)  # threshold the image\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/6 and B[2]<col_size/6*5):\n            good_labels.append(prop.label)\n    mask = np.ndarray([row_size,col_size],dtype=np.int8)\n    mask[:] = 0\n\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    #\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(first_patient[50].pixel_array, display=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.animation as animation\n\nfig = plt.figure()\n\nims = []\nfor scan in first_patient:\n    im = plt.imshow(scan.pixel_array, animated=True, cmap=\"Greys\")\n    plt.axis(\"off\")\n    ims.append([im])\n\nani = animation.ArtistAnimation(fig, ims, interval=100, blit=False,\n                                repeat_delay=1000)\n\nHTML(ani.to_jshtml())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(5,5))\n\nims = []\nfor scan in first_patient:\n    im = ax.imshow(make_lungmask(scan.pixel_array), animated=True, cmap=\"gray\")\n    ax.axis(\"off\")\n    ims.append([im])\n\nani = animation.ArtistAnimation(fig, ims, interval=150, blit=False,\n                                repeat_delay=1000)\n\nHTML(ani.to_jshtml())","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}