Kaimo Hu
email: hukaimo02@gmail.com.

Dong-Ming Yan
email: yandongming@gmail.com.

David Bommes
email: bommes@aices.rwth-aachen.de.

Pierre Alliez
email: pierre.alliez@inria.fr.

Bedrich Benes
email: bbenes@purdue.edu

• K Hu

D.-M Yan

Error-Bounded and Feature Preserving Surface Remeshing with Minimal Angle Improvement

Keywords: error-bounded, feature preserving, minimal angle improvement, feature intensity. !

The typical goal of surface remeshing consists in finding a mesh that is (1) geometrically faithful to the original geometry, (2) as coarse as possible to obtain a low-complexity representation and (3) free of bad elements that would hamper the desired application. In this paper, we design an algorithm to address all three optimization goals simultaneously. The user specifies desired bounds on approximation error δ, minimal interior angle θ and maximum mesh complexity N (number of vertices). Since such a desired mesh might not even exist, our optimization framework treats only the approximation error bound δ as a hard constraint and the other two criteria as optimization goals. More specifically, we iteratively perform carefully prioritized local operators, whenever they do not violate the approximation error bound and improve the mesh otherwise. Our optimization framework greedily searches for the coarsest mesh with minimal interior angle above θ and approximation error bounded by δ. Fast runtime is enabled by a local approximation error estimation, while implicit feature preservation is obtained by specifically designed vertex relocation operators. Experiments show that our approach delivers high-quality meshes with implicitly preserved features and better balances between geometric fidelity, mesh complexity and element quality than the state-of-the-art.

Index Terms-surface remeshing,

INTRODUCTION

S Urface remeshing is a key component in many geom- etry processing applications such as simulation, deformation, or parametrization [START_REF] Botsch | Polygon Mesh Processing[END_REF]. While many remeshing techniques are goal-specified, a common goal is to find a satisfactory balance between the following three criteria:

•

The output mesh should be a good approximation of the input, making the geometric fidelity, usually measured as the approximation error, a key requirement for most applications.

•

The quality of mesh elements is crucial for robust geometry processing and numerical stability of simulations, which requires fairly regular meshes in terms of both geometry and connectivity. Particularly, a lower bound on the minimal angle is vital for many simulation applications [START_REF] Shewchuk | What is a good linear element? interpolation, conditioning, and quality measures[END_REF]. methods that generate meshes with high element quality often require high mesh complexity or introduce high approximation error [START_REF] Alliez | Centroidal Voronoi diagrams for isotropic surface remeshing[END_REF], [START_REF] Valette | Generic remeshing of 3D triangular meshes with metric-dependent discrete Voronoi diagrams[END_REF]. The error-driven methods, while preserving the results in controllable geometric fidelity and low mesh complexity, do not deal with the element quality [START_REF] Cohen-Steiner | Variational shape approximation[END_REF], [START_REF] Mandad | Isotopic approximation within a tolerance volume[END_REF]. In addition, many approaches require the sharp features to be specified or detected in advance [START_REF] Alliez | Interactive geometry remeshing[END_REF]- [START_REF] Yan | Gap processing for adaptive maximal Poisson-disk sampling[END_REF], which is usually difficult and error-prone.

We propose an approach that controls both the approximation error and element quality simultaneously. Our approach only requires the user to specify the error-bound threshold δ, the desired minimal angle θ and an upper bound on the mesh complexity N , measured as the number of vertices. In an initial phase our algorithm coarsens the mesh as much as possible while respecting the error-bound δ. It then iteratively improves the minimal angle of the mesh until the desired bound for θ or N is met. Since all atomic operations respect the error-bound δ, the result is guaranteed to satisfy geometric fidelity. In contrast, low mesh complexity and high quality of mesh elements are optimization goals, which depending on the model and the desired bounds might or might not be met. However, experiments show that given a reasonable vertex budget N and an angle bound θ ≤ 35 • our algorithm is usually able to reach all three goals simultaneously. Moreover, our method preserves geometric features like sharp creases and ridges without explicitly specifying or detecting them.

Our method is inspired by the remeshing methods that proceed by applying a series of local operators such as edge split, edge collapse, edge flip, and vertex relocate [START_REF] Hoppe | Mesh optimization[END_REF], [START_REF] Botsch | A remeshing approach to multiresolution modeling[END_REF]. Contrary to the existing methods that iteratively apply these operations sequentially and globally, our core algorithm employs a dynamic priority queue to maintain all angles which are smaller than the user specified threshold, and then greedily applies local operators whenever an improvement of the mesh is possible. After the initial coarsening our method only improves the mesh in local regions where the element quality is poor, and thus modifies the model as little as necessary with respect to the minimal angle threshold θ and the error-bound constraint δ.

Since we apply local operators directly on the mesh, without using any surface parameterizations or density functions, our algorithm is naturally suitable for high genus, irregular, and multi-component inputs, as well as meshes with very high/low resolutions. In summary, we claim the following contributions:

• A surface remeshing algorithm with minimal angle improvement based on applying local operators, which bounds the "worst element quality" (Sec. 3).

• A reliable and efficient local error update scheme based on approximated symmetric Hausdorff distance, which bounds the geometric fidelity (Sec. 4).

• Two feature intensity functions designed for vertex position relocation, which enable implicit feature preservation and support geometric fidelity (Sec. 5).

RELATED WORK

The variety of applications leads to a large number of different remeshing techniques. We restrict the discussion to the most relevant aspects of our algorithm, i.e. high-quality remeshing, error-driven remeshing and feature-preserving remeshing. For a more complete discussion we refer the reader to the survey [START_REF] Alliez | Recent advances in remeshing of surfaces[END_REF]. High quality remeshing is typically based on sampling and Centroidal Voronoi Tessellation (CVT) optimization [START_REF] Du | Centroidal Voronoi tessellations: applications and algorithms[END_REF]. Early approaches apply 2D CVT in a parametric domain [START_REF] Alliez | Centroidal Voronoi diagrams for isotropic surface remeshing[END_REF], [START_REF] Alliez | Interactive geometry remeshing[END_REF], [START_REF] Alliez | Isotropic surface remeshing[END_REF]- [START_REF] Fuhrmann | Direct resampling for isotropic surface remeshing[END_REF]. Instead of CVT optimization, Vorsatz et al. [START_REF] Vorsatz | Dynamic remeshing and applications[END_REF] utilize a partial system approach in the parametric domain. In general, parametrization-based methods suffer from the additional distortion of the map and the need to stitch parameterized charts for high genus surfaces. Valette et al. [START_REF] Valette | Generic remeshing of 3D triangular meshes with metric-dependent discrete Voronoi diagrams[END_REF] perform a discrete version of CVT directly on the input surface. However, the resulting mesh quality can be poor due to the inexact computation. Yan et al. [START_REF] Yan | Isotropic remeshing with fast and exact computation of restricted Voronoi diagram[END_REF], [START_REF] Yan | Low-resolution remeshing using the localized restricted Voronoi diagram[END_REF], [START_REF] Yan | Non-obtuse remeshing with centroidal Voronoi tessellation[END_REF] avoid the parameterization by computing the 3D CVT restricted to the surface. Additionally, they proposed blue-noise remeshing techniques using adaptive maximal Poisson-disk sampling [START_REF] Yan | Gap processing for adaptive maximal Poisson-disk sampling[END_REF], [START_REF] Guo | Efficient maximal Poisson-disk sampling and remeshing on surfaces[END_REF], farthest point optimization [START_REF] Yan | Blue-noise remeshing with farthest point optimization[END_REF], and push-pull operations [START_REF] Ahmed | A simple push-pull algorithm for blue-noise sampling[END_REF], which improve the element quality as well as introducing blue-noise properties. However, these approaches still suffer from common limitations, e.g., geometric fidelity and the minimal angle cannot be explicitly bounded. Moreover, sharp features must be specified in advance.

Another way to avoid the stitching problem is to operate directly on the surface mesh [START_REF] Frey | Mesh generation, application to finite elements[END_REF], [START_REF] Narain | Adaptive anisotropic remeshing for cloth simulation[END_REF]. An efficient isotropic approach proposed by Botsch and Kobbelt [START_REF] Botsch | A remeshing approach to multiresolution modeling[END_REF] takes an edge length as input and repeatedly splits long edges, collapses short edges, equalizes vertex valences and relocates vertex positions until all edges are approximately of the specified target edge length. To extend this work to an adaptive version, Dunyach et al. [START_REF] Dunyach | Adaptive remeshing for real-time mesh deformation[END_REF] replace the constant target edge length with an adaptive sizing field that is sensitive to local curvatures. Since this kind of methods requires neither surface parameterization nor density functions, they are easy to implement, robust for high genus inputs, and efficient for real-time applications. Our method falls into this category. However, we enrich the local operators and apply them in a more selective manner in order to obtain guarantees on the geometric fidelity as well as higherquality results and implicit feature preservation. Error-driven remeshing amounts to generating a mesh that optimizes the tradeoff between geometric fidelity and mesh complexity. Cohen-Steiner et al. [START_REF] Cohen-Steiner | Variational shape approximation[END_REF] propose an error-driven clustering method to coarsen the input mesh. They formulate the approximation problem as a variational geometric partitioning problem, and optimize a set of planes iteratively using Lloyd's iteration [START_REF] Lloyd | Least squares quantization in PCM[END_REF] to minimize a predefined approximation error.

