A Matlab toolbox implementing Level Set Methods

Baris Sumengen, Vision Research Lab at UC Santa Barbara

• Download
• Purpose
• Capabilities
• Tutorial
• Examples

This set of Matlab files implements Level Set Methods and follows Osher and Fedkiw's book. A combination of curvature-based forces, vector field-based forces and forces in  the normal direction can be used.

Download: Download by clicking here. This package consist of a set of m-files compressed into a zip file. Unzip them to a folder and then add this folder to Matlab's path.
Any type of comments, bug reports, benchmarks, profiling efforts are highly appreciated. Please send your comments to  .

- I would be happy if you acknowledge this toolbox in your papers if you end up using it for your research.

Note: See also Ian Mitchell's Matlab toolbox for Level Set Methods for a different implementation of the same problem.
Note 2: See also  Gabriel Peyré's Matlab toolbox for Fast Marching Methods. This uses a C implementation and .mex files. Currently fast marching methods are not implemented in my toolbox but will be soon.

Purpose: Creating a simple to use and understand implementation of  the Level Set Methods. Speed was not the main motivating point as much as readability of the code. On the other hand, given Matlab's capabilities, this toolbox should be reasonably fast (unintuitive Matlab tricks are avoided intentionally to prevent confusion). It  can be a good learning tool to accompany Osher and Fedkiw's book on Level Set Methods. Also it would be a good template for a C or Java based implementation.
This toolbox in its current state can be used for implementing a 2-D curve evolution or a diffusion of a 2-D function phi(x,y), e.g. anisotropic diffusion on a gray-scale image.

Capabilities:

• This toolbox currently implements only a 2-D interface (curve) evolution. Extending to 3-D interface (surface) evolution is trivial because of its modular structure based on 1-D functions (I'll do this soon).
• 1st, 2nd and 3rd order accurate Essentially Non-Oscillatory (ENO) and 5th order accurate Weighted ENO (WENO) are supported throughout the toolbox for derivative calculations.
• A function for re-initialization of a signed distance function is provided (reinit_SD).
• A high-level function (evolve2D) is provided that can utilize any combination of: 1) a force in the normal direction to the curve, 2) a curvature-based force, 3) evolution under the influence of an external vector field.

Tutorial

Also see my power point presentation: Overview of Level Set Methods and Image Segmentation

There are three m-files that needs to be talked about first: reinit_SD(), evolve2D() and test.m

test.m include a fully working program that implements curve evolution: starts with a square shaped curve and shrinks it under constant normal force.

The function evolve2D() is a high level function that takes an input, evolves it N iterations and returns the result. Here is the prototype and documentation for this function (taken from the m-file):

The function reinit_SD() can be used for re-initialization of a signed distance function. This function would be especially effective if the input is not too far away from a signed distance function.

evolve2D() does not call reinit_SD() within itself (and as a result may not preserve signed distance-ness of phi). On the other hand, if desired, in an outside loop you can call evolve2D() for 5-10 iterations followed by a call to reinit_SD(), then call evolve2D() again for 5-10 more iterations, etc.

Examples

Let us first start with a rectangular curve. This curve can be defined as the boundary of the following binary function (or image, matrix) consisting of 1's and -1's:

phi = ones(100,100);
phi(30:60,30:60) = -1;

