Parameterization and Flattening Gabriel Peyré CEREMADE, Université Paris-Dauphine www.numerical-tours.com

parameterization Mesh Parameterization - Overview 2 texture mapping M R3

parameterization D R2 Mesh Parameterization - Overview 2 1 re-sampling zoom texture mapping M R3

Overview •Dirichlet Energy on Meshes • Harmonic Parameterization • Spectral Flattening • Barycentric Coordinates for Warping • Geodesic Flattening • High Dimensional Data Analysis 3

Local Averaging 4 Local operator: W = (wij ) i,j V where wij = > 0 if j Vi, 0 otherwise. (Wf) i = (i,j) E wijfj. Examples: for i j, wij = 1 combinatorial wij = 1 ||xj xi ||2 distance wij = cot( ij ) + cot(⇥ij ) conformal Local averaging operator ˜ W = ( ˜ wij ) i,j V : ⇥ (i, j) E, ˜ wij = wij (i,j) E wij . ˜ W = D 1W with D = diag i (di ) where di = (i,j)⇥E wij. Averaging: ˜ W1 = 1. (explanations later)

Voronoi and Dual Mesh 5 Definition for a planar triangulation M of a mesh M R2. Voronoi for vertices: ⇧ i ⇤ V, Ei = {x ⇤ M \ ⇧ j ⌅= i, ||x xi || ⇥ ||x xj ||} Voronoi for edges: ⌅ e = (i, j) ⇥ E, Ee = {x ⇥ M \ ⌅ e ⇤= e, d(x, e) d(x, e )} Partition of the mesh: M = i V Ei = e E Ee. i j Ai cf i j A(ij) cf Dual mesh 1:3 subdivided mesh

Approximating Integrals on Meshes Approximation of integrals on vertices and edges: ⇥ M f(x)dx i V Ai f(xi ) e=(i,j) E Ae f([xi, xj ]). Theorem : ⇥ e = (i, j) E, Ae = Area(Ee ) = 1 2 ||xi xj ||2 (cot( ij ) + cot(⇥ij )) i j A(ij) cf

Approximating Integrals on Meshes Approximation of integrals on vertices and edges: ⇥ M f(x)dx i V Ai f(xi ) e=(i,j) E Ae f([xi, xj ]). Theorem : ⇥ e = (i, j) E, Ae = Area(Ee ) = 1 2 ||xi xj ||2 (cot( ij ) + cot(⇥ij )) i j A(ij) cf A B C O h + + = 2 A(ABO) = ||AB|| h = ||AB|| ||AB|| 2 tan( ) A(ABO) = ||AB||2 2 tan ⇤ 2 ( + ⇥) Proof:

Cotangent Weights Sobolev norm (Dirichlet energy): J(f) = M || f(x)||2dx

