Lung Segmentation with Machine Learning

In this article, I will introduce you to the application of Machine Learning in healthcare. I will show you how we can work on the lung segmentation task with machine learning using python.

Lung Segmentation with Machine Learning

Lung segmentation is one of the most useful tasks of machine learning in healthcare. Lung CT image segmentation is an initial step necessary for lung image analysis, it is a preliminary step to provide accurate lung CT image analysis such as detection of lung cancer.

Also, Read – Cross-Validation in Machine Learning.

Now let’s see how we can use machine learning for the lung segmentation task. Before we start, I’ll import a few packages and determine the available patients:

import numpy as np
import pandas as pd 
import dicom
import os
import scipy.ndimage
import matplotlib.pyplot as plt

from skimage import measure, morphology
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

# Some constants 
INPUT_FOLDER = 'path to sample_images'
patients = os.listdir(INPUT_FOLDER)
patients.sort()Code language: PHP (php)

Loading Data for Lung Segmentation

Dicom is the de-facto repository in medical imaging. This is my first time working with it, but it seems pretty straightforward. These files contain a lot of metadata. This analysis pixel size/coarseness differs from analysis to analysis, which can adversely affect the performance of CNN approaches.

Below is the code to load an analysis, which consists of multiple slices, which we simply save in a Python list. Each record in the dataset is an analysis. A metadata field is missing, the size of the pixels in the Z direction, which is the thickness of the slice. Fortunately, we can infer it and add it to the metadata:

def load_scan(path):
    slices = [dicom.read_file(path + '/' + s) for s in os.listdir(path)]
    slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))
    try:
        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])
    except:
        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)
        
    for s in slices:
        s.SliceThickness = slice_thickness
        
    return slicesCode language: PHP (php)

The unit of measure for CT scans is the Hounsfield Unit (HU), which is a measure of radiodensity. CT scanners are carefully calibrated to measure this accurately. From Wikipedia:

By default, the returned values are not in this unit. We first need to fix this.

Some scanners have cylindrical scan limits, but the output image is square. Pixels that go outside these limits get the fixed value -2000. The first step is to set these values to 0, which is currently air. Then back to HU units, multiplying by the rescaling slope and adding the intercept:

def get_pixels_hu(slices):
    image = np.stack([s.pixel_array for s in slices])
    # Convert to int16 (from sometimes int16), 
    # should be possible as values should always be low enough (&amp;lt;32k)
    image = image.astype(np.int16)

    # Set outside-of-scan pixels to 0
    # The intercept is usually -1024, so air is approximately 0
    image[image == -2000] = 0
    
    # Convert to Hounsfield units (HU)
    for slice_number in range(len(slices)):
        
        intercept = slices[slice_number].RescaleIntercept
        slope = slices[slice_number].RescaleSlope
        
        if slope != 1:
            image[slice_number] = slope * image[slice_number].astype(np.float64)
            image[slice_number] = image[slice_number].astype(np.int16)
            
        image[slice_number] += np.int16(intercept)
    
    return np.array(image, dtype=np.int16)Code language: PHP (php)

Now let’s take a look at one of the patients:

first_patient = load_scan(INPUT_FOLDER + patients[0])
first_patient_pixels = get_pixels_hu(first_patient)
plt.hist(first_patient_pixels.flatten(), bins=80, color='c')
plt.xlabel("Hounsfield Units (HU)")
plt.ylabel("Frequency")
plt.show()

# Show some slice in the middle
plt.imshow(first_patient_pixels[80], cmap=plt.cm.gray)
plt.show()Code language: PHP (php)

By looking at the information of Lung CT measurements from Wikipedia and the histogram above, we can see which pixels are air and which are tissue. We will use this for the lung segmentation task later.

Resampling

A CT scan normally has a pixel spacing of [2.5, 0.5, 0.5], which means that the distance between the slices is 2.5 millimetres. For a different analysis, this can be [1,5, 0,725, 0,725], it can be problematic for an automatic analysis (eg using ConvNets).

A common method of solving this problem is to resample the entire data set to a certain isotropic resolution:

def resample(image, scan, new_spacing=[1,1,1]):
    # Determine current pixel spacing
    spacing = np.array([scan[0].SliceThickness] + scan[0].PixelSpacing, dtype=np.float32)

    resize_factor = spacing / new_spacing
    new_real_shape = image.shape * resize_factor
    new_shape = np.round(new_real_shape)
    real_resize_factor = new_shape / image.shape
    new_spacing = spacing / real_resize_factor
    
    image = scipy.ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest')
    
    return image, new_spacingCode language: PHP (php)

