Level set segmentation on GPUs using OpenCL

Brain segmented from synthetic MR images (generated at BrainWeb) on the GPU using OpenCL and the Level Set method

The actual contour, is defined by the zero level set, which are the coordinates (x,y,z) where the level set function is zero:

To move the contour, the level set function is derivated in respect to time:

F is called the speed function and defines how fast and in which direction the contour moves. The speed function can be tailored for any problem. In image segmentation it is usual to model the speed function to be high at coordinates where the image has a desired intensity and visa versa. To make the contour smooth and avoid leaking into surrounding regions a curvature term (\(\kappa = \nabla \cdot \frac{\nabla \phi}{|\nabla \phi|}\)) is often included in the speed function. A popular choice of speed function for image segmentation is:

Here \(\alpha \in [0,1]\) is a weighting parameter between the intensity and the curvature term. The parameters T and \(\epsilon\) are used to drive to contour toward voxels with intensity in the range \(I \in [T-\epsilon,T+\epsilon]\).

Level set surface moving in the image plane. The red circles show the zero level set at various time steps. As time goes, the surface is moved down through the image plane and the zero level set change according to the shape of the surface.

The level set method is very computationally expensive because each voxel has to be updated for each iteration. However, each voxel can be updated in parallel using the same instructions, making level sets ideal for GPUs (see [2,3,4] for details on different GPU implementations). I have created a simple GPU accelerated version of level set volume segmentation using OpenCL. The implementation uses 3D textures on the GPU to reduce memory access latency. Read more on textures in OpenCL my previous post on Gaussian Blur using OpenCL. If you want to look into further optimizing the level set computation you should look into the narrow band, sparse field or fast marching methods (see [1] for more details).

The level set gradient \(\nabla \phi\) and the curvature \(\kappa\) has to approximated numerically. This can be done using the upwinding scheme.

The level set function has to be initialized. It is common to initialize it to the distance transform which calculates the distance from each voxel to the initial contour. The signed distance is negative for voxels inside the initial contour and positive outside. If we use a spherical initial contour the signed distance transform can be easily calculated in parallel for each voxel using the following equation \(d = |\vec x – \vec c| – r\) where \(\vec x\) is the coordinate of the voxel, \(\vec c\) is the position of the center and r is the radius.

Download and run the example

The code is available on GitHub for download.

The program uses the Simple Image Processing Library (SIPL) for loading, storing and displaying the volumes. This library is dependent on GTK 2.

Below are a set of commands for downloading, compiling and running the example on Ubuntu.

# Install dependencies (OpenCL has to be installed manually)
sudo apt-get install libgtk2.0-dev
 
# Download
git clone git://github.com/smistad/OpenCL-Level-Set-Segmentation.git
cd OpenCL-Level-Set-Segmentation
git submodule init
git submodule update
 
# Compile and run
cmake .
make
./levelSetSeg example_data/mr_brain.mhd result.mhd 100 100 100 10 2000 125 40 0.05 125 255

References

1. Level Set Methods and Fast Marching Methods by J.A. Sethian. Cambridge University Press
2. Rumpf, M., Strzodka, R. Level set segmentation in graphics hardware. Proceedings 2001 International Conference on Image Processing 1103–1106
3. Lefohn, A., Cates, J., & Whitaker, R. . Interactive, gpu-based level sets for 3d segmentation. Medical Image Computing and Computer-Assisted Intervention – MICCAI 2003. 564–572
4. Roberts, M., Packer, J., Sousa, M. C., & Mitchell, J. R. (2010). A Work-Efficient GPU Algorithm for Level Set Segmentation. Proceedings of the Conference on High Performance Graphics. 123–132.
5. BrainWeb. http://brainweb.bic.mni.mcgill.ca/brainweb/