{"cells":[{"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)\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#for 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":{},"cell_type":"markdown","source":"## Libraries","execution_count":null},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"%matplotlib inline\n\nimport numpy as np\nimport pydicom\nimport os\nimport matplotlib.pyplot as plt\nfrom glob import glob\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nimport scipy.ndimage\nfrom skimage import morphology\nfrom skimage import measure\nfrom skimage.transform import resize\nfrom sklearn.cluster import KMeans\nfrom plotly import __version__\nfrom plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\nfrom plotly.tools import FigureFactory as FF\nfrom plotly.graph_objs import *\nimport plotly.express as px\ninit_notebook_mode(connected=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Data Loading","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"data_path = '/kaggle/input/osic-pulmonary-fibrosis-progression/train/'\ntrain_data = pd.read_csv('/kaggle/input/osic-pulmonary-fibrosis-progression/train.csv')\npatient_list = train_data.Patient.unique()\nprint(\"Number of Patients :\",len(patient_list) )","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Data Exploration","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"patient = pd.DataFrame()\npid = []\ncount = []\npath = []\nfor pat in patient_list: \n    g = glob(data_path +pat +'/*.dcm')\n    pid.append(pat)\n    path.append(g)\n    count.append(len(g))\npatient['pid'] = pid\npatient['scan_count'] = count\npatient['path'] = path\nfig = px.bar(patient, x = 'pid', y = 'scan_count')\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#      \n# Loop over the image files and store everything into a list.\n# \n\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\ndef 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":"# ID00007637202177411956430\npatient = load_scan(data_path+'ID00007637202177411956430/')\nimgs = get_pixels_hu(patient)\nout_path = '/kaggle/working/'\nid = 0\nnp.save(out_path + \"fullimages_%d.npy\" % (id), imgs)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(imgs[25], cmap=plt.cm.bone)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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\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","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get masks \nimport copy\nsegmented_lungs = segment_lung_mask(imgs, fill_lung_structures=False)\nsegmented_lungs_fill = segment_lung_mask(imgs,fill_lung_structures=True)\ninternal_structures = segmented_lungs_fill - segmented_lungs\n# isolate lung from chest\ncopied_pixels = copy.deepcopy(imgs)\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\nplt.figure(figsize = (20,15))\nplt.imshow(seg_lung_pixels[25], cmap=plt.cm.bone)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"file_used=out_path+\"fullimages_%d.npy\" % id\nimgs_to_process = np.load(file_used).astype(np.float64) \nplt.figure(figsize = (15,12))\nplt.hist(imgs_to_process.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":"id = 0\nimgs_to_process = np.load(out_path+'fullimages_{}.npy'.format(id))\n\ndef sample_stack(stack, rows=5, cols=5, start_with=1, show_every=1):\n    fig,ax = plt.subplots(rows,cols,figsize=[12,12])\n    for i in range(rows*cols):\n        ind = start_with + i*show_every\n        ax[int(i/rows),int(i % rows)].set_title('slice %d' % ind)\n        ax[int(i/rows),int(i % rows)].imshow(stack[ind],cmap='gray')\n        ax[int(i/rows),int(i % rows)].axis('off')\n    plt.show()\n\nsample_stack(imgs_to_process)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Resampling","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":{"trusted":true},"cell_type":"code","source":"id = 0\nimgs_to_process = np.load(out_path+'fullimages_{}.npy'.format(id))\ndef 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)\n    \n    return image, new_spacing\n\nprint(\"Shape before resampling\\t\", imgs_to_process.shape)\nimgs_after_resamp, spacing = resample(imgs_to_process, patient, [1,1,1])\nprint(\"Shape after resampling\\t\", imgs_after_resamp.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3D Plotting","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_mesh(image, threshold=-300, step_size=1):\n\n    print(\"Transposing surface\")\n    p = image.transpose(2,1,0)\n    \n    print(\"Calculating surface\")\n    verts, faces, norm, val = measure.marching_cubes_lewiner(p, threshold, step_size=step_size, allow_degenerate=True) \n    return verts, faces\n\ndef plotly_3d(verts, faces):\n    x,y,z = zip(*verts) \n    \n    print(\"Drawing\") \n    \n    # Make the colormap single color since the axes are positional not intensity. \n#    colormap=['rgb(255,105,180)','rgb(255,255,51)','rgb(0,191,255)']\n    colormap=['rgb(236, 236, 212)','rgb(236, 236, 212)']\n    \n    fig = FF.create_trisurf(x=x,\n                        y=y, \n                        z=z, \n                        plot_edges=False,\n                        colormap=colormap,\n                        simplices=faces,\n                        backgroundcolor='rgb(64, 64, 64)',\n                        title=\"Interactive Visualization\")\n    iplot(fig)\n\ndef plt_3d(verts, faces):\n    print(\"Drawing\") \n    x,y,z = zip(*verts) \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], linewidths=0.05, alpha=1)\n    face_color = [1, 1, 0.9]\n    mesh.set_facecolor(face_color)\n    ax.add_collection3d(mesh)\n\n    ax.set_xlim(0, max(x))\n    ax.set_ylim(0, max(y))\n    ax.set_zlim(0, max(z))\n    #ax.set_axis_bgcolor((0.7, 0.7, 0.7))\n    plt.show()\nv, f = make_mesh(imgs_after_resamp, 350)\nplt_3d(v, f)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Segmentation","execution_count":null},{"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    \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/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    #\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":"img = imgs_after_resamp[260]\nmake_lungmask(img, display=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"masked_lung = []\n\nfor img in imgs_after_resamp:\n    masked_lung.append(make_lungmask(img))\n\nsample_stack(masked_lung, show_every=10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Source blog : https://www.raddq.com/dicom-processing-segmentation-visualization-in-python/","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Processing with SimpleITK","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def sitk_show(img, title=None, margin=0.05, dpi=40 ):\n    nda = SimpleITK.GetArrayFromImage(img)\n    spacing = img.GetSpacing()\n    figsize = (1 + margin) * nda.shape[0] / dpi, (1 + margin) * nda.shape[1] / dpi\n    extent = (0, nda.shape[1]*spacing[1], nda.shape[0]*spacing[0], 0)\n    fig = plt.figure(figsize=figsize, dpi=dpi)\n    ax = fig.add_axes([margin, margin, 1 - 2*margin, 1 - 2*margin])\n\n    plt.set_cmap(\"gray\")\n    ax.imshow(nda,extent=extent,interpolation=None)\n    \n    if title:\n        plt.title(title)\n    \n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import SimpleITK","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"reader = SimpleITK.ImageSeriesReader()\nfilenamesDICOM = reader.GetGDCMSeriesFileNames(data_path+'ID00007637202177411956430/')\nreader.SetFileNames(filenamesDICOM)\nimgOriginal = reader.Execute()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"idxSlice = 25\nimgOriginal = imgOriginal[:,:,idxSlice]\nsitk_show(imgOriginal)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Image Smoothing","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"imgSmooth = SimpleITK.CurvatureFlow(image1=imgOriginal,\n                                    timeStep=0.125,\n                                    numberOfIterations=5)\nsitk_show(imgSmooth)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## If you like this kernel plz <font color='red'>UPVOTE</font>.\n## If you have suggestions to improve this kernel plz <font color='red'>COMMENT</font>.\n\n<font color='Blue'>************************ Notebook is under construction ************************</font>","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}