Abstract :The goals of this lecture is to use the level set framework in order to do curve evolution. The mean curvature motion is the basic tool, and it can be extended into edge-based (geodesic active contours) and region-based (Chan-Vese) snakes.

- First download the Matlab toolbox toolbox_fast_marching.zip. Unzip it into your working directory. You should have a directory toolbox_fast_marching/ in your path.
- The first thing to do is to install this toolbox in your path.
- Recompile the mex file for your machine (this can produce some warning). If it does not work, either use the already compiled mex file (they should be available in toolbox_fast_marching/ for MacOs and Unix) or try to set up matlab with a C compiler (e.g. gcc) using 'mex -setup'.

path(path, 'toolbox_fast_marching/');

path(path, 'toolbox_fast_marching/data/');

path(path, 'toolbox_fast_marching/toolbox/');

cd toolbox_fast_marching

compile_mex;

cd ..

- In order to perform curve evolution, we will deal with a distance function stored in a 2D image D. The curve will be embedded in the level set D=0. This curve will be evolved by modifying the image D. A curve evolution ODE can be replaced by an PDE on D. This allows to deal with topological changes when a curve split or two curves merge.
- During the curve evolution, the image D might become far from being a distance function. In order to stabilize the algorithm, one needs to re-compute this distance function.

n = 200; % size of the image

% load a distance function

D0 = compute_levelset_shape('circlerect2', n);

% type 'help compute_levelset_shape' to see other

% basic curve you can load.

% display the curve

clf; hold on;

imagesc(D); axis image; axis off; axis([1 n 1 n]);

[c,h] = contour(D,[0 0], 'r');

set(h, 'LineWidth', 2);

hold off;

colormap gray(256);

% do the union of two curves

options.center = [0.15 0.15]*n;

options.radius = 0.1*n;

D1 = compute_levelset_shape('circle', n,options);

imagesc(min(D0,D1)<0);

% here we simulate a modification of the distance function

[Y,X] = meshgrid(1:n,1:n);

D = (D0.^3) .* (X+n/3);

D1 = perform_redistancing(D);

% display both the original and the new,

% redistanced, curve (should be very close)

...

- In order to compute differential quantities (tangent, normal, curvature, etc) on the curve, you can compute derivatives of the image D.
- The mean curvature motion of the level sets of some image u is driven be
the following equation.

Implement this evolution explicitly in time using finite differences.

% the gradient

g0 = divgrad(D);

% display the gradient (as arrow field with 'quiver', ...)

...

% the normalized gradient

d = max(eps, sqrt(sum(g0.^2,3)) );

g = g0 ./ repmat( d, [1 1 2] );

% display

...

% the curvature

K = d .* divgrad( g );

% display

...

Tmax = 1000; % maximum time of evolution

dt = 0.4; % time step (should be small)

niter = round(Tmax/dt); % number of iterations

D = D0; % initialization

for i=1:niter

% compute the right hand size of the PDE

...

% update the distance field

D = ...;

% redistance the function from time to time

if mod(i,30)==0

D = perform_redistancing(D);

end

% display from time to time

if mod(i,30)==1

% display here

...

end

end

Curve evolution under the mean curvature motion (the background is the distance function D). |

- Given a background image M to segment, one needs to compute an edge-stopping function E. It should be small in area of high gradient, and high in area of large gradient.
- The geodesic active contour evolution is a mean curvature motion modulated
by the edge-stopping function:

Implement this evolution explicitly in time using finite differences [warning: the g in the formula is the edge stopping function E].

% load an image

name = 'brain';

M = rescale( sum( load_image(name, n), 3) );

% display it

...

% compute a smoothed gradient

sigma = 4; % blurring size

G = divgrad( perform_blurring(M,sigma) );

% compute the norm of the gradient

d = ...

% compute the edge-stopping function

E = 1./max(d, 1e-3);

% rescale it so that it is in realistic ranges

E = rescale(E,0.3,1);

% initialize the level set

D0 = compute_levelset_shape('square', n);

% precompute the gradient of E [which is g in the formula]

...

% do the iterations

...

Segmentation with geodesic active contours. |

- The geodesic active contour uses an edge-based energy. It has lots of local
minima and is very sensitive to initialization. In order to circumvent these
drawbacks, one can use a region based energy like the Mumford-Shah functional.
Re-casted into the level set framework, it reads:

The corresponding gradient descent is the Chan-Vese active contour method:

Implement this evolution explicitly in time using finite differences, when c1 and c2 are known in advanced. - In the case that one does not know in advance the constants c0 and c1, how to update them automatically during the evolution ? Implement this method.

% initialize with a complex distance function

D0 = compute_levelset_shape('small-disks', n);

% set parameters

lambda = 0.8;

c1 = ...; % black

c2 = ...; % gray

...

Segmentation with Chan-Vese active contour without edges.

Copyright © 2006 Gabriel Peyré