When you apply this, to the whole dataset due to rounding, this may be slightly off from the desired spacing. We need to resample our patient’s pixels to an isomorphic resolution of 1 by 1 by 1 mm:

pix_resampled, spacing = resample(first_patient_pixels, first_patient, [1,1,1])
print("Shape before resampling\t", first_patient_pixels.shape)
print("Shape after resampling\t", pix_resampled.shape)Code language: PHP (php)

Shape before resampling	 (134, 512, 512)
Shape after resampling	 (335, 306, 306)

3D Plotting of Lungs

For visualization, it is useful to be able to display a 3D image of the scan. So we’re going to use walking cubes to create a rough mesh for our 3D object, and the plot that with matplotlib:

def plot_3d(image, threshold=-300):
    
    # Position the scan upright, 
    # so the head of the patient would be at the top facing the camera
    p = image.transpose(2,1,0)
    
    verts, faces = measure.marching_cubes(p, threshold)

    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111, projection='3d')

    # Fancy indexing: `verts[faces]` to generate a collection of triangles
    mesh = Poly3DCollection(verts[faces], alpha=0.70)
    face_color = [0.45, 0.45, 0.75]
    mesh.set_facecolor(face_color)
    ax.add_collection3d(mesh)

    ax.set_xlim(0, p.shape[0])
    ax.set_ylim(0, p.shape[1])
    ax.set_zlim(0, p.shape[2])

    plt.show()
Code language: PHP (php)

Our plot function takes an argument which we can use to plot structures, such as all tissue or only bone. 400 is a good threshold to show only bones:

plot_3d(pix_resampled, 400)

Lung segmentation

In order to reduce the problem space, we can segment the lungs (and usually certain tissues around it). The method that my fellow students and I developed was quite effective.

It involves several smart steps. It consists of a series of regional growth applications and morphological operations. In this case, we will only use the analysis of connected components.

The steps:

Image Threshold (-320 HU is a good threshold, but whatever for this approach)
Make connected components, determine the label of the air around the person, fill it with 1s in the binary image
Optional: For each axial section of the scan, determine the largest connected solid component (the body + air around the person) and set the others to 0. This fills the lung structures of the mask.
Keep only the largest air pocket (the human body has other air pockets here and there).

def largest_label_volume(im, bg=-1):
    vals, counts = np.unique(im, return_counts=True)

    counts = counts[vals != bg]
    vals = vals[vals != bg]

    if len(counts) &amp;gt; 0:
        return vals[np.argmax(counts)]
    else:
        return None

def segment_lung_mask(image, fill_lung_structures=True):
    
    # not actually binary, but 1 and 2. 
    # 0 is treated as background, which we do not want
    binary_image = np.array(image &amp;gt; -320, dtype=np.int8)+1
    labels = measure.label(binary_image)
    
    # Pick the pixel in the very corner to determine which label is air.
    #   Improvement: Pick multiple background labels from around the patient
    #   More resistant to "trays" on which the patient lays cutting the air 
    #   around the person in half
    background_label = labels[0,0,0]
    
    #Fill the air around the person
    binary_image[background_label == labels] = 2
    
    
    # Method of filling the lung structures (that is superior to something like 
    # morphological closing)
    if fill_lung_structures:
        # For every slice we determine the largest solid structure
        for i, axial_slice in enumerate(binary_image):
            axial_slice = axial_slice - 1
            labeling = measure.label(axial_slice)
            l_max = largest_label_volume(labeling, bg=0)
            
            if l_max is not None: #This slice contains some lung
                binary_image[i][labeling != l_max] = 1

    
    binary_image -= 1 #Make the image actual binary
    binary_image = 1-binary_image # Invert it, lungs are now 1
    
    # Remove other air pockets insided body
    labels = measure.label(binary_image, background=0)
    l_max = largest_label_volume(labels, bg=0)
    if l_max is not None: # There are air pockets
        binary_image[labels != l_max] = 0
 
    return binary_imageCode language: PHP (php)

But there is one thing we can fix, it’s probably a good idea to include structures in the lungs (like the nodules are solid), we don’t just want to ventilate in the lungs:

segmented_lungs = segment_lung_mask(pix_resampled, False)
segmented_lungs_fill = segment_lung_mask(pix_resampled, True)
plot_3d(segmented_lungs, 0)
Code language: PHP (php)

It’s better. Let’s also visualize the difference between the two:

plot_3d(segmented_lungs_fill - segmented_lungs, 0)

Also, Read – Data Leakage in Machine Learning.

I hope you liked this article on the Lung Segmentation as an application of Machine Learning on Healthcare. Feel free to ask your valuable questions in the comments section below. You can also follow me on Medium to learn every topic of Machine Learning.