Gabriel Peyré
email: peyre@cmapx.polytechnique.frcohen@ceremade.dauphine.fr

Laurent Cohen

Surface Segmentation Using Geodesic Centroidal Tesselation

Keywords:

In this paper, we solve the problem of mesh partition using intrinsic computations on the 3D surface. The key concept is the notion of centroidal tesselation that is widely used in an eucidan settings. Using the Fast Marching algorithm, we are able to recast this powerful tool in the language of mesh processing. This method naturally fits into a framework for 3D geometry modelling and processing that uses only fast geodesic computations. With the use of classical geodesic-based building blocks, we are able to take into account any available information or requirement such as a 2D texture or the curvature of the surface.

Introduction

The applications of 3D geometry processing abound nowadays. They range from finite element computation to computer graphics, including solving all kinds of surface reconstruction problems. The most common representation of 3D objects is the triangle mesh, and the need for fast algorithms to handle this kind of geometry is obvious. Classical 3D triangulated manifold processing methods have several well identified shortcomings: mainly, their high complexity when dealing with large meshes, and their numerical instabilities.

To overcome these difficulties, we propose a geometry processing pipeline that relies on intrinsic information of the surface and not on its underlying triangulation. Borrowing from well established ideas in different fields (including image processing and perceptual learning) we are able to process very large meshes efficiently.

Overview

In section 2 we introduce some concepts we use in our geodesic computations. This includes basic facts about the Fast Marching algorithm, the computations of Voronoi diagrams on surfaces, and a recently proposed greedy algorithm for manifold sampling.

Then, in section 3, we will introduce the notion of geodesic centroidal tessellation in order to compute a segmentation of the manifold. There, a fast algorithm based on a geodesic gradient descent is derived, and we will explain how to take into account special features of the mesh, or user defined constraints.

In the conclusion, we will apply this partioning algorithm to texture each patch of the mesh. We will then give a complete study of the timings of each part of our algorithm, including a comparison with classical methods.

Related Work

Geodesic Computations. Distances computation on manifolds is a complex topic, and a lot of algorithms have been proposed such as Chen and Han shortest path method [START_REF] Chen | Shortest Path on a Polyhedron[END_REF] which is of quadratic complexity. Kimmel and Sethian's Fast Marching algorithm [START_REF] Kimmel | Computing Geodesic Paths on Manifolds[END_REF] allows finding numerically the geodesic distance from a given point on the manifold, in O(n log(n)) in the number of vertices. They deduce minimal geodesics between two given points. Some direct applications of geodesic computations on manifolds have been proposed, such as in [START_REF] Kimmel | Fast Voronoi Diagrams on Triangulated Surfaces[END_REF], which applies the Fast Marching algorithm to obtain Voronoi diagram and offset curves on a manifold. Mesh Segmentation. The segmentation of a complex surface is the first step for constructing a parameterization of a whole mesh. In [START_REF] Levy | Least Squares Conformal Maps for Automatic Texture Atlas Generation[END_REF] such a segmentation is performed to build an atlas for mesh texturing. In [START_REF] Zhang | Feature-based surface parameterization and texture mapping[END_REF] the segmentation borrows ideas from Morse theory. The work of [START_REF] Katz | Hierarchical mesh decomposition using fuzzy clustering and cuts[END_REF] decomposes the object into meaningful components using fuzzy clustering. The importance of texture seams reduction is studied in [START_REF] Sheffer | Seamster: Inconspicuous Low-Distortion Texture Seam Layout[END_REF]. The geometry images paradigm [START_REF] Gu | Geometry Images[END_REF] uses a single chart to parametrize the mesh on a 2D image. In [START_REF] Sander | Multi-chart Geometry Images[END_REF], the authors also use a segmentation procedure to build geometry images for complex meshes. This segmentation is computed using a form of combinatorial Lloyd iteration. In section 3 we will describe a more powerful approach that uses real geodesic computations. The main issue of these techniques is the discontinuities of the parameterization along patch boundaries. These discontinuities can be nearly removed as explained in [START_REF] Khodakovsky | Globally Smooth Parameterizations with Low Distortion[END_REF] at the cost of additional complexity.

Geodesic-Based Building Blocks

Fast Marching Algorithm

The classical Fast Marching algorithm is presented in [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF], and a similar algorithm was also proposed in [START_REF] Tsitsiklis | Efficient Algorithms for Globally Optimal Trajectories[END_REF]. This algorithm is used intensively in computer vision, for instance it has been applied to solve global minimization problems for deformable models [START_REF] Cohen | Global Minimum for Active Contour models: A Minimal Path Approach[END_REF].

This algorithm is formulated as follows. Suppose we are given a metric P(s)ds on some manifold S such that P > 0. If we have two points x 0 , x 1 ∈ S , the weighted geodesic distance between x 0 and x 1 is defined as

d(x 0 , x 1 ) def. = min γ 1 0 ||γ ′ (t)||P(γ(t))dt , ( 1
)
where γ is a piecewise regular curve with γ(0) = x 0 and γ(1) = x 1 . When P = 1, the integral in (1) corresponds to the length of the curve γ and d is the classical geodesic distance. To compute the distance function U(x)

def.

= d(x 0 , x) with an accurate and fast algorithm, this minimization can -→ n x , where -→ n x is the exterior unit vector normal to the curve at x, and the function U satisfies the nonlinear Eikonal equation:

||∇U(x)|| = P(x). ( 2
)
The function F = 1/P > 0 can be interpreted as the propagation speed of the front C t . The Fast Marching algorithm on an orthogonal grid makes use of an upwind finite difference scheme to compute the value u of U at a given point x i, j of a grid:

max(u -U(x i-1, j ), u -U(x i+1, j ), 0) 2 + max(u -U(x i, j-1 ), u -U(x i, j+1 ), 0) 2 = h 2 P(x i, j ) 2 .
This is a second order equation that is solved as detailed for example in [START_REF] Cohen | Multiple Contour Finding and Perceptual Grouping Using Minimal Paths[END_REF]. An optimal ordering of the grid points is chosen so that the whole computation only takes O(N log(N)), where N is the number of points.

