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_geodesic_interpolation.m
function [M,G] = perform_geodesic_interpolation(W,points,f,options)

% perform_geodesic_interpolation - interpolate function values
%
%   [M,G] = perform_geodesic_interpolation(W,points,f,options);
%
%   f(:,i) is the value of the function at the point points(:,i).
%
%   options.method can be 'powerlaw' or 'gaussian'.
%
%   Copyright (c) 2007 Gabriel Peyr

n = size(W,1);
s = size(f,1); % dimenstionality
npoints = size(points,2);
points = round(points);
points = clamp(points, 1,n);

options.null = 1;
if isfield(options,'method')
method = options.method;
else
method = 'powerlaw';
end
if isfield(options,'sigma')
sigma = options.sigma;
else
sigma = 5/n;
end
if isfield(options,'sigma')
alpha = options.alpha;
else
alpha = 3;
end

% compute geodesic distances
G = zeros(n,n,1,npoints);

for i=1:npoints
progressbar(i,npoints);
[G(:,:,1,i),Z,Q] = perform_fast_marching(W, points(:,i), options);
end

% perform the interpolation
switch lower(method)
case 'powerlaw'
S = 1 ./ ( G.^alpha + sigma^alpha );
case 'gaussian'
S = exp( -(G/sigma).^2 );
otherwise
error('Unknown method.');
end

D = repmat( sum(S,4), [1 1 1 npoints] );
D(D<eps) = 1;
S = S ./ D;
S = repmat(S, [1 1 s 1]);
T = repmat( reshape(f, [1 1 s npoints]), [n n 1 1] );
M = S .* T;
M = sum( M, 4 );