Example a/segmentation: Segmentation of image based fields

This example demonstrates how to segment image based computed fields. It also shows how to use the results of the segmentation to generate data points to describe the boundary of an anatomical feature.

Segmentation refers to the process of partitioning an image into various regions. These regions are identified by sets of pixels with the same intensity values. Segmentation of medical images is typically used to identify the boundaries between different types of tissue and pinpoint the position of objects such as bones and organs.

The human brain is very good at segmenting images but it is quite a difficult task to automate. Many different segmentation algorithms have been developed and for a given image one approach may produce better results than another.

Cmgui implements segmentation via computed fields, using the Insight Tookit (ITK) library code. For background information on segmentation see chapter 9 of the ITK software guide. The latest version of the guide can be found at www.itk.org/ItkSoftwareGuide.pdf

Currently cmgui supports the following segmentation algorithms:

This example demonstrates the use of both methods to segment out different areas from a sample brain image

Connected threshold segmentation

The filter pipeline used for the connected threshold segmentation is as follows:
Input image
Curvature
anisotropic
diffusion
filter
Smoothed
Connected
threshold
filter
 
Output image
Input image Smoothed Segmented image

The images show the stages in segmenting out the left and right ventricle. Two seed points were used for the connected threshold filter, one within the left ventricle and one within the right.

Fast marching segmentation

The filter pipeline used for the fast marching segmentation is as follows:
Input image
Curvature
anisotropic
diffusion
filter
 
Smoothed
Gradient
magnitude
recursive
gaussian
filter
Gradient
Sigma
filter
 
 
 
Sigma
Fast
marching
filter
 
 
Time cross map
Binary
threshold
filter
 
 
Output image
Input image Smoothed Gradient Speed term Time-crossing map Segmented image

The images show the stages in segmenting out the right ventricle and the gray matter. Two seed points were used for the fast marching filter, one within the right ventricle and one within the left half of the gray matter.

Here the speed term is calculated using a sigmoid function, but it could also be calculated using the negative exponential exp(-x) or the reciprocal 1/(1+x).

Created by Peter Bier, July 2007.

Screenshot of example a/segmentation


The comfile run by this example is as follows:

# This example demonstrates the use of two methods which can be used to
# segment out different regions from a sample brain image.
# The segmentation algorithms are described in more depth in the
# ITK software guides, www.itk.org/ItkSoftwareGuide.pdf
# Once the image has been segmented, data points can be generated
# on the boundaries of the segmented regions, as shown at the
# end of this example.
#

# CONNECTED THRESHOLD SEGMENTATION

# Connected threshold segmentation is a very simple region growing method.
# One or more seed points are specified for the image and the region grows 
# from those points, adding neighbouring pixels to the region if the pixel
# intensity falls within a specific interval.  The regions grow until no 
# more pixels can be added.

# Noise present in an input image can reduce the effectivness of this filter, 
# so it is often useful to remove noise by first passing the image through
# an edge preserving smoothing filter.

# The filter pipeline is as follows (inputs to the filters shown in brackets):

#     Input image
#        |
#        V
# Curvature anisotropic diffusion filter (iterations, time step, conductance)
#        |
#        V
#     Smoothed image
#        |
#        V
# Connected threshold filter (lower threshold, upper threshold, seed points)
#        |
#        V
#     Output image

# We will use the connected threshold filter to segment out the left and
# right ventricle.

# Read in a mesh to hold image which is currently required for cmiss.
gfx read nodes $example/square.exnode
gfx read elements $example/square.exelem

# read in the image to a texture and define a field based on it
gfx create texture input_image image png:$example/BrainProtonDensitySlice.png nearest width 1 height 1;
gfx create material input_image texture input_image

# Set up a dummy image of the correct size for the output field
# This will be used to write out filtered images
gfx create texture output_image image png:$example/BrainProtonDensitySlice.png specify_format i nearest width 1 height 1;
gfx create material output_image texture output_image

# create a spectrum
gfx create spectrum gray_spectrum clear overwrite_colour
gfx modify spectrum gray_spectrum linear range 0 1 extend_above extend_below monochrome colour_range 0 1 ambient diffuse component 1

# Set up the pipeline for the connected threshold filter

