Mesh Processing Course: Geodesic Mesh Processing

Slide 1

Slide 1 text

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

Slide 2

Slide 2 text

Local vs. Global Processing 2 Local Processing Differential Computations Global Processing Geodesic Computations Surface ﬁltering Fourier on Meshes Front Propagation on Meshes Surface Remeshing

Slide 3

Slide 3 text

Overview •Metrics and Riemannian Surfaces. • Geodesic Computation - Iterative Scheme • Geodesic Computation - Fast Marching • Shape Recognition with Geodesic Statistics • Geodesic Meshing 3

Slide 4

Slide 4 text

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

Slide 5

Slide 5 text

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

Slide 6

Slide 6 text

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 ¯ ¯

Slide 7

Slide 7 text

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 ¯ ¯

Slide 8

Slide 8 text

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

Slide 9

Slide 9 text

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:

Slide 10

Slide 10 text

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.

Slide 11

Slide 11 text

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.

Slide 12

Slide 12 text

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.

Slide 13

Slide 13 text

Slide 14

Slide 14 text

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 .

Slide 15

Slide 15 text

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 .

Slide 16

Slide 16 text

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 .

Slide 17

Slide 17 text

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( )

Slide 18

Slide 18 text

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

Slide 19

Slide 19 text

Slide 20

Slide 20 text

Slide 21

Slide 21 text

Slide 22

Slide 22 text

Slide 23

Slide 23 text

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}

Slide 24

Slide 24 text

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).

Slide 25

Slide 25 text

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).

Slide 26

Slide 26 text

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|)

Slide 27

Slide 27 text

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||)

Slide 28

Slide 28 text

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

Slide 29

Slide 29 text

Isotropic Metric Design: Vessels 10 f ˜ f W = ( + max( ˜ f, 0)) 3D Volumetric datasets: Remove background: ˜ f = G ⇥ f f, ⇥vessel width.

Slide 30

Slide 30 text

Overview • Metrics and Riemannian Surfaces. •Geodesic Computation - Iterative Scheme • Geodesic Computation - Fast Marching • Shape Recognition with Geodesic Statistics • Geodesic Meshing 11

Slide 31

Slide 31 text

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

Slide 32

Slide 32 text

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

Slide 33

Slide 33 text

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

Slide 34

Slide 34 text

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

Slide 35

Slide 35 text

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

Slide 36

Slide 36 text

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

Slide 37

Slide 37 text

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

Slide 38

Slide 38 text

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

Slide 39

Slide 39 text

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 ,

Slide 40

Slide 40 text

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)

Slide 41

Slide 41 text

Slide 42

Slide 42 text

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)

Slide 43

Slide 43 text

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:

Slide 44

Slide 44 text

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

Slide 45

Slide 45 text

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)

Slide 46

Slide 46 text

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 .

Slide 47

Slide 47 text

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:

Slide 48

Slide 48 text

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

Slide 49

Slide 49 text

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

Slide 50

Slide 50 text

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 < +

Slide 51

Slide 51 text

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 < +

Slide 52

Slide 52 text

Numerical Examples on Meshes 19

Slide 53

Slide 53 text

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).

Slide 54

Slide 54 text

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).

Slide 55

Slide 55 text

Overview • Metrics and Riemannian Surfaces. • Geodesic Computation - Iterative Scheme •Geodesic Computation - Fast Marching • Shape Recognition with Geodesic Statistics • Geodesic Meshing 21

Slide 56

Slide 56 text

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.

Slide 57

Slide 57 text

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

Slide 58

Slide 58 text

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

Slide 59

Slide 59 text

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

Slide 60

Slide 60 text

Fast Marching on an Image 24

Slide 61

Slide 61 text

Fast Marching on Shapes and Surfaces 25

Slide 62

Slide 62 text

Volumetric Datasets 26

Slide 63

Slide 63 text

Propagation in 3D 27

Slide 64

Slide 64 text

Overview • Metrics and Riemannian Surfaces. • Geodesic Computation - Iterative Scheme • Geodesic Computation - Fast Marching •Shape Recognition with Geodesic Statistics • Geodesic Meshing 28

Slide 65

Slide 65 text

Bending Invariant Recognition 29 [Zoopraxiscope, 1876] Shape articulations:

Slide 66

Slide 66 text

Bending Invariant Recognition 29 [Zoopraxiscope, 1876] Shape articulations: ˜ x1 ˜ x2 M Surface bendings: [Elad, Kimmel, 2003]. [Bronstein et al., 2005].

Slide 67

Slide 67 text

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

Slide 68

Slide 68 text

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