The mesh simplification techniques are similar to errordriven remeshing in some way. Garland and Heckbert [START_REF] Garland | Surface simplification using quadric error metrics[END_REF] use iterative contractions of vertex pairs to simplify models and maintain surface approximation error based on quadric error metrics. Borouchaki and Frey [START_REF] Borouchaki | Simplification of surface mesh using Hausdorff envelope[END_REF] define a fidelity metric named Hausdorff envelope, and simplify and optimize the reference mesh that stays inside the tolerance volume. While they consider the geometric fidelity and element quality simultaneously, nothing is done to improve the worst element quality. Our method guarantees the worst element quality by improving the minimal angle and keeping the mesh complexity as low as possible, with respect to a given error-bound. Based on the concept of tolerance vol-ume, Mandad et al. [START_REF] Mandad | Isotopic approximation within a tolerance volume[END_REF] propose an isotopic approximation method. Their algorithm generates a surface triangle mesh guaranteed to be within a given tolerance volume. However, they generate the mesh from scratch, and mainly focuse on mesh complexity rather than element quality. Instead, we strive for good balances between mesh complexity and element quality with a given input. Feature preservation is crucial in remeshing. However, automatically identifying sharp features on a surface mesh is a difficult problem that depends both on the local shape and on the global context and semantic information of the model. This makes feature detection a strongly ill-posed problem. A wide range of approaches address this problem [START_REF] Sun | Triangle mesh-based edge detection and its application to surface segmentation and adaptive surface smoothing[END_REF]- [START_REF] Chen | Bilateral blue noise sampling[END_REF]. However, none of them works reliably for all kinds of meshes. Most remeshing algorithms avoid this problem by assuming that the features have been specified in advance [START_REF] Yan | Isotropic remeshing with fast and exact computation of restricted Voronoi diagram[END_REF], [START_REF] Yan | Gap processing for adaptive maximal Poisson-disk sampling[END_REF], [START_REF] Yan | Low-resolution remeshing using the localized restricted Voronoi diagram[END_REF], [START_REF] Guo | Efficient maximal Poisson-disk sampling and remeshing on surfaces[END_REF], [START_REF] Yan | Blue-noise remeshing with farthest point optimization[END_REF], [START_REF] Zhong | Particle-based anisotropic surface meshing[END_REF]. Some remeshing techniques try to preserve features implicitly [START_REF] Kobbelt | Feature sensitive mesh processing[END_REF], [START_REF] Lévy | Lp centroidal Voronoi tessellation and its applications[END_REF]. Vorsatz et al. [START_REF] Vorsatz | Feature sensitive remeshing[END_REF] first apply a relaxation in the parameter domain, and then snap the vertices to feature edges and corners. Since they separate remeshing and feature-snapping, the resulting mesh quality near sharp features might be poor. Valette [START_REF] Valette | Generic remeshing of 3D triangular meshes with metric-dependent discrete Voronoi diagrams[END_REF] alleviates this issue by embedding the Quadric Error Metric approximation (QEM) criterion inside the CVT optimization. However, the performance of their clusterbased method is highly dependent on the quality of the input, and the sharp features might not be well preserved when users specify a small vertex budget. Jakob et al. [START_REF] Jakob | Instant field-aligned meshes[END_REF] propose a general framework for isotropic triangular/quaddominant remeshing using a unified local smoothing operator, in which the edges naturally align to sharp features. However, little attention is paid on the approximation error and element quality. We address this problem by optimizing the element quality explicitly in combination with implicit feature preservation based on the new defined feature intensity functions.

ALGORITHM OVERVIEW

Given a 2-manifold triangular mesh M I , the goal of our algorithm consists in finding an improved surface mesh M R with approximation error below δ, minimal interior angle above θ and mesh complexity below N . The main idea is to transform the mesh by a series of discrete local operators as illustrated below: In our algorithm, only edge collapse, edge split and vertex relocation operators are employed. Edge flips are implicitly performed as a combination of an edge split followed by an edge collapse (cf. Fig. 2). This convention not only lowers the combinatorial complexity of all operators but is also advantageous for our implicit feature preservation since edge flips tend to destroy sharp creases and would require additional nontrivial checks.

The remeshing algorithm is designed to perform local operators in a conservative manner. More specifically, a local operator is only executed if it respects the approximation error bound δ, does not introduce new interior angles below current minimal angle θ min and maintains the 2manifoldness of the mesh. This behavior can be interpreted as a conservative greedy algorithm, which in each step identifies the most promising operation that improves the result (either coarsens or improves angles), while never leaving the feasible set of meshes with approximation error below δ. Note that since pure topological edge collapses and edge splits improve the mesh quality very rarely, these operations are always combined with vertex position optimization.

The most crucial design choices of the algorithm are the scheduling of different operators and the efficient modeling and handling of approximation error queries (cf. Sec.4). The algorithm passes through three different stages:

The initial phase concentrates on coarsening and runs a mesh simplification, which solely performs edge collapses.

The second phase then tries to lift the minimal interior angle above the user-provided bound θ. For this task we devised three different processes, which are tried subsequently, as illustrated in Fig. 3. Since an edge collapse reduces the complexity, it would be the best way to improve the minimal angle. If edge collapse is not possible, we try to improve the minimal angle by relocating one of the triangle's vertices. If vertex relocation also fails, edge splits are considered as the last option as they increase the mesh complexity. While edge splits do not directly improve the minimal angle, they are crucial to enrich the local mesh connectivity in order to enable improvements in subsequent steps. We apply an approach similar to the longest-side propagation path [START_REF] Rivara | New mathematical tools and techniques for the refinement and/or improvement of unstructured triangulations[END_REF], described in more detail in Sec. 3.2.

The third stage of our algorithm freezes the mesh connectivity and concentrates on global vertex relocation. This stage is designed to improve the average quality of the mesh elements, while not violating the user-specified bounds. The pseudocode of our remeshing algorithm is shown in Alg. 1. update Q and #V ; 9: end while 10: F inalV ertexRelocation(M I , M R , δ, θ);

In the following we provide more details on the three phases of initial mesh simplification, greedy angle improvement and vertex relocation.

Initial Mesh Simplification

The input mesh is often densely sampled such that the mesh complexity can be significantly reduced without violating the approximation error bound. While such a mesh simplification could potentially also be done on the fly, a specialized pre-process turns out to be more efficient. The goal of this step consists in finding a mesh, which is significantly coarser and offers a good starting point for the second phase of element quality improvement. Simplification is achieved through iteratively collapsing edges and relocating related vertices as long as the approximation error does not exceed δ. Short edges, or those opposite to small angles are likely to lower the mesh quality and are consequently collapsed first. This is achieved through using a priority queue sorted by increasing values of the following function:

l(h) • (θ 1 (h) + θ 2 (h))/2
, where l(h) denotes the length of halfedge h and θ i (h) denote the two angles opposite to h and h's opposite halfedge. The pseudocode of the corresponding function is in Alg. 2.

Alg. 2 InitialM eshSimplif ication(M I , M R , δ)

1: fill Q with all halfedges;{Q is a dynamic priority queue} 2: while Q = ∅ do 3:

h ← P op(Q); 4: if CollapseAndRelocateImproves(h, δ, M I , M R ) then 5: CollapseAndRelocate(h, δ, M I , M R ); 6: end if 7:
update Q; 8: end while

Greedy Improvement of Angles

The second phase is designed to improve the mesh quality by iteratively increasing the smallest angle in the mesh, until the desired angle bound θ is satisfied or the complexity limit N is reached. Our approach repeats the following process: We simulate a potential operation, test whether the resulting mesh improves and only in this case perform the candidate operation. The mesh improvement test additionally contains several important validity checks.

Mesh Improvement Test: we simulate each potential operation and measure if the following constraints are satisfied:

• Topology. For edge collapses topology changes are prevented by checking the link condition [START_REF] Edelsbrunner | Geometry and topology for mesh generation[END_REF].

• Geometry. The operator should not create fold-overs by flipping the orientation facets (cf. Fig. 4).

• Fidelity. The approximation error between M R and M I should remain below δ (cf. Fig. 6).

(a) Edge collapse.

(b) Vertex relocate. (b) Vertex P will be relocated. ST P creates fold-overs if P before and P after relocation are located on the two sides of line ST . In 3D, we say that a triangle creates fold-overs if its normals before and after applying a local operator have opposite orientations.

In order to improve the minimal angle, we greedily apply operators that pass the mesh improvement test. The operators are tested in the following order: edge collapses, vertex relocations and edge splits. The pseudocode of the greedy improvement is shown in Alg. 3.

Alg. 3 GreedyImproveAngle(θ min , M I , M R , δ, θ) Input: θ min {Angle requiring improvement} M I {Input mesh} δ > 0 {Error-bound threshold} θ ≥ 0 {Minimal angle threshold} Output: M R {Locally improved mesh} 1: h ← halfedge opposite to θ min ; 2: if CollapseAndRelocateImproves(h, θ min , θ, δ, M I , M R ) then 3: CollapseAndRelocate(h, θ min , θ, δ, M I , M R ); 4:
return; 5: end if 6: let v o , v s and v e be h's opposite, start and end vertex;

7: for v ∈ {v o , v s , v e } do 8: if RelocateImproves(v, θ min , θ, δ, M I , M R ) then 9:
Relocate(v, θ min , θ, δ, M I , M R ); SplitAndRelocate(h l , θ, δ, M I , M R ); 16: end if Longest-Side Propagation Path: If neither edge collapse nor vertex relocation is possible, we search the longest edge h l in the neighborhood [START_REF] Rivara | New mathematical tools and techniques for the refinement and/or improvement of unstructured triangulations[END_REF] and then split it. Starting from an edge e, we iteratively move on to the longest edge of the neighboring two triangles until no further enlargement is possible or we hit the boundary. Though the edge split operator does not increase θ min directly, it modifies the local connections such that θ min can be improved in later iterations (cf. Fig. 3(c)). While there might exist other ways for local connectivity modification, experiments show that this strategy works well in our algorithm.

Final vertex relocation

To achieve a better overall element quality, in the third phase we perform a series of vertex relocations until the angles can no longer be significantly improved. In contrast to the second phase which only focuses on specific regions where angle improvement is required, the third stage optimizes all vertex locations. We maintain the relocation candidate set in a queue, which is initialized with all vertices. Whenever a vertex is relocated, all its neighbors are added to the queue, since a change in the neighborhood might enable further improvements. The pseudocode of the corresponding function is shown in Alg. 4.

Alg. 4 F inalV ertexRelocation(M

I , M R , δ, θ) 1: fill Q with all vertices of M R ;{Q is a FIFO queue} 2: while Q = ∅ do 3: v ← P op(Q); 4: if RelocateImproves(v, θ min , θ, δ, M I , M R ) ≥ ∆θ then 5: Relocate(v, θ min , θ, δ, M I , M R ); 6: add neighbors of v to Q 7:
end if 8: end while Optimal vertex positions would be found if no vertex could be relocated to a better position anymore. However, we empirically restrict the angle improvement ∆θ = 0.1 • (cf. Sec. 6.4), since afterward there are usually no significant improvements anymore.

ERROR METRIC

Hausdorff Distance

We use the Hausdorff distance to measure the approximation error between M R and M I . Let d(p, q) denote the Euclidean distance between two points p and q in 3D space. The distance of a point p to a surface M is defined as the shortest distance between p and any point of

M d(p, M ) = min q∈M d(p, q).
(

) 1
The one-sided Hausdorff distance from a source surface M to a target surface T is defined as the maximum of all such point to surface distances :

d h (M, T ) = max p∈M d(p, T ). (2)
The one-sided Hausdorff distance is in general not symmetric, i.e.

d h (M, T ) = d h (T, M ). It is easily possible to construct counter-intuitive cases where d h (M, T ) = 0 but d h (T, M ) is arbitrarily large 1 .
The two-sided Hausdorff distance [START_REF] Henrikson | Completeness and total boundedness of the hausdorff metric[END_REF] between M and T resolves this issue by symmetrization

d H (M, T ) = max{d h (M, T ), d h (T, M )}. (3)

Approximating d H with Stratified Sampling

The exact evaluation of the two-sided Hausdorff distance is computationally very expensive [START_REF] Tang | Interactive Hausdorff distance computation for general polygonal models[END_REF]. However, by careful surface sampling in combination with local updates of shortest point-to-surface links it is possible to obtain an efficient yet sufficiently accurate approximation, as discussed next.

Assume that M is sampled by a point set S M ⊂ M . Then the one-sided Hausdorff distance can be approximated by

d h (M, T ) ≈ max a∈S M d(a, T ). (4)
By additionally sampling T we obtain an approximation of the two-sided Hausdorff distance

d H (M, T ) ≈ max max a∈S M d(a, T ), max b∈S T d(b, M ) , (5)
with S T being a set of point samples on T . Note that our approximation still measures the exact distance from sample points to the complete surface, which provides significantly higher accuracy than a point cloud distance d H (S M , S T ).

Moreover, it ensures that we strictly underestimate the real distance. Following the triangle inequality, the approximation error of our sampled Hausdorff distance is bounded by max{d h (M, S M ), d h (T, S T )}, i.e. the maximum gap between sample points. Consequently, in order to guarantee a good approximation, we target a uniform sampling of the surfaces. However, for piecewise linear triangle meshes the maximal distance often occurs at creases or corners, i.e. at mesh edges or vertices. This two observations motivate our stratified sampling approach, which is uniform on faces but additionally adds samples on edges and vertices, as a kind of greedy importance sampling.

Sampling the facets. Instead of uniformly sampling the complete surface, we empirically found that better results can be obtained by uniformly sampling per triangle. In this way we obtain an automatic adaptivity, since more samples are available in structurally complex and highly curved areas, where also more triangles are necessary to describe such shape. Since our meshes are often highly non-uniform, we add a slight local smoothing of the sampling density by

n(f i ) = n f • 1 + |N fi | 1 + fj ∈N f i Aj Ai , (6)
where n f is the average number of samples per facet specified by the user, N fi are the neighbor facets that share vertices with f i , and A i is the area of facet f i . Choosing a large 1. Choose T as a sphere of radius r and M as a hemisphere subset. Then with r → ∞ also

d h (T, M ) → ∞.
This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.

The final version of record is available at http://dx.doi.org/10.1109/TVCG.2016.2632720

Copyright (c) 2016 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

n f offers a tighter approximation of the Hausdorff distance, however, resulting in high computational complexity. In our experiments, we found that n f = 10 gives a good tradeoff between efficiency and effectiveness (cf. Fig. 12).

To distribute the samples evenly on a triangle f i , we first generate n(f i ) samples randomly, and then perform a Lloyd relaxation process onto a bounded Voronoi Diagram (BVD) [START_REF] Tournois | 2D Centroidal Voronoi Tessellations with Constraints[END_REF] By setting S I and S R as the stratified samples on M I and M R , the Hausdorff distance between M I and M R is approximated using Eq. ( 5), as illustrated by Fig. 6.

Local Update Scheme

Each operator of our remeshing algorithm only changes a local area L ⊂ M R of the target mesh M R (cf. Fig. 3 green area). Hence, it is possible to rapidly compute the Hausdorff distance for the updated mesh. In general both directions of the Hausdorff distance change and require an update. Identifying the correct partitioning of samples S I would in general require global tests and is thus time consuming. Therefore, in order to enable a rapid update we only approximate this decomposition by tracking the history of samples in S I . The key idea is to store for each sample s j ∈ S I a link to its closest triangle t j ∈ M R (Fig. 6). After a local modification of M R these links will typically change in the vicinity of the modified area L. Thus, in order to avoid global checks, we only update the links of samples connecting to a region L + , which enlarges L by an additional ring of neighboring triangles (cf. Fig. 3 gray area). These new links can be efficiently found by constructing and querying an the axis aligned bounding box tree of the local area L + . Note that the resulting approximated Hausdorff distance is a strict overestimator because shorter links to L+ might exist that are not investigated by our approximation. A positive side effect of our localized update scheme is that links cannot jump between geodesically far away regions, which could cause topologically inconsistent links.

IMPLICIT FEATURE PRESERVATION

A clean representation of geometric features as for instance sharp creases is important in many applications ranging from visualization to simulation. In a polygonal mesh, we distinguish between vertex and edge features. Vertex features include tips, darts, cusps and corners, while edge features are either creases or boundaries. Typical examples are depicted in Fig. 7.

Instead of requiring an explicit tagging of features as most remeshing approaches, e.g. [START_REF] Yan | Isotropic remeshing with fast and exact computation of restricted Voronoi diagram[END_REF], [START_REF] Yan | Gap processing for adaptive maximal Poisson-disk sampling[END_REF], [START_REF] Yan | Low-resolution remeshing using the localized restricted Voronoi diagram[END_REF], [START_REF] Yan | Blue-noise remeshing with farthest point optimization[END_REF], [START_REF] Zhong | Particle-based anisotropic surface meshing[END_REF], our goal is to implicitly preserve features. This not only releases users from the time-consuming burden of manual tagging but moreover often enables the recovery of features that were lost, e.g. through inappropriate sampling by a 3D laser scanner.

In principle feature preservation could be simply a byproduct of approximation error minimization since incorrectly meshing features usually induces a large approximation error. Nevertheless, there are several reasons why special care is still required. First of all, the minimization of the approximation error is a non-convex problem such that bad initializations might lead to low-quality local minima, not well representing features. Moreover, anticipating feature locations and placing vertices accordingly will speed up the overall method. In our method special feature handling is done whenever vertices are either newly placed or relocated, i.e. during (i) edge collapse, (ii) edge split and (iii) vertex relocation. Since robust and automatic feature detection is a difficult yet unsolved problem, we rely on a softer identification of features by means of a feature intensity function defined at vertices of the mesh.

Feature Intensity

Feature vertices can be characterized by large Gaussian curvature K(v), which in the discrete setting is identical to the angle defect, i.e.

K(v) = π -θ sum (v) v is on a boundary, 2π -θ sum (v) otherwise, (7)
where θ sum (v) is the sum of interior angles adjacent to v.

