A Novel Algorithm for Volume-Preserving Parameterizations of 3-Manifolds

(1)

Mei-Heng Yueh^†, Tiexiang Li^‡, Wen-Wei Lin^§, and Shing-Tung Yau^¶

Abstract. Manifold parameterizations have been applied to various fields of commercial industries. Several efficient algorithms for the computation of triangular surface mesh parameterizations have been proposed in recent years. However, the computation of tetrahedral volumetric mesh parameterizations is more challenging due to the fact that the number of mesh points would become enormously large when the higher resolution mesh is considered and the bijectivity of parameterizations is more difficult to guarantee. In this paper, we develop a novel volumetric stretch energy minimization algorithm for volume-preserving parameterizations of simply connected 3-manifolds with a single boundary under the restriction that the boundary is a spherical area-preserving mapping. In addition, our algorithm can also be applied to compute spherical angle- and area-preserving parameterizations of genus-zero closed surfaces, respectively. Several numerical experiments indicate that the developed algorithms are more efficient and reliable compared to other existing algorithms. Numerical results on applications of the manifold partition and the mesh processing for three-dimensional printing are demonstrated thereafter to show the robustness of the proposed algorithm.

Key words. volumetric stretch energy minimization, volume-preserving parameterization, simply connected 3-manifold, manifold parameterization, genus-zero closed surface

AMS subject classifications. 15B48, 52C26, 65F05, 68U05, 65D18

1. Introduction. A manifold parameterization refers to a bijective mapping between the manifold and the domain of a simple canonical shape. The mapping induces a canonical coordinate system on the manifold, which can be used to simplify the issues arising from geometry processing and computer graphics. In particular, surface (2-manifold) parameterizations have been widely studied and applied in various tasks of computer vision, such as surface registration, remeshing, morphing, and texture mapping. Several numerical algorithms for the computation of surface parameterizations have been developed by different research groups in recent years. Most of the parameterization algorithms for surfaces are based on minimizing the distortions of angle or area, or balancing between them. An angle-preserving parameterization is called a conformal parameterization, which has been applied in the texture mapping of surfaces [40, 50] and the image analysis of proximal femur surfaces [55]. On the other hand, an area-preserving parameterization is called an authalic or equiareal parameterization, which has been applied in the surface registration and the remeshing [18, 76]. More details

∗Author’s preprint; last modified August 16, 2022.

Funding: The work of the second author was partially supported by National Natural Science Foundation of China grant 11471074. The work of the authors was partially supported by the Ministry of Science and Technology (MOST), the National Center for Theoretical Sciences (NCTS), the Taida Institute for Mathematical Sciences, the ST Yau Center in Taiwan, the Shing-Tung Yau Center at Southeast University, and the Center of Mathematical Sciences and Applications at Harvard University.

†Department of Mathematics, National Taiwan Normal University, Taipei, Taiwan (yue@ntnu.edu.tw).

‡School of Mathematics, Southeast University, Nanjing, China (txli@seu.edu.cn).

§Department of Applied Mathematics, National Yang Ming Chiao Tung University, Hsinchu, Taiwan (wwlin@math.nctu.edu.tw).

¶Department of Mathematics, Harvard University, Cambridge, MA 02138 (yau@math.harvard.edu).

1

(2)

for methods and applications of surface parameterizations can be found in the survey papers [31,64,42,11,33,39].

In the past, most of the related works merely consider the surface data in three-dimensional (3D) space, which is mathematically a 2D object that can be bijectively mapped onto a planar domain. With the development of 3D imaging technology, such as computed tomography scans and magnetic resonance imaging, 3D images in the real world can be obtained easily. Such a 3D image is mathematically equivalent to a simply connected 3-manifold in R³ that can be bijectively mapped into a unit solid ball while preserving the local volume. To realize the volume-preserving parameterization, Su et al. [67] proposed a useful algorithm based on the discrete optimal mass transportation (OMT) [34]. However, it is usually very time-consuming to compute a required parameterization for 3-manifolds, especially when high-resolution mesh data are considered, e.g., it would cost more than 15 hours¹ using the OMT algorithm to compute a spherical volume-preserving mapping of a volumetric mesh of 290K tetrahedrons.

In addition, the bijectivity of a volumetric mapping is more difficult to guarantee [30,32].

1.1. Contributions. In this paper, we develop a novel volumetric stretch energy minimization (VSEM) for the computation of a volume-preserving parameterization that maps a simply connected volumetric mesh with a single boundary into a unit solid ball. First, the boundary of the 3-manifold is mapped conformally into a unit sphere by minimizing the Dirichlet energy. Then, the boundary mapping is deformed to an equiareal mapping by minimizing the stretch energy. Finally, the volume-preserving parameterization is computed by minimizing the volumetric stretch energy. The contribution of this paper can be separated into three areas.

• Universality: The proposed energy minimization algorithm can be used to compute (i) angle-preserving (conformal) and area-preserving (equiareal) parameterizations, respectively, of genus-zero closed surfaces, and (ii) volume-preserving parameterizations of simply connected 3-manifolds with a single boundary.

• Improved effectiveness and accuracy: The effectiveness and the accuracy of parameterizations computed by the proposed volume-preserving parameterization algorithm are significantly improved compared to the OMT algorithm [67].

• Applications: Applications on the manifold partition for 3-manifolds and the mesh processing for 3D printing can be performed robustly by the proposed parameterization algorithm.

1.2. Notation and overview. The following notation is used in this paper. Other notation will be clearly defined whenever they appear.

• Bold letters, e.g., u, v, w, denote (complex) vectors.

• Capital letters, e.g., A, B, C, denote matrices.

• Typewriter letters, e.g., I, J, K, denote ordered sets of indices.

• v_i denotes the ith entry of the vector v.

• v_I denotes the subvector of v composed of vi for i ∈ I.

• |v| denotes the vector with the ith entry being |v_i|.

• diag(v) denotes the diagonal matrix with the (i, i)th entry being v_i.

1See the computational cost of the OMT algorithm for the Bimba mesh model inTable 6.3.

(3)

• A_i,j denotes the (i, j)th entry of the matrix A.

• A_I,J denotes the submatrix of A composed of A_i,j for i ∈ I and j ∈ J.

