610

# 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
69. None
70. None
71. None
72. None
73. None
74. None
75. None
76. None
77. None
78. None
79. None
80. None
81. None
82. None
83. None
84. None
85. None
86. None
87. None
88. None
89. None