$\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\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}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$
This tour explores 2D flattening of 3D surfaces using spectral methods.
addpath('toolbox_signal')
addpath('toolbox_general')
addpath('toolbox_graph')
addpath('solutions/meshdeform_3_flattening')
Mesh flattening finds 2D locations that minimize a variational energy with a non-degenracy constraint (for instance maximal variance).
For the Dirichlet energy (Sobolev norm), the resulting location are described by the first eigenvectors of the Laplacian.
This method is refered to as "Laplacian eigenmaps" in manifold learning, see:
The advantage over fixed boundary harmonic parameterization is that the boundary of the flattened domain is not fixed, but the drawback is that the parameterization is not guaranteed to be valid (bijective).
First load a mesh.
name = 'nefertiti';
options.name = name;
[vertex,faces] = read_mesh(name);
n = size(vertex,2);
Display it.
clf;
plot_mesh(vertex,faces, options);
shading faceted;
Compute the mesh Laplacian matrix.
options.symmetrize = 1;
options.normalize = 0;
L = compute_mesh_laplacian(vertex,faces,'conformal',options);
Compute the eigenvalues and eigenvectors
[U,S] = eig(full(L)); S = diag(S);
[S,I] = sort(S,'ascend'); U = U(:,I);
The vertex positions are the eigenvectors 2 and 3.
vertexF = U(:,2:3)';
Use translation / rotation to align the parameterization.
icenter = 88;
irotate = 154;
vertexF = vertexF - repmat(vertexF(:,icenter), [1 n]);
theta = -pi/2+atan2(vertexF(2,irotate),vertexF(1,irotate));
vertexF = [vertexF(1,:)*cos(theta)+vertexF(2,:)*sin(theta); ...
-vertexF(1,:)*sin(theta)+vertexF(2,:)*cos(theta)];
Display the flattened mesh.
clf;
plot_mesh(vertexF,faces);
Exercise 1
Perform the same flattening, but with the combinatorial Laplacian.
exo1()
%% Insert your code here.
Another (nonlinear) embedding can be computed by minimizing the geodesic distortion between points on the surface and points over the parameterized domain.
First we compute the geodesic distance on the mesh using the Fast Marching algorithm.
D = zeros(n);
for i=1:n
D(:,i) = perform_fast_marching_mesh(vertex,faces,i);
end
Enforce symmetry.
D = (D+D')/2;
Compute the centered matrix.
J = eye(n) - ones(n)/n;
W = -J*(D.^2)*J;
Diagonalize the centered matrix.
[U,S] = eig(W);
S = diag(S);
[S,I] = sort(S,'descend'); U = U(:,I);
Display the decay of the eigenvalues. If the mesh was isometric to the plane, then only the two largest eigenvalues would be non zero.
clf;
hh = plot(S(1:30), '.-'); axis('tight');
set(hh, 'LineWidth', 2);
Isomap embedding is defined from the two largest eigenvalues.
vertexF = U(:,1:2)' .* repmat(sqrt(S(1:2)), [1 n]);
Align the parameters.
vertexF = vertexF - repmat(vertexF(:,icenter), [1 n]);
theta = -pi/2+atan2(vertexF(2,irotate),vertexF(1,irotate));
vertexF = [vertexF(1,:)*cos(theta)+vertexF(2,:)*sin(theta); ...
-vertexF(1,:)*sin(theta)+vertexF(2,:)*cos(theta)];
Display and compare with Laplacian embeddedding.
clf;
plot_mesh(vertexF,faces,options);
Exercise 2
Compute the embedding using Stress minimization with SMACOF. See the numerical tours on bending invariants for more details.
exo2()
%% Insert your code here.
Plot stress evolution during minimization.
clf;
plot(stress(2:end), '.-');
axis('tight');
Exercise 3
Compute mesh parameterization using a circle as boundary. otan weights oundary ixed positions ystem et up the right hand sizes with the fixed position. olve lign isplay
exo3()
%% Insert your code here.