In [START_REF] Kimmel | Computing Geodesic Paths on Manifolds[END_REF], a generalization to an arbitrary triangulation is proposed. This allows performing front propagations on a triangulated manifold, and computing geodesic distances with a fast and accurate algorithm. The only issue arises when the triangulation contains obtuse angles. The numerical scheme presented above is not monotone anymore, which can lead to numerical instabilities. To solve this problem, we follow [START_REF] Kimmel | Computing Geodesic Paths on Manifolds[END_REF] who proposes to "unfold" the triangles in a zone where we are sure that the update step will work [START_REF] Kimmel | Computing Geodesic Paths on Manifolds[END_REF]. To get more accurate geodesic distance on meshes of bad quality, one can use higher order approximations, e.g. [START_REF] Manay | Second-order Models for Computing Distance Transforms[END_REF], which can be extended to triangulations using a local unfolding of each 1-ring.

Figure 2 shows the propagation of a front and the calculation of a geodesic path computed using a gradient descent of the distance function.

We should emphasize the fact that the Fast Marching computations are consistent with the continuous setting, which is not the case for traditional graph-based methods. When the triangulation becomes denser, geodesics will converge to the exact geodesics of the underlying continuous manifold. For that reason, we will work with this continuous setting in mind, almost without any reference to the underlying triangulation that supports the computations.

Extraction of Voronoi Regions

It is possible to start several fronts from points {x 1 , . . . , x n } and make them evolve together, as shown on figure 3. The areas with various colours define the Voronoi diagram of the starting points, namely the tessellation into the regions, for i ∈ {1, . . . , n}

V i def. = {x ∈ S \ ∀ j = i, d(x, x j ) > d(x, x i )}.
To accurately compute the boundaries of the Voronoi regions, we allow an overlap of the front on one vertex. Suppose a front a arrives at a vertex v 1 with time arrival t a

λ )v 1 + λ v 2 where λ = d a 2 -d a 1 +d b 1 -d b 2 d b 1 -d a 1 .

A Greedy Algorithm for Sampling a Manifold

A new method for sampling a 3D mesh was recently proposed in [START_REF] Peyré | Geodesic Remeshing Using Front Propagation[END_REF] that follows a farthest point strategy based on the weighted distance obtained through Fast Marching on the initial triangulation. This is related to the method introduce in [START_REF] Cohen | Multiple Contour Finding and Perceptual Grouping Using Minimal Paths[END_REF]. A similar approach was proposed independently and simultaneously in [START_REF] Moenning | Fast Marching Farthest Point Sampling[END_REF]. It follows the farthest point strategy, introduced with success for image processing in [START_REF] Eldar | The Farthest Point Strategy for Progressive Image Sampling[END_REF] and related to the remeshing procedure of [START_REF] Chew | Guaranteed-Quality Mesh Generation for Curved Surfaces[END_REF]. This greedy solution for sampling has been used with success in other fields such as computer vision (component grouping, [START_REF] Cohen | Multiple Contour Finding and Perceptual Grouping Using Minimal Paths[END_REF]), halftoning (void-and-cluster, [START_REF] Ulichney | The Void-and-Cluster Method for Generating Dither Arrays[END_REF]) and remeshing (Delaunay refinement, [START_REF] Ruppert | A Delaunay Refinement Algorithm for Quality 2-Dimensional Mesh Generation[END_REF]).

This approach iteratively adds new vertices based on the geodesic distance on the surface. Figure 4 shows the first steps of our algorithm on a square. The result of the algorithm gives a set of vertices uniformly distributed on the surface according to the geodesic distance.

Once we have found enough points, we can link them together to form a geodesic Delaunay triangulation. This is done incrementally during the algorithm, and leads to a powerful remeshing method that can be either uniform (figure 5) or adaptive.

Adaptive Sampling Using Texture Information

To introduce some adaptivity in the sampling performed by this algorithm, we use a speed function F = 1/P (which is the right hand side of the Eikonal equation ( 2)) that is not constant across the surface.

When a mesh is obtained from range scanning, a picture I of the model can be mapped onto the 3D mesh. Using a function F of the form

F(x) = 1 1+µ|grad(I(x))| ,
where µ is a user-defined constant, one can refine regions with high variations in intensity. Figure 6 shows such a 3D model. On figure 1, one can see a 3D head remeshed with various µ ranging from µ = 0 (uniform) to µ = 20/ max(|grad(I(x))|) (highly adaptive). Local density based on curvature information can be used, see [START_REF] Peyré | Geodesic Remeshing Using Front Propagation[END_REF].

In the following section, we will make heavy use of adaptive sampling based on the curvature of the surface. Let us denote by τ(x) def.

= |λ 1 | + |λ 2 | the total curvature at a given point x of the surface, where λ i are the eigenvalues of the second fundamental form. We can introduce two speed functions

F 1 (x) def. = 1 + ετ(x) and F 2 (x) def. = 1 1 + µτ(x) , (3)
where ε and µ are two user-defined parameters. Figure 7 (a) shows that by using function F 1 , we avoid putting more vertices in regions of the surface with high curvature. The speed function F 1 can be interpreted as an "edge repulsive" function. On the other hand, function F 2 could be called "edge attractive" function, since it forces the sampling to put vertices in region with high curvature such as mesh corners and edges. Figure 7 (b) shows that this speed function leads to very good results for the remeshing of a surface with sharp features, which is obviously not the case for the "edge repulsive" speed function (figure 7 (a)).

Mesh Segmentation Using Centroidal Tessellation

In this section, we present an automatic mesh segmentation method which is the first part of our surface parameterization pipeline. This method is adapted to the tessellation of a complex manifold in elementary domains topologically equivalent to a disk.

In this method the mesh is cut into regions that best satisfy the following properties: (C1) Boundaries of the regions agree with sharp features of the surface. (C2) Regions are as compact as possible (the ratio area/perimeter should be large), enclosing equal areas. Condition (C1) is natural because the resulting parameterization will have discontinuities on the boundaries of the patches, and it has no consequences if these singularities coincide with singularities of the surface. Condition (C2) prevents degenerated patches and minimizes the length of the singularities which occur along the boundary. These two conditions are contradictory, and a trade-off can be found using an iterative approach.

To achieve this trade-off, we first make an initial choice of regions using our greedy sampling algorithm to satisfy (C1). Then we use these regions to seed a geodesic-based centroidal tesselation to match requirement (C2).

Initial Choice of Regions

The goal of the algorithm is to build a segmentation S = n i=1 V i of a triangulated manifold S . The V i will be the Voronoi regions associated with a given set of points {v 1 , . . . , v n }.

These points are chosen using the sampling algorithm of section 2.3. To avoid the clustering of base points v i near regions of high curvature, we choose the "edge repulsive" speed function F 1 of equation ( 3) for the sampling.