• Sⁿ:= {x ∈ Rⁿ⁺¹| ∥x∥ = 1} denotes the n-sphere in Rⁿ⁺¹.

• Bⁿ:= {x ∈ Rⁿ| ∥x∥ ≤ 1} denotes the solid n-ball in Rⁿ.

• [v₀, . . . , v_m] denotes the m-simplex determined by the points v₀, . . . , v_m.

• |[v₀, . . . , vm]| denotes the volume of the m-simplex [v0, . . . , vm].

• i denotes the imaginary unit√

−1.

• I_n denotes the identity matrix of size n × n.

• 0 denotes the zero vectors and matrices of appropriate sizes.

This paper is organized as follows. First, we introduce the related previous works and background in sections 2 and 3, respectively. Then, we propose a novel VSEM algorithm in section 4. The bijectivity of volume-preserving parameterizations is discussed in section 5.

Numerical comparisons between our algorithm and the OMT algorithm are presented in section 6. Applications on the manifold partition and the mesh processing for 3D printing are demonstrated in section 7. Concluding remarks are given in section 8.

2. Previous works. In this section, we briefly review the related previous works on computational algorithms for surface and volumetric parameterizations, respectively.

2.1. Surface parameterizations. The major categories of surface parameterizations con- tain conformal parameterizations, equiareal parameterizations, and balancing parameterizations between angle and area distortions.

A conformal parameterization aims to minimize the angular distortion. Numerical methods for surface conformal parameterizations have been widely developed because the conformal parameterization preserves the local shape well. The existence of conformal mappings between surfaces is guaranteed by the Poincar´e–Klein–Koebe uniformization theorem of 100 years ago. In particular, for a genus-zero closed surface, some feasible numerical methods including the linear Laplace–Beltrami equation [12], the nonlinear heat diffusion [35,43], and the fast landmark aligned spherical harmonic (FLASH) algorithm [19] have been proposed to map it conformally to a unit sphere. In addition, for a simply connected open surface, some efficient numerical methods including the fast disk mapping [20], the linear disk mapping [17], and the conformal energy minimization (CEM) [75] have been recently developed to map it conformally to a unit disk. Furthermore, classical methods including discrete conformal parameterization [22], least-squares conformal mapping [50], angle-based flattening [62, 63,77], discrete conformal equivalence [66], spectral conformal parameterization [53,44], and orbifold Tutte embedding [9,10,8] have been proposed for the conformal parameterization of topological disks with free or mild shape constraints of boundaries. For conformal parameterizations of higher genus surfaces, some well-known methods have been developed, including the holo- morphic one-form [36, 37, 46], the discrete Ricci flow [45, 78], and the discrete Calabi flow [79].

On the other hand, an equiareal parameterization aims to minimize the area distortion.

Numerical methods for computing surface equiareal parameterizations have been gradually getting much more attention in recent years because equiareal mappings preserve the density of the vertices well. Some efficient numerical methods, including the stretch-minimizing method [60,74], the Lie advection method [81], discrete OMT [80,68], the density-equalizing mapping

(4)

[18], and the stretch energy minimization (SEM) [76] have been developed to achieve this purpose.

Additionally, for parameterizations of balancing distortions, some well-developed numerical methods, including the as-rigid-as-possible surface parameterization [51, 72], the most isometric parametrization [41, 21], the isometric distortion minimization [59], and boundary first flattening [61], have been proposed to reach a trade-off between minimizing the angle and area distortions.

2.2. Volumetric parameterizations. The computation of 3-manifold parameterizations is more challenging due to the fact that the number of vertices would become enormously large when the high-resolution volumetric mesh is considered. In addition, the bijectivity of the volumetric mapping is more difficult to guarantee because convex combination mappings in 3D space do not need to be bijective [32]. Note that conformal mappings between 3-manifolds, in general, do not exist. Still, volumetric mappings with small angular distortion are frequently adopted. Paillé and Poulin [56] and Chern, Pinkall, and Schröder [16] proposed conformal-based volumetric mapping methods by applying the Cauchy–Riemann equation to each canonical orthogonal plane in R³ and a low shear distortion on the decoupling of scaling and rotation, respectively. Kovalsky et al. [48, 49] developed methods to deform a given volumetric mapping into a bijective one with bounded distortion. Paillé et al. [57] proposed a spectral method based on the dihedral angle representation for computing the locally injective mapping of tetrahedral meshes. Jin et al. [47] proposed a method for computing the volumetric parameterization balancing between angle and volume distortions by minimizing the stretch-distortion energy. Naitsat, Saucan, and Zeevi [54] and Rabinovic et al. [59] proposed deformation methods for volumetric meshes based on the quasi-conformal homeomorphism and minimizing a linear combination of local isometric distortion measures [65], respectively.

The above methods have mainly been developed for the computation of bijective volumetric mappings of 3-manifolds by minimizing angular or isometric distortions so that the qualities on angles and local shapes of tetrahedral meshes are well-preserved, but in general they are not volume-preserving.

On the other hand, it is worth noting that only a few of the existing algorithms consider the volumetric parameterizations for prescribed shapes that minimize the volume distortion.

On the basis of OMT, Su et al. [67] developed a pioneering algorithm of volume-preserving parameterizations for tetrahedral meshes that maps a simply connected 3-manifold with a single boundary into a unit solid ball, which has advantages on manifold resampling and registration.

3. Discrete manifolds and parameterizations. A 3-manifold refers to a 3D topological space in which each point of the manifold has a neighborhood homeomorphic to a subset of R³. In this paper, we consider simply connected discrete 3-manifolds with a single boundary that are embedded in R³. In practice, a discrete 3-manifold is represented as a tetrahedral mesh M composed of n vertices with coordinates in R³

V(M) =n

v_s≡ v_s¹, v_s², v_s³⊤

∈ R³on s=1

(5)

and tetrahedrons

T (M) =[v_i, v_j, v_k, v_ℓ] ⊂ R³ for some vertices {v_i, v_j, v_k, v_ℓ} ⊂ V(M) .

Here the bracket [v_i, v_j, v_k, v_ℓ] denotes the convex hull (3-simplex) of the affinely independent points {vi, vj, vk, vℓ}. Furthermore, we denote the set of triangular faces and edges of the mesh M by

F (M) = {[v_i, v_j, v_k] | [v_i, v_j, v_k, v_ℓ] ∈ T (M) for some v_ℓ∈ V(M)}

