Part 3: Segmentation of an aorta

In the previous notebook, we presented the most fundamental operations of image processing. Let us now showcase their use on a simple toy problem: the segmentation of an aorta from axial slices.

As always, we first need to re-import the necessary modules and define our custom display routine:

In [1]:
%matplotlib inline
import center_images             # Center our images
import matplotlib.pyplot as plt  # Graphical display
import numpy as np               # Numerical computations
from imageio import imread       # Load .png and .jpg images

def display(im):  # Define a new Python routine
    """
    Displays an image using the methods of the 'matplotlib' library.
    """
    plt.figure(figsize=(8,8))                    # Create a square blackboard
    plt.imshow(im, cmap="gray", vmin=0, vmax=1)  # Display 'im' using gray colors,
                                                 #     from 0 (black) to 1 (white)

We will showcase our method on images kindly provided by Tom Boeken from the HEGP hospital - strictly following legal requirements and consent:

In [2]:
# Import our images as grayscale, normalized from 0 (black) to 1 (white)
Images = [ imread("data/aortes/{}.jpg".format(i), as_gray=True) / 255 
           for i in [1,2,3,4,5,6] ]
# 'Images' is a list of six 2d arrays, denoted by 'Images[0]', ..., 'Images[5]'

plt.figure(figsize=(15,10)) # Rectangular blackboard
for i in [1,2,3,4,5,6] :    # Loop on indices i from 1 to 6
    plt.subplot(2,3,i)      # Sub-plot "i" in a 2x3 waffle plot
    plt.imshow( Images[i-1], cmap="gray", vmin=0, vmax=1) # Images[0] is top-left, etc.
    plt.title("Images[{}]".format(i-1))

First, we focus on the top-left slice and re-name it I:

In [3]:
I = Images[0]
display(I)

How can we extract the aorta from the surrounding background? An acute human eye performs this "trivial" pre-processing step in a fraction of a second... But our computer cannot benefit from the machinery hardcoded in the human brain! One way or another, we must find a way to decompose this operation into simple, basic lines of code.

Fortunately, in this toy example, the color of the aorta is clearly different from that of the other tissues. Using a simple threshold, we should thus be able to make a huge step in the right direction:

In [4]:
mask = (.6 < I) & (I < .9)  # Only keep pixel whose intensity lies in the (.6,.9) range
display(mask)

Great! Our aorta is now the only structure that is both thick and bright. To extract it, we will use simply operations from mathematical morphology:

In [5]:
from scipy.ndimage.morphology import binary_erosion  # "Erode" white structures
J = binary_erosion(mask)  # by removing pixels which have a black neighbour
display(J)
In [6]:
# Remove pixels which are close to the black background:
core = binary_erosion(mask, iterations=5)  # Iterate "erosion" five times
display(core)

A succession of thresholdings and erosions has allowed us to locate the core of our aorta. We can now let it grow to delimit a reasonable trust region:

In [7]:
from scipy.ndimage.morphology import binary_dilation  # Invert "erosion":
region = binary_dilation(core, iterations=10)  # Let the "white region" grow 10 times
display(region)

And, finally, intersect this mask with the set of pixels "which are bright enough":

In [8]:
aorta  = (I > .6) & region  # Intersection between "image is not too dark" 
                            #                        and the trust region

# Fancy display in a 3x2 waffle grid ---------------------------------------------
plt.figure(figsize=(16,24))
plt.subplot(3,2,1); plt.imshow(I,         cmap="gray"); plt.title("Raw data")
plt.subplot(3,2,2); plt.imshow(mask,      cmap="gray"); plt.title("Thresholded")
plt.subplot(3,2,3); plt.imshow(core,      cmap="gray"); plt.title("Eroded")
plt.subplot(3,2,4); plt.imshow(region,    cmap="gray"); plt.title("Inflated")
plt.subplot(3,2,5); plt.imshow(aorta,     cmap="gray"); plt.title("Intersection with the bright pixels")
plt.subplot(3,2,6); plt.imshow(aorta+.5*I,cmap="gray"); plt.title("Segmented aorta");

