{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.7.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceType":"competition","sourceId":20604,"databundleVersionId":1357052}],"dockerImageVersionId":29980,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# UPDATE:\n\nThe notebook has been updated and **works much faster (8x less processing time).**\n\n### Each slice now takes approx. 120-220 ms, compared to the 2.5 - 5 seconds processing time.","metadata":{}},{"cell_type":"markdown","source":"## Update 2: I've been working on finding lung volumes and chest circumferences, I'll keep you updated once it's done!","metadata":{}},{"cell_type":"markdown","source":"# Lung Segmentation using Marker-Controlled Watershed Transformation\n\nPreviously-suggested segmentation methods use image thresholding based on HU (Hounsfield Value) or other methods which are susceptible for picking up regions which aren't of our interest. **Watershed Transform** is a really powerful segmentation algorithm which is based on [watersheds](https://science.howstuffworks.com/environmental/conservation/issues/watershed1.htm) where we think the image as a surface.\n\n\n![Blobs](https://www.mathworks.com/company/newsletters/articles/the-watershed-transform-strategies-for-image-segmentation/_jcr_content/mainParsys/image_1.adapt.1200.high.jpg/1542750811892.jpg)![Catchment Basins](https://www.mathworks.com/company/newsletters/articles/the-watershed-transform-strategies-for-image-segmentation/_jcr_content/mainParsys/image_2.adapt.1200.high.gif/1542750811908.gif)\n\n## Watershed Transformation\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)","metadata":{}},{"cell_type":"markdown","source":"We'll be using `pydicom` for dealing with the scans, feel free to use any available library:","metadata":{}},{"cell_type":"code","source":"!pip install pydicom","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T14:07:29.739702Z","iopub.execute_input":"2026-04-27T14:07:29.740128Z","iopub.status.idle":"2026-04-27T14:07:36.886858Z","shell.execute_reply.started":"2026-04-27T14:07:29.740091Z","shell.execute_reply":"2026-04-27T14:07:36.885586Z"}},"outputs":[],"execution_count":null},{"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","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T12:52:52.892021Z","iopub.execute_input":"2026-04-27T12:52:52.892460Z","iopub.status.idle":"2026-04-27T12:52:54.606214Z","shell.execute_reply.started":"2026-04-27T12:52:52.892419Z","shell.execute_reply":"2026-04-27T12:52:54.605061Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Let's load the patients' data:","metadata":{}},{"cell_type":"code","source":"INPUT_FOLDER = '/kaggle/input/osic-pulmonary-fibrosis-progression/train/'\n\npatients = os.listdir(INPUT_FOLDER)\npatients.sort()\n\nprint(\"Some examples of patient IDs:\")\nprint(\",\\n\".join(patients[:5]))","metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"execution":{"iopub.status.busy":"2026-04-27T12:52:58.610649Z","iopub.execute_input":"2026-04-27T12:52:58.611010Z","iopub.status.idle":"2026-04-27T12:52:58.623100Z","shell.execute_reply.started":"2026-04-27T12:52:58.610980Z","shell.execute_reply":"2026-04-27T12:52:58.622277Z"}},"outputs":[],"execution_count":null},{"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)","metadata":{}},{"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","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T12:53:33.073525Z","iopub.execute_input":"2026-04-27T12:53:33.073878Z","iopub.status.idle":"2026-04-27T12:53:33.082576Z","shell.execute_reply.started":"2026-04-27T12:53:33.073848Z","shell.execute_reply":"2026-04-27T12:53:33.081405Z"}},"outputs":[],"execution_count":null},{"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\n![HU Table](http://patentimages.storage.googleapis.com/WO2005055806A2/imgf000011_0001.png)\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","metadata":{}},{"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)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T12:53:38.649614Z","iopub.execute_input":"2026-04-27T12:53:38.650023Z","iopub.status.idle":"2026-04-27T12:53:38.658085Z","shell.execute_reply.started":"2026-04-27T12:53:38.649987Z","shell.execute_reply":"2026-04-27T12:53:38.656892Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Let's store store the slices and the images:","metadata":{}},{"cell_type":"code","source":"test_patient_scans = load_scan(INPUT_FOLDER + patients[24])\ntest_patient_images = get_pixels_hu(test_patient_scans)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"We'll be taking a random slice to perform segmentation:","metadata":{}},{"cell_type":"code","source":"plt.imshow(test_patient_images[12], cmap='gray')\nplt.title(\"Original Slice\")\nplt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Marker-Controlled Watershed Transformation\n\n**Watershed Transform** is a really powerful segmentation algorithm, but has a drawback:\n\n- **Over Segmentation:** Oversegmentation occurs because every regional minimum forms its own catchment basin. Here is an example where steel grains are over-segmented by watershed transformation:\n\n![Steel Grains](https://www.mathworks.com/company/newsletters/articles/the-watershed-transform-strategies-for-image-segmentation/_jcr_content/mainParsys/image_9.adapt.1200.high.gif/1542750812181.gif)![Oversegmented](https://www.mathworks.com/company/newsletters/articles/the-watershed-transform-strategies-for-image-segmentation/_jcr_content/mainParsys/image_10.adapt.1200.high.gif/1542750812206.gif)\n\n> **Left:** Steel Grains, **Right:** Oversegmented image as a result of using normal watershed transformation.\n\nTo 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)","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","metadata":{}},{"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","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-27T13:03:40.960600Z","iopub.execute_input":"2026-04-27T13:03:40.960962Z","iopub.status.idle":"2026-04-27T13:03:40.971409Z","shell.execute_reply.started":"2026-04-27T13:03:40.960930Z","shell.execute_reply":"2026-04-27T13:03:40.970240Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Let's get our markers for the sample slice:","metadata":{}},{"cell_type":"code","source":"test_patient_internal, test_patient_external, test_patient_watershed = generate_markers(test_patient_images[12])\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()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"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`.","metadata":{}},{"cell_type":"code","source":"# Lists to store computation times and iterations\ncomputation_times = []\niteration_titles = []","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"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","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"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.","metadata":{}},{"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[12], itrs)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"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","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"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()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"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()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"f, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize = (12, 12))\n\nax1.imshow(test_patient_images[12], 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()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Conclusion:\n\nIn this kernel, we discussed about:\n\n- Idea behind Watershed Transformation\n- Loading Scans\n- Getting HU values from pixel data\n- Drawback of Watershed Transforamtion\n- Moving to Marker-Controlled Watershed Transformation\n- Creating Internal, External and Watershed Markers\n- Using Sobel Gradient\n- Performing Black Top-hat Morphology\n- Getting the segmented lung images\n\nI think that there's still a long way to go in this competition. These segmented lung images can be used for determination of lung volume or for any other purpose.\n\nWatershed transformation is comparatively better than algorithms. We use a marker-based approach which reduces the chances of over-segmentation and preserves the original lung border very accurately.","metadata":{}}]}