Segmentation: Region Growing

In this notebook we use one of the simplest segmentation approaches, region growing. We illustrate the use of three variants of this family of algorithms. The common theme for all algorithms is that a voxel's neighbor is considered to be in the same class if its intensities are similar to the current voxel. The definition of similar is what varies:

  • ConnectedThreshold: The neighboring voxel's intensity is within explicitly specified thresholds.
  • ConfidenceConnected: The neighboring voxel's intensity is within the implicitly specified bounds $\mu\pm c\sigma$, where $\mu$ is the mean intensity of the seed points, $\sigma$ their standard deviation and $c$ a user specified constant.
  • VectorConfidenceConnected: A generalization of the previous approach to vector valued images, for instance multi-spectral images or multi-parametric MRI. The neighboring voxel's intensity vector is within the implicitly specified bounds using the Mahalanobis distance $\sqrt{(\mathbf{x}-\mathbf{\mu})^T\Sigma^{-1}(\mathbf{x}-\mathbf{\mu})}<c$, where $\mathbf{\mu}$ is the mean of the vectors at the seed points, $\Sigma$ is the covariance matrix and $c$ is a user specified constant.

We will illustrate the usage of these three filters using a cranial MRI scan (T1 and T2) and attempt to segment one of the ventricles.

In [1]:
# To use interactive plots (mouse clicks, zooming, panning) we use the notebook back end. We want our graphs 
# to be embedded in the notebook, inline mode, this combination is defined by the magic "%matplotlib notebook".
%matplotlib notebook

import SimpleITK as sitk

%run update_path_to_download_script
from downloaddata import fetch_data as fdata
import gui

# Using an external viewer (ITK-SNAP or 3D Slicer) we identified a visually appealing window-level setting
T1_WINDOW_LEVEL = (1050,500)

Read Data and Select Seed Point(s)

We first load a T1 MRI brain scan and select our seed point(s). If you are unfamiliar with the anatomy you can use the preselected seed point specified below, just uncomment the line.

In [2]:
img_T1 = sitk.ReadImage(fdata("nac-hncma-atlas2013-Slicer4Version/Data/A1_grayT1.nrrd"))
# Rescale the intensities and map them to [0,255], these are the default values for the output
# We will use this image to display the results of segmentation
img_T1_255 = sitk.Cast(sitk.IntensityWindowing(img_T1, 
                                               windowMinimum=T1_WINDOW_LEVEL[1]-T1_WINDOW_LEVEL[0]/2.0, 
                                               windowMaximum=T1_WINDOW_LEVEL[1]+T1_WINDOW_LEVEL[0]/2.0), 
                       sitk.sitkUInt8)

Fetching nac-hncma-atlas2013-Slicer4Version/Data/A1_grayT1.nrrd
In [3]:
point_acquisition_interface = gui.PointDataAquisition(img_T1, window_level=(1050,500))

#preselected seed point in the left ventricle  
point_acquisition_interface.set_point_indexes([(132,142,96)])
In [4]:
initial_seed_point_indexes = point_acquisition_interface.get_point_indexes()

ConnectedThreshold

We start by using explicitly specified thresholds, you should modify these (lower/upper) to see the effects on the resulting segmentation.

In [5]:
seg_explicit_thresholds = sitk.ConnectedThreshold(img_T1, seedList=initial_seed_point_indexes, lower=100, upper=170)
# Overlay the segmentation onto the T1 image
gui.MultiImageDisplay(image_list = [sitk.LabelOverlay(img_T1_255, seg_explicit_thresholds)],                   
                      title_list = ['connected threshold result'])
Out[5]:
<gui.MultiImageDisplay at 0x12f6281d0>

ConfidenceConnected

This region growing algorithm allows the user to implicitly specify the threshold bounds based on the statistics estimated from the seed points, $\mu\pm c\sigma$. This algorithm has some flexibility which you should familiarize yourself with:

  • The "multiplier" parameter is the constant $c$ from the formula above.
  • You can specify a region around each seed point "initialNeighborhoodRadius" from which the statistics are estimated, see what happens when you set it to zero.
  • The "numberOfIterations" allows you to rerun the algorithm. In the first run the bounds are defined by the seed voxels you specified, in the following iterations $\mu$ and $\sigma$ are estimated from the segmented points and the region growing is updated accordingly.
In [6]:
seg_implicit_thresholds = sitk.ConfidenceConnected(img_T1, seedList=initial_seed_point_indexes,
                                                   numberOfIterations=0,
                                                   multiplier=2,
                                                   initialNeighborhoodRadius=1,
                                                   replaceValue=1)

gui.MultiImageDisplay(image_list = [sitk.LabelOverlay(img_T1_255, seg_implicit_thresholds)],                   
                      title_list = ['confidence connected result'])
Out[6]:
<gui.MultiImageDisplay at 0x1352c2400>

VectorConfidenceConnected

We first load a T2 image from the same person and combine it with the T1 image to create a vector image. This region growing algorithm is similar to the previous one, ConfidenceConnected, and allows the user to implicitly specify the threshold bounds based on the statistics estimated from the seed points. The main difference is that in this case we are using the Mahalanobis and not the intensity difference.

In [7]:
img_T2 = sitk.ReadImage(fdata("nac-hncma-atlas2013-Slicer4Version/Data/A1_grayT2.nrrd"))
img_multi = sitk.Compose(img_T1, img_T2)

Fetching nac-hncma-atlas2013-Slicer4Version/Data/A1_grayT2.nrrd
In [8]:
seg_implicit_threshold_vector = sitk.VectorConfidenceConnected(img_multi, 
                                                               initial_seed_point_indexes, 
                                                               numberOfIterations=2, 
                                                               multiplier=4)

gui.MultiImageDisplay(image_list = [sitk.LabelOverlay(img_T1_255, seg_implicit_threshold_vector)],                   
                      title_list = ['vector confidence connected result'])
Out[8]:
<gui.MultiImageDisplay at 0x1352df4a8>

Clean up, Clean up...

Use of low level segmentation algorithms such as region growing is often followed by a clean up step. In this step we fill holes and remove small connected components. Both of these operations are achieved by using binary morphological operations, opening (BinaryMorphologicalOpening) to remove small connected components and closing (BinaryMorphologicalClosing) to fill holes.

SimpleITK supports several shapes for the structuring elements (kernels) including:

  • sitkAnnulus
  • sitkBall
  • sitkBox
  • sitkCross

The size of the kernel can be specified as a scalar (same for all dimensions) or as a vector of values, size per dimension.

The following code cell illustrates the results of such a clean up, using closing to remove holes in the original segmentation.

In [9]:
vectorRadius=(1,1,1)
kernel=sitk.sitkBall
seg_implicit_thresholds_clean = sitk.BinaryMorphologicalClosing(seg_implicit_thresholds, 
                                                                vectorRadius,
                                                                kernel)

And now we compare the original segmentation to the segmentation after clean up (using the GUI you can zoom in on the region of interest for a closer look).

In [10]:
gui.MultiImageDisplay(image_list = [sitk.LabelOverlay(img_T1_255, seg_implicit_thresholds), 
                                sitk.LabelOverlay(img_T1_255, seg_implicit_thresholds_clean)], 
                  shared_slider=True,
                  title_list = ['before morphological closing', 'after morphological closing'])
Out[10]:
<gui.MultiImageDisplay at 0x13a368b00>
In [ ]: