Gabriel Peyré
January 01, 2011
920

# Mesh Processing Course: Geodesic Mesh Processing

January 01, 2011

## Transcript

2. ### Local vs. Global Processing 2 Local Processing Differential Computations Global

Processing Geodesic Computations Surface ﬁltering Fourier on Meshes Front Propagation on Meshes Surface Remeshing
3. ### Overview •Metrics and Riemannian Surfaces. • Geodesic Computation - Iterative

Scheme • Geodesic Computation - Fast Marching • Shape Recognition with Geodesic Statistics • Geodesic Meshing 3
4. ### Parametric Surfaces 4 Parameterized surface: u ⇥ R2 ⇤ (u)

⇥ M. u1 u2 ⇥ ⇥u1 ⇥ ⇥u2
5. ### Parametric Surfaces 4 Parameterized surface: u ⇥ R2 ⇤ (u)

⇥ M. Curve in parameter domain: t ⇥ [0, 1] ⇤ (t) ⇥ D. u1 u2 ⇥ ⇥u1 ⇥ ⇥u2
6. ### Parametric Surfaces 4 Parameterized surface: u ⇥ R2 ⇤ (u)

⇥ M. Curve in parameter domain: t ⇥ [0, 1] ⇤ (t) ⇥ D. Geometric realization: ¯(t) def. = ⇥( (t)) M. u1 u2 ⇥ ⇥u1 ⇥ ⇥u2 ¯ ¯
7. ### Parametric Surfaces 4 Parameterized surface: u ⇥ R2 ⇤ (u)

⇥ M. Curve in parameter domain: t ⇥ [0, 1] ⇤ (t) ⇥ D. Geometric realization: ¯(t) def. = ⇥( (t)) M. For an embedded manifold M Rn: First fundamental form: I = ⇥ ⇥ui , ⇥ ⇥uj ⇥ ⇥ i,j=1,2 . u1 u2 ⇥ ⇥u1 ⇥ ⇥u2 L( ) def. = 1 0 ||¯ (t)||dt = 1 0 ⇥ (t)I (t) (t)dt. Length of a curve ¯ ¯
8. ### Isometric and Conformal M is locally isometric to the plane:

I = Id. Exemple: M =cylinder. Surface not homeomorphic to a disk:
9. ### Isometric and Conformal M is locally isometric to the plane:

I = Id. Exemple: M =cylinder. ⇥ is conformal: I (u) = (u)Id. Exemple: stereographic mapping plane sphere. Surface not homeomorphic to a disk:
10. ### Riemannian Manifold 6 Length of a curve (t) M: L(

) def. = 1 0 ⇥ (t)TH( (t)) (t)dt. Riemannian manifold: M Rn (locally) Riemannian metric: H(x) Rn n, symmetric, positive deﬁnite.
11. ### Riemannian Manifold 6 Length of a curve (t) M: L(

) def. = 1 0 ⇥ (t)TH( (t)) (t)dt. W(x) Euclidean space: M = Rn, H(x) = Id n . Riemannian manifold: M Rn (locally) Riemannian metric: H(x) Rn n, symmetric, positive deﬁnite.
12. ### Riemannian Manifold 6 Length of a curve (t) M: L(

) def. = 1 0 ⇥ (t)TH( (t)) (t)dt. W(x) Euclidean space: M = Rn, H(x) = Id n . 2-D shape: M R2, H(x) = Id 2 . Riemannian manifold: M Rn (locally) Riemannian metric: H(x) Rn n, symmetric, positive deﬁnite.
13. ### Riemannian Manifold 6 Length of a curve (t) M: L(

) def. = 1 0 ⇥ (t)TH( (t)) (t)dt. W(x) Euclidean space: M = Rn, H(x) = Id n . 2-D shape: M R2, H(x) = Id 2 . Riemannian manifold: M Rn (locally) Riemannian metric: H(x) Rn n, symmetric, positive deﬁnite. Isotropic metric: H(x) = W(x)2Id n .
14. ### Riemannian Manifold 6 Length of a curve (t) M: L(

) def. = 1 0 ⇥ (t)TH( (t)) (t)dt. W(x) Euclidean space: M = Rn, H(x) = Id n . 2-D shape: M R2, H(x) = Id 2 . Image processing: image I, W(x)2 = ( + || I(x)||) 1. Riemannian manifold: M Rn (locally) Riemannian metric: H(x) Rn n, symmetric, positive deﬁnite. Isotropic metric: H(x) = W(x)2Id n .
15. ### Riemannian Manifold 6 Length of a curve (t) M: L(

) def. = 1 0 ⇥ (t)TH( (t)) (t)dt. W(x) Euclidean space: M = Rn, H(x) = Id n . 2-D shape: M R2, H(x) = Id 2 . Parametric surface: H(x) = Ix (1st fundamental form). Image processing: image I, W(x)2 = ( + || I(x)||) 1. Riemannian manifold: M Rn (locally) Riemannian metric: H(x) Rn n, symmetric, positive deﬁnite. Isotropic metric: H(x) = W(x)2Id n .
16. ### Riemannian Manifold 6 Length of a curve (t) M: L(

) def. = 1 0 ⇥ (t)TH( (t)) (t)dt. W(x) Euclidean space: M = Rn, H(x) = Id n . 2-D shape: M R2, H(x) = Id 2 . Parametric surface: H(x) = Ix (1st fundamental form). Image processing: image I, W(x)2 = ( + || I(x)||) 1. DTI imaging: M = [0, 1]3, H(x)=di usion tensor. Riemannian manifold: M Rn (locally) Riemannian metric: H(x) Rn n, symmetric, positive deﬁnite. Isotropic metric: H(x) = W(x)2Id n .
17. ### Geodesic Distances 7 Geodesic distance metric over M Rn Geodesic

curve: (t) such that L( ) = dM (x, y). Distance map to a starting point x0 M: Ux0 (x) def. = dM (x0, x). dM (x, y) = min (0)=x, (1)=y L( )
18. ### Geodesic Distances 7 Geodesic distance metric over M Rn Geodesic

curve: (t) such that L( ) = dM (x, y). Distance map to a starting point x0 M: Ux0 (x) def. = dM (x0, x). metric geodesics dM (x, y) = min (0)=x, (1)=y L( ) Euclidean
19. ### Geodesic Distances 7 Geodesic distance metric over M Rn Geodesic

curve: (t) such that L( ) = dM (x, y). Distance map to a starting point x0 M: Ux0 (x) def. = dM (x0, x). metric geodesics dM (x, y) = min (0)=x, (1)=y L( ) Euclidean Shape
20. ### Geodesic Distances 7 Geodesic distance metric over M Rn Geodesic

curve: (t) such that L( ) = dM (x, y). Distance map to a starting point x0 M: Ux0 (x) def. = dM (x0, x). metric geodesics dM (x, y) = min (0)=x, (1)=y L( ) Euclidean Shape Isotropic
21. ### Geodesic Distances 7 Geodesic distance metric over M Rn Geodesic

curve: (t) such that L( ) = dM (x, y). Distance map to a starting point x0 M: Ux0 (x) def. = dM (x0, x). metric geodesics dM (x, y) = min (0)=x, (1)=y L( ) Euclidean Shape Isotropic Anisotropic
22. ### Geodesic Distances 7 Geodesic distance metric over M Rn Geodesic

curve: (t) such that L( ) = dM (x, y). Distance map to a starting point x0 M: Ux0 (x) def. = dM (x0, x). metric geodesics dM (x, y) = min (0)=x, (1)=y L( ) Euclidean Shape Isotropic Anisotropic Surface
23. ### Anisotropy and Geodesics 8 H(x) = 1 (x)e1 (x)e1 (x)T

+ 2 (x)e2 (x)e2 (x)T with 0 < 1 2, Tensor eigen-decomposition: x e1 (x) M e2 (x) 2 (x) 1 2 1 (x) 1 2 { \ H(x) 1}
24. ### Anisotropy and Geodesics 8 H(x) = 1 (x)e1 (x)e1 (x)T

+ 2 (x)e2 (x)e2 (x)T with 0 < 1 2, Tensor eigen-decomposition: x e1 (x) M e2 (x) 2 (x) 1 2 1 (x) 1 2 { \ H(x) 1} Geodesics tend to follow e1 (x).
25. ### Anisotropy and Geodesics 8 H(x) = 1 (x)e1 (x)e1 (x)T

+ 2 (x)e2 (x)e2 (x)T with 0 < 1 2, Tensor eigen-decomposition: Local anisotropy of the metric: 4 ECCV-08 submission ID 1057 Figure 2 shows examples of geodesic curves computed from a single starting point S = {x1 } in the center of the image = [0, 1]2 and a set of points on the boundary of . The geodesics are computed for a metric H(x) whose anisotropy ⇥(x) (deﬁned in equation (2)) is increasing, thus making the Riemannian space progressively closer to the Euclidean space. Image f = .1 = .2 = .5 = 1 Image f = .5 = 0 = .95 = .7 (x) = ⇥1 (x) ⇥2 (x) ⇥1 (x) + ⇥2 (x) [0, 1] x e1 (x) M e2 (x) 2 (x) 1 2 1 (x) 1 2 { \ H(x) 1} Geodesics tend to follow e1 (x).
26. ### Isotropic Metric Design 9 Image f Metric W(x) Distance Ux0

(x) Geodesic curve (t) Image-based potential: H(x) = W(x)2Id 2 , W(x) = ( + |f(x) c|)
27. ### Isotropic Metric Design 9 Image f Metric W(x) Distance Ux0

(x) Geodesic curve (t) Image f Metric W(x) U{x0,x1} Image-based potential: H(x) = W(x)2Id 2 , W(x) = ( + |f(x) c|) Geodesics Gradient-based potential: W(x) = ( + || xf||)
28. ### Isotropic Metric Design: Vessels 10 f ˜ f W =

( + max( ˜ f, 0)) Remove background: ˜ f = G ⇥ f f, ⇥vessel width.
29. ### Isotropic Metric Design: Vessels 10 f ˜ f W =

( + max( ˜ f, 0)) 3D Volumetric datasets: Remove background: ˜ f = G ⇥ f f, ⇥vessel width.
30. ### Overview • Metrics and Riemannian Surfaces. •Geodesic Computation - Iterative

Scheme • Geodesic Computation - Fast Marching • Shape Recognition with Geodesic Statistics • Geodesic Meshing 11
31. ### Eikonal Equation and Viscosity Solution 12 U(x) = d(x0, x)

Distance map: Theorem: U is the unique viscosity solution of || U(x)||H(x) 1 = 1 with U(x0 ) = 0 where ||v||A = v Av
32. ### Eikonal Equation and Viscosity Solution 12 Geodesic curve between x1

and x0 solves (t) = ⇥t H( (t)) 1 Ux0 ( (t)) (0) = x1 t > 0 with U(x) = d(x0, x) Distance map: Theorem: U is the unique viscosity solution of || U(x)||H(x) 1 = 1 with U(x0 ) = 0 where ||v||A = v Av
33. ### Eikonal Equation and Viscosity Solution 12 Geodesic curve between x1

and x0 solves Example: isotropic metric H(x) = W(x)2Id n , (t) = ⇥t H( (t)) 1 Ux0 ( (t)) (0) = x1 t > 0 || U(x)|| = W(x) (t) = ⇥t U( (t)) and with U(x) = d(x0, x) Distance map: Theorem: U is the unique viscosity solution of || U(x)||H(x) 1 = 1 with U(x0 ) = 0 where ||v||A = v Av
34. ### Simpliﬁed Proof V solving ||⇤V (x)||2 H 1 = H

1(x)⇤V (x), ⇤V (x)⇥ = 1, V (x0 ) = 0. U(x) = min :x0 x L( ) = 1 0 H( (t)) (t), (t) dt
35. ### Simpliﬁed Proof V solving ||⇤V (x)||2 H 1 = H

1(x)⇤V (x), ⇤V (x)⇥ = 1, V (x0 ) = 0. U V , ⇤V ⇥ = H1/2 , H 1/2⇤V ⇥ ||H1/2 ||||H 1/2⇤V || C.S. = 1 If V is smooth on : Let : x0 x be any smooth curve. U(x) = min :x0 x L( ) = 1 0 H( (t)) (t), (t) dt
36. ### Simpliﬁed Proof V solving ||⇤V (x)||2 H 1 = H

1(x)⇤V (x), ⇤V (x)⇥ = 1, V (x0 ) = 0. U V , ⇤V ⇥ = H1/2 , H 1/2⇤V ⇥ ||H1/2 ||||H 1/2⇤V || = U(x) = min L( ) V (x) C.S. = 1 = 0 L( ) = 1 0 ||H1/2 || 1 0 ⇥ , ⌅V ⇤ = V ( (1)) V ( (0)) = V (x) If V is smooth on : Let : x0 x be any smooth curve. U(x) = min :x0 x L( ) = 1 0 H( (t)) (t), (t) dt
37. ### Simpliﬁed Proof (cont.) U V Deﬁne: (t) = H 1(

(t)) V ( (t)) x0 x (0) = x Let x be arbitrary.
38. ### Simpliﬁed Proof (cont.) U V Deﬁne: (t) = H 1(

(t)) V ( (t)) x0 x (0) = x Let x be arbitrary. dV ( (t)) dt = (t), V ( (t)) = 1 If V is smooth on ([0, tmax )), then = (tmax ) = x0
39. ### Simpliﬁed Proof (cont.) 14 13 U V Deﬁne: = 1

= 0 (t) = H 1( (t)) V ( (t)) x0 H , = H 1 V, V = 1 x (0) = x Let x be arbitrary. dV ( (t)) dt = (t), V ( (t)) = 1 If V is smooth on ([0, tmax )), then = (tmax ) = x0 = tmax 0 , V = V ( (tmax )) + V ( (0)) = V (x) One has: U(x) L( ) = tmax 0 H , = tmax 0 H ,
40. ### Discretization 15 x x0 B(x) Control (derivative-free) formulation: U(x) =

d(x0, x) is the unique solution of y U(x) = (U)(x) = min y B(x) U(y) + d(x, y)
41. ### Discretization 15 x x0 B(x) xj xi xk B(x) Control

(derivative-free) formulation: U(x) = d(x0, x) is the unique solution of Manifold discretization: triangular mesh. U discretization: linear ﬁnite elements. H discretization: constant on each triangle. y U(x) = (U)(x) = min y B(x) U(y) + d(x, y)
42. ### Discretization 15 x x0 B(x) xj xi xk B(x) Control

(derivative-free) formulation: U(x) = d(x0, x) is the unique solution of Manifold discretization: triangular mesh. U discretization: linear ﬁnite elements. H discretization: constant on each triangle. Ui = (U) i = min f=(i,j,k) Vi,j,k Vi,j,k = min 0 t 1 tUj + (1 t)Uk on regular grid: equivalent to upwind FD. explicit solution (solving quadratic equation). xj xk xi txj + (1 t)xk y +||txj + (1 t)xk xi ||Hijk U(x) = (U)(x) = min y B(x) U(y) + d(x, y)
43. ### Update Step on a triangulation 16 Vi,j,k = min 0

t 1 tUj + (1 t)Uk (U) i = min f=(i,j,k) Vi,j,k xi xj xk +||txj + (1 t)xk xi ||Hijk Discrete Eikonal equation:
44. ### Update Step on a triangulation 16 Vi,j,k = min 0

t 1 tUj + (1 t)Uk Distance function in (i, j, k): (U) i = min f=(i,j,k) Vi,j,k xi xj xk +||txj + (1 t)xk xi ||Hijk g Unknowns: Discrete Eikonal equation: = Vi,j,k U(x) = x xi, g + d gradient
45. ### Update Step on a triangulation 16 Vi,j,k = min 0

t 1 tUj + (1 t)Uk Distance function in (i, j, k): (U) i = min f=(i,j,k) Vi,j,k X = (xj xi, xk xi ) Rd 2 u = (Uj, Uk ) R2 S = (X X) 1 R2 2 I = (1, 1) R2 xi xj xk +||txj + (1 t)xk xi ||Hijk g Notations: Unknowns: Discrete Eikonal equation: = Vi,j,k U(x) = x xi, g + d gradient Hi,j,k = w2Id 3 (for simpliﬁty)
46. ### Update Step on a triangulation (cont.) 17 xi xj xk

0 X g + dI = u = = S(u dI) Find g = X , R2 and d = Vi,j,k .
47. ### Update Step on a triangulation (cont.) 17 xi xj xk

0 X g + dI = u = = S(u dI) Find g = X , R2 and d = Vi,j,k . || U(xi )|| = ||g|| = w Discrete Eikonal equation:
48. ### Update Step on a triangulation (cont.) 17 xi xj xk

0 X g + dI = u = d2 2bd + c = 0 a = SI, I b = SI, u c = Su, u w2 = = S(u dI) Quadratic equation: Find g = X , R2 and d = Vi,j,k . || U(xi )|| = ||g|| = w Discrete Eikonal equation: ||XS(u dI)||2 = w2
49. ### Update Step on a triangulation (cont.) 17 = b2 ac

