Pro Yearly is on sale from $80 to $50! »

Mesh Processing Course: Geodesic Mesh Processing

Mesh Processing Course: Geodesic Mesh Processing

E34ded36efe4b7abb12510d4e525fee8?s=128

Gabriel Peyré

January 01, 2011
Tweet

Transcript

  1. Geodesic Data Processing Gabriel Peyré CEREMADE, Université Paris-Dauphine www.numerical-tours.com

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

    Processing Geodesic Computations Surface filtering 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 definite.
  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 definite.
  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 definite.
  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 definite. 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 definite. 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 definite. 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 definite. 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) (defined 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. Simplified 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. Simplified 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. Simplified 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. Simplified Proof (cont.) U V Define: (t) = H 1(

    (t)) V ( (t)) x0 x (0) = x Let x be arbitrary.
  38. Simplified Proof (cont.) U V Define: (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. Simplified Proof (cont.) 14 13 U V Define: = 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 finite 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 finite 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 simplifity)
  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 < +
  52. Numerical Examples on Meshes 19

  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 (first 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
  60. Fast Marching on an Image 24

  61. Fast Marching on Shapes and Surfaces 25

  62. Volumetric Datasets 26

  63. Propagation in 3D 27

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

    Iterative Scheme • Geodesic Computation - Fast Marching •Shape Recognition with Geodesic Statistics • Geodesic Meshing 28
  65. Bending Invariant Recognition 29 [Zoopraxiscope, 1876] Shape articulations:

  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