# Smooth the image using anisotropic diffusion with 5 iterations
# Note that if trying to duplicate the itk segmentation tests that the 
# conductance and time step do not require dividing by 255
gfx define field smoothed curvature_anisotropic_diffusion_filter field input_image num_iterations 5 time_step 0.125 conductance 9.0;

# Segment the image using the connected threshold filter
# The seed points are defined in terms of their position within the texture
# described by normalised coordinates.
# In terms of pixel value the seed points are (85,107) and (96,107)
# The image pixel dimensions are 181x217 which gives the normallised seed 
# points as (0.47,0.5) and (0.53,0.5)
# Pixels will be added to the region if the fall within the intensity range
# 0.85 to 1.0
gfx define field connected_threshold_segmentation connected_threshold_filter field smoothed num_seed_points 2 dimension 2 seed_points 0.47 0.5 0.53 0.5 lower_threshold 0.85 upper_Threshold 1.0

# write out the segmented image to a file
gfx modify texture output_image nearest width 1 height 1 evaluate_image field connected_threshold_segmentation spectrum gray_spectrum texture_coord xi element_group square;
gfx write texture output_image file connected_threshold_segmentation.png;

# FAST MARCHING SEGMENTATION

# Fast marching segmentation is a level set segmentation method which can be 
# used if the equation governing the level set evolution is very simple. 
# Level set segmentation is a very powerful technique that uses a numerical 
# method for tracking the evolution of contours and surfaces, based on terms
# derived from an input image.

# In the case of the fast marching segmentation filter, the "speed term"
# required to calculate the segmentation is typically some function of the
# image gradient.  As with the connected threshold example, the input image is
# usually smoothed first, before computing the gradient and then some
# function of the gradient to give the speed image.  The fast marching
# filter then produces a time map, indicating how long it takes to get from
# a seed point to any given point on the image.  By using a threshold filter
# on the time map a segemented image can be obtained.

# The fast marching filter segmentation method demonstrated here uses the same 
# filters outlined in the collaboration diagram in the ITK software guide 
# (See Figure 9.15)
#
# The filter pipeline is as follows (inputs to the filters shown in brackets):

#      Input image
#         |
#         V
# Curvature anisotropic diffusion filter (iterations, time step, conductance)
#         |
#         V
#      Smoothed image
#         |
#         V
# Gradient magnitude recursive gaussian (sigma)
#         |
#         V
#      Gradient image
#         |
#         V
# Sigmoid filter (alpha, beta)
#         |
#         V
#      Speed term
#         |
#         V
# Fast marching filter (seed points, stopping value)
#         |
#         V
#      Time-crossing map
#         |
#         V
# Binary Threshold filter (lower threshold, upper threshold)
#         |
#         V
#      Segmented image


# We will use the fast marching image filter to segment out the white
# matter and the right ventricle.

# Assemble the pipeline for the fast marching image filter segmentation

# Smooth the image using anisotropic diffusion with 5 iterations
# Note that if trying to duplicate the itk segmentation tests that the 
# conductance and time step do not require dividing by 255
gfx define field smoothed curvature_anisotropic_diffusion_filter field input_image num_iterations 5 time_step 0.125 conductance 9.0;


# Calculate the gradient of the image using a sigma value of 1.0/255
gfx define field gradient gradient_magnitude_recursive_gaussian_filter field smoothed sigma 0.0039216;

# Caclulate the speed term using a sigmoid with alpha -0.3/255 and 2.0/255 
gfx define field sigma sigmoid_filter field gradient alpha -0.0011765 beta  0.0078431;

# Calculate a time-crossing map indicating how long it takes to reach
# pixels from the starting seed points, using the fast marching image filter
# The seed points are defined in terms of their position within the texture
# described by normalised coordinates.
# In terms of pixel value the seed points are (56,92) and (96,107)
# The image pixel dimensions are 181x217 which gives the normallised seed 
# points as (0.30939,0.423396) and (0.53,0.5)

gfx define field time_cross_map fast_marching_filter field sigma num_seed_points  2 seed_points 0.31 0.42 0.53 0.5 stopping_value 1000000
 
