# 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

Level sets is a mathematical method of evolving contours in Cartesian grids such as images. The method works by considering a function $$\phi$$, called the level set function, which has one more dimension than the Cartesian grid we want to evolve the contour on. Thus, for a 2D image the level set function defines a 3D surface, while for a 3D volume the level set function is a 4D hypersurface. For each point on the grid (x, y, z), it defines the height h from the surface to the grid at a given time t: $$h = \phi(x,y,z,t)$$

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

$$\phi(x,y,z,t) = 0$$

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

$$\frac{\partial \phi}{\partial t} = -F|\nabla \phi|$$

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:

$$F = -\alpha (\epsilon – |T – I(x,y,z)|) + (1-\alpha)\kappa(x,y,z)$$

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.

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/

### 4 Responses

1. Thyago Maia says:

Hi, Erik

When I try to execute the same code in the CPU (Intel), the program fails creating the context and launches this error: Failed to create an OpenCL context!

I have already installed the Intel OpenCL SDK and I’ve done the following change (main.cpp, Line 227): ocl.context = createCLContext ();

Do I need to do something more?

Thanks a lot for your help!

• Erik Smistad says:

Hi!

Most often this happens because OpenCL is not properly installed. If you are on linux, make sure that the .icd file is located in /etc/OpenCL/vendors, otherwise try to reinstall or at least verify that your installation works.

• Thyago Maia says:

Hi, Erik! Thanks for the quick response!

The intel64.icd file is in the folder (I’m on Ubuntu). The example of the vector sum (available on this blog) worked normally. The error only occurs on this project.

• Thyago Maia says:

Another detail: I removed the “typedef struct OpenCL” declared in the main.cpp file, due an error in make (Error: previous definition of typedef struct OpenCL in file OpenCLUtilities.hpp)