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.

test_geodesic_interpolation.m
% test for data interpolation using geodesic distances

rep  = 'data/';
rep  = '';

name = 'mm';
name = 'cavern';
name = 'disk';
name = 'toto';
name = 'boat';
name = 'cartoon';

clear options;
options.null = 0;

%% compute a binary shape
n = 300;
if strcmp(name, 'toto')
M = sum(M0,3);
M = 1-(M==M(1));
% compute the constraint map
CM = zeros(n) + Inf; CM(M==0) = -Inf;
options.constraint_map = CM; % constraint the propagation inside the shape
W = ones(n);
mode = 'interpolation';
else
n0 = [];
if strcmp(name, 'disk')
n0 = n;
end
if strcmp(name, 'cartoon')
n0 = 400;
end
M = rescale(crop(M,n)); M0 = M;
M = rescale(sum(M,3));
W = compute_edge_energy(M, 1.2);
W = 1./rescale(W,0.001,1);
mode = 'colorization';
end

% RGB2COL matrix
CM = rand(3); CM(:,1) = 1;
[CM,R] = qr(CM); CM(:,1) = 1/3; CM = CM';

%% pick some points
b = 1;
points = [];
while b==1
clf;
hold on;
imagesc(rescale(M)); axis image; axis off;
if not(isempty(points))
plot_scattered(points(end:-1:1,:));
end
colormap gray(256); axis ij;
hold off;
[y,x,b] = ginput(1);
if b==1
points(:,end+1) = [x;y];
end
end
npoints = size(points,2);
% random colors (sqrt to make it more colorfull
f = rand(3,npoints);
f = ( f ./ repmat(sum(f,1),[3, 1]) ).^(1);
% use fixed palette
fc = [1 0 0; 0 1 0; 0 0 1; 1 1 0; 1 0 1; 0 1 1; 0.2 0.7 1; 1 0.2 0.7; 0.7 1 0.2; 0.2 1 0.7]';
sel =  floor(rand(npoints,1)*6)+1;
f = fc(:,sel);

if strcmp(mode, 'colorization')
% random CbCr
f = CM*f; f = f(2:end,:);
end

%% perform interpolation
if strcmp(mode, 'colorization')
sigma_list = 0.001/n;
alpha_list = 3;
else
% the smaller, the more interpolation is performed
% the higher, the more diffusion it is
sigma_list = [25 25 25 0.01 0.01 0.01]/n;
% the higher, the more "voronoi"-like it is
alpha_list = [1 3 50 1 3 50];
end

nsigma = length(sigma_list);
nrows = 1;
if nsigma>3
nrows = 2;
end

clf;
for i=1:nsigma
options.sigma = sigma_list(i);
options.alpha = alpha_list(i);
[A,G] = perform_geodesic_interpolation(W,points,f,options);

if strcmp(mode, 'colorization')
% cat the luminosity component
A = cat(3, M,A);
A = reshape(A, [n^2 3])';
A = ( (CM^(-1)) * A )';
A = reshape( A, [n n 3] );
end
A = clamp(A);
subplot(nrows,ceil(nsigma/nrows), i);
hold on;
imagesc(rescale(A)); axis image; axis off; axis ij;
plot_scattered(points(end:-1:1,:));
hold off;
title(['\alpha=' num2str(options.alpha) ', \sigma=' num2str(options.sigma)]);
end