⇧ (i, j) ⇤ E, wij = 1. Distance weights: they depends both on the geometry and the to require faces information, ⇧ (i, j) ⇤ E, wij = 1 ||xj xi ||2 . Conformal weights: they depends on the full geometrical realiza require the face information ⇧ (i, j) ⇤ E, wij = cot( ij ) + cot(⇥ Figure 1.2 shows the geometrical meaning of the angles ij and ij = ⇥(xi, xj, xk1 ) and ⇥ij = ⇥(xi where (i, j, k1 ) ⇤ F and (i, j, k2 ) ⇤ F are the two faces adjacent in the next section the explanation of these celebrated cotangent xi xj xk1 xk2 ij ij Cotangent Weights 7 Approximation of Dirichelet energy: where wij = cot( ij ) + cot(⇥ij ). ⇥ M ||⇤xf||2dx ⇥ e E Ae |(Gf) e |2 = (i,j) E Ae |f(xj ) f(xi )|2 ||xj xi ||2 = (i,j) E wij |f(xj ) f(xi )|2 Sobolev norm (Dirichlet energy): J(f) = M || f(x)||2dx

⇧ (i, j) ⇤ E, wij = 1. Distance weights: they depends both on the geometry and the to require faces information, ⇧ (i, j) ⇤ E, wij = 1 ||xj xi ||2 . Conformal weights: they depends on the full geometrical realiza require the face information ⇧ (i, j) ⇤ E, wij = cot( ij ) + cot(⇥ Figure 1.2 shows the geometrical meaning of the angles ij and ij = ⇥(xi, xj, xk1 ) and ⇥ij = ⇥(xi where (i, j, k1 ) ⇤ F and (i, j, k2 ) ⇤ F are the two faces adjacent in the next section the explanation of these celebrated cotangent xi xj xk1 xk2 ij ij Cotangent Weights 7 Approximation of Dirichelet energy: Theorem : wij > 0 ⇥ ij + ⇥ij < ⇤ where wij = cot( ij ) + cot(⇥ij ). ⇥ M ||⇤xf||2dx ⇥ e E Ae |(Gf) e |2 = (i,j) E Ae |f(xj ) f(xi )|2 ||xj xi ||2 = (i,j) E wij |f(xj ) f(xi )|2 Sobolev norm (Dirichlet energy): J(f) = M || f(x)||2dx

Overview • Dirichlet Energy on Meshes •Harmonic Parameterization • Spectral Flattening • Barycentric Coordinates for Warping • Geodesic Flattening • High Dimensional Data Analysis 8

Mesh Parameterization 3D space (x,y,z) 2D parameter domain (u,v) boundary boundary 9 Parameterization: bijection : M ⇤ D ⇥ R2. Hypothesis: = ( 1, 2 ) is smooth, minimizes ⇥ xi ⇥M, (i) = 0(i) ⇥D. min 0 (i,j) E wi,j (| 1 (i) 1 (j)|2 + | 2 (i) 2 (j)|2) With boundary conditions 0:

Mesh Parameterization 3D space (x,y,z) 2D parameter domain (u,v) boundary boundary 9 Parameterization: bijection : M ⇤ D ⇥ R2. ⇥ sparse linear system to solve. Hypothesis: = ( 1, 2 ) is smooth, minimizes ⇥ xi ⇥M, (i) = 0(i) ⇥D. min 0 (i,j) E wi,j (| 1 (i) 1 (j)|2 + | 2 (i) 2 (j)|2) ⇥ i / ⇥M, (L 1 )(i) = (L 2 )(i) = 0 ⇥ i ⇥M, (i) = 0(i) ⇥D. Optimality conditions: With boundary conditions 0:

Mesh Parameterization 3D space (x,y,z) 2D parameter domain (u,v) boundary boundary 9 Parameterization: bijection : M ⇤ D ⇥ R2. ⇥ sparse linear system to solve. Theorem: (Tutte) if i, j, wij > 0, then is a bijection. ⇥ i / ⇥M, (i) = (i,j) E ˜ wi,j (j). Hypothesis: = ( 1, 2 ) is smooth, minimizes ⇥ xi ⇥M, (i) = 0(i) ⇥D. min 0 (i,j) E wi,j (| 1 (i) 1 (j)|2 + | 2 (i) 2 (j)|2) ⇥ i / ⇥M, (L 1 )(i) = (L 2 )(i) = 0 ⇥ i ⇥M, (i) = 0(i) ⇥D. Optimality conditions: Remark: each point is the average of its neighbors: With boundary conditions 0:

Examples of Parameterization Combinatorial Conformal Mesh 10

Examples of Parameterization Combinatorial Conformal Mesh 11

Application to Remeshing 12 parameterization 1 zoom re-sampling P. Alliez et al., Isotropic Surface Remeshing, 2003.

Application to Texture Mapping 13 param eterization texture g(u) color g( (x))

Mesh Parameterization #1 14 S S0 : S0 ⇥ S S0 ⇥ S ⇥ x S0 \⇥S0, = 0

Mesh Parameterization #1 14 W = make_sparse(n,n); for i=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 a = 1 ./ tan( acos(sum(pp.*qq,1)) ); a = max(a, 1e-2); % avoid degeneracy W = W + make_sparse(faces(i2,:),faces(i3,:), a, n, n ); W = W + make_sparse(faces(i3,:),faces(i2,:), a, n, n ); end S S0 : S0 ⇥ S S0 ⇥ S ⇥ x S0 \⇥S0, = 0 Discretization of :

Mesh Parameterization #1 14 W = make_sparse(n,n); for i=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 a = 1 ./ tan( acos(sum(pp.*qq,1)) ); a = max(a, 1e-2); % avoid degeneracy W = W + make_sparse(faces(i2,:),faces(i3,:), a, n, n ); W = W + make_sparse(faces(i3,:),faces(i2,:), a, n, n ); end D = spdiags(full( sum(W,1) ), 0, n,n); L = D - W; L1 = L; L1(boundary,:) = 0; L1(boundary + (boundary-1)*n) = 1; S S0 : S0 ⇥ S S0 ⇥ S ⇥ x S0 \⇥S0, = 0 Formation of the linear system: Discretization of :

Mesh Parameterization #1 14 W = make_sparse(n,n); for i=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 a = 1 ./ tan( acos(sum(pp.*qq,1)) ); a = max(a, 1e-2); % avoid degeneracy W = W + make_sparse(faces(i2,:),faces(i3,:), a, n, n ); W = W + make_sparse(faces(i3,:),faces(i2,:), a, n, n ); end D = spdiags(full( sum(W,1) ), 0, n,n); L = D - W; L1 = L; L1(boundary,:) = 0; L1(boundary + (boundary-1)*n) = 1; Rx = zeros(n,1); Rx(boundary) = x0; Ry = zeros(n,1); Ry(boundary) = y0; x = L1 \ Rx; y = L1 \ Ry; S S0 : S0 ⇥ S S0 ⇥ S ⇥ x S0 \⇥S0, = 0 Formation of the linear system: Formation of the RHS and resolution: Discretization of :

Mesh Parameterization #2 15 Exercise: perform the linear interpolation of the parameterization. Exercise: display the geometry image using a checkboard texture. Geometry image: re-sample X/Y/Z coordinates of on a grid. store the surface as a color (R/G/B) image.

Mesh Parameterization #3 16 Exercise: Locate the position of the eyes / the mouth Exercise: Compute an a ne transformation to re-align the texture. in the texture and on the mesh.

Mesh Deformations 17 Initial position: xi R3. Displacement of anchors: i I, xi xi = xi + i R3 Linear deformation: xi xi I i I, (i) = i i / I, (i) = 0 xi = xi + (i)

Mesh Deformations 17 Initial position: xi R3. Displacement of anchors: i I, xi xi = xi + i R3 Linear deformation: xi xi I % modify Laplacian L1 = L; L1(I,:) = 0; L1(I + (I-1)*n) = 1; % displace vertices vertex = vertex + ( L1 \ Delta0' )'; i I, (i) = i i / I, (i) = 0 xi = xi + (i)

Mesh Deformations 17 Initial position: xi R3. Displacement of anchors: i I, xi xi = xi + i R3 Linear deformation: xi xi xi Non-linear deformation: xi = ˜ xi + i ˜ xi coarse scale details ˜ xi ˜ xi Linear deformation: xi = xi + i, ni ˜ ni Extrusion along normals: ˜ ni I % modify Laplacian L1 = L; L1(I,:) = 0; L1(I + (I-1)*n) = 1; % displace vertices vertex = vertex + ( L1 \ Delta0' )'; i I, (i) = i i / I, (i) = 0 xi = xi + (i)

Overview • Dirichlet Energy on Meshes • Harmonic Parameterization •Spectral Flattening • Barycentric Coordinates for Warping • Geodesic Flattening • High Dimensional Data Analysis 18

Mesh Flattening 19 min 1⇥Rn ||G 1 ||2 = ⇥ i j wi,j | 1 (i) 1 (j)|2 with || 1 || = 1, ⇥ 1, 1⇤ = 0. min 2⇥Rn ||G 2 ||2 = ⇧ i j wi,j | 2 (i) 2 (j)|2 with ⌅ ⇤ ⌅ ⇥ || 2 || = 1, ⇥ 2, 1 ⇤ = 0, ⇥ 2, 1⇤ = 0. No boundary condition, minimize:

Mesh Flattening 19 min 1⇥Rn ||G 1 ||2 = ⇥ i j wi,j | 1 (i) 1 (j)|2 with || 1 || = 1, ⇥ 1, 1⇤ = 0. min 2⇥Rn ||G 2 ||2 = ⇧ i j wi,j | 2 (i) 2 (j)|2 with ⌅ ⇤ ⌅ ⇥ || 2 || = 1, ⇥ 2, 1 ⇤ = 0, ⇥ 2, 1⇤ = 0. ( 1 (i), 2 (i)) R2 Theorem: ⇥i = i+1L⇥i , where 0 = 0 1 2 . . . n 1 are eigenvalues of L = G⇥G. No boundary condition, minimize: conformal combinatorial

Proof 20 L = G G = D W = U U Spectral decomposition: = diag( i ) where 0 = 1 < 2 . . . n U = (ui )n i=1 orthonormal basis of Rn. u1 = 1

Proof 20 L = G G = D W = U U Spectral decomposition: = diag( i ) where 0 = 1 < 2 . . . n U = (ui )n i=1 orthonormal basis of Rn. If , 1 = 0, then u1 = 1 E(⇥) = ||G⇥||2 = n i=1 i | ⇥, ui ⇥|2 E(⇥) = n i=2 iai where ai = | ⇥, ui ⇥|2

Proof 20 L = G G = D W = U U Spectral decomposition: = diag( i ) where 0 = 1 < 2 . . . n U = (ui )n i=1 orthonormal basis of Rn. If , 1 = 0, then u1 = 1 E(⇥) = ||G⇥||2 = n i=1 i | ⇥, ui ⇥|2 Constrained minimization: linear program: minimum reached at a = i . min P n i=2 ai=1 n i=2 iai ±u2 = argmin , 1⇥=0,|| || E( ) E(⇥) = n i=2 iai where ai = | ⇥, ui ⇥|2

Flattening Examples 21 Main issue: No guarantee of being valid (bijective). conformal combinatorial

Overview • Dirichlet Energy on Meshes • Harmonic Parameterization • Spectral Flattening •Barycentric Coordinates for Warping • Geodesic Flattening • High Dimensional Data Analysis 22

Barycentric Coordinates 23 x1 x2 x3 x ˜ i (x) = A(x, xi+1, xi+2 ) A(x1, x2, x3 )

Barycentric Coordinates 23 x1 x2 x3 x1 x3 x2 x x Barycentric coordinates: { i (x)}i I Positivity: i (x) 0. Interpolation: Reproduction of a ne functions: ˜ i (x) = A(x, xi+1, xi+2 ) A(x1, x2, x3 ) Normalized: ˜ i (x) = i (x) j j (x) i I ˜ i (x)xi = x ˜ ⇥i (xj ) = i,j

Barycentric Coordinates 23 x1 x2 x3 x1 x3 x2 x x Barycentric coordinates: { i (x)}i I Positivity: i (x) 0. Interpolation: Reproduction of a ne functions: Application: interpolation of data {fi }i I ˜ i (x) = A(x, xi+1, xi+2 ) A(x1, x2, x3 ) Normalized: ˜ i (x) = i (x) j j (x) i I ˜ i (x)xi = x f(x) = i I ˜ i (x)fi ˜ ⇥i (xj ) = i,j

Barycentric Coordinates 23 x1 x2 x3 x1 x3 x2 x x Barycentric coordinates: { i (x)}i I Positivity: i (x) 0. Interpolation: Reproduction of a ne functions: Application: interpolation of data {fi }i I Application: mesh parameterization: xi xj ˜ i (x) = A(x, xi+1, xi+2 ) A(x1, x2, x3 ) Normalized: ˜ i (x) = i (x) j j (x) i I ˜ i (x)xi = x f(x) = i I ˜ i (x)fi wi,j = i (xj ) ˜ ⇥i (xj ) = i,j

Mean Value Coordinates 24 x xi+1 xi xi i ˜ i ˜ i i Conformal Laplacian weights: ⇥i (x) = cotan( i (x)) + cotan(˜ i (x)) Mean-value coordinates: not necessarily positive. valid coordinates. extend to non-convex coordinates (oriented angles). ⇥i (x) = tan( i (x)/2) + tan(˜ i (x)/2) ||x xi || ˜ 1 (x) ˜ 2 (x)

Barycentric Coordinates for Warping 25 Example: textured grid, 3D model, etc. Data points: {yj }j J C. Cage C: polygon with vertices {xi }i I . Initialization: data anchoring, compute j J, i I, i,j = i (yj ). Data warping: yj yj = i I i,jxi Satisfies: yi = i I i,jxi Cage warping: xi xi xi xi yi yi

Harmonic Coordinates 26 Mean value Harmonic Mean value coordinates: “non-physical” behavior, passes “through” the cage. ⇥ x C, i (x) = 0. ⇥ x ⇥C, i (x) = 0 i (x). Harmonic mapping: Boundary conditions:

Warping Comparison 27 Mean value Harmonic Initial shape

Overview • Dirichlet Energy on Meshes • Harmonic Parameterization • Spectral Flattening • Barycentric Coordinates for Warping •Geodesic Flattening • High Dimensional Data Analysis 28

Geodesic Distances 29 Length of a curve (t) M: L( ) def. = 1 0 W( (t))|| (t)||dt. Geodesic distance: dM (x, y) = min (0)=x, (1)=y L( ) Euclidean Shape Isotropic W = 1 Surface dM (x, y) = L( ) Geodesic curve :

Computation of Geodesic Distances 30 Non-linear PDE: Distance map to a point: Ux0 (x) = dM (x0, x). Ux0 x0 || Ux0 (x)|| = W(x) Ux0 (x0 ) = 0, (viscosity) Upwind finite di erences approximation. Fast Marching: front propagation in O(N log(N)) operations.

Manifold Flattening 31 Input manifold M, dM geodesic distance on M. Input geodesic distance matrix: ˜ D = (dM (˜ xi, ˜ xj )2) i,j , for ˜ xi M. x1 ˜ x1 x2 ˜ x2 M R3 M R2 ˜ x1 ˜ x2 x1 x2 Flattening: find X = (xi )p i=1 ⇤ Rn p such that ||xi xj || ⇥ dM (˜ xi, ˜ xj ). Surface parameterization Bending invariant

Stress Minimization 32 Geodesic stress: di,j = dM (˜ xi, ˜ xj ) S(X) = i,j |||xi xj || di,j |2, SMACOF algorithm: X( +1) = 1 N X( )B(X( )) where B(X) i,j = di,j ||xi xj || Non-convex functional: X( ) X local minimizer of S.

Projection on Distance Matrices 33 D(X) i,j = ||xi xj ||2 Di,j = d2 i,j ||xi xj ||2 = ||xi ||2 + ||xj ||2 2⇥xi, xj ⇤ min X1=0 i,j |||xi xj ||2 d2 i,j |2 = ||D(X) D||2, =⇥ D(X) = d1T + 1d 2XTX where d = (||xi ||2) i ⇤ Rn

Projection on Distance Matrices 33 D(X) i,j = ||xi xj ||2 Di,j = d2 i,j ||xi xj ||2 = ||xi ||2 + ||xj ||2 2⇥xi, xj ⇤ min X1=0 i,j |||xi xj ||2 d2 i,j |2 = ||D(X) D||2, =⇥ D(X) = d1T + 1d 2XTX where d = (||xi ||2) i ⇤ Rn For centered points: Centering matrix: J = Id n 11T /N 1 2JD(X)J = XTX JX = X J1 = 0

Projection on Distance Matrices 33 D(X) i,j = ||xi xj ||2 Di,j = d2 i,j ||xi xj ||2 = ||xi ||2 + ||xj ||2 2⇥xi, xj ⇤ min X1=0 i,j |||xi xj ||2 d2 i,j |2 = ||D(X) D||2, =⇥ D(X) = d1T + 1d 2XTX where d = (||xi ||2) i ⇤ Rn For centered points: Centering matrix: J = Id n 11T /N 1 2JD(X)J = XTX JX = X J1 = 0 Replace ||D(X) D|| by Explicit solution: diagonalize 1 2 JDJ = U UT i i 1 k = diag( 0, . . . , k 1 ), Uk = (u0, . . . , uk 1 )T, X = kUk min X || J(D(X) D)J/2|| = ||XTX + JDJ/2||

Isomap vs. Laplacian 34 Flattening: f = (f1, f2 ) R2. Isomap: global constraints: ||f(x) f(y)|| ⇥ dM (x, y). ⇥ (f1, f2 ) eigenvectors (#2,#3) of L = GTG. Laplacian: local smoothness: fi = argmin ||Gf|| subj. to ||f|| = 1. Mesh Lapl. combin. Lapl. conformal Isomap Bijective Not bijective ⇥ (f1, f2 ) eigenvectors (#1,#2) of J(dM (xi, xj )2) ijJ.

Bending Invariants of Surfaces 35 Bending invariants: [Elad, Kimmel, 2003]. Surface M, Isomap dimension reduction: x ⇥ M ⇤ IM (x) ⇥ R3. [Elad, Kimmel, 2003]. M IM

Bending Invariants of Surfaces 35 Bending invariants: [Elad, Kimmel, 2003]. Surface M, Isomap dimension reduction: x ⇥ M ⇤ IM (x) ⇥ R3. Geodesic isometry : M M : dM (x, y) = dM ( (x), (y)). Theorem: up to rigid motion, IM is invariant to geodesic isometries: IM (x) = v + UIM ( (x)) where U O(3) and v R3. [Elad, Kimmel, 2003]. M IM

Bending Invariants of Surfaces 35 Bending invariants: [Elad, Kimmel, 2003]. Surface M, Isomap dimension reduction: x ⇥ M ⇤ IM (x) ⇥ R3. Geodesic isometry : M M : dM (x, y) = dM ( (x), (y)). Theorem: up to rigid motion, IM is invariant to geodesic isometries: IM (x) = v + UIM ( (x)) where U O(3) and v R3. [Elad, Kimmel, 2003]. [Bronstein et al., 2005]. M IM

Bending Invariants 36 IM M

Face Recognition 37 Rigid similarity Non-rigid similarity Alex A. M. Bronstein et al., IJCV, 2005

Overview • Dirichlet Energy on Meshes • Harmonic Parameterization • Spectral Flattening • Barycentric Coordinates for Warping • Geodesic Flattening •High Dimensional Data Analysis 38

High Dimensional Data Sets 39

Graph and Geodesics 40

Isomap Dimension Reduction 41

Isomap vs PCA Flattening 42

Laplacian Spectral Dimension Reduction 43

Parameterization of Image Datasets 44

Library of Images 45

When Does it Works? 46

Local patches in images 47