d = b + a Admissible solution: dj = Uj + Wi ||xi xj || (ui ) = d if 0 min(dj, dk ) otherwise. xi xj xk 1 0 0 X g + dI = u = d2 2bd + c = 0 a = SI, I b = SI, u c = Su, u w2 = = S(u dI) Quadratic equation: Find g = X , R2 and d = Vi,j,k . || U(xi )|| = ||g|| = w Discrete Eikonal equation: ||XS(u dI)||2 = w2
50. ### Numerical Schemes 18 Fixed point equation: U = (U) is

monotone: U V = (U) (V ) U( +1) = (U( )) U( ) U solving (U) = U U(0) = 0, U( ) Iterative schemes: || (U( )) U( )|| = U( +1) U( ) C < +
51. ### Numerical Schemes 18 Fixed point equation: U = (U) is

monotone: U V = (U) (V ) U( +1) = (U( )) U( ) U solving (U) = U U(0) = 0, U( ) Iterative schemes: || (U( )) U( )|| Minimal path extraction: ( +1) = ( ) ⇥ H( ( )) 1 U( ( )) = U( +1) U( ) C < +

53. ### Discretization Errors 20 For a mesh with N points: U[N]

RN solution of (U[N]) = U[N] Linear interpolation: ˜ U[N](x) = i U[N] i i (x) Uniform convergence: || ˜ U[N] U|| N + ⇥ 0 Continuous geodesic distance U(x).
54. ### Discretization Errors 20 1 N i |UN i U(xi )|2

For a mesh with N points: U[N] RN solution of (U[N]) = U[N] Linear interpolation: ˜ U[N](x) = i U[N] i i (x) Uniform convergence: || ˜ U[N] U|| N + ⇥ 0 Numerical evaluation: Continuous geodesic distance U(x).
55. ### Overview • Metrics and Riemannian Surfaces. • Geodesic Computation -

Iterative Scheme •Geodesic Computation - Fast Marching • Shape Recognition with Geodesic Statistics • Geodesic Meshing 21
56. ### Causal Updates 22 j i, (U) i Uj Causality condition:

The value of Ui depends on {Uj }j with Uj Ui . Compute (U) i using an optimal ordering. Front propagation, O(N log(N)) operations.
57. ### Causal Updates 22 j i, (U) i Uj Causality condition:

The value of Ui depends on {Uj }j with Uj Ui . Compute (U) i using an optimal ordering. Front propagation, O(N log(N)) operations. u = (U) i is the solution of max(u Ui 1,j, u Ui+1,j, 0)2+ max(u Ui,j 1, u Ui,j+1, 0)2 = h2W2 i,j (upwind derivatives) Isotropic H(x) = W(x)2Id, square grid. xi+1,j xi,j+1 xi,j
58. ### Causal Updates 22 j i, (U) i Uj Causality condition:

The value of Ui depends on {Uj }j with Uj Ui . Compute (U) i using an optimal ordering. Front propagation, O(N log(N)) operations. triangulation with no obtuse angles. Bad Good u = (U) i is the solution of max(u Ui 1,j, u Ui+1,j, 0)2+ max(u Ui,j 1, u Ui,j+1, 0)2 = h2W2 i,j (upwind derivatives) Isotropic H(x) = W(x)2Id, square grid. Surface (ﬁrst fundamental form) xi xj xk Good Bad xi+1,j xi,j+1 xi,j
59. ### Front Propagation 23 x0 Algorithm: Far Front Computed. 2) Move

from Front to Computed . Iteration Front Ft , Ft = {i \ Ui t} Ft State Si {Computed, Front, Far} 3) Update Uj = (U) j for neighbors 1) Select front point with minimum Ui and

64. ### Overview • Metrics and Riemannian Surfaces. • Geodesic Computation -

Iterative Scheme • Geodesic Computation - Fast Marching •Shape Recognition with Geodesic Statistics • Geodesic Meshing 28

66. ### Bending Invariant Recognition 29 [Zoopraxiscope, 1876] Shape articulations: ˜ x1

˜ x2 M Surface bendings: [Elad, Kimmel, 2003]. [Bronstein et al., 2005].
67. ### 2D Shapes 30 2D shape: connected, closed compact set S

R2. Piecewise-smooth boundary S. Geodesic distance in S for uniform metric: dS (x, y) def. = min ⇥P(x,y) L( ) where L( ) def. = 1 0 | (t)|dt, Shape S Geodesics
68. ### Distribution of Geodesic Distances 31 0 20 40 60 80

0 20 40 60 80 0 20 40 60 80 Distribution of distances to a point x: {dM (x, y)}y M