Gabriel Peyré
January 01, 2011
730

Mesh Processing Course: Geodesic Sampling

January 01, 2011

Transcript

1. Geodesic Sampling and Meshing
Gabriel Peyré

2. Overview
•Riemannian Metrics.
•Riemannian Voronoi and Delaunay.
•Farthest Point Sampling.
•Anisotropic Triangulations.
•Anisotropic Delaunay Reﬁnement.
2

3. Riemannian Manifold
3
Riemaniann manifold: abstract parametric space M Rn.
Metric: M equiped x ⇥ M ⇤ H(x) ⇥ Rn n positive deﬁnite.
Length of a curve (t) M: L( ) def.
=
1
0

(t)TH( (t)) (t)dt.
W(x)
Examples:
Euclidean space: M = Rn and H(x) = Id
n
.
2D shape: M R2 and H(x) = Id
2
.
Parametric surface: H(x) = Ix
ﬁrst fundamental form.
Isotropic metric: H(x) = W(x)Id
n
, W(x) > 0 weight function.
Image processing: I : [0, 1]2 ⇥ R, W(x) = ( + ||⇤xI||) 1.
DTI imaging: M = [0, 1]3, H(x) =di usion tensor.

4. Geodesic Distances
4
Geodesic distance metric over M Rn
⇥ (x, y) M2, dM
(x, y) def.
= min
T >0, PT (x,y)
L( )
where PT
(x, y) def.
= { \ (0) = x and (T) = y} .
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
2 ECCV-08 submission ID 1057
Euclidean Shape Isotropic Anisotropic Surface

5. Overview
•Riemannian Metrics.
•Riemannian Voronoi and Delaunay.
•Farthest Point Sampling.
•Anisotropic Triangulations.
•Anisotropic Delaunay Reﬁnement.
5

6. Euclidean Delaunay Triangulation
6
Delaunay graph DS
: (xi, xj
) ⇤ DS
⇥ Ci
⌃ Cj
⌅= ⇧.

7. Euclidean Voronoi Segmentation
7

8. Euclidean Delaunay Triangulation
8
Delaunay graph DS
: (xi, xj
) ⇤ DS
⇥ Ci
⌃ Cj
⌅= ⇧.

9. Geodesic Voronoi and Delaunay
9
Voronoi segmentation:
Outer cell: C0
= Closure( c).
Delaunay graph DS
: (xi, xj
) ⇤ DS
⇥ Ci
⌃ Cj
⌅= ⇧.

10. Voronoi Diagrams on Surfaces
10

11. Double and Triple Points
11
S large enough = DS
is a planar triangulation.
Tripple point: si,j,k
Ci
⇥ Cj
⇥ Ck.
Double point: wi,j
= argmin
x Ci
⇥Cj
d(x, xi
).
Geometric realization of DS
:
xj
xk
xi
si,j,k
wi,j
Union of two geodesics starting from wi,j
.
Special case: boundary curves if xi
and xj
.
Fast computation in O(n log(n))
with Fast Marching.

12. Anisotropic Voronoi Segmentation
12
ECCV-08 submission ID1057 7
f
= .1 = .2 = .5 = 1
Fig. 4. Examples of anisotropic distances (top row) and Voronoi diagrams (bottom
row) with an decreasing anisotropy . The metric tensor is computed using the structure
tensor, equation (8).
= .5 = 0
= .95 = .7
Distances Voronoi

13. Overview
•Riemannian Metrics.
•Riemannian Voronoi and Delaunay.
•Farthest Point Sampling.
•Anisotropic Triangulations.
•Anisotropic Delaunay Reﬁnement.
13

14. Coverings and Packings
14

15. Geodesic Delaunay and Voronoi
15
xi
xi,j
xk
x
xi,j,k
xi,j,
xj

