{"cells":[{"metadata":{},"cell_type":"markdown","source":"Most of cooler stuff in here is inspired by the following:<br>\nhttps://www.kaggle.com/allunia/pulmonary-fibrosis-dicom-preprocessing\nhttps://www.kaggle.com/ankasor/improved-lung-segmentation-using-watershed\nhttps://www.kaggle.com/arnavkj95/candidate-generation-and-luna16-preprocessing\n\nAll the ugly bits are from me :)\n\nI'm also somewhat new to this, so all commentary should be deemed speculative from this moment forward\n\nWe're given CT scan's in this competition and are asked to help predict the presence of Pulmonary Embolism (https://en.wikipedia.org/wiki/Pulmonary_embolism)\n\nI thought before the OSIC competition, I was not aware of the stark difference between CT scans and X-Rays. Basically an X-Ray is projection of the density of 3-Dimensional object into 2-D. The CT scan retains this third dimension of the information. In a slice of CT Scan, we can view each item as a pixel, but when they are stacked together to encompass some volume, we refer to it as voxel. Each voxel has a value telling the average mass density of the matter at that particular point.\n\nDid you say density?? It's technically Radio Density, which is a function of both the mass' density & the atomic number of the material in question. \n\nIn this notebook we'll do the following: <br>\n1) Briefly explore the data provided in the training dataframe <br>\n2) Explore Dicoms <br>\n3) Watershed segmentation(!) <br>\n4) Another segmentation technique (to be decided) <br>\n5) Compare the two techniques\n6) Sprinkle in (potentially) incorrect and useless musings :) <br>\n\nThere are some things to keep in mind when dealing with CT scans:<br>\n1) There does not seem to be a standard across scanners -- some manufacturers have specific settings in their scanners and the people who use them have the capacity to also fiddle with these settings as they see fit (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6115360/) -- this increases the complexity of the problem because as we will see, some scans are (way) bigger than others in terms of slices and the overall thickness of slice.<br>\n2) CT scanners are like really expensive, so the amount of CT scanning data floating around and thus the amount of documenation on tackling problems related to them is a bit spare when compared to X-Rays (well, maybe this competition and OSIC will put a dent into that)"},{"metadata":{"trusted":true,"_kg_hide-output":true,"_kg_hide-input":true,"collapsed":true},"cell_type":"code","source":"!conda install -c conda-forge gdcm -y","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os\nfrom pathlib import Path\nimport matplotlib.pyplot as plt\nimport pydicom\nimport cv2\nimport seaborn as sns\nimport gdcm\nfrom skimage import measure, segmentation, morphology\nfrom skimage.morphology import disk, opening, closing\nfrom scipy import ndimage","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"input_path = Path('../input/rsna-str-pulmonary-embolism-detection')\nos.listdir(input_path)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"train_df = pd.read_csv(input_path/'train.csv')\ntest_df = pd.read_csv(input_path/'test.csv')\nsub_df = pd.read_csv(input_path/'sample_submission.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.shape, test_df.shape, sub_df.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So we have 1,790,594 slices available for training and validation<br>\n146,853 slices we're predicting on<br>\nThere are more rows in submission file than slices in the test set -- so be mindful of how you are creating submission file!"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"list(train_df.columns)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Not entirely sure what the negative_exam_for_pe (exam level) is indicating yet, will come back around for it"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#let's just compute the % of scans with a particular attribute\npositive = train_df['pe_present_on_image'].value_counts()[1] / len(train_df['pe_present_on_image'])\nprint(\"{0:.2f}% of the training data shows Pulmonary Embolism visually\".format(positive * 100))\n\nmotion_issue = (train_df['qa_motion'].value_counts()[1] / len(train_df))\nprint(\"{0:.2f}% of the scans are noted that motion may have caused issues\".format(motion_issue*100))\n\ncontrast_issue = (train_df['qa_contrast'].value_counts()[1]) / len(train_df)\nprint(\"{0:.2f}% of the scans are noted with contrast issues\".format(contrast_issue*100))\n\n\nleft_pe = (train_df['leftsided_pe'].value_counts())\nleft_pe_pct = left_pe[1] / len(train_df['leftsided_pe'])\nprint(\"{0:.2f}% of the scans are noted with PE on left side\".format(left_pe_pct*100))\n\n\nright_pe = (train_df['rightsided_pe'].value_counts())\nright_pe_pct = right_pe[1] / len(train_df['rightsided_pe'])\nprint(\"{0:.2f}% of the scans are noted with PE on right side\".format(right_pe_pct*100))\n\ncentral_pe = (train_df['central_pe'].value_counts())\ncentral_pe_pct = central_pe[1] / len(train_df['central_pe'])\nprint(\"{0:.2f}% of the scans are noted with PE on right side\".format(central_pe_pct*100))\n\nchronic_pe = train_df['chronic_pe'].value_counts()\nchronic_pe_pct = chronic_pe[1] / len(train_df['chronic_pe'])\nprint(\"{0:.2f}% of the scans feature Chronic PE\".format(chronic_pe_pct * 100))\n\nacute_and_chronic = train_df['acute_and_chronic_pe'].value_counts()\nacute_chr_pct = acute_and_chronic[1] / len(train_df['acute_and_chronic_pe'])\nprint(\"{0:.2f}% of PE present are both acute AND Chronic :(\".format(chronic_pe_pct * 100))\n\nindeterminate = train_df['indeterminate'].value_counts()\nindeterminate_pct = indeterminate[1] / len(train_df['indeterminate'])\nprint(\"{0:.2f}% of scans had QA issues\".format(chronic_pe_pct * 100))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#no information on this in competition description\ntrain_df['flow_artifact'].value_counts() ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Should we remove the indeterminate ones for modeling?"},{"metadata":{"trusted":true},"cell_type":"code","source":"training_path = input_path/'train'\nlen(os.listdir(training_path))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ok so we have 7,729 folders inside of our training folder which represent unique 7,729 scans -- each of which has a specific number of slices. What's the distribution of these slices?\n\nInside each folder of the training folder there is another folder.... inside of this folder will you find the dicoms...."},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nscans_per_folder = []\nfor x in os.listdir(training_path):\n    path = Path(str(training_path) + '/' + str(x))\n    scans_per_folder.append(len(os.listdir(path)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"len(scans_per_folder), pd.Series(scans_per_folder).unique()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The above two lines just confirm that we only have one folder inside of each folder in the training set"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nslices_per_scan = []\nfor x in os.listdir(training_path):\n    path = Path(str(training_path) + '/' + str(x))\n    for folder in os.listdir(path):\n        scan_path = Path(str(path)+ '/' + str(folder))\n        slices_per_scan.append(len(os.listdir(scan_path)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The annoying folder structure makes it so that an easy to read list comprehension is less feasible to express cleanly"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.title('Distribution of Slices per Scan')\nplt.xlabel('Number of slices')\nplt.ylabel('Frequency')\nplt.hist(slices_per_scan, bins=50);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ok so most scans have in the range of 200-300 slices(!) -- each slices is 2-Dimensional array of numbers. So when these slices are stacked up adjacent to one another - we can volumetric information!\n\nLet's create a function that extracts the dicoms from one folder, we'll reverse the ordering because the slices towards the feet are at the start"},{"metadata":{"trusted":true},"cell_type":"code","source":"#returns a list with the dicoms in order\ndef dcm_sort(scan_path, scan_folder):\n    #a list comprehension create distill our exact file paths -- ugh\n    dcm_paths = [(str(scan_path) + '/' + file) for file in scan_folder]\n    #list comprehension that runs through each slice in the folder\n    dcm_stacked = [pydicom.dcmread(dcm) for dcm in dcm_paths]\n    dcm_stacked.sort(key=lambda x: int(x.InstanceNumber), reverse=True)\n    #returning a python list of dicoms sorted\n    return dcm_stacked","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scan_path = training_path/'858a11d72ad0/7829612362e8'\nscan_folder = os.listdir(scan_path)\nprint(\"There are {} slices in the selected scan\".format(len(scan_folder)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nsorted_scan = dcm_sort(scan_path, scan_folder)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sorted_scan[0]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A whole lot of info in one file!\n\nAs a side note, this was developed in the 80s and just like almost everything else developed by and for the medical community, it's kind of filled with a ton of garbage :)\n\nThe following resource helps interpret some of it: http://dicom.nema.org/medical/dicom/2017d/output/chtml/part03/sect_C.7.6.2.html\n\n(0020, 0013) Instance Number --> IS: \"18\" Corresponds to the fact that this is the 18th dicom for this patient (are they ordered in a spatial way?)\n\nThese two determine voxel size:\n\n(0018, 0050) Slice Thickness --> DS: \"5.0\" this is expressed in millimeters --> z-axis\n(0028, 0030) Pixel Spacing --> DS: [0.683, 0.683] physical distance between center of each pixel. The pair of values indicates adjacent row spacing and adjacent column spacing --> x/y plane\n\n(0020, 1041) Slice Location --> DS: \"82.0\" relative position of the image plane expressed in millimeters\n\n(0020, 0032) Image Position (Patient) --> DS: [-174.2187, -175.0000, 1773.500] tells the x,y,z coordinates of the upper left hand corner (center of the firxt voxel transmitted) of image in millimeters\n\nWe will also take into account the Rescale slope & intercept below in order to convert the pixel values into the Hounsfield scale shortly: https://en.wikipedia.org/wiki/Hounsfield_scale"},{"metadata":{"trusted":true},"cell_type":"code","source":"sorted_scan[0].PixelData[0:100]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Yikes, that's not something we can work with... luckily we can access the pixel information with the following:"},{"metadata":{"trusted":true},"cell_type":"code","source":"sorted_scan[0].pixel_array[0:5]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Take note of that -2000 at all four corners of this slice -- we'll come back to that shortly\n\nWe know there are 261 slices in this scan and now that we have them in order, let's pick something near the middle so that we can see more of the lungs(!)"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(sorted_scan[100].pixel_array);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"That's not very easy to see, change the lens we're interpreting the pixel information with by choosing a better color map"},{"metadata":{"trusted":true},"cell_type":"code","source":"#let's concentrate on the a section of 60 slices in this scan\nmiddle_scan = sorted_scan[80:140]\n\nfig,ax = plt.subplots(4,5, figsize=(12,8))\nfor n in range(4):\n    for m in range(5):\n        ax[n,m].imshow(middle_scan[n*5+m].pixel_array, cmap='Blues_r')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#let's take a look at a single slice's pixel distribution\nplt.hist(middle_scan[20].pixel_array);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we mentioned -- there's a ton of pixels at -2000.... why?<br>\nThe scan represents areas outside of the body as these extremely low values -- let's set them to 0, which represents water"},{"metadata":{"trusted":true},"cell_type":"code","source":"one_slice = middle_scan[20].pixel_array\none_slice[one_slice <= -1000] = 0\nplt.imshow(one_slice, cmap='Blues_r');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ok, that looks better in terms of the border but what's up with that weird thing at the bottom of the image? Is it the surface the patient is laying on? We'll have to ensure the segmentation technique we use can filter this out\n\nLet's take into account the rescale intercept and slope(!)"},{"metadata":{"trusted":true},"cell_type":"code","source":"one_slice = middle_scan[20]\none_slice.RescaleIntercept, one_slice.RescaleSlope","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def scan_transformed_hu(dcm_sorted, threshold=-1000, replace=-1000):\n    intercept = dcm_sorted[0].RescaleIntercept\n    slices_stacked = np.stack([dcm.pixel_array for dcm in dcm_sorted])\n    slices_stacked = slices_stacked.astype(float)\n    \n    #converts the unknown values to desired replacement\n    slices_stacked[slices_stacked <= threshold] = replace\n    \n    #turn into hounsfield scale\n    slices_stacked += np.int16(intercept)\n    \n    return np.array(slices_stacked, dtype=np.int16)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"middle_slices_hu = scan_transformed_hu(middle_scan, replace=0)\n\nfig,ax = plt.subplots(12,5, figsize=(20,20))\nfor n in range(12):\n    for m in range(5):\n        ax[n,m].imshow(middle_slices_hu[n*5+m], cmap='Blues_r')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ok these look better but you can still see the difference in contrast between slices\n\nLet's start some segmentation -- the notebooks above list the concise functions of segmentations we'll explore. But we'll try to step through the segmentations one line at a time"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"test_slice = middle_slices_hu[10]\nfig, ax = plt.subplots(1,2, figsize=(12,3))\nax[0].imshow(test_slice)\nax[0].set_title('Slice #90') #middle scan was 80-120 and this is 10th one....\nax[1].set_title('Pixel Distribution of Slice')\nax[1].hist(test_slice);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"First we'll want to threshold the image -- this means we want everything below or above a certain condition to be set as true and everything else is set as false -- we know that lung tissue should be around -400 and air is -1000 so let's pick something a bit higher to be safe -- this will represent our internal marker"},{"metadata":{"trusted":true},"cell_type":"code","source":"internal_marker = test_slice < -300\ninternal_marker[203:207, 203:220] #just to show the discrepeny in middle somewhere","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#this represents the region we know definitely features lung tissue\nplt.title('prelimary internal marker')\nplt.imshow(segmentation.clear_border(internal_marker), cmap='gray');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"internal_marker_labels = measure.label(segmentation.clear_border(internal_marker))\nplt.imshow(internal_marker_labels, cmap='gray');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#explicating the next list comprehension\nmeasure.regionprops(internal_marker_labels)[0:3]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"areas = [x.area for x in measure.regionprops(internal_marker_labels)]\nareas.sort()\nareas","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for region in measure.regionprops(internal_marker_labels):\n    if region.area < areas[-2]:\n        for coordinates in region.coords:\n            internal_marker_labels[coordinates[0], coordinates[1]] = 0","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"marker_internal = internal_marker_labels > 0\nplt.title('Internal marker')\nplt.imshow(marker_internal, cmap='gray');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"If you look closely this is NOT the same as the preliminary marker from above -- specifically those tiny white spots below the lungs are no longer visible\n\nNow we want to generate the external marker - which is the area we know is outside our region of interest. <br>\n\nHow is this done??<br> A morphological dilation of the internal marker, with two iterations done and then we'll find the difference<br>\n\nWhat does that mean?<br>\nLet's take a look first"},{"metadata":{"trusted":true},"cell_type":"code","source":"external_a = ndimage.binary_dilation(marker_internal, iterations=10)\nexternal_b = ndimage.binary_dilation(marker_internal, iterations=50)\nmarker_external = external_b ^ external_a\n#since they're set to binary values - finding sum will tell you how much white\nexternal_a.sum(), external_b.sum() ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1, 3, figsize=(20,8))\nax[0].imshow(external_a, cmap='gray')\nax[0].set_title('dilation of internal marker - 10 iterations')\nax[1].imshow(external_b, cmap='gray')\nax[1].set_title('dilation of internal marker - 50 iterations')\nax[2].set_title('external marker')\nax[2].imshow(marker_external, cmap='gray');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ok, conceptually -- what just happened?\n\nFrom the internal marker that we set -- we used the dilation operation to gradually enlarge the boundaries of regions of forground pixels (white pixels).<br>\n\n1) Areas of foreground pixels grows in size -- so the boundary of the lungs are expanding<br>\n\n2) Holes within those regions become smaller --> notice how in the first iteration you already no longer see all those pockets of black inside of the lungs<br>\n\n3) by taking the difference, you find an outline that is beyond the lungs themselves to ensure that you are only pulling that regions\n\nThis was originally implemented on a dataset/competition where you are looking for nodules in the lungs\n\nWhich would end up on the lining of the lungs themselves, so they wanted to ensure they captured that boundary point fully -- unsure if that's needed here\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"watershed_marker = np.zeros((512, 512), dtype=np.int)\nwatershed_marker += marker_internal * 255 #high intensity\nwatershed_marker += marker_external * 128 #medium intensity","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The above items do the following:<br>\n1) initialze an empty array<br>\n2) set values known to be lungs as the highest intensity<br>\n3) set values of the external marker to be of medium intensity<br>"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.title('Watershed marker!')\nplt.imshow(watershed_marker, cmap='gray');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we'll be using Sobel Kernels (https://en.wikipedia.org/wiki/Sobel_operator) which is just two convolution kernels with set weights that we use to compute over the image to determine edges in both the x and y plane"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1, 2, figsize=(20,8))\nax[0].imshow(ndimage.sobel(test_slice, 0), cmap='gray')\nax[0].set_title('vertical edges')\nax[1].imshow(ndimage.sobel(test_slice, 1), cmap='gray')\nax[1].set_title('horizontal edges');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"x_edges = ndimage.sobel(test_slice, 1)\ny_edges = ndimage.sobel(test_slice, 0)\nsobel_grad = np.hypot(x_edges, y_edges)\nsobel_grad *= 255.0 / np.max(sobel_grad)\nplt.title('sobel gradient')\nplt.imshow(sobel_grad, cmap='gray');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So we used sobel algorithm to find edges in our image and now we'll take the output of this and apply the watershed segmentation algorithm:\nhttps://scikit-image.org/docs/dev/auto_examples/segmentation/plot_watershed.html\n\nThe algorithm takes the watershed markers we found earlier -- it then treats the image topographically: pixel values are viewed as having a certain elevation based on their intensity. It floors basins from these markers until a basin from another make contact"},{"metadata":{"trusted":true},"cell_type":"code","source":"img_watershed = segmentation.watershed(test_slice, watershed_marker)\nwatershed = segmentation.watershed(sobel_grad, watershed_marker)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2, figsize=(20,8))\nax[0].imshow(img_watershed, cmap='gray')\nax[0].set_title('watershed seg w/ original img')\nax[1].set_title('watershed seg w/ sobel gradient')\nax[1].imshow(watershed, cmap='gray');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see the output via the sobel gradient image gives us smoother edges of the lungs vs the more jagged texture as seen with the original image\n\nNext we reduce what we have to an outline!\nThis is done with the morphological_gradient fxn in scipy, which seems to take a morphological dilation and a morphological erosion of the input and then finds the difference of the two: https://scipy.github.io/devdocs/generated/scipy.ndimage.morphological_gradient.html#scipy.ndimage.morphological_gradient\n\nBased on that link from earlier -- the erosion takes the binary image and erodes away at the pixels in the boundary regions -- shrinking the foreground in size. Since the dilation is doing the opposite, by finding the difference between the two we can outline the lungs :)"},{"metadata":{"trusted":true},"cell_type":"code","source":"#let's try out different kernel sizes :)\nfig, ax = plt.subplots(1, 3, figsize=(20,8))\nax[0].imshow(ndimage.morphological_gradient(watershed, size=(2,2)))\nax[0].set_title('outline derived from 2x2 kernel')\nax[1].imshow(ndimage.morphological_gradient(watershed, size=(3,3)))\nax[1].set_title('outline derived from 3x3 kernel')\nax[2].set_title('outline derived from 7x7 kernel')\nax[2].imshow(ndimage.morphological_gradient(watershed, size=(7, 7)));","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We definitely don't want to select the 2x2 kernel -- it can't even trace the bottom right section of this slice -- the other implementations of this all have a 3x3 kernel, but it might be interesting to see what happens if we were to use a larger kernel -- that outline is bold but perhaps unneccessary"},{"metadata":{"trusted":true},"cell_type":"code","source":"outline = ndimage.morphological_gradient(watershed, size=(3,3))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now onto black (top) hat morphology..... the documentation for scipy provides no explanation. openCV documentation to the rescue :) <br>\nhttps://docs.opencv.org/trunk/d9/d61/tutorial_py_morphological_ops.html <br>\nBasically it's the difference between closing of the input image and the input image. What does closing mean?<br>\n\nClosing is a dilation followed by an erosion -- it allows you to close small holes inside the foreground objects.<br>\n\nDilation: you use a kernel still and a pixel element is 1 if at least ONE pixel under the kernel is 1. Meaning, it increases the white region in the image.<br>\n\nErosion: a pixel is consdiered a 1 only if all the pixels under that kernel are 1 -- otherwise it's eroded."},{"metadata":{"trusted":true},"cell_type":"code","source":"#openCV has kernel fxns - not so with scipy, hmm\nblackhat_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\nblackhat_kernel = ndimage.iterate_structure(blackhat_struct, 8)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"blackhat_outline = outline + ndimage.black_tophat(outline,\n                                structure=blackhat_kernel)\n\nplt.title('Blackhat outline')\nplt.imshow(blackhat_outline, cmap='gray');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"lung_filter = np.bitwise_or(marker_internal, blackhat_outline)\nplt.title('lung filter')\nplt.imshow(lung_filter, cmap='gray');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"lung_filter = ndimage.morphology.binary_closing(lung_filter,\n                structure=np.ones((5,5)), iterations=3)\nplt.title('Lung Filter via internal marker and blackhat outline')\nplt.imshow(lung_filter, cmap='gray');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#where you see a 1 in lung filter - put the actual pixel value\n#everywhere else include -2000\nplt.title('Our segmented slice!!')\nplt.imshow(np.where(lung_filter == 1, test_slice, -2000), cmap='gray');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"That looks pretty good! Of course it isn't exactly lining up along the border of lungs, but that's totally fine for the high level purposes we're trying for here. \n\nOk let's wrap all of that up into a function that operates over each slice in a scan one by one -- of just applies it to a single slice if that's all that's fed to it\n\nNo commentary below - if confused, check the cells ran above :)"},{"metadata":{"trusted":true},"cell_type":"code","source":"def gen_internal_marker(slices_s, threshold= -300):\n    internal_marker = slices_s < threshold\n    internal_marker_labels = measure.label(segmentation.clear_border(internal_marker))\n    areas = [x.area for x in measure.regionprops(internal_marker_labels)]\n    areas.sort()\n    for region in measure.regionprops(internal_marker_labels):\n        if region.area < areas[-2]:\n            for coordinates in region.coords:\n                internal_marker_labels[coordinates[0], coordinates[1]] = 0\n    marker_internal = internal_marker_labels > 0\n                \n    return marker_internal\n\ndef gen_external_marker(internal_marker, iter_1 = 10, iter_2 = 50):\n    external_a = ndimage.binary_dilation(internal_marker, \n                                         iterations=iter_1)\n    external_b = ndimage.binary_dilation(internal_marker, \n                                         iterations=iter_2)\n    external_marker = external_b ^ external_a\n    return external_marker\n\ndef gen_watershed_marker(internal_marker, external_marker):\n    watershed_marker = np.zeros((512, 512), dtype=np.int)\n    watershed_marker += internal_marker * 255\n    watershed_marker += external_marker * 128\n    return watershed_marker\n\ndef gen_sobel_grad(one_slice):\n    x_edges = ndimage.sobel(one_slice, 1)\n    y_edges = ndimage.sobel(one_slice, 0)\n    sobel_grad = np.hypot(x_edges, y_edges)\n    sobel_grad *= 255.0 / np.max(sobel_grad)\n    return sobel_grad\n\ndef gen_blackhat_outline(watershed, blackhat_struct, b_hat_iters=1):\n    outline = ndimage.morphological_gradient(watershed, size=(3,3))\n    blackhat_kernel = ndimage.iterate_structure(blackhat_struct, \n                                               b_hat_iters)\n    blackhat_outline = outline + ndimage.black_tophat(outline, \n                            structure=blackhat_kernel)\n    return blackhat_outline\n\n\ndef gen_lung_filter(internal_marker, blackhat_outline,\n                   kernel_size=(5,5), iterations=3):\n    pre_filter = np.bitwise_or(internal_marker, blackhat_outline)\n    lung_filter = ndimage.morphology.binary_closing(pre_filter,\n                            structure=np.ones(kernel_size),\n                            iterations=iterations)\n    return lung_filter\n\n\ndef watershed_seg(slice_s, blackhat_struct, threshold=-350,\n                  b_hat_iters=1, iter_1=10, iter_2=50):\n    \n    scan = [] #initialize an empty list\n    for one_slice in slice_s:\n        internal_marker = gen_internal_marker(one_slice)\n\n        external_marker = gen_external_marker(internal_marker)\n        \n        watershed_marker = gen_watershed_marker(internal_marker,\n                                               external_marker)\n        \n        sobel_grad = gen_sobel_grad(one_slice)\n       \n        watershed = segmentation.watershed(sobel_grad, \n                                           watershed_marker)\n        \n        blackhat_outline = gen_blackhat_outline(watershed,\n                                blackhat_struct, b_hat_iters)\n        \n        lung_filter = gen_lung_filter(internal_marker,\n                                      blackhat_outline)\n        \n        segmented_slice = np.where(lung_filter == 1, one_slice, -2000)\n        scan.append(segmented_slice)\n        \n    return np.array(scan)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#just a reminder of the shape of the middle section of the scan we pulled earlier\nmiddle_slices_hu.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's look at the differences between 1 and 6 iterations of the blackhat fxn -- we'll consider both the time and segmentation results"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nsegmented_scan_1 = watershed_seg(middle_slices_hu, blackhat_struct,\n                                b_hat_iters=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nsegmented_scan_6 = watershed_seg(middle_slices_hu, blackhat_struct,\n                              b_hat_iters=6)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1, 2, figsize=(20,8))\nax[0].imshow(segmented_scan_1[32], cmap='Blues_r')\nax[1].imshow(segmented_scan_6[32], cmap='Blues_r');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"7,279 scans in train folder * 112 seconds means this would take 226 hours to segment all that with 6 iterations, that's not feasible. At one iteration it'll still take just over a whole day on a kaggle kernel to complete(!) \n\nThe number of iterations done in the blackhat function expands the segmentation mask a bit but at a cost of over 9 times the time it takes for a single iteration."},{"metadata":{"trusted":true},"cell_type":"code","source":"fig,ax = plt.subplots(12,5, figsize=(20,20))\nfor n in range(12):\n    for m in range(5):\n        ax[n,m].imshow(segmented_scan_1[n*5+m], cmap='Blues_r')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A note to keep in mind. We did this on the middle 60 slices of the scan above, had we chosed one of the extremes, towards either the feet or head -- we would end up not seeing much or any lung at all and the segmentation technique would might end up picking up things we don't want, like the table/surface that the patient is laying on. This is NOT ideal and probably something we would want to select against. However in this competition, the segmentation seems to work better than in the OSIC one, not sure why. \n\n**To-Do List:**<br>\n1) implement another segmentation method - most likely the one suggested by @[allunia](https://www.kaggle.com/allunia) in her awesome OSIC notebook!\n\n2) compare the differences in the segmentations"},{"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}