Feature edges are characterized by large dihedral angles. Accordingly, for a vertex we define the feature edge intensity E(v) to be the maximal unsigned dihedral angle of an edge adjacent to v

E(v) = max e∈Ne(v) |D(e)| (8)
where N e (v) are the edges adjacent to v and D(e) is the dihedral angle at e.

Finally, the feature intensity F(v) is defined as the combination

F(v) = (τ (|K(v)|) + 1) • (τ (E(v)) + 1) -1
with the transfer function τ (x) = min{π, 2 • x}, which rescales values by a factor of 2 and clamps them at π. Thus, the feature intensity is a value between 0 and ((π + 1) 2 -1) and corresponds to a logical or whenever one of the individual intensities vanishes. Fig. 8 shows an examples of the three fields |K|, E and F.

|K(v)| E(v) F(v)

Fig. 8. From left to right Gaussian curvature |K|, feature edge intensity E and combined feature intensity F , with higher intensities in red.

Feature-Sensitive Vertex Relocation

In our remeshing algorithm each local operation is combined with a subsequent relocation of the modified vertex in order to minimize the approximation error. This relocation of a single vertex is done in two stages. First a feature-sensitive initialization, specifically adapted to the local operation, and then a nonlinear minimization of the two-sided Hausdorff distance. A careful initialization of vertex positions is important to avoid poor local minima of the non-convex Hausdorff energy and additionally increases performance.

Edge Collapse Initialization: During an edge collapse of edge e ij the two vertices v i and v j merge into a new vertex v m and a new position for v m must be specified. One feature-sensitive standard technique uses error quadrics [START_REF] Valette | Generic remeshing of 3D triangular meshes with metric-dependent discrete Voronoi diagrams[END_REF], [START_REF] Garland | Surface simplification using quadric error metrics[END_REF]. However, since we anyway perform a more accurate approximation error-driven relocation subsequently, a simpler and faster position initialization based on our feature intensity is sufficient. We distinguish two cases. If the edge is incident to a single strong feature we want that v m snaps onto this feature, meaning that the default behavior is to snap to the vertex v k with higher feature intensity, i.e. v k = arg max v∈{vi,vj } F(v). An unclear situation arises, when both feature intensities are similar, i.e. F(v i ) ≈ F(v j ). This happens either in regions without features or for edges along a crease, where it is reasonable to initialize v m to the edge midpoint. We decide for the midpoint initialization based on a parameter ω, whenever

|F(v i ) -F(v j )| < ω • max(F(v i ), F(v j ))
with ω = 0.15 in all our examples.

Edge Split Initialization: Since an edge split does not change the geometry we simply always initialize the new point as the edge midpoint. Further improvement is achieved through a nonlinear minimization of the twosided Hausdorff distance, as discussed next.

Vertex Relocation Initialization: Without topological mesh modifications (edge collapse or split), a vertex relocation is initialized in the following way. The goal is to anticipate a good relocation position for v while preserving feature corners and creases. Therefore, we first classify the vertex v as either being (i) a feature vertex, (ii) a crease vertex or (iii) a smooth vertex. If v is a feature vertex it simply remains at its position. If it is on a crease, we relocate v to the CVT barycenter of its two neighboring crease vertices.

And only if it is a smooth vertex we relocate it to the CVT barycenter of all one-ring neighbors. The classification is done by counting how many neighbors v i are of similar or higher importance as v, i. This optimization is highly nonlinear, since changing the position of v changes the samples S R of the modified mesh M R as well as the links from the input mesh samples S I to M R . We perform an optimization similar to the Hausdorff distance minimization technique proposed by Winkler et al. [START_REF] Winkler | Mesh massage[END_REF]. However, we improve the feature-sensitivity by adjusting the weighting scheme with our feature intensity function F, as detailed next.

We want to relocate v in order to minimize the twosided Hausdorff distance between the local area L ⊂ M R consisting of all one-ring triangles of v and the subregion of M I with links into L, referred to as M I→L . Assume that the corresponding subsets of samples for our Hausdorff distance approximation are S L ⊂ S R and S I→L ⊂ S I . Then according to Motzkin and Walsh's theorem [START_REF] Motzkin | Polynomials of best approximation on a real finite point set[END_REF] there exist weights w i and ŵi such that minimizing Eq. ( 5) w.r.t. the center vertex v is equivalent to minimizing

ai∈S L w i |a i (v) -âi | 2 + bi∈S I→L ŵi | bi (v) -b i | 2 (9)
where âi ∈ M I→L and bi (v) ∈ L are the closest points to a i and b i respectively. Freezing the closest point pairs and assuming a linear barycentric relation

a i (v) -â = α i v + β i d i + γ i e i -â = α i v -p i
with constant p i = â -β i d i -γ i e i and similarly expressed bi (v) -b i = αi v -pi , the optimal position v * can be computed analytically

v * = ai w i α i p i + bi ŵi αi pi ai w i α 2 i + bi ŵi α2 i ( 10
)
To find optimal weights w i and ŵi , Winkler et al. [44] use Lawson's algorithm [START_REF] Lawson | Contributions to the Theory of Linear Least Maximum Approximation[END_REF] with iterative updates of the form

w (k+1) i = w (k) i • d(a (k) i , â(k) i ), ( 11
)
where d(a

(k) i , â(k) i )
is the Euclidean distance of the closestpoint pair (a i , âi ) after the k-th iteration, and initialization w (0) i = 1. The idea behind this scheme is that samples with larger distances get a higher weight in the next iteration. Based on our feature intensity function and a sample density estimation, the weight update can be improved to

w (k+1) i = w (k) i • d(a (k) i , â(k) i ) • V(a (k) i ) • F(a (k) i ) ( 12
)
where V(a i ) is the Voronoi cell area of sample a i , and F(a i ) is the linearly interpolated feature intensity value of a i . Fig. 11 illustrates the additional weighting factors. Advantages are a better feature preservation, improved robustness w.r.t. non-uniform sampling and faster convergence.

Each time when a local operator is applied, we iteratively optimize the new position of the effected vertex by minimizing Eq. ( 9) using the new weighting defined in Eq. ( 12). The optimization procedure is similar to the Expectationmaximization (EM) algorithm: in each iteration, we first calculate the optimal position v * of vertex v using Eq. ( 10), and then move v to v + λ(v * -v), λ ∈ (0, 1] with default λ = 0.9 and update the closest point pairs. In practice, a few iterations usually suffice to get very close to the optimum.

V(ai) F(ai) V(ai) • F (ai)

Fig. 11. Sample weights of a cube corner, with the integration of V(a i ) and F (a i ). Higher intensity of red highlights larger weight values.

EXPERIMENT RESULTS AND DISCUSSIONS

We implemented our approach in C++ and tested on a 64bit Windows 8.1 operating system. The CGAL library [START_REF]CGAL, Computational Geometry Algorithms Library[END_REF] provided most of the basic data structures and operations. Timings for all the examples were conducted on a single 3.40GHz Intel(R) Core(TM) i7-2600K CPU with a 16GB RAM. We provide next a complete evaluation of our algorithm and comparisons with state-of-the-art approaches.

Evaluation of the Local Error Update Scheme

For efficient local error update we use the axis aligned bounding box tree of CGAL. As global error update is too compute-intensive we only verify how the inner patch sizes and sampling densities affect the effectiveness and efficiency of the local error update scheme (Fig. 12). We find that even when the inner patch size is set to one-ring (cf. Fig. 3), our local error update approach catches more than 99.9% of the global nearest points. The accuracy increases little with higher sampling density and larger inner patch size. We set the inner patch sizes as one-ring facets and the average sampling number per facet as ten in all experiments. 2) as input.

Evaluation of Vertex Position Optimization

Solely applying the vertex position initialization makes most vertices of M R stay near or on the surface of M I . Though this quasi-interpolation preserves sharp features, the optimal geometric fidelity is usually achieved when M R is an approximation of M I . We visually demonstrate the difference between interpolation and approximation in Fig. 13. In the interpolation case, the Hausdorff distance (Hdist) [START_REF] Cignoni | Metro: measuring error on simplified surfaces[END_REF] between the sphere and the icosahedron is 6.48(%bb), and the root mean square (RMS) distance is 4.84(%bb); while in the approximation case, the Hdist and RMS distances between them are 5.00(%bb) and 1.09(%bb), respectively. However, the approximation might destroy features when minimizing the local sample pair distances [START_REF] Winkler | Mesh massage[END_REF]. Fig. 14 compares the average distance and RMS distance based on Lawson's weighting scheme (Eq. ( 11)) and our improved weighting scheme (Eq. ( 12)). Generally, our weighting scheme reduces the average distance and RMS distance about 2.3% and 3.1% respectively. However, it reduces the approximation error of vertices on sharp features about 11.8% and 12.8%, respectively. Therefore, applying the improved weighting scheme does not only reduce the approximation error, but also better preserves sharp features. Fig. 14. Comparison of geometric fidelity between Lawson's weighing scheme [START_REF] Winkler | Mesh massage[END_REF] and ours(Eq. 12). We sort the vertices in M I according to their saliency function values in ascending order, and compute their average distance and RMS distance in each bin. In this experiment, δ = 0.2(%bb), and θ = 30 • . The above data are the averages of 10 consecutive executions with the Hand model (Fig. 10) as input.

