Gabriel Peyré
January 01, 2011
820

# Mesh Processing Course: Geodesic Mesh Processing

January 01, 2011

## Transcript

1. Geodesic Data Processing
Gabriel Peyré
www.numerical-tours.com

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}
(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}
(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.
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

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

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.

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

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