and

E(M) = {[v_i, v_j] | [v_i, v_j, v_k] ∈ F (M) for some v_k∈ V(M)} ,

respectively. The union of T (M), F (M), E (M), and V(M) forms a homogeneous simplicial 3-complex. Similarly, a 2-manifold is called a surface, which is represented as a triangular mesh composed of vertices and triangular faces.

A piecewise affine mapping f : M → R³ can be represented as a matrix (3.1) f = f(v1) · · · f (vn)⊤

≡

f1 · · · fn

⊤

∈ R^n×3.

Note that for the point v ∈ M which is not a vertex, v must belong to a tetrahedron [v_i, v_j, v_k, v_ℓ] ∈ T (M). Then, the value f (v) can be represented as a linear combination of fi, fj, fk, and fℓ with coefficients being the barycentric coordinates of v in [vi, vj, vk, vℓ], that is,

f |_[v_i_,v_j_,v_k_,v_ℓ_](v) = λi(v)fi+ λj(v)fj+ λk(v)fk+ λℓ(v)fℓ, where λ_i(v) = _|[v^|[v,v^j^,v^k^,v^ℓ^]|

i,vj,vk,vℓ]|, λ_j(v) = _|[v^|[vⁱ^,v,v^k^,v^ℓ^]|

i,vj,vk,vℓ]|, λ_k(v) = _|[v^|[vⁱ^,v^j^,v,v^ℓ^]|

i,vj,vk,vℓ]|, and λ_ℓ(v) = _|[v^|[vⁱ^,v^j^,v^k^,v]|

Remark. The matrixf represents the unique piecewise affine mapping satisfying f(vs) =fs

for every v_s∈ V(M).

In this paper, the considered parameterization of a 3-manifold M ⊂ R³ is a bijective piecewise affine mapping f that maps M to the unit solid ball B³. A parameterization f : M → B³is said to be volume-preserving if the Jacobian matrix J_f⁻¹ =h

∂f⁻¹

∂u¹

∂f⁻¹

∂u²

∂f⁻¹

∂u³

i satisfies

det J_f⁻¹(u¹, u², u³) = 1.

In other words, f preserves the local volume.

In addition, for a genus-zero closed triangular mesh S, a parameterization f : S → S² is said to be conformal or angle-preserving if the first fundamental form I_f⁻¹ satisfies

I_f⁻¹(u¹, u²) = λ(u¹, u²)I₂

for some positive scaling function λ, i.e., f preserves local angles. A parameterization f : S → S² is said to be equiareal or area-preserving if the first fundamental form I_f⁻¹ satisfies det I_f⁻¹(u¹, u²) = 1, i.e., f preserves the local area.

(6)

4. The novel algorithm. In this section, we develop a novel energy minimization algorithm for the computation of a volume-preserving parameterization that maps a simply connected 3-manifold with a single boundary to a unit solid ball. Given a simply connected tetrahedral mesh M with the boundary ∂M = S being a genus-zero closed surface, our algorithm is composed of three stages as follows. First, the initial boundary mapping is chosen to be a spherical conformal parameterization f : S → S², which is computed by minimizing the discrete Dirichlet energy [23]. Then, the boundary mapping is deformed into an equiareal parameterization, which is iteratively computed by minimizing the stretch energy [76]. Finally, the volume-preserving parameterization with the prescribed boundary mapping is computed by minimizing the volumetric stretch energy.

4.1. Initial conformal boundary parameterizations. The initial spherical conformal parameterization is carried out by applying the CEM algorithm [75]. The original CEM algorithm aims to compute disk-shaped conformal parameterizations of simply connected open surfaces by minimizing the discrete Dirichlet energy functional [23,38]

(4.1) E_D(f ) = 1

2tracef^⊤LDf , where LD is the Laplacian matrix with

(4.2) [LD]i,j =







−w_i,j ≡ −¹₂(cot α_i,j+ cot α_j,i) if [v_i, v_j] ∈ E (S), P

k̸=iwi,k if j = i,

0 otherwise,

in which α_i,jand α_j,iare two angles opposite to the edge [v_i, v_j]. The modified CEM algorithm for genus-zero closed surfaces is stated as follows. First, an initial mapping h⁽⁰⁾ : S → C is obtained by solving the Laplace–Beltrami equation

(4.3) △_Sh⁽⁰⁾ = ∂

∂x− i ∂

∂y

δp,

where δ_p is the Dirac delta function at a selected point p on S, and (x, y) are the local coordinates defined on a neighborhood of p. This method was originally proposed by Angenent et al. [12] and has a benefit that the resulting spherical parameterization is bijective. A straightforward proof for the bijectivity of the method [12] is given in Theorem 5.3. The algorithm for solving (4.3)is summarized in Algorithm 4.1.

Unsatisfactorily, the angular distortion of the mapping computed byAlgorithm 4.1would be relatively large at the neighborhood of p. To remedy this drawback, we improve the conformality of the mapping by iteratively minimizing the Dirichlet energy

(4.6)

(△_Sh^(k)(v) = 0 if

Inv ◦ h^(k−1)(v) < r, h^(k)(v) = Inv ◦ h^(k−1)(v) otherwise,

where Inv denotes the inversion Inv(z) = ¹_z. Once the iteration (4.6) converges, the desired spherical conformal parameterization is obtained by the inverse stereographic projection Π⁻¹

S² :

(7)

Algorithm 4.1 Spherical conformal parameterizations [12]

Input: A genus-zero closed mesh S.

Output: A spherical conformal parameterizationf.

1: Let n be the number of vertices of S.

2: Find the most regular triangular face by

[v_a, v_b, v_c] = argmin

[vi,vj,vk]∈F (S)





|[v_i, v_j]| − ¹₃(|[v_i, v_j]| + |[v_j, v_k]| + |[v_k, v_i]|)

|[v_j, v_k]| − ¹₃(|[v_i, v_j]| + |[v_j, v_k]| + |[v_k, v_i]|)

|[v_k, vi]| − ¹₃(|[vi, vj]| + |[vj, vk]| + |[vk, vi]|)



 .

3: Set B = {a, b, c} and I = {1, . . . , n}\B.

4: Set α = ^(v^c^−v_∥v^a⁾^⊤^(v^b^−v^a⁾

b−va∥²₂ and

(4.4) h_B =







−1

∥v_b−v_a∥₂ 1

∥vb−va∥2

0





+ i







1−α

∥vc−(va+α(vb−va))∥2

α

−1





.

5: Compute h by solving the linear system

(4.5) [LD]I,IhI = −[LD]I,BhB, where LD is defined as in(4.2).

6: The spherical parameterizationf is obtained by Π⁻¹_S2(h_ℓ) in (4.7), ℓ = 1, . . . , n.

C → S² defined by

(4.7) Π⁻¹

S²(z) =

2 Re z

|z|²+ 1, 2 Im z

|z|²+ 1,|z|²− 1

|z|²+ 1

.

The improved algorithm derived by (4.6)for the computation of spherical conformal parameterizations is summarized inAlgorithm 4.2.

Remark. In fact, the FLASH algorithm [19] for spherical conformal parameterization is satisfactory with high accuracy and effectiveness. However, the accuracy and effectiveness of Algorithm 4.2are slightly better than those of the FLASH algorithm [19]. Numerical comparisons between the FLASH algorithm [19] andAlgorithm 4.2are demonstrated inAppendix C.

4.2. Equiareal boundary parameterizations. The boundary spherical equiareal parameterization is carried out by applying the SEM algorithm [76]. The original SEM algorithm [76]

aims to compute disk-shaped equiareal parameterizations of simply connected open surfaces by minimizing the stretch energy functional

(4.8) E_S(f ) = 1

2tracef^⊤L_S(f )f ,

(8)

Algorithm 4.2 CEM for spherical conformal parameterizations

Input: A genus-zero closed mesh S, a tolerance ε, a radius r (e.g., ε = 10⁻⁶, r = 1.2).

Output: A spherical conformal parameterizationf.

2: Compute a spherical conformal parameterizationg usingAlgorithm 4.1.

3: Perform the stereographic projection h_ℓ ← ₁₋^g^ℓ,1_g

ℓ,3 + i₁₋^g^ℓ,2_g

ℓ,3, ℓ = 1, . . . , n.

4: Let δ ← ∞.

5: while δ > ε do

6: Perform the inversion h ← diag(|h|)⁻²h.

7: Update the index sets I = {i | |hi| < r } and B = {1, . . . , n}\I.

8: Update h by solving the linear system

[LD]I,IhI= −[LD]I,BhB, where LD is defined as in(4.2).

9: Setfℓ ← Π⁻¹

S² (h_ℓ) as in(4.7), ℓ = 1, . . . , n.

10: Update δ ← E_D(g) − ED(f).

11: Updateg ← f.

12: end while

13: return The spherical parameterizationf.

where LS(f ) is the stretch Laplacian matrix with

(4.9) [LS(f )]_i,j =











−w_i,j(f ) ≡ −¹₂

_cot(α

i,j(f ))

σ_{f −1}([vi,vj,vk]) +_σ^cot(α^j,i^{(f ))}

f −1([vj,vi,vℓ])

if [vi, vj] ∈ E (S), P

ℓ̸=iwi,ℓ(f ) if j = i,

0 otherwise

in which σ_f⁻¹([v_i, v_j, v_k]) = _{|f ([v}^|[vⁱ^,v^j^,v^k^]|

i,vj,vk])| is the stretch factor of f on the triangular face [vi, vj, vk]. The modified SEM algorithm for genus-zero closed surfaces is stated as follows.

First, the initial mapping f⁽⁰⁾ is a spherical conformal parameterization computed by Al- gorithm 4.2 and let h⁽⁰⁾_ℓ = ^f

(0) ℓ,1

1−fℓ,3⁽⁰⁾

+ i ^f

(0) ℓ,2

1−fℓ,3⁽⁰⁾

, ℓ = 1, . . . , n. Then, the stretch energy (4.8) is iteratively minimized by

(4.10) [LS(f^(k))]_I(k),I^(k)h^(k+1)

I^(k) = −[LS(f^(k))]_I(k),B^(k)h^(k)

B^(k), where LS is defined in (4.9),

I^(k)= n

ℓ

|(diag(|h^(k)|)⁻²h^(k))_ℓ| < ro

and B^(k)= {1, . . . , n}\I^(k).

In practice, the radius r is chosen to be 1.2. Once h^(k) in the iteration (4.10) converges to h^(∗), the desired spherical equiareal parameterizationf is obtained by the inverse stereographic projection operates on h^(∗).

The algorithm for the computation of spherical equiareal parameterizations by (4.10) is summarized inAlgorithm 4.3.

(9)

Algorithm 4.3 SEM for spherical equiareal parameterizations

Input: A genus-zero closed mesh S, a tolerance ε, a radius r (e.g., ε = 10⁻⁶, r = 1.2).

Output: A spherical equiareal parameterization f.

2: Compute a spherical conformal parameterizationg usingAlgorithm 4.2.

3: Perform the stereographic projection h_ℓ = ₁₋^g^ℓ,1_g

ℓ,3+ i₁₋^g^ℓ,2_g

ℓ,3, ℓ = 1, . . . , n.

4: Let δ ← ∞.

5: while δ > ε do

6: Update the matrix A ← L_S(f ), where L_S(f ) is defined as in (4.9).

7: Perform the inversion h ← diag(|h|)⁻²h.

8: Update the index sets I = {i | |h_i| < r } and B = {1, . . . , n}\I.

9: Update h by solving the linear system A_I,Ih_I= −A_I,Bh_B.

10: Setfℓ ← Π⁻¹

S² (hℓ) as in(4.7), ℓ = 1, . . . , n.

11: Update δ ← E_S(g) − ES(f).

12: Updateg ← f.

13: end while

14: return The spherical parameterizationf.

4.3. Volume-preserving parameterizations. Let M be a tetrahedral mesh of n vertices in R³ and let f : M → R³ be a piecewise affine mapping on M. The volumetric Dirichlet energy functional [69,70,71] is defined by

(4.11) E_D(f ) = 1

2tracef^⊤L_Df , where L_D∈ R^n×n is the volumetric Laplacian matrix with

(4.12) [L_D]i,j =







−w_i,j if [v_i, v_j] ∈ E (M), P

k̸=iwi,k if j = i,

0 otherwise,

in which wi,j is the cotangent weight on the edge [vi, vj] given by

(4.13) w_i,j = 1

6

X

τ ∈T (M) [vi,vj]∪[v_k,v_ℓ]⊂τ [vi,vj]∩[vk,vℓ]=∅

|[v_k, v_ℓ]| cot θ_i,j^k,ℓ,

where θ^k,ℓ_i,j is the dihedral angle between [v_i, v_k, v_ℓ] and [v_j, v_ℓ, v_k] in the tetrahedron τ on the edge [v_k, v_ℓ], as illustrated inFigure 4.1.

Similarly as [58,73, 76], the cotangent formula (4.13)can be fairly changed according to the image of a mapping f by

(4.14) wi,j(f ) = 1 6

X

τ ∈T (M) [vi,vj]∪[vk,vℓ]⊂τ [vi,vj]∩[v_k,v_ℓ]=∅

|f ([v_k, v_ℓ])| cot θ_i,j^k,ℓ(f ),

(10)

θ^k,ℓ_i,j p v_i

vℓ

v_k vj

Figure 4.1. An illustration for the dihedral angle between [vi, vk, vℓ] and [vj, vℓ, vk] in the tetrahedron [vi, vj, vk, vℓ].

where θ_i,j^k,ℓ(f ) is the dihedral angle between the triangular faces f ([v_i, v_k, v_ℓ]) and f ([v_j, v_ℓ, v_k]) in the tetrahedron f (τ ). We now define the stretch factor of f on a tetrahedron τ ∈ T (M) as

(4.15) σ_f⁻¹(τ ) = |τ |

|f (τ )|.

In addition, for the purpose of the volume-preserving parameterization, the condition |f (τ )| =

|τ |, for each τ ∈ T (M), is required to be satisfied. We now impose the stretch factor σ_f⁻¹(τ ) into wi,j(f ) and modify the formula as

(4.16) wi,j(f ) = 1 6

X

τ ∈T (M) [vi,vj]∪[vk,vℓ]⊂τ [vi,vj]∩[vk,vℓ]=∅

|f ([v_k, v_ℓ])| cot θ_i,j^k,ℓ(f ) σ_f⁻¹(τ ) ,

where σ_f⁻¹ is the stretch factor as in(4.15), which is equivalent to the Jacobian determinant (see Appendix A).

The volumetric stretch energy functional on M is defined as

(4.17) E_S(f ) = 1

2tracef^⊤L_S(f )f , where L_S(f ) is the stretch volumetric Laplacian matrix with

(4.18) [L_S(f )]_i,j =







−w_i,j(f ) if [v_i, v_j] ∈ E (M), P

ℓ̸=iwi,ℓ(f ) if j = i,

0 otherwise,

in which wi,j(f ) is the modified weight given in(4.16).

Suppose an equiareal spherical boundary mapping fB⁽⁰⁾ is computed by Algorithm 4.3.

Then, the volume-preserving parameterization is computed by minimizing the volumetric

(11)

stretch energy (4.17)via the iteration

(4.19) [L_S(f^(k))]I,IfI^(k+1) = −[L_S(f^(k))]I,BfB⁽⁰⁾

for solving the sequential quadratic programmings, where the matrix L_S is defined in (4.18), B = {s | v_s∈ ∂M} and I = {1, . . . , n}\B.

The VSEM algorithm for the computation of volume-preserving parameterizations by (4.19)is summarized in Algorithm 4.4.

Remark. A mapping f is volume-preserving if and only if the stretch factor σ_f⁻¹(τ ) = 1 for every τ ∈ T (M). As a result, f is a fix point of the iteration(4.19), i.e., a critical point of the energy functional E_Sin(4.17). However, it is not obvious the other way around since the energy functional(4.17)is highly nonlinear. So, it is reasonable to numerically check whether the resulting stretch factor σ_f⁻¹(τ ) is very close to 1 for every τ ∈ T (M). The distributions of σ_f⁻¹ for various benchmark mesh models are demonstrated later in Figure 6.4.

Algorithm 4.4 VSEM for volume-preserving parameterizations

Input: A simply connected tetrahedral mesh M, a tolerance ε (e.g. ε = 10⁻⁶).

Output: A volume-preserving parameterization f.

1: Let n be the number of vertices of M.

2: Let B = {s | v_s∈ ∂M} and I = {1, . . . , n}\B.

3: Compute a spherical equiareal parameterization gB by Algorithm 4.3.

4: Compute g by solving the linear system

[L_D]_I,IgI = −[L_D]_I,BgB, where L_D is defined as in(4.12).

5: Let δ ← ∞.

6: LetfB←gB.

7: while δ > ε do

8: Update A ← L_S(g), where L_S(g) is defined as in (4.18).

9: Updatef by solving the linear system AI,IfI = −AI,BfB.

10: Update δ ← E_S(g) − E_S(f).

11: Updateg ← f.

12: end while

13: return The volumetric mappingf.

5. Bijectivity of the parameterizations. The bijectivity is one of the most important in the parameterization process. In the following, we first show that the spherical conformal parameterization of a genus-zero closed mesh produced by Algorithm 4.1 is bijective under the Delaunay assumption.

For convenience, we give the definition of an M-matrix [13] and a related lemma.

Definition 5.1. A matrix A ∈ R^m×n is nonnegative (positive) if all entries of A are nonnegative (positive). A squared matrix A ∈ R^n×n is an M-matrix if A = sI − B with B being nonnegative and s ≥ ρ(B), where ρ(B) is the spectral radius of B.

(12)

Lemma 5.2 (Theorems 1.4.7 and 1.4.10 in [52]). Suppose A ∈ R^n×n is a singular, irreducible M-matrix. Then each principal submatrix ˆA of A other than A itself is a nonsingular M-matrix and ˆA⁻¹ is nonnegative.

Theorem 5.3. Given a genus-zero closed Delaunay triangular mesh M of n vertices, the spherical conformal parameterization of M produced by Algorithm 4.1 is bijective.

Proof. The Delaunay property guarantees that the weights wi,j of L_D in(4.2)are positive.

Let D = diag (P_n

j=1w_1,j, . . . ,P_n

j=1w_n,j)^⊤ . Then

D⁻¹L_Dh

i= X

vj∈N (v_i)

w_i,j Pn

j=1wi,j

(h_j− h_i) , i = 1, . . . , n,

where N (v_i) is the set of neighboring vertices of v_i. SinceP

vj∈N (vi) wi,j

Pn

j=1wi,j = 1, D⁻¹L_Dh = 0 is equivalent to

(5.1) hi = X

vj∈N (vi)

w_i,j Pn

j=1wi,j

hj, i = 1, . . . , n.

It is known that L_D is a singular, irreducible M-matrix [75]. From Lemma 5.2, the matrix [LD]I,I in (4.5) is invertible and −[LD]⁻¹_I,I[LD]I,B is nonnegative. Therefore, from (5.1), the unique mapping h_I= [L_D]⁻¹_I,I[L_D]_I,Bh_B obtained by solving the linear system(4.5)is a convex combination mapping with the boundary hB in (4.4) being a triangle. From the result in [30], the interior mapping hI in (4.5) is bijective. Hence, the bijectivity of the spherical conformal parameterization f in step 6 of Algorithm 4.1 follows from the bijectivity of the inverse stereographic projection.

On the other hand, the bijectivity of the spherical conformal parameterization computed by the CEM algorithm, Algorithm 4.2, can be theoretically guaranteed if the boundary mapping h_B in step 8 forms a convex polygon. Unfortunately, the polygon formed by h_B is not necessarily convex, so the result in [30] cannot be applied. However, from the construction of the index set B in step 7 of Algorithm 4.2, the polygon formed by hB is nearly convex, i.e., there exists a circle γ : S¹→ C with radius r satisfying

(5.2) sup

x∈S¹,j∈B

|γ(x) − h_j| < sup

[vi,vj]∈E(S) i /∈B,j∈B

|h_i− h_j|.

It is reasonable to demonstrate the bijectivity of the spherical conformal parameterization via numerical checking of the resulting mapping computed by Algorithm 4.2. The method for checking the bijectivity of a discrete spherical mapping f : S → S² can be found in Appendix B.1.

As we will see in Table C.2, the spherical conformal parameterization computed by the CEM algorithm, Algorithm 4.2, is bijective, for each tested benchmark mesh model.

In addition, for the SEM algorithm,Algorithm 4.3, the resulting mapping seldom fails to be bijective, e.g., for the boundary mapping of the Bunny model,² the image of the resulting

2Among all 17 tested surface mesh models, the Bunny model is the only one for which the resulting spherical mapping by the SEM algorithm,Algorithm 4.3, is not bijective.

(13)

spherical mapping f contains three overlapped triangular faces. It is not surprising since the stretch cotangent weights in (4.9)are not necessarily positive so that the resulting mapping is not inevitably a convex combination mapping. To remedy this drawback, we unfold the overlapped triangular faces on the image of f by the procedures inAppendix B.2.

Furthermore, for the simply connected tetrahedral mesh with a single boundary, the convex combination property of the mapping follows similarly to Theorem 5.3.

Corollary 5.4. Given a simply connected tetrahedral mesh M of n vertices with the boundary and interior indices being B and I, respectively, suppose L is a Laplacian matrix as in (4.18) with positive weights and f : M → R³ is a piecewise affine mapping satisfying that the boundary mapping {fs:= f (vs) | s ∈ B } forms a convex polyhedron. Then LI,I is invertible and the mapping

(5.3) fI= −L⁻¹_I,IL_I,BfB

is a convex combination mapping.

Remark. The bijectivity of the volume-preserving parameterization of tetrahedral meshes, in general, is not guaranteed even if the mapping is a convex combination mapping with the boundary mapping being convex. An elegant counterexample has been given in [32].

6. Numerical experiments. In this section, we demonstrate numerical experiments of the VSEM algorithm for volume-preserving parameterizations of simply connected tetrahedral meshes. All linear systems in our algorithms are solved by the backslash operator (\) in MATLAB. Some of surface and tetrahedral mesh models are obtained from TurboSquid [7], the AIM@SHAPE shape repository [3], the Stanford 3D scanning repository [6], a project page of ALICE [1], and Gu’s website [2]. Some of tetrahedral mesh models are generated using JIGSAW mesh generators [24,27,25,28,26].

In order to better understand the shape of the mesh models and the image of the volume- preserving parameterizations computed by the VSEM algorithm,Algorithm 4.4, we first show the tetrahedral meshes of David Head and Human Brain and their volume-preserving parameterizations inFigures 6.1 and 6.2, respectively.

Then, we compare the effectiveness and accuracy of the VSEM algorithm,Algorithm 4.4, to the state-of-the-art OMT algorithm [67]. The executable program files of the OMT algorithm are obtained from Gu’s website [2].

We now introduce the volume distortion to measure the accuracy of a volume-preserving parameterization by the total volume distortion as well as the mean and standard deviation (SD) of local volume ratios. The total volume distortion of a mapping f on a mesh M is defined as

(6.1) D_M(f ) = 1

4 X

v∈V(M)

P

τ ∈N (v)|τ |

|M| −

P

τ ∈N (v)|f (τ )|

|f (M)|

,

where N (v) = {τ ∈ T (M) | v ⊂ τ } is the set of neighboring tetrahedrons of the vertex v, and |M| and |f (M)| denote volumes of M and its image, respectively. A mapping f is

(14)

volume-preserving if DM(f ) = 0. The local volume ratio R_f on a vertex v is defined as

(6.2) R_f(v) =

P

τ ∈N (v)|τ |/|M|

P

τ ∈N (v)|f (τ )|/^{|f (M)|}.

The mean and SD of R_f(v), for all v ∈ V(M), are used to measure the local volume distortion.

A mapping f is volume-preserving if the mean is 1 and the SD is 0.

In Tables 6.1 and 6.2 and Figure 6.3, we show the total volume distortion as well as the mean and SD of local volume ratios of volume-preserving parameterizations of various tested benchmark mesh models obtained by OMT and the VSEM algorithm, Algorithm 4.4, respectively. We see that the VSEM has much smaller total volume distortions and local volume ratios than those of OMT. In addition, we show the execution time of OMT and the VSEM inTable 6.3, which indicates that the effectiveness of the VSEM outperforms the OMT by several hundreds of times.

Furthermore, in Figure 6.4, we show histograms of the stretch factors σ_f⁻¹(τ ), for τ ∈ T (M), of volume-preserving parameterizations f computed by the VSEM algorithm, Algo- rithm 4.4. It indicates that the stretch factors of most tetrahedrons are closed to 1, as we desired.

Recall that the bijectivity of convex combination mappings on 3-manifolds is, in general, not guaranteed [32]. It is reasonable to numerically check the bijectivity of volume-preserving parameterizations computed by the VSEM algorithm,Algorithm 4.4. The method for checking the bijectivity of a volumetric mapping can be found in Appendix B.3. InTable 6.4, we show the percentages of bijective tetrahedrons of volume-preserving parameterizations by OMT and VSEM, respectively. We observe that most of tetrahedrons are mapped bijectively by OMT and VSEM, and the bijectivity of the VSEM is much better than that of OMT.

Remark. Thanks to the large-scale bounded distortion mapping [4, 49], for each tested tetrahedral volumetric mesh model listed in Table 6.4, the overlapped tetrahedrons can be unfolded into a 100% bijective volume-preserving parameterization, slightly sacrificing the total volume distortion and the spherical boundary.

In summary, the proposed VSEM algorithm, Algorithm 4.4, has better accuracy and effectiveness than OMT that would cost less than 10 seconds for computing volume-preserving parameterizations of meshes of 100,000 tetrahedrons. It should be quite satisfactory in prac- tical applications.

7. Applications. In this section, we present sample applications of volume-preserving parameterizations for 3-manifolds, namely, the manifold partition and the mesh processing for 3D printing.

7.1. Manifold partitions. The aim of the manifold partition problem is to find an optimal partition of a manifold in terms of minimizing a certain energy functional [15, 14]. In the previous works, the considered manifolds are usually of simple shapes, e.g., a disk, a sphere, a solid ball, and a cube. It is worth noting that the optimal partition usually has the property that each part has a similar volume. With aid of the VSEM algorithm, Algorithm 4.4, for volume-preserving parameterizations, a manifold can be easily partitioned into several simply

(15)

Table 6.1

The total volume distortion (6.1) of volume-preserving parameterizations by OMT and VSEM algorithm, Algorithm 4.4. Here ”—” means the executable program file of OMT does not work.

Model # Tetrahedrons OMT [67] VSEM

Bunny 84,787 0.0627 0.0322

Fandisk 88,374 0.0645 0.0506

David Head 90,575 0.0627 0.0142

Lion 97,131 0.0666 0.0262

Max Planck 115,767 0.0400 0.0086

Venus 130,429 0.0346 0.0052

Apple 169,888 0.0300 0.0030

Bimba 290,353 0.0500 0.0184

Human Brain 852,565 — 0.0790

Table 6.2

The mean and SD of local volume ratios of volume-preserving parameterizations by OMT and VSEM algorithm,Algorithm 4.4. Here ”—” means the executable program file of OMT does not work.

Mean SD Mean SD

Bunny 84,787 1.5015 42.5488 1.1044 1.4497

Fandisk 88,374 1.3520 6.3902 1.2468 2.8231

David Head 90,575 1.0621 0.6465 1.0323 0.6634

Lion 97,131 1.0417 3.1432 1.0246 0.1215

Max Planck 115,767 1.0568 1.6567 1.0235 0.7186

Venus 130,429 1.0183 0.4585 1.0066 0.3662

Apple 169,888 1.0048 0.0564 1.0002 0.0086

Bimba 290,353 1.0342 1.6458 1.0066 0.4764

Human Brain 852,565 — — 1.0266 0.8850

Table 6.3

The computational cost (second) of volume-preserving parameterizations by OMT and VSEM algorithm, Algorithm 4.4. Here ”—” means the executable program file of OMT does not work.

Time #Iter. Time #Iter.

Bunny 84,787 24438.40 231 6.26 10

Fandisk 88,374 8542.10 61 6.35 10

David Head 90,575 6911.32 24 6.13 10

Lion 97,131 13716.90 50 6.95 10

Max Planck 115,767 6679.72 47 9.81 10

Venus 130,429 5558.35 40 13.35 10

Apple 169,888 5599.84 37 20.08 10

Bimba 290,353 54856.90 193 33.09 10

Human Brain 852,565 — — 157.79 10

(16)

Figure 6.1. The tetrahedral mesh of David Head (left) and its volume-preserving parameterization (right) by VSEM algorithm,Algorithm 4.4.

connected submanifolds with a similar volume that could be good initially for solving the partition problem [15,14].

First, a 3-manifold M is mapped into a unit solid ball B³ by f using VSEM. Then, the partition of M can be achieved by slicing B³ with suitable cutting planes. The uniform sampling of B³ can be controlled by using the spherical coordinate system, and the tetrahedral mesh can be constructed by the Delaunay triangulation algorithm (delaunayTriangulation) in MATLAB due to the fact that B³ is convex.

Figure 7.1 shows the equal-volume partition of the David Head model M. The cutting planes are chosen to be x = 0, y = 0 and z = 0 so that the uniform tetrahedral mesh of a ball B³ is partitioned into eight equal-volume partitions. Then, each partition of B³ is mapped to M via f⁻¹ and the corresponding volume is shown inTable 7.1. Note that, here, the volume of M is normalized to be 1. From Table 7.1, we see that each partition of the manifold has a volume around ¹/8with SD less than 2.73 × 10⁻⁴, which is quite satisfactory.

7.2. Mesh processing for 3D printing. Nowadays, 3D printing is a popular topic. Issues arising from the data processing for 3D printing are getting much attention. In order to save printing material, it is common to hollow out the interior of the model when a solid object is

(17)

Figure 6.2.The tetrahedral mesh of Human Brain (left) and its volume-preserving parameterization (right) by VSEM algorithm,Algorithm 4.4.

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07

Total Volume Distortion

Bunny Fandisk David Head

Lion Max Planck

Venus Apple Bimba

OMT VSEM

Figure 6.3. The total volume distortion (6.1) of volume-preserving parameterizations by OMT and VSEM algorithm,Algorithm 4.4.

(18)

0.5 1 1.5 Local Volume Ratios 0

500 1000 1500 2000 2500

0.5 1 1.5

Stretch Factors 0

1000 2000 3000 4000 5000

0.5 1 1.5

Stretch Factors 0

500 1000 1500 2000

0.5 1 1.5

Stretch Factors 0

1000 2000 3000 4000 5000

(a) Bunny (b) Fandisk (c) David Head (d) Lion

0.5 1 1.5

Stretch Factors 0

500 1000 1500 2000

0.5 1 1.5

Stretch Factors 0

5000 10000 15000

0.5 1 1.5

Stretch Factors 0

5000 10000 15000

0 1 2

Stretch Factors 0

5 10 ×10⁴

(e) Max Planck (f) Venus (g) Bimba (h) Human Brain

Figure 6.4. Histograms of stretch factors σ_f−1(τ ) of volume-preserving parameterizations obtained by VSEM algorithm,Algorithm 4.4.

Table 6.4

The percentages of bijective tetrahedrons for volume-preserving parameterizations by OMT and VSEM algorithm,Algorithm 4.4. Here ”—” means the executable program file of OMT does not work.

Bunny 84,787 99.19% 99.90%

Fandisk 88,374 98.79% 99.94%

David Head 90,575 99.90% 99.99%

Lion 97,131 99.58% 99.97%

Max Planck 115,767 99.91% 99.99%

Venus 130,429 99.96% 99.99%

Apple 169,888 99.96% 100.00%

Bimba 290,353 99.29% 99.96%

Human Brain 852,565 — 99.81%

Front View Left Side View Back View Right Side View

Figure 7.1. The partition of the David Head model in four views.

(19)

Table 7.1

The partition of the David Head model and their volumes. The total volume is normalized to be 1.

Part Volume Part Volume Part Volume Part Volume

0.1249 0.1247 0.1257 0.1247

0.1248 0.1250 0.1255 0.1247

SD of volumes: 2.7227 × 10⁻⁴

printed. With aid of volume-preserving parameterizations, the thickness of the model can be easily controlled. First, a 3-manifold M is mapped into the ball B³ with VSEM. Then, the uniform tetrahedral mesh of B³ is constructed by the spherical coordinate system as before.

In Figure 7.2, we show the hollowed David Head and Max Planck models in three views.

The capacity of the interior is roughly 63% of the volume of the original solid manifold. From Figure 7.2, we see that the thickness of each hollowed model is almost equally distributed with a smooth interior surface.

8. Concluding remarks. In this paper, we develop a novel algorithm by minimizing the volumetric stretch energy, which can be used to compute volume-preserving parameterizations for simply connected tetrahedral meshes with a single boundary. In addition, we generalize the CEM [75] and SEM [76] algorithms to compute the spherical conformal and equiareal parameterizations of genus-zero closed triangular meshes, respectively. Numerical experiments indicate that the proposed VSEM algorithm for volume-preserving parameterizations outperforms the state-of-the-art OMT algorithm [67] with better effectiveness and accuracy. Appli- cations of the manifold partition and the mesh processing for 3D printing are demonstrated to show the practicality of our algorithm.

Acknowledgments. The authors want to thank Prof. Xianfeng David Gu, Prof. Lok Ming Ronald Lui and Mr. Gary Pui-Tung Choi for useful discussions and the executable program files of the OMT and the FLASH algorithms, respectively.

REFERENCES

[1] ALICE. http://alice.loria.fr/(2016).

[2] David Xianfeng Gu’s Home Page. http://www3.cs.stonybrook.edu/^∼gu/(2017).

[3] Digital Shape Workbench - Shape Repository. http://visionair.ge.imati.cnr.it/ontologies/shapes/(2016).

[4] Large-scale bounded distortion mappings. https://shaharkov.github.io/LargeScaleBD.html(2018).

[5] Spherical Conformal Map - File Exchange - MATLAB Central. https://www.mathworks.com/

matlabcentral/fileexchange/65551-spherical-conformal-map(2018).

[6] The Stanford 3D Scanning Repository. http://graphics.stanford.edu/data/3Dscanrep/(2016).

[7] TurboSquid. http://www.turbosquid.com/(2016).

[8] N. Aigerman, S. Z. Kovalsky, and Y. Lipman, Spherical orbifold tutte embeddings, ACM Trans.

Graph., 36 (2017), pp. 90:1–90:13,https://doi.org/10.1145/3072959.3073615.

(20)

Front View Side View Back View

Figure 7.2. The hollowed David Head and Max Planck models in three views.

[9] N. Aigerman and Y. Lipman, Orbifold tutte embeddings, ACM Trans. Graph., 34 (2015), pp. 190:1–

190:12,https://doi.org/10.1145/2816795.2818099.

[10] N. Aigerman and Y. Lipman, Hyperbolic orbifold tutte embeddings, ACM Trans. Graph., 35 (2016), pp. 217:1–217:14,https://doi.org/10.1145/2980179.2982412.

[11] P. Alliez, G. Ucelli, C. Gotsman, and M. Attene, Recent Advances in Remeshing of Sur- faces, Springer Berlin Heidelberg, Berlin, Heidelberg, 2008, pp. 53–82, https://doi.org/10.1007/

978-3-540-33265-7 2.

[12] S. Angenent, S. Haker, A. Tannenbaum, and R. Kikinis, On the Laplace-Beltrami operator and brain surface flattening, IEEE Trans. Med. Imaging, 18 (1999), pp. 700–711,https://doi.org/10.1109/42.

796283.

[13] A. Berman and R. Plemmons, Nonnegative Matrices in the Mathematical Sciences, Society for Indus- trial and Applied Mathematics, 1994,https://doi.org/10.1137/1.9781611971262.

[14] B. Bogosel, Efficient algorithm for optimizing spectral partitions, Appl. Math. Comput., 333 (2018), pp. 61–75,https://doi.org/10.1016/j.amc.2018.03.087.

[15] B. Bourdin, D. Bucur, and E. Oudet, Optimal partitions for eigenvalues, SIAM J. Sci. Comput., 31 (2010), pp. 4100–4114,https://doi.org/10.1137/090747087.

[16] A. Chern, U. Pinkall, and P. Schr¨oder, Close-to-conformal deformations of volumes, ACM Trans.

Graph., 34 (2015), pp. 56:1–56:13,https://doi.org/10.1145/2766916.

[17] G. P.-T. Choi and L. M. Lui, A linear formulation for disk conformal parameterization of simply-connected open surfaces, Adv. Comput. Math., (2017), pp. 1–28, https://doi.org/10.1007/