16. Farthest Point Sampling
16
Surface sampling {x1, . . . , xn
} M.
Parametric domain sampling: : [0, 1]2 ⇥ M, xi,j
= (i/

n, j/

n).
-covering:
i
B (xi
) = M, where B (x) def.
= {y \ dM
(x, y) }.
Theorem: The sampling {x1, . . . , xn
} is an -covering that is -separated for
= max
i=1,...,n
min
j=1,...,n
dM
(xi, xj
).
1. Initialization: x1
random, d(x) dM
(x0, x), set i = 1.
2. Select point: xi+1
= argmax
x
d(x), = d(xi+1
).
3. Local update of the distance: d(x) min(d(x), dM
(xi+1, x)).
4. Stop: If i < n or > 0
, set i i + 1 and go back to 2.
Farthest point
-separated: min(dM
(xi, xj
)) .

17. Farthest Point Sampling
17
# samples
Metric
W(x)
W(x)
W(x) small = front moves slowly,
= denser sampling.
W(x)

18. Farthest Point Triangulation
18
# samples

19. Meshing a Planar Shape
19
Uniform

20. Uniform Remeshing
20

21. Remeshing of the David
21
Original:
100k vert.
Remeshed:
10k vert.

22
# samples

23. Density Given by a Texture
23
Parameterized surface: : D [0, 1]2 ⇥ S R3.
Texture: I = [0, 1]2 R.
|| 1(p)
I||
⇥ p S, W(p) = ( + ||⇤ 1(p)
I||) 1
W = 1 large small
I( 1(p))

24. Curvature of a Height ﬁeld
24
Height ﬁeld: x ⇥ [0, 1]2 (x, f(x)) ⇥ R3.
f(x + h) = f(x) + ⇤xf, h⇥ + Hx
(f)h, h⇥ + O(||h||2)
f(x + h) (f(x) + ⇥⌅xf, h⇤) = µ1
⇥h, e1
⇤2 + µ2
⇥h, e2
⇤2 + O(h2).
Hx
= µ1e1e1
T + µ2e2e2
T with |µ1
| > |µ2
|.
Gaussian curvature: µ1
(x)µ2
(x) = det(Hx
(f)).
Mean curvature: (µ1
(x) + µ2
(x))/2 = tr(Hx
(f))/2.
Total curvature: |µ1
(x)| + |µ2
(x)|.
Intrinsic
Extrinsic
e1
(x)
e2
(x)
x

25. Curvature of a Surface
25
Numerical estimation:
• polynomial ﬁt.
• local covariance analysis.
Hx
= ⇥2 (x)
⇥xi⇥xj
, n(x)⇥

i,j=1,2
Second fundamental form:
Normal: n = ⇥
⇥x1

⇥x2
/ ⇥
⇥x1

⇥x2
.
Parametric surface: x ⇥ R2 ⇤ (x) ⇥ R3.
|µ1
(x)| + |µ2
(x)|.

26. Edge-aware Remeshing
26
Original mesh W = 1 Small
W(x) = ( + |µ1
(x)| + |µ2
(x)|) 1
Curvature-driven metric:

27. Overview
•Riemannian Metrics.
•Riemannian Voronoi and Delaunay.
•Farthest Point Sampling.
•Anisotropic Triangulations.
•Anisotropic Delaunay Reﬁnement.
27

