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.
path(path, 'toolbox_graph/');
path(path, 'toolbox_graph/off/');
path(path, 'toolbox_graph/toolbox/');
path(path, 'meshes/');
% for large meshes, you can try 'skull', 'fandisk', 'lionhead', 'bunny'
% for small meshes, you can try 'mushroom', 'nefertiti', 'hammerheadtriang'
name = 'elephant50kv';
% load from file
[vertex,face] = read_mesh(name);
% display the mesh
clf;
plot_mesh(vertex,face);
% add the display of the triangle
shading faceted;
% when you rotate the camera, focuss the light using
camlight;


Exemple of a rendered triangle meshe (interp and faceted). 
L0 = D + W  (symmetric, nonnormalized)  
L1 = D^{1}*L0 = Id + D^{1}*W  (nonsymmetric, normalized)  
L2 = D^{1/2}*L0*D^{1/2} = Id + D^{1/2}*W*D^{1/2}  (symmetric, normalized) 
W(i,j)=1  (combinatorial)  
W(i,j)=1/vivj^2  (distance)  
W(i,j)=cot(alpha_ij)+cot(beta_ij)  (harmonic) 
% 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);
G0 = compute_mesh_gradient(vertex,face,laplacian_type,options);
disp(['Error (should be 0): ' num2str(norm(L0G0'*G0, 'fro')) '.']);
options.normalize = 1;
L1 = compute_mesh_laplacian(vertex,face,laplacian_type,options);
G1 = compute_mesh_gradient(vertex,face,laplacian_type,options);
disp(['Error (should be 0): ' num2str(norm(L1G1'*G1, 'fro')) '.']);
% these matrices are stored as sparse matrix
spy(L0);
% compute the unsymmetric 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. 
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. 
% solve the equation
vertex1 = ...;
% display
...
% 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(:,end10); % 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. 
% 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. 
% load a mesh
name = 'nefertiti';
[vertex,face] = read_mesh([name '.off']);
A = triangulation2adjacency(face);
% 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. 
% 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
...
% 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, 