Once the v i are computed, we can build the Voronoi regions V i , using the method described in section 2.2. In order to force the boundaries of the regions V i to follow the discontinuities of the surface, we use the "edge attractive" speed function F 2 (v) defined in equation ( 3) for the Voronoi cells extraction. This will allow us to "freeze" the front in regions with high curvature. This way the resulting Voronoi regions will have boundaries aligned with sharp features of the surface, and condition (C1) will be satisfied.

We also put a topological constraint on the shape of the Voronoi regions so that each cell is a topological disk. In order to maintain this constraint, we keep track of the shape of the front during the Fast Marching procedure. When this front fails to enclose a topological disk (e.g. on a sphere with a single base, the front will collapse onto itself), we stop the process and add a new base point at the location of the failure.

Centroidal Voronoi Tessellation

The segmentation provided by the previous section satisfies condition (C1), but does not match the requirements of (C2). In order to make the size of regions more homogeneous, we use a clustering technique to refine the Voronoi segmentation. This method is well known in the case of a planar segmentation, and allows the building of centroidal Voronoi diagrams [START_REF] Du | Centroidal Voronoi Tessellations: Applications and Algorithms[END_REF].

In a formal fashion, the planar segmentation corresponds to minimizing the energy function

E(v i , V i ) def. = n ∑ i=1 V i ρ(y)d(y, v i ) 2 dy, (4)
where ρ represents a density function, and d is a distance on the plane. If one uses d(x, y) = ||x -y|| the Euclidean distance, with ρ = 1, and we can prove that:

• The V i regions must be the Voronoi cells of the v i points. • The v i points must be the centers of mass of the V i cells.

One can show that such a segmentation exists (although it need not be unique), and it is called a centroidal Voronoi tessellation. This segmentation can be computed using the Lloyd algorithm as follows. We are given an integer n and an initial configuration {(v i , V i )} n i=1 . The algorithm iteratively updates the V i and the v i using these steps:

(i) Compute the centers of mass v * i of V i . (ii) Replace v i with v * i . (iii) Compute the Voronoi regions V * i , associated with v i , for the metric d. (iv) Replace V i with V * i . If ∀i, V i = V * i ,
terminate the algorithm. Otherwise, go back to (i). Figure 8 shows some iterations of the algorithm, restricted to a square. If we omit the problem on the boundaries, the centroidal segmentation corresponds to the well-known bee-hive structure which gives both an optimal compactness and an equal repartition of area among the regions. It is important to note that we have started the algorithm with a random distribution of points. If we had used the output of the previous section (using our greedy sampling method), the number of iterations would have been reduced significantly (approximately 5 iterations instead of 100 for the same result).

The distribution of the points v i really corresponds to the notion of isotropic sampling. For example, suppose we want to construct a Delaunay triangulation with the size of the triangles controlled by a function H(x). Then, one can show that we may solve this problem asymptotically using a centroidal tessellation by choosing ρ(x) ∝ 1/H(x) 4 [START_REF] Du | Centroidal Voronoi Tessellations: Applications and Algorithms[END_REF].

Geodesic Lloyd Algorithm

In this paragraph, we explain how to generalize the centroidal Voronoi segmentation to surfaces. To that end we consider, on a triangulated manifold S , the weighted geodesic distance d defined by equation [START_REF] Alliez | Isotropic Surface Remeshing[END_REF].

The centroidal tessellation paradigm has already been recognized as a powerful tool for remeshing a triangulated manifold. In [START_REF] Alliez | Isotropic Surface Remeshing[END_REF], an isotropic remeshing is performed using a Lloyd algorithm in parameter space. We, however, propose to compute the Lloyd iterations directly on the mesh. This allows us to handle meshes with arbitrary genus and without the distortion introduced by the parameterization of the surface. In [START_REF] Surazhsky | Isotropic Remeshing of Surfaces: a Local Parameterization Approach[END_REF] such a parameterization is avoided by using a local Lloyd scheme. However, due to the scaling induced by the conformal parameterization, this is not equivalent to the computation of a geodesic centroidal tesselation. In [START_REF] Sander | Multi-chart Geometry Images[END_REF] a heuristic is used (back-propagation from the boundary) which does not solve the real center of mass equation.

To perform the Lloyd algorithm, we need to define the intrinsic center of mass v * i of V i , which is the minimizer of the energy function

E i (w) def. = x∈V i d(x, w) 2 ds.
where ds is the area element on the surface. On a Riemannian surface, such a minimizer exists and is unique under some mild assumptions on the local curvature, which are given in details in [START_REF] Jost | Riemannian Geometry and Geometric Analysis[END_REF]. It is important to note that this intrinsic center of mass is different from the classical Euclidean center of mass (which in general does not even lie on the surface). The computation of the intrinsic center of mass has been introduced in the computer vision community in [START_REF] Srivastava | Geometric Analysis of Constrained Curves for Image Understanding[END_REF] and is closely related to the generalised procrustes analysis, see [START_REF] Le | Mean Size-and-shape and Mean Shapes: a Geometric Point of View[END_REF]. We note, however, that this center of mass has never been explicitly used for mesh processing, for example in [START_REF] Du | A Constrained Centroidal Voronoi Tessellations for Surfaces[END_REF], the authors uses a projection of the center of mass, and not an intrinsic barycenter.

In order to compute v * i , we perform a gradient descent of the energy function using

--→ ∇E i (w) = 1 2 x∈V i d(w, x) ---→ n w (x) ds,
where ---→ n w (x) stands for the unit vector tangent at w to the geodesic path joining x and w. To estimate --→ ∇E i (w), we first compute the distance map from w to all points in V i .

Tracing back the geodesics, we are able to compute ---→ n w (x) for each vertex x of V i . At last, we use a first order numerical integration formula to approximate the integral. We compute the boundary of V i by interpolation as explained in section 2.2.

Figure 9 shows the results we obtain with a uniform speed function (so d is just the classical geodesic distance). We obtain a kind of geodesic bee-hive structure. The initialization of the seed point is done at random. With our greedy sampling algorithm for initialization, only one Lloyd iteration instead of three is enough to get the same centroidal tessellation.

But the most powerful feature of our method is that we can use a varying speed function to define the geodesic distance. In the following, we will use the "edge attractive" speed function F 2 , that we have already used in subsection 3.1 to construct the initial regions. This way we ensure that condition (C1) will still be satisfied while the Lloyd iterations improves the segmentation with respect to condition (C2).

To perform a case study of our method, the figure 10 takes the example of a cube with 10 base points:

