Important: Please read the installation page for details about how to install the toolboxes. $\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 parameterization of 3D surfaces onto a sphere.
We use a simple minimization of the Dirichlet energy under spherical constraints. There is no theoritical guarantee, but for some meshes, it seems to work correctly.
from __future__ import division
import nt_toolbox as nt
from nt_solutions import meshdeform_2_parameterization_sphere as solutions
%matplotlib inline
%load_ext autoreload
%autoreload 2
name = 'horse'
[vertex, faces] = read_mesh(name)
n = size(vertex, 2)
m = size(faces, 2)
clear options; options.name = name
Display the mesh.
options.lighting = 1
plot_mesh(vertex, faces, options)
shading faceted
Compute the weights. The weights should be positive for the method to work.
weight = 'conformal'
weight = 'combinatorial'
switch weight
case 'conformal'
W = make_sparse(n, n)
for i in 1: 3:
i1 = mod(i-1, 3) + 1
i2 = mod(i , 3) + 1
i3 = mod(i + 1, 3) + 1
pp = vertex(: , faces(i2, : )) - vertex(: , faces(i1, : ))
qq = vertex(: , faces(i3, : )) - vertex(: , faces(i1, : ))
% normalize the vectors
pp = pp ./ repmat(sqrt(sum(pp.^2, 1)), [3 1])
qq = qq ./ repmat(sqrt(sum(qq.^2, 1)), [3 1])
% compute angles
ang = acos(sum(pp.*qq, 1))
a = max(1 ./ tan(ang), 1e-1); % this is *very* important
W = W + make_sparse(faces(i2, : ), faces(i3, : ), a, n, n)
W = W + make_sparse(faces(i3, : ), faces(i2, : ), a, n, n)
case 'combinatorial'
E = [faces([1 2], : ) faces([2 3], : ) faces([3 1], : )]
E = unique_rows([E E(2: -1: 1, : )]')'
W = make_sparse(E(1, : ), E(2, : ), ones(size(E, 2), 1))
Compute the normalized weight matrix |tW| such that its rows sums to 1.
d = full(sum(W, 1))
D = spdiags(d(: ), 0, n, n)
iD = spdiags(d(: ).^(-1), 0, n, n)
tW = iD * W
It is possible to smooth the positions of the mesh on the sphere by performing an averaging according to |W|, and projecting back on the sphere.
Compute an initial mapping on the sphere. This simply a radial projection.
vertex1 = vertex
vertex1 = vertex1 - repmat(mean(vertex1, 2), [1 n])
vertex1 = vertex1 ./ repmat(sqrt(sum(vertex1.^2, 1)), [3 1])
Check which faces have the correct orientation.
normal to faces
[normal, normalf] = compute_normal(vertex1, faces)
center of faces
C = squeeze(mean(reshape(vertex1(: , faces), [3 3 m]), 2))
inner product
I = sum(C.*normalf)
Ratio of inverted triangles. For a bijective mapping, there should not be any inverted triangle.
disp(['Ratio of inverted triangles: ' num2str(sum(I(: ) <0)/ m, 3) '%'])
Display on the sphere.
options.name = 'none'
options.face_vertex_color = double(I(: ) >0)
plot_mesh(vertex1, faces, options)
colormap gray(256); axis tight
shading faceted
Perform smoothing and projection.
vertex1 = vertex1*tW'
vertex1 = vertex1 ./ repmat(sqrt(sum(vertex1.^2, 1)), [3 1])
Exercise 1
Perform iterative smoothing and projection. Record the evolution of the number of inverted triangle in |ninvert|. Record also the evolution of the Dirichlet energy in |Edir|.
solutions.exo1()
## Insert your code here.
Display the decay of the evolution of the Dirichlet energy.
plot(Edir/ Edir(1))
axis('tight')
Display final spherical configuration.
plot_mesh(vertex1, faces)
colormap gray(256); axis tight
shading faceted
Display the evolution of the number of inverted triangles.
plot(ninvert/ m); axis tight
Using this spherical parameterization, one maps the surface on a sphere, then on an octahedron, and finally on a square. This allows to map the surface on a 2D image, thus creating a geometry image.
The method is originaly described in
Spherical Parameterization and Remeshing E. Praun, H. Hoppe Proceedings of SIGGRAPH 2003
q = 128
options.verb = 0
M = perform_sgim_sampling(vertex, vertex1, faces, q, options)
Display the spherical geometry image.
plot_geometry_image(M, 1, 1)
axis equal
colormap gray(256)
view(134, -61)
By mapping two meshes on the same sphere, one computes a bijection between two meshes.
By linearly interpolating the positions of the points that are in correspondance, one performs a warp of a mesh onto another one.
Exercise 2
Implement the mesh morphing.
solutions.exo2()
## Insert your code here.
Spherical relaxation leads to an uncontrolled evolution because triangle are not constrained in size.
To avoid this, it is possible to penalize the size of large triangle.
This is similar to the method proposed in:
Unconstrained Spherical Parameterization Ilja Friedel, Peter Schr der, and Mathieu Desbrun Journal of Graphics Tools, 12(1), pp. 17-26, 2007.
First initialize the gradient descent.
vertex1 = vertex
vertex1 = vertex1 - repmat(mean(vertex1, 2), [1 n])
vertex1 = vertex1 ./ repmat(sqrt(sum(vertex1.^2, 1)), [3 1])
Step size for the gradient descent.
eta = .5
Compute the center of the faces.
A = squeeze(mean(reshape(vertex1(: , faces), [3 3 m]), 2))
Compute the Dirichlet energy of each face.
E = zeros(1, m)
for i in 1: 3:
i1 = mod(i, 3) + 1
% directed edge
u = vertex(: , faces(i, : )) - vertex(: , faces(i1, : ))
% norm squared
u = sum(u.^2)
% weights between the vertices
w = W(faces(i, : ) + (faces(i1, : )-1)*n)
E = E + w.*u
Compute gradient direction.
G = zeros(3, n)
for j in 1: m:
f = faces(: , j)
Alpha = A(: , j)
alpha = norm(Alpha)
for i in 1: 3:
i1 = mod(i , 3) + 1
i2 = mod(i + 1, 3) + 1
% directed edges
u1 = vertex(: , f(i)) - vertex(: , f(i1))
u2 = vertex(: , f(i)) - vertex(: , f(i2))
% weights between the vertices
w1 = W(f(i) + (f(i1)-1)*n)
w2 = W(f(i) + (f(i2)-1)*n)
G(: , f(i)) = G(: , f(i)) + (w1*u1 + w2*u2) ./ alpha^2 - Alpha/ alpha^4 * E(j)
Perform the gradient descent step and the projection.
vertex1 = vertex1 - eta*G
vertex1 = vertex1 ./ repmat(sqrt(sum(vertex1.^2, 1)), [3 1])
Exercise 3
Perform the full descent. Record the decay of the energy in |Edir|.
solutions.exo3()
## Insert your code here.
Plot the decay of the energy.
plot(Edir/ Edir(1))