{"cells":[{"metadata":{},"cell_type":"markdown","source":"# OSIC: Pulmonary Fibrosis Progression\n\n**In this notebook, I will try to perform some basic operations (EDA, Preprocessing & Augmentation) on this dataset. As I am a beginner myself, I will try to explain the findings as much as possible. Please give your valuable opinions and suggestions in the comments.**\n","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Part 1: EDA","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"**At first, Let's import the modules.**","execution_count":null},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport plotly.express as px\nimport cv2\n\nfrom matplotlib import pyplot as plt\n\n# Set Color Palettes for the notebook (https://color.adobe.com/)\ncolors_nude = ['#FFE61A','#B2125F','#FF007B','#14B4CC','#099CB3']\nsns.palplot(sns.color_palette(colors_nude))\n\n# Set Style\nsns.set_style(\"whitegrid\")\nsns.despine(left=True, bottom=True)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Load the training and test data and look at the first 5 entries of the Train data**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data = pd.read_csv(\"../input/osic-pulmonary-fibrosis-progression/train.csv\")\ntest_data = pd.read_csv(\"../input/osic-pulmonary-fibrosis-progression/test.csv\")\n\nprint(train_data.head())\nprint()\nprint(\"The number of total records in the Train Data: \",len(train_data))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Now, let's check if there's any missing data**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"print(train_data.isna().sum())\nprint()\nprint('Here, we can see there is no missing data in any of the columns.')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Now, let's see the summary of the columns in the train data**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data.describe(percentiles=[.20,.40,.60,.80])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Let's see the number of unique patients in the dataset.**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"The number of unique patients: \",len(train_data[\"Patient\"].unique()))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Let's see a sunburst plot to get a general idea on the data**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = px.sunburst(data_frame=train_data,\n                  path=['Age', 'Sex', 'SmokingStatus'],\n                  color='Sex',\n                  maxdepth=-1,\n                  title='Sunburst Chart')\n\nfig.update_traces(textinfo='label+percent parent')\nfig.update_layout(margin=dict(t=0, l=0, r=0, b=0))\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Now, let's see some distribution plots.**","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"**VARIABLE: AGE**","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"**First, let's see the age distribution of the patients**","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"sns.distplot(train_data['Age'],rug=True)\n\nprint(\"This plot shows most of the patients are from 65 to 75 years old.\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Now let's see the most common ages among the patients.**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Occurance of landmark_id in decreasing order(Top categories)\ntemp = pd.DataFrame(train_data.Age.value_counts().head(10))\ntemp.reset_index(inplace=True)\ntemp.columns = ['Age','Number of Patients']\n\n# Plot the most frequent landmark_ids\nplt.figure(figsize = (9, 10))\nplt.title('The most common ages of the patients')\nsns.set_color_codes(\"deep\")\nsns.barplot(x=\"Age\", y=\"Number of Patients\", data=temp,\n            label=\"Count\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**VARIABLE: SEX**","execution_count":null},{"metadata":{"_kg_hide-input":false,"trusted":true},"cell_type":"code","source":"df= train_data\ndf= df.drop_duplicates(subset='Patient',keep='first') \n\n\ntemp = pd.DataFrame(df.Sex.value_counts().head(10))\ntemp.reset_index(inplace=True)\ntemp.columns = ['Sex','Number of Patients']\n\n\n# Plot the most frequent landmark_ids\nplt.figure(figsize = (9, 10))\nplt.title('Sex distribution of the patients')\nsns.set_color_codes(\"deep\")\nsns.barplot(x=\"Sex\", y=\"Number of Patients\", data=temp,\n            label=\"Count\")\nplt.show()\n\nprint(\"This plot shows the number of male patients are about 3.5 times higher than female patients.\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Age distribution of the patients with respect to their gender**","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"sns.catplot(x=\"Patient\", y=\"Age\", hue=\"Sex\", kind=\"swarm\", data=df);\nprint(\"From this plot we can see that the oldest patient is Male and the youngest patient is a Female\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**VARIABLE: SmokingStatus**","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"df= train_data\ndf= df.drop_duplicates(subset='Patient',keep='first') \n\n\ntemp = pd.DataFrame(df.SmokingStatus.value_counts().head(10))\ntemp.reset_index(inplace=True)\ntemp.columns = ['Smoking Status','Number of Patients']\n\n\n# Plot the most frequent landmark_ids\nplt.figure(figsize = (9, 10))\nplt.title('Smoking Status distribution of the patients')\nsns.set_color_codes(\"deep\")\nsns.barplot(x=\"Smoking Status\", y=\"Number of Patients\", data=temp,\n            label=\"Count\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Smoking condition of the patients with respect to their age**","execution_count":null},{"metadata":{"_kg_hide-input":false,"trusted":true},"cell_type":"code","source":"sns.catplot(x=\"Patient\", y=\"Age\", hue=\"SmokingStatus\", kind=\"swarm\", data=df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Smoking condition of the patients with respect to their age and sex**","execution_count":null},{"metadata":{"_kg_hide-input":false,"trusted":true},"cell_type":"code","source":"sns.catplot(x=\"SmokingStatus\", y=\"Age\", hue=\"Sex\", kind=\"swarm\", data=df);\nprint(\"This plot shows that most of the Male patients are Ex-smoker and most the Frmale patients have Never Smoked\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"VARIABLE: FVC & Percent","execution_count":null},{"metadata":{"_kg_hide-input":false,"trusted":true},"cell_type":"code","source":"# Figure\nf, (ax1, ax2) = plt.subplots(1, 2, figsize = (16, 6))\n\na = sns.distplot(train_data[\"FVC\"], ax=ax1, hist=True)\nb = sns.distplot(train_data[\"Percent\"], ax=ax2, hist=True)\n\na.set_title(\"FVC Distribution\", fontsize=16)\nb.set_title(\"Percent Distribution\", fontsize=16);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Smoking condition of the patients with respect to their respective FVC and Gender**","execution_count":null},{"metadata":{"_kg_hide-input":false,"trusted":true},"cell_type":"code","source":"sns.catplot(x=\"SmokingStatus\", y=\"FVC\", hue=\"Sex\", kind=\"box\", data=df);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Smoking condition of the patients with respect to their respective Percent and Gender**","execution_count":null},{"metadata":{"_kg_hide-input":false,"trusted":true},"cell_type":"code","source":"sns.catplot(x=\"SmokingStatus\", y=\"Percent\", hue=\"Sex\", kind=\"box\", data=df);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Part 2: DICOM DATA\n\n**Digital Imaging and Communications in Medicine (DICOM) is the standard for the communication and management of medical imaging information and related data.**","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"**Import the necessary modules**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import pydicom\nfrom ipywidgets.widgets import * \nimport ipywidgets as widgets\n\nimport re\nfrom PIL import Image\nfrom IPython.display import Image as show_gif\nimport scipy.misc\nimport matplotlib\nfrom skimage import measure\nfrom skimage import morphology\nfrom skimage.transform import resize\nfrom sklearn.cluster import KMeans","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Let's look at some of the images.**","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Some Utility Function for loading the images taken from [here](https://medium.com/@hengloose/a-comprehensive-starter-guide-to-visualizing-and-analyzing-dicom-images-in-python-7a8430fcb7ed)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_scan(path):\n    slices = [pydicom.dcmread(path + '/' + s) for s in               \n              os.listdir(path)]\n    slices = [s for s in slices if 'SliceLocation' in s]\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    for s in slices:\n        s.SliceThickness = slice_thickness\n    return slices\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)\n\ndef dicom_animation(x):\n    plt.figure(figsize = (15,5))\n    plt.imshow(patient_pixels[x],cmap='gray')\n    return x","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Let's see a patient's CT scans in an interactive manner.**","execution_count":null},{"metadata":{"_kg_hide-input":false,"trusted":true},"cell_type":"code","source":"path = '../input/osic-pulmonary-fibrosis-progression/train/ID00010637202177584971671'\npatient_dicom = load_scan(path)\npatient_pixels = get_pixels_hu(patient_dicom)\n\nprint('There are {} images in this scan.'.format(len(patient_pixels)))\nprint()\nprint('Pateint: ', path.split('/')[-1])\ninteract(dicom_animation, x=(0, len(patient_pixels)-1))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Now let's see the number of scans per patient**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"path= '../input/osic-pulmonary-fibrosis-progression/train/'\npatients = os.listdir(path)\nimage_counts = []\nfor p in patients:\n    number_of_images= len(os.listdir(path+p))\n    #print(p, number_of_images)\n    image_counts.append(number_of_images)\n\ndata={'Patients':patients,'Number of Scans':image_counts}\n\nplt.figure(figsize = (20, 8))\np=sns.barplot(x='Patients', y='Number of Scans',data=data,palette=\"Blues_d\")\nplt.xlabel('Patient ID', fontsize=14)\nplt.ylabel('Frequency', fontsize=14)\n\nplt.title(\"Number of CT Scans per Patient\", fontsize = 15)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Now let's extract the meta datas from the DICOM files**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"path='../input/osic-pulmonary-fibrosis-progression/train/ID00019637202178323708467/1.dcm'\nimage= pydicom.dcmread(path)\nprint(image)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Creating a GIF from a patients CT scans**\n\nRef: https://www.kaggle.com/andradaolteanu/pulmonary-fibrosis-competition-eda-dicom-prep#To-Be-Continued-","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"**Function for creating GIF from a patients CT scans**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def create_gif(number_of_CT = 87):\n    \"\"\"Picks a patient at random and creates a GIF with their CT scans.\"\"\"\n    \n    # Select one of the patients\n    # patient = \"ID00007637202177411956430\"\n    patient = patients[image_counts.index(number_of_CT)]\n    \n    print('Patient: ',patient)\n    \n    # === READ IN .dcm FILES ===\n    patient_dir = \"../input/osic-pulmonary-fibrosis-progression/train/\" + patient\n    datasets = []\n\n    # First Order the files in the dataset\n    files = []\n    for dcm in list(os.listdir(patient_dir)):\n        files.append(dcm) \n    files.sort(key=lambda f: int(re.sub('\\D', '', f)))\n\n    # Read in the Dataset from the Patient path\n    for dcm in files:\n        path = patient_dir + \"/\" + dcm\n        datasets.append(pydicom.dcmread(path))\n        \n        \n    # === SAVE AS .png ===\n    # Create directory to save the png files\n    if os.path.isdir(f\"png_{patient}\") == False:\n        os.mkdir(f\"png_{patient}\")\n\n    # Save images to PNG\n    for i in range(len(datasets)):\n        img = datasets[i].pixel_array\n        matplotlib.image.imsave(f'png_{patient}/img_{i}.png', img)\n        \n        \n    # === CREATE GIF ===\n    # First Order the files in the dataset (again)\n    files = []\n    for png in list(os.listdir(f\"../working/png_{patient}\")):\n        files.append(png) \n    files.sort(key=lambda f: int(re.sub('\\D', '', f)))\n\n    # Create the frames\n    frames = []\n\n    # Create frames\n    for file in files:\n    #     print(\"../working/png_images/\" + name)\n        new_frame = Image.open(f\"../working/png_{patient}/\" + file)\n        frames.append(new_frame)\n\n    # Save into a GIF file that loops forever\n    frames[0].save(f'gif_{patient}.gif', format='GIF',\n                   append_images=frames[1:],\n                   save_all=True,\n                   duration=200, loop=0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Creating and showing the GIF**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"create_gif()\nshow_gif(filename=\"./gif_ID00199637202248141386743.gif\", format='png', width=400, height=400)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Part 3: Preprocessing","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"**Let's start by applying various types masks on the Input image**","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"**Masking type: 1**\n\nReference: https://www.kaggle.com/azaemon/starter-keras-implementation","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def split_lung_parenchyma(target,size,thr):\n    img=cv2.imdecode(np.fromfile(target,dtype=np.uint8),cv2.IMREAD_GRAYSCALE)\n    try:\n        img_thr= cv2.adaptiveThreshold(img, 255, cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY_INV,size,thr).astype(np.uint8)\n    except:\n        img_thr= cv2.adaptiveThreshold(img, 255, cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY_INV,999,thr).astype(np.uint8)\n    img_thr=255-img_thr\n    img_test=measure.label(img_thr, connectivity = 1)\n    props = measure.regionprops(img_test)\n    img_test.max()\n    areas=[prop.area for prop in props]\n    ind_max_area=np.argmax(areas)+1\n    del_array = np.zeros(img_test.max()+1)\n    del_array[ind_max_area]=1\n    del_mask=del_array[img_test]\n    img_new = img_thr*del_mask\n    mask_fill=fill_water(img_new)\n    img_new[mask_fill==1]=255\n    img_out=img*~img_new.astype(bool)\n    return img_out\n\ndef fill_water(img):\n    copyimg = img.copy()\n    copyimg.astype(np.float32)\n    height, width = img.shape\n    img_exp=np.zeros((height+20,width+20))\n    height_exp, width_exp = img_exp.shape\n    img_exp[10:-10, 10:-10]=copyimg\n    mask1 = np.zeros([height+22, width+22],np.uint8)   \n    mask2 = mask1.copy()\n    mask3 = mask1.copy()\n    mask4 = mask1.copy()\n    cv2.floodFill(np.float32(img_exp), mask1, (0, 0), 1) \n    cv2.floodFill(np.float32(img_exp), mask2, (height_exp-1, width_exp-1), 1) \n    cv2.floodFill(np.float32(img_exp), mask3, (height_exp-1, 0), 1) \n    cv2.floodFill(np.float32(img_exp), mask4, (0, width_exp-1), 1)\n    mask = mask1 | mask2 | mask3 | mask4\n    output = mask[1:-1, 1:-1][10:-10, 10:-10]\n    return output","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Let's apply this masking**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fig= plt.figure(figsize=(16,6))\nindex= './png_ID00199637202248141386743/img_30.png'\na= fig.add_subplot(1,2,1)\na.set_title('Original Image')\nplt.imshow(plt.imread(index))\nplt.grid(None)\n\na= fig.add_subplot(1,2,2)\na.set_title('Masked Image')\nimg_split=split_lung_parenchyma('./png_ID00199637202248141386743/img_30.png',15599,-66)\nplt.imshow(img_split)\nplt.grid(None)\n\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Masking type: 2**\n\nReference: https://www.raddq.com/dicom-processing-segmentation-visualization-in-python/","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    \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    \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_split=make_lungmask(cv2.imdecode(np.fromfile(index,dtype=np.uint8),cv2.IMREAD_GRAYSCALE), display=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Now let's apply Watershed Algorithm on an image**\n\n\n\nReference: https://www.kaggle.com/ankasor/improved-lung-segmentation-using-watershed","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"**First, we need to extract the internal (Lung tissue) and external (Region of interest) markers from the image**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from skimage import measure, morphology, segmentation\nimport scipy.ndimage as ndimage\n\ndef generate_markers(image):\n    #Creation of the internal Marker\n    marker_internal = image < -400\n    marker_internal = segmentation.clear_border(marker_internal)\n    marker_internal_labels = measure.label(marker_internal)\n    areas = [r.area for r in measure.regionprops(marker_internal_labels)]\n    areas.sort()\n    if len(areas) > 2:\n        for region in measure.regionprops(marker_internal_labels):\n            if region.area < areas[-2]:\n                for coordinates in region.coords:                \n                       marker_internal_labels[coordinates[0], coordinates[1]] = 0\n    marker_internal = marker_internal_labels > 0\n    #Creation of the external Marker\n    external_a = ndimage.binary_dilation(marker_internal, iterations=10)\n    external_b = ndimage.binary_dilation(marker_internal, iterations=55)\n    marker_external = external_b ^ external_a\n    #Creation of the Watershed Marker matrix\n    marker_watershed = np.zeros((512, 512), dtype=np.int)\n    marker_watershed += marker_internal * 255\n    marker_watershed += marker_external * 128\n    \n    return marker_internal, marker_external, marker_watershed\n\n#Show some example markers from the middle   \npath = '/kaggle/input/osic-pulmonary-fibrosis-progression/train/ID00010637202177584971671'\npatient_dicom = load_scan(path)\npatient_pixels = get_pixels_hu(patient_dicom)\n\ntest_patient_internal, test_patient_external, test_patient_watershed = generate_markers(patient_pixels[50])\n\nfig= plt.figure(figsize=(12,6))\na= fig.add_subplot(1,3,1)\na.set_title('Internal Marker')\nplt.imshow(test_patient_internal, cmap='gray')\nplt.grid(None)\n\na= fig.add_subplot(1,3,2)\na.set_title('External Marker')\nplt.imshow(test_patient_external, cmap='gray')\nplt.grid(None)\n\na= fig.add_subplot(1,3,3)\na.set_title('External Marker')\nplt.imshow(test_patient_watershed, cmap='gray')\nplt.grid(None)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Now we apply the marker based Watershed algorithm to find the precise border of the Lung located in the Black strip of the Watershed marker shown above.**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def seperate_lungs(image):\n    #Creation of the markers as shown above:\n    marker_internal, marker_external, marker_watershed = generate_markers(image)\n    \n    #Creation of the Sobel-Gradient\n    sobel_filtered_dx = ndimage.sobel(image, 1)\n    sobel_filtered_dy = ndimage.sobel(image, 0)\n    sobel_gradient = np.hypot(sobel_filtered_dx, sobel_filtered_dy)\n    sobel_gradient *= 255.0 / np.max(sobel_gradient)\n    \n    #Watershed algorithm\n    watershed = morphology.watershed(sobel_gradient, marker_watershed)\n    \n    #Reducing the image created by the Watershed algorithm to its outline\n    outline = ndimage.morphological_gradient(watershed, size=(3,3))\n    outline = outline.astype(bool)\n    \n    #Performing Black-Tophat Morphology for reinclusion\n    #Creation of the disk-kernel and increasing its size a bit\n    blackhat_struct = [[0, 0, 1, 1, 1, 0, 0],\n                       [0, 1, 1, 1, 1, 1, 0],\n                       [1, 1, 1, 1, 1, 1, 1],\n                       [1, 1, 1, 1, 1, 1, 1],\n                       [1, 1, 1, 1, 1, 1, 1],\n                       [0, 1, 1, 1, 1, 1, 0],\n                       [0, 0, 1, 1, 1, 0, 0]]\n    blackhat_struct = ndimage.iterate_structure(blackhat_struct, 8)\n    #Perform the Black-Hat\n    outline += ndimage.black_tophat(outline, structure=blackhat_struct)\n    \n    #Use the internal marker and the Outline that was just created to generate the lungfilter\n    lungfilter = np.bitwise_or(marker_internal, outline)\n    #Close holes in the lungfilter\n    #fill_holes is not used here, since in some slices the heart would be reincluded by accident\n    lungfilter = ndimage.morphology.binary_closing(lungfilter, structure=np.ones((5,5)), iterations=3)\n    \n    #Apply the lungfilter (note the filtered areas being assigned -2000 HU)\n    segmented = np.where(lungfilter == 1, image, -2000*np.ones((512, 512)))\n    \n    return segmented, lungfilter, outline, watershed, sobel_gradient, marker_internal, marker_external, marker_watershed\n\n#Some Testcode:\ntest_segmented, test_lungfilter, test_outline, test_watershed, test_sobel_gradient, test_marker_internal, test_marker_external, test_marker_watershed = seperate_lungs(patient_pixels[50])\n\nprint (\"Sobel Gradient\")\nplt.imshow(test_sobel_gradient, cmap='gray')\nplt.grid(None)\nplt.show()\nprint (\"Watershed Image\")\nplt.imshow(test_watershed, cmap='gray')\nplt.grid(None)\nplt.show()\n\nprint (\"Outline after reinclusion\")\nplt.imshow(test_outline, cmap='gray')\nplt.grid(None)\nplt.show()\nprint (\"Lungfilter after closing\")\nplt.imshow(test_lungfilter, cmap='gray')\nplt.grid(None)\nplt.show()\nprint (\"Segmented Lung\")\nplt.imshow(test_segmented, cmap='gray')\nplt.grid(None)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Part 4: Data Augmentation","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"**Data augmentation is a strategy that enables practitioners to significantly increase the diversity of data available for training models, without actually collecting new data. Now let's apply various types of augmentation techniques on one of the images**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from albumentations import (\n    HorizontalFlip, IAAPerspective, ShiftScaleRotate, CLAHE, RandomRotate90,\n    Transpose, ShiftScaleRotate, Blur, OpticalDistortion, GridDistortion, HueSaturationValue,\n    IAAAdditiveGaussianNoise, GaussNoise, MotionBlur, MedianBlur, IAAPiecewiseAffine,\n    IAASharpen, IAAEmboss, RandomContrast, RandomBrightness, Flip, OneOf, Compose, RandomGamma, ElasticTransform, ChannelShuffle,RGBShift, Rotate\n)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**There are various types of augmentation techniques available. In this case, we applied:**\n\n1. Horizontal Flip\n2. Random contrast\n3. Random Gamma\n4. Random Brightness\n5. Shift Scale Rotate\n6. Channel Shuffle\n7. Elastic Transform","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"albumentation_list =  [\n    HorizontalFlip(p=0.5),\n    RandomContrast(limit=0.5, p=1),\n    RandomGamma(gamma_limit=(80, 120), p=1),\n    RandomBrightness(limit=0.5, p=0.5),\n    ShiftScaleRotate(\n        shift_limit=0.0625, scale_limit=0.1, \n        rotate_limit=15, border_mode=cv2.BORDER_REFLECT_101, p=0.8), \n    ChannelShuffle(p=1),\n    ElasticTransform(p=1,border_mode=cv2.BORDER_REFLECT_101,alpha_affine=60)\n]\n\n\nchosen_image= plt.imread('/kaggle/working/png_ID00199637202248141386743/img_30.png')\nimg_matrix_list = []\nbboxes_list = []\nfor aug_type in albumentation_list:\n    img = aug_type(image = chosen_image)['image']\n    img_matrix_list.append(img)\n    \nimg_matrix_list.insert(0,chosen_image)    \n\ntitles_list = [\"Original\",\"Horizontal Flip\",\"Random Contrast\",\"Random Gamma\",\"RandomBrightness\",\n               \"Shift Scale Rotate\",\"Channel Shuffle\", \"Elastic Transform\"]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Function for plotting the augmeted images**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_multiple_img(img_matrix_list, title_list, ncols, main_title=\"Data Augmentation\"):\n    fig, myaxes = plt.subplots(figsize=(20, 15), nrows=2, ncols=ncols, squeeze=True)\n    fig.suptitle(main_title, fontsize = 30)\n    #fig.subplots_adjust(wspace=0.3)\n    #fig.subplots_adjust(hspace=0.3)\n    for i, (img, title) in enumerate(zip(img_matrix_list, title_list)):\n        myaxes[i // ncols][i % ncols].imshow(img)\n        myaxes[i // ncols][i % ncols].grid(None)\n        myaxes[i // ncols][i % ncols].set_title(title, fontsize=15)\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Now Let's see the augmented images**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_multiple_img(img_matrix_list, titles_list, ncols = 4)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Now, Let's try applying composite of augmentations on an image**\n\nReference: https://github.com/albumentations-team/albumentations/blob/a5bbfedfc500e9cac1c3689b3720769e4fa727b3/notebooks/example.ipynb","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def strong_aug(p=1):\n    return Compose([\n        RandomRotate90(),\n        Flip(),\n        Transpose(),\n        OneOf([\n            IAAAdditiveGaussianNoise(),\n            GaussNoise(),\n        ], p=0.2),\n        OneOf([\n            MotionBlur(p=.2),\n            MedianBlur(blur_limit=3, p=.1),\n            Blur(blur_limit=3, p=.1),\n        ], p=0.2),\n        ShiftScaleRotate(shift_limit=0.0625, scale_limit=0.2, rotate_limit=45, p=.2),\n        OneOf([\n            OpticalDistortion(p=0.3),\n            GridDistortion(p=.1),\n            IAAPiecewiseAffine(p=0.3),\n        ], p=0.2),\n        OneOf([\n            CLAHE(clip_limit=2),\n            IAASharpen(),\n            IAAEmboss(),\n            RandomContrast(),\n            RandomBrightness(),\n        ], p=0.3),\n        #HueSaturationValue(p=0.3),\n    ], p=p)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"aug = strong_aug(p=1)\nimg = aug(image = chosen_image)['image']\nplt.imshow(img)\nplt.grid(None)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Thank you\n\nDon't forget to upvote if you like this notebook.","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}