• On the top row, we use a random initialization. After 2

Lloyd iterations, we get a tessellation that respects the faces of the cube, so it is compliant with condition (C1). This is due to the fact that we use the "edge attractive" speed function. After 3 more iterations, we get a perfectly centered segmentation, and condition (C2) is respected on each face. • With a random initialization, we do not get an optimal tessellation, since the distribution of base points is as follow: 1 face with 3 points, 2 faces with 2 points, and 3 faces with 1 point. On the bottom row, we use our greedy initialization. Thanks to the use of the "edge repulsive" speed function, the seeds are nearly optimally placed (no faces have 3 seeds anymore). We then perform the segmentation and Lloyd iterations using the "edge attractive" function. Only one iteration instead of three is enough to get a perfect segmentation.

Figure 11 shows the segmentation we obtain on more complex models. In the close-up we can see that the cell boundaries try to follow the edges of the mesh whenever it is possible.

Segmentation using texture information

Following the ideas of the section 2.4, we can use a texture function to modulate the speed functions F 1 and F 2 . The resulting segmentation will take in account both the texture intensity and the curvature information, according to the user will. Figure 12 shows the segmentation of a textured earth, and the resulting segmentation of the texture.

On such a simple model (a sphere), one could perform a segmentation directly on the 2D image (with periodic boundary conditions), with a special speed function that takes into acount both the gradient of the image and the length distortion due to the cylindrical coordinates. However, this cannot be easily extended to models with complex topology such as the one shown on figure 1. One could also use both the texture and the curvature information, as shown on figure 13.

Results and Discussion

Texturing of a Complex Model.

Once we have performed the segmentation of a 3D mesh, we can flatten each patch and compute texture mapping on this atlas. Although the study of this step is outside the scope of this paper, we note that most of these classical methods come from graph-drawing theory, [START_REF] Floater | Parameterization of Manifold Triangulations[END_REF] gives a survey of these techniques. The boundary-free formulation of [START_REF] Desbrun | Intrinsic Parameterizations of Surface Meshes[END_REF] and [START_REF] Levy | Least Squares Conformal Maps for Automatic Texture Atlas Generation[END_REF] (which are mathematically equivalent) is very interesting for texture mapping. We rather choose the boundary-free formulation of [START_REF] Peyré | Geodesic Computations for Fast and Accurate Surface Flattening[END_REF] that is not conformal but seems to better preserve area across the parameterization. This flattening scheme also uses only geodesic computations and naturally fits into an intrinsic framework for mesh processing together with our segmentation procedure.

On figure 14 one can see the whole pipeline in action. This includes first a centroidal tessellation of the mesh, then the extraction and flattening of each cell, and lastly the texturing of the model.

Figure 15 shows a bunny and a shark that have undergone the same texturing process. On the bunny, left image, all the base points for all patches have been depicted (30 point per patch). On the shark, left image, we have shown in blue the movement of some seed points during the Lloyd iteration process (4 iterations). This clearly shows one of the main features of our segmentation method, which is able to drive its base points to the most relevant areas with respect to the tessellation quality.

To store the texture into memory, we pack all the patches into a single square image, as shown on figure 16. We use a simple strategy to minimize the lost space, as suggested in [START_REF] Levy | Least Squares Conformal Maps for Automatic Texture Atlas Generation[END_REF]. It is important to note that the textures we use are mapped directly on the parameter space, without taking care of the junctions between boundaries. Although it is out of the scope of this article, a more complex texturing process could take the output of our pipeline and use it to perform real time painting directly on the mesh. Discussion. Our tests on various 3D models enlight the Our future works include a theoretical study of the convergence of the geodesic Lloyd algorithm. We also would like to analyze experimentally the quality of the parameterization. A good way of evaluating the efficiency of such a scheme is to use its output to perform mesh compression. The mesh atlas provided by our algorithm is an ideal pre-processing step for performing wavelet transform in parameter space, in a fashion similar to [START_REF] Sander | Multi-chart Geometry Images[END_REF].

Conclusion

We have described a new algorithm to perform the segmentation of a triangulated manifold. The main tool that allows to have a fast algorithm is the fast marching on a triangulated mesh, together with some improvements we added. This segmentation step is the first stage of most mesh processing pipelines. Our contribution includes a geodesic extension of the Lloyd algorithm that is able to construct a geodesic centroidal tessellation. This iterative algorithm takes into account curvature of the surface as well as texture information and is very well suited to building a set of base domains for mesh flattening.

Figure 2 .

 2 Figure 2. Front Propagation on the David, level sets of the distance function and geodesic path. be reformulated as follows. The level set curve C t def. = {x \ U(x) = t} propagates following the evolution equation ∂ C t ∂t (x) = 1 P(x)

Figure 3 .

 3 Figure 3. Progression of the fronts, Voronoi diagram, and resulting tessellation.