We extensively tested how the iteration count and relocate ratio (λ) affect the approximation error and execution time in the vertex position optimization procedure (in Fig. 15). We found that setting the iteration count to two and the relocate ratio to 0.9 achieved the best compromise between effectiveness and efficiency. This configuration is used in all our later experiments.

Evaluation of Initial Mesh Simplification

In general, applying Alg. 2 further reduces 20% vertices on average. However, the execution time would be 2-3 times slower than not applying it. Fig. 16 shows the remeshing results of the Egea model with and without applying initial mesh simplification, and Tab. 2 further compares the differences (OUR vs. OUR*). Usually, if users care more about the execution time than the mesh complexity, Alg. 2 can be disabled. However, we enable Alg. 2 by default for a better compromise among mesh complexity, element quality and approximation error.

Evaluation of Final Vertex Relocation

To measure how the overall element quality is improved by applying Alg. 4, we introduce two new measurements: the first is the average minimal angles of all triangles in M R and the second is the average value of triangle qualities defined as

Q t = 2 √ 3S t /(p t h t ) [49]
, where S t is the area of triangle t, p t the in-radius of t and h t the length of the longest edge in t. We tested the statistic element quality and the execution time with varying ∆θ in Alg. 4 (Fig. 17), and found when ∆θ < 0.1 • , the quality improvement is not significant. In our experiments, we set the default ∆θ as 0.1 • in Alg. 4. . 2) as input.

Influence of the Minimal Angle Threshold and the Mesh Complexity

Our algorithm produces results with either desired minimal angle threshold or desired mesh complexity, depending on the user specified parameters. In order to get the desired minimal angle, N in Alg. 1 should be set very large; otherwise, θ should be set very large (e.g. 60 • ). Usually, the larger the value θ or N is specified, the better element quality is achieved. We tested the remeshing results of the Fandisk model with θ varying from 0 • to 40 • and N varying from 0.15k to 2.8k, and show the results in Fig. 1. The complete attributes are listed in Tab. 1. Note that in Fig. 1(g), one vertex has been relocated a little away from the crease after refinement, such that the minimal angle is improved. This happens when optimizing a mesh with sharp features up to a high minimal angle threshold.

TABLE 1

Influence of θ and σ. V 567 is the percent of vertices with valences 5, 6 and 7. We set δ = 0.2(%bb), and the Fandisk (Fig. 1(a)) is the input. We compare our results with the state-of-the-art methods, and find only in a small portion of results presented in [START_REF] Yan | Low-resolution remeshing using the localized restricted Voronoi diagram[END_REF] and [START_REF] Yan | Gap processing for adaptive maximal Poisson-disk sampling[END_REF], the minimal angles exceed 35 • (with highest record 38 • ). For most other methods, the minimal angles vary between [25 • , 35 • ]. Contrary to the previous work, our method is able to generate results with minimal angles higher than 35 • in all test cases. The complete comparison with the state-of-the-art methods is shown in Sec. 6.8.

θ #V Q min θmax RMS(%bb) V 567 (%)

Influence of the Approximation Error Threshold

We demonstrate the influence of δ in Fig. 18. The results indicate that δ does not influence the mesh complexity and the execution time significantly with a fixed θ value. However, two interesting phenomena are observed: 1) when θ ≤ 35 • , the larger the θ value is, the lower the mesh complexity we achieve; however, when θ is set to 40 • , both mesh complexity and execution time increase dramatically. This is because when θ is small, the edge collapse operator is preferentially applied, which increases the minimal angles while reduces the mesh complexity. However, when θ is large, more edge split operators are applied to modify local connections. 2) Within a fixed θ, when δ increases, the mesh complexity decreases slightly and smoothly, since the higher δ is set, the more edge collapse operators are triggered. We compare our results with those provided by Yan et al. [START_REF] Yan | Gap processing for adaptive maximal Poisson-disk sampling[END_REF], [START_REF] Yan | Low-resolution remeshing using the localized restricted Voronoi diagram[END_REF], [START_REF] Yan | Blue-noise remeshing with farthest point optimization[END_REF]. The best record of the approximation error is 0.10(%bb), which is generated by [START_REF] Yan | Isotropic remeshing with fast and exact computation of restricted Voronoi diagram[END_REF] with the Homer model as input. However, it can not be explicitly controlled. For other methods, the approximation error is between 0.3-0.5(%bb), and still cannot be strictly bounded. In striking contrast, our algorithm is able to explicitly and strictly control the approximation error, and achieve the approximation error as low as 0.07(%bb) with the same input. More complete comparisons are in Sec. 6.8.

Robustness

Since in our local error update scheme, the closest point pairs of stratified samples are reliably initialized and locally updated, our method is robust to models with complex topology, holes, and intersections of surfaces. For example, the Close spheres model (Fig. 19(c)) is composed of two rounded half spheres that are positioned very close to each other, and the Klein bottle model (Fig. 19(h)) exhibits nonorientable surfaces that are self-intersecting. Our method generates the right results for both of them. By integrating the feature intensity function, our method successfully handles models with and without clear sharp features. The Helmet (Fig. 19(a)), Fertility (Fig. 19(b)) and Close spheres models are smooth, whereas the U-part (Fig. 19(e)) and the Smooth crease models (Fig. 19(f)) have sharp features. In addition, since the feature intensity also captures boundaries, our method is capable of remeshing models with boundaries (Fig. 19(i)) and smooth features (Fig. 19(f)).

Our method requires no surface parameterization, making it suitable for high-genus models (Fig. 19 and Fig. 19(b)) as well as models with multiple components (Fig. 19(g)). Note that adaptivity is created automatically by the approximation error parameter δ if reasonable without requiring an a priori estimation of a density function. The Helmet, Fertility, U-part and Lion head models in Fig. 19 clearly illustrate this advantage.

Our method is suitable to process models with very high/low resolutions and/or badly shaped triangles. To tackle dense models, typically the initial mesh simplification strongly reduces the mesh complexity. For very coarse models, users can optionally increase the sampling density for better error control. We present two typical examples: the input Fertility model has 13k vertices, whereas our remeshing result (Fig. 19(b)) has only 2.9k vertices, thanks to the effectiveness of Alg. 2. The input U-part model possesses only 86 vertices. By sampling 50 points in each facet averagely, we get the result with 347 vertices, and the geometric fidelity is well-controlled.

Comparisons

We compare our approach to the state-of-the-art techniques in terms of efficiency, geometric fidelity, element quality and mesh complexity. For simplicity, only the most efficient methods (RAR [START_REF] Dunyach | Adaptive remeshing for real-time mesh deformation[END_REF] and MMGS, an improvement of YAMS [START_REF] Frey | Mesh generation, application to finite elements[END_REF]) and the methods that produce the best results in Yan and Wonka's conduction [START_REF] Yan | Gap processing for adaptive maximal Poisson-disk sampling[END_REF] (CVT [START_REF] Yan | Isotropic remeshing with fast and exact computation of restricted Voronoi diagram[END_REF] (100 iterations in our experiments), CVD [START_REF] Valette | Generic remeshing of 3D triangular meshes with metric-dependent discrete Voronoi diagrams[END_REF] and MPS [START_REF] Yan | Gap processing for adaptive maximal Poisson-disk sampling[END_REF]) are compared with identical inputs. Among all the compared methods, CVT and CVD require the number of vertices to be specified, and the Hdist is required in RAR and MMGS. To make the results comparable, we set the Hdist of RAR to the same value as that of our method, and carefully adjusted the Hdist parameter for MMGS, such that it generates results with the same complexity as CVT, MPS and CVD. In our method, OUR* means Alg. 2 is disabled while OUR means it is enabled. A detailed comparison is listed in Tab. 2, and Fig. 20 shows a close-up comparison. More visual comparisons are provided in the supplemental materials, available online.

From all compared methods, RAR performs most efficiently and introduces the lowest mesh complexity. However, the geometric fidelity cannot be guaranteed. MMGS is also efficient, yet it introduces much higher Hdist distances. Our method is almost at the same level of efficiency as CVT when θ = 35 • , but is much slower when θ = 40 • .

According to [START_REF] Yan | Gap processing for adaptive maximal Poisson-disk sampling[END_REF], CVT performs best in keeping high geometric fidelity, but cannot explicitly bound the approximation error. By setting δ comparable to CVT's best results, our method consistently produces results with strictly bounded Hdist, and produces lower Hdist and RMS distances in most cases. We also find that our method consistently generates results with higher minimal angle and triangle quality than MPS, due to the fact that our method provides a better improvement of "worst element quality" measured by the minimal angle. However, since our method does not optimize the global connectivity of the input, our results have lower regularity (measured as V 567 ) than MPS.

From Tab. 2, we see that by setting the desired resolution lower than MMGS, CVT, MPS and CVD, we still get results with higher geometric fidelity and better bounds of minimal angle. For the very dense Bunny model, when Alg. 2 is enabled, our method even reduces the complexity of the Comparison with the state-of-the-art methods. For all the method, the input parameters are highlighted with italic fonts, and the best results are highlighted in bold. In column Hdist(%bb), the values before "/" are the input, and the values after "/" are the real HDist measured by Metro [START_REF] Cignoni | Metro: measuring error on simplified surfaces[END_REF].

