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.
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()
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)) try: slice_thickness = np.abs(slices.ImagePositionPatient - slices.ImagePositionPatient) except: slice_thickness = np.abs(slices.SliceLocation - slices.SliceLocation) for s in slices: s.SliceThickness = slice_thickness return slices
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 (<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)
Now let’s take a look at one of the patients:
first_patient = load_scan(INPUT_FOLDER + patients) 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, cmap=plt.cm.gray) plt.show()
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.
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.SliceThickness] + scan.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_spacing
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)
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) ax.set_ylim(0, p.shape) ax.set_zlim(0, p.shape) plt.show()
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:
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.
- 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) > 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 > -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_image
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)
It’s better. Let’s also visualize the difference between the two:
plot_3d(segmented_lungs_fill - segmented_lungs, 0)
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.