Code covered by the BSD License

# Toolbox Fast Marching

### Gabriel Peyre (view profile)

24 Oct 2004 (Updated )

A toolbox for the computation of the Fast Marching algorithm in 2D and 3D.

perform_active_contour(D0, motion, options)
```function D = perform_active_contour(D0, motion, options)

% perform_active_contour - perform active contour resolution
%
%   D = perform_active_contour(D0, motion, options);
%
%   Level set implementation of various active contour, all related to
%   motion by mean curvature.
%
%   The resolution is done by an explicit euler, the time sep is options.dt
%   and should be quite small for stability. The maximum time is
%   options.Tmax.
%
%   In the following, one denotes the curvature of the level sets of D as
%   and D' = d(D)/dt the time derivative.
%
%   motion can be:
%       'mean': D' = |grad(D)| * Curv(D)
%       'affine': D' = |grad(D)| * Curv(D)^(1/3)
%       'errosion': D' = |grad(D)| * max(Curv(D),0)^(1/3)
%       'chan-vese': D' = |grad(D)| * ( Curv(D) - lambda*(D-c1)^2 + lambda*(D-c2)^2 )
%
%   You can provide the parameters in options.E, options.c1, options.c2,
%   options.lambda. You can turn on the automatic update of c1/c2 using
%   options.update_c.
%
%   The distance function D is redistanced every options.redistance_freq
%   iterations.
%
%   To turn off the display, options.do_display=0. You can provide a
%   background image in options.M.
%
%
%   Copyright (c) 2007 Gabriel eyre

dt              = getoptions(options, 'dt', 1000);
Tmax            = getoptions(options, 'Tmax', 1000);
do_display      = getoptions(options, 'display', 1);
display_freq    = getoptions(options, 'display_freq', 20);
redistance_freq = getoptions(options, 'redistance_freq', 30);
lambda          = getoptions(options, 'lambda', 0.8);
update_c        = getoptions(options, 'update_c', 1);
c1              = getoptions(options, 'c1', 0.1);
c2              = getoptions(options, 'c2', 0.8);
nb_svg          = getoptions(options, 'nb_svg', 0);
solver          = getoptions(options, 'solver', 'cg');
svg_path        = getoptions(options, 'svg_path', ''); % ['results/active-contour/' motion '/']);
if strcmp(motion,'chan-vese') || strcmp(motion, 'snake') || strcmp(motion, 'chan-vese-user')
E           = getoptions(options, 'E', 0, 1);
else
E = [];
end
M               = getoptions(options, 'M', E(:,:,1));
options.bound = 'per';

if not(isempty(svg_path)) && not(exist(svg_path))
mkdir(svg_path);
end

nb_phase = 1;
if not(isempty(E))
nb_phase = 2 - (size(E,3)==1);
end

if strcmp(motion, 'snake')
% pre compute gradient of the energy
end

D = D0;
n = size(D0,1);
epsilon = 1e-3;
niter = round(Tmax/dt);

if nb_svg>0
% distribute more svg at the begining of the iterations
svg_list = round( 1 + linspace(0,1,nb_svg+2).^3*niter ); svg_list(1) = [];
% delta_svg = niter/nb_svg;
if not(exist(svg_path))
mkdir(svg_path);
end
else
delta_svg = Inf;
end
next_svg = 0;
nb = 0;

for i=1:niter
progressbar(i,niter);
for k=1:nb_phase
d(:,:,k) = max(epsilon, sqrt(sum(g0(:,:,:,k).^2,3)) );
g(:,:,:,k) = g0(:,:,:,k) ./ repmat( d(:,:,k), [1 1 2] );
end

switch lower(motion)
case 'mean'
dD = d .* divgrad( g,options );
case 'affine'
dD = d .* sign(gg) .* abs(gg).^(1/3);
case 'errosion'
case 'snake'
dD = E .* d .* divgrad( g,options ) + sum(dE.*g0, 3);
case 'chan-vese'
% compute the inner / outer constant
if update_c
c1 = mean( E(D>=0) );
c2 = mean( E(D<0) );
end
dD = d .* ( divgrad( g,options ) - lambda*(E-c1).^2 + lambda*(E-c2).^2 );
case 'chan-vese-user'
if nb_phase==1
% single phase
dD = d .* ( divgrad( g,options ) + E );
else
for k=1:nb_phase
dD(:,:,k) = d(:,:,k) .* ( dG - E(:,:,2*k-1).*(D(:,:,3-k)>0) - E(:,:,2*k).*(D(:,:,3-k)<=0)  );
end
end
otherwise
error('Unknown motion');
end
D = D + dt*dD;
else
options.E = [];
switch lower(motion)
case 'mean'
C = zeros(n);
case 'snake'
options.E = E;
C = sum(dE.*g0, 3);
case 'chan-vese'
% compute the inner / outer constant
if update_c
c1 = mean( E(D>=0) );
c2 = mean( E(D<0) );
end
C = - lambda*(E-c1).^2 + lambda*(E-c2).^2;
case 'chan-vese-user'
C = E;
otherwise
error('Unknown motion');
end
% right hand side
y = D + dt*C;
options.d = d;
options.dt = dt;
options.n = n; options.ncols = n^2;
options.niter_max = 20;
ptions.epsilon = 1e-8;
options.x = D(:);
D = reshape(D,n,n);
end

if mod(i,redistance_freq)==0
for k=1:nb_phase
D(:,:,k) = perform_redistancing(D(:,:,k), options);
end
end
if do_display && ( mod(i,display_freq)==1 || i>=next_svg )
A = prod(D,3);
if strcmp(motion, 'snake') || strcmp(motion, 'chan-vese') || strcmp(motion, 'chan-vese-user')
A = M;
end
clf;
display_segmentation(D,A);
drawnow;

if 0
clf;
hold on;
if strcmp(motion, 'snake') || strcmp(motion, 'chan-vese') || strcmp(motion, 'chan-vese-user')
imagesc(M); axis image; axis off; axis([1 n 1 n]);
else
imagesc(D); axis image; axis off; axis([1 n 1 n]);
end
%        h = plot(c(1,:), c(2,:), 'r');
[c,h] = contour(D,[0 0], 'r');
set(h, 'LineWidth', 2);
drawnow;
hold off;
colormap gray(256);
end

if i>=next_svg
nb = nb+1;
% save image
saveas(gcf, [svg_path motion '-' num2string_fixeddigit(nb, 2) '.png'], 'png');
next_svg = svg_list(nb); % next_svg+delta_svg;
end
end
end```