Input Methods #V Q min θ min ( • ) θmax( • )
Hdist (%bb) RMS (%bb) θ < 30 Bunny model to 5% of the input without violating the errorbound constraint. When it is disabled, the resolutions are still lower than the inputs in most cases, since the edge collapse operator has high priority in Alg. 1.

To the best of our knowledge, RAR, MMGS, CVT and MPS require sharp features to be specified or detected in advance, which may be time-consuming or error-prone. Though CVD is able to preserve features implicitly, it leads to results with lower geometric fidelity and element quality than our method.

Since our method does not optimize the element quality globally, the average element quality is not superior to remeshing based on farthest point optimization [START_REF] Yan | Blue-noise remeshing with farthest point optimization[END_REF] (FPO). However, we consistently produce results with better geometric fidelity and larger minimal angle than FPO.

Limitations

Although practically the minimal angles can be improved to values above 35 • in all our test cases, we do not have any theoretical guarantee for the convergence with a specified θ. To challenge our algorithm, we set θ to its theoretical upper bound and show the best results that our algorithm achieved in Tab. 3. We found the approximation error can still be bounded. However, the algorithm runs into infinite loops or generates degenerated edges while refining.

Another limitation is that we can only tackle 2-manifold meshes, for the reason that our local operators highly rely on the topology information of local regions. Finally, our method cannot remesh noisy models robustly, since the feature intensity function will interpret the noise as some kind of features and thus tries to preserve it.

CONCLUSION

We presented a novel surface remeshing algorithm based on minimal angle improvement. In this framework, the minimal angle of the input model is sequentially improved by applying local operators. Furthermore, an efficient and reliable local error update scheme was designed and embedded for explicitly bounding the approximation error, and two novel feature intensity functions were defined and integrated in the vertex relocation, in order to preserve features implicitly. Compared to the state-of-the-art, our method consistently generates results with higher element quality and lower mesh complexity but satisfied error-bounds. The resulting meshes are well suited for robust numerical simulations, since they offer bounded approximation error and large minimal angles. However, there are still some limitations (Sec. 6.9), which motivate future work in the following aspects: 1) providing a theoretical convergence guarantee with a specified θ; 2) extending the current implementation such that this framework can be applied to triangle soups or even point clouds, for error-bounded and feature preserving mesh reconstruction; and 3) exploring more robust feature intensity functions for measuring sharp features on noisy models.

  (a) Input model. (b) θ = 0 • . (c) θ = 10 • . (d) θ = 20 • . (e) θ = 30 • . (f) θ = 35 • . (g) θ = 40 • .

Fig. 1 .

 1 Fig. 1. Examples of surface meshes generated with our approach. The input model (a) has 7.2k vertices. From (b) to (g) are the results with different minimal angle threshold θ. The angle distributions and their corresponding approximation error are shown below each result, in blue bars and red curves respectively. The error-bound threshold δ is set to 0.2% of the diagonal length of the input's bounding box (%bb).

Fig. 2 .

 2 Fig. 2. Common local operators [1].

Alg. 1 2 :

 12 Remeshing(M I , M R , δ, θ, N ) Input: M I {Input mesh} δ > 0 {Error-bound threshold} θ ≥ 0 {Minimal angle threshold} N > 0 {Desired mesh complexity} Output: M R {Remeshing result} 1: M R ← M I ; InitialM eshSimplif ication(M I , M R , δ); 3: fill Q with angles of M R smaller than θ; {Q is a dynamic priority queue} 4: #V ← the number of vertices in M R ; 5: while Q = ∅ and #V < N do 6: θ min ← P op(Q);

7 :

 7 GreedyImproveAngle(θ min , M I , M R , δ, θ); 8:

  (a) Edge collapse. (b) Vertex relocation. (c) Edge split.

Fig. 3 .

 3 Fig. 3. Local operators. The local patches before/after applying the operators are shown in the first/second rows, respectively. The inner patch (green) contains the facets that will be directly affected by the local operators, and the outer patch (gray) includes the facets who share vertices with facets of the inner patch. In each sub figure, the left shows the inner case while the right shows the boundary case, in which the boundary edges and vertices are depicted in blue. We depict the current minimal angle θ min in yellow.

Fig. 4 .

 4 Fig. 4. Creation of fold-overs in 2D. (a) Edge P Q will be collapsed into Q. ST Q creates fold-overs if P and Q lie on the two sides of line ST ;(b) Vertex P will be relocated. ST P creates fold-overs if P before and P after relocation are located on the two sides of line ST . In 3D, we say that a triangle creates fold-overs if its normals before and after applying a local operator have opposite orientations.

Fig. 5 .

 5 Fig. 5. Stratified sampling process. (a) Initial facet samples; (b) Optimized facet samples with BVD; (c) Edge samples. The facet, edge, and vertex samples are rendered in green, yellow and red, respectively. Sampling the edges. By counting the number of incident Voronoi cells to an edge we first estimate the local sampling density. The resulting number of samples is then evenly distributed along the edge (Fig 5(c)). Sampling the vertices. The position of a vertex is simply regarded as its own vertex sample.

Fig. 6 .

 6 Fig. 6. Hausdorff distance approximation with stratified samples. Samples on vertices, edges, and faces are rendered in red, yellow, and green respectively. We use the circles centered with crosses to indicate the samples on one mesh, and the circles centered with dots to indicate their closest points on the other mesh. The approximated one-side Hausdorff distances are shown as the blue links, and the overall approximated Hausdorff distance is simply the length of the longer one.

  Updating d h (M R , M I ): First of all notice that the one-sided Hausdorff distance can be decomposed w.r.t. to a modified local areaL ⊂ M R into d h (M R , M I ) = max{d h (L, M I ), d h ( L, M I )} with M R = L ∪ L being a partition of M R . Since L is unchanged, checking d h (L, M I ) ≤ δis sufficient to verify that the modification does not violate the approximation error bound δ. Thus, we first re-sample the modified local area L with stratified samples S L , and then efficiently evaluate the approximated Hausdorff distance d h (S L , M I ) with a pre-computed axis aligned bounding box tree[START_REF] Alliez | 3D fast intersection and distance computation[END_REF] for the static mesh M I .Updating d h (M I , M R ): Checking the opposite direction is more difficult for two reasons. Firstly, since M R changes, we cannot simply pre-compute a static search tree. Secondly, decomposing the Hausdorff distance in the second argument is more intricate. In addition to M R , M I must also be decomposed correctly. More specifically, we haved h (M I , M R ) ≈ d h (S I , M R ) = max{d h (S L I , L), d h (S L I , L)}where S I = S L I ∪ S L I is a partitioning of samples on M I into those closer to L and L (the rest of M R ) respectively.

Fig. 7 .

 7 Fig. 7. Vertex features: (a) Tips, (b) Darts, (c) Cusps and (d) Corners are feature vertices who are adjacent to zero, one, two and three sharp creases respectively.

Fig. 9 .

 9 Fig. 9. Visualization of F (v) in (a). The color of v's Voronoi cell represents the number of important neighbors, with red=0, blue=2, purple=3 and gray=degree(v) is shown in (b) and (c). A wrong classification based on solely feature intensity in (b) is corrected by additionally checking importance of edges in (c).The parameter ζ controls the feature classification. For higher values of ζ more vertices are implicitly treated as features and thus prevented from free movement. Fig.10shows an example classification for ζ varying between 0.3 and 0.7. In our experiments the default is ζ = 0.5

Fig. 10 .

 10 Fig. 10. Influence of the classification tolerance ζ. The Voronoi cells of vertices color code the classification with red=feature vertex, blue=crease vertex and white=smooth vertex. We show the Anchor model with obvious features (top) and the Hand model without obvious features (bottom).

Fig. 12 .

 12 Fig. 12. The nearest point searching accuracy and execution time of our local error scheme, with respect to different sampling densities and inner patch sizes. In this experiment, δ = 0.2% of the bounding box's diagonal length(%bb) and θ = 30 • . The above data are the averages of 10 consecutive executions with the Rockerarm model (Table2) as input.

Fig. 13 .

 13 Fig.[START_REF] Du | Centroidal Voronoi tessellations: applications and algorithms[END_REF]. Demonstration of interpolation and approximation. We use an icosahedron (a) to interpolate and approximate a sphere ((b) and (c)). In interpolation, the vertices of M R are guaranteed to be on the surface of M R , while in approximation, the error is minimized, regardless of whether the vertices of M R are on the surface of M I .

Fig. 15 .

 15 Fig. 15. Effectiveness of the iteration count and relocate ratio. For separability, RMS and Hdist bars are plotted in different scales. We set δ = 0.2(%bb), and θ = 30 • . The above data are the averages of 10 consecutive executions with the Homer model (Table.2) as input.

  Fig. 15. Effectiveness of the iteration count and relocate ratio. For separability, RMS and Hdist bars are plotted in different scales. We set δ = 0.2(%bb), and θ = 30 • . The above data are the averages of 10 consecutive executions with the Homer model (Table.2) as input.

