Geodesic Sampling and Meshing Gabriel Peyré CEREMADE, Université Paris Dauphine http://www.ceremade.dauphine.fr/~peyre/

Overview •Riemannian Metrics. •Riemannian Voronoi and Delaunay. •Farthest Point Sampling. •Anisotropic Triangulations. •Anisotropic Delaunay Refinement. 2

Riemannian Manifold 3 Riemaniann manifold: abstract parametric space M Rn. Metric: M equiped x ⇥ M ⇤ H(x) ⇥ Rn n positive definite. 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 first 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.

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

Overview •Riemannian Metrics. •Riemannian Voronoi and Delaunay. •Farthest Point Sampling. •Anisotropic Triangulations. •Anisotropic Delaunay Refinement. 5

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

Euclidean Voronoi Segmentation 7

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

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

Voronoi Diagrams on Surfaces 10

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.

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

Overview •Riemannian Metrics. •Riemannian Voronoi and Delaunay. •Farthest Point Sampling. •Anisotropic Triangulations. •Anisotropic Delaunay Refinement. 13

Coverings and Packings 14

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

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

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

Farthest Point Triangulation 18 # samples

Meshing a Planar Shape 19 Uniform Adaptive

Uniform Remeshing 20

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

Adaptive Remeshing 22 # samples

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

Curvature of a Height field 24 Height field: 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

Curvature of a Surface 25 Numerical estimation: • polynomial fit. • 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)|.

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

Overview •Riemannian Metrics. •Riemannian Voronoi and Delaunay. •Farthest Point Sampling. •Anisotropic Triangulations. •Anisotropic Delaunay Refinement. 27

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 flat, 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 fitness for interpolation. Table 2 gives two upper bounds on ∇f − ∇g ∞ over a triangle. The first 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, bad shape good shape 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

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 significantly 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 finite elements on a triangulation for a function without and with additional blurring. 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 satisfies ||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 define a thin band of width M 2 around all the edge curves, see figure ??. 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 find 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 significantly 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 finite elements on a triangulation for a function without and with additional blurring. 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 satisfies ||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 figure ??. 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 figure ??. 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 satisfies ||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 defined on M such higher order Sharp edges Blurred edges Sharp edges Blurred edges

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

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)

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

Overview •Riemannian Metrics. •Riemannian Voronoi and Delaunay. •Farthest Point Sampling. •Anisotropic Triangulations. •Anisotropic Delaunay Refinement. 33

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 define graphs and trian- 107 gulations that reflect 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 defined 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 figure 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 subdivided instead). 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 . Good situation Bad: x3 encroaches 1,2

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

Riemannian Delaunay Refinement 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 circumradius to shortest edge ratio ⇥(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 defined by H(x). 248 = circumradius 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 refinement to reach some quality criterion, while a bound U⇥ enforces a uniform refinement 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

Isotropic vs. Anisotropic Meshing 37 Isotropic Anisotropic

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 refi ning an aniso- nforce the conditions that guarantee that nd to remove any poor-quality elements

Bonus • Centroidal Tesselations 39

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

Centroidal Tesselations of the Plane 41 = Delaunay = Voronoi uniform adaptive #iterations

Centroidal Tesselations of Surfaces 42 #iterations