Mesh Processing Course: Geodesic Mesh Processing

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