28. Anisotropy is Important
28
Better respect of features.
Better approximation of functions.
8 Jonathan Richard Shewchuk
Figure 2: A visual illustration of how large angles, but not small angles, can cause the error ∇f − ∇g to
explode. In each triangulation, 200 triangles are used to render a paraboloid.
35
50
65
40 20
20
40
40
40
20
Figure 3: As the large angle of the triangle approaches 180◦, or the sliver tetrahedron becomes arbitrarily
ﬂat, the magnitude of the vertical component of ∇g becomes arbitrarily large.
sphere, then perturbing one of the vertices just off the equator so that the sliver has some (but not much)
volume.
Because of this sensitivity, mesh generators usually choose the shapes of elements to control ∇f −
∇g ∞, and not f − g ∞, which can be reduced simply by using smaller elements. Section 6.1 presents
quality measures that judge the shape of elements based on their ﬁtness for interpolation.
Table 2 gives two upper bounds on ∇f − ∇g ∞ over a triangle. The ﬁrst upper bound is almost tight,
to within a factor of two. The “weaker but simpler upper bound” of 3ctrcirc is not as good an indicator as
the stronger upper bound, but it has the advantages of being smooth almost everywhere (and therefore more
Isotropic Anisotropic, Anisotropic,
f(x + h) = f(x) + ⇤xf, h⇥ + Hx
(f)h, h⇥ + O(||h||2)
f(x + h) (f(x) + ⇥⌅xf, h⇤) = µ1
⇥h, e1
⇤2 + µ2
⇥h, e2
⇤2 + O(h2).
= Local optimal shape:
Hx
= µ1e1e1
T + µ2e2e2
T with |µ1
| > |µ2
|.
width/length = |µ2
|/|µ1
|
x
e1

µ2

µ1

29. Approximation of Images with Edges
29
For a geometrically regular image f, that is C2 outside C2 edges, one can however built an
adapted triangulation that enhances signiﬁcantly the approximation performance of a wavelet
basis. The following theorem sketch the construction of such a triangulation.
f ∗ h
f
Figure 9.2: Approximation with ﬁnite elements on a triangulation for a function without and with
Theorem 2. If f is C2 outside a set of C2 contours, then it exists a constant C such that for
any M, it exists an adapted triangulation (V, T) over which the piecewise-linear interpolation fM
satisﬁes
||f fM
||2 C M 2. (9.2)
Proof. We give here only a sketch of the proof. [todo : show an image of the tubes] Near an edges,
the lengths of the triangles should be of order M 1 and their widths should be of order M 2. We
deﬁne a thin band of width M 2 around all the edge curves, see ﬁgure ??. Since the edge curve are
C2, this band can be sub-divided in elongated tubes of length proportional to M 1, each of witch is
decomposed in two elongated triangles. The number of such triangle is 2LM where L is the total length
of the edges in the image. A special care should be taken near edge corner or edge crossings but we
ignore these technicalities here. Over these tubes, the function is bounded and thus the approximation
error ||f e
f||2
L2( )
in the tubes is of the order of area( )||f||2

= LM 2||f||2

. In the complementary
of the tubes c, one packs M large equilateral triangles of edge lengths approximatively ⌃ ⇥ M 1/2.
Since f is C2 outside the edge curve, the approximation error of a linear interpolation e
f over such
The piecewise linear interpolation ⇥
fM
of f on (V, T) is such that ⇥
fM
(xi
) = f(xi
) and ⇥
fM
is
linear on each triangular face tj
. This is a two-dimensional extension of the spline approximations
studied in Section ??. This piecewise linear approximation depends on the position of the vertices
V and the connectivity T of the triangulation. In order to e⇥ciently approximate a given function
f with M triangles, one needs to ﬁnd the optimal shape of the triangles.
For a geometrically regular image f, that is C2 outside C2 edges, one can however built an
adapted triangulation that enhances signiﬁcantly the approximation performance of a wavelet
basis. The following theorem sketch the construction of such a triangulation.
f ∗ h
f
Figure 9.2: Approximation with ﬁnite elements on a triangulation for a function without and with
Theorem 2. If f is C2 outside a set of C2 contours, then it exists a constant C such that for
any M, it exists an adapted triangulation (V, T) over which the piecewise-linear interpolation ⇥
fM
satisﬁes
||f ⇥
fM
||2 C M 2. (9.2)
Proof. We give here only a sketch of the proof. [todo : show an image of the tubes] Near an edges,
9.4. Geometric Image Approximation
As already seen in the proof of theorem ??, one can design triangles of width and legnth M
the tube so that the error in this area is also of the order of M 2, which concludes the p
The proof of this theorem shows that when an image is smoothed by an unknow
width ⇤, the triangulation should depend on ⇤ in order to get a fast decay of the app
error. To reach an error decay of O(M 2), in the neighborhood of a contour smoothed
of width ⇤, the triangle should have a length of order ⇤1/4M 1/2 and a width of order ⇤
as shown on ﬁgure ??. The scale ⇤ is most of the time unknown and one thus needs an
algorithm to devise the size of the triangles.
s1/ 4M1/ 2
triangles
s1/ 4M−1/ 2
s3/ 4M−1/ 2
Figure 9.4: Aspect ratio of triangles for the approximation of a blurred contou
This shows that locally, an optimal triangle should be aligned with the direction µ
function is the most regular. The local behavior of a smooth C2 image is well descr
second order by a quadratic approximation
f(x + h) = f(x) + ⌅⌥xf, h⇧ + ⌅Hx
(f)h, h⇧ + O(||h||2)
2⇥2
||f ⇥
fM
||2 C M 2. (9.2)
. We give here only a sketch of the proof. [todo : show an image of the tubes] Near an edges,
ngths of the triangles should be of order M 1 and their widths should be of order M 2. We
e a thin band of width M 2 around all the edge curves, see ﬁgure ??. Since the edge curve are
his band can be sub-divided in elongated tubes of length proportional to M 1, each of witch is
mposed in two elongated triangles. The number of such triangle is 2LM where L is the total length
e edges in the image. A special care should be taken near edge corner or edge crossings but we
e these technicalities here. Over these tubes, the function is bounded and thus the approximation
||f e
f||2
L2( )
in the tubes is of the order of area( )||f||2

= LM 2||f||2

. In the complementary
e tubes c, one packs M large equilateral triangles of edge lengths approximatively ⇤ M 1/2.
f is C2 outside the edge curve, the approximation error of a linear interpolation e
f over such
gles is ||f e
f||⇤( c
||f||C
2 and thus the error satisﬁes
||f e
f||2 = O(LM 2||f||2

+ ||f||2
C
M 2).
M−1
M−2
M−1
M−α
Figure 9.3: Finite elements for the approximation around a singularity curve.
theorem shows that an adapted triangulation should balance the approximation error
d outside the singularity curves. The triangles located far away from these singularities
e large and isotropic, while triangles that cover the edges should be stretched along the
y curves. This construction can be generalized by replacing triangles by higher order
c primitives whose boundaries are polynomial curves of degree , as shown on the right
gure ??. The adapted approximation using polynomials deﬁned on M such higher order
Sharp edges
Blurred edges
Sharp edges Blurred edges

30. Metric Design for Anisotropic Meshing
30
Bossen−Heckbert [1996] George−Borouchaki [1998]
Generating Anisotropic Meshes
Heuristic Algorithms for
x
2
(x)
1
(x)
e1
(x)
Local orientation: e1
(x).
Local shape:
1
(x)/ 2
(x).
Idea: design the metric H(x).

31. Image Approximation
31
Garland & Heckbert, Approximation of Terrains 19
Figure 16: Mandrill original, a 200×200 raster image. Figure 17: Mandrill approximated with Gouraud
shaded triangles created by subsampling on a uniform
20×20 grid (400 vertices).
Figure 16: Mandrill original, a 200×200 raster image. Figure 17: Mandrill approximated with Gouraud
shaded triangles created by subsampling on a uniform
20×20 grid (400 vertices).
Figure 18: Mandrill approximated with Gouraud Figure 19: Mesh for the image to the left.
Original image Isotropic mesh (400 vert)
Anisotropic mesh (400 vert)

32. Image Compression
32
16 L. Demaret, N. Dyn, M.S. Floater, A. Iske
(a) (b)
(c) (d)
Fig. 8. Reflex. (a) Original image of size 128 × 128. Compression at 0.251 bpp and
reconstruction by (b) SPIHT with PSNR 30.42 db, (d) AT∗
2
with PSNR 41.73 db.
16 L. Demaret, N. Dyn, M.S. Floater, A. Iske
(a) (b)
(c) (d)
Fig. 8. Reflex. (a) Original image of size 128 × 128. Compression at 0.251 bpp and
Original image Wavelets @.25bpp
Triangulation @.25bpp

33. Overview
•Riemannian Metrics.
•Riemannian Voronoi and Delaunay.
•Farthest Point Sampling.
•Anisotropic Triangulations.
•Anisotropic Delaunay Reﬁnement.
33

34. Planar Domains
34
This section, shows how several tools from computational geometry extend
105
to the setting of a Riemannian metric.
106
Starting from a set of points S = {xi
}m
i=1
, one can deﬁne graphs and trian-
107
gulations that reﬂect the geometry of the Riemannian manifold. These points
108
and the corresponding graphs are the basic building blocks of the algorithms for
109
perceptual grouping and planar domain meshing.
110
In the following, the boundary ⇧ of is
assumed to be a set of closed smooth curves.
At least one point of S is located on each
curve, and these boundary points segment
⇧ as a set of sub-curves
⇧ =

i,j
⇥i,j
with ⇥i,j
⇤ P(xi, xj
)
⇥i,j
↵ S = {xi, xj
}.
(one can have xi
= xj
if there is only one
point on a curve).

Ωc
x1
x2
x3
x4
111
2.1 Delaunay and Voronoi Graphs
112
The segmentation of the domain in Riemannian Voronoi cells is
113
= C0

xi
S
Ci
where Ci
= {x ⇤ \ ⇧ j ⌅= i, d(xi, x) d(xj, x)} . (7)
The outer Voronoi cell is deﬁned as C0
= Closure( c).
114
The Delaunay graph DS
of S is a graph where two points are connected if
their respective Voronoi cells are adjacent
(xi, xj
) ⇤ DS
⇥ Ci
↵ Cj
⌅= ⌃.
To each Delaunay edge (xi, xj
) ⇤ DS
corresponds a double point wi,j
, which is
the closest point to xi
and xj
on the common Voronoi cell boundary
wi,j
= argmin
x Ci
⌅Cj
d(x, xi
).
ECCV-08 submission ID1057 11
e to the boundary. A boundary sub-curve
i,j
is 245
S if it exists a triple point wi,k,0 i,j
, see ﬁgure 246
nnot be part of the Delaunay triangulation, and 247
orithm by inserting a mid point. Similarly, triple 248
croach any boundary sub-curve (the sub-curve is 249
250
he Delaunay graph DS
of S is not necessarily a 251
ng S is not dense enough, see [7]. This is because 252
connected to only one other point of S in DS
. 253
dd points on the Voronoi cell boundary of such a 254
255
x2 x1 x2
x3
aches the boundary curve ⇥1,2. Right: the vertex x3
e (x1, x2) is a Delaunay edge.
ECCV
if a point xk
is located too close to the boundary.
245
said to be encroached by xk
S if it exists a triple
246
7. Such an encroached edge cannot be part of the
247
is automatically split by the algorithm by inserting
248
points are not added if they encroach any boundar
249
250
Another di⇥culty is that the Delaunay graph
251
valid triangulation if the sampling S is not dense en
252
of some isolated point, that is connected to only
253
The algorithm automatically add points on the Vo
254
point.
255
x1 x2
x3
x1
Fig. 7. Left: the vertex x3 encroaches the boundary
does not encroach anymore because (x1, x2) is a Delau
2D Riemannian manifold: [0, 1]2 equipped with H(x) ⇥ R2 2.
Boundaries:
⇥ =

i,j
i,j
with i,j
P(xi, xj
)
i,j
⇥ S = {xi, xj
}.
Requirement: DS
.
encroaches
1,2

35. Face and Edge Splitting
35
Voronoi edge
Delaunay edge
Missing boundary x3 x4
Bad quality triangle (x1, x2, x3
)

36. Riemannian Delaunay Reﬁnement
36
It extends the isotropic farthest point seeding strategy of [?] w
242
metrics and domains with arbitrary boundaries.
243
Our anisotropic meshing algorithm proceeds by iteratively i
point si,j,k
⇥ S
to an already computed set of points S. In or
an anisotropic mesh with triangles of high quality with resp
metric, one inserts si,j,k
for a Delaunay triangle (xi, xj, xk
) w
⇥(si,j,k
) = d(si,j,k, xi
)
min(d(xi, xj
), d(xj, xk
), d(xk, xi
)),
which is a quantity computed for each triple point in parallel to
244
ing propagation. In the Euclidean domain, a triangle (xi, xj, xk
)
245
of ⇥(si,j,k
) is badly shaped since its smallest angle is close to 0 .
246
[?], this property extends to an anisotropic metric H(x) if angl
247
using the inner product deﬁned by H(x).
248
shortest edge
Fig. 7. Left: the vertex x3 encroaches the boundary curve ⇥1,2. Right: the vertex x3
does not encroach anymore because (x1, x2) is a Delaunay edge.
Table 2 details this algorithm. A bound ⇥⇥ on ⇥ enforces the reﬁnement to
reach some quality criterion, while a bound U⇥ enforces a uniform reﬁnement to
match some desired triangle density.
1. Initialization: set S with at least one point on each curve of , compute US with
a Fast Marching.
2. Boundary enforcement: while it exists ⇥i,j ⌃ encroached by some xk
⇤ S,
subdivide: S ⇥ S ⌃ argmax
x⇤⇤i,j
US (x). Update US with a local Fast Marching.
3. Triangulation enforcement: while it exists (xi, xj) ⇤ DS with xi or xj isolated,
insert w⇥ = argmax
w⇤Ci⇧Cj
d(xi, w).
4. Select point: s⇧ ⇥ argmin
s⇤ S ⌃⇥
⇤(s).
– If in S ⌃ {s⇧}, s⇧ encroaches some ⇥i,j ⌃ , subdivide: S ⇥ S ⌃ argmax
x⇤⇤i,j
US (x).
– Otherwise, add it: S ⇥ S ⌃ s⇧.
Update US with a local Fast Marching.
5. Stop: while ⇤(s⇧) > ⇤⇧ or US (s⇧) > U⇧, go back to 2.
x1
x3
s1,2,3
s1,2,3
x3
x2
x2
x1

37. Isotropic vs. Anisotropic Meshing
37
Isotropic
Anisotropic

38. Examples of Anisotropic Meshing
38
M,
Given polygonal domain and metric tensor field
generate anisotropic mesh.
nisotropic mesh.
If metric tensor
derivatives, no triangle has angle < 20° as
measured by any point in the triangle.
Main Result
M is smooth with bounded
etric dual of an anisotropic Voronoi di-
angulation. We describe conditions in
guaranteed to be entirely visible from
ction 5. For the special case of two di-
ons also guarantee that the planar dual
with no inverted triangles.
ible an algorithm that generates high-
isotropic meshes by reﬁ ning an aniso-
nforce the conditions that guarantee that
nd to remove any poor-quality elements

39. Bonus
• Centroidal Tesselations
39

40. Centroidal Tesselation
40
A ⇥ Rn, Euclidean center of gravity: g(A) = argmin
x
A
||x y||2dy.
A M, geodesic center of gravity: gM
(A) = argmin
x
A
dM
(x, y)2dy
Sampling {xi
}m
i=1
, Voronoi tessellation Voronoi
M
({xi
}i
) = {Vi
}m
i=1
.
Centroidal tessellation i, xi
= gM
(Vi
), local minimizer of
E({xi
}, {Vi
}i
) =
m
i=1

Vi
dM
(xi, y)2dy.
1. Initialize: {xi
}i
random.
2. Update regions: ⇥ i, {Vi
}i
Voronoi
M
({xi
}i
).
3. Update centers: ⇥ i, xi gM
(Vi
).
4. Stop: while not converged, go back to 2.
Lloyd-Max
x
y

nx
(y)
EA
(x) =
A
dM
(x, y)2dy =⇤ ⌅x
EA
=
A
dM
(x, y)

nx
(y)dy

41. Centroidal Tesselations of the Plane
41
= Delaunay = Voronoi
uniform