Geodesic Bending Invariants for Surfaces
This tour explores the computation of bending invariants of surfaces.
Contents
Installing toolboxes and setting up the path.
You need to download the following files: signal toolbox, general toolbox and graph toolbox.
You need to unzip these toolboxes in your working directory, so that you have toolbox_signal, toolbox_general and toolbox_graph in your directory.
For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.
Recommandation: You should create a text file named for instance numericaltour.sce (in Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab command you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.
Execute this line only if you are using Matlab.
getd = @(p)path(p,path); % scilab users must *not* execute this
Then you can add the toolboxes to the path.
getd('toolbox_signal/'); getd('toolbox_general/'); getd('toolbox_graph/');
Bending Invariants
Bending invariant replace the position of the vertices in a mesh by new position that are insensitive to isometric deformation of the mesh. This defines a bending invariant signature that can be used for surface matching.
Bending invariant were introduced in:
On bending invariant signatures for surfaces, A. Elad and R. Kimmel, IEEE Transactions onPattern Analysis and Machine Intelligence, Vol. 25(10), p. 1285-1295, 2003.
A related method was developped for brain flattening in:
E.L. Schwartz and A. Shaw and E. Wolfson, _A Numerical Solution to the Generalized Mapmaker's Problem: Flattening Nonconvex Polyhedral Surfaces, IEEE Transactions on Pattern Analysis and Machine Intelligence, 11(9), p. 1005-1008, 1989.
This method is related to the Isomap algorithm for manifold learning:
A Global Geometric Framework for Nonlinear Dimensionality Reduction, J. B. Tenenbaum, V. de Silva and J. C. Langford, Science 290 (5500): 2319-2323, 22 December 2000
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.
Bending Invariant by Strain Minimization
A first method to compute bending invariant is by strain minimization, that requires the eigendecomposition of a kernel matrix. This corresponds to the Isomap algorithm method.
Load a mesh.
name = 'camel';
options.name = name;
[vertex,faces] = read_mesh(name);
n = size(vertex,2);
Display it.
clf; plot_mesh(vertex,faces, options);
Exercice 1: (the solution is exo1.m) Compute the geodesic distance matrix D between all samplg points on the mesh. It is going to take some of time.
exo1;
Compute a centered kernel for the Landmarks, that should be approximately a matrix of inner products.
J = eye(n) - ones(n)/n; K = -1/2 * J*(D.^2)*J;
Perform classical MDS on the reduced set of points, to obtain new positions in 3D.
opt.disp = 0;
[Xstrain, val] = eigs(K, 3, 'LR', opt);
Xstrain = Xstrain .* repmat(sqrt(diag(val))', [n 1]);
Xstrain = Xstrain';
Display the bending invariant surface.
clf; plot_mesh(Xstrain,faces, options);
Bending Invariant with Stress Minimization using SMACOF
The spectral localization X computed before is minimizing a weird functional, the strein. In practice, it is much better to minimize the stress, that is a non-convex functional.
Initialize the positions for the algorithm.
Xstress = vertex;
Compute the distance matrix.
D1 = repmat(sum(Xstress.^2,1),n,1); D1 = sqrt(D1 + D1' - 2*Xstress'*Xstress);
Compute the value of the stress function given by the configuration Xstress.
s = sqrt( sum( abs(D(:)-D1(:)).^2 ) / n^2 );
disp(['Stress = ' num2str(s) ]);
Stress = 0.24905
Compute the scaling matrix for the SMACOF update.
B = -D./max(D1,1e-10); B = B - diag(sum(B));
Update the positions.
Xstress = (B*Xstress')' / n;
Exercice 2: (the solution is exo2.m) Perform the SMACOF iterative algorithm. Plot the decay of the stress.
exo2;
Plot the optimized invariant shape.
clf; plot_mesh(Xstress,faces, options);
Plot stress evolution during minimization.
clf; plot(stress(1:end), '.-'); axis('tight');
Value of the stress for the SMACOF minimizer.
D1 = repmat(sum(Xstress.^2,1),n,1);
D1 = sqrt(D1 + D1' - 2*Xstress'*Xstress);
s = sqrt( sum( abs(D(:)-D1(:)).^2 ) / n^2 );
disp(['Stress (for SMACOF minimizer) = ' num2str(s) ]);
Stress (for SMACOF minimizer) = 0.047388
Surface Retrieval with Bending Invariant.
One can compute a bending invariant signature for each mesh in a library of 3D surface.
Isometry-invariant retrival is then perform by matching the invariant signature.
Exercice 3: (the solution is exo3.m) Implement a surface retrival algorithm based on these bending invariants.
exo3;