That's it! Using super-simple operations provided by the scipy library, we segmented our aorta. We can package our routine as a Python function, and provide users with a new keyword segment:

In [9]:
def segment(I, vmin=.6, vmax=.9, erosion=5, dilation=15) :
    """
    Tries to segment an aorta using the default parameters presented above.
    """
    mask   = (vmin < I) & (I < vmax)                     # Threshold the image values
    core   = binary_erosion( mask, iterations=erosion)   # Remove the thin components 
    region = binary_dilation(core, iterations=dilation)  # Core growth -> trust region
    aorta  = (I > .6) * region   # Intersection between "image is not too dark" 
                                 #                        and the trust region
    
    # Fancy display in a 3x2 waffle grid ---------------------------------------------
    plt.figure(figsize=(16,24))
    plt.subplot(3,2,1); plt.imshow(I,         cmap="gray"); plt.title("Raw data")
    plt.subplot(3,2,2); plt.imshow(mask,      cmap="gray"); plt.title("Thresholded")
    plt.subplot(3,2,3); plt.imshow(core,      cmap="gray"); plt.title("Eroded")
    plt.subplot(3,2,4); plt.imshow(region,    cmap="gray"); plt.title("Inflated")
    plt.subplot(3,2,5); plt.imshow(aorta,     cmap="gray"); plt.title("Intersection with the bright pixels")
    plt.subplot(3,2,6); plt.imshow(aorta+.5*I,cmap="gray"); plt.title("Segmented aorta");

Exercise : Let us display our "dataset" of axial slices once again:

In [10]:
plt.figure(figsize=(15,10)) # Rectangular blackboard
for i in [1,2,3,4,5,6] :    # Loop on indices i from 1 to 6
    plt.subplot(2,3,i)      # Sub-plot "i" in a 2x3 waffle plot
    plt.imshow( Images[i-1], cmap="gray", vmin=0, vmax=1) # Images[0] is top-left, etc.
    plt.title("Images[{}]".format(i-1))

What can you say about them? As far as you can tell, what are the images that should be "easy to handle"?

Playing with the parameters of our brand new program, try to segment the five remaining aortas Images[1], ..., Images[5]. Comment.

In [11]:
segment(Images[1], dilation=15)
In [12]:
segment(Images[2], erosion=7)
In [13]:
segment(Images[3])
In [14]:
segment(Images[4], erosion=7, vmax=.85)
In [15]:
segment(Images[5], erosion=7, dilation=15, vmax=.85)

Solution + to remember: "AI" or not, all computer programs rely on simple operations defined on numerical arrays. In specific settings (e.g. when one color stands out), some problems may be surprisingly easy to solve. But performances plummet as soon as programs leave their comfort zones: Images[3] and Images[4] are nearly impossible to segment with the simple tools presented in this notebook.

Here, the fact that kidneys are opacified by the constrast agent just as well as the aorta (a characteristic of "healthy", well-vascularized patients) prevents our simplistic color thresholding from working as intended: the aorta's color is too similar to that of the kidneys and vertebrae. CT-scan reconstruction artifacts worsen the problem, as they induce a "bubbling" noise inside the aorta that completely messes up the erosion step... As the simplistic hypotheses encoded in our program became invalid, we ended up with crappy results.

Today, Convolutional "Neural Networks" allow us to segment biomedical images by leveraging finely tuned texture descriptors, fitted on expert-annotated databases. This is a huge breakthrough, that will probably impact your clinical practice. But you should never fall in the traps of "magical thinking" and "Pygmalion syndrome": success in a specific task cannot be extrapolated to other fields - with the exception of sci-fi novels, of course.

Keep in mind that chess, the king of Western games, was "cracked" in 1997!