{"cells":[{"metadata":{},"cell_type":"markdown","source":"
\n

Pulmonary Embolism

\n \n \n
","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## About Pulmonary Embolism:\n\nPulmonary embolism (PE) is a sudden blockage in a lung artery. It usually happens when a blood clot breaks loose and travels through the bloodstream to the lungs. \n\nPE is a serious condition that can cause:\n\n- Permanent damage to the lungs\n- Low oxygen levels in your blood\n- Damage to other organs in your body from not getting enough oxygen\n\nPE can be life-threatening, especially if a clot is large, or if there are many clots.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## How can we help?\n\nWe can use ML/DL techniques which could help to more accurately identify PE cases, making management and treatment more effective for patients.\n\nThe competition data consists of CT scans, which consist of hundreds of images that require detailed review to identify clots within the pulmonary arteries.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Table of Contents:\n\n- ## 1.1 [Lung Segmentation](https://www.kaggle.com/aadhavvignesh/making-use-of-ct-scans-pulmonary-embolism#1.1-Lung-Segmentation)\n- ## 1.2 [Calculating Lung Volume](https://www.kaggle.com/aadhavvignesh/making-use-of-ct-scans-pulmonary-embolism#1.2-Calculating-Lung-Volumes-from-CT-Scans)\n- ## 1.3 [Inspecting Tabular Data](https://www.kaggle.com/aadhavvignesh/making-use-of-ct-scans-pulmonary-embolism#1.3-Checking-the-tabular-data)","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### This work is based on my [previous work](https://www.kaggle.com/aadhavvignesh/lung-segmentation-by-marker-controlled-watershed).","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# 1.1 Lung Segmentation\nThe basic idea behind watershed segmentation is that any grayscale can be considered as a topographic surface.\nIf we flood the surface from its minima, and successfully prevent merging of waters, we partition the image into two different sets: the catchment basins and the watershed lines.\n\n![Watershed](http://www.cmm.mines-paristech.fr/~beucher/lpe1.gif)![Final Watersheds](http://www.cmm.mines-paristech.fr/~beucher/ima3.gif)\n\n> Image Source: [CMM](http://www.cmm.mines-paristech.fr/~beucher/wtshed.html)","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"We'll be using `pydicom` for dealing with the scans, feel free to use any available library:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install pydicom","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%matplotlib inline\n\nimport numpy as np\nimport pandas as pd\nimport pydicom\nimport os\nimport scipy.ndimage as ndimage\nfrom skimage import measure, morphology, segmentation\nimport matplotlib.pyplot as plt\n\nimport time","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's load the patients' data:","execution_count":null},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"INPUT_FOLDER = '/kaggle/input/rsna-str-pulmonary-embolism-detection/train/'\n\npatients_studyuid = os.listdir(INPUT_FOLDER)\npatients_studyuid.sort()\n\nprint(\"Some examples of patient Study UIDs:\")\nprint(\",\\n\".join(patients_studyuid[:5]))\n\n\nseries_uids = []\nfor study_uid in patients_studyuid[:5]:\n series_uid = os.listdir(INPUT_FOLDER + study_uid)[0]\n# print(\"Series UID:\", series_uid)\n series_uids.append(INPUT_FOLDER + study_uid + '/' + series_uid)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now let's load the scans:\n\nThis code for loading the scans is based on Franklin Heng's [Medium article.](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 \"\"\"\n Loads scans from a folder and into a list.\n \n Parameters: path (Folder path)\n \n Returns: slices (List of slices)\n \"\"\"\n \n slices = [pydicom.read_file(path + '/' + s) for s in os.listdir(path)]\n slices.sort(key = lambda x: int(x.InstanceNumber))\n \n try:\n slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\n except:\n slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n \n for s in slices:\n s.SliceThickness = slice_thickness\n \n return slices","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Hounsfield Units\n\nThe unit of measurement in CT scans is the Hounsfield Unit (HU), which is a measure of radiodensity. \n\n**Hounsfield units (HU)** are a dimensionless unit universally used in computed tomography (CT) scanning to express CT numbers in a standardized and convenient form. Hounsfield units are obtained from a linear transformation of the measured attenuation coefficients.\n\nHUs can be calculated from the pixel data with a DICOM Image using the following formula:\n\n$\\large HU = m*P + b$\n\nwhere,\n\n$m$ = `RescaleSlope` attribute of the DICOM image,\n\n$b$ = `RescaleIntercept` attribute of the DICOM image,\n\n$P$ = Pixel Array","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_pixels_hu(scans):\n \"\"\"\n Converts raw images to Hounsfield Units (HU).\n \n Parameters: scans (Raw images)\n \n Returns: image (NumPy array)\n \"\"\"\n \n image = np.stack([s.pixel_array for s in scans])\n image = image.astype(np.int16)\n\n # Since the scanning equipment is cylindrical in nature and image output is square,\n # we set the out-of-scan pixels to 0\n image[image == -2000] = 0\n \n \n # HU = m*P + b\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":{},"cell_type":"markdown","source":"Let's store store the slices and the images:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"test_patient_scans = load_scan(series_uids[0])\ntest_patient_images = get_pixels_hu(test_patient_scans)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We'll be taking a random slice to perform segmentation:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(test_patient_images[100], cmap='gray')\nplt.title(\"Original Slice\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Marker-Controlled Watershed Transformation\n\n**Watershed Transform** is a really powerful segmentation algorithm, but has a drawback where every regional minimum forms its own catchment basin. To overcome this drawback, we use a marker-controlled watershed transformation, where we manually create markers where we start the flooding process.\n\n## About the Algorithm:\n\nThe image is seen as a topographical surface where grey values are deemed as elevation of the surface at that location. Then, flooding process starts in which water effuses out of the minimum grey value or the marker. When flooding across two minimum converges then a dam is built to identify the boundary across them.\n\n\n![Markers](http://www.cmm.mines-paristech.fr/~beucher/ima4.gif)\n![Flood](http://www.cmm.mines-paristech.fr/~beucher/lpe2.gif)","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Marker Generation:\n\nFor using marker-controlled watershed segmentation, we'll need to identify markers. Internal marker, which is our region of interest, i.e lung tissue and an external marker, which is the region outside of our region of interest.\n\nWe create the external marker is created by morphological dilation of the internal marker, by iterating twice and subtracting the results. The watershed marker is created by superimposing both the markers.\n\nSome of the code is based from @arnavkj95's kernel: https://www.kaggle.com/arnavkj95/candidate-generation-and-luna16-preprocessing","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def generate_markers(image):\n \"\"\"\n Generates markers for a given image.\n \n Parameters: image\n \n Returns: Internal Marker, External Marker, Watershed Marker\n \"\"\"\n \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 \n areas = [r.area for r in measure.regionprops(marker_internal_labels)]\n areas.sort()\n \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 \n marker_internal = marker_internal_labels > 0\n \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 \n # Creation of the Watershed Marker\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","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's get our markers for the sample slice:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"test_patient_internal, test_patient_external, test_patient_watershed = generate_markers(test_patient_images[100])\n\nf, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True, figsize=(15,15))\n\nax1.imshow(test_patient_internal, cmap='gray')\nax1.set_title(\"Internal Marker\")\nax1.axis('off')\n\nax2.imshow(test_patient_external, cmap='gray')\nax2.set_title(\"External Marker\")\nax2.axis('off')\n\nax3.imshow(test_patient_watershed, cmap='gray')\nax3.set_title(\"Watershed Marker\")\nax3.axis('off')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Sobel Gradient and Edge Outlining\n\nThe Sobel operator performs a 2D spatial gradient measurement on an image and so emphasizes regions of high spatial frequency that correspond to edges.\n\nIt consists of a pair of 3×3 convolution kernels.\n\n![Conv Filters](http://homepages.inf.ed.ac.uk/rbf/HIPR2/figs/sobmasks.gif)\n\nThese kernels can then be combined together to find the absolute magnitude of the gradient at each point and the orientation of that gradient.\n\nThe gradient magnitude is given by:\n$G = sqrt(Gx^2 + Gy^2)$\n\nSobel gradient can be calculated by `scipy.ndimage`.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Lists to store computation times and iterations\ncomputation_times = []\niteration_titles = []","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def seperate_lungs(image, iterations = 1):\n \"\"\"\n Segments lungs using various techniques.\n \n Parameters: image (Scan image), iterations (more iterations, more accurate mask)\n \n Returns: \n - Segmented Lung\n - Lung Filter\n - Outline Lung\n - Watershed Lung\n - Sobel Gradient\n \"\"\"\n \n # Store the start time\n start = time.time()\n \n marker_internal, marker_external, marker_watershed = generate_markers(image)\n \n \n '''\n Creation of Sobel Gradient\n '''\n \n # 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 \n '''\n Using the watershed algorithm\n \n \n We pass the image convoluted by sobel operator and the watershed marker\n to morphology.watershed and get a matrix matrix labeled using the \n watershed segmentation algorithm.\n '''\n watershed = morphology.watershed(sobel_gradient, marker_watershed)\n \n '''\n Reducing the image to outlines after Watershed algorithm\n '''\n outline = ndimage.morphological_gradient(watershed, size=(3,3))\n outline = outline.astype(bool)\n \n \n '''\n Black Top-hat Morphology:\n \n The black top hat of an image is defined as its morphological closing\n minus the original image. This operation returns the dark spots of the\n image that are smaller than the structuring element. Note that dark \n spots in the original image are bright spots after the black top hat.\n '''\n \n # Structuring element used for the filter\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 \n blackhat_struct = ndimage.iterate_structure(blackhat_struct, iterations)\n \n # Perform Black Top-hat filter\n outline += ndimage.black_tophat(outline, structure=blackhat_struct)\n \n '''\n Generate lung filter using internal marker and outline.\n '''\n lungfilter = np.bitwise_or(marker_internal, outline)\n lungfilter = ndimage.morphology.binary_closing(lungfilter, structure=np.ones((5,5)), iterations=3)\n \n '''\n Segment lung using lungfilter and the image.\n '''\n segmented = np.where(lungfilter == 1, image, -2000*np.ones((512, 512)))\n \n # Append computation time\n end = time.time()\n computation_times.append(end - start)\n iteration_titles.append(\"{num} iterations\".format(num = iterations))\n \n \n return segmented, lungfilter, outline, watershed, sobel_gradient","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Comparison of iterations with time\n\nWe'll be checking for iterations in the range of 1-8. `iterations = 1` is the default for the `seperate_lungs` function.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"for itrs in range(1, 9):\n test_segmented, test_lungfilter, test_outline, test_watershed, test_sobel_gradient = seperate_lungs(test_patient_images[100], itrs)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"itr_dict = {'Iterations' : iteration_titles, 'Computation Times (in seconds)': computation_times}\n\ncolors = ['#30336b',] * 8\ncolors[0] = '#eb4d4b'\n\nimport plotly.express as px\nimport plotly.graph_objects as go\n\n\nfig = go.Figure(data=[go.Bar(\n x=itr_dict['Iterations'],\n y=itr_dict['Computation Times (in seconds)'],\n marker_color = colors\n)])\nfig.update_traces(texttemplate='%{y:.3s}', textposition='outside')\n\n\nfig.update_layout(\n title = 'Iterations vs Computation Times',\n yaxis=dict(\n title='Computation Times (in seconds)',\n titlefont_size=16,\n tickfont_size=14,\n ),\n autosize=False,\n width=800,\n height=800)\n\nfig.show()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"f, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize = (12, 12))\n\nax1.imshow(test_sobel_gradient, cmap='gray')\nax1.set_title(\"Sobel Gradient\")\nax1.axis('off')\n\nax2.imshow(test_watershed, cmap='gray')\nax2.set_title(\"Watershed\")\nax2.axis('off')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"f, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True, figsize = (12, 12))\n\nax1.imshow(test_outline, cmap='gray')\nax1.set_title(\"Lung Outline\")\nax1.axis('off')\n\nax2.imshow(test_lungfilter, cmap='gray')\nax2.set_title(\"Lung filter\")\nax2.axis('off')\n\nax3.imshow(test_segmented, cmap='gray')\nax3.set_title(\"Segmented Lung\")\nax3.axis('off')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"f, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize = (12, 12))\n\nax1.imshow(test_patient_images[100], cmap='gray')\nax1.set_title(\"Original Lung\")\nax1.axis('off')\n\nax2.imshow(test_segmented, cmap='gray')\nax2.set_title(\"Segmented Lung\")\nax2.axis('off')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 1.2 Calculating Lung Volumes from CT Scans","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def generate_markers_set(images):\n return np.array([generate_markers(img)[0] for img in images])\n\ndef calculate_lungvolume(patient_scans, patient_images):\n lung_volume = 0\n patient_masks = generate_markers_set(patient_images)\n \n for idx,patient_scan in enumerate(patient_scans):\n pixel_spacing = patient_scan.PixelSpacing\n slice_thickness = patient_scan.SliceThickness\n lung_volume += np.count_nonzero(patient_masks[idx]) * pixel_spacing[0] * pixel_spacing[0] * slice_thickness\n \n # Returns volume in cm3\n return lung_volume*0.001","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"lung_vol = calculate_lungvolume(test_patient_scans, test_patient_images)\nprint(\"The lung volume of the first patient is:\", lung_vol, \"cm^3\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 1.3 Checking the tabular data","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = pd.read_csv(\"../input/rsna-str-pulmonary-embolism-detection/train.csv\")\ntest_df = pd.read_csv(\"../input/rsna-str-pulmonary-embolism-detection/test.csv\")\nsubm_df = pd.read_csv(\"../input/rsna-str-pulmonary-embolism-detection/sample_submission.csv\")\n\ntrain_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":"subm_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# More to come!","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Conclusion:\n\nIn this kernel, we discussed about:\n\n- What is Watershed Transformation?\n- Loading Scans\n- Creating Internal, External and Watershed Markers\n- Getting the segmented lung images\n- Using the segmented lung images for calculating volume.\n- Inspecting tabular data.","execution_count":null},{"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}