# Lecture 4 - Mesh Processing

Abstract : The goal of this lecture is to manipulate a 3D mesh. This includes the loading and display of a 3D mesh and then the processing of the mesh. This processing is based on computations involving various kinds of Laplacians. These Laplacians are extensions of the classical second order derivatives to 3D meshes. They can be used to perform heat diffusion (smoothing), compression and parameterization of the mesh.

## Setting up Matlab.

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

• You can load a mesh from a file and display it. Remember that a mesh is given by two matrix: vertex store the 3D vertex positions and face store 3-tuples giving the number of the vertices forming faces.
• % for large meshes, you can try 'skull', 'fandisk', 'lion-head', 'bunny'
% for small meshes, you can try 'mushroom', 'nefertiti', 'hammerheadtriang'
name = 'elephant-50kv';
% display the mesh
clf;
plot_mesh(vertex,face);
% add the display of the triangle
% when you rotate the camera, focuss the light using
camlight;

Exemple of a rendered triangle meshe (interp and faceted).

## Heat Diffusion on 3D Meshes.

• A Laplacian on a 3D mesh is a (n,n) matrix L, where n is the number of vertices, that generalize the classical laplacian of image processing to a 3D mesh. There are several forms of laplacian depending whether it is symmetric (L'=L) and normalized (1 on the diagonal) :
 L0 = D + W (symmetric, non-normalized) L1 = D^{-1}*L0 = Id + D^{-1}*W (non-symmetric, normalized) L2 = D^{-1/2}*L0*D^{-1/2} = Id + D^{-1/2}*W*D^{-1/2} (symmetric, normalized)
Where W is a weight matrix, W(i,j)=0 if (i,j) is not an edge of the graph. There are several kinds of such weights
 W(i,j)=1 (combinatorial) W(i,j)=1/|vi-vj|^2 (distance) W(i,j)=cot(alpha_ij)+cot(beta_ij) (harmonic)
where {vi}_i are the vertex of the mesh, i.e. vi=vertex(:,i), and (alpha_ij,beta_ij) is the pair of adjacent angles to the edge (vi,vj). A gradient matrix G associated to the laplacian L is an (m,n) matrix where m is the number of edges in the mesh, that satisfies L=G'*G. It can be computed as G((i,j),k)=+sqrt(W(i,j)) if k=j and G((i,j),k)=-sqrt(W(i,j)) if k=i, and G((i,j),k)=0 otherwise. In the following, you can compute gradient, weights and laplacian using the compute_mesh_xxx functions.
• % kind of laplacian, can be 'combinatorial', 'distance' or 'conformal' (slow)
laplacian_type = ...;
% load two different kind of laplacian and check the gradient factorization
options.symmetrize = 1;
options.normalize = 0;
L0 = compute_mesh_laplacian(vertex,face,laplacian_type,options);
disp(['Error (should be 0): ' num2str(norm(L0-G0'*G0, 'fro')) '.']);
options.normalize = 1;
L1 = compute_mesh_laplacian(vertex,face,laplacian_type,options);
disp(['Error (should be 0): ' num2str(norm(L1-G1'*G1, 'fro')) '.']);
% these matrices are stored as sparse matrix
spy(L0);

• In order to smooth a vector f of size n (i.e. a function defined on each vertex of the mesh), one can perform a heat diffusion by solving the following PDE
d F / dt = -L*F       with       F(x,t=0)=f
until some stopping time t.
When this diffusion is applied to each component of the positions of the vertices f=vertex(i,:), this smoothes the 3D mesh. Implement this PDE using an explicit discretization in time.
• % compute the un-symmetric laplacian
options.symmetrize = 0;
options.normalize = 1;
L = compute_mesh_laplacian(vertex,face,laplacian_type,options);
% the time step should be small enough
dt = 0.1;
% stopping time
Tmax = 6;
% number of steps
niter = round(Tmax/dt);
% initialize the 3 vectors at time t=0
vertex1 = vertex;
% solve the diffusion
for i=1:niter
% update the position by solving the PDE
vertex1 = vertex1 + ...;
end

 Heat diffusion on a 3D mesh, at times t=0, t=10, t=40, t=200.

• The heat equation can be solved to perform denoising of a damaged 3D mesh. In order to simulate this, one can add some noise to a known 3D mesh. This noise is added in the normal direction to each component of the vertices. In order to select the best denoising, one can test different stopping time Tmax and select the one that gives the smallest denoising error.
• nvert = size(vertex,2);
% compute the normals
normals = compute_normal(vertex,face)';
% add some noise to the vertices
sigma = .01; % noise level
rho = randn(1,nvert)*sigma;
vertex1 = vertex + repmat(rho,3,1).*normals;
% solve the heat equation
for i=1:niter
vertex1 = ...;
% record the denoising error
error(i) = ...;
end
% display
...
% compute the optimal diffusion time t as a function of the noise sigma

[errmin,topt] = min(error);
% compute the evolution of topt as a function of the noise level sigma
...

 Mesh denoising with several stopping times Tmax.

• Another way to smooth a mesh is to perform the following quadratic regularization for each componnent f=vertex(i,:)
F = argmin_g |f-g|^2 + t * |G*g|^2
The solution of this optimization is given in closed form using the Laplacian L=G'*G as the solution of the following linear system:
(Id+t*L)*F = f
Solve this problem for various t on a 3D mesh. You can use the operator \ to solve a system. How does this method compares with the heat diffusion ?
• % solve the equation
vertex1 = ...;
% display
...

## Combinatorial Laplacian: Spectral Decomposition and Compression.

• The combinatorial laplacian is a linear operator (thus a NxN matrix where N is the number of vertices). It depends only on the connectivity of the mesh, thus on face only. The eigenvector of this matrix forms an orthogonal basis of the vector space of signal of NxN values (one real value per vertex). Those functions are the extension of the Fourier oscillating functions to surfaces. For a small mesh (less than 1000) vertices, one can compute this set of vectors using the svd functions. For large meshes, one can compute only a small (e.g. 50) number of low pass eigenvectors using the sparse SVD function svds.
• % combinatorial laplacian computation
options.symmetrize = 1;
options.normalize = 0;
L0 = compute_mesh_laplacian(vertex,face,'combinatorial',options);
%% Performing eigendecomposition
if size(vertex,2)<1000
[U,D] = eig(full(L0));
else
opts.disp = 0;
[U,D] = eigs(L0,50,'SM',opts);
U = real(U(:,end:-1:1));
end

% extract one of the eigenvectors
c = U(:,end-10); % you can try with other vector with higher frequencies
% assign a color to each vertex
options.face_vertex_color = rescale(c, 0,255);
% display
clf;
plot_mesh(vertex,face, options);

 Eigenvectors of the combinatorial laplacian with increasing frequencies from left to right.

• Like the Fourier basis, the laplacian eigenvector basis can be used to perform an orthogonal decomposition of a function. In order to perform mesh compression, we decompose each coordinate X/Y/Z of the mesh into this basis. Once this decomposition has been performed, a compression is achieved by keeping only the low pass coefficients. To perform this compression, you need to load a small mesh, like 'venus' or 'mushroom'.
• % load a small mesh (like 'venus')
...
% warning : vertex should be of size 3 x nvert
keep = 5; % number of percents of coefficients kept
pf = (U'*vertex'); % projection of the vector in the laplacian basis
% display the spectrum pf
...
% set to zero high pass coefficients
q = round(keep/100*nvert); % number of zeros
pf(q+1:end,:) = 0;
vertex1 = (U*pf)';
% display the reconstructed mesh
...

 Spectral mesh compression performed by decomposition on the eigenvectors of the laplacian.

## Mesh Flattening and Parameterization.

• In order to do Mesh parameterization, one looks for 3 harmonic vectors vertex1(i,:) that satisfy L*vertex1(i,:)=0 inside the mesh, and with the position of the boundary vertices fixed to some closed convex curve. The flattening can then be proved to be 1:1 (no triangle flip) and long as the curve is convex and weights of the laplacian are positive. Mesh parameterization involve the solution of a sparse linear system, which is quite fast.
name = 'nefertiti';
% compute the parameterization
xy = compute_parameterization(vertex,face,options);
options.boundary = 'circle'; % the other choices are 'square' and 'triangle'
options.laplacian = 'combinatorial'; % the other choice is 'conformal'
% display
clf;
plot_graph(A,xy,'k.-');
axis tight; axis square; axis off;
% test with other laplacian and boundary conditions
...

 The first row shows parameterization using the combinatorial laplacian (with various boundary conditions). This assumes that the edge lengths is 1. To take into account the geometry of the mesh, the second row uses the conformal laplacian. The last row shows another example of parameterization with a conformal laplacian.

• The eigenvector of the combinatorial laplacian can also be used to perform mesh flattening. Flattening means finding a 2D location for each vertex of the mesh. These two coordinates are composed by the eigenvectors n°2 and n°3 of the laplacian (the first eigenvector being constant). In order to have a better flattening, you can try the same thing but with the conformal laplacian, who approximate the Laplace Beltrami operator of the continuous underlying surface.
• % compute the laplacian (must be symmetric)
L = ...;
% compute eigendecomposition
[U,S] = eig(full(L));
% exctact the correct component from the eigenvectors
vertex1 = ...;
% display the flattened mesh
plot_graph(A,vertex1);
axis tight; axis square; axis off;
% do the flattening with the conformal laplcian
...

• Another way to compute the flattening is to use the Isomap algorithm. This algorithm is not based on local differential operator such as laplacian. Instead, the geodesic distance between points on the mesh graph is first computed (see the course 4 for example of geodesic computations on graphs). Then the 2D layout of point is computed in order to match the geodesic distance with the distance in the plane.
• % the embedding is now computed with isomap
options.method = 'isomap';
xy = compute_parameterization(vertex,face, options);
% display
plot_graph(A,xy,'k.-');
axis tight; axis square; axis off;

 Comparison of the flattening obtained with the combinatorial laplacian, the conformal laplacian and isomap. Top row shows a model for with the flattening is valid (1:1) and bottom rows shows a model with high iso-perimetric distortion with leads to a non-valid flattening.