1 and

 1 another front b arrives at a vertex v 2 (connected to v 1 ) with time t b 2 . Allowing an overlap of the fronts, we record the time arrival t a 2 of a at v 2 , and t b 1 of b at v 1 . Then the two fronts collapse at (1 -

Figure 4 .

 4 Figure 4. A n overview of the greedy sampling algorithm. 300 points 5,000 points 20,000 points

Figure 5 .

 5 Figure 5. Geodesic remeshing with an increasing number of points.

Figure 6 .

 6 Figure 6. Remeshing of a 3D model using increasing weight for the speed function.

Figure 7 .

 7 Figure 7. Uniform versus curvature-based sampling and remeshing.

Figure 8 .

 8 Figure 8. Iterations of Lloyd algorithm on a square.

Figure 9 .Figure 10 .

 910 Figure 9. Lloyd iterations on various models.

Figure 11 .

 11 Figure 11. Segmentation of two complex models.

Figure 12 .Figure 13 .

 1213 Figure 12. Speed function, segmentation of the earth and resulting regions on the texture map (right).

Figure 15 .

 15 Figure 15. Texturing a bunny and a shark.

Figure 16 .

 16 Figure 16. Example of chart packing. strengths of our approach: • The boundary of our patches are smooth, almost straight in regions with low curvature variations. We can see that in parameter space, the base domain of a patch located in the middle of the David is a convex polyhedron with straight edges. • The Voronoi cells of the tessellation are very compact, with almost equal areas. In regions with low curvature variations, this leads to convex base domains, and in regions with sharp features, the boundary of the cells follows the discontinuities. • Using our greedy initialization, the convergence of the geodesic Lloyd algorithm is very fast (3 iterations).Our future works include a theoretical study of the convergence of the geodesic Lloyd algorithm. We also would like to analyze experimentally the quality of the parameterization. A good way of evaluating the efficiency of such a scheme is to use its output to perform mesh compression. The mesh atlas provided by our algorithm is an ideal pre-processing step for performing wavelet transform in parameter space, in a fashion similar to[START_REF] Sander | Multi-chart Geometry Images[END_REF].

Figure 14 .

 14 Figure 14. An overview of our pipeline. The mesh is first segmented using a weighted geodesic centroidal tessellation. Each resulting patch is then flattened using the Geodesic LLE procedure. At last, we can perform texture mapping on each base domain.


Gabriel Peyré
email: peyre@cmapx.polytechnique.frcohen@ceremade.dauphine.fr

Laurent Cohen

Geodesic Remeshing Using Front Propagation

In this paper, we present a method for remeshing triangulated manifolds by using geodesic path calculations and distance maps. Our work builds on the Fast Marching algorithm, which has been extended to arbitrary meshes by Sethian and Kimmel in [17]. First, a set of points that are evenly spaced across the surface is automatically found. A geodesic Delaunay triangulation of the set of points is then created, using a Voronoi diagram construction based on Fast Marching. At last, we use the distance information to find a simple parameterization of the manifold. Marching algorithm makes this method computationally inexpensive, and gives very good results. Examples are shown for synthetic and real surfaces.

Introduction

The applications of 3D manifold sampling are nowadays almost everywhere. It ranges from finite element computations to computer graphics, including all kinds of questions related to surface reconstruction techniques used in 3D image segmentation, like deformable models. The most common representation of 3D objects is the triangle mesh, and the need for the construction of a nice triangulation of a given surface is obvious. Thus it is often needed to enhance a given 3D structure (for example obtained from a 3D medical image or artist modeling) using a so-called remeshing algorithm.

The manifold parameterization problem is important as it is the basis for building "good" mesh representations. For example, obtaining a parameterization that minimizes cer-tain distortion measures is the first step used in texture mapping or semi-regular remeshing. Also, once a semi-regular representation of the mesh is built, it is easy to compute wavelet transforms [START_REF] Sweldens | The Lifting Scheme: A construction of second generation wavelets[END_REF][START_REF] Schröder | Spherical wavelets: efficiently representing functions on the sphere[END_REF] and perform data compression [START_REF] Guskov | Normal meshes[END_REF].

Overview

In this paper, we introduce in section 3 an isotropic remeshing algorithm that provides automatically a set of points that are spaced on the surface either uniformly or adaptively with a given density. Thiese points are defined iteratively, adding each point according to geodesic (weighted) distance map to the current set of points. A fast algorithm is proposed that updates by Fast Marching the distance map at each iteration. Basic facts about Fast Marching and geodesic computations are recalled in section 2. At last, in section 4, we propose a novel parameterization method based on geodesic information.

In order to parameterize a manifold, we first calculate a coarse mesh by defining a Voronoi diagram on the surface. The parameterization is then interpolated inside each geodesic triangle, by using appropriate distance data.

Related work

Geodesic calculations on manifolds is a large topic, and numerous algorithms have been proposed. Some of the most interesting exact methods are: -Chen and Han's shortest path algorithm [START_REF] Chen | Shortest path on a polyhedron[END_REF] is of quadratic order, and difficult to implement. -Polthier and Schmies's straightest geodesic algorithm [START_REF] Akeley | Straightest geodesics on polyhedral surfaces[END_REF] allows solving the problem with initial conditions (a starting point, and direction). This problem is well defined (uniqueness), but cannot be used for remeshing and building surface parameterizations. -Sethian and Kimmel's Fast Marching algorithm [START_REF] Sethian | Computing geodesic paths on manifolds[END_REF] allows solving the problem with boundary conditions (start and end points), which has not necessarily a unique solution. The algorithm is fast, of order O(n log(n)) in the number of vertices on the manifold. Some applications of geodesic computations on manifolds have been proposed, such as [START_REF] Kimmel | Fast Voronoi diagrams on triangulated surfaces[END_REF], which applies the Fast Marching algorithm to obtain Voronoi diagram and offset curves on a manifold. This idea is used in our algorithm to calculate Delaunay triangulations.

Remeshing methods roughly fall into two categories: -Isotropic remeshing: a surface density of points is defined, and the algorithm tries to position the new vertices to match this density. For example the algorithm of Terzopoulos and Vasilescu [START_REF] Terzopoulos | Adaptive meshes and shells: Irregular triangulation, discontinuities, and hierarchical subdivision[END_REF] uses dynamic models to perform the remeshing. -Anisotropic remeshing: the algorithm takes into account the principal directions of the surface to align locally the newly created triangles and/or rectangles. The algorithm proposed in [START_REF] Alliez | Anisotropic polygonal remeshing[END_REF] used lines of curvature to build a quaddominant mesh. Finite element methods make heavy use of such remeshing algorithms [START_REF] Kunert | Towards anisotropic mesh construction and error estimation in the finite element method[END_REF]. Our method relies on an isotropic distribution of points so that it can be cast into an Eikonal problem.

The last point presented in this paper is mesh parameterization. Some of the most effective methods available are: -Energy based approach: first introduced by Floater [START_REF] Floater | Parameterization of manifold triangulations[END_REF], these methods consist in solving a linear system in order to minimize a certain energy measure. Geodesic paths are often chosen as boundaries to partition the the mesh. -Geodesic surface flattening: in [START_REF] Zigelman | Texture mapping using surface flattening via multi-dimensional scaling[END_REF] Zigelman et al use geodesic information. They use multidimensional scaling to flatten the manifold while attempting to maintain the accuracy of distance information.

Geodesic Calculations

Fast Marching on an orthogonal grid

The classical Fast Marching algorithm is presented for example in [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF][START_REF] Cohen | Global minimum for active contour models: A minimal path approach[END_REF] for finding 2D paths or in [START_REF] Deschamps | Fast extraction of minimal paths in 3D images and applications to virtual endoscopy[END_REF][START_REF] Cohen | Grouping connected components using minimal path techniques. Application to reconstruction of vessels in 2D and 3D images[END_REF] for 3D extension and improvements. A similar algorithm was also proposed in [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF]. Fast Marching enables solving the nonlinear Eikonal equation:

∇U (x) = P (x), (1)
where U is the weighted distance function to a given set of points in the plane, and F = 1/P ≥ 0 is the speed of front propagation. In this section, we consider the orthogonal triangulation as represented on the left of figure 2. The Fast Marching algorithm makes use of an upwind finite differences scheme to compute the value u at a given point x i,j of the grid:

x x 1 x 2 x x 1 x 2
max(u -U (x i-1,j ), u -U (x i+1,j ), 0) 2 + max(u -U (x i,j-1 ), u -U (x i,j+1 ), 0) 2 = h 2 P (x i,j ) 2 .
This is a second order equation that is solved as detailed in [START_REF] Cohen | Multiple contour finding and perceptual grouping using minimal paths[END_REF]. An optimal ordering of the grid points is chosen so that the whole computation only takes O(N log(N )), where N is the number of points.

Study over a more general triangulation

Following [START_REF] Sethian | Computing geodesic paths on manifolds[END_REF], to generalize the above construction to an arbitrary triangulation, we consider the neighborhood of a point x represented by right side of figure 2. More specifically, we try to calculate a value for U (x) in the triangle x, x 1 , x 2 , if possible. Since several triangles around x may yield valid solutions, we consider only the smallest one. The problem that we face in approximating the gradient of U over the triangle x, x 1 , x 2 is that there is no "natural" coordinate system that we can use (unlike the case of an orthogonal grid). In [START_REF] Sethian | Computing geodesic paths on manifolds[END_REF] such a system is judiciously chosen, and trigonometric calculations are performed.

In a more general manner, [START_REF] Sethian | Fast methods for the Eikonal and related Hamilton-Jacobi equations on unstructured meshes[END_REF] proposes a gradient calculation method in arbitrary dimension. It consists in considering the matrix M of size 2 × 2 whose lines are constituted of the vectors --→ xx 1 and --→ xx 2 . We can then calculate the directional derivatives of U in x according to the vectors --→ xx 1 and --→ xx 2 with the first order formula:

ν := U (x) -U (x 1 ) U (x) -U (x 2 ) = U (x)a + b, where a = 1 1 and b = - U (x 1 ) U (x 2 )
. Now, the derivative in the direction of --→ xx k is by definition equal to ∇U, --→ xx k , hence the formula ∇U M -1 ν. By plugging this into the Eikonal equation (1), and denoting Q := (MM T ) -1 , we get:

(a T Qa)U (x) 2 -(2a T Qb)U (x) + b T Qb = P (x) 2 .
This is a second order equation, its discriminant being positive by the Cauchy-Schwartz formula. The upwind propagation condition we impose, as is the case in the orthogonal grid, is that -∇U (x) points towards the interior of the triangle. We can show that if the angle in x is acute, then the scheme is monotone.

Figure 3 shows the propagation of a front and the calculation of a geodesic path (see 2.3).

Figure 3. Front Propagation and geodesic path (on the right).

If, on the other hand, the triangulation contains obtuse angles, then the numerical scheme presented above is not monotone anymore, which can lead to numerical instabilities. To solve this problem, we follow Sethian who proposes to "unfold" the triangles in a zone where we are sure that the update step will work [START_REF] Sethian | Computing geodesic paths on manifolds[END_REF].

Geodesic Extraction

Once we have calculated the function U , which is the geodesic distance to a vertex v, we need to determine the geodesic path starting from a point v 0 and reaching the point v. To that effect, we look for the parametric curve C(t) that verifies back propagation equation [START_REF] Cohen | Global minimum for active contour models: A minimal path approach[END_REF]:

dC dt = - --→ ∇U C(0) = v 0 .
We solve this equation by a numerical method such as Runge Kutta 4 [START_REF] Press | Numerical Recipes in C : the art of computer programming[END_REF]. To that end, we need to calculate the value of --→ ∇U at any given point on the surface. To obtain a continuous variation of --→ ∇U over the surface, we introduce a novel interpolation scheme that differs from the one used in [START_REF] Sethian | Computing geodesic paths on manifolds[END_REF]. First, we calculate the value of the gradient of U at every vertex v on the surface. Then, we linearly interpolate the value of --→ ∇U on each face. To estimate the value of --→ ∇U at a vertex v, we consider the 1-ring represented in figure 4 (left), where the v i 's are the neighbors of v. By replacing each angle α i by α i in order to have a sum of 2π, we flatten the neighborhood of v on a plane (figure 4, on the right). We can thus approximate the

⇒ v v 1 v 2 v n v 1 v 2 v v n α 1 α 1
∇U (i) = U (v i ) -U (v) |v i -v| .
The gradient at v i satisfies the n equations:

∇U, e i = ∇U (i) ,

where e i := v-vi | v-vi| .
We solve this over-determined system by least squares, which correspond to solving the 2 × 2 system:

e i0 e i0
e i0 e i1 e i0 e i1 e i1 e i1

(∇U ) 0 (∇U ) 1 = e i0 ∇U (i) e i1 ∇U (i) ,
where (e i0 , e i1 ) represent the coordinates of e i in a local frame of the flattening plane.

Figure 5 shows the geodesic path calculation for different models. We can also start several fronts and make them

Isotropic remeshing of a triangulation

Our approach iteratively adds new vertices based on the geodesic distance on the surface. The result of the algorithm gives a set of vertices uniformly distributed on the surface according to the geodesic distance. Taking into account a local density of vertices will be done in sections 3.3 and 3.4.

This is an extension of an idea of the algorithm of [5] that updated a front propagation from a set of points in order to iteratively find pairs of points to link together. In our context, there are two main differences. The first is that we find points on a triangulated surface instead of an image. The second is that we give a precise approximation of the point where the maximal geodesic distance is reached.

Iterative choice of basis points

We now describe how to automatically build an evenly spaced set of points on a triangulated surface. A first point x 1 is chosen at random on the mesh and its geodesic distance map U 1 computed by fast marching. A more elaborate choice consists in replacing this random point by the point with maximum distance from it.

Then we assume we have already computed a set of points S n = {x 1 , . . . , x n }, together with U n the geodesic distance map to S n . To add a new point x n+1 , we simply select a point on the manifold that is furthest away from S n . To compute the new distance map U n+1 , we use the fact that U n+1 = min(U n , U xn+1 ), where we have noted U xn+1 the distance map to x n+1 . So we simply need to update U n by starting a front from x n+1 (using the Fast Marching algorithm exposed in section 2) and to confine it on the set {x ; U xn+1 (x) ≤ U n (x)}. This assures that the whole remeshing process roughly takes less than O(N log(N ) 2 ) operations (which would be the case if we desire as many points as the original mesh has, i.e. N ).

At each iteration, the new point x n+1 needs not to be a vertex of the original mesh. It can be positioned accurately by interpolating the distance map. To be more precise, we add a new vertex v = x n+1 in one of the triangles v 1 , v 2 , v 3 around the vertex of greatest distance v 1 . We choose v 2 the point around v 1 of greatest geodesic distance, and similarly for v 3 but on the neighbors of the edge [v 1 , v 2 ]. It happens most often that the three vertices v i are reached by three different fronts (the other cases are trivial). We denote t 1 , t 2 , and t 3 the arrival times of the three fronts at the three vertices. We calculate the point w 3 at which the fronts 1 and 2 meet by linearly interpolating the distance map along edge [v 1 , v 2 ] (here we assume a front speed of 1):

w 3 := λv 1 + (1 -λ)v 2 with: λ := |v 1 -v 2 | + t 2 -t 1 2|v 1 -v 2 |
We then choose for v the center of mass of the three intersection points w i , as shown in figure 7. Note that v is also added to the set of vertices of the original mesh. We choose to stop the algorithm either when the last added point x n+1 satisfy U n (x n+1 ) ≤ δ, where δ is a given threshold, or when a given number of points have been distributed.

Figure 8 shows three stages of the process of inserting points on a square. (in addition to the triangulation obtained, which is explained in the next section). The intensity reflects the geodesic distance to the nearest basis point (white for value 0).

Calculation of the geodesic triangles

Once we have found the complete set S n0 , we must determine which vertices to link together to obtain our basis triangulation T which is built incrementally during the algorithm. To that end, during the point distribution process we keep track of saddle points, which are vertices v that satisfy these two criterions: -When the value of U (v) is set by the Fast Marching algorithm, two fronts coming from different basis points x i and x j must meet for the first time at v (see [START_REF] Cohen | Multiple contour finding and perceptual grouping using minimal paths[END_REF]). -Adding edge [x i , x j ] to the basis triangulation T must keep the triangulation valid (e.g. the edge must have less than two adjacent faces). Note that when we update a distance map U n+1 , a previously found saddle point v can disappear (if v is reached by the front coming from x n+1 ), and of course new saddle points can be created.

The set of saddle points tells us which vertices should be linked together to obtain a valid triangulation T . We can also trace on the original mesh the geodesic path corresponding to the edges of T . If we have a saddle point v where two fronts from x i and x j meet, we simply perform the gradient descent described in section 2.3 from v in both the direction of x i and the direction of x j .

Figure 9 shows progressive remeshing of the bunny. At last, we can note that some holes might still be found in the resulting triangulation, for instance if there are holes in the original mesh. If we want to avoid such holes (e.g. if we want to use the resulting mesh for parameterizing the manifold), we must fill in the gaps. This can be performed by first unfolding and projecting the hole on a plane, and then computing a constrained Delaunay triangulation. This method is illustrated in figure 10.

Adaptive remeshing

In the algorithm presented in sections 3.1 and 3.2, the fronts propagate at a constant speed which results in uniformly spaced mesh. Our approach can be extended to work with an arbitrary speed function F > 0 for the front propagation, thus computing a geodesic distance weighted by P = 1/F . Since vertices are added at maximal values of the geodesic weighted distance, the resulting triangulation will be dense in regions with smaller F , and in regions with higher F the triangulation will be sparse. This is due to the fact that the algorithm distributes points in such a way that their weighted geodesic distances to neighbors are almost equal. The geodesic distance to vertices in a region with higher value of P is thus smaller. Function F can reflect the need of the user to refine some specific regions with more vertices.

Figure 11 shows a uniform distribution of points on the head and the distribution of points with a split of the mesh into two regions, one with high F and the other with low F . Similar results are shown for the bunny on the left and middle images of figure 13.

When a mesh is obtained from range scanning, a picture of the model can be mapped onto the 3D mesh. Using a function F that is inversely proportional to the norm of the gradient of the image, the user can refine regions with high variations in intensity. Figure 12 shows such a 3D model.

Curvature-based remeshing

The local density of vertices can also reflect some geometric proprerties of the surface. The most natural choice is to adapt the mesh in order to be finer in regions where the local curvature is larger. To that end we choose the speed F to be inversely proportional to max(|λ 1 |, |λ 2 |), where λ i are the eigenvalues of the local curvature tensor. This ten- sor is computed in a pre-processing step using least square fitting of a quadric. The evaluation of the curvature tensor is a vast topic we used a robust construction proposed recently [START_REF] Cohen-Steiner | Restricted Delaunay triangulations and normal cycles[END_REF]. Figure 13 shows on the right a curvaturebased distribution of points on the bunny. Figure 14 shows remeshing of two models with sharp features. Clearly the curvature-based method gives better reconstruction results.

Application to mesh parameterization 4.1 Problem statement

A parameterization f of a 3D triangulated manifold V corresponds to a set of mappings f i : U i → V, where the U i ⊂ R 2 , and the f i are homeomorphisms whose images cover all of V. To simplify, we will consider a triangulated manifold V 0 that contains very few triangles, but that has the same topology as V. The maps f i will then be piecewise linear, and the U i are the triangles of V 0 . Figure 15 shows such a parameterization. The functions f -1 i are affine on every small triangle of V, and the image of each triangle U i of V 0 by f i delimits a region of V. As basis domain for the mapping f we choose the basis mesh constructed with our algorithm in section 3. So the borders of each region f i (U i ) are geodesic paths, and these paths need to be added to the original mesh. Figure 16 shows how to perform this inclusion. Since inside each small triangle of a region f i (U i ) the function f -1 i is affine, for every vertex x of V, we only need to know the barycentric coordinates of f -1 i (x) in the triangle U i . The remaining question is how to calculate those barycentric coordinates judiciously.

Parameterization interpolation

Last, we need to calculate the barycentric coordinates (ν 1 , ν 2 , ν 3 ) of f -1 (v) inside a given triangle x 1 , x 2 , x 3 , for each vertex v. Note, the points v 1 = f (x 1 ), v 2 = f (x 2 ), and v 3 = f (x 3 ) are the basis points determined in paragraph 3.1. To that end, we will use the geodesic distances (t 1 , t 2 , t 3 ) between v and each of the v i . We will in fact place the point f -1 (v) inside the triangle x 1 , x 2 , x 3 whith conservation of the calculated geodesic distances.

Therefore, knowing the distances (t 1 , t 2 , t 3 ) of a point x to the three vertices (x 1 , x 2 , x 3 ) of the triangle, we need to calculate the barycentric coordinates (ν 1 , ν 2 , ν 3 ) of x in a manner analogous to the case of a planar triangle. We denote (l 1 , l 2 , l 3 ) the lengths of the geodesic edge opposite to each point (for example, l 1 is given by the arrival time in v 3 of the front departing from v 2 ). Now, in a plane, Heron's formula says that the area A 3 of the triangle x, x 1 , x 2 is:

A 3 = p(p -l 3 )(p -t 1 )(p -t 2 ),
with p := (t 1 + t 2 + l 3 )/2. Figure 17 shows the triangle considered. Similarly, we can calculate A 1 and A 2 , which

t 1 t 2 t 3 l 1 l 2 l 3 x 1 x 2 x 3 A 3
ν 1 = A 1 A 1 + A 2 + A 3 ν 2 = A 2 A 1 + A 2 + A 3 ν 3 = A 3 A 1 + A 2 + A 3
Those formulas give us a parameterization of the point x, as long as there is no other point that has the same values for (t 1 , t 2 , t 3 ). This condition is met in practice as long as V 0 is topologically equivalent to V, and that the basis meshing is fine enough.

To compute the three distances (t 1 , t 2 , t 3 ) for each vertex, we simply start a front and confine it to the geodesic neighborhood of each basis vertex, as depicted on the left of figure 18. Doing this, we are sure that each point will be reached by three fronts and only three. On the right of figure 18, we can see the geodesic parameterization obtained. We simply assign a random color to each vertex, and interpolate the color using the parameterization.

Application to multiresolution constructions

We can sample each basis triangle according to a regular subdivision scheme such as those depicted in figure These kinds of triangulations are of primary importance, both because the connectivity is known in advance (which is a great advantage for compression), and because the subdivision process gives a natural multiresolution representation of the mesh. Wavelet transforms can then be built using the lifting scheme [START_REF] Sweldens | The Lifting Scheme: A construction of second generation wavelets[END_REF]. Wavelet coefficients can be encoded using a zerotree coder [START_REF] Guskov | Normal meshes[END_REF], which gives good distortion results.

Figure 20 show the process of subdivision. The parameterization of the mesh is shown in figure 18 (bottom right).

Final remarks and future work

The running times of our algorithm are very reasonable, and it takes about 10 seconds to distribute 500 points on a model of about 5000 vertices. The distribution of more points is almost immediate, thanks to the incremental structure of our method. The coarse mesh construction and the parameterization steps take about the same time. So our current implementation allows us to process a whole model in less than one minute. The isotropic remeshing algorithm proposed here gives very good results, similar to those given in [START_REF] Terzopoulos | Adaptive meshes and shells: Irregular triangulation, discontinuities, and hierarchical subdivision[END_REF] for a 2D mesh. The advantage of our method is that it is conceptually simple and very fast. The main component is the Fast Marching algorithm, which is pretty well understood.

The mesh parameterization scheme proposed in this paper gives a smooth map in regions of the manifold that are smooth. Tests and distortion measures remain to be made to compare this method with classical ones, such as Floater's scheme [START_REF] Floater | Parameterization of manifold triangulations[END_REF].

Our future directions of research include taking into account special features (such as sharp edges) during the remeshing and the parameterization. One could define a linear density of vertices on these 1D features. An interesting question is finding theoretical bounds for the parameterization distortion in case of smooth surfaces.

Conclusion

We have described an iterative algorithm for distributing a set of points on a surface according to a local density function. This algorithm is fast and efficient, and provides a simple way to automatically find a geodesic Delaunay triangulation of the original surface. We have applyed this construction to build a coarse basis domain for a parameterization of the triangulated manifold. This method provides a smooth mapping from the basis mesh to the original triangulation, allowing multiresolution constructions and wavelet compression.

  Figure 1 shows the different steps involved in the parameterization process: -Original mesh. -Automatic determination of basis points. -Determination of associated Delaunay triangulation. -Calculation of the corresponding geodesic triangles. -Distance data calculation. -Parameterization interpolation. -Semi-regular remeshing.

Figure 2 .

 2 Figure 2. Orthogonal and arbitrary triangulations.

Figure 4 .

 4 Figure 4. Flattening of a neighborhood. gradient along each edge (v, v i+1 ):

Figure 5 .

 5 Figure 5. Geodesic paths. evolve at the same time, as it is shown in figure 6. The different colored regions form a Voronoi diagram of the starting points, and the intensity represents the geodesic distance to the closest point.

Figure 6 .

 6 Figure 6. Front propagation (top), geodesic distance map (bottom, left), and Voronoi diagram (bottom, right).

3 t 1 = 1 t 2 = 0. 9 t 3 = 0. 8 vFigure 7 .

 31129387 Figure 7. Determination of intersection points of three propagating fronts.

Figure 8 .

 8 Figure 8. Insertion of points in a square (top), and corresponding triangulation (bottom).

Figure 9 .

 9 Figure 9. Geodesic remeshing with 100, 300, 800 and 1500 points.

Figure 10 .

 10 Figure 10. Original hole, flattened hole, and constrained Delaunay triangulation.

Figure 11 .

 11 Figure 11. Uniform (left) and user-defined remeshing (right).

Figure 12 .

 12 Figure 12. Original model, F function, adaptive distribution of points, and triangulation. Thanks to authors of [22] for providing data.

Figure 13 .

 13 Figure 13. Uniform, user defined, and curvature-based distribution of points.

Figure 14 .

 14 Figure 14. Original model (left), uniform remeshing (middle) and curvature-based remeshing (right).

Figure 15 .

 15 Figure 15. Piecewise linear parameterizationAs basis domain for the mapping f we choose the basis mesh constructed with our algorithm in section 3. So the borders of each region f i (U i ) are geodesic paths, and these paths need to be added to the original mesh. Figure16shows how to perform this inclusion.

Figure 16 .

 16 Figure 16. Triangulation before inclusion of the path (left) and after (right).

Figure 17 .

 17 Figure 17. Determination of the barycentric coordinates of a point x gives us the following barycentric coordinates:

Figure 18 .

 18 Figure 18. Examples of parameterizations.19. The resulting triangulation is a so-called semi-regular mesh.

Figure 19 .

 19 Figure 19. Triangle-based and quad-based subdivision schemes.

Figure 20 .

 20 Figure 20. Four steps of regular subdivision.