Quasi-Isometric Mesh Parameterization Using Heat-Based Geodesics and Poisson Surface Fills

In the context of CAD, CAM, CAE, and reverse engineering, the problem of mesh parameterization is a central process. Mesh parameterization implies the computation of a bijective map φ from the original mesh M ∈ R3 to the planar domain φ(M) ∈ R2. The mapping may preserve angles, areas, or distances. Distance-preserving parameterizations (i.e., isometries) are obviously attractive. However, geodesic-based isometries present limitations when the mesh has concave or disconnected boundary (i.e., holes). Recent advances in computing geodesic maps using the heat equation in 2-manifolds motivate us to revisit mesh parameterization with geodesic maps. We devise a Poisson surface underlying, extending, and filling the holes of the mesh M. We compute a near-isometric mapping for quasi-developable meshes by using geodesic maps based on heat propagation. Our method: (1) Precomputes a set of temperature maps (heat kernels) on the mesh; (2) estimates the geodesic distances along the piecewise linear surface by using the temperature maps; and (3) uses multidimensional scaling (MDS) to acquire the 2D coordinates that minimize the difference between geodesic distances on M and Euclidean distances on R2. This novel heat-geodesic parameterization is successfully tested with several concave and/or punctured surfaces, obtaining bijective low-distortion parameterizations. Failures are registered in nonsegmented, highly nondevelopable meshes (such as seam meshes). These cases are the goal of future endeavors.


Introduction
Mesh parameterization is the process by which a piecewise linear surface (i.e., triangular mesh) M is mapped with the least possible distortion onto a planar (R 2 ) region, via a bijective continuous function φ : M → R 2 . The mesh M is supposed to be a connected 2-manifold with border (and possibly holes).
Mesh parameterization is central in tool path generation for surface machining, texture mapping, thermoforming of thin layers (metal, leather, plastic, etc.), reverse engineering, finite element remeshing, facial expressions, morphing, etc.
A geodesic curve between two points of a continuous surface is the shortest path within the surface joining the two points. The length of such a path is the geodesic (shortest) distance, embedded in the surface, between those two points.
Given any two points x i , x j ∈ M, ideal parameterizations of such a surface seek to map them to φ(x i ), φ(x j ) ∈ R 2 so that the geodesic distance between x i and x j in M matches the Euclidean 2D distance between φ(x i ) and φ(x j ) in R 2 . On the rare occasions that this goal is possible, M is called a developable surface and φ is an isometric map. When the distortion in such distances is small, one qualifies M as a quasi-developable surface. This case is sufficiently frequent, since large triangular meshes can be segmented with a goal being that the resulting submeshes are quasi-developable or developable.
It is not convenient, when the mesh has holes or concavities in its boundary, to parameterize the mesh via geodesics. The reason is that mesh points being close neighbors in the surface may be far apart via geodesic curves due to mesh gaps or concavities.
Mesh parameterization algorithms can be classified depending on the type of distortion being minimized, as follows: (a) Area-preserving (authalic) algorithms; (b) angle-preserving (conformal) algorithms; and (c) distance-preserving (isometric) algorithms. The remainder of this section discusses a summary of recent mesh parameterization algorithms already present in the literature (Detailed surveys on mesh parameterization algorithms are presented in [1][2][3]).

Area-Preserving Mesh Parameterization
Area-preserving (authalic) parameterization algorithms rely on the minimization of an area preserving continuous cost function. Zou et al. [4] solved a Lie advection problem on the mesh M. The gradient of the scalar Lie advection field was then added to an initial parameterization φ 0 of M, resulting in an authalic parameterization. Zhao et al. [5] solved an optimal mass transport problem from the mesh M to its parameterization φ(M). The optimal mass transport poses a partial differential equation in which the parameterization φ(x) locally preserves the area at every point x ∈ M. Since most optimal transport methods only parameterize meshes with a connected boundary (i.e., without holes), Su et al. [6] introduced additional boundary conditions to allow authalic parameterizations of meshes with more complex topologies.

Angle-Preserving Mesh Parameterization
Angle-preserving (conformal) optimization aims to minimize the parameterization angle distortion. Since this objective can be achieved by collapsing all triangles to a single point, these algorithms rely on constraining the parameterized boundary to a region in R 2 . Disk geometries are usually used in this context [7,8] however, other geometries such as squared domains have also been proposed [9,10]. The imposed boundary restrictions in these constrained optimization algorithms induce additional distortion in the resulting parameterization.
Free boundary algorithms allow unrestricted boundary parameterizations, producing less distorted maps. Sawhney and Crane [11] presented an algorithm in which the mesh boundary is mapped to R 2 according to its shape. The parameterized boundary is then used as a constraint to produce a boundary-free parameterization. Starting from a disk parameterization, Bright et al. [12] performed nonlinear optimization while unconstraining the boundary edges of the parameterized mesh. The resulting parameterization allows the (initially mapped to disk) boundary to move freely in the parameter space. Smith and Schaefer [13] presented a mesh parameterization method which introduced a barrier function in its optimization process. The introduced barrier function penalizes nonadjacent triangle overlaps, which guarantees global bijectiveness in the resulting parameterization.
Dimensionality reduction is a superset of mesh parameterization, in which a d-manifold embedded in a higher dimensional space R D , is parameterized to its corresponding R d domain. As a consequence, these algorithms have been applied successfully in mesh parameterization applications. Such algorithms include Laplacian Eigenmaps [14] and Hessian Locally Linear Embedding [15].

Distance-Preserving Mesh Parameterization
By definition, a distance-preserving (isometric) mapping is a function that simultaneously preserves areas and angles. Mejia et al. [16] presented a nonlinear minimization algorithm for area-angle (isometry) preservation. The minimization function is a linear combination of area and angle distortion, and the weighting parameters for each distortion term are adjusted by the user. The authors pointed out that the algorithm performs better when the angle-preserving term is preponderant over the area-preserving one. Similarly, Yu et al. [17] used polar factorization to introduce area-angle preserving mappings, in which area distortion increases as angle distortion decreases. ARAP (As Rigid As Possible) algorithms divide the parameterization into two optimization steps-local parameterization and global parameterization-performing these steps iteratively one after another until convergence [18][19][20]. These algorithms produce different bijective parameterizations, but since the weighting parameters are nonoptimized (as they are user-defined), the resulting parameterization is rarely the optimal distance-preserving map.
Ruiz et al. [14] used a dimensionality-reduction geodesic-based algorithm (Isomap) for the computation of isometric parameterization of quasi-developable meshes. However, in addition to the classic nonconvex parameterization problems, such algorithms estimate geodesics using shortest-path graph algorithms which introduce additional distortions in the resulting parameterization. Li et al. [21] presented a geodesic approximation approach in which cutting planes are intersected with the mesh to estimate geodesic paths. This approach solves the problem of nonconvex surfaces and distortion errors induced by mesh graph approximation. However, the method is limited to geodesic curves embedded in R 2 (i.e., the cutting plane).

Conclusions of the Literature Review
Most of the distance-preserving parameterization algorithms rely on weighting angle vs. area distortion. Such a weighting is defined by the user and drastically changes the resulting parameterization, not providing the optimal isometric mapping. Geodesic-based parameterization algorithms solve this problem by directly minimizing the distance distortion. However, current geodesic-based algorithms rely on shortest-path graph algorithms for geodesics estimation, introducing unnecessary distortion in the resulting parameterization. In addition, estimation of geodesics fails when the surface is nonconvex (such as surfaces with holes and boundary concavities).
To overcome these problems, this manuscript presents a novel heat-geodesic mesh parameterization algorithm. Our algorithm computes a set of temperature maps (heat kernels) on the mesh M, which are then used to retrieve the set of point-to-point geodesic distances on the discrete mesh. Afterwards, a near-isometric parameterization is obtained by minimizing the difference between the parameterization Euclidean distances and their corresponding geodesics. Since our method relies on finite element mesh discretization to estimate the temperature maps and geodesics, our geodesics estimation is unaffected by mesh-graph topology (as opposed to shortest-path graph algorithms). To overcome the nonconvexity problem, our algorithm uses Poisson surface reconstruction [22], in which the surface holes and boundary concavities are temporarily filled for parameterization. The resulting parameterization for the Poisson reconstructed surface is trimmed with the original boundary of M, producing a trimmed surface. The implementation and integration of these different techniques provide a novel geodesic-based mesh parameterization algorithm which is (1) unaffected by mesh holes and/or concavities and (2) less sensitive to mesh graph topology.

Methodology
Consider M = (X, T ) (with point set X = {x 1 , x 2 , · · · , x n } and triangle set T = {t 1 , t 2 , · · · , t m }), a connected 2-manifold with border (and possibly holes) embedded in R 3 . The problem of mesh parameterization consists of finding a set of points Φ = {φ 1 , φ 2 , · · · , φ n } ⊂ R 2 such that φ i is the image of x i ∈ M under the image of a homeomorphism φ : M → R 2 (i.e., φ i = φ(x i )). The function φ must satisfy the following conditions: 1. Continuity: If t i ∈ T and t j ∈ T (t i = t j ) are adjacent triangles in M, then φ(t i ) and φ(t j ) are adjacent in φ(M). 2. Local bijectiveness: All mapped triangles φ(M) share the same orientation in R 2 . 3. Global bijectiveness: Triangles in φ(M) do not overlap each other. This can happen even if all triangles share the same orientation as illustrated in [13].
In addition, define g : , then φ is an isometric mapping (i.e., φ preserves geodesic distances) and M is a developable surface.
As most of the surfaces are not developable in practice, we aim to find the discrete mapping Φ that minimizes the difference between these two distances as follows: where the restriction ∑ n i=1 φ i = 0 indicates that the parameterization is mean centered, i.e., the center of mass of the parameterization points is 0 ∈ R 2 . Such a restriction is introduced to obviate all the possible translations of the same solution.

Algorithm Overview
Our mesh parameterization algorithm aims to retrieve a parameterization φ(M) which highly preserves the geodesic distances of M as Euclidean distances. In order to estimate the geodesic distances g in M, the heat-based geodesics algorithm presented in [23] is implemented. Afterwards, we use classic multidimensional scaling to retrieve the 2D coordinates of the parameterization φ from the computed geodesic distances. In the case that M presents holes or concavities, our algorithm applies Poisson surface reconstruction [22] and computes a parameterization on a trimmed surface instead. A summary of the algorithm is presented in Figure 1.
The remainder of this section details the steps to solve Equation (1), and the details of our mesh parameterization algorithm. The algorithm was implemented in MATLAB R [24], except for the Poisson Surface Fills which were implemented in C++ using the PCL library [25]. Figures including triangular meshes, scalar fields, vector fields, and 2D parameterizations were produced in MATLAB R . Figures including texture maps on the surface were produced using MeshLab [26].

Mesh Heat Kernels
A heat kernel of a point x i ∈ M is a function u i : M × (0, T] → R that satisfies the following partial differential equation [27]: where ∆ is the Laplace-Beltrami operator, u i is the temperature distribution (heat kernel) associated to the source point x i , x ∈ M, t ∈ (0, T] are the spatial and time coordinates, respectively, and T > 0 is the simulation time. The term ∂u(x,t) ∂n ∂M = 0 corresponds to the Neumann boundary condition (no boundary heat flux). Finally, the term u i (x, 0) = δ x i (x) corresponds to Dirac delta initial conditions, i.e., The above initial conditions dictate initial infinite temperature at vertex x i and 0 everywhere else. After some timem t > 0 has passed, heat dissipates from x i as illustrated in Figure 2. Equation (2) can be solved using a Finite Element discretization scheme, as follows:

Source point
where ∆t is the time step, U t)), and L and B are the n × n Laplace-Beltrami and mass (sparse and symmetric) matrices, respectively. For a given edge e ij = (x i , x j ), the Laplace-Beltrami matrix L is defined as follows: where α ij , β ij ∈ (0, π) are the two opposite angles to the edge e ij , and S * i = {k|e ik ∈ Edges(M)} is the index set of all incident edges to the vertex x i ∈ X. The entries L ij of the matrix L are known as cotangent weights [28].
Similarly, the mass matrix B is defined as follows: where t 1 , t 2 ∈ T are the pair of triangles adjacent to the edge e ij and |t l | is the area of the triangle t l (l = 1, 2). For each vertex x i ∈ M, Equation (4) is solved for U i using a sparse Cholesky linear solver. It is worth noting that for every x i and t ∈ (0, T], the matrices ∆tL and B are the same. As a consequence, these matrices are prefactored only once using Cholesky factorization, which speeds up the computation of the heat kernels.
Finally, the simulation parameters ∆t, T are chosen according to [23]: with ∆t computed as the magnitude of the largest edge in M, and T is equal to ∆t. These values have provided better results in our parameterization experiments than other values.

Heat-Based Geodesic Distance
The vector field −∇u i (x, t) (∇: gradient operator on M) describes the heat flux on M for the respective heat kernel u i . The normalized heat flux vector field H i is defined as follows: It is worth noting that the magnitude of the vector field H i is 1 everywhere ( H i (x, t) = 1). In addition, as illustrated in Figure 3, the temperature contours are perpendicular to the geodesic paths from x i ∈ M, and the corresponding vector field points in the same direction as such paths.
The geodesic field g(x i , x j ) satisfies the following differential equation [23]: where ∇ · H i (x, T) is the divergence field of the normalized heat flux. Similar to Equation (2), Equation (9) is discretized using the same Finite Element scheme, as follows: LG where G i = {g i1 , g i2 , · · · , g in } is the vector of geodesic distances g ij = g(x i , x j ) and K in } is the divergence field of the normalized gradient k Figure 4 plots the estimated geodesic distance field g(x i , x) for the vertex x i .

Multidimensional Scaling (MDS)
After the geodesic field g ij = g(x i , x j ) has been estimated on M, the minimization problem in Equation (1) can be solved. Classic multidimensional scaling poses an equivalent minimization problem [29]: Let C be the symmetric, semidefinite positive matrix whose entries contain the mean centered squared geodesics (i.e., C ij = − 1 2 [g 2 ij − 1 n (∑ kl g 2 kl )]). In matrix form, C is equivalent to where I n , J n are the n × n identity and all-ones matrices, respectively, and G 2 is the n × n symmetric matrix whose entries contain the squared geodesic distances in M, i.e., G 2 ij = g 2 ij . Then, Equation (11) becomes: ij as the (squared) Frobenius norm of A. Finally, Equation (13) can be solved by an eigendecomposition of C as follows [29]. Let λ 1 and λ 2 be the largest positive eigenvalues of C, with respective n × 1 eigenvectors V 1 and V 2 . The near-isometric parameterization of M becomes where √ λ 1 V 1 are the discrete u-coordinates and √ λ 2 V 2 are the discrete v-coordinates of the parameterization Φ = {φ 1 , φ 2 , . . . , φ n } ⊂ R 2 . Figure 5 plots the resulting parameterization using MDS on the estimated geodesic fields G.
from the estimated geodesic distances.

Poisson Mesh Reconstruction
In the case that M is nonconvex (i.e., it has holes or concavities), we seek to compute an underlying extending surface M * . Such a surface contains the points in M, and fills the holes and concavities by extending M in such areas (M ⊂ M * ). As an example, a geodesic path in a nonconvex M circles a given hole (Figure 6a). On the other hand, the same geodesic path in the extended surface M * crosses through the hole (Figure 6b). To compute the surface M * , our parameterization algorithm uses Poisson surface reconstruction [22] from the PCL library [25]. Define χ : R 3 → R as an indicator function such that χ(x) = 1 if x ∈ R 3 is "inside" the solid defined by M * ; and χ(x) = 0 if x is "outside" such solid. The surface M * is composed by the points in between, as follows [22]: The indicator function χ is computed by solving the following partial differential equation in R 3 [22]: where ∆ and ∇ are the R 3 Euclidean Laplacian and gradient operators, respectively. It is worth noting that these Laplacian and gradient operators are different from the 2-manifold version presented in Sections 2.2 and 2.3. N : R 3 → R 3 is a vector field in R 3 whose value N(x) is the normal vector to the original surface M if x ∈ M, and N(x) = 0 everywhere else. To solve Equation (16), the PCL library uses a hierarchical 3D spatial discretization and a Finite Differences approach [25].

Results
To test our mesh parameterization algorithm, we run tests with several parameterizable surfaces. Section 3.1 presents a comparison of our mesh parameterization algorithm without Poisson surface filling vs. Poisson surface filling, for quasi-developable nonconvex meshes. Section 3.2 presents parameterization results for some challenging, strongly nondevelopable data sets. Finally, Section 3.3 presents the application of our parameterization algorithm for the reverse engineering of a scanned cow vertebra. Figure 11 plots the parameterization results (2D Φ coordinates and 3D texture map) for the Mask data set. Without using Poisson surface reconstruction, the resulting parameterization is bijective with relatively low distortion, except for the eye holes (Figure 11a). However, such a parameterization is improved by applying the Poisson reconstruction, reducing the distortion close to mesh holes and boundary concavities (Figure 11b). A more extreme case is illustrated in Figure 12, with the S-trimmed-on-cone data set from [14]. Since the S shape is trimmed on a cone, M is a fully developable surface. Yet, the heat-geodesic parameterization on M fails to compute a bijective parameterization (Figure 12a

Parameterization of Strongly Nondevelopable Meshes
In this section, the public data sets Foot, Gargoyle, and Cow are parameterized with our heat-based geodesics algorithm. These benchmark datasets contain an artificially introduced border [30], which allows their parameterization. Figure 13a,b plot the parameterization results for the seam Foot and Gargoyle, respectively. The resulting parameterization is bijective, with some noticeable distortion (e.g., in the Gargoyle head). Figure 13c plots our parameterization results for the seamed Cow, which is not bijective in the head, legs, and tail. It is worth noting that although parameterizable, these benchmark data sets are strongly nondevelopable without a proper mesh presegmentation. This fact is illustrated in the next section.

Reverse Engineering of Cow Vertebra
For this section, the vertebra of a cow is scanned using a 3D optical scanner. The scanned mesh is closed and therefore accepts no (bijective) parameterization. The closed mesh is segmented into quasi-developable meshes using a heat-based segmentation approach [31]. Figure 14 plots the parameterization results for the segmented mesh. Each submesh bijective parameterization presents low distortion, enabling further reverse engineering operations such as NURBs reparameterization [14], finite element analysis, structural optimization, and/or dimensional inspection [31].

Conclusions
This article presents the implementation of a novel application of heat propagation in 2-manifolds used for mesh parameterization. The temperature contours for the heat kernels computed on the mesh are perpendicular at each point to the geodesic curves in the surface. This principle permits determination of geodesic maps in the mesh and in particular vertex-to-vertex geodesic distances. Although Finite Element methods produce better results as the mesh resolution increases (higher polygon count), the geodesics estimation method still produces robust results for low polygon count meshes as each geodesic path traverses across the mesh faces (contrary to graph algorithms which traverse the mesh graph). A quasi-isometric bijective function (i.e., the parameterization) is synthesized to map the 3D mesh to the parameter (2D) space. This parameterization is near isometric in the sense that geodesic distance on the mesh between any two points on the mesh approximates the Euclidean distance between their images in the parametric space. This approach is obviously limited to meshes that are quasi-developable or developable, since mesh developability is necessary for the existence of an isometric parameterization for it.
Our approach circumvents the weakness of geodesic maps in the presence of mesh interruptions (boundary concavities and mesh holes) by devising an underlying continuous Poisson surface that approximates the input mesh but contains no such interruptions. This underlying surface allows for geodesic maps to be computed on it, which are also valid in the input mesh. In this manner, the parameterization computed for the Poisson surface is valid for the input mesh. Finally, the boundaries of the input mesh are explicitly marked on the parameterization to obtain a trimmed surface or FACE in the Boundary Representation jargon.
Future work is required in these aspects-(a) to eliminate redundant computation that is present in the construction of heat-based geodesic maps and (b) to use failures in the bijectiveness of the computed parameterizations to force mesh segmentation, when the input mesh is strongly nondevelopable. φ: Continuous and bijective function (homeomorphism) φ : M → R 2 that maps M to a planar region in R 2 . In this manuscript, φ is a nearly-isometric map (i.e., highly preserves distances). g: Geodesic distance function g : M × M → R + defined on M. g ij = g(x i , x j ). δ x i Dirac delta (temperature) distribution δ x i : M → [0, ∞] associated to the source point x i , such that the temperature at x i is infinite and 0 everywhere else. u i : Heat kernel function u i : M × (0, T] → R associated to vertex x i ∈ M. u i (x, t) is the temperature at the point x, t, due to an initial infinite-heat point-source δ x i .

χ:
Continuous function χ : R 3 → that indicates if a point p ∈ R 3 is "inside" (χ(p) = 1), "outside" (χ(p) = 0), or "in-between" (χ(p) = 1 n ∑ x∈M χ(x)) a solid defined by M * . U (t) : n × n matrix with the discrete heat kernels maps associated to each vertex in M, i.e., U (t) ij = u i (x j , t). K (t) : n × n matrix of the divergence fields of all heat kernels U (t) . K (t) ij = ∇ · H i (x j , t). G, G 2 : n × n matrices of geodesic and squared geodesic distances, respectively, defined on M. G ij = g(x i , x j ), G 2 ij = g(x i , x j ) 2 . C: n × n symmetric matrix whose entries contain the mean centered squared geodesic distances in M, i.e., C ij = − 1 2 [g 2 ij − 1 n (∑ kl g 2 kl )]. C is semidefinite positive. I n , J n : n × n identity I n and all-ones J n matrices. λ 1 , V 1 : Largest (positive) eigenvalue λ 1 of the semidefinite positive matrix C and its corresponding n × 1 eigenvector V 1 . √ λ 1 V 1 are the discrete u-coordinates of the parameterization Φ ⊂ R 2 . λ 2 , V 2 : Second largest (positive) eigenvalue λ 2 of the semidefinite positive matrix C and its corresponding n × 1 eigenvector V 2 . √ λ 2 V 2 are the discrete v-coordinates of the parameterization Φ ⊂ R 2 .