{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Overview\n\nThe purpose of this notebook is to improve lung segmentation shown here [here](https://www.kaggle.com/andradaolteanu/pulmonary-fibrosis-competition-eda-dicom-prep) and demostrate tissue segmentation.\n"},{"metadata":{},"cell_type":"markdown","source":"# Load Libraries"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# Regular Imports\nimport numpy as np\nimport pandas as pd\nimport os\nimport random\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm\nfrom PIL import Image\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.model_selection import GroupKFold, KFold, StratifiedKFold\nimport seaborn as sns\nimport math\nimport cv2\nimport pydicom\nimport os\nfrom tqdm import tqdm\nimport glob\nimport pickle as pkl\nimport matplotlib.image as mpimg\nfrom tabulate import tabulate\nimport missingno as msno \nfrom IPython.display import display_html\nfrom PIL import Image\nimport gc\nfrom skimage.transform import resize\nimport copy\nimport re\nfrom scipy.stats import pearsonr\n%matplotlib inline\n\n# Segmentation\nimport 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 *\ninit_notebook_mode(connected=True) \n\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\npd.set_option(\"display.max_columns\", 100)\n\ncustom_colors = ['#74a09e','#86c1b2','#98e2c6','#f3c969','#f2a553', '#d96548', '#c14953']\nsns.palplot(sns.color_palette(custom_colors))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"N_ROWS = 512\nN_COLS = 512\npath = '../input/osic-pulmonary-fibrosis-progression/'\nBATCH_SIZE=128","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Helper Functions"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"# Areas with the same number of pixels on the edges are not required. Crop it.\ndef crop_image(img: np.ndarray):\n    edge_pixel_value = img[0, 0]\n    mask = img != edge_pixel_value\n    return img[np.ix_(mask.any(1),mask.any(0))]\n\n# Load images, crop thick borders(if any) and resize\ndef load_image(path):\n    dataset = pydicom.dcmread(path)\n    img = dataset.pixel_array\n    img = crop_image(img)\n#     img = cv2.resize(img, (512,512))\n    return img\n\n# Get Nth percentile image\ndef get_img(perc, patient_id, data):\n    \n    l = glob.glob(path+'{0}/{1}/*.dcm'.format(data, patient_id))\n    img_ids = []\n    for x in l:\n        y = x.split('/')[-1]\n        z = int(y.split('.')[0])\n        img_ids.append(z)\n\n    img_ids.sort()\n    \n    return img_ids[math.ceil(perc*(len(img_ids)))-1]\n\n# Get num of slices bw two percentiles\ndef num_img_bw_perc(p1, p2, patient_id, data):\n    \n    l = glob.glob(path+'{0}/{1}/*.dcm'.format(data, patient_id))\n    img_ids = []\n    for x in l:\n        y = x.split('/')[-1]\n        z = int(y.split('.')[0])\n        img_ids.append(z)\n\n    img_ids.sort()\n    \n    return len(img_ids[math.ceil(p1*(len(img_ids)))-1:math.ceil(p2*(len(img_ids)))])-1\n    \n\n# Get number of images per patient\ndef get_num_images(patient_id, data):\n    \n    return len(glob.glob(path+'{0}/{1}/*.dcm'.format(data, patient_id)))\n\n# Get the lung area in the image slice\ndef lung_seg_pixel_ratio(img_array):\n    \n    c = 0\n    for i in range(img_array.shape[0]):\n        for j in range(img_array.shape[1]):\n            if img_array[i][j] != 0:\n                c+=1\n    \n    return c, round(c/(img_array.shape[0]*img_array.shape[1]),4)\n\n# Get dicom meta data\ndef get_dicom_meta(path):\n    \n    '''Get information from the .dcm files.\n    path: complete path to the .dcm file'''\n\n    image_data = pydicom.read_file(path)\n\n    # Dictionary to store the information from the image\n    observation_data = {\n        \"SliceThickness\" : float(image_data.SliceThickness),\n        \"PixelSpacing\" : float(image_data.PixelSpacing[0]),\n    }\n    \n    return observation_data\n\n# Get tissue mask\n# To extract the tissues from the segmented lung all we need to do is get rid of the border parts from the segmented lung \n# Grey pixels present within the border of the lung is assumed to be tissue.\n# Inorder to get rid of the border pixels of the lung we slightly perturb the segmented lung to the right, left, top and bottom\n# The intersection of all the perturbed images gets rid of the border lung pixels \n# This resultant image serves as the mask for the tissue segmentation\ndef tissue_mask(img, mask, shift_perc):\n    \n    r_dim, c_dim = img.shape[0], img.shape[1]\n    \n    # Move the image by shift_perc to the left\n    del_left_cols = int(shift_perc*c_dim)\n    \n    mask1, mask2 = mask.copy(), np.zeros((r_dim, c_dim)).astype(int)\n    mask1 = mask1[:,del_left_cols:]\n    mask2[:,:c_dim-del_left_cols] = mask1\n    \n    # Move the image by shift_perc to the right\n    del_right_cols = int(shift_perc*c_dim)\n    \n    mask3, mask4 = mask.copy(), np.zeros((r_dim, c_dim)).astype(int)\n    mask3 = mask3[:,:c_dim-del_right_cols]\n    mask4[:,del_right_cols:] = mask3\n    \n    # Move the image by shift_perc to the top\n    del_top_rows = int(shift_perc*c_dim)\n    \n    mask5, mask6 = mask.copy(), np.zeros((r_dim, c_dim)).astype(int)\n    mask5 = mask5[del_top_rows:,:]\n    mask6[:r_dim-del_top_rows,:] = mask5\n    \n    # Move the image by shift_perc to the bottom\n    del_bottom_rows = int(shift_perc*r_dim)\n    \n    mask7, mask8 = mask.copy(), np.zeros((r_dim, c_dim)).astype(int)\n    mask7 = mask7[:r_dim-del_bottom_rows,:]\n    mask8[del_bottom_rows:,:] = mask7\n\n    #Obtain the final mask\n    final_mask = ((mask2==1) & (mask4==1) & (mask6==1) & (mask8==1)).astype(int)\n    \n    return final_mask\n    \n# Get tissue features\ndef tissue_features(tissue_mask, img, thresh = 0.35):\n    \n    final_img = tissue_mask*img\n\n    checker = np.zeros((final_img.shape[0], final_img.shape[1]))\n    counter, other_counter = 0, 0\n    for i in range(final_img.shape[0]):\n        for j in range(final_img.shape[1]):\n            if final_img[i][j]>=thresh:\n                checker[i][j] = 1\n                counter+=1\n            else:\n                checker[i][j] = 0\n                other_counter+=1\n                \n    tissue_by_total = counter/(final_img.shape[0]**2)\n    tissue_by_lung = counter/((tissue_mask==1).sum())\n   \n    return counter, tissue_by_total, tissue_by_lung, checker\n\n# https://www.raddq.com/dicom-processing-segmentation-visualization-in-python/\ndef make_lungmask(img, display=False):\n    \n    img = img.astype(float)\n    \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    \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    \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    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*0.9) and (B[3]-B[1]<col_size*0.9) and (B[2]-B[0]>row_size*0.20) \n            and (B[3]-B[1]>col_size*0.10) and (B[0]>row_size*0.03) and (B[2]<row_size*0.97) \n            and (B[1]>col_size*0.03) and (B[3]<col_size*0.97)):\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    \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    # Compute Lung Area in the slice\n    lung_pixels, slice_lung_area = lung_seg_pixel_ratio(mask)\n    \n    # Tissue Mask\n    t_mask = tissue_mask(img, mask, shift_perc = 0.02)\n    \n    # Extract tissue features\n    num_t_pixels, tissue_by_total, tissue_by_lung, checker = tissue_features(t_mask, img, thresh = 0.35)\n    \n    if (display):\n        fig, ax = plt.subplots(4, 2, figsize=[18, 18])\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        ax[3, 0].set_title(\"Inner Lung Mask\")\n        ax[3, 0].imshow(t_mask, cmap='gray')\n        ax[3, 0].axis('off')\n        ax[3, 1].set_title(\"Segmented Tissue\")\n        ax[3, 1].imshow(checker, cmap='gray')\n        ax[3, 1].axis('off')\n        \n        plt.show()\n        \n    return lung_pixels, slice_lung_area, num_t_pixels, tissue_by_total, tissue_by_lung","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"make_lungmask(load_image(path+'train/ID00007637202177411956430/15.dcm'), display = True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"make_lungmask(load_image(path+'train/ID00048637202185016727717/11.dcm'), display = True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"make_lungmask(load_image(path+'train/ID00072637202198161894406/12.dcm'), display = True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# What next? - Feature Extraction"},{"metadata":{},"cell_type":"markdown","source":"We choose 11 images from the middle 30-60 percentile for each CT scan. Using thise images we extract the following features:\n- Lung volume\n- Average of number of tissue pixels present across all selected images for a patient\n- Average tissue area\n- Average tissue volume\n- Average of (number of tissue pixels)/(Total pixels on image)\n- Average of (tissue area)/(lung area)\n\nFeature extraction, data processing and final submission notebook can be found [here](https://www.kaggle.com/abhishekgbhat/quantreg-linear-decay-efficientnet-b1-su)"},{"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}