Abstract : The goals of this lecture is to manipulate the fast marching algorithm on triangulated meshes. Application to geodesic extraction, Voronoi segmentation, remeshing and bending invariants are presented.
path(path, 'toolbox_fast_marching/');
path(path, 'toolbox_fast_marching/toolbox/');
path(path, 'toolbox_graph/');
path(path, 'toolbox_graph/off/');
path(path, 'meshes/');
cd toolbox_fast_marching
compile_mex;
cd ..
% read the 3D mesh
name = 'elephant50kv'; % other choices includes 'skull' or 'bunny'
[vertex,faces] = read_mesh([name '.off']);
% select random starting points for the propagation
nverts = max(size(vertex)); % number of vertices
nstart = 8; % number of starting points
start_points = floor(rand(nstart,1)*nverts)+1;
% perform the front propagation, Q contains an approximate segementation
[D,S,Q] = perform_fast_marching_mesh(vertex, faces, start_points);
% display the result using the distance function
options.start_points = start_points;
plot_fast_marching_mesh(vertex,faces, D, [], options);
% extract precisely the voronoi regions, and display it
[Qexact,DQ, voronoi_edges] = compute_voronoi_mesh(vertex,faces, start_points, options);
options.voronoi_edges = voronoi_edges;
plot_fast_marching_mesh(vertex,faces, D, [], options);

Examples of voronoi segmentation on two different 3D meshes for an increasing number of starting points. 
options.nb_iter_max = ...; % trying with a varying number of iterations
[D,S,Q] = perform_fast_marching_mesh(vertex, faces, start_points, options);
% display the distance
...

Example of front propagation. 
% select an initial starting point
x = ...;
% compute the vertices 1ring
vring = compute_vertex_ring(faces);
% initiate the path
path = x;
% build the path
while true
% find next vertex
path{end+1} = ...;
...
end
% display the path
...
% in order to extract a smooth path, one needs to use a gradient descent
x = = floor(rand(20,1)*nverts)+1; % select random ending points for the geodesics
options.method = 'continuous';
paths = compute_geodesic_mesh(D, vertex, faces, x, options);
plot_fast_marching_mesh(vertex,faces, Q, paths, options);

Geodesics between several ending points and an increasing
number of starting points (from left to right: 1, 10, 20, 50 and 100 points). The color corresponds to the geodesic distance to the closest starting point. 
nbr_landmarks = ...; % number of points, eg. 400
W = ones(nverts,1); % the speed function, for now constant speed to perform uniform remeshing
% perform the sampling of the surface
landmark = farthest_point_sampling_mesh( vertex,faces, [], nbr_landmarks, options );
% compute the associated triangulation
[D,Z,Q] = perform_fast_marching_mesh(vertex, faces, landmark);
[vertex_voronoi,faces_voronoi] = compute_voronoi_triangulation_mesh(Q,vertex,faces);
% display the distance function (same as before)
...
% display the remeshed triangulation (same but with vertex_voronoi and faces_voronoi)
...
% you need in fact to reorient the mesh
faces_voronoi = perform_faces_reorientation(vertex_voronoi,faces_voronoi);

Farthest point sampling with an increasing number of
points (upper row) and corresponding remeshing (bottom row). 
% first kind of speed: low on the left, high on the right
v = rescale(vertex(1,:));
options.W = rescale(v>0.5, 1, 3); options.W = options.W(:);
% do the remeshing
...
% second kind of speed: continuously increasing
v = rescale(vertex(1,:),1,8);
options.W = v(:);
% do the remeshing
...

Rows 1&2: uniform sampling and remeshing. Rows 3: adapted (split left/right) remeshing. Rows 4: adapted (continuously increasing) remeshing. 
nlandmarks = 300; % number of landmarks (low to speed up)
landmarks = floor(rand(nlandmarks,1)*nverts)+1; % samples landmarks at random
% compute the distance between landmarks and the rest of the vertices
Dland = zeros(nverts,nlandmarks);
for i=1:nlandmarks
fprintf('.');
[d,S,Q] = perform_fast_marching_mesh(vertex, faces, landmarks(i));
Dland(:,i) = d(:).^2;
end
fprintf('\n');
% perform isomap on the reduced set of points
D1 = Dland(landmarks,:); % reduced pairwise distances
D1 = (D1+D1')/2; % force symmetry
J = eye(nlandmarks)  ones(nlandmarks)/nlandmarks; % centering matrix
K = 1/2 * J*D1*J; % inner product matrix
% compute the rank3 approximation of the inner product to compute embedding
opt.disp = 0;
[xy, val] = eigs(K, 3, 'LR', opt);
xy = xy .* repmat(1./sqrt(diag(val))', [nlandmarks 1]);% interpolation on the full set of points
% extend the embeding using geodesic interpolation
vertex1 = zeros(nverts,3);
deltan = mean(Dland,1);
for x=1:nverts
deltax = Dland(x,:);
vertex1(x,:) = 1/2 * ( xy' * ( deltandeltax )' )';
end
% display both the original mesh and the embedding
...
Original mesh (top row) and bending invariant (bottom row).