We can think of the curve as the zero level set of this function (in continuous domain): phi(x,y) == 0. As can be seen, any function phi(x,y) whose zero level set corresponds to this curve can be a representation of this curve. Since a signed distance function will make life easier for us due to its smoothness, we would like to create a function phi(x,y) such that |phi(x,y)| = d, where d is the distance of (x,y) from the curve itself. The sign is taken as (+) if (x,y) is outside the curve, (-) if (x,y) is inside the curve. One way to create this signed distance function is calling reinit_SD(), which will smooth phi using a reinitialization PDE (see equations 7.4 and 7.5 in Osher&Fedkiw's book).

phi = reinit_SD(phi, dx, dy, 0.5, 'WENO', 200);

The following images show the progress of this function. As can be seen, the corners of the curve are slightly smoothed but except that the zero level set of phi is preserved during the diffusion.

Evolution of a binary function towards a signed distance function under re-initialization PDE

phi after 0 iterations	curve (zero level set) after 0 iterations
phi after 30 iterations	curve (zero level set) after 30 iterations
phi after 200 iterations	curve (zero level set) after 200 iterations

Since this process takes longer than necessary and smoothes the curve slightly, we'll use for the rest of this discussion a different way of creating this signed distance function. We will utilize Matlab's bwdist() function.

phi = ones(100,100);
phi(30:60,30:60) = -1;
phi = double((phi > 0).*(bwdist(phi < 0)-0.5) - (phi < 0).*(bwdist(phi > 0)-0.5));
figure;surf(phi);
figure;contour(phi, [0 0], 'r');

Now we can start evolving this initial curve using the evolve2D() function. Note that, phi does not have to be a signed distance function, but if it is this increases the stability and quality of the evolution (especially if a vector field-based force and a force in normal direction are combined).

1) Shrinking the curve in the normal direction:

% Vn: force in the normal direction
Vn = -0.2*ones(100,100);

dx=1;
dy=1;

figure;contour(phi, [0 0], 'b');
phi = evolve2D(phi,dx,dy,0.5,15,'ENO2',[],1,Vn,0,[],[],0,[]);
hold on;contour(phi, [0 0], 'r');hold off;axis equal;

2) Translating the curve using an external vector field

% External vector field
u = -0.5*(1/sqrt(2)).*ones(100,100); %x component of the vector field
v = -0.5*(1/sqrt(2)).*ones(100,100); %y component of the vector field

dx=1;
dy=1;

figure;contour(phi,[0 0],'b');
phi = evolve2D(phi,dx,dy,0.5,25,'ENO3',[],0,[],1,u,v,0,[]);
hold on;contour(phi,[0 0],'r'); hold off; axis equal;

3) Smoothing a curve under its curvature

% b: weighting for curvature-based force. b needs to be positive
b = 0.3*ones(100, 100);

dx=1;
dy=1;

figure;contour(phi,[0 0],'b');
phi = evolve2D(phi,dx,dy,0.5,25,[],[],0,[],0,[],[],1,b);
hold on; contour(phi,[0 0],'r'); hold off; axis equal;

As can be seen, curvature of a straight line is 0, so only the corners are smoothed.

4) Simultanous translation and expansion of the curve

% External vector field
u = -0.5*(1/sqrt(2)).*ones(100, 100); %x component of the vector field
v = -0.5*(1/sqrt(2)).*ones(100, 100); %y component of the vector field

% Vn: force in the normal direction
Vn = 0.2*ones(100, 100);

dx=1;
dy=1;

figure;contour(phi,[0 0],'b');
phi = evolve2D(phi,dx,dy,0.5,60,'ENO2',1,1,Vn,1,u,v,0,[]);
hold on; contour(phi,[0 0],'r'); hold off; axis equal;
 

5) Smoothing an image under the curvature of its level sets (Mean curvature flow)

In this example, all level sets of the image are smoothed under curvature-based forces.
Note that Matlab's contour function selects the levels automatically, so the level sets displayed might not correspond to each other.

dx=1;
dy=1;

im = imread('retina.bmp');
phi = double(im(:,:,1));
figure; imagesc(phi); axis image; colormap gray;
figure; contour(phi); axis image;

% b:weighting for curvature-based force. b needs to be positive
b = 0.3*ones(size(phi));

phi = evolve2D(phi,dx,dy,0.5,25,[],[],0,[],0,[],[],1,b);
figure; imagesc(phi); axis image; colormap gray;
figure; contour(phi); axis image;

Original Image

Level Sets of the original image

After 5 iterations of smoothing

Corresponding level sets (5 iterations)

After 25 iterations of smoothing

Corresponding level sets (25 iterations)