\[ \newcommand{\NN}{\mathbb{N}} \newcommand{\CC}{\mathbb{C}} \newcommand{\GG}{\mathbb{G}} \newcommand{\LL}{\mathbb{L}} \newcommand{\PP}{\mathbb{P}} \newcommand{\QQ}{\mathbb{Q}} \newcommand{\RR}{\mathbb{R}} \newcommand{\VV}{\mathbb{V}} \newcommand{\ZZ}{\mathbb{Z}} \newcommand{\FF}{\mathbb{F}} \newcommand{\KK}{\mathbb{K}} \newcommand{\UU}{\mathbb{U}} \newcommand{\EE}{\mathbb{E}} \newcommand{\Aa}{\mathcal{A}} \newcommand{\Bb}{\mathcal{B}} \newcommand{\Cc}{\mathcal{C}} \newcommand{\Dd}{\mathcal{D}} \newcommand{\Ee}{\mathcal{E}} \newcommand{\Ff}{\mathcal{F}} \newcommand{\Gg}{\mathcal{G}} \newcommand{\Hh}{\mathcal{H}} \newcommand{\Ii}{\mathcal{I}} \newcommand{\Jj}{\mathcal{J}} \newcommand{\Kk}{\mathcal{K}} \newcommand{\Ll}{\mathcal{L}} \newcommand{\Mm}{\mathcal{M}} \newcommand{\Nn}{\mathcal{N}} \newcommand{\Oo}{\mathcal{O}} \newcommand{\Pp}{\mathcal{P}} \newcommand{\Qq}{\mathcal{Q}} \newcommand{\Rr}{\mathcal{R}} \newcommand{\Ss}{\mathcal{S}} \newcommand{\Tt}{\mathcal{T}} \newcommand{\Uu}{\mathcal{U}} \newcommand{\Vv}{\mathcal{V}} \newcommand{\Ww}{\mathcal{W}} \newcommand{\Xx}{\mathcal{X}} \newcommand{\Yy}{\mathcal{Y}} \newcommand{\Zz}{\mathcal{Z}} \newcommand{\al}{\alpha} \newcommand{\la}{\lambda} \newcommand{\ga}{\gamma} \newcommand{\Ga}{\Gamma} \newcommand{\La}{\Lambda} \newcommand{\Si}{\Sigma} \newcommand{\si}{\sigma} \newcommand{\be}{\beta} \newcommand{\de}{\delta} \newcommand{\De}{\Delta} \renewcommand{\phi}{\varphi} \renewcommand{\th}{\theta} \newcommand{\om}{\omega} \newcommand{\Om}{\Omega} \renewcommand{\epsilon}{\varepsilon} \newcommand{\Calpha}{\mathrm{C}^\al} \newcommand{\Cbeta}{\mathrm{C}^\be} \newcommand{\Cal}{\text{C}^\al} \newcommand{\Cdeux}{\text{C}^{2}} \newcommand{\Cun}{\text{C}^{1}} \newcommand{\Calt}[1]{\text{C}^{#1}} \newcommand{\lun}{\ell^1} \newcommand{\ldeux}{\ell^2} \newcommand{\linf}{\ell^\infty} \newcommand{\ldeuxj}{{\ldeux_j}} \newcommand{\Lun}{\text{\upshape L}^1} \newcommand{\Ldeux}{\text{\upshape L}^2} \newcommand{\Lp}{\text{\upshape L}^p} \newcommand{\Lq}{\text{\upshape L}^q} \newcommand{\Linf}{\text{\upshape L}^\infty} \newcommand{\lzero}{\ell^0} \newcommand{\lp}{\ell^p} \renewcommand{\d}{\ins{d}} \newcommand{\Grad}{\text{Grad}} \newcommand{\grad}{\text{grad}} \renewcommand{\div}{\text{div}} \newcommand{\diag}{\text{diag}} \newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} } \newcommand{\pdd}[2]{ \frac{ \partial^2 #1}{\partial #2^2} } \newcommand{\dotp}[2]{\langle #1,\,#2\rangle} \newcommand{\norm}[1]{|\!| #1 |\!|} \newcommand{\normi}[1]{\norm{#1}_{\infty}} \newcommand{\normu}[1]{\norm{#1}_{1}} \newcommand{\normz}[1]{\norm{#1}_{0}} \newcommand{\abs}[1]{\vert #1 \vert} \newcommand{\argmin}{\text{argmin}} \newcommand{\argmax}{\text{argmax}} \newcommand{\uargmin}[1]{\underset{#1}{\argmin}\;} \newcommand{\uargmax}[1]{\underset{#1}{\argmax}\;} \newcommand{\umin}[1]{\underset{#1}{\min}\;} \newcommand{\umax}[1]{\underset{#1}{\max}\;} \newcommand{\pa}[1]{\left( #1 \right)} \newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. } \newcommand{\enscond}[2]{ \left\{ #1 \;:\; #2 \right\} } \newcommand{\qandq}{ \quad \text{and} \quad } \newcommand{\qqandqq}{ \qquad \text{and} \qquad } \newcommand{\qifq}{ \quad \text{if} \quad } \newcommand{\qqifqq}{ \qquad \text{if} \qquad } \newcommand{\qwhereq}{ \quad \text{where} \quad } \newcommand{\qqwhereqq}{ \qquad \text{where} \qquad } \newcommand{\qwithq}{ \quad \text{with} \quad } \newcommand{\qqwithqq}{ \qquad \text{with} \qquad } \newcommand{\qforq}{ \quad \text{for} \quad } \newcommand{\qqforqq}{ \qquad \text{for} \qquad } \newcommand{\qqsinceqq}{ \qquad \text{since} \qquad } \newcommand{\qsinceq}{ \quad \text{since} \quad } \newcommand{\qarrq}{\quad\Longrightarrow\quad} \newcommand{\qqarrqq}{\quad\Longrightarrow\quad} \newcommand{\qiffq}{\quad\Longleftrightarrow\quad} \newcommand{\qqiffqq}{\qquad\Longleftrightarrow\qquad} \newcommand{\qsubjq}{ \quad \text{subject to} \quad } \newcommand{\qqsubjqq}{ \qquad \text{subject to} \qquad } \]

Geodesic Bending Invariants with Landmarks

Geodesic Bending Invariants with Landmarks

This tour explores the use of farthest point sampling to compute bending invariant with classical MDS (strain minimization).

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/');

Farthest Points Landmarks Seeding

For large mesh, computing all the pairwise distances is intractable. It is possible to speed up the computation by restricting the computation to a small subset of landmarks.

This seeding strategy was used for surface remeshing in:

Geodesic Remeshing Using Front Propagation, Gabriel Peyré and Laurent Cohen, International Journal on Computer Vision, Vol. 69(1), p.145-156, Aug. 2006.

Load a mesh.

name = 'elephant-50kv';
options.name = name;
[vertex,faces] = read_mesh(name);
nverts = size(vertex,2);

Display it.

clf;
plot_mesh(vertex,faces, options);

Compute a sparse set of landmarks to speed up the geodesic computations. The landmarks are computed using farthest point sampling.

First landmarks, at random.

landmarks = 23057;
Dland = [];

Perform Fast Marching to compute the geodesic distance, and record it.

[Dland(:,end+1),S,Q] = perform_fast_marching_mesh(vertex, faces, landmarks(end));

Select farthest point. Here, min(Dland,[],2) is the distance to the set of seed points.

[tmp,landmarks(end+1)] = max( min(Dland,[],2) );

Update distance function.

[Dland(:,end+1),S,Q] = perform_fast_marching_mesh(vertex, faces, landmarks(end));

Display distances.

clf;
options.start_points = landmarks;
plot_fast_marching_mesh(vertex,faces, min(Dland,[],2) , [], options);

Exercice 1: (check the solution) Compute a set of n = 300 vertex by iterating this farthest point sampling. Display the progression of the sampling.

exo1;

Compute the distance matrix restricted to the landmarks.

D = Dland(landmarks,:);
D = (D+D')/2;

Bending Invariant by Strain Minimization and Nistrom Interpolation

One can compute the bending invariant of the set of landmarks, and then apply it to the whole mesh using interpolation.

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';

Interpolate the locations to the whole mesh by Nystrom eigen-extrapolation, as detailed in

Sparse multidimensional scaling using landmark points V. de Silva, J.B. Tenenbaum, Preprint.

vertex1 = zeros(nverts,3);
deltan = mean(Dland.^2,1);
for i=1:nverts
    deltax = Dland(i,:).^2;
    vertex1(i,:) = 1/2 * ( Xstrain * ( deltan-deltax )' )';
end
vertex1 = vertex1';

Display the bending invariant mesh.

clf;
plot_mesh(vertex1,faces,options);

Farthest Point for Stress Minimization

The proposed interpolation method is valid only for the Strain minimizer (spectral Nistrom interpolation). One thus needs to use another interpolation method.

See for instance this work for a method to do such an interpolation:

A. M. Bronstein, M. M. Bronstein, R. Kimmel, Efficient computation of isometry-invariant distances between surfaces, SIAM J. Scientific Computing, Vol. 28/5, pp. 1812-1836, 2006.

Exercice 2: (check the solution) Create an interpolation scheme to interpolate the result of MDS dimensionality reduction with Stree minimization (SMACOF algorithm).

exo2;