Spectral Mesh Flattening

$\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.

In [2]:

Spectral Mesh 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.

In [3]:
name = 'nefertiti'; 
options.name = name;
[vertex,faces] = read_mesh(name);
n = size(vertex,2);

Display it.

In [4]:
plot_mesh(vertex,faces, options);
shading faceted;

Compute the mesh Laplacian matrix.

In [5]:
options.symmetrize = 1;
options.normalize = 0;
L = compute_mesh_laplacian(vertex,faces,'conformal',options);

Compute the eigenvalues and eigenvectors

In [6]:
[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.

In [7]:
vertexF = U(:,2:3)';

Use translation / rotation to align the parameterization.

In [8]:
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); ...

Display the flattened mesh.

In [9]:

Exercise 1

Perform the same flattening, but with the combinatorial Laplacian.

In [10]:
In [11]:
%% Insert your code here.

Geodesic Embedding (Isomap)

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.

In [12]:
D = zeros(n);
for i=1:n
    D(:,i) = perform_fast_marching_mesh(vertex,faces,i);

Enforce symmetry.

In [13]:
D = (D+D')/2;

Compute the centered matrix.

In [14]:
J = eye(n) - ones(n)/n;
W = -J*(D.^2)*J;

Diagonalize the centered matrix.

In [15]:
[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.

In [16]:
hh = plot(S(1:30), '.-'); axis('tight');
set(hh, 'LineWidth', 2);