Zichun Zhong
email: zichunzhong@wayne.edu

Xiaohu Guo

Wenping Wang
email: wenping@cs.hku.hk

Bruno Lévy

Jing Hua
email: jinghua@wayne.edu

Bruno Levy
email: bruno.levy@inria.fr

Computing a High-Dimensional Euclidean Embedding from an Arbitrary Smooth Riemannian Metric

This article presents a new method to compute a self-intersection free high-dimensional Euclidean embedding (SIFHDE) for surfaces and volumes equipped with an arbitrary Riemannian metric. It is already known that given a high-dimensional (high-d) embedding, one can easily compute an anisotropic Voronoi diagram by back-mapping it to 3D space. We show here how to solve the inverse problem, i.e., given an input metric, compute a smooth intersection-free high-d embedding of the input such that the pullback metric of the embedding matches the input metric. Our numerical solution mechanism matches the deformation gradient of the 3D → higher-d mapping with the given Riemannian metric. We demonstrate applications of the method, by being used to construct anisotropic Restricted Voronoi Diagram (RVD) and anisotropic meshing, that are otherwise extremely difficult to compute. In the SIFHDE-space constructed by our algorithm, difficult 3D anisotropic computations are replaced with simple Euclidean computations, resulting in an isotropic RVD and its dual mesh on this high-d embedding. The results are compared with the state-ofthe-art in anisotropic surface and volume meshings using several examples and evaluation metrics.

Introduction

Anisotropic meshes are important for not only improving the accuracy of the numerical simulations [START_REF] Alauzet | High-order sonic boom modeling based on adaptive methods[END_REF][START_REF] Narain | Adaptive anisotropic remeshing for cloth simulation[END_REF] but also better approximating the shapes [START_REF] Simpson | Anisotropic mesh transformations and optimal error control[END_REF].

Figure 1: Illustration of self-intersection free high-dimensional embedding (SIFHDE) for anisotropic Restricted Voronoi Diagram and meshing on Riemannian 2-and 3-manifolds. Left: For a simple example of a 2D Riemannian manifold with a metric field, we can compute a SIFHDE of it in a 3D Euclidean space. The unit circles and bean-shaped curves are iso-distance contours in metric domains. Based on the proposed SIFHDE framework, a high-dimensional embedding for an arbitrary 2-or 3-manifold can be computed. Then, we can generate the Voronoi diagrams and meshes in such embedding with Euclidean metric. Finally, when they are mapped to the original domains, the anisotropic Voronoi diagrams (Top Right) and the anisotropic meshes (Bottom Right) will exhibit the desired anisotropy. Note that Right images show anisotropic Voronoi diagrams and meshes generated by our method in different scenarios, such as metrics defined by a 2D analytic function, 3D surface curvature, and 3D analytic function.

Unlike isotropic meshes, where the elements are chosen to be as regular and uniform as possible, anisotropic meshes are designed with elongated mesh elements that follow user specified orientations and aspect ratios.

Anisotropic meshes are notoriously difficult to compute. A possible strategy to overcome the difficulty is to map the (complicated) anisotropic 3D space onto a higher-d space, where geometric computations are made simpler (Euclidean / isotropic). In a certain sense, "anisotropy is traded for additional dimensions". This idea is inspired by the celebrated Nash embedding theorem, that states that every Riemannian manifold can be isometrically embedded into some high-dimensional (high-d) Euclidean space [START_REF] Nash | C 1 -isometric embeddings[END_REF][START_REF] Kuiper | On C 1 -isometric embeddings I[END_REF]. In such highd embedding space, the metric is uniform and isotropic. When the isotropic mesh is computed in this embedding space and back-projected, the mesh elements on the original manifold will exhibit the desired anisotropic property. Recently, a few methods investigated this problem on surface meshing. Zhong et al. [ZGW + 13] introduced a particle-based meshing approach that conceptualizes the inter-particle energy optimization in a high-d "embedding space". This embedding space was used in their energy and force formulation without computing such a space explicitly. However, their method has some limitations due to lack of explicit embedding. Specifically, without an explicit embedding, K-Nearest Neighbors (K-NN) search of particles and anisotropic Restricted Voronoi Diagram (RVD) computation under a Riemannian metric space are very challenging and time-consuming. Another group of research focused on computing explicit embeddings, such as 3D embedding [START_REF] Panozzo | Frame fields: Anisotropic and non-orthogonal cross fields[END_REF], 5D or 6D embedding [START_REF] Lévy | Variational anisotropic surface meshing with Voronoi parallel linear enumeration[END_REF][START_REF] Dassi | Curvature-adapted remeshing of CAD surfaces[END_REF][START_REF] Dassi | Anisotropic finite element mesh adaptation via higher dimensional embedding[END_REF]. However, these methods have their own limitations, such as exhibiting self-intersections [START_REF] Panozzo | Frame fields: Anisotropic and non-orthogonal cross fields[END_REF], working only on specific dimensions [PPTSH14, LB12, DMS14, DSPS15], or using only specific embeddings (i.e., not taking into account arbitrary input Riemannian metrics) [START_REF] Lévy | Variational anisotropic surface meshing with Voronoi parallel linear enumeration[END_REF][START_REF] Dassi | Curvature-adapted remeshing of CAD surfaces[END_REF][START_REF] Dassi | Anisotropic finite element mesh adaptation via higher dimensional embedding[END_REF].

The novelty introduced by our method is that our input is a smooth arbitrary user-defined metric, from which the embedding is deduced, whereas previous work observed that some predefined shapes of the embedding result in some specific anisotropies when back-mapped to 3D space. Our main contribution is to provide a numerical method that solves the "inverse problem". It makes our method much more general in terms of supported metrics. Besides that, previous work considered the surfacic case, whereas we also experiment with the (significantly more challenging) volume meshing problem. Our algorithm works with arbitrary embedding dimension and avoids self-intersections. On the computed embedded surface or volume, we can directly optimize for a uniform, isotropic particle distribution and generate a high-quality isotropic RVD and its dual mesh on it. When mapped back to the original surface or volume, an anisotropic RVD as well as a simplicial mesh are obtained. Fig. 1 illustrates the method in the flat 2D setting, 3D surface setting and volumetric setting.

The embedding problem is formulated as an optimization that minimizes the deviation between the given metric and the deformation gradient of a map from the original surface / volume to the high-d embedded surface / volume. Our empirical experiments detailed further demonstrate that our method constructs an embedding with a small metric deviation error. The benefits of the high-d embedding are as follows:

the so-constructed high-d space can be used to compute anisotropic Voronoi diagrams for surfaces and volumes with a prescribed anisotropy. To our knowledge, this is the first time that an effective computational algorithm is proposed to construct such objects in the volumetric case with arbitrary metric; it introduces new computational strategies for some recent anisotropic meshing algorithms, such as particle optimization, RVD, and dual meshing in the high-d Euclidean embedding space. Our empirical evaluations in several scenarios demonstrate a practical utility, using a 8D embedding space.

2 Related Work

Anisotropic Triangulation

Anisotropic Centroidal Voronoi Tessellation: In order to compute well sampled anisotropic meshes, Du et al. [START_REF] Du | Anisotropic centroidal Voronoi tessellations and their applications[END_REF] further generalized the concept of Centroidal Voronoi Tessellation (CVT) to the anisotropic case -Anisotropic CVT (ACVT), and then computed its dual mesh graph. An Anisotropic Voronoi Diagram (AVD) with the given Riemannian metric needs to be constructed in each Lloyd iteration [START_REF] Lloyd | Least squares quantization in PCM[END_REF], which is a time-consuming operation. To make the computation faster, Valette et al. [START_REF] Valette | Generic remeshing of 3D triangular meshes with metric-dependent discrete Voronoi diagrams[END_REF] described a discrete approximation of ACVT by clustering the vertices of a dense pre-triangulation of the domain, at the expense of degraded mesh quality. Recently, Zhong et al. [START_REF] Zhong | Anisotropic surface meshing with conformal embedding[END_REF] provided a method to solve the anisotropic meshing by conformally mapping the metric surface to an appropriate 2D parametric domain and then compute CVT on it. But the limitation of the conformal embedding method is that it cannot handle surfaces with complicated topologies.

Surface Meshing in Higher Dimensional Space: Several solutions were proposed for certain classes of surface meshing problems by embedding them in high-d spaces [START_REF] Cañas | Surface remeshing in arbitrary codimensions[END_REF][START_REF] Boissonnat | Anisotropic diagrams: Labelle shewchuk approach revisited[END_REF][START_REF] Kovacs | Anisotropic quadrangulation[END_REF][START_REF] Lévy | Variational anisotropic surface meshing with Voronoi parallel linear enumeration[END_REF]. Lévy and Bonneel [START_REF] Lévy | Variational anisotropic surface meshing with Voronoi parallel linear enumeration[END_REF] extended the computation of CVT to a 6D space in order to achieve the curvatureadaptation. The main idea of their method is to transform the anisotropic meshing problem on a 3D surface to an isotropic one embedded in 6D space described by vertex positions and normals. Recently, Dassi et showed that by formulating the inter-particle energy optimization in a high-d space, such optimizations have very fast convergence without any need for the explicit control of particle population. But they still needed to compute the anisotropic RVD on the surface in 3D space in their final step of mesh generation, which is less robust and efficient. Another bottleneck in their framework is the search speed of neighboring particles. When the anisotropic stretching ratio is high, it uses a large set of K-NN candidates in the original space that contains many false positives. A comparison with our approach is shown in Sec. [START_REF] Boissonnat | Anisotropic Delaunay mesh generation[END_REF][START_REF] Boissonnat | Anisotropic Delaunay meshes of surfaces[END_REF] introduced a Delaunay refinement framework, which makes the star around each vertex x i to be consisting of the triangles that are exactly Delaunay for the metric associated with x i . In order to "stitch" the stars of neighboring vertices, refinement algorithms add new vertices gradually to reach the final anisotropic meshes. Recently, Rouxel-Labbé et al. [START_REF] Rouxel-Labbé | Discretized Riemannian Delaunay triangulations[END_REF] introduced an algorithm to compute discrete approximations of Riemannian Voronoi diagrams on 2-manifolds by using the numerical computation of geodesic distances. The main strength of these methods is that they come with theoretical proofs, both on termination (finite number of inserted points) and existence of a straight dual. However, since they use local criteria to insert points, the final result is often of lower quality than when using a global optimization [ZGW + 13, FLSG14]. Since they use the same representation of the metric as in [START_REF] Lévy | Variational anisotropic surface meshing with Voronoi parallel linear enumeration[END_REF] and as us (piecewise constant in the simplices), their method can be used as a post-processing after a global optimization, in order to benefit from both advantages (quality of global optimization and robustness / theoretical guarantees of Delaunay refinement).