Fig. 16 .

 16 Fig. 16. Egea models with/without executing Alg. 2. we set θ = 40 • and δ = 0.20(%bb). The input (a) has 8.3k vertices. The result with Alg. 2 enabled has 4.2k vertices (b), and spends 330 seconds; the result with Alg. 2 disabled has 7.6k vertices (c), and spends 133 seconds.

Fig. 17

 17 Fig. 17. The effectiveness of ∆θ in Alg. 4. In each sub figure, the first bar indicates the value when no final vertex relocation is applied. In this experiment, δ = and θ = 30 • . The above data are the averages of 10 consecutive executions with the Helmet model (Table.2) as input.

  Fig. 17. The effectiveness of ∆θ in Alg. 4. In each sub figure, the first bar indicates the value when no final vertex relocation is applied. In this experiment, δ = and θ = 30 • . The above data are the averages of 10 consecutive executions with the Helmet model (Table.2) as input.

Fig. 18 .

 18 Fig. 18. Influence of the error-bound threshold δ. The above data are the averages of 10 consecutive executions with the Elephant model (Table2) as input.

Fig. 19 .

 19 Fig. 19. Selected results. In this experiment, δ is set to 0.2%(bb) and θ is set to 35 • . In (h) and (i), the boundaries are rendered in red. The dark part of (h) means the triangle normals are inside.

Fig. 20 .

 20 Fig. 20. A close-up comparison of results with state-of-the-art approaches. In our method, δ is set to 0.20(%bb) and θ is set to 40 • .

  This indicates the execution time of Alg. 2.

							Time
	0	0.15k	0.040	174.4	0.043	76.5	1:42*
	10 0.20k	0.188	155.1	0.038	71.3	2:18
	20 0.26k	0.340	135.5	0.036	78.6	2:41
	30 0.40k	0.482	117.7	0.031	81.7	3:01
	35 0.73k	0.552	109.3	0.029	85.5	3:44
	40	2.8k	0.640	98.8	0.022	98.5	6:14

*

TABLE 2

 2

TABLE 3

 3 Cases when too high θ is specified. We set θ = 60 • , run Alg. 1 until it fails, and record the best results it achieved. Alg. 2 is enabled here.

	Input	#V	θ min	Hdist(%bb)	Fail Types
	Rockerarm	4.2k	41.22	0.20/0.19	Degenerated edges
	Homer	6.1k	41.13	0.20/0.14	Infinite loops
	Triceratops 4.8k	40.02	0.20/0.19	Degenerated edges
	Elephant	4.8k	41.23	0.20/0.13	Infinite loops
	Bunny	2.4k	41.42	0.20/0.19	Infinite loops

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.The final version of record is available at http://dx.doi.org/10.1109/TVCG.2016.2632720Copyright (c) 2016 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

ACKNOWLEDGMENTS

