Toolbox Fast Marching - A toolbox for Fast Marching and level sets computations
Copyright (c) 2008 Gabriel Peyre
Contents
The toolbox can be downloaded from Matlab Central http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=6110&objectType=FILE
Introduction to Fast Marching and Level Sets
The Fast Marching algorithm, introduced by J. Sethian (1996) is a numerical algorithm that is able to catch the viscosity solution of the Eikonal equation norm(grad(D))=P. The level set {x \ F(x)=t} can be seen as a front advancing with speed P(x).
The resulting function D is a distance function, and if the speed P is constant, it can be seen as the distance function to a set of starting points.
The Fast Marching is very similar to the Dijkstra algorithm that finds shortest paths on graphs. Using a gradient descent of the distance function D, one is able to extract a good approximation of the shortest path (geodesic) in various settings (Euclidean for P constant, and a weighted Riemanian manifold with P varying).
The main reference about the Fast Marching algorithm is the book
- Level Set Methods and Fast Marching Methods Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science, J.A. Sethian, Cambridge University Press, 1999, Cambridge Monograph on Applied and Computational Mathematics
A good review of the Fast Marching in 3D together with some applications can be found in
- Fast extraction of minimal paths in 3D images and application to virtual endoscopy, T.Deschamps and L.D. Cohen, Medical Image Analysis, volume 5, Issue 4, December 2001.
Setting-up the Path and Compiling the Mex Files.
Befor begining, we add some directories to the path.
path(path, 'toolbox/'); path(path, 'data/');
The main computation are done in a mex file so it is quite fast (using a Fibonacci heap structure). To compile the mex files, simply run the following lines.
% set up a C compiler % mex -setup % uncomment this line % compile % compile_mex; % uncomment this line
2D Fast Marching Computations
The toolbox allows to compute a distance map from a set of starting point using an arbitrary isotropic metric. This distance map is then used to compute geodesic paths from any point that joins the closest starting point, according to the geodesic distance.
We first load a 2D map, and ask to the user for the starting point.
n = 400; [M,W] = load_potential_map('road2', n); start_point = [16;219]; end_point = [394;192]; % You can use instead the function % [start_point,end_point] = pick_start_end_point(M);
We then perform the front propagation
clear options; options.nb_iter_max = Inf; options.end_points = end_point; % stop propagation when end point reached [D,S] = perform_fast_marching(W, start_point, options); % nicde color for display A = convert_distance_color(D); imageplot({W A}, {'Potential map' 'Distance to starting point'}); colormap gray(256);

At last we can extract the geodesic curve, and display the result.
gpath = compute_geodesic(D,end_point);
plot_fast_marching_2d(W,S,gpath,start_point,end_point, options);
title('Shortest path');

One can also use several starting and ending points and perform progressive propagation
[M,W] = load_potential_map('mountain', n); % random seeding of the points nstart = 8; nend = 30; start_points = round( rand(2,nstart)*(n-1)+1 ); end_points = round( rand(2,nend)*(n-1)+1 ); % FM computation options.end_points = []; A = {}; nlist = round( linspace(0.05,1,8)*n^2 ); for i=1:length(nlist) options.nb_iter_max = nlist(i); [D,S] = perform_fast_marching(W, start_points, options); mask = repmat(S==-1,[1 1 3]); A{i} = convert_distance_color(D); A{i} = mask.*A{i} + (1-mask) .* repmat(rescale(W),[1 1 3]); end clf; imageplot(A, '', 2,4); % paths extractions gpath = compute_geodesic(D,end_points); % display clf; subplot(1,2,1); imageplot(W,'Potential map'); subplot(1,2,2); plot_fast_marching_2d(A{end},[],gpath,start_points,end_points, options); title('Shortest paths');

3D Volumetric Datasets
You can load a 3D volume of data.
% load the whole volume load ../toolbox_fast_marching_data/brain1-crop-256.mat % crop to retain only the central part n = 100; M = rescale( crop(M,n) ); % display some horizontal slices slices = round(linspace(10,n-10,6)); Mlist = mat2cell( M(:,:,slices), n, n, ones(6,1)); clf; imageplot(Mlist);

You can also perform a volumetric display. By changing the value of options.center, you can display various level sets.
clf;
options.center = .35;
imageplot(M,options);
title('center=.35');

3D Fast Marching
The procedure for volumetric FM is the same as for 2D FM.
We ask to the user for some input points.
delta = 5;
start_point = [91;15;delta];
% You can use instead the function pick_start_end_point
We compute a potential that is high only very close to the value of M at the selected point.
W = abs(M-M(start_point(1),start_point(2),start_point(3))); W = rescale(-W,.001,1); % perform the front propagation options.nb_iter_max = Inf; [D,S] = perform_fast_marching(W, start_point, options); % display the distance map D1 = rescale(D); D1(D>.7) = 0; clf; imageplot(D1,options); alphamap('rampup'); colormap jet(256);

Fast Marching on 3D Meshes
The FM algorithm also works on triangulated meshes (see Kimmel and Sethian paper). The metric is the metric induced by the embedding of the 3D surface in R^3, but one can also modulate this metric by a scalar metric map.
The triangulated mesh is loaded from a .off file. You can see my graph toolbox for more functions for mesh processing.
% load the mesh [vertex,faces] = read_mesh('elephant-50kv'); % display the mesh clf; plot_mesh(vertex, faces); shading interp;