Optimal Delaunay Triangulation: Another idea for anisotropic meshing is through Delaunay-based variational approach, by minimizing the error between a lifted quadratic target function and the linear interpolation of the constructed mesh. Chen et al. [START_REF] Chen | Optimal Delaunay triangulations[END_REF][START_REF] Chen | Optimal anisotropic meshes for minimizing interpolation errors in L p -norm[END_REF] proposed Optimal Delaunay Triangulation (ODT) method to minimize the error function for both isotropic and anisotropic meshing. Alliez et al. [START_REF] Alliez | Variational tetrahedral meshing[END_REF] presented a robust approach for isotropic tetrahedral meshing by exploiting the ODT formalism. In the context of varying densities, the objective function was carefully analyzed and improved in [CWL + 14]. ODT can be considered as lifting the points onto a paraboloid and estimating volume integrals on this paraboloid. Budninskiy et al. [BLdG + 16] presented a variational method to take into account an anisotropy by using a convex function instead of a paraboloid. While elegant and interesting from a theoretical point of view, their method is limited to a small class of anisotropies stemming from convex functions. The constraint is not reasonable in practical applications (for instance, it could not account for the anisotropy in Fig. 1-Right, nor with any non-convex geometry). Fu et al. [START_REF] Fu | Anisotropic simplicial meshing using local convex functions[END_REF] proposed a Locally Convex Triangulation (LCT) method to compute the anisotropic meshes based on constructing convex functions that locally match the specified Riemannian metric. We compare the quality of our generated mesh with Fu et al.'s method in Sec. 6.

Anisotropic Tetrahedralization

Theoretically, it is possible to extend ACVT-based method for anisotropic meshing from 2D domain or surface [DW05a, VCP08, ZSJG14] into 3D volume, but this group of methods needs to compute AVD iteratively in the ambient 3D space with specified metrics, which makes the computation very complicated and impractical. In computer graphics and geometric modeling areas, there is few literature in anisotropic tetrahedral meshing: Dobrzynski and Frey's local mesh adaptation method based on an extension of Delaunay kernel [START_REF] Dobrzynski | Anisotropic Delaunay mesh adaptation for unsteady simulations[END_REF], Klingner et al.'s aggressive optimization [START_REF] Matthew | Agressive tetrahedral mesh improvement[END_REF][START_REF] Matthew | Tetrahedral Mesh Improvement[END_REF] and Fu et al.'s LCT method [START_REF] Fu | Anisotropic simplicial meshing using local convex functions[END_REF] on some simple 3D models. Yamakawa and Shimada [START_REF] Yamakawa | High quality anisotropic tetrahedral mesh generation via packing ellipsoidal bubbles[END_REF] proposed an ellipsoidal bubble packing method, but inserting or deleting particles / bubbles is necessary in the computation. In mechanical and biomedical engineering, anisotropic meshing is widely used in computational fluid dynamics (CFD) [FA05, AL16, LL16] and blood flow simulations [START_REF] Marchandise | Cardiovascular and lung mesh generation based on centerlines[END_REF][START_REF] Sauvage | Metric field construction for anisotropic mesh adaptation with application to blood flow simulations[END_REF]. However, their methods are numerical and error estimation-based approaches. For instance, they are very time-consuming, since partial differential equations (PDEs) are solved at each solver time step, which make it difficult to efficiently compute complicated 3D anisotropic volume meshes.

Other Related Embedding Work

Bronstein et al. [START_REF] Bronstein | Multigrid multidimensional scaling[END_REF][START_REF] Bronstein | A multigrid approach for multi-dimensional scaling[END_REF] demonstrated the performance of multidimensional scaling (MDS). Instead of raising the dimensions, they mapped the original metric domain into a parameter domain (similar to a conformal embedding). Verma [START_REF] Verma | Distance preserving embeddings for general ndimensional manifolds[END_REF] proposed the distance preserving embeddings for general n-dimensional manifolds in the machine learning field. It applies "spiraling perturbations / corrections" for globally isometric mapping. This cannot compute a smooth high-d embedding for our computer graphics research and applications. Recently, Panozzo et al. [START_REF] Panozzo | Frame fields: Anisotropic and non-orthogonal cross fields[END_REF] proposed a 3D embedding framework with self-intersections to compute a surface deformation that warps a frame field into a cross field, and used in the adaptive quad meshing. In their method, the parameterization is used to handle the final meshing with self-intersections. However, our proposed method does not need an extra step to parameterize the embedding for meshing, since we can directly produce the Voronoi diagrams and meshes on the computed SIFHDE. Another limitation of their method is that the frame-driven deformation for 3D embedding computation may deteriorate the quality of the input triangulation, especially when the stretching ratio is high, such as ≥ 10, which leads to failure. However, our method theoretically can compute the embedding for any high-stretching ratio case, since we provide enough degrees of freedom in higher dimensions. We compare the embedding quality of our method with Panozzo et al.'s method in Sec. 6. Besides surface cases, our proposed high-d embedding approach can work on volume cases, which is an unexplored topic in the previous work.

3 High-D Embedding

Anisotropy and High-D Embedding

We consider that a metric M(.), i.e., a symmetric positive definite (SPD) bilinear form, is defined over the surface or volume domain Ω ⊂ R m . In other words, at a given point x ∈ Ω, the dot product between two vectors a and b is given by a, b M(x) . The metric can be represented by a symmetric m × m matrix M(x), in which case the dot product becomes:

a, b M(x) = a T M(x)b. (1)
Through Singular Value Decomposition (SVD), the metric matrix M(x) can be decomposed into:

M(x) = R(x) T S(x) 2 R(x), (2)
where the diagonal matrix S(x) 2 contains its ordered eigenvalues (i.e., a scaling field), and the orthogonal matrix R(x) contains its eigenvectors (i.e., a rotation field). In the surface case [DW05a, ZGW + 13], they are defined on the tangent spaces of the surface. s2 s1 is defined as stretching ratio in the surface metric field M(x), where s 1 and s 2 are the two diagonal items in S(x) and s 1 ≤ s 2 . Similarly, in volume case [START_REF] Yamakawa | High quality anisotropic tetrahedral mesh generation via packing ellipsoidal bubbles[END_REF], they represent the scaling and rotation in a 3D volume space.

Metric through High-D Embedding: For an arbitrary metric field M(x) defined on the surface or volume Ω ⊂ R m (i.e., Riemannian 2-or 3-manifold), Nash theorem [START_REF] Nash | C 1 -isometric embeddings[END_REF] states that there exists a high-d space R m (i.e., m ≤ m), in which the surface or volume can be embedded with Euclidean metric as Ω ⊂ R m . In this article, Ω denotes the "embedded surface or volume". Consider a mapping denoted by φ : Ω → Ω.

From the definition of anisotropy above, it can be seen that introducing anisotropy means changing the definition of the dot product. If we consider two vectors a and b from a given location x ∈ Ω, then they are transformed into a = J(x)a and b = J(x)b on the high-d embedded surface or volume Ω, where J(x) denotes the Jacobian matrix of mapping φ at x. The dot product between a and b is given by the pullback metric of Φ, defined as:

a, b = a T J(x) T J(x)b = a T M(x)b.
(3)

Importance of High-D Embedding: (1) In a high-d space, more degrees of freedom are available to deform and embed the given surface or volume, so that better embedding quality is obtained. Experiments are shown in Sec. 6. (2) When using the high-d space, we can avoid self-intersections of the embedded surface or volume as discussed in Sec. 3.4.2, instead of embedding them in the original space. (3) By generating a high-d Euclidean embedding without selfintersections, we are able to drastically simplify several Riemannian geometric applications, such as computing anisotropic RVD and meshing on surface and volume with high-quality elements using only isotropic Euclidean computations. This is both more efficient and more accurate than the computations in metric space. Details are presented in Sec. 4.

High-D Embedding Transformation

In this work, we assume the surfaces are all represented as triangle meshes and the volumes are all represented as tetrahedron meshes. The two surfaces / volumes Ω (in the original space) and Ω (in the high-d space) share the same number of vertices and the same connectivity between vertices. The only difference between them is their vertex coordinates and dimensions; thus, we assume that the two surfaces / volumes have the same indices of vertices and triangles / tetrahedrons. In the following subsection, we discuss the surface and volume cases, respectively.

Surface: For a triangle j on Ω and Ω, let {v j1 ,v j2 , v j3 } and {v j1 , v j2 , v j3 } denote their vertices, respectively. Their corresponding edge vectors can be represented as the following two matrices W j and W j :

W j = [v j2 -v j1 , v j3 -v j1 ] W j = [v j2 -v j1 , v j3 -v j1 ]. (4)
Volume: For a tetrahedron j on Ω and Ω, let {v j1 , v j2 , v j3 , v j4 } and {v j1 , v j2 , v j3 , v j4 } denote their vertices, respectively. Their corresponding edge vectors can be represented as the following two matrices W j and W j :

W j = [v j2 -v j1 , v j3 -v j1 , v j4 -v j1 ] W j = [v j2 -v j1 , v j3 -v j1 , v j4 -v j1 ].
(5)

Their relationship can be represented as:

W j = J j W j , (6)
where J j is the Jacobian transformation matrix for triangle or tetrahedron j.

In what follows, we assume that the original metric is smooth over the surface or the volume, and sampled as a constant symmetric matrix M j attached to each mesh simplex j. Note that in practice, it may be specified in the input at each vertex. In this case, several techniques can be used to interpolate the metric and deduce a reasonable constant value on each simplex (see [START_REF] Courty | Continuous metric and mesh adaptation[END_REF] for a discussion).

We now elaborate on the relationship between the Jacobian matrix of the mapping J j and the metric M j .

To better understand J j , it is necessary to explore the relationship between the original triangle or tetrahedron j on Ω and the embedded triangle or tetrahedron j on Ω. J j is an m × m matrix, and it can be represented as the product of a rotation in the high-d embedding space, and a scaling and rotation in the original space. That is:

W j = U j Sj Rj 0 W j , (7)
where U j is an m × m unitary matrix (i.e., a rotation matrix in R m ); 0 is a (m -m) × m block of zeros; S j and R j are the diagonal scaling matrix and rotation matrix extracted from SVD of the metric M j as shown in Eq. ( 2). Thus, the embedding transformation J j is:

J j = U j Sj Rj 0 . ( 8
)
By denoting Q j = Sj Rj 0 , we can simply represent J j as:

J j = U j Q j . (9)

High-D Deformation Gradient

Intuitively, the transformation between the original surface or volume Ω and its high-d embedded one Ω can be considered as the deformation of Ω. It is represented by the field of the deformation gradient over the surface or volume as being intuitively by Sumner and Popović [START_REF] Sumner | Deformation transfer for triangle meshes[END_REF]. This field is defined per triangle face or tetrahedron volume. We extend their idea from 3D surface to our high-d embedding as well as volume cases. Given a manifold with no deformation (i.e., the original surface or volume Ω) and its corresponding deformed manifold (i.e., the embedded surface or volume Ω), the deformation gradient T j for triangle or tetrahedron j (T j is an m × m matrix) can be defined by:

T j W j = W j . (10)
In order to obtain T j , we need to compute the inverse of W j . Since the surface and volume have different representations of edge vectors, there are two cases as follows:

(1) Surface: As for the surface, W j is a 3 × 2 matrix, so we cannot directly compute its inverse matrix. With QR factorization of W j [START_REF] Golub | Matrix Computations[END_REF]:

W j = O j Pj 0 = [ OjαO jβ ] Pj 0 = O jα P j , (11)
where P j is a 2 × 2 upper triangular matrix, O j is an m × m orthogonal matrix, O jα is an m × 2 matrix, and O jβ is an m × (m -2) matrix, where m = 3. Then we have the minimum norm solution of the deformation gradient for triangle j:

T j = W j P -1 j O T jα . (12)
(2) Volume: As for the volume, W j is a 3 × 3 square matrix, so the inverse of it can be directly computed as:

T j = W j W -1 j . (13)

Embedding Optimization

From Eq. ( 9) and Eq. ( 12) (or Eq. ( 13)), we can clearly see that in essence, the high-d embedding transformation J j and the high-d deformation gradient T j are the same. For J j in Eq. ( 9), Q j is coming from the given metric M j on the surface or volume, while U j is the high-d rotation matrix that is unknown. For T j in Eq. ( 12), P j and O jα are from the edge vector W j of the original surface, or in Eq. ( 13), W -1 j is from the edge vector W j of the original volume, while W j is the edge vector of the high-d embedded triangle or tetrahedron that is unknown.

Thus, knowing either of these two matrices, i.e., J j and T j , guides us in computing the other. We can formulate an expression to minimize the function E em defined as the difference between these two matrices:

E em (v 1 , • • • , v nv ) = min n ele j=1 T j -J j 2 F , (14)
where n v is the number of vertices, n ele is the number of mesh elements, i.e., triangles or tetrahedrons, and {v 1 , • • • , v nv } are the vertices of the high-d embedded surface or volume Ω. The matrix norm ||.|| F is the Frobenius norm. By substituting Eq. ( 9) and Eq. ( 12) (or Eq. ( 13)) into Eq. ( 14), we can get the more detailed expressions for this objective function :

Surface:

E em (v 1 , • • • , v nv ) = min n ele j=1 W j P -1 j O T jα -U j Q j 2 F . (15)
Volume:

E em (v 1 , • • • , v nv ) = min n ele j=1 W j W -1 j -U j Q j 2 F . (16)

Regularization Term

It is well known that solutions of Nash equations often present wrinkles, known as corrugations, that can form fractal patterns (see [START_REF] Borelli | Flat tori in three-dimensional space and convex integration[END_REF] and references herein on convex integration). Thus, a direct minimization of the criterion in Eq. (15) or Eq. ( 16) may lead to an oscillating high-d embedding, with wrinkles, that may have an influence on the stability / robustness of the subsequent geometric operations: in the high-d embedding space, we would like to use the Euclidean distance to approximate the metric distance on the original surface or volume, in order to efficiently compute: (1) K-NN for particle optimization in the high-d space; and (2) high-d RVD and dual mesh based on uniformly optimized particles. Thus the embedded surface or volume has to be smooth enough, in order to ensure that the K-NN and RVD computed with Euclidean distances in the high-d space are sufficiently accurate.

Thus, to regularize the embedding, we add the following term:

E reg (v 1 , • • • , v nv ) = n v-b i=1 m d=1 ( k∈N (i) (v d k -v d i ) |N (i)| ) 2 , ( 17
)
where n v-b is the total number of vertices excluding those on the open boundaries in surface case or boundary surface in volume case, N (i) is the set of one-ring neighbors of vertex i, |N (i)| is the size of set N (i), and v d i , v d k are the d-th dimensional coordinates of vertices v i and v k . The regularity term E reg is a summation of the square of Graph Laplacian operations over every vertex in the embedding space except those vertices on the open boundaries in surface case or boundary surface in volume case. Now our regularized objective function for the embedding optimization includes two terms: the similarity between two transformations and the regularity used to achieve smoothness of the embedding:

E total = E em + µE reg , ( 18
)
where µ is a weighting factor to balance the similarity and regularity terms during optimization. In our experiments, µ is set to be 100.

Avoiding Intersections and Choosing Target Dimension

According to Nash embedding theorem, using the mapping φ : Ω → Ω given by

v 1:m → (v 1:m , v m+1 , • • • , v m
), one obtains an embedding without self-intersections.

In our case, to avoid self-intersections that may come from numerical approximations or from the regularization term, in addition we keep the original 3D coordinates, thus, in a certain sense, the additional nD coordinates that we compute act as "correcting terms" to solve for the metric. Clearly, if the original 3D surface is free of self-intersections, then it is also the case of our (n+3)D embedding. Compared with methods that avoid self-intersections a-posteriori in a 3D embedding [START_REF] Panozzo | Frame fields: Anisotropic and non-orthogonal cross fields[END_REF], this is both simpler, and leaves more degrees of freedom to fit the user-specified metric. Generally speaking, the higher the dimension is, the smaller the embedding error is, since more degrees of freedom are provided to deform and embed the given surface or volume in such an embedding space; it is necessary to choose the dimension in a way that balances the computational efficiency and the embedding errors. Our empirical results showed that 8D is a reasonable target dimension for embedding most general metric surfaces and volumes in small errors and without any self-intersection. More details and empirical validation will be given in Sec. 6 (particularly in Sec. 6.1.2).

Numerical Solution Mechanism

The optimization defined above is a non-linear problem with two unknown parameters (i.e., W j and U j ). We use an iterative method to obtain the optimal solution of W j . It should be noted that the optimized vertex coordinates v 1 , • • • , v nv of the embedded surface or volume are included in the edge vectors W j . Initially, for each vertex i, we set its high-d embedding coordinates as

v i = (v 1:m i , v m+1 i , • • • , v m i )
, where the first m dimensions are the same as the original surface or volume's coordinates, and use a small random perturbation to each higher dimensional coordinate (from dimension m + 1 to dimension m), i.e., not all equal to zero. The insightful reason is that at first we fix W j to compute the rotation matrix U j in high-d space. U j is computed from Eq. ( 19) or Eq. (20) as follows, which includes W j . If the higher dimensional coordinates are all zeros, W j will be degraded in the original space R m , leading to the high-d rotation U j restricted in the original space as well, so that we use this trick to add the small perturbation value to let the optimization "jump" out of the original space.

When the vertex coordinates

{v i |i = 1, • • • , n v } are fixed (i.
e., all W j s are fixed), the optimal U j can be computed as the orthogonal factor of the following Polar decomposition: Surface:

U j = P olar(W j P -1 j O T jα Q T j ). (19)
Volume:

U j = P olar(W j W -1 j Q T j ). (20)
The detailed derivations are given in the supplementary material (Appendix A). When all U j s are fixed, we can compute the optimal vertex coordinates

{v i |i = 1, • • • , n v }, where v 1:3 i = v 1:3
i , by solving a linear system, since the objective function E total is a quadratic form of the vertex coordinates. This is similar to the method used by Sumner and Popović [START_REF] Sumner | Deformation transfer for triangle meshes[END_REF]. Deformation gradients are invariant to translation, so there are infinite solutions to this optimization problem: all translations of one optimal solution are also optimal. Thus we fix the coordinates of the first vertex in the mesh in order to make the solution unique.

Our goal is to ensure that the regularized objective function always decreases during optimization. The optimization strategy is similar to [START_REF] Sorkine | As-rigid-as-possible surface modeling[END_REF][START_REF] Panozzo | Frame fields: Anisotropic and non-orthogonal cross fields[END_REF]. In summary, as shown in the following Alg. 1, we can optimize the vertex coordinates v 1 , • • • , v nv of the embedded surface or volume iteratively, until convergence by satisfying a specified stopping condition, e.g., the magnitude of the gradient is smaller than a threshold or the desired number of iterations is reached. In our experiments, we use the number of iterations as 50 for the stopping condition. One example of the evolution of the objective function during optimization is given in Fig. 2 (b), which can reach a satisfactory small embedding error.

Anisotropic RVD and Meshing

In order to demonstrate the importance and usefulness of the proposed SIFHDE, we demonstrate two applications, i.e., computing anisotropic RVD and anisotropic Since the high-d Euclidean embedding without self-intersections has been computed in the previous section (Sec. 3), in the following section, under the isotropic metric space, we use the particle optimization technique (Sec. 4.1) to compute isotropic RVD and its dual Delaunay triangulation on the high-d embedded surface or volume (Sec. 4.2), due to its efficiency to generate regular hexagonal patterns of particles on 2-manifold and regular body-centered-cubic (BCC) lattices on 3-manifold, which are similar to the results of CVT [DFG99, LWL + 09, DW05b]. Moreover, the particle-based method is easily formulated in the high-d space. Finally the generated RVD and triangulation / tetrahedralization can be mapped from the high-d embedding space to the original space, on which the RVD and triangulation / tetrahedralization will exhibit the desired anisotropic property.

Particle Optimization with Efficient K-NN in High-D Embedding

In this section, we discuss the isotropic particle optimization on the explicitlycomputed high-d embedded surface or volume, which is the main difference (2) during the particle optimization, they had to have a very large range for searching neighboring particles in the original metric space when the anisotropic stretching ratio is high, which is less efficient. Besides that, they only worked on surface case without considering volume case. However, based on our method, once the high-d embedding is computed as discussed in Sec. 3, we can directly define the inter-particle energy and force in such embedding space. The goal is to let particles reach isotropic uniform distribution at their equilibrium state, where the particle distribution on the original surface or volume will exhibit the desired anisotropic property.

In our method, we formulate the isotropic inter-particle energy and force with explicit coordinates of the particles in the high-d embedding space, which is based on [ZGW + 13]. The energy E pq between particles p and q in the high-d embedding space is defined as:

E pq = e -xp -xq 2 4σ 2
, where σ is kernel width, and x p and x q are particle positions in the embedding space. Note that our goal is to let the particles uniformly and isotropically sampled on Ω, so σ = 0.3 d |Ω|/n, as suggested in [ZGW + 13], where d = 2 in surface and d = 3 in volume, |Ω| denotes the area or the volume of Ω in the embedding space, n is the number of particles. The force applied from particle p to particle q in the embedding space is defined as the gradient of the inter-particle energy. To sum up all inter-particle energies, the computational complexity of particle optimization is O(cn p ), where c is a constant, i.e., the average number of K-NN for each particle, and n p is total number of particles (i.e., the number of RVD cells / final mesh vertices). The k-d tree data structure is used to compute K-NN. Since n p is provided by user, c is the key factor to affect the efficiency of the computation. With the computed high-d Euclidean embedded surface or volume, the particle positions are defined and optimized on such embedding, which helps us to determine the c around 20∼30, i.e., nearest neighbors in Euclidean metric space within five standard deviations of the Gaussian energy kernel (5σ). If without the high-d Euclidean embedding, hundreds or thousands of neighbors need to be checked and then prune the spurious neighbors when the stretching ratio in the given metric field is high, such as ≥ 10. The K-NN comparison experiments are given in Sec. 6.2. Finally, we use the L-BFGS algorithm [START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF] to optimize the particle positions constrained on the high-d embedded surface or inside the high-d embedded volume.

Restricted Voronoi Diagram and Mesh Generation

High-D RVD and Meshing: The input of high-d RVD computation are the triangle / tetrahedral mesh of the high-d embedded surface / volume Ω and the sites (i.e., optimized particle positions in high-d space) of high-d Voronoi cells. The key task is to identify the high-d Voronoi cells that overlap each triangle / tetrahedron of the embedded surface / volume Ω and compute their intersections. The RVD computation on surface is based on Yan et al.'s method [YLL + 09] and the volume version is based on Lévy and Bonneel's method [START_REF] Lévy | Variational anisotropic surface meshing with Voronoi parallel linear enumeration[END_REF] with the exact geometric predicates from [START_REF] Lévy | Robustness and efficiency of geometric programs: The predicate construction kit (PCK)[END_REF], and then extended in high-d space. All these computations are done under the Euclidean metric, which is easy and efficient, even in the high-d space. In particular, note that the geometric predicates in [START_REF] Lévy | Robustness and efficiency of geometric programs: The predicate construction kit (PCK)[END_REF] compute barycentric coordinates, that depend on the intrinsic dimension (2 for surfaces, 3 for volumes) rather than ambient dimension (highd). Without using the high-d Euclidean embedding, one would need to compute anisotropic RVD and its dual anisotropic mesh on the original surface with the given metric, requiring a high-resolution tessellation of the input surface / volume and additional costly computations, such as computation of anisotropic Voronoi cell boundaries, search for neighboring Voronoi sites in anisotropic metric space, . . . . Once the RVD is obtained, we can easily compute its dual graph, i.e., Restricted Delaunay Triangulation (RDT) as follows.

Anisotropic RVD and Mesh: Using the barycentric coordinates of each output site or vertex, we can back-project the RVD and RDT from the highd embedding space onto the original space, and generate the final anisotropic RVD and mesh. As we mentioned before, if the given Riemannian metrics are smooth, the proposed high-d embedding framework can compute a smooth and small-error embedding in higher dimensions. The algorithm robustly computes the RVD using filtered geometric predicates and symbolic perturbation to resolve degeneracies [START_REF] Lévy | Robustness and efficiency of geometric programs: The predicate construction kit (PCK)[END_REF]. Note that in general, there is no guarantee that the dual mesh (i.e., the associated RDT) will not have inverted elements (negative Jacobian) when back-projected to 3D. Whenever such an inverted element is detected, our implementation inserts additional points using the provably terminating algorithm in [START_REF] Rouxel-Labbé | Discretized Riemannian Delaunay triangulations[END_REF], that uses the same discretization of the metric as us (also used in [START_REF] Lévy | Variational anisotropic surface meshing with Voronoi parallel linear enumeration[END_REF]). In practice, we did not observe such inverted elements in our empirical experiments (the refinement step was never triggered).

Evaluations

Embedding Quality

In order to evaluate the optimized high-d embedding of the surface or volume equipped with a given metric (in practice, the given metric is specified in the input mesh at each vertex), we consider the edge lengths of the input triangular surface or tetrahedral volume under the given metric, and use them as the ground truth to evaluate computed embedding results.

For each edge with vertices v a and v b , we use the average metric of two vertices Q ab = (Q(v a ) + Q(v b ))/2 to approximate the metric for this edge. Then the length of each edge can be computed according to its own Q ab and we use it as the ground truth. The absolute edge length error is the absolute difference between the edge length of computed embedding and the ground truth. The relative edge length error is the percentage of the absolute difference with respect to the ground truth. In our experiments, we only use relative edge length error to measure the computed embedding, because different models may have different sizes and scales. The criteria for evaluating the embedding quality are: L rel max is the maximal value of relative edge length errors of all embedded triangles / tetrahedrons; and L rel avg is the average value of relative edge length errors of all embedded triangles / tetrahedrons. The color-coded diagram of the relative edge length errors is rendered.

K-NN Efficiency

The following criteria of K-NN efficiency are used: T knn is computational time on K-NN searching for all particles; N avg is the average number of nearest neighbors to be searched; N max is the maximal number of nearest neighbors to be searched. We use the above criteria to compare Zhong et al.'s anisotropic particle optimization [ZGW + 13] with our proposed isotropic particle optimization on the high-d Euclidean embedding in Sec 6.2.

Anisotropic Mesh Quality

Since RVD is not straightforward to be evaluated as compared with its dual triangular and tetrahedral meshes, in the experiments, we will focus on the mesh quality for evaluations.

Surface: In the anisotropic triangular meshing, for each triangle abc in the final mesh, we use its approximated metric

Q( abc ) = (Q(x a ) + Q(x b ) + Q(x c
))/3 to affine-transform it from the original anisotropic space into the Euclidean space. After that, we employ the following isotropic triangular criteria, as suggested by Frey and Borouchaki [START_REF] Frey | Surface mesh evaluation[END_REF] and used in Zhong et al. [ZGW + 13] and Fu et al. [START_REF] Fu | Anisotropic simplicial meshing using local convex functions[END_REF]'s recent work, to evaluate the quality of generated anisotropic triangular mesh: The quality of a triangle is measured by G = 2 √ 3 S ph , where S is the triangle area, p is its half-perimeter, and h is the length of its longest edge; G min is the minimal quality of all triangles; G avg is the average quality of all triangles; θ min is the smallest angle of the minimal angles of all triangles; θ avg is the average angle of the minimal angles of all triangles; % <30 • is the percentage of triangles with their minimal angles smaller than 30 • ; The angle histogram is also provided to show the distribution of angles for all triangles. θ avg should be 60 • if it is a regular triangle. G is between 0 and 1, where 0 denotes a skinny triangle and 1 denotes a regular triangle.

Volume: In the anisotropic tetrahedral mesh, for each tetrahedron tet abcd , we have its approximated metric

Q(tet abcd ) = (Q(x a ) + Q(x b ) + Q(x c ) + Q(x d ))/4.
Then we use the corresponding Q(tet abcd ) to affine-transform the tetrahedron in the Euclidean space, and employ the following isotropic tetrahedral mesh quality measurements [DVS + 09, YWLL13] to evaluate the quality of generated anisotropic tetrahedral mesh: The quality of a tetrahedron is measured by

G = 12 3 √ 9V 2 l 2 i,j
, where V is the volume of the tetrahedron, and l i,j is the length of the edge which connects vertices v i and v j . G min , G avg are the minimal and average qualities of all tetrahedrons. θ min , θ avg are the smallest and average angles of the minimal dihedral angles of all tetrahedrons. % <15 • is the percentage of tetrahedrons with their dihedral angles smaller than 15 • , which are considered as slivers. The angle histogram is also provided to show the distribution of minimal dihedral angles for all tetrahedrons. θ avg should be ≈ 70.53 • if it is a regular tetrahedron. G is between 0 and 1, where 0 denotes a sliver and 1 denotes a regular tetrahedron.

Results

We implement our algorithms by using Microsoft Visual C++ 2013. The embedding and mesh quality evaluations are implemented with Matlab R2015a. For the hardware platform, the experiments are run on a desktop computer with Intel(R) Core(TM) i7-6700 CPU with 3.40GHz, 32GB DDR4 RAM.

Riemannian Metric Design: In our experiments, users firstly design a smooth scaling field S(x) and a rotation field R(x) that is smooth in regions other than some singularities, and then compose them to M(x) = R(x) T S(x) 2 R(x), which is similar to the approach mentioned in [START_REF] Du | Anisotropic centroidal Voronoi tessellations and their applications[END_REF][START_REF] Yamakawa | High quality anisotropic tetrahedral mesh generation via packing ellipsoidal bubbles[END_REF]. Since this article focuses on both 2-and 3-manifolds, the metric design has some slight differences as follows.

For the anisotropic meshing in 2D domains, we use the following 2×2 metric tensor:

M = [v 1 , v 2 ]diag(s 1 2 , s 2 2 )[v 1 , v 2 ] T , (21)
where v 1 and v 2 are two orthogonal principal directions, composing a rotation field. s 1 and s 2 are two user-specified stretching factors along principal directions.

For the anisotropic meshing on 3D surfaces, we use the following 3×3 metric tensor:

M = [v min , v max , n]diag(s 1 2 , s 2 2 , 0)[v min , v max , n] T , (22)
where v min and v max are the directions of the principal curvatures, n is the unit surface normal. s 1 and s 2 are two user-specified stretching factors along principal curvature directions. Since the surface models are computed by curvaturebased metric tensor fields. We use the metric of Eq. ( 22) with s 1 = √ K min and s 2 = √ K max , where K min and K max are the principal curvatures. We set small thresholds to preserve both K min and K max not vanishing. As suggested by Alliez et al. [ACSD + 03], Laplacian smoothing is applied to both the stretching factors and directions, to ensure smoothness of the input metric field on the surface.

For the anisotropic meshing in 3D volumes, we use the following 3 × 3 metric tensor:

M = [v 1 , v 2 , v 3 ]diag(s 1 2 , s 2 2 , s 3 2 )[v 1 , v 2 , v 3 ] T , ( 23
)
where v 1 , v 2 , and v 3 are three orthogonal principal directions. s 1 , s 2 , and s 3 are three user-specified stretching factors along principal directions.

SIFHDE Computation

Importance of Higher Dimensions

In order to demonstrate the importance of the high-d embedding, instead of embedding computed in the original space. We use two examples in detail to explain its advantages. One example is a 2D domain with a Riemannian metric field. Fig. 2 shows a simple 3D self-intersection free embedding computed by our proposed method for a 2D domain equipped with an anisotropic metric, which can "perfectly" embed the original 2D domain into a Gaussian surface in 3D space. Since we have both the original domain with the metric and its ground truth of the embedded surface (the target dimension is known, i.e., 3D), we can explicitly evaluate our high-d embedding optimization framework both quantitatively and visually.

The pipeline of the proposed work is given in Fig. 2, i.e., to compute the self-intersection free high-d embedding, isotropic particle optimization on such embedding, high-d isotropic RVD and dual meshing, final anisotropic RVD and mesh. The anisotropic metric field is M(x) = R(x) T diag{Stretch(x) 2 , 1}R(x), where the stretching field Stretch(x) ∈ [1, 3.85] with the rotation field R(x). Our proposed embedding optimization framework converges superbly fast, i.e., 20 iterations in this 3D embedding. The embedding quality is quite satisfactory, i.e., L rel avg = 0.92% and L rel max = 10.58%. Besides the high quality, there is no selfintersection, due to the strategy we discuss in Sec. 3.4.2, i.e., fixing the original 2D coordinates and deforming the new added 3rd coordinates in the embedding computation. Visually, the shape of the embedding result at 20 iterations (Fig. 2 b) looks very similar with the shape of the ground truth Gaussian surface (Fig. 2 a). By using this 3D embedding, we can compute the uniform distribution of particles (Fig. 2 c) for RVD and dual meshing on it (Fig. 2 d), and finally generate high-quality anisotropic RVD and triangular mesh with 2000 output samples as shown in Fig. 2 (e). Note that without using an extra higher dimension (i.e., the 3rd coordinate), the output embedding will have self-intersections, since the input metric controls the deformations of each triangle in the original 2D space, which will very easily cause self-intersections, especially when the metric field has large stretching ratios and directions. Regularizing the embedding to avoid self-intersections is not a good strategy, since it would dramatically increase the error of the embedding either.

Another example is a 3D Cyclide surface with a Riemannian metric field, which is designed according to the surface curvature. The stretching ratio is s2 s1 ∈ [1.6, 9.4]. Without a high-d embedding, i.e., using a 3D embedding as suggested by Panozzo et al. [START_REF] Panozzo | Frame fields: Anisotropic and non-orthogonal cross fields[END_REF], in Fig. 3, though the embedding errors are relatively small, i.e., L rel avg = 2.36% and L rel max = 30.62% (still larger than our 8D embedding errors as shown in Fig. 4), there are many self-intersecting faces in the computed 3D embedding (i.e., 1146 self-intersecting faces in this Cyclide model), which are harmful for the particle optimization, RVD and dual mesh computations. The major fundamental difference between their framework and ours is that they deform the embedding directly in the original 3D coordinates, which may lead self-intersections; however, we deform the embedding in highd space with fixing the first three coordinates to avoid self-intersections (as discussed in Sec. 3.4.2). Furthermore, we also investigate some other surface models and Tab. 1 shows that self-intersecting faces do widely exist in the 3D embedding (i.e., Panozzo et al.'s method [START_REF] Panozzo | Frame fields: Anisotropic and non-orthogonal cross fields[END_REF]); however, our high-d embedding framework does not have such problem. Due to the page limit, we do not provide the quantitative statistics of self-intersecting faces in 2D domain and 3D volume models, but it is clear that the same conclusion can be drawn.

Choosing the Dimension of the Embedding

We elaborate on how to choose the dimension of the high-d target for computing the SIFHDE with small relative edge length errors. In general, we do not know beforehand the required dimension for the high-d embedding. From Nash theorem, it is known that for constructing a C 1 isometric embedding of a nD metric, 2n dimensions suffice (Nash-Kuiper). For higher smoothness (C k , k ≥ 2), the bound on embedding dimension becomes n(3n + 11)/2 (Nash-Moser). With our regularization term, we constrain the embedding to be C 2 and thus avoid the oscillations of C 1 embeddings, but the Nash-Moser theoretical bound becomes 17D for surfaces and 30D for volumes. We now evaluate the required embedding dimension in practice for surfaces and volumes and show empirically that in both cases 8D suffices to accurately account for the input metric:

Surface Case: We conduct different experiments with surfaces. The associated metric deviation errors of a Cyclide surface are plotted in Fig. 4. It is noted that in our case, when using embedding dimensions larger than 8D, the relative edge length errors are quite stable and small. Furthermore, we also investigate some other surface models (as shown in Fig. 5) to conclude that 8D space is a good target dimension for computing SIFHDE with smaller relative edge length errors, e.g., Nefertiti surface: L rel avg = 2.98% and L rel max = 46.33%, Kitten surface: L rel avg = 3.82% and L rel max = 139.21%.

Volume Case: As for the volume models, we also investigate the best target dimension for computing the self-intersection free embedding with smaller relative edge length errors. Similar to surface case, volume models (as shown in Fig. 6) also conclude that 8D is a good target dimension for computing SIFHDE with smaller relative edge length errors, e.g., Cube volume: L rel avg = 1.65% and L rel max = 17.54%, Sphere volume: L rel avg = 3.21% and L rel max = 65.01%. As in the surface case, we also observe that beyond 8D, the relative edge length errors remain stable and small.

K-NN Efficiency on High-D Embedding

In order to compare the efficiency of the K-NN between the high-d embedding space and the original space, we use the same search radius 5σ, where σ is the kernel width of the inter-particle energy in the Euclidean embedding space 21 (as mentioned in Sec. 4.1). However, when the K-NN computation happens in the original space, we need to transform this search radius into the original anisotropic space, resulting in a large range once the anisotropic stretching ratio is high. After that, we need to further check and prune the spurious neighbors under the given metric. Based on the optimized particles on the self-intersection free 8D embedding, we can compute the RVD and its dual mesh, and finally map them back to the original 3D space to generate the anisotropic RVD and its dual anisotropic mesh. The results of the Cyclide model are shown in Fig. 8.

Besides Cyclide model, we also measure the K-NN efficiency on the high-d embeddings of other models, such as Kitten, Vase, Knot, Club, etc., and the average K-NN is quite stable and small, i.e., N avg = 20. However, without using Euclidean embedding, the average K-NN highly depends on the stretching ratio in the Riemannian metric. The higher the stretching ratio is, the larger the average K-NN is, such as N avg = 232 on another Cyclide model (as shown in Fig. 12) with the larger stretching ratio s2 s1 ∈ [2, 29].

2D RVD and Meshing

A 2D domain can be considered as a 3D flat surface, so 8D is a good dimension for the embedding computation. Besides the previous 2D domain with Gaussian embedding metric in Fig. 2, we also evaluate our method on the largely varying anisotropic field. Fig. 9 shows the 8D embedding errors (L rel avg = 4.04% and L rel max = 76.79%), the anisotropic RVD, and meshing results of a 2D square domain with 4000 output samples, with a highly-varying metric field M(x) = R(x) T diag{Stretch(x) 2 , 1}R(x), where the stretching field Stretch(x) ∈ [1, 40] with the rotation field R(x), which is defined by the Hessian of the non-convex analytic functions u(x, y) = tanh(10(sin(5y) -2x)) + x 2 y + y 3 .

3D Surface RVD and Meshing

Our Method Results: Fig. 10 shows the anisotropic 3D surface RVD and meshing results of Kitten, Vase, Knot, Club, Rocker Arm, and Hand models with the anisotropic metrics designed by the surfaces' curvature tensors. The statistics and timings for our 8D embedding computation and surface meshing are shown in Tab. 2. It is noted that our embedding computation in 8D space is quite efficient, and all surface models provided in this article need only dozens of seconds. Comparison with Other Embedding Methods: Fig. 11 shows the comparison with 2D conformal embedding method [START_REF] Zhong | Anisotropic surface meshing with conformal embedding[END_REF] and Lévy and Bonneel's 6D embedding [START_REF] Lévy | Variational anisotropic surface meshing with Voronoi parallel linear enumeration[END_REF], and our method, on a Kitten surface model with 5000 output samples. The 2D conformal embedding method produces poor meshing result, since it computes the embedded surface into 2D space, which has very limited degrees of freedom. However, the 6D embedding has more degrees of freedom, but it applies a 6D metric in terms of vertex positions and normals, without considering the curvature tensors, so its meshing result is even worse than 2D conformal embedding. Our high-d embedding method, i.e., 8D, obtains better mesh quality result and it also better matches the input curvature-based metric field. We also use Hausdorff distance to measure the faithfulness accuracy to the original input surface of these three methods. The Hausdorff distance of our result is 0.004964% of the bounding-box's diagonal length, while the Hausdorff distances of 2D conformal embedding and 6D embedding results are 0.013415% and 0.005110%, respectively. It is noted that the meshing fidelity of 2D embedding is worse than both those of 6D and our 8D embeddings, though it has a reasonable mesh quality result to match the designed metric. Essentially, when using dimension reduction strategy, such as mapping a 3D shape into a 2D domain, it may lose some 3D shape fidelity information. We also try to compare with the 3D embedding method [START_REF] Panozzo | Frame fields: Anisotropic and non-orthogonal cross fields[END_REF], but there are 1001 self-intersecting faces in the Kitten model as shown in Tab. 1. Then, if we would like to compute RVD and meshing on it, some parameterization methods are needed, which is beyond the scope of this article. Due to this limitation, we do not provide the quantitative comparison between the 3D embedding and our method for meshing computation.

Comparison with Other Anisotropic Surface Meshing Methods: In this subsection, we show further comparative analysis and experiments with other state-of-the-art anisotropic meshing approaches. Fig. 12 compares our method with both particle-based method [ZGW + 13] and LCT method [START_REF] Fu | Anisotropic simplicial meshing using local convex functions[END_REF] on a Cyclide model with the large stretching ratio s2 s1 ∈ [2, 29]. Our method improves the mesh quality significantly as shown in Fig. 12 compared with the original particle-based method [ZGW + 13], since we have computed the particle optimization directly on the high-d embedding. Compared with LCT method (one of the optimal high-quality anisotropic meshing methods), our method 13 visualizes the final mesh quality errors (i.e., 1-G, where G is the quality of triangle as mentioned in Sec. 5.3) of the above three methods. We can see that our errors are more smoothly distributed, whereas both particle-based method and LCT method have more poor quality triangles.

3D Volume RVD and Meshing

The anisotropic metric fields in 3D volume models are designed by different analytic functions in our experiments. The statistics and timings for our 8D embedding computation and volume meshing are shown in Tab. 3. Fig. 14 shows the anisotropic 3D volume RVD and meshing results of a Cube model with domain [0, 1] 3 of two different Riemannian metric fields M(x) = diag{Stretch(x) 2 , 1, 1}. The first one is designed by a linear function with stretching ratios Stretch(x) = 10x, where Stretch(x) ∈ [1, 10], and the second one is designed by a highly nonlinear function with large stretching ratios Stretch(x) = 0.2(0.0025 + 0.2(1e -|x-0.5| )) -0.8 , ∈ [1,25].

A unit Sphere model with a cylindrical anisotropic metric field is given in Fig. 1. The Riemannian metric is designed as M(x)=R(x) T diag{Stretch(x) 2 , 1, 1}R(x), where Stretch(x) ∈ [1, 5], and R(x)'s three columns are (x/ x 2 + y 2 , y/ x 2 + y 2 , 0) T , (-y/ x 2 + y 2 , x/ x 2 + y 2 , 0) T , and (0, 0, 1) T . Fig. 15 shows one volumetric RVD result on a Cube model with domain [1, 11] 3 . The targeted cylindrical Riemannian metric field is specified via a highly nonlinear analytic function as M(x) = R(x) T diag{Stretch(x) 2 , 1, 1}R(x), where Stretch(x) = 2(0.1 + 2(1 -e -0.01|x 2 +y 2 -49| )) -1 , and R(x) is a cylindrical rotation field as defined previously. The stretching ratio is Stretch(x) ∈ [1, 20]. By using our method, we can generate highly curved 3D Voronoi cells in regions where the metric varies quickly, while in low stretching regions, the regular 3D Voronoi cells are obtained. Its dual tetrahedral mesh is given in the following section.

Fig. 16 shows the meshing comparison with our method and LCT method [FLSG14] on a Cube model with the previous cylindrical Riemannian metric field. By using almost the same number of vertices, i.e., 6338 vertices in LCT method and 6300 vertices in our method, we achieve better mesh quality results in average value measurements as shown in Fig. 16. Visually, our result achieves better mesh stretching ratios and directions in both surface and interior of the volume. The reason why we have worse results in minimal mesh quality and smallest angle is we do not apply any sliver elimination strategy in our meshing method. Since the sliver elimination is a post-processing in both methods, it is more straightforward to compare the performance of two algorithms without doing

Discussion and Future Work

To the best of our knowledge, this is the first article in the literature for computing the self-intersection free Euclidean embedding in arbitrary dimensions and using it in Voronoi Diagram, surfaces and volume meshing equipped with Riemannian metrics. In particular, this is the first practical algorithm for computing volumetric anisotropic Voronoi diagrams. While there is a direct application to anisotropic meshing, we feel that there are many unexplored applications as limitations and future work: (1) We have not yet explored the cases where the metric has sharp discontinuities, in such case, finding a practical embedding may be challenging. (2) In this work, we would like to emphasize that all the costs of computations in a high-d space are justified since we can obtain high-quality anisotropic RVD and meshes as compared with previous methods. In the future, we will improve the computational speed by using GPU-based parallel algorithm and implementation. (3) In order to generate the sliver-free tetrahedral meshes, we will apply some sliver suppressing strategies, such as perturbations [START_REF] Tournois | Perturbing slivers in 3D Delaunay meshes[END_REF], sliver elimination [START_REF] Fu | Anisotropic simplicial meshing using local convex functions[END_REF] or gradient-based shape matching [NZL + 17]. (4) Besides anisotropic RVD and meshing, we will explore other important and interesting applications by using the proposed SIFHDE.

  al. [DMS14, DSPS15] used Lévy and Bonneel's framework to compute remeshing, but instead of optimizing the CVT in 6D, they used local operations (i.e., edge flips). It should be noted that both Lévy and Bonneel's and Dassi et al.'s work does not provide user's flexibility to control the anisotropy through an input arbitrary metric field. Particle-Based Anisotropic Meshing: Particle-based approaches for anisotropic meshing have been proposed during the past two decades. Bossen and Heckbert [BH96] simulated the repulsion and attraction forces between particles based on a distance function with the metric tensor. Shimada et al. [SG95, SYI97] physically modeled the inter-bubble forces by a bounded cubic function of the distance, and further extended it to anisotropic meshing by converting spherical bubbles to ellipsoidal ones. Both Bossen et al. and Shimada et al.'s work requires dynamic population control schemes, that is to adaptively insert or delete particles / bubbles in certain regions. Zhong et al. [ZGW + 13]

  Self-intersection free high-d embedding optimization (d) RVD and its dual mesh on high-d embedded surface (e) Mapping the RVD and its dual mesh to the original anisotropic metric domain (a) Input: 2D domain with anisotropic metric errors of 3D embedding at 20 iterations (c) Uniform particle distribution on high-d embedded surface G min = 0.81 G avg = 0.98 θ min = 40.

Figure 2 :

 2 Figure 2: Illustration of the proposed high-d embedding for anisotropic RVD and surface meshing. (a) A simple example of a 2D domain equipped with anisotropic metric, which has an ideal 3D embedding -the Gaussian surface as the ground truth. (b) Our embedding optimization iteratively deforms the original domain into the embedded surface. After embedding we can (c) uniformly and isotropically sample particles and (d) compute the isotropic RVD and its dual mesh on the embedded surface. (e) When they are mapped to the original 2D domain, the RVD and the triangular mesh will exhibit the desired anisotropy. triangular / tetrahedral meshing.Since the high-d Euclidean embedding without self-intersections has been computed in the previous section (Sec. 3), in the following section, under the isotropic metric space, we use the particle optimization technique (Sec. 4.1) to compute isotropic RVD and its dual Delaunay triangulation on the high-d embedded surface or volume (Sec. 4.2), due to its efficiency to generate regular hexagonal patterns of particles on 2-manifold and regular body-centered-cubic (BCC) lattices on 3-manifold, which are similar to the results of CVT [DFG99, LWL + 09, DW05b]. Moreover, the particle-based method is easily formulated in the high-d space. Finally the generated RVD and triangulation / tetrahedralization can be mapped from the high-d embedding space to the original space, on which the RVD and triangulation / tetrahedralization will exhibit the desired anisotropic property.

Data:

  Original surface or volume Ω with user-specified metric M and target dimension m Result: Vertex coordinates {v i |i = 1, • • • , n v } of the embedded surface or volume Ω in m space Initialize vertex coordinates in m space; while stopping condition not satisfied do for each triangle or tetrahedron j do Compute U j by fixing {v i |i = 1, • • • , n v } and using Eq. (19) or Eq. (20); end Compute {v i |i = 1, • • • , n v } by fixing all U j s and fixing v 1 , i.e., solve a linear system to minimize the regularized objective function E total of Eq. (18); end Algorithm 1: High-D Embedding Optimization Framework. between ours and Zhong et al.'s method [ZGW + 13]. The limitations of their methods are: (1) they made use of the concept of embedding space for their energy and force formulations, without computing such high-d embedding explicitly, which may lead to more approximation error;

Figure 3 :

 3 Figure 3: A 3D embedding result of a Cyclide surface with the stretching ratio s2 s1 ∈ [1.6, 9.4]. There are 1146 self-intersecting faces out of total 21,600 faces as shown in green color.

Figure 4 :

 4 Figure 4: The color-coded relative edge length error images of high-d embeddings (i.e., 4D ∼ 20D) of the Cyclide surface.

Figure 5 :Figure 6 :

 56 Figure 5: The relative edge length errors of high-d embeddings (i.e., 4D ∼ 20D) of the Nefertiti and Kitten surfaces.

Fig. 7 Figure 7 :Figure 8 :

 778 Figure 7: Comparison of K-NN searching of 2000 particles optimization on the original anisotropic 3D Cyclide surface and its 8D Euclidean embedded surface.

  Stretch θ G min = 0.34 G avg = 0.89 θ min = 18.

Figure 9 :

 9 Figure 9: The embedding errors, anisotropic RVD, and anisotropic meshing results of a 2D domain with a highly-varying metric field M(x) = diag{Stretch(x) 2 , 1}, where Stretch(x) ∈ [1, 40].

Figure 10 :

 10 Figure 10: The anisotropic 3D surface RVD and meshes (Kitten, Vase, Knot, Club, Rocker Arm and Hand models) generated by our 8D embedding method.

G

  Figure 11: Comparison on anisotropic meshing results with 2D conformal embedding method [ZSJG14], Lévy and Bonneel's 6D embedding [LB12], and our method on the Kitten surface model.

Figure 12 :

 12 Figure 12: Comparison on anisotropic surface meshing results (8000 output samples) with particle-based method [ZGW + 13], LCT method [FLSG14], and our method on the Cyclide surface model with the stretching ratio s2 s1 ∈ [2, 29]. increases the minimal angle from 26.94 • to 29.07 • and % <30 • = 0.03% versus % <30 • = 0.009%, meanwhile we have similar quality of average values. Fig.13visualizes the final mesh quality errors (i.e., 1-G, where G is the quality of triangle as mentioned in Sec. 5.3) of the above three methods. We can see that our errors are more smoothly distributed, whereas both particle-based method and LCT method have more poor quality triangles.

Figure 13 :

 13 Figure 13: Comparison on final mesh quality errors (i.e., 1-G) with particlebased method [ZGW + 13], LCT method [FLSG14], and our method on the Cyclide surface model with the stretching ratio s2 s1 ∈ [2, 29].

Figure 14 :

 14 Figure 14: The anisotropic 3D volume RVD and meshes on two Cube models (with different Riemannian metric fields) generated by our 8D embedding method. Both the surface and interior of the volume RVD and tetrahedral meshes are shown. In order to better visualize the 3D RVD results, we shrink each Voronoi cell (i.e., polyhedron) a little bit. The blue color represents the outside of each 3D Voronoi cell, and the red color represents the inside of each 3D Voronoi cell.

Figure 15 :

 15 Figure 15: The anisotropic 3D volume RVD on the Cube model with the large cylindrical anisotropic metric field generated by our 8D embedding method. Both the surface and interior of the volume RVD are shown.

Figure 17 :

 17 Figure 17: The anisotropic 3D volume RVD on the Sphere model (with a cosine wave anisotropy) generated by our 8D embedding method. Both the surface and interior of the volume RVD are shown.

  6. Refinement-Based Anisotropic Delaunay Triangulation: Delaunay refinement-based approaches involve iterative point insertion in a Delaunay triangulation. Many researchers have extended these approaches to work with anisotropic triangulations, and there are some practical applications [BGH + 97, BGM97, DF08]. Cheng et al. [CDES01, CDRR04, CDR06] applied Delaunay refinement to anisotropic surface meshing. Boissonnat et al. [BWY08b,

Table 1 :

 1 Statistics of self-intersecting faces for embeddings in 3D and high-d spaces on different surface models.

	Model Cyclide1 Cyclide2 Kitten Vase Knot Club Rocker Arm Hand
	3D	1146	1751	1001 502 6552 4444	7545	616
	High-D	0	0	0	0	0	0	0	0

Table 2 :

 2 Statistics and timings for embedding computation and surface meshing.

	Model	In # V.	Stretch	L rel avg	L rel max	Embed. (s) Out # V. G min Gavg	θ min	θavg	% <30 • Mesh (s)
	Cyclide1 (Fig. 8)	10, 800 [1.6, 9.4] 1.83% 25.81%	9.14	2000	0.57	0.92	28.88 • 54.47 •	0.05%	6.22
	Cyclide2 (Fig. 12) 25, 920	[2, 29]	2.18% 19.39%	23.63	8000	0.63	0.92	29.07 • 53.70 • 0.009%	42.83
	Kitten	20, 000	[1, 6]	3.82% 139.21%	17.56	2000	0.33	0.85	20.69 • 48.60 •	0.34%	7.87
	Vase	20, 000	[1, 5]	4.18% 157.97%	17.09	2000	0.36	0.86	21.33 • 49.52 •	0.24%	7.48
	Knot	24, 392	[2, 8]	4.49% 168.46%	21.16	5000	0.47	0.93	25.37 • 53.79 •	0.05%	26.41
	Club	30, 000	[1, 6]	5.44% 143.44%	27.82	5000	0.46	0.91	24.10 • 51.25 •	0.07%	29.94
	Rocker Arm	35, 840	[1, 7]	5.95% 138.31%	34.77	5000	0.37	0.86	19.59 • 47.65 •	0.16%	33.39
	Hand	36, 619	[1, 5]	2.68% 121.57%	35.78	10, 000	0.28	0.90	16.28 • 50.14 •	0.23%	55.50
	Note: Embed. Time: timing for embedding computation with 50 iterations. Mesh Time: timing for both
		particle optimization with 50 iterations and RVD computation.

Our method LCT

sliver elimination. As mentioned in [START_REF] Fu | Anisotropic simplicial meshing using local convex functions[END_REF], they have 517 slivers (i.e., tetrahedrons with the dihedral angle < 15 • ) before applying sliver elimination, and we have a smaller number of slivers, i.e., 311. It is noted that the LCT result in Fig. 16 is after applying sliver elimination, however, our result is not.

Lastly, we show one more volumetric RVD result on a unit Sphere model with a smaller number of samples, i.e., 1000, in Fig. 17. The Riemannian metric is defined as a cosine wave anisotropy:

where Stretch(x) = 10, and R(x)'s three columns are (2cos(4x), 1, 0) T , (1, -2cos(4x), 0) T , and (0, 0, 1) T . A closeup shows clearly how the 3D Voronoi cells are extremely curved. Both the surface and interior of the volumetric RVD are illustrated. This example is to demonstrate the advantage of using our proposed SIFHDE as a tool to compute arbitrary anisotropic Voronoi diagrams. Without using SIFHDE, such as Budninskiy et al. [BLdG + 16]'s work by using a convex function to represent the anisotropy, it is very difficult to obtain our result.

Through the above experiments, it is noted that our proposed SIFHDE and high-d RVD are effective computational algorithms to generate the high-quality anisotropic 3D Voronoi diagrams and tetrahedral meshes with arbitrary smooth Riemannian metric fields.


Leonid V Averyanov

B B Thinh

Nguyen Thi Lieu

L H Truong

Quach Van Hoi
email: quachvanhoi@gmail.com

Luong Van Dung

Roman V Doudkin

Bao Bui

Thinh

& Luu

Hong Truong

Camellia proensis (Theaceae, sect. Piquetia), a new species from Southern Vietnam

Keywords: Camellia proensis, Lam Dong, taxonomy, Theaceae, Vietnam

  L'archive ouverte pluridisciplinaire

Introduction

Camellia L. is the largest genus in the Theaceae family, with about 280 species [START_REF] Chang | Camellias[END_REF][START_REF] Ming | Theaceae. In Flora of China[END_REF]. Distribution of this species ranges from Nepal, Bhutan, Northeastern India, Myanmar, China, Japan, to Southeast Asia [START_REF] Ho | Theaceae: In Cay co Viet Nam [An illustrated flora of Vietnam[END_REF][START_REF] Ninh | Biodiversity of the genus Camellia of Vietnam[END_REF][START_REF] Ming | Theaceae. In Flora of China[END_REF]. According to [START_REF] Ho | Theaceae: In Cay co Viet Nam [An illustrated flora of Vietnam[END_REF], about 30 species of Camellia were recorded in Vietnam. Later, Ninh's statistics showed that there were 50 species [START_REF] Ninh | Biodiversity of the genus Camellia of Vietnam[END_REF]. So far, more than 95 species of Camellia have been reported in Vietnam (Le et al. 2020). Vietnam is the center of diversity for Camellia species. Recently, many more new species of Camellia were discovered and described in Vietnam [START_REF] Orel | A new species of Camellia section Piquetia (Theaceae) from Vietnam[END_REF][START_REF] Tran | Camellia dalatensis: a new species and precious gene should be conserved[END_REF][START_REF] Ninh | Camellia dilinhensis: A new yellow species from Vietnam[END_REF][START_REF] Orel | Four new species and two new sections of Camellia (Theaceae) from Vietnam[END_REF][START_REF] Dung | Camellia thuongiana -a new yellow Camellia species from Vietnam[END_REF][START_REF] Truong | A new species of the family Theaceae from central Vietnam[END_REF][START_REF] Nguyen | Two new species of Camellia (Theaceae) from Vietnam[END_REF][START_REF] Pham | Camellia velutina (Theaceae, sect. Chrysantha), a new species from northern Vietnam[END_REF][START_REF] Do | A new yellow Camellia (Theaceae) from central Vietnam[END_REF], Nguyen et al. 2020, Do et al. 2020).

In a field trip to Pro Commune, Don Duong District, Lam Dong Province, Vietnam, in 2015, the authors discovered a Camellia species with yellow flowers but did not see the fruit. After many field trips, we collected the flowers and fruits. We then compared the morphological characteristics with previously published species and conducted a thorough review of specimens from the following herbaria: HN, VNM, and DLU (herbarium codes follow [START_REF] Thiers | Index Herbariorum: A global directory of public herbaria and associated staff[END_REF]). We also used digital images on the web from the Muséum National d'Histoire Naturelle (https://science.mnhn.fr/), Chinese Virtual Herbarium (http://www.cvh.ac.cn/) and Jstor Global Plant (https://plants.jstor.org/). After an exhaustive review, we were unable to find any matches to this Camellia. Particular attention compared the collected samples with other closely related species. We will name this new species Camellia proensis, based on the classification system of [START_REF] Chang | Camellias[END_REF] and confirm its placement in Section Piquetia (Pierre 1877[START_REF] Sealy | A revision of the genus Camellia[END_REF].

Taxonomic description

Camellia proensis L.V. Dung,Doudkin & Q.V.Hoi,sp. nov. (Figs. 1,[START_REF]Flowers orange or yellow[END_REF] Type:-VIETNAM. Lam Dong Province, Don Duong District, Pro commune, 800-1100 m elev., 8 March 2015, Quach Van Hoi, Luong Van Dung DUNG.150301 (holotype, isotypes, DLU!).

Description:-Small tree, 5-6 m tall. Young branches glabrous, light brown buds. Leaves oblong-elliptic, apex acuminate, base rounded or acute, adaxial and abaxial surfaces glabrous, 7.5-8.5 cm wide, 33-35 cm long; petioles 1-2 cm long, glabrous; lateral veins almost parallel, midrib and lateral veins clearly show in the abaxial surfaces, margin denticulate. Flowers 1-3 in the axillary, flower buds globose, 1.3-1.7 cm wide, flower 4.0-5.5 cm diameter, peduncle branched or not, 3.0-5.5 cm long, glabrous. Bractlets 2-3, light green, triangular, dropping quickly before the fruiting season, glabrous, 1-2 mm wide 1.5-3.0 mm tall. Sepals 5, hemisphere, persistent, light green, pubescent on the outside, glabrous on the inside, 6-12 mm tall, and 9-15 mm wide. Petals 5-6, almost round to obovate, light yellow, pubescent on the outside, glabrous on the inside, 1.7-1.8 wide and 1.7-2.4 cm tall. Stamens numerous, in 4-5 whorls, 1.5-2.0 cm long, outer filaments basally connate into a tube, filament tube 5-7 mm tall and adnate to the base of petals, pubescent at the base of inner filaments. Ovary globose, tomentose, 3.5-4.0 mm tall and 4-5 mm wide, has grooves, the ovary has 4-5 locular. Styles ripped 4-5 to base, tomentose to ½ from the base, 1.5-1.8 cm long. Fruit flat-globose, 4.0-5.5 cm tall and 8.5-10.0 cm wide, sparse hair, capsule five locular, 1-locular with two seeds, pericarp 1-1.3 cm thick. Seeds spherical, hemispherical or wedge-shaped, glabrous, 1.6-2.0 cm wide, 2.3-2.6 cm long.

Affinity:-Camellia proensisis morphologically similar to C. sonthaiensis in glabrous young branches and leaves; sepals 5, hairy on the outside; petal light yellow, hairy on the outside, from which it is distinguished by the size of mature fruit, 3.0-3.5 cm in diameter, 1.5-2.0 cm deep (C. sonthaiensis), 4.0-5.5 cm tall, 8.5-10.0 cm wide (C. proensis).

Etymology:-The species is named after the location where the authors first found this plant. Pro is the name of Pro commune, Don Duong District, Lam Dong Province, Vietnam.

Habitat:-This species is distributed in the evergreen broadleaf forest, the foot of the mountain where the slope is low, along the sides of the stream with relatively high humidity, or near the edge of the forest, where there is plenty of light. The species is distributed at elevations 800-1100 m, Pro commune, Don Duong District, Ta Nang Commune, Duc Trong District, Lam Dong Province.

Phenology:-These species produce buds and start to flower from January to March. Flowers bloom in late February, bearing fruits from August to October.

Note:-Camellia proensis has characteristics flowers 1-3 on a short bracteate shoot in axils of old leaves; flowers pedicellate, nodding. Bractlets 2-3, dropping quickly before the fruiting season; sepals 5, persistent; petals 5-6, inter petals united with outermost filaments. Gynoecium superior, ovary 4-5 locular, tomentose; styles 4-5, free. The new species has many characteristics of the section Piquetia (Pierre) Sealy as flowers are borne on a short shoots in the leaves' axils, usually, 3-5 flowers on a shoot but sometimes one only. Flowers pedicellate, nodding, pedicel very stout, thickened upwards; bracteoles 2-3, persistent; sepals 5, persistent; petals eight or more; stamen free above the union with the petals, puberulous inside; gynoecium densely hairy; styles 5 or 6 (Pierre 1877[START_REF] Sealy | A revision of the genus Camellia[END_REF]. Therefore Camellia proensis is classified into section Piquetia. Currently, the species in the sect. Piquetia

  occurring in Vietnam are: C. dalatensis,C. dongnaiensis, C. honbaensis, C. longii, C. piquetiana and C. sonthaiensis (Pierre 1877[START_REF] Sealy | A revision of the genus Camellia[END_REF][START_REF] Orel | A new species of Camellia section Piquetia (Theaceae) from Vietnam[END_REF][START_REF] Tran | Camellia dalatensis: a new species and precious gene should be conserved[END_REF][START_REF] Orel | Four new species and two new sections of Camellia (Theaceae) from Vietnam[END_REF][START_REF] Truong | A new species of the family Theaceae from central Vietnam[END_REF]. The new taxon is easily distinguishable from C. honbaensis, C. piquetiana and C. longii by the yellow flowers; from C. dalatensis by the glabrous twigs and leaves; from C. dongnaiensis by light yellow flowers; from C. sonthaiensis by mature fruit flat-globose, 4.0-5.5 cm tall and 8.5-10.0 cm wide.

FIGURE 1 .

 1 FIGURE 1. Camellia proensis. A, B. Young twig. C, E. Flowering shoot and arrangement of petals. D. Open flower. F, G. Fruits. H. Seeds. Photos by Quach Van Hoi.

FIGURE 2 .

 2 FIGURE 2. Camellia proensis. A. Leaf. B. Venation of leaf (adaxial surface). C. Venation of leaf (abaxial surface). D. Flower bud. E, F. Flower. G. Sepals. H. Petals. I. Stamens. K. Styles and gynoecium. L. Fruit. M. Seed. Drawings by Luong Van Dung.

Acknowledgments

We would like to thank Stephen Maciejewski, The Gesneriad Society/American Begonia Society, and Michael LoFurno, Adjunct Professor, Temple University, Philadelphia PA, USA, for their editorial assistance.