We wish to thank Mario Botsch for providing the source code and results of their method, and Pascal Frey for providing the MMGS software. This work was partially supported by the European Research Council (ERC Starting Grant Robust Geometry Processing #257474), the National Science Foundation of China (61373071, 61372168, 61620106003) and the German Research Foundation (DFG, grant GSC 111, Aachen Institute for Advanced Study in Computational Engineering Science).


Christophe Chalons

Régis Duvigneau
email: regis.duvigneau@inria.fr

Camilla Fiorini

Sensitivity analysis for the Euler equations in Lagrangian coordinates

Keywords: Sensitivity analysis, p-system, Dirac delta function, Finite Volume. MSC (2010): 35L60, 65M08, 49Q12

Sensitivity analysis (SA) is the study of how changes in the inputs of a model affect the outputs. SA has many applications, among which uncertainty quantification, quick evaluation of close solutions, and optimization. Standard SA techniques for PDEs, such as the continuous sensitivity equation method, call for the differentiation of the state variable. However, if the governing equations are hyperbolic PDEs, the state can be discontinuous and this generates Dirac delta functions in the sensitivity. The aim of this work is to define and approximate numerically a system of sensitivity equations which is valid also when the state is discontinuous: to do that, one can define a correction term to be added to the sensitivity equations starting from the Rankine-Hugoniot conditions, which govern the state across a shock. We show how this procedure can be applied to the Euler barotropic system with different finite volumes methods.

Introduction

Sensitivity analysis (SA) concerns the quantification of changes in Partial Differential Equations (PDEs) solution due to perturbations in the model input. It has been a

Christophe Chalons and Camilla Fiorini Laboratoire de Mathématiques de Versailles, UVSQ, CNRS, Université Paris-Saclay, 45 avenue des États-Unis, 78035 Versailles, France e-mail: christophe.chalons@uvsq.fr,camilla.fiorini@uvsq.fr topic of active research for the last years, due to its many applications, for instance in uncertainty quantification, quick evaluation of close solutions [START_REF] Duvigneau | A sensitivity equation method for fast evaluation of nearby flows and uncertainty analysis for shape parameters[END_REF], and optimization [START_REF] Borggaard | A PDE sensitivity equation method for optimal aerodynamic design[END_REF], to name but a few. Note that SA approaches differ from adjoint methods, which are restricted to the evaluation of functional derivatives. Standard SA methods work only under certain hypotheses of regularity of the solution U [START_REF] Bardos | A formalism for the differentiation of conservation laws[END_REF]. These assumptions are not verified in the case of hyperbolic systems of the general form

∂ t U + ∂ x F(U) = 0, x ∈ R, t > 0, U(x, 0) = U 0 (x),
due to possible discontinuities, which can occur even when the initial condition is smooth. If the state U is discontinuous, Dirac delta functions will appear in the sensitivity U a = ∂ a U. Here and throughout this work, a denotes the input parameter of the model which may vary and induce a non trivial sensitivity U a of the state solution U.

In this work, we consider the Euler equations in Lagrangian coordinates in a barotropic case, i.e. the p-system:

∂ t τ -∂ x u = 0, ∂ t u + ∂ x p(τ) = 0, (1)
where τ > 0 is the covolume (i.e. τ = 1 ρ , where ρ is the density of the fluid), u is the velocity and the pressure p(τ) is a function of τ such that p (τ) < 0 and p (τ) > 0. Under these assumptions, (1) is strictly hyperbolic with eigenvalues λ ± = ±c where c = -p (τ) is the Lagrangian sound speed. In this work we will consider p(τ) = τ -γ , where γ = 1.4 is the heat capacity ratio. By differentiating the system (1) with respect to the parameter of interest a and considering smooth solutions of (1), we obtain the sensitivity equations:

∂ t τ a -∂ x u a = 0, ∂ t u a + ∂ x (p (τ)τ a ) = 0. (2)
One can define the following vectors:

U = τ u , F(U) = -u p(τ) , U a = τ a u a , F a (U, U a ) = -u a p (τ)τ a ,
and rewrite the systems (1) and ( 2) in a vectorial form:

∂ t U + ∂ x F(U) = 0, ∂ t U a + ∂ x F a (U, U a ) = 0. (3)
At this stage, it is easy but important to remark that the global system (3) admits the same real eigenvalues λ ± (both with multiplicity 2) as the original system (1) but it is only weakly hyperbolic as soon as τ a = 0. We recall that weak hyperbol-icity means that the Jacobian matrix of the system admits real eigenvalues but is not R-diagonalizable. As a consequence and without any modification of (3), discontinuous weak solutions of the state variable U will generally induce Dirac delta functions in the sensitivity variable U a , in addition to the usual discontinuity, so that the solutions of (3) have to be understood in the sense of measures. However, and as already stated above, we are not interested in considering Dirac delta functions. Instead, we would like to introduce a modification in the sensitivity equations in order to make the system (3) valid in the usual sense of weak solutions also for discontinuous state variables (as done in [START_REF] Guinot | Upwind finite volume solution of sensitivity equations for hyperbolic systems of conservation laws with discontinuous solutions[END_REF]). This is the aim of the next section.

Source term

In order to remove the Dirac delta functions that are naturally present in the solutions of (3), we suggest to add to (2) a source term S, which is of the following form:

S = N s ∑ k=1 δ k ρ ρ ρ k , (4)
where N s is to be associated to the number of discontinuities in the state solution U, ρ ρ ρ k is the amplitude of the k-th correction (to be computed), and δ k is the Dirac delta function δ k = δ (xx s,k ), where x s,k is the position of the k-th discontinuity. Let us consider then a control volume (x 1 , x 2 ) × (t 1 ,t 2 ) containing a single discontinuity indexed by k and propagating at speed σ k . By integrating the equations ( 2) with the additional source term (4) over the control volume, when the size of the control volume tends to zero one has:

ρ ρ ρ k (t) = (U - a -U + a )σ k + F + a -F - a , (5)
where the plus (respectively minus) indicates the value of the sensitivity U a and of the flux F a to the right (respectively left) of the discontinuity. In other words, [START_REF] Guinot | Upwind finite volume solution of sensitivity equations for hyperbolic systems of conservation laws with discontinuous solutions[END_REF] gives a natural meaning of ρ ρ ρ k in terms of a defect measure of the Rankine-Hugoniot relations applied to [START_REF] Borggaard | A PDE sensitivity equation method for optimal aerodynamic design[END_REF]. It is now a matter of defining ρ ρ ρ k in such a way that the new model including the source term is also valid for discontinuities of the state variable (recall that (2) was obtained by differentiating with respect to a the smooth solutions of ( 1)). Considering the Rankine-Hugoniot conditions across a discontinuity of the state variable, namely (U --U + )σ k = F --F + , we differentiate with respect to the parameter a to obtain

(U - a -U + a )σ k + (U --U + )σ k,a = F - a -F + a , with σ k,a = ∂ a σ k .
Comparing the latter with (5), one is thus led to set

ρ ρ ρ k (t) = σ k,a (U + -U -). ( 6
)
3 Numerical schemes

We introduce a constant space step ∆ x and a varying time step ∆t n . We define the mesh interfaces x j+1/2 = j∆ x, the cells C j = [x j-1/2 , x j+1/2 ], the cell centres x j and the intermediate times t n+1 = t n + ∆t n , where ∆t n is chosen according to the usual CFL condition. The Godunov method.

In this paragraph, we present the usual Godunov method based on the exact resolution of the Riemann problem including the source term, and associated with the initial data U(x, 0

) = U L 1 (x<x c ) + U R 1 (x>x c ) and U a (x, 0) = U a,L 1 (x<x c ) + U a,R 1 (x>x c ) .
The details are not reported here but one has been able to prove that the analytical solution is known and that its structure is resumed in Figures 12. In particular, the solution for the state consists of two waves, which can be either shocks or rarefaction waves, and whose speed can be computed exactly. On the other hand, the sensitivity has the same two-wave structure, however both of the waves are discontinuities. This simplification for the sensitivity is due to the fact that we are considering a reduced Euler system, under barotropic conditions (cf. [START_REF] Guinot | Sensitivity equations for hyperbolic conservation law-based flow models[END_REF]).

Since the state equations ( 1) are conservative, the Godunov method can be written with the classic update formula

U n+1 j = U n j - ∆t ∆ x (F(U * j+1/2 ) -F(U * j-1/2 )), (7)
where U * j-1/2 denotes the exact intermediate state variable in the Riemann solution associated to left and right states in the j -1-th and j-th cells. Due to the presence of the source term in the sensitivity equations, a conservative update formula like (7) cannot be obtained. However, recall that the structure of the sensitivity at each interface is very simple and made of two shocks. As a consequence, one can easily perform the average on the cells, provided that the slopes κ 1 and κ 2 , i.e. the slope of the red lines and the blue solid line in Figure 2, are known at each interface j -1/2. More precisely, we easily get

U n+1 a, j = U n a, j + ∆t ∆ x κ 2, j-1/2 (U * a, j-1/2 -U n a, j ) -κ 1, j+1/2 (U * a, j+1/2 -U n a, j ) , (8)
where the intermediate states U * a, j-1/2 and U * a, j+1/2 are known analytically and depend on U * j-1/2 and U * j+1/2 respectively. We remark that κ , j±1/2 depends on the solution structure: in case of shock it is given by the Rankine-Hugoniot conditions, in case of rarefaction wave it is the eigenvalue λ evaluated in the intermediate state U * j±1/2 . Let us observe that (8) already encompasses the source term to remove the Dirac delta function, therefore we do not need to discretise it. However, this correction is taken into account even in a numerical rarefaction profile (we recall that two points in a rarefaction wave are not necessarily linked by a rarefaction), and this is the cause of the failure of Godunov's method. In the next section, we present a Roe's method which uses shock detectors in order to overcome this problem.

x t x c U L U * U R (a) 1-shock-2-shock. x t x c U L U * U R U (b) 1-shock-2-rarefaction x t x c U L U * U R U (c) 1-rarefaction-2-shock. x t x c U L U * U R U 1 U 2 (d) 1-rarefaction-2-rarefaction.
Fig. 1: Configurations for the state variable U. A Roe-type method.

x t x c U a,L U * a U a,R (a) 1-shock-2-shock. x t x c U a,L U * a U a,R (b) 1-shock-2-shock. x t x c U a,L U * a U a,R (c) 1-shock-2-shock. x t x c U a,L U * a U a,R (
The proposed approximate Riemann solver of Roe-type consists of three constant states (say U L , U * and U R for the state, and U a,L , U * a and U a,R for the sensitivity) separated by two shock waves propagating at velocities

λ ROE L, j-1/2 = --(p(τ n j-1 ) -p(τ n j ))/(τ n j-1 -τ n j ) and λ ROE R, j-1/2 = -λ ROE L, j-1/2 if τ n j-1 = τ n j (and
of course ∓ -p (τ n j ) otherwise). In the following, we will use the notation

λ ROE j-1/2 = λ ROE R, j-1/2 = -λ ROE L, j-1/2 .
The fact that the velocities at each interface are equal and opposite in sign allows us to write at the interface j -1/2 the Harten, Lax and van Leer consistency relations for the state in the following way:

U * j-1/2 = 1 2 (U n j-1 + U n j ) - F(U n j ) -F(U n j-1 ) 2λ ROE j-1/2 . ( 9
)
Knowing U * j-1/2 and the velocity λ ROE j-1/2 at each interface one can average the solution value on the cells, obtaining the following update formula for the state:

U n+1 j = U n j - ∆t ∆ x (Φ n j+1/2 -Φ n j-1/2 ), (10)
where Φ n is the numerical flux and is defined as follows:

Φ n j-1/2 = F(U n j ) + F(U n j-1 ) 2 - λ ROE j-1/2 2 (U n j -U n j-1 )
.

We now consider the sensitivity equation, with the following source term, defined according to equations ( 4) to ( 6) and given by

∆ xS j-1/2 = -∂ a λ ROE j-1/2 (U * j-1/2 -U n j-1 )d 1, j-1 + ∂ a λ ROE j-1/2 (U n j -U * j-1/2 )d 2, j ,
where d , j is a shock detector to be defined, and d , j = 1 if there is an -shock in the j-th cell, it is zero elsewhere. The shock detector used in this work is based on the fact that the velocity u is always decreasing across a shock, and the covolume τ is decreasing across a 1-shock, whilst it is increasing across a 2-shock. Finally, the update formula for the sensitivity is the following:

U n+1 a, j = U n a, j - ∆t ∆ x (Φ n a, j+1/2 -Φ n a, j-1/2 ) + ∆t 2 (S j-1/2 + S j+1/2 ), (11)
where Φ n a is defined as follows:

Φ n a, j-1/2 = F a (U n j , U n a, j ) + F a (U n j-1 , U n a, j-1 ) 2 - λ ROE j-1/2
2 (U n a, j -U n a, j-1 ).

We remark that in (11), we add the source term S j-1/2 to both cells j-th and ( j -1)-th: indeed, it is defined starting from an integral balance done on both cells. Finally, we extended this scheme to the second order: we used a standard twostep Runge-Kutta method in time, whereas in space we used a MUSCL scheme adapted to take into account a second order discretization of the source term.

Numerical results

In this section we present some numerical results. The spatial domain is (0, 1), x c = 0.5, and final time T = 0.03.

The test case considered is a 1-shock-2-rarefaction, with the following initial conditions: U L = (0.7, 0) T , U R = (0.2, 0) T , U a,L = (0, 1) T , U a,R = (0, 0) T , and the parameter of interest is a = u L . Figure 3 and u a have a similar behaviour) at the final time T . As one can see, all the methods approximate well the state solution and the modified sensitivity formulation succeeds in removing the Dirac delta function located at point x ≈ 0.4 and noticeable without source term (scheme labelled "No source term" in Figure 3), however we remark that there are two main problems in the sensitivity solution: the shock corresponding to the state rarefaction is not captured and the value in the star zone is not correct. The first problem can be solved by refining the mesh or by using higher-order schemes (one can observe that the second order Roe method captures the discontinuity better); whilst the second problem, in our opinion, is due to the numerical diffusion in the shock. Figure 4 shows the convergence of the schemes: for the state we have the expected convergence; concerning the sensitivity, for coarser meshes the error is decreasing because its main part is in the rarefaction zone, however when this contribution becomes comparable with the error in the star zone, a plateau is reached.

Conclusion and discussion

The numerical results show that the modified sensitivity system here proposed is well defined and it allows us to achieve the main goal of this work, i.e. to have a sensitivity without Dirac delta function. However, the proposed modified formulation yields an incorrect value of the sensitivity in the star zone. Interestingly, numerical results show that even the exact Godunov method does not provide a correct solution, neither does a higher order scheme. In our opinion, this problem is due to the numerical diffusion in the shock. To illustrate this, we briefly show in Figure 5 the results obtained with a modified Godunov method based on sampling techniques, introduced in [START_REF] Chalons | Godunov scheme and sampling technique for computing phase transitions in traffic flow modeling[END_REF], which does not have any numerical diffusion and seems to provide a correct sensitivity in the star region.

  d) 1-shock-2-shock.

Fig. 2 :

 2 Fig. 2: Corresponding configurations for the sensitivity U a .

Fig. 3 :

 3 Fig. 3: Test case: 1-shock -2-rarefaction, ∆ x = 10 -3

Fig. 4 :Fig. 5 :

 45 Fig. 4: Convergence results for test case 1-shock -2-rarefaction.