Select some starting point and do the propagation.
start_points = 20361; [D,S,Q] = perform_fast_marching_mesh(vertex, faces, start_points, options);
Extract some geodesics and display the result.
npaths = 30; nverts = size(vertex,2); % select some points that are far enough from the starting point [tmp,I] = sort( D(:) ); I = I(end:-1:1); I = I(1:round(nverts*1)); end_points = floor( rand(npaths,1)*(length(I)-1) )+1; end_points = I(end_points); % precompute some usefull information about the mesh options.v2v = compute_vertex_ring(faces); options.e2f = compute_edge_face_ring(faces); % extract the geodesics options.method = 'continuous'; options.verb = 0; paths = compute_geodesic_mesh(D,vertex,faces, end_points, options); % display options.colorfx = 'equalize'; plot_fast_marching_mesh(vertex,faces, D, paths, options); shading interp;

Fast Marching Inside a 2D Shape
It is possible to retrict the propagation inside a 2D shape.
First we load a 2D shape, which is a BW image
n = 128;
M = rescale( load_image('chicken',n), 0,1 );
M = double(M>0.5);
Compute geodesic distance to a point inside the shape
start_points = [65;65]; % use constant metric W = ones(n); % create a mask to restrict propagation L = zeros(n)-Inf; L(M==1) = +Inf; options.constraint_map = L; % do the FM computation [D,S,Q] = perform_fast_marching(W, start_points, options);
Select a set of end points, here we locate them on the boundary of the shape.
bound = compute_shape_boundary(M); nbound = size(bound,1); npaths = 40; sel = round(linspace(1,nbound+1,npaths+1)); sel(end) = []; end_points = bound(sel,:)';
Extract the geodesics from these ending points
D1 = D; D1(M==0) = 1e9; paths = compute_geodesic(D1,end_points);
We display the result.
ms = 30; lw = 3; % size for plot % display A = convert_distance_color(D); clf; hold on; imageplot(A); axis image; axis off; for i=1:npaths h = plot( paths{i}(2,:), paths{i}(1,:), 'k' ); set(h, 'LineWidth', lw); h = plot(end_points(2,i),end_points(1,i), '.b'); set(h, 'MarkerSize', ms); end h = plot(start_points(2),start_points(1), '.r'); set(h, 'MarkerSize', ms); hold off; colormap jet(256); axis ij;

Farthest Point Sampling
One can sample an image using a greedy algorithm that distribute, at each step, a new point that is located the farthest away from the previously sampled set of points. See the paper of Peyre and Cohen, IJCV.
clear options; n = 400; [M,W] = load_potential_map('mountain', n); npoints_list = round(linspace(20,200,6)); landmark = []; options.verb = 0; ms = 15; clf; for i=1:length(npoints_list) nbr_landmarks = npoints_list(i); landmark = perform_farthest_point_sampling( W, landmark, nbr_landmarks-size(landmark,2), options ); % compute the associated triangulation [D,Z,Q] = perform_fast_marching(W, landmark); % display sampling and distance function D = perform_histogram_equalization(D,linspace(0,1,n^2)); subplot(2,3,i); hold on; imageplot(D'); plot(landmark(1,:), landmark(2,:), 'r.', 'MarkerSize', ms); title([num2str(nbr_landmarks) ' points']); hold off; colormap jet(256); end

Anisotropic Fast Marching
The toolbox also implements anisotropic FM, where at each location, the metric is given by a tensor field.
We first create a 2D vector field. The tensor field is aligned with this tensor field.
n = 200;
U = randn(n,n,2);
options.bound = 'per';
U = perform_vf_normalization( perform_blurring(U, 30,options) );
Create a set of ending points, located on the boundary of the image.
start_points = [n;n]/2;
% end points on the boundary
t = 1:n;
x = [t t*0+n t(end-1:-1:1) t*0+1];
y = [t*0+1 t t*0+n t(end-1:-1:1)];
npaths = 50;
s = round(linspace( 1,length(x), npaths+1) ); s(end) = [];
end_points = cat(1, x(s),y(s));
We build a metric with an increasing level of anisotropy, and perform the FM each time.
% test for various degree of anisotropy aniso_list = [.01 .05 .1 .2 .5 1]; ms = 30; lw = 3; % display params clf; for ianiso = 1:length(aniso_list) % build the tensor field aniso = aniso_list(ianiso); V = cat(3, -U(:,:,2), U(:,:,1)); % orthogonal vector T = perform_tensor_recomp(U,V, ones(n),ones(n)*aniso ); % propagation [D,S,Q] = perform_fast_marching(T, start_points); % for sexy display D1 = perform_histogram_equalization(D, linspace(0,1,n^2)); % extract tons of geodesics paths = compute_geodesic(D,end_points, options); % display subplot(2,3,ianiso); hold on; imageplot(D1); axis image; axis off; colormap jet(256); title(['Anisotropy=' num2str(1-aniso)]); for i=1:npaths end_point = end_points(:,i); h = plot( paths{i}(2,:), paths{i}(1,:), 'k' ); set(h, 'LineWidth', lw); h = plot(end_point(2),end_point(1), '.b'); set(h, 'MarkerSize', ms); end h = plot(start_points(2),start_points(1), '.r'); set(h, 'MarkerSize', ms); hold off; colormap jet(256); axis ij; end