# Caclulate the segmented image by threshold the image.  
# Pixels with intensities between 0 (black) and 0.95 (almost white) 
# will be included in the regions
gfx define field fast_marching_segmentation binary_threshold_filter field time_cross_map lower_threshold 0 upper_threshold 0.95

# write out the output images for each stage of filtering to check that 
# the filter pipeline has worked

gfx modify texture output_image nearest width 1 height 1 evaluate_image field smoothed spectrum gray_spectrum texture_coord xi element_group square;
gfx write texture output_image file smoothed.png;

gfx modify texture output_image nearest width 1 height 1 evaluate_image field gradient spectrum gray_spectrum texture_coord xi element_group square;
gfx write texture output_image file gradient.png;

gfx modify texture output_image nearest width 1 height 1 evaluate_image field sigma spectrum gray_spectrum texture_coord xi element_group square;
gfx write texture output_image file sigma.png;

gfx modify texture output_image nearest width 1 height 1 evaluate_image field time_cross_map spectrum gray_spectrum texture_coord xi element_group square;
gfx write texture output_image file time_cross_map.png;

gfx modify texture output_image nearest width 1 height 1 evaluate_image field fast_marching_segmentation spectrum gray_spectrum texture_coord xi element_group square;
gfx write texture output_image file fast_marching_segmentation.png;


# generate node points around the segmented image

# Use a fine discretization with for generating the isolines.
gfx modify g_element square general element_discretization "512*512*512"

# Define the density field used to control the number of points generated 
# Here we want a constant number of points per segment length.
gfx define field point_density constant 50;

# Create the isolines
gfx modify g_element square iso_surfaces as isolines select_on iso_scalar fast_marching_segmentation iso_value 0.5 use_elements select_on material red selected_material default_selected render_shaded data point_density;


if (!$TESTING)
  {
	# set up a window to view the output of the fast marching segmentation

	gfx create window
	gfx modify window 1 set perturb_lines;
  gfx modify window 1 layout 2d
  gfx mod g_element square surfaces coordinate coordinates select_on material output_image texture_coordinates xi selected_material default_selected render_shaded

  }

# Generate points on the isolines calculated above.
gfx create ngroup data;
gfx define graphics_filter isolines match_graphics_name isolines;
gfx convert graphics filter isolines render_surface_node_cloud region data coordinate coordinates;

# Write out these points
gfx write nodes group data data_cloud.exnode;

if (!$TESTING)
  {
	 gfx modify g_element data node_points glyph sphere general size ".05*.05*.05" centre 0,0,0 font default select_on material green selected_material default_selected;
  }

Files used by this example are:

Name                         Modified     Size

segmentation.com 17-Mar-2014 9.6k BrainProtonDensitySlice.png 17-Mar-2014 33k COPYRIGHT 17-Mar-2014 504 square.exelem 17-Mar-2014 1.5k square.exnode 17-Mar-2014 400

Download the entire example:

Name                            Modified     Size

examples_a_segmentation.tar.gz 09-Mar-2016 420k

Testing status by version:

StatusTestedReal time (s)
i686-linux
cmgui-wxFailureSun Mar 6 00:16:08 20164
last breakTue Feb 24 03:18:00 20153
cmgui-wx-debugFailureSun Mar 6 00:15:34 20165
last breakTue Feb 24 03:19:00 20154
cmgui-wx-debug-memorycheckFailureSun Mar 6 00:19:10 20165
last breakTue Feb 24 03:21:00 20155
cmgui-wx-debug-valgrindFailureSun Mar 6 01:34:09 2016133
last breakSun Mar 6 01:31:00 2016133
x86_64-linux
cmgui-wxFailureSun Mar 6 00:01:32 20160
last breakSun Mar 6 00:01:00 20160
cmgui-wx-debugFailureSun Mar 6 00:01:33 20161
last breakSun Mar 6 00:01:00 20161
cmgui-wx-debug-memorycheckFailureSun Mar 6 00:01:33 20161
last breakSun Mar 6 00:01:00 20161
cmgui-wx-debug-valgrindFailureSun Mar 6 00:02:56 20169
last breakSun Mar 6 00:02:00 20169

Testing status by file:


Html last generated: Wed Mar 9 16:02:36 2016

Input last modified: Wed Mar 9 15:49:40 2016


CMISS Help / Examples / a / segmentation