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.

Note : this lecture is a Scilab variant of my original Matlab lectures. If you want to have access to more functionalities for mesh processing (including parameterization and vizualisation), you should check the Matlab lectures.

Setting up Scilab.

Download the Scilab toolbox toolbox_graph_sci.zip. Unzip it into your working directory. You should have a directory toolbox_graph_sci/ in your path.
Start Scilab, and change of directory (using cd or using the File menu). To load all the function from the directory, execute 'getd()' from the command line.
You should write each part of the lecture in a separate file (for instance 'myfile.sce'). To execute your code, simply type 'exec(myfile.sce);' in the command line. If you need to create a function, then write it into a file (e.g. 'myfunc.sci'), but do not forget to reload the set of function using 'getd()'.

Mesh Loading and Displaying.

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. Beware that Scilab display is much slower than Matlab.

// for large meshes, you can try 'bunny', 'elephant-50kv'
// for small meshes, you can try 'mushroom', 'nefertiti'
name = 'mushroom';
// load from file
[vertex,face] = read_off(strcat([name '.off']));
// display the mesh
plot_mesh(vertex,face);
// add the light to the display
options.lighting = 1;
options.faceted = 1;
plot_mesh(vertex,face, options);
// set a color for display
options.face_vertex_color = vertex(2,:);
options.lighting = 0; // you can try also with 1
options.faceted = 0;
plot_mesh(vertex,face, options);
// try to rotate the camera, change the color of display, etc
...
// try to morph the mesh (using handcraft deformation of vertex, and using the normal)
normal = compute_normal(vertex,face);
...

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}L0D^{-1/2} = Id + D^{-1/2}WD^{-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);
// check that the laplacian is symmetric, test if L*1=0 or not
...

Using the weight matrix W that is non-symmetric but normalized, you can smooth each coordinate of a mesh

options.symmetrize=0;
options.normalize=1;
L = compute_mesh_weight(vertex,face,'combinatorial',options);
vertex1 = (W*vertex')'; // do you understand the role of ' ?
// display the mesh
...
// iterate the smoothing and display - relate this with the theory (convergence, etc)
...
// try with other kind of weights
...

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. Remember that F is in fact 3 functions (X/Y/Z of the mesh), but you can solve for all the three functions together !

// 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 = 5;
// 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 + (ones(3,1)*rho).*normals;
// solve the heat equation
for i=1:niter
        vertex1 = ...;
        // record the denoising error
        error(i) = ...;
end
// display the error decay
...
// 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
...
// display the denoised mesh - comments ?
...

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. Remember that f is in fact 3 functions (X/Y/Z of the mesh), but you can solve for all the three functions together ! 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.

// Warning: use this code on a SMALL mesh (e.g. mushroom or nefertiti)
// combinatorial laplacian computation
options.symmetrize = 1;
options.normalize = 0;
L0 = compute_mesh_laplacian(vertex,face,'combinatorial',options);
// Performing eigendecomposition
[U,lambda] = spec(full(L));
// display the spectrum lambda - comments ?
plot(lambda);
// 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 = c;
// display
clf;
plot_mesh(vertex,face, options);
// try for other eigenvectors - comments ?

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 - comments ?
...

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

Mesh Flattening.

The eigenvector of the 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.

// try first this method on 'nefertiti' mesh
// compute the laplacian (must be symmetric)
L = ...;
// compute eigendecomposition
pf = (U'*vertex');
// exctact the correct component from the eigenvectors
vertex1 = ...;
// display the flattened mesh
plot_graph(triangulation2adjacency(face),vertex1);
// do the flattening with the conformal laplacian - comment
...
// try on manequin
...
// try on a sphere with a hole of increasing size
...

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.

Bonus questions.

Instead of solving for the Heat equation (diffusion), one can solve for the wave equation on the mesh (propagation), which reads
d^2 F / dt^2 = -L*F with F(x,t=0)=f and F'(x,t=0)=g
Solve this problem for various t on a 3D mesh. Note that this time, F is just 1 function (not anymore the three coordinates X/Y/Z !).
Start with f being some randomly +1/-1 at some locations, and g=0. Comments ? How do you explain this ?
Modify f in order to be a set of little gaussian bump (positive and negative). Comments ?

// load a mesh, e.g. bunny, compute normalized laplacian
...
// the time step should be small enough
dt = 0.2;
// stopping time
Tmax = 20;
// number of steps
niter = round(Tmax/dt);
// initialize the vectors at time t=0
n = size(vertex,2);
f = ...; // set random +1/-1
// solve the equation
fprev = f;
for i=1:niter
        // update the position by solving the PDE
        fsvg = f;
        f = 2*f - ...;
        fprev = f;
end
// display
...
// use another f (gaussian bumps)
...

Wave propagation with several stopping times Tmax.
The initial condition f is a collection of little gaussian bumps.

Implement the parameterization using the Tutte method.

This requires :
1) identifying the vertex of the mesh that are on the boundary.
2) computing for these vertex a 2D position (x0,y0) on a circle (or any convex curve).
3) computing the laplacian.
4) for each index i on the boundary, set L(i,:)=0 and then L(i,i)=1.
5) compute a right hand size vector b=0 with b(i)=x0(i) for those i on the boundary. Solve for x=L\b.
6) same thing but for y.

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.