namespace CGAL {
/*!
\mainpage User Manual
\anchor Chapter_Triangulated_Surface_Mesh_Approximation

\cgalAutoToc
\authors Pierre Alliez, David Cohen-Steiner, Lingjie Zhu

\section sma_introduction Introduction

This package implements the <em>Variational Shape Approximation</em> \cgalCite{cgal:cad-vsa-04} (%VSA) method to approximate an input surface mesh by a simpler surface triangle mesh.
The input of the algorithm must be:
- Triangulated
- Combinatorially 2-manifold

The output is a triangle soup and can be built into a polygon surface mesh.

Given an input surface triangle mesh, %VSA leverages a discrete clustering algorithm to approximate it by a set of local simple shapes referred to as proxies. Each cluster is represented as a connected set of triangles of the input mesh, and the output mesh is constructed by generating a surface triangle mesh which approximates the clusters.
The approximation error is one-sided, defined between the clusters and their associated proxies. Two error metrics (\f$ \mathcal{L}^2 \f$, \f$ \mathcal{L}^{2,1} \f$) for planar proxies are provided via the classes `CGAL::Surface_mesh_approximation::L2_metric_plane_proxy` and `CGAL::Surface_mesh_approximation::L21_metric_plane_proxy`, and the algorithm design is generic to other user-defined metrics. The current proxies are planes or vectors, however the algorithm design is generic for future extensions to non-planar proxies \cgalCite{cgal:ywly-vmsqs-12}\cgalCite{cgal:wk-srhvs-05}. The default \f$ \mathcal{L}^{2,1} \f$ metric is recommended in terms of computation and visual perception \cgalCite{cgal:cad-vsa-04}. A brief background about <em>%Proxy</em> and <em>%ErrorMetric</em> is provided in Section \ref sma_background.

\cgalFigureBegin{Approximation_teaser, teaser.jpg}
Variational shape approximation on two models with the \f$ \mathcal{L}^{2,1} \f$ error metric and planar proxies. From left to right: partition of the input surface triangle mesh, anchor vertices and edges, and output triangle mesh. The partition is optimized via discrete clustering of the input triangles, so as to minimize the approximation error from the clusters to the planar proxies (not shown).
\cgalFigureEnd

This package offers both the approximation and mesh construction functionalities, through the free function `CGAL::Surface_mesh_approximation::approximate_triangle_mesh()` which runs a fully automated version of the algorithm:
\cgalExample{Surface_mesh_approximation/vsa_simple_approximation_example.cpp}

A class interface is also provided for advanced users, in which a series of pliant operators offer interactive capabilities during clustering and customization in terms of error and proxies.

\section sma_overview Overview

The package contains 3 main components: approximation algorithm, pliant operators and meshing as shown in Figure \cgalFigureRef{workflow}.

\cgalFigureBegin{workflow, workflow.svg}
From left to right are 3 components of the approximation package: approximation algorithm (left), optional pliant operations (middle) and meshing (right).
\cgalFigureEnd

\subsection sma_approximation Approximation

The left part of Figure \cgalFigureRef{workflow} depicts the workflow of the approximation algorithm.

\subsubsection sma_clustering Clustering Iteration

Figure \cgalFigureRef{iterations} depicts several Lloyd \cgalCite{cgal:l-lsqp-82} clustering iterations on the plane-sphere model with planar proxies and the \f$ \mathcal{L}^{2,1} \f$ metric. We plot the fitting error against each iteration. After 8 iterations, the error barely changes. Based on this observation, we consider that the clustering converges if the error change between the current and previous iteration is lower than a user-specified threshold (indicated by two green dash lines).

\cgalFigureBegin{iterations, iterations.jpg}
Discrete Lloyd iterations on the plane-sphere model with planar proxies and the \f$ \mathcal{L}^{2,1} \f$ metric: (left) random seeding of 6 proxies; (center) after one iteration; (right) after 8 iterations, the regions settle. The red lines depict the proxy normals drawn at the seed faces.
\cgalFigureEnd


\subsubsection sma_seeding Seeding

Each proxy is always associated to a <em>seed</em> triangle face in the input surface mesh. While the proxies may be viewed as centers (or best representative) in a geometric error sense, the seed of each proxy is used as the starting point in the clustering process. Seeding is the processing of deciding how to select a seed face where a new proxy/partition can be initialized from.

To start the clustering iterations, we need an initial set of proxies. The default (hierarchical) approach generates one proxy per connected component, seeded at arbitrarily chosen faces. It then adds more proxies in batches, in order to drive the error down. After each batch of proxies added, it performs several inner clustering iterations, which is referred to as <em>relaxation</em> in the seeding step.

Assuming a clustering partition of \f$n\f$ regions with errors \f$ \{E_k\}_{k=1\cdots n} \f$, and we wish to add \f$m\f$ proxies. We provide 3 different seeding methods:
- <b>Random</b>. \f$m\f$ seed faces are picked randomly on the surface, excluding the current seed faces.
- <b>Incremental</b>. Each new proxy is initialized from a face of the region with the largest approximation error. The face itself is chosen as the one realizing the largest error in its region.
- <b>Hierarchical</b>. \f$m\f$ seed faces are dispatched within the current partition, where each region is refined with a number of proxies chosen in accordance to its fitting error:
  - calculate total error \f$ E_{total} \f$, then average error \f$ E_{avg} = E_{total} / m \f$ (assuming that each new proxy shares the same amount of error)
  - sort errors \f$ \{E_{min},\cdots,E_{max}\} \f$
  - from \f$ E_{min} \f$ to \f$ E_{max} \f$, we diffuse the error <em>hierarchically</em> one after another. More specifically, the number of proxies \f$N_k\f$ added to the \f$k\f$th region is proportional to its error:
\f[ N_k = \lfloor E_k / E_{avg} + 0.5 \rfloor, \f]
and the remaining error is added to the next proxy error in order to keep the total error unchanged:
\f[ E'_{k+1} = (E_k - N_k * E_{avg}) + E_{k+1}. \f]


Figure \cgalFigureRef{seeding_method} depicts different seeding methods. Random initialization randomly selects a set of input triangle faces as proxy seeds. While it is very fast, the subsequent clustering process can be entrapped in a bad local minimum, especially on shapes with regions surrounded by sharp creases (left closeup). Incremental initialization adds the proxies one by one at the most distorted region. It can thus be slow due to the large number of interleaved relaxation iterations. Hierarchical initialization (selected by default) repeatedly doubles the number of proxies in a hierarchical refinement sequence, so as to generate clustering regions with evenly distributed fitting errors. Time consumption is typically in-between the former two. Statistics and comparisons are available in Section \ref sma_perf.

\cgalFigureBegin{seeding_method, seeding_method.jpg}
Different seeding methods on the sphere-cube model. From left to right: initial partition (\f$ \mathcal{L}^{2,1} \f$ metrics and 20 proxies), add 5 proxy seeds (red faces) with random, incremental and hierarchical methods respectively.
\cgalFigureEnd


\subsubsection sma_stop Stop Criteria

To determine when to stop adding more proxies, we can specify either the maximum number of proxies required to approximate the geometry or the minimum error drop percentage with respect to the very first partition. More specifically, we can decide:
- <b>Maximum number of proxies</b>. Adding proxies until the specified number is met.
- <b>Minimum error drop</b>. Starting from the very first partition (with one proxy per connected component) with fitting error \f$ \hat E \f$, the algorithm adds proxies until the approximation error drops below the specified percentage \f$ target\_drop * \hat E \f$.
As depicted by Figure \cgalFigureRef{nb_proxies}, specifying a minimum error drop of 10% (yellow dash lines) as stopping criterion, yields 12 proxies on the plane-sphere model. When both criteria are provided, the first criterion met stops the seeding. Different seeding examples are depicted by Figure \cgalFigureRef{meshing}.
To balance between performance and speed, we recommend using hierarchical seeding and specifying both stopping criteria.


\cgalFigureBegin{nb_proxies, nb_proxies.jpg}
Using different number of proxies to approximate the plane-sphere model. From left to right: 8, 14, 20 proxies. We plot right the error against the number of proxies.
\cgalFigureEnd


\subsection sma_operations Pliant Operators

For interactive use, the approach can yield better approximations of the geometry via adding/removing proxies and tunneling out of local minima via additional operations:
- <b>Merging</b>. Merging two adjacent regions.
- <b>Splitting</b>. Splitting one specified region into smaller ones to reduce the error. By default bisection is performed but N-section is also possible. We first select the request number of face seeds from the specified region, then perform the refitting process confined to the region.
- <b>Adding</b>. Adding one or more proxies to further reduce the approximation error. As for the seeding process, addition can be performed incrementally or hierarchically, as shown in Figure \cgalFigureRef{seeding_method}.
- <b>Teleporting</b>. A teleportation operator is a combination of merging and adding proxies: merging the pair of adjacent regions and adding a proxy seed to the worst region. More specifically, the pair of regions whose merging realizes the smallest error after merging and local re-fitting, is selected for merging. In practice, the teleport operation can temporarily either decrease or increase the total approximation error. We provide an optional heuristic to evaluate if the teleportation is worth it by verifying whether the error increase induced by a (simulated) deletion is smaller than half of the error of the worst region. If this test fails, no teleportation is judged necessary.

\cgalFigureBegin{operations, operations.jpg}
Operations on the sphere-cube union model. Upper row: merging (9 proxies) and the reverse bisection splitting (10 partitions) on the confined region. Lower row: one teleportation operation of merging and adding a face seed, the sphere is approximated with one more proxy.
\cgalFigureEnd

As depicted in Figure \cgalFigureRef{operations}, teleportation provides a means to relocate a local minimum region entrapped in the planar part (left) to the most needed regions on the sphere (right). In \ref sma_example3, the class interface is used to control the approximation process through the aforementioned operations.

\subsection sma_meshing Meshing

This package implements the meshing algorithm described in \cgalCite{cgal:cad-vsa-04} by generating a triangle mesh approximation of the clustering partition. The meshing algorithm has two major steps:
- <b>Finding anchors</b>. <em>Anchors</em> are just a subset of the vertices of the input mesh on the region boundaries.
- <b>Discrete constrained 2D Delaunay triangulation</b>. Connecting the anchors to generate the final approximated triangle mesh.

\subsubsection sma_anchors_basic Basic Anchors

A vertex is selected as a basic anchor if it is:
- not a mesh boundary vertex and adjacent to at least three regions,
- or a mesh boundary vertex adjacent to at least two regions.

\subsubsection sma_anchors_subdivision Subdivision Anchors

Walking along the boundary of a proxy region  (counterclockwise), a <em>chord</em> is a sequence of halfedges connecting two anchors. One cluster boundary cycle may consist of several chords. A connected region with holes may yield several boundary cycles (Figure \cgalFigureRef{operations}, planar part before teleportation).

In order to approximate complex boundaries well, more anchors are generated by recursive chord subdivision (Figure \cgalFigureRef{chord}). An anchor \f$ \mathbf{c} \f$ is added at the furthest vertex of a chord \f$ (\mathbf{a}, \mathbf{b}) \f$, split it into \f$ (\mathbf{a}, \mathbf{c}) \f$ and \f$ (\mathbf{c}, \mathbf{b}) \f$, if the distance \f$ d = \Vert \mathbf{c} , (\mathbf{a}, \mathbf{b}) \Vert \f$ is beyond certain threshold. To make \f$ d \f$ independent to the scale of the input:
\f[ d = d / input\_mesh\_average\_edge\_length. \f]
Optionally, \f$ d \f$ can be measured as the ratio of the chord length:
\f[ d = d / \Vert(\mathbf{a}, \mathbf{b})\Vert. \f]
Also, we can add a dihedral angle weight \f$ sin(\mathbf{N}_i,\mathbf{N}_j) \f$ to the distance measurement, where \f$ \mathbf{N}_i,\mathbf{N}_j \f$ are the normals of the proxies separated by the chord \f$ (\mathbf{a}, \mathbf{b}) \f$. If the angle between proxy \f$ P_i \f$ and \f$ P_j \f$ is rather small, then a coarse approximation will do as it does not add geometric information on the shape. Trivial chords (less than 4 edges) are not subdivided if they are non-circular. In case of circular chords, additional anchors maybe added to maintain the topology even if they are trivial, as detailed in Section \ref sma_anchors_additional.

\cgalFigureBegin{chord, chord.jpg}
Varying the chord error. From left to right: clustering partition, and meshing with decreasing absolute chord error 5, 3 and 1 without dihedral angle weight. The boundaries of the partition (red lines) are approximated with increasing accuracy.
\cgalFigureEnd

\subsubsection sma_anchors_additional Additional Anchors

For a boundary cycle without any anchor such as the hole depicted Figure \cgalFigureRef{operations}, we first add a starting anchor to the boundary. We then subdivide this circular chord to ensure that every boundary cycle has at least 2 anchors (i.e., every chord is connecting 2 different anchors, Figure \cgalFigureRef{anchors}). Finally, we add additional anchors to ensure that at least three anchor vertices are generated on every boundary cycle.

\cgalFigureBegin{anchors, anchors.jpg}
Adding anchors. From left to right: starting from a partition (grey) with a hole (while) and two encircled regions (green and blue), we add a starting anchor (orange disk) to the boundary cycle (red dash line) without any anchor (2nd), subdivide the circular chord (3rd, the number indicates the level of recursion) and add anchors to the boundary cycle with less than 2 anchors (4th, red dash lines).
\cgalFigureEnd

\subsubsection sma_triangulation Discrete Triangulation

With the anchors defined, their chord connection graph forms a general polygon mesh. Because of non-flat, concave polygon or polygons with holes, we need to triangulate this initial polygon mesh. The triangulation is generated by computing a discrete variant of a constrained 2D Delaunay triangulation, where distances are measured on the input triangle mesh.

The first image of Figure \cgalFigureRef{triangulation} depicts how the Delaunay triangulation of set of points (colored disks) is deduced from its dual Voronoi diagram (colored regions separated by blue lines) by connecting the points indicated by the vertices (red circles) where 3 Voronoi cells meet. In an analogous manner, we construct discrete Voronoi cells from which the triangulation is extracted.


In a first step, we start a flooding of the interior of the region, coloring the vertices according to their closest anchor vertex. We then only flood the boundary of a region so that every vertex on it is colored depending on the closest anchor vertex. This enforces the constrained edges by forcing the boundary to be in it.

\cgalFigureBegin{triangulation, triangulation.jpg}
Discrete constrained triangulation on a sphere model. The triangulation process first floods the inner vertices (red arrows, 2nd) to simulate the Voronoi diagram. It then constructs constrained edges between anchor vertices, by flooding along the boundary edges (red arrows, 3rd). Finally, triangles (red hollow triangle, 4th) are formed by connecting the source anchors (black arrows, 4th) of the faces where 3 Voronoi cells meet (red solid triangle, 4th).
\cgalFigureEnd

After each clustering region is triangulated, the final anchor positions are recomputed by averaging the projections of an anchor on its incident proxies.

In \cgalFigureRef{meshing}, the bear model is approximated through \f$ \mathcal{L}^{2,1} \f$ metric and the final number of proxies is determined by monitoring the error drop.
The anchor points (black) are attached to the corresponding vertex on the mesh (white). The red lines connecting the anchor points approximate the boundary of each region.

\cgalFigureBegin{meshing, meshing.jpg}
Meshing the bear model with decreasing target error drop.
From left to right, the target error drop are 6\%, 4\% and 2\% to the initial error respectively, the output mesh densifies. Notice the boundary subdivision in the black rectangle area.
\cgalFigureEnd

As there is no guarantee that the output mesh is 2-manifold and oriented, the main input is an indexed triangle set. We can use \ref PolygonSoups to build the triangle soup into a valid polygon mesh.

\subsection sma_api API

This package can be used with any class model of the concept `FaceListGraph` described in \ref PkgBGL "CGAL and the Boost Graph Library".

Free function with \ref vsa_namedparameters options.
- `CGAL::Surface_mesh_approximation::approximate_triangle_mesh()`: given a triangle mesh, approximate the geometry with default \f$ \mathcal{L}^{2,1} \f$ metric.

Class interface:
- `CGAL::Variational_shape_approximation`: allowing more customization of the proxy, metric and approximation process. As shown in Figure \cgalFigureRef{workflow}, typical calling order of the approximation and meshing process is:
  - \link CGAL::Variational_shape_approximation::initialize_seeds initialize seeds \endlink
  - \link CGAL::Variational_shape_approximation::run run clustering iterations \endlink
  - \link CGAL::Variational_shape_approximation::extract_mesh extract mesh \endlink
  - \link CGAL::Variational_shape_approximation::output take outputs \endlink

One thing to note is that some parameters depend heavily on the input, like the number of proxies. Although we can approximate a geometry with any number of proxies regardless of the quality, it is not recommended to use all the defaults without any consideration of the input.

\section sma_examples Examples

\subsection sma_example1 Free Function Approximation

The following example calls the free function `CGAL::Surface_mesh_approximation::approximate_triangle_mesh()` on the input triangle mesh with default `CGAL::Surface_mesh_approximation::L21_metric_plane_proxy`.

\cgalExample{Surface_mesh_approximation/vsa_approximation_example.cpp}

The function parameters are provided through \ref vsa_namedparameters. Setting the non-default parameter values requires calling the functions with the required parameters, connected by a dot and in an arbitrary order. The following example shows a different configuration of parameters and retrieves the indexed face proxy map and the proxies:

\cgalExample{Surface_mesh_approximation/vsa_approximation_2_example.cpp}

The face proxy index map and the output proxies provide a means to access the partition and the proxy parameters as illustrated by \cgalFigureRef{iterations}.

\subsection sma_example2 Free Function Segmentation

The free function can be used for retrieving the segmentation via proxy ids output into the proxy output iterator:

\cgalExample{Surface_mesh_approximation/vsa_segmentation_example.cpp}

\subsection sma_example3 Class Interface

The class interface `CGAL::Variational_shape_approximation` offers a flexible means to control of the algorithm. The following example uses the `CGAL::Surface_mesh_approximation::L2_metric_plane_proxy` to approximate the shape.

\cgalExample{Surface_mesh_approximation/vsa_class_interface_example.cpp}

\subsection sma_example4 Customized Proxy and Error Metric

\cgalFigureBegin{vsa_metric_comparison, vsa_metric_comparison_200_30.jpg}
Comparison of different error metrics on the bear model, with 200 proxies and hierarchical seeding. From left to right: \f$ \mathcal{L}^{2,1} \f$ metric, \f$ \mathcal{L}^2 \f$ metric and custom compact metric.
\cgalFigureEnd

The following example defines a point-wise proxy to yield an isotropic approximation. The output mesh is depicted in Figure \cgalFigureRef{vsa_metric_comparison}.

\cgalExample{Surface_mesh_approximation/vsa_isotropic_metric_example.cpp}

\section sma_perf Performances

We provide some performance comparisons with the free function API `CGAL::Surface_mesh_approximation::approximate_triangle_mesh`.
Timings are recorded on a PC running Windows10 X64 with an Intel Xeon E5-1620 clocked at 3.70 GHz with 32GB of RAM.
The program has been optimized with the O2 option with Visual Studio 2015. By default the kernel used is `Exact_predicates_inexact_constructions_kernel` (`EPICK`).

Runtime in seconds with target number of proxies of different seeding method:

<center>
    Model    | \#Triangles | \#Proxies | Random | Incremental | Hierarchical
----------:  | ---------:  | -------:  | ----:  | ---------:  | -----------:
plane-sphere |    6,826    |     20    |     0  |     0.87    | 0.17
bear         |   20,188    |    200    |     0  |   36.749    | 1.194
masque       |   62,467    |    200    | 0.002  |  133.901    | 4.308
</center>

Runtime in seconds with target error drop of different seeding method. The benchmark is running on the bear model with 20,188 faces. Each column records the time and the resulting number of proxies:

<center>
 Target Error Drop |   Random   | Incremental | Hierarchical
----------------:  | --------:  | ---------:  | -----------:
        0.06       |   1.03/64  |   9.053/53  | 1.017/64
        0.04       |  1.207/128 |  15.422/88  | 1.2/128
        0.02       |  1.415/256 |  35.171/192 | 1.428/256
</center>

Runtime of the 3 phases of the algorithm in seconds: seeding, clustering iteration and meshing. The seeding method is hierarchical with target number of proxies.

<center>
    Model    | \#Triangles | \#Proxies | \#Iterations | Seeding | Clustering | Meshing | Total
----------:  | ---------:  | -------:  | ----------:  | -----:  | --------:  | -----:  | ----:
plane-sphere |    6,826    |     20    |       20     |   0.17  |    0.228   | 0.044   | 0.442
bear         |   20,188    |    200    |       20     |  1.194  |    0.784   | 0.128   | 2.006
masque       |   62,467    |    200    |       20     |  4.308  |    2.974   | 0.349   | 7.631
</center>

<!-- With different metrics, we record the running time and complexity of the output mesh.

<center>
    Model    | \#Triangles | \f$ \mathcal{L}^{2,1} \f$ | \f$ \mathcal{L}^2 \f$ | %Compact Metric
----------:  | ---------:  | -----------------------:  | -------------------:  | -------------:
plane-sphere |    6,826    |           tbd             |           tbd         | tbd
bear         |   20,188    |           tbd             |           tbd         | tbd
masque       |   62,467    |           tbd             |           tbd         | tbd
</center> -->

\section sma_background Background

The %VSA method has two key geometric concepts:

- <b>Proxy</b> \f$ P \f$. The parameterized best-fit surrogate of a piece of surface geometry.

- <b>Error metric</b> \f$ E \f$. To measure how well a proxy approximates the corresponding geometry.

Given an error metric \f$ E \f$, a desired number of \f$ k \f$ proxies,
and an input surface \f$ S \f$, we denote by <em>optimal shape proxies</em> a set \f$ P \f$ of proxies \f$ P_i \f$ associated to the regions \f$ R_i\f$ of a partition \f$ \mathcal{R} \f$ of \f$ S \f$ that minimizes the total fitting error \cgalCite{cgal:cad-vsa-04} :

\f[ E(\mathcal{R}, P) = \sum_{i = 1..k} E(\mathcal{R}_i, P_i). \f]

By casting the approximation problem into an optimal discrete clustering one, the algorithm leverages the effective Lloyd algorithm \cgalCite{cgal:cad-vsa-04} to drive the total error down iteratively.
More specifically, during each iteration two different steps are conducted, for the \f$ m \f$th iteration:
- <b>Partition process</b>. Firstly, all triangle faces are partitioned into \f$ k \f$ connected regions \f$ \{ R^{m}_1, \cdots, R^{m}_k \} \f$ by assigning each face to its nearest proxy \f$ P^{m-1}_i \f$.
- <b>Fitting process</b>. Then, the algorithm fits a proxy and update the parameters \f$ P^{m}_i \f$ from the corresponding region \f$ R^{m}_i \f$.

For a sequence of iterations with the fitting error \f$ \{ E^1, \cdots, E^m \} \f$, the iteration is repeated until one of the stopping criteria is met:
- The maximum number of iterations is reached: \f$ m >= max\_iterations \f$.
- No significant error change between two iterations: \f$ (E^{m-1} - E^{m}) / E^{m-1} < converge\_threshold \f$.

Intuitively, each region \f$ \mathcal{R}_i \f$ of a partition \f$ \mathcal{R} \f$ can be summarily represented to first order as an "average" point \f$ X_i \f$ and an "average" normal \f$ N_i \f$.
We denote such local representative pair \f$ P_i = (X_i, N_i) \f$, a <em>planar proxy</em> of the associated region.

Defining an appropriate error metric is the key ingredient for the algorithm.
The \f$ \mathcal{L}^2 \f$ metric is defined as:
\f[ \mathcal{L}^2(\mathcal{R}_i, P_i) = \iint_{x \in \mathcal{R}_i}\Vert x - \Pi_i(x)\Vert^2 dx. \f] where \f$ \Pi_i(\cdot) \f$ denotes the orthogonal projection of the argument onto the proxy plane passing through \f$ X_i \f$ and normal to \f$ N_i \f$. The \f$ \mathcal{L}^2 \f$ metric tries to match the input shape through approximation of the geometric position.

In the original paper \cgalCite{cgal:cad-vsa-04} the author proposed the \f$ \mathcal{L}^{2,1} \f$ metrics, arguing that the normals are important to the visual interpretation of the shape. The \f$ \mathcal{L}^{2,1} \f$ is defined as: \f[ \mathcal{L}^{2,1}(\mathcal{R}_i, P_i) = \iint_{x \in \mathcal{R}_i}\Vert \mathbf{n}(x) - \mathbf{n}_i\Vert^2 dx. \f]

The \f$ \mathcal{L}^{2,1} \f$ is numerically superior to \f$ \mathcal{L}^2 \f$ in several ways:
- The anisotropy of the surface is better captured.
- Finding the best normal proxy is as simple as averaging the normals.


\section sma_history Design and Implementation History

This package is the result of the work of Lingjie Zhu during the 2017 season of the Google Summer of Code, mentored by Pierre Alliez. The code is based on an initial research code written by Pierre Alliez at Inria in 2003, for a paper published at the ACM SIGGRAPH conference in 2004, co-authored by David Cohen-Steiner, Pierre Alliez and Mathieu Desbrun \cgalCite{cgal:cad-vsa-04}.

*/

}
