# Lecture 3 - Geodesic Computations on 3D Meshes

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.

## Setting up Matlab.

• The first thing to do is to install these toolboxes in your path.
• path(path, 'toolbox_fast_marching/');
path(path, 'toolbox_fast_marching/toolbox/');
path(path, 'toolbox_graph/');
path(path, 'toolbox_graph/off/');
path(path, 'meshes/');

• Recompile the mex files for your machine (this can produce some warning). If it does not work, either use the already compiled mex file (they should be available in toolbox_fast_marching/ for MacOs and Unix) or try to set up matlab with a C compiler (e.g. gcc) using 'mex -setup'.
• cd toolbox_fast_marching
compile_mex;
cd ..

## Distance Computation on Triangulated Surface.

• Using the fast marching on a triangulated surface, one can compute the distance from a set of input points. This function also returns the segmentation of the surface into geodesic Voronoi cells.
• % read the 3D mesh
name = 'elephant-50kv'; % other choices includes 'skull' or 'bunny'

% 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.

• You can even stop the propagation after a fixed number of steps and see the resulting partially propagated front.
• 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.

• In order to extract a geodesic from a given ending point x, you need to fill a vector path such that path(1)=x and D(path(i+1))<D(path(i)) with [path(i), path(i+1)] being an edge of the mesh. In order to do so, first compute the edge 1-ring (list of adjacent vertices to all the vertices in the mesh) and then iterate to construct path. If you do not manage to program this, you can use the function compute_discrete_geodesic_mesh that compute the geodesic.
• % select an initial starting point
x = ...;
% compute the vertices 1-ring
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.

## Geodesic Remeshing.

• A regular sampling of the surface can be performed by seeding in a greedy farthest fashion samples. These points can be linked according to the Voronoi cells adjacency which gives a powerful but yet simple remeshing scheme.
• 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 re-orient 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).

• One can use a non-constant speed function in order to have an adapted remeshing. The sampling will be denser in area where the speed function is low.
• % 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.

## Bending Invariants of a Surface.

• One can use the Isomap procedure in order to modify the surface to get bending invariant signature, useful to perform articulation-invariant recognition of 3D surfaces. One first need to compute the pairwise geodesic distances between points on the surfaces. One then look for new 3D positions of the vertices so that the new Euclidean distances matches closely the geodesic distances. In order to speed up computation, the geodesic distances are computed only on a reduced set of landmarks points. The 3D new locations are then interpolated.

• 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 rank-3 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' * ( deltan-deltax )' )';
end
% display both the original mesh and the embedding

...
Original mesh (top row) and bending invariant (bottom row).