Human-Centered Data Visualization
Angle-Preserving Quadrilateral Mesh Parameterization Wenyong Gong ■ Jilin University, China Xiaohua Xie ■ Sun Yat-Sen University, China Rui Ma and Tieru Wu ■ Jilin University, China
M
esh parameterization, which is a classic issue in computer graphics, is widely used in texture mapping, texture synthesis, surface fitting, remeshing, model repairing, architecture surface design, and industrial garment design. Triangle meshes have been popularly used to represent a surface, so existing parameterization methods are almost all based on triangular meshes.1–4 However, in many situations, quadrilateral meshes are directly used in 3D modeling. Compared with triangular meshes, quadrilateral meshes have many excellent properties. For instance, the quadrilateral has a unique opposite edge for each edge, which makes cutting a quadrilateral mesh more convenient. Furthermore, the edge of quadrilateral mesh can better reflect the edge flow direction of the object surface. Because they are easy to manipulate, quadrilateral meshes are supported by many industrial 3D software applications such as Autodesk’s Maya and 3ds Max and Pixologic’s ZBrush, especially in sculpture modeling applications. Consequently, with the growing use of quadrilateral meshes in practical applications, it is necessary to develop a parameterization method for quadrilateral meshes. (See the “Related Work in Surface Conformal Parameterization” sidebar for a description of earlier works.) To achieve the parameterization for a quadrilateral mesh surface, a straightforward approach is to convert the meshes to triangular meshes and then Published by the IEEE Computer Society
g6xie.indd 51
apply the triangular mesh based parameterization algorithm. This method, however, is fussy and may ignore some of the quadrilateral’s properties. In this article, we present two algorithms for the direct parameterization of quadrilateral meshes with minimal angle distortion: one is for a topological disk surface while the other is for a topological sphere surface. The In response to the growing use parameterization for topological disk surfaces involves devising a of quadrilateral meshes in realdiscrete conformal energy func- world applications, a method tion with a length-preserving for the direct parameterization boundary condition to flatten of quadrilateral meshes is now the quadrilateral meshes. For necessary. The two proposed topological sphere surfaces, we algorithms map a topological derive a Tuette energy function disk surface onto a Euclidean for the initialization of spherical plane and map a topological conformal parameterization on sphere surface onto a unit a quadrilateral mesh and then sphere. define a harmonic energy function to approximate the quadrilateral meshes’ harmonic energy. Finally, we can obtain the parametrization result by minimizing the harmonic energy function.
Planar Conformal Mapping In this section, we devise an angle system of quadrilateral meshes to construct an energy function. By minimizing this energy function with a given length-preserving boundary condition, we obtain a planar discrete conformal mapping for topological disk quadrilateral meshes.
0272-1716/15/$31.00 © 2015 IEEE
IEEE Computer Graphics and Applications
51
10/22/15 1:38 PM
Human-Centered Data Visualization
Related Work in Surface Conformal Parameterization
T
he computer graphics community has extensively studied surface conformal parameterization. Based on the surface topologies, conformal parameterization algorithms fall into three categories: a conformal map of disk-like surfaces, conformal map of sphere-like surfaces, and conformal map of high genus surfaces. In the area of conformal maps of disk-like surfaces, Alla Sheffer and Eric de Sturler proposed unfolding triangular meshes by using angle-based flattening.1 In practice, they minimize a point-wise energy function that measures the angle differences between the original and corresponding planar meshes. Mathieu Desbrun and his colleagues2 and Bruno Levy and his colleagues 3 independently developed a method to linearize the Cauchy-Riemann equation and compute the discrete conformal mappings. Liliya Kharevych and her colleagues presented a method called Circle Patten to compute discrete conformal mappings.4 Li Liu and his colleagues recently presented a method to directly parameterize quadrilateral meshes by minimizing a combinatorial energy.5 For conformal maps of sphere-like surfaces, Steven Haker and his colleagues proposed computing the global conformal parameterization from a topological sphere surface to a unit sphere by solving a linear system.6 Their work is based on the Laplacian-Beltrami operator. Because a harmonic mapping is equivalent to a conformal mapping for genus zero closed surfaces, Xianfeng Gu and his colleagues proposed a nonlinear method and iteratively computed conformal mappings by minimizing the harmonic energy.7,8 Lastly, for conformal maps of high genus surfaces, it is impossible to generalize the harmonic mapping method to high genus surfaces because the harmonic mapping is not equivalent to the conformal mapping for high genus cases. Based on the Hodge theory, Xianfeng Gu and Shing-Tung Yau proposed using the simplicial cohomology to approximate the De Rham cohomology and to com-
pute a basis of holomophic 1-forms.8 They then computed the conformal structures for high genus meshes. Andrei Khodakovsky and his colleagues proposed a method to compute a globally smooth parameterization with low distortion.9
References 1. A. Sheffer and E. de Sturler, “Parameterization of Faceted Surfaces for Meshing Using Angle Based Flattening,” Engineering with Computers, vol. 17, no. 3, 2001, pp. 326–337. 2. M. Desbrun, M. Meyer, and P. Alliez, “Intrinsic Parameterizations of Surface Meshes,” Computer Graphics Forum, vol. 21, no. 3, 2002, pp. 209–218. 3. B. Levy et al., “Least Squares Conformal Maps for Automatic Texture Atlas Generation,” ACM Trans. Graphics, vol. 21, no. 3, 2002, pp. 362–371. 4. L. Kharevych, B. Springborn, and P. Schröder, “Discrete Conformal Mappings via Circle Patterns,” ACM Trans. Graphics, vol. 25, no. 2, 2006, pp. 412–438. 5. L. Liu, C.M. Zhang, and F. Cheng, “Parametrization of Quadrilateral Meshes,” Proc. 7th Int’l Conf. Computational Science (ICCS), Part II, LNCS 4488, Springer, 2007, pp. 17–24. 6. S. Haker et al., “Conformal Surface Parameterization for Texture Mapping,” IEEE Trans. Visualization and Computer Graphics, vol. 6, no. 2, 2000, pp. 181–189. 7. X. Gu et al., “Genus Zero Surface Conformal Mapping and Its Application to Brain Surface Mapping,” IEEE Trans. Medical Imaging, vol. 23, no. 8, 2004, pp. 949–958. 8. X. Gu and S.T. Yau, “Computing Conformal Structures of Surfaces,” Comm. Information and Systems, vol. 2, no. 2, 2002, pp. 121–146. 9. A. Khodakovsky, N. Litke, and P. Schröder, “Globally Smooth Parameterizations with Low Distortion,” ACM Trans. Graphics, vol. 22, no. 3, 2003, pp. 350–357.
Angle System of Quadrilateral Meshes For the mesh parameterization, the ideal solution is isometric (area-preserving and angle-preserving), but this is impossible in most cases. Motivated by circle patterns5 here we establish an angle system that approximates the angles of a quadrilateral mesh as much as possible. Consider a planar quadrilateral mesh M = {V, E, Q} of finitely many points U = {ui} in the plane, where V = {v i}, E = {e i}, and Q = {qi} represent the sets of vertices, edges, and quadrilaterals, respectively, and ui is the point position of vertex v i. is a quadrilateral mesh in R3, we look forIf M ward to seeking a corresponding planar quadrilat52
g6xie.indd 52
eral mesh M with angles that approximate those as close as possible. Therefore, we minimize of M the following energy function: E (θ) =
∑θ
k j , j+1, j+2
− θ jk, j+1, j+2
2
,
which is subject to the following constraints:
■■
positivity, ∀θjk, j+1, j+2 : θjk, j+1, j+2 > 0 ; quadrilateral sum condition, ∀q = { vk v j v j+1v j+2 } ∈ Q :
■■
θjk, j+1, j+2 + θkj , j+1, j+2 + θkj+, j ,1j+2 + θkj+, j ,2j+1 = 2π ; and
■■
■■
interior vertex sum condition, ∀vi ∈ Vint : vk ∈qk θi = 2π ,
∑
November/December 2015
10/22/15 1:38 PM
vj vk
where the notations refer to Figure 1. Obviously, the first two constraints imply that q < 2p for arbitrary angle q. The boundary condition is important because it determines the shape of the parameterization mesh. In this article, we only consider the natural boundary conditions, thus the following constraint on the boundary vertex (see Figure 1) is needed:
vj+1 vj+m vj+2 vj+6
vi vj+3
∀vk ∈ Vbdy :
∑
k vk ∈qk , j , j +1, j +2 q j , j+1, j+2
< 2p vj+4
The angle system is therefore reduced an energy minimizing problem with 4|Q| variables, |Q| + |Vint| equality constraints and 4|Q| + |Vbdy| inequality constraints. As a quadric constraint optimization problem, this problem can be solved by many optimization software applications, such as MOSEK.
Figure 1. Boundary vertex vk and interior vertex vi on a quadrilateral mesh. The constraint is on the boundary vertex because we only consider the natural boundary conditions in this work. vj+1
vi
Free Boundary Planar Discrete Conformal Mapping
l r l r cot αij + cot αij + cot βij + cot βij u − u Ep = j i 4 4 4 4 [ vi v j ]∈E
∑
γij′ η ′ 2 + + cot ij ui − u j′ , cot 4 4 vi v j ′ ∈Ed
g6xie.indd 53
vi vj ′
β lij
ηij ′
β rij
vj
2
Figure 2. Edges, diagonals, and angles in the discrete conformal energy calculation. Both edges and diagonals are considered in our energy function, forming a new angle system.
where E and Ed represent the edge and diagonal sets, respectively. Figure 2 describes cot angles in the energy function. Because coefficients of the energy function are always positive, the local injectivity of the discrete conformal mapping can be guaranteed. Our energy function is similar to the one used in previous work,6 but we modified it in two ways: we consider the diagonals and use the angles calculated in the last section. An obvious advantage of using these angles is we can efficiently avoid edges overlapping each other. The quadric energy function only depends on the angles, so the conformality of the discrete mapping can be preserved better by using the angles we described in the last section instead of the angles from the original surface. To flatten a mesh satisfying minimal angle distortions, the minimum discrete conformal energy is required. Thus, we just need a vanishing derivative of Ep with respect to ui:
r ij
α lij
∑
γij ′ α
Conformal mapping is angle-preserving. To unfold a mesh into a plane with minimal angle distortions between the original mesh and the corresponding planar mesh, we devise a discrete conformal energy function that is motivated by the previous version defined on triangular meshes.1
vj+5
r l r αl ∂E p cot ij + cot αij + cot βij + cot βij u − u =2 j) ( i ∂ui 4 4 4 4 v j ∈dN(i)
∑
η ′ γij′ , +2 + cot ij (ui − u j′ ) cot 4 4 v j ′ ∈iN(i)
∑
(1)
where angles have the same meaning as before and dN(i) represents the points that directly connect with vi. Thus we also call dN(i) the direct onering neighborhood of vi. On the other hand, iN(i) denotes the points that share a common quadrilateral with vi but do not connect with vi. We call iN(i) the indirect one-ring neighborhood of vi. From Equation 1, we get alinear equations with a coefficient matrix that is sparse and coefficients that are related to angles we established earlier. To solve this equation, we still need to prescribe the boundary condition. The natural boundary condition is used when we calculate the angle system, so we map the boundary of a mesh into a plane by using a length-preserving free boundary manner. IEEE Computer Graphics and Applications
53
10/22/15 1:38 PM
Human-Centered Data Visualization Pi
Pj+1
Pj –1
Pi
Pj+1
Pi Pj ′
Pj ′ (a)
(b)
Pj (c)
βj1
Pj+1
α j4
αj2
Pj ′
Pj
Pj
α j1
Pj
Pi
Pi
βj2
Pj
α j3
Pj ′ (d)
(e)
Pj+1
Pj ′
Figure 3. Delaunay and anti-Delaunay triangulation. (a) One-ring neighbors of a point Pi. (b) Delaunay triangulation of one-ring neighbors of Pi on a quadrilateralmesh. (c) Anti-Delaunay triangulation of one-ring neighbors of Pi on a quadrilateral mesh. (d) Anti-Delaunay triangulation of a quadrilateral. (e) Delaunay triangulation of a quadrilateral.
Spherical Discrete Conformal Mapping Here we focus on the genus zero closed quadrilateral meshes, whose natural parameterization domain is the unit sphere in R3.
Laplacian of Quadrilateral Meshes In this subsection, we investigate the Laplacian on quadrilateral meshes by taking advantage of the calculation formulation on triangular meshes.7 Let M ∈ R3 be a two-manifold triangular mesh. Given an arbitrary vertex P i of M, its discrete Laplacian can be formulated as D M Pi = 2
∑w
ij
(Pi − Pj ) ,
j∈N(i)
where N(i) is the index set of the one-ring neighboring points of P i. Mathieu Desbrun and his colleagues used the COT scheme to calculate wij.7 To calculate the discrete differential operator, we first introduce the Delaunay triangulation and anti-Delaunay triangulation of the one-ring neighbors of a point P i on a quadrilateral mesh. The one-ring neighbors of a point P i on a quadrilateral mesh form a set of quadrilaterals that share the same vertex P i. Figure 3a provides an illustration. Delaunay triangulation means every triangle in a triangulation contains no interior points, and it maximizes the minimum angles. A single quadrilateral Delaunay triangulation satisfies the summation of opposite angles is less than 180 degrees (aj1 + aj4 < 180 degrees as in Figure 3d). For simplicity, we denote the Delaunay triangulation and anti-Delaunay triangulation of a one-ring neighbors of a point on a quadrilateral mesh as 1-DTQ and 1-aDTQ, respectively. Now we give the formal definitions of both. Definition 1 (1-DTQ). 1-DTQ refers to every one-ring neighboring quadrilateral of a point P that is a single quadrilateral Delaunay triangulation. 54
g6xie.indd 54
Definition 2 (1-aDTQ). 1-aDTQ refers to every onering neighboring quadrilateral of a point P that is not a single quadrilateral Delaunay triangulation. 1-DTQ and 1-aDTQ represent the good and bad triangulation. Furthermore, the Laplacian is estimated with the aid of 1-DTQ and 1-aDTQ. Figures 3b and 3c illustrate the 1-DTQ and 1-aDTQ, respectively. According to the COT scheme,7 we now derive the Laplacian of P i on quadrilateral meshes. In 1-DTQ and 1-aDTQ, some triangles (such as DP jP j′P j+1 in Figure 3c) are not used when computing the Laplacian. Therefore, vertices such as P j′ are only used once, while vertices such as P j are used D twice exactly. Let AM (Pi ) and AMa (Pi ) represent the area of 1-DTQ and 1-aDTQ, respectively. Then we give the area estimation Q AM (Pi ) =
1 D 1 a AM (Pi ) + AM (Pi ). 2 2
Consequently, the Laplacian of P i can be formulated as DQM (Pi ) =
(Pi )
1 Q 2 AM
∑w j
ij
(Pi − Pj ) − ∑ w ij′ (Pi − Pj′ ) , j′
(2) where wij = (cotaj2 + cotaj–1,3 + cotbj2 + cotbj–1,1)/2, wij′ = (cotaj1 + cota2)/2. Figure 3d provides descriptions of aik and bik.
Harmonic Energy of Quadrilateral Meshes Because conformal mapping is angle-preserving, it has extensive applications in computer graphics and medical image processing. 2 According to the harmonic analysis theory,8 a mapping is harmonic if it minimizes the harmonic energy. Moreover, a degree one harmonic map h: M → S2 is conformal. Suppose f is a function defined on a surface M. Then the harmonic energy of f is
November/December 2015
10/22/15 1:38 PM
Pk
Pi
E(f) = ∫M ||∆f||2 ds,
Pi
where ∇ is the gradient operator and ds is the area element. When M is a triangular mesh and f is a piecewise linear function defined on M, then E(f) has the following discretization:
∑
E( f ) =
2
g i1 = (1 – u)(1 – v)f i + (1 – u)vf j + u(1 – v)f j′ + uvf j+1 , g i2 = (1 – u)(1 – v)f j + (1 – u)vf j′ + u(1 – v)f j+1 + uvf i , g i3 = (1 – u)(1 – v)f j′ + (1 – u)vf j+1 + u(1 – v)f i + uvf j , and g i4 = (1 – u)(1 – v)f j+1 + (1 – u)vf i + u(1 – v)f j + uvf j′ , where u, v ∈ [0, 1]. According to the definition of the harmonic energy, we have the following energy of g i1 defined on Qi: 1 2 ( f i − f j′ ) 3 1 1 1 2 2 2 + ( f i − f j ) + ( f j − f j+1 ) + ( f j+1 − f j′ ) 3 3 3 1 1 + ( f i − f j′ )( f j − f j+1 ) − ( f i − f j )( f j+1 − f j′ ) . 3 3 1
∫ ∫ 0
1
0
2
∇g i1 du dv =
The E(g i2), E(g i3), and E(g i4) have similar representations. Let f=
1 ( g i1 + g i2 + g i3 + g i 4 ) . 4
Then the energy of f defined on Qi can be formulated as Ei ( f ) =
(
1 2 2 2 4 ( f i − f j′ ) + 4 ( f j − f j+1 ) + 2( f i − f j ) 12 2
2
2
+2( f j+1 − f j′ ) + 2( f i − f j′ ) + 2( f i − f j+1 )
)
∑ 6 E ( f ) = ∑ (2 ( f − f ) i
i
j′
2
2
2
)
+ ( f i − f j ) + ( f j+1 − f j′ ) + ( f i − f j′ ) + ( f i − f j+1 ) .
g6xie.indd 55
(b)
Figure 4. Two kinds of considered edges in the harmonic energy calculation. (a) The edge [PiPj] shared by two adjacent quadrilaterals. (b) The diagonal edge [PiPj′] dividing a quadrilateral into two quadrilaterals Ta, Tb. Their weights for harmonic energy calculation are the summation of the cotangents of related angles.
According to the Tuette energy, we also could define a Laplacian, and we call it the Tuette-Laplacian, which is defined as the derivative of ET(f) with respect to f. ∆T f (Pi ) =
∑
4 ( fi − f j ) +
[ Pj Pj ]∈E
∑
4 ( f i − f j′ ) ,(4)
Pj Pj ′ ∈Ed
where Ed is the set of diagonal edges. With reference to the relationship between the Tuette energy and the harmonic energy defined in a triangular mesh, we define the harmonic energy of quadrilateral meshes as follows:
∑
2
∑
2
EH ( f ) = kPi ,Pj f i − f j + kPi ,Pj ′ f i − f j′ . Pi Pj ′ ∈Ed [ PiPj ]∈E (5) Suppose [P iP j] has two adjacent quadrilaterals Qa = {P iP jP j′P j+1} and Qb = {P iP jPk′Pk} (see Figure 4a). Let ∠PiPj Pj′ ∠Pi Pj′ Pj+1 + cot 4 4 ∠Pj′ Pj+1Pi ∠Pj+1PiPj + cot + cot . 4 4
k aPiPj = cot
Similarly, kPbi ,Pj is defined. Thus, kPi ,Pj = kPai ,Pj + kPbi ,Pj . In addition, a diagonal edge [P iP j′] divides a quadrilateral into two triangles T a, T b (see Figure 4b), then we define ∠PiPj Pj′ ∠Pi Pi+1Pj′ + cot . 4 4
Spherical Discrete Conformal Mapping
2
+ 2( f j − f j+1 ) 2
(a)
kPi ,Pj ′ = cot
2
2
Pj ′
Pj ′
.
Therefore, we define the Tuette energy of f on M as follows: ET ( f ) =
Pj
Pj
where u and v are any two vertices of M, [u, v] is an edge of M, and ku,v is the coefficient that is related to the two adjacent triangles of [u, v]. The energy is known as the Tuette energy, especially when ku,v = 1. Suppose f is a function defined on a quadrilateral mesh M, and f i is the function value of vertex P i. For given quads Qi = {P iP jP j′P j+1} (see Figure 3d), we use the following formulations to approximate f on the quads Qi:
E ( g i1 ) =
Tα Qα
ku ,v f (u) − f ( v ) ,
[ u ,v ]∈M
Tβ
Pi+1
Qβ
Pk ′
Pi+1
(3)
For a genus zero closed mesh M1 and a unit spherical mesh S2, there are many conformal mappings between them. That is, the conformal mapping between two meshes is not unique. The basis of IEEE Computer Graphics and Applications
55
10/22/15 1:38 PM
Human-Centered Data Visualization
Input: A quadrilateral mesh M1, time step dt, threshold e Output: A spherical Tuette mapping h 1. Compute the Guassian map n: M1 → S2 . Let h = n, compute the Tuette energy ET (h). 2. Update the map: h = h – dt × DTh. 3. Compute the Tuette energy ET′ (h) . 4. If δE = ET′ (h) − ET (h) < ε , return h. Otherwise let ET (h) = ET′ (h) , and repeat steps 2–4.
component DN(h) and the tangential component DT(h)—that is, D(h) = DN(h) + DT(h). To ensure the unique solution, an extra constraint is required. The constraint is that the mass center of S2 is coincided with the origin:
∫
S2
h dσ M1 = 0 ,(9)
Figure 5. The algorithm flow chart for building Spherical Tuette mapping. The key of the algorithm is to compute the homeomorphism h by minimizing the Tuette energy in Equation 3.
where σ M1 is the area element. Then we can construct the geometric partial differential equation
Input: A quadrilateral mesh M1, time step dt, threshold e Output: A spherical conformal mapping h 1. Compute the spherical Tuette map ht: M1 → S2 . Let h = ht, compute the harmonic energy Eh = EH (h). 2. Update the map: h = h – dt × DTh. 3. Compute a Möbius transformation T, such that T ° h satisfies the constraint in Equation 9 of the mass center. 4. Compute the harmonic energy Eh′ = EH (h) . 5. If δE = Eh′ − E h < ε , return h. Otherwise let , and repeat steps 2–5.
with the constraint of Equation 9. The steady state solution h of the geometric flow is the conformal mapping from M1 to S2. Equation 10 can be efficiently solved by the steepest descendent algorithm. The basic idea of computing the discrete spherical conformal mapping has two steps: First, we compute a homeomorphism h from M1 to S2, and then we diffuse h to be harmonic. To compute the homeomorphism h between M1 and S2, we construct the spherical Tuette mapping as in Figure 5, which minimizes the Tuette energy defined in Equation 3. Viewing the spherical Tuette mapping as the initial homeomorphism, we construct the spherical conformal mapping for genus zero surfaces as in Figure 6.
Figure 6. The algorithm flow chart for building Spherical harmonic mapping. The spherical Tuette mapping is used as the initialization.
computing conformal mappings is that degree one harmonic mappings are conformal for genus zero surfaces.9 For a map h : M1 → S2 where h = (h0, h1, h2), the Tuette energy Er(h) of h can be defined as ET ( f ) =
2
∑E
T
(hi ) ,(6)
i=0
where ET(hi) is the same as in Equation 3. Also, the harmonic energy EH(h) of h is given by EH (h) =
2
∑E
H
(hi ) ,(7)
i=0
where EH(hi) is defined as in Equation 5. The Laplacian of h is D(h) = (D(h0), D(h1), D(h2)),(8) where D is the Laplacian operator in Equations 2 or 4. In practice, we use Equation 4 for spherical Tuette mapping (see Figure 5), whereas we use Equation 2 for spherical harmonic mapping (see Figure 6). Because a map h : M1 → S2 is harmonic if and only if D(h) has vanishing tangential component, h could be written as the sum of the normal 56
g6xie.indd 56
∂h = −∆T h (10) ∂t
Experimental Results In this section, we demonstrate the proposed quadrilateral mesh parameterization methods on three 3D models—Nicolo, Venus, and Moai—by displaying the texture mapping results based on specific parameterization manners. We also report the quantification evaluation based on the angle distortion. Some related parameterization methods were also carried out for comparison. Figure 7 illustrates the parameterization results on the Nicolo model when using different disk parameterization methods based on triangular meshes: least squares conformal Maps (LSCM),10 intrinsic,1 and low-stretch methods.11 To execute the triangular mesh based parameterization algorithms, we converted the quadrilateral Nicolo model (left image in Figure 7a) into the triangle mesh model (left image in Figure 7c) by connecting diagonals. As shown, our method can preserve the angles better than other methods. Furthermore, because our method adopts the free boundary condition, it is in a fair way to attain a very small change in the boundary length. In the experiments, the maximum boundary length ratios between the parameterized meshes and the original meshes are below 3 percent.
November/December 2015
10/22/15 1:38 PM
(a)
(b)
(c) Figure 7. The disk parameterization on the Nicolo model. (a) The quadrilateral Nicolo model and its texture mapping result achieved by our parameterization method. (b) The corresponding disk parameterization domain of the Nicolo model. (c) The triangulated Nicolo model and the parameterization results generated by least squares conformal maps (LSCM),10 intrinsic parameterization,1 and low-stretch parameterization methods,11 respectively, from left to right.
The spherical parameterization was implemented on the Venus and Moai models. Figure 8 shows the texture mapping results on Venus model based on different parameterization methods. For the sake of comparison, the figure shows the result on the quadrangular mesh when using the Maya software and our method and the result on the triangulated mesh when using the conformal spherical parameterization algorithm. 2 Our method attains better angle preservation
(a)
(b)
than the Maya software (especially in the neck region), and our result are comparable with the Gu parameterization method.2 The quadrangular Venus and Moai models were also parameterized onto unit spheres using several methods. As shown in Figures 8 and 9, our method attains slightly better result than the Maya software, but its results are worse than the Gu triangular mesh based parameterization method. 2
(c)
Figure 8. The spherical parameterization on the Venus model. (a) The Venus model and the parameterization result by our method. (b) The texture mapping results when using Maya software (left) and our method (right). (c) The triangulated Venus model. (d) The texture mapping based on the Gu parameterization method2 on the triangular mesh.
g6xie.indd 57
IEEE Computer Graphics and Applications
57
10/22/15 1:38 PM
Human-Centered Data Visualization
(a)
(b)
(c)
Figure 9. The spherical parameterization on the Moai model. (a) The Moai model. (b) The spherical parameterization domain of the Moai model (left), the texture mapping results based on the spherical parameterization by using Maya software (middle), and our method (right). (c) The triangulated Moai model. (d) The texture mapping result based on the Gu parameterization method on the triangular mesh.2
Table 1. Angle distortion of different parameterization methods. Model
Planar parameterization methods
Spherical parameterization methods
Low stretch
Intrinsic
LSCM
Ours
Gu
Nicolo
0.9312
0.3230
0.2995
0.0879
–
Venus
–
–
–
–
0.0025
0.02175
Moai
–
–
–
–
0.074
0.24
To quantify the parameterization distortion, we use the following function to measure the angle distortion: TA =
n
k
∑∑ f =1 i=1
σp 2 ρ f io − 1 ,(11) σi
where n is the face number and k = 3(4) when the mesh is triangular (quadrilateral). The σip and σio are the angle after parameterization and the original angle, respectively. The weight rf is defined as rf = Af/Sf A f, where Af is the area of the triangle or quadrilateral f. Table 1 gives the angle distortions of different models with respect to various methods, including the low-stretch,11 intrinsic parameterization,1 LSCM,10 Gu conformal spherical parameterization,2 and our planar and spherical parameterization. The table shows that our spherical parameterization method gets larger distortion values than the Gu method,2 but our planar parameterization method attains the smallest distortion value among all the compared planar parameterization algorithms. These examples all have demonstrated the feasibility and practicability of our parameterization methods for quadrilateral meshes. For the planar case, unlike existing parameterization methods that 58
g6xie.indd 58
Ours –
simultaneously determine the angles and boundaries,1,10,11 our method specifies the angle system in advance so that it has better conformal capability. For the spherical case, because the harmonic mapping on a quadrilateral mesh is not well-defined in mathematics, our method in fact attains slightly worse performance compared with the triangle mesh based parameterization methods.2
W
ith the growing use of quadrilateral meshes in 3D modeling, direct parameterization methods for quadrilateral meshes are required. In the future, we will investigate the angle-preserving parameterization globally for high genus quadrilateral meshes. Furthermore, we would like to extend the proposed scheme to the parameterization of a point cloud.
Acknowledgments We thank Joel Daniels II from the New York University Polytechnic Institute for providing quad meshes (Moai and Venus quad models). This work is supported by the National Nature Science Foundation of China (grants 61202223 and 61373003) and the GuangZhou Program (grant 201508010032).
November/December 2015
10/22/15 1:38 PM
References 1. M. Desbrun, M. Meyer, and P. Alliez, “Intrinsic Parameterizations of Surface Meshes,” Computer Graphics Forum, vol. 21, no. 3, 2002, pp. 209–218. 2. X. Gu et al., “Genus Zero Surface Conformal Mapping and Its Application to Brain Surface Mapping,” IEEE Trans. Medical Imaging, vol. 23, no. 8, 2004, pp. 949–958. 3. A. Khodakovsky, N. Litke, and P. Schröder, “Globally Smooth Parameterizations with Low Distortion,” ACM Trans. Graphics, vol. 22, no. 3, 2003, pp. 350–357. 4. A. Sheffer and E. de Sturler, “Parameterization of Faceted Surfaces for Meshing Using Angle Based Flattening,” Engineering with Computers, vol. 17, no. 3, 2001, pp. 326–337. 5. L. Kharevych, B. Springborn, and P. Schröder, “Discrete Conformal Mappings via Circle Patterns,” ACM Trans. Graphics, vol. 25, no. 2, 2006, pp. 412–438. 6. L. Liu, C.M. Zhang, and F. Cheng, “Parametrization of Quadrilateral Meshes,” Proc. 7th Int’l Conf. Computational Science (ICCS), Part II, LNCS 4488, Springer, 2007, pp. 17–24. 7. M. Desbrun et al., “Implicit Fairing of Irregular Meshes Using Diffusion and Curvature Flow,” Proc. 26th Ann. Conf. Computer Graphics and Interactive Techniques (SIGGRAPH), 1999, pp. 317–324. 8. R. Schoen and S.T. Yau, Lectures on Harmonic Maps, Harvard Univ. Press, 1997. 9. X. Gu and S.T. Yau, “Computing Conformal Structures of Surfaces,” Comm. Information and Systems, vol. 2, no. 2, 2002, pp. 121–146. 10. B. Levy et al., “Least Squares Conformal Maps for Automatic Texture Atlas Generation,” ACM Trans. Graphics, vol. 21, no. 3, 2002, pp. 362–371.
11. S. Yoshizawa, A. Belyaev, and H.P. Seidel, “A Fast and Simple Stretch-Minimizing Mesh Parameterization,” Proc. Shape Modeling and Application, 2004, pp. 200–208. Wenyong Gong is a senior algorithm engineer at AEE in Shenzhen, China. His research interests include computer graphics, image reconstruction, and image recognition. Gong has a PhD from the Institute of Mathematics at Jilin University. Contact him at [email protected]. Xiaohua Xie (corresponding author) is a research professor at Sun Yat-Sen University. His research interests include image processing, computer vision, pattern recognition, and computer graphics, especially focusing on image understanding and object modeling. Xie has a PhD in applied mathematics from Sun Yat-Sen University. He received the Overseas High-Caliber Personnel Award in Shenzhen in 2012. Contact him at [email protected]. Rui Ma is a PhD candidate in the School of Computing Science at Simon Fraser University. His research interests include computer graphics, high-level geometry processing, and shape analysis. Ma has an MS in mathematics from Jilin University. Contact him at [email protected]. Tieru Wu is a professor in the Institute of Mathematics at Jilin University, China. He research interests include developing various techniques of computer graphics, geometry processing, machine learning in image processing, and image retrieval. His research is supported in part by the National Science Foundation of Chinese. Wu has a PhD in quantitative economics from Jilin University. Contact him at [email protected].
EXECUTIVE STAFF PURPOSE: The IEEE Computer Society is the world’s largest association of computing professionals and is the leading provider of technical information in the field. MEMBERSHIP: Members receive the monthly magazine Computer, discounts, and opportunities to serve (all activities are led by volunteer members). Membership is open to all IEEE members, affiliate society members, and others interested in the computer field. COMPUTER SOCIETY WEBSITE: www.computer.org
Next Board Meeting: 15–16 November 2015, New Brunswick, NJ, USA EXECUTIVE COMMITTEE President: Thomas M. Conte President-Elect: Roger U. Fujii; Past President: Dejan S. Milojicic; Secretary: Cecilia Metra; Treasurer, 2nd VP: David S. Ebert; 1st VP, Member & Geographic Activities: Elizabeth L. Burd; VP, Publications: Jean-Luc Gaudiot; VP, Professional & Educational Activities: Charlene (Chuck) Walrad; VP, Standards Activities: Don Wright; VP, Technical & Conference Activities: Phillip A. Laplante; 2015–2016 IEEE Director & Delegate Division VIII: John W. Walz; 2014–2015 IEEE Director & Delegate Division V: Susan K. (Kathy) Land; 2015 IEEE Director-Elect & Delegate Division V: Harold Javid
BOARD OF GOVERNORS Term Expiring 2015: Ann DeMarle, Cecilia Metra, Nita Patel, Diomidis Spinellis, Phillip A. Laplante, Jean-Luc Gaudiot, Stefano Zanero Term Expriring 2016: David A. Bader, Pierre Bourque, Dennis J. Frailey, Jill I. Gostin, Atsuhiro Goto, Rob Reilly, Christina M. Schober Term Expiring 2017: David Lomet, Ming C. Lin, Gregory T. Byrd, Alfredo Benso, Forrest Shull, Fabrizio Lombardi, Hausi A. Muller
g6xie.indd 59
Executive Director: Angela R. Burgess; Director, Governance & Associate Executive Director: Anne Marie Kelly; Director, Finance & Accounting: Sunny Hwang; Director, Information Technology Services: Ray Kahn; Director, Membership: Eric Berkowitz; Director, Products & Services: Evan M. Butterfield; Director, Sales & Marketing: Chris Jensen
COMPUTER SOCIETY OFFICES Washington, D.C.: 2001 L St., Ste. 700, Washington, D.C. 20036-4928 Phone: +1 202 371 0101 • Fax: +1 202 728 9614 • Email: [email protected] Los Alamitos: 10662 Los Vaqueros Circle, Los Alamitos, CA 90720 Phone: +1 714 821 8380 • Email: [email protected] MEMBERSHIP & PUBLICATION ORDERS Phone: +1 800 272 6657 • Fax: +1 714 821 4641 • Email: [email protected] Asia/Pacific: Watanabe Building, 1-4-2 Minami-Aoyama, Minato-ku, Tokyo 107-0062, Japan • Phone: +81 3 3408 3118 • Fax: +81 3 3408 3553 • Email: [email protected]
IEEE BOARD OF DIRECTORS President & CEO: Howard E. Michel; President-Elect: Barry L. Shoop; Past President: J. Roberto de Marca; Director & Secretary: Parviz Famouri; Director & Treasurer: Jerry Hudgins; Director & President, IEEE-USA: James A. Jefferies; Director & President, Standards Association: Bruce P. Kraemer; Director & VP, Educational Activities: Saurabh Sinha; Director & VP, Membership and Geographic Activities: Wai-Choong Wong; Director & VP, Publication Services and Products: Sheila Hemami; Director & VP, Technical Activities: Vincenzo Piuri; Director & Delegate Division V: Susan K. (Kathy) Land; Director & Delegate Division VIII: John W. Walz revised 5 June 2015
IEEE Computer Graphics and Applications
59
10/22/15 1:38 PM
Angle-Preserving Quadrilateral Mesh Parameterization Wenyong Gong ■ Jilin University, China Xiaohua Xie ■ Sun Yat-Sen University, China Rui Ma and Tieru Wu ■ Jilin University, China
M
esh parameterization, which is a classic issue in computer graphics, is widely used in texture mapping, texture synthesis, surface fitting, remeshing, model repairing, architecture surface design, and industrial garment design. Triangle meshes have been popularly used to represent a surface, so existing parameterization methods are almost all based on triangular meshes.1–4 However, in many situations, quadrilateral meshes are directly used in 3D modeling. Compared with triangular meshes, quadrilateral meshes have many excellent properties. For instance, the quadrilateral has a unique opposite edge for each edge, which makes cutting a quadrilateral mesh more convenient. Furthermore, the edge of quadrilateral mesh can better reflect the edge flow direction of the object surface. Because they are easy to manipulate, quadrilateral meshes are supported by many industrial 3D software applications such as Autodesk’s Maya and 3ds Max and Pixologic’s ZBrush, especially in sculpture modeling applications. Consequently, with the growing use of quadrilateral meshes in practical applications, it is necessary to develop a parameterization method for quadrilateral meshes. (See the “Related Work in Surface Conformal Parameterization” sidebar for a description of earlier works.) To achieve the parameterization for a quadrilateral mesh surface, a straightforward approach is to convert the meshes to triangular meshes and then Published by the IEEE Computer Society
g6xie.indd 51
apply the triangular mesh based parameterization algorithm. This method, however, is fussy and may ignore some of the quadrilateral’s properties. In this article, we present two algorithms for the direct parameterization of quadrilateral meshes with minimal angle distortion: one is for a topological disk surface while the other is for a topological sphere surface. The In response to the growing use parameterization for topological disk surfaces involves devising a of quadrilateral meshes in realdiscrete conformal energy func- world applications, a method tion with a length-preserving for the direct parameterization boundary condition to flatten of quadrilateral meshes is now the quadrilateral meshes. For necessary. The two proposed topological sphere surfaces, we algorithms map a topological derive a Tuette energy function disk surface onto a Euclidean for the initialization of spherical plane and map a topological conformal parameterization on sphere surface onto a unit a quadrilateral mesh and then sphere. define a harmonic energy function to approximate the quadrilateral meshes’ harmonic energy. Finally, we can obtain the parametrization result by minimizing the harmonic energy function.
Planar Conformal Mapping In this section, we devise an angle system of quadrilateral meshes to construct an energy function. By minimizing this energy function with a given length-preserving boundary condition, we obtain a planar discrete conformal mapping for topological disk quadrilateral meshes.
0272-1716/15/$31.00 © 2015 IEEE
IEEE Computer Graphics and Applications
51
10/22/15 1:38 PM
Human-Centered Data Visualization
Related Work in Surface Conformal Parameterization
T
he computer graphics community has extensively studied surface conformal parameterization. Based on the surface topologies, conformal parameterization algorithms fall into three categories: a conformal map of disk-like surfaces, conformal map of sphere-like surfaces, and conformal map of high genus surfaces. In the area of conformal maps of disk-like surfaces, Alla Sheffer and Eric de Sturler proposed unfolding triangular meshes by using angle-based flattening.1 In practice, they minimize a point-wise energy function that measures the angle differences between the original and corresponding planar meshes. Mathieu Desbrun and his colleagues2 and Bruno Levy and his colleagues 3 independently developed a method to linearize the Cauchy-Riemann equation and compute the discrete conformal mappings. Liliya Kharevych and her colleagues presented a method called Circle Patten to compute discrete conformal mappings.4 Li Liu and his colleagues recently presented a method to directly parameterize quadrilateral meshes by minimizing a combinatorial energy.5 For conformal maps of sphere-like surfaces, Steven Haker and his colleagues proposed computing the global conformal parameterization from a topological sphere surface to a unit sphere by solving a linear system.6 Their work is based on the Laplacian-Beltrami operator. Because a harmonic mapping is equivalent to a conformal mapping for genus zero closed surfaces, Xianfeng Gu and his colleagues proposed a nonlinear method and iteratively computed conformal mappings by minimizing the harmonic energy.7,8 Lastly, for conformal maps of high genus surfaces, it is impossible to generalize the harmonic mapping method to high genus surfaces because the harmonic mapping is not equivalent to the conformal mapping for high genus cases. Based on the Hodge theory, Xianfeng Gu and Shing-Tung Yau proposed using the simplicial cohomology to approximate the De Rham cohomology and to com-
pute a basis of holomophic 1-forms.8 They then computed the conformal structures for high genus meshes. Andrei Khodakovsky and his colleagues proposed a method to compute a globally smooth parameterization with low distortion.9
References 1. A. Sheffer and E. de Sturler, “Parameterization of Faceted Surfaces for Meshing Using Angle Based Flattening,” Engineering with Computers, vol. 17, no. 3, 2001, pp. 326–337. 2. M. Desbrun, M. Meyer, and P. Alliez, “Intrinsic Parameterizations of Surface Meshes,” Computer Graphics Forum, vol. 21, no. 3, 2002, pp. 209–218. 3. B. Levy et al., “Least Squares Conformal Maps for Automatic Texture Atlas Generation,” ACM Trans. Graphics, vol. 21, no. 3, 2002, pp. 362–371. 4. L. Kharevych, B. Springborn, and P. Schröder, “Discrete Conformal Mappings via Circle Patterns,” ACM Trans. Graphics, vol. 25, no. 2, 2006, pp. 412–438. 5. L. Liu, C.M. Zhang, and F. Cheng, “Parametrization of Quadrilateral Meshes,” Proc. 7th Int’l Conf. Computational Science (ICCS), Part II, LNCS 4488, Springer, 2007, pp. 17–24. 6. S. Haker et al., “Conformal Surface Parameterization for Texture Mapping,” IEEE Trans. Visualization and Computer Graphics, vol. 6, no. 2, 2000, pp. 181–189. 7. X. Gu et al., “Genus Zero Surface Conformal Mapping and Its Application to Brain Surface Mapping,” IEEE Trans. Medical Imaging, vol. 23, no. 8, 2004, pp. 949–958. 8. X. Gu and S.T. Yau, “Computing Conformal Structures of Surfaces,” Comm. Information and Systems, vol. 2, no. 2, 2002, pp. 121–146. 9. A. Khodakovsky, N. Litke, and P. Schröder, “Globally Smooth Parameterizations with Low Distortion,” ACM Trans. Graphics, vol. 22, no. 3, 2003, pp. 350–357.
Angle System of Quadrilateral Meshes For the mesh parameterization, the ideal solution is isometric (area-preserving and angle-preserving), but this is impossible in most cases. Motivated by circle patterns5 here we establish an angle system that approximates the angles of a quadrilateral mesh as much as possible. Consider a planar quadrilateral mesh M = {V, E, Q} of finitely many points U = {ui} in the plane, where V = {v i}, E = {e i}, and Q = {qi} represent the sets of vertices, edges, and quadrilaterals, respectively, and ui is the point position of vertex v i. is a quadrilateral mesh in R3, we look forIf M ward to seeking a corresponding planar quadrilat52
g6xie.indd 52
eral mesh M with angles that approximate those as close as possible. Therefore, we minimize of M the following energy function: E (θ) =
∑θ
k j , j+1, j+2
− θ jk, j+1, j+2
2
,
which is subject to the following constraints:
■■
positivity, ∀θjk, j+1, j+2 : θjk, j+1, j+2 > 0 ; quadrilateral sum condition, ∀q = { vk v j v j+1v j+2 } ∈ Q :
■■
θjk, j+1, j+2 + θkj , j+1, j+2 + θkj+, j ,1j+2 + θkj+, j ,2j+1 = 2π ; and
■■
■■
interior vertex sum condition, ∀vi ∈ Vint : vk ∈qk θi = 2π ,
∑
November/December 2015
10/22/15 1:38 PM
vj vk
where the notations refer to Figure 1. Obviously, the first two constraints imply that q < 2p for arbitrary angle q. The boundary condition is important because it determines the shape of the parameterization mesh. In this article, we only consider the natural boundary conditions, thus the following constraint on the boundary vertex (see Figure 1) is needed:
vj+1 vj+m vj+2 vj+6
vi vj+3
∀vk ∈ Vbdy :
∑
k vk ∈qk , j , j +1, j +2 q j , j+1, j+2
< 2p vj+4
The angle system is therefore reduced an energy minimizing problem with 4|Q| variables, |Q| + |Vint| equality constraints and 4|Q| + |Vbdy| inequality constraints. As a quadric constraint optimization problem, this problem can be solved by many optimization software applications, such as MOSEK.
Figure 1. Boundary vertex vk and interior vertex vi on a quadrilateral mesh. The constraint is on the boundary vertex because we only consider the natural boundary conditions in this work. vj+1
vi
Free Boundary Planar Discrete Conformal Mapping
l r l r cot αij + cot αij + cot βij + cot βij u − u Ep = j i 4 4 4 4 [ vi v j ]∈E
∑
γij′ η ′ 2 + + cot ij ui − u j′ , cot 4 4 vi v j ′ ∈Ed
g6xie.indd 53
vi vj ′
β lij
ηij ′
β rij
vj
2
Figure 2. Edges, diagonals, and angles in the discrete conformal energy calculation. Both edges and diagonals are considered in our energy function, forming a new angle system.
where E and Ed represent the edge and diagonal sets, respectively. Figure 2 describes cot angles in the energy function. Because coefficients of the energy function are always positive, the local injectivity of the discrete conformal mapping can be guaranteed. Our energy function is similar to the one used in previous work,6 but we modified it in two ways: we consider the diagonals and use the angles calculated in the last section. An obvious advantage of using these angles is we can efficiently avoid edges overlapping each other. The quadric energy function only depends on the angles, so the conformality of the discrete mapping can be preserved better by using the angles we described in the last section instead of the angles from the original surface. To flatten a mesh satisfying minimal angle distortions, the minimum discrete conformal energy is required. Thus, we just need a vanishing derivative of Ep with respect to ui:
r ij
α lij
∑
γij ′ α
Conformal mapping is angle-preserving. To unfold a mesh into a plane with minimal angle distortions between the original mesh and the corresponding planar mesh, we devise a discrete conformal energy function that is motivated by the previous version defined on triangular meshes.1
vj+5
r l r αl ∂E p cot ij + cot αij + cot βij + cot βij u − u =2 j) ( i ∂ui 4 4 4 4 v j ∈dN(i)
∑
η ′ γij′ , +2 + cot ij (ui − u j′ ) cot 4 4 v j ′ ∈iN(i)
∑
(1)
where angles have the same meaning as before and dN(i) represents the points that directly connect with vi. Thus we also call dN(i) the direct onering neighborhood of vi. On the other hand, iN(i) denotes the points that share a common quadrilateral with vi but do not connect with vi. We call iN(i) the indirect one-ring neighborhood of vi. From Equation 1, we get alinear equations with a coefficient matrix that is sparse and coefficients that are related to angles we established earlier. To solve this equation, we still need to prescribe the boundary condition. The natural boundary condition is used when we calculate the angle system, so we map the boundary of a mesh into a plane by using a length-preserving free boundary manner. IEEE Computer Graphics and Applications
53
10/22/15 1:38 PM
Human-Centered Data Visualization Pi
Pj+1
Pj –1
Pi
Pj+1
Pi Pj ′
Pj ′ (a)
(b)
Pj (c)
βj1
Pj+1
α j4
αj2
Pj ′
Pj
Pj
α j1
Pj
Pi
Pi
βj2
Pj
α j3
Pj ′ (d)
(e)
Pj+1
Pj ′
Figure 3. Delaunay and anti-Delaunay triangulation. (a) One-ring neighbors of a point Pi. (b) Delaunay triangulation of one-ring neighbors of Pi on a quadrilateralmesh. (c) Anti-Delaunay triangulation of one-ring neighbors of Pi on a quadrilateral mesh. (d) Anti-Delaunay triangulation of a quadrilateral. (e) Delaunay triangulation of a quadrilateral.
Spherical Discrete Conformal Mapping Here we focus on the genus zero closed quadrilateral meshes, whose natural parameterization domain is the unit sphere in R3.
Laplacian of Quadrilateral Meshes In this subsection, we investigate the Laplacian on quadrilateral meshes by taking advantage of the calculation formulation on triangular meshes.7 Let M ∈ R3 be a two-manifold triangular mesh. Given an arbitrary vertex P i of M, its discrete Laplacian can be formulated as D M Pi = 2
∑w
ij
(Pi − Pj ) ,
j∈N(i)
where N(i) is the index set of the one-ring neighboring points of P i. Mathieu Desbrun and his colleagues used the COT scheme to calculate wij.7 To calculate the discrete differential operator, we first introduce the Delaunay triangulation and anti-Delaunay triangulation of the one-ring neighbors of a point P i on a quadrilateral mesh. The one-ring neighbors of a point P i on a quadrilateral mesh form a set of quadrilaterals that share the same vertex P i. Figure 3a provides an illustration. Delaunay triangulation means every triangle in a triangulation contains no interior points, and it maximizes the minimum angles. A single quadrilateral Delaunay triangulation satisfies the summation of opposite angles is less than 180 degrees (aj1 + aj4 < 180 degrees as in Figure 3d). For simplicity, we denote the Delaunay triangulation and anti-Delaunay triangulation of a one-ring neighbors of a point on a quadrilateral mesh as 1-DTQ and 1-aDTQ, respectively. Now we give the formal definitions of both. Definition 1 (1-DTQ). 1-DTQ refers to every one-ring neighboring quadrilateral of a point P that is a single quadrilateral Delaunay triangulation. 54
g6xie.indd 54
Definition 2 (1-aDTQ). 1-aDTQ refers to every onering neighboring quadrilateral of a point P that is not a single quadrilateral Delaunay triangulation. 1-DTQ and 1-aDTQ represent the good and bad triangulation. Furthermore, the Laplacian is estimated with the aid of 1-DTQ and 1-aDTQ. Figures 3b and 3c illustrate the 1-DTQ and 1-aDTQ, respectively. According to the COT scheme,7 we now derive the Laplacian of P i on quadrilateral meshes. In 1-DTQ and 1-aDTQ, some triangles (such as DP jP j′P j+1 in Figure 3c) are not used when computing the Laplacian. Therefore, vertices such as P j′ are only used once, while vertices such as P j are used D twice exactly. Let AM (Pi ) and AMa (Pi ) represent the area of 1-DTQ and 1-aDTQ, respectively. Then we give the area estimation Q AM (Pi ) =
1 D 1 a AM (Pi ) + AM (Pi ). 2 2
Consequently, the Laplacian of P i can be formulated as DQM (Pi ) =
(Pi )
1 Q 2 AM
∑w j
ij
(Pi − Pj ) − ∑ w ij′ (Pi − Pj′ ) , j′
(2) where wij = (cotaj2 + cotaj–1,3 + cotbj2 + cotbj–1,1)/2, wij′ = (cotaj1 + cota2)/2. Figure 3d provides descriptions of aik and bik.
Harmonic Energy of Quadrilateral Meshes Because conformal mapping is angle-preserving, it has extensive applications in computer graphics and medical image processing. 2 According to the harmonic analysis theory,8 a mapping is harmonic if it minimizes the harmonic energy. Moreover, a degree one harmonic map h: M → S2 is conformal. Suppose f is a function defined on a surface M. Then the harmonic energy of f is
November/December 2015
10/22/15 1:38 PM
Pk
Pi
E(f) = ∫M ||∆f||2 ds,
Pi
where ∇ is the gradient operator and ds is the area element. When M is a triangular mesh and f is a piecewise linear function defined on M, then E(f) has the following discretization:
∑
E( f ) =
2
g i1 = (1 – u)(1 – v)f i + (1 – u)vf j + u(1 – v)f j′ + uvf j+1 , g i2 = (1 – u)(1 – v)f j + (1 – u)vf j′ + u(1 – v)f j+1 + uvf i , g i3 = (1 – u)(1 – v)f j′ + (1 – u)vf j+1 + u(1 – v)f i + uvf j , and g i4 = (1 – u)(1 – v)f j+1 + (1 – u)vf i + u(1 – v)f j + uvf j′ , where u, v ∈ [0, 1]. According to the definition of the harmonic energy, we have the following energy of g i1 defined on Qi: 1 2 ( f i − f j′ ) 3 1 1 1 2 2 2 + ( f i − f j ) + ( f j − f j+1 ) + ( f j+1 − f j′ ) 3 3 3 1 1 + ( f i − f j′ )( f j − f j+1 ) − ( f i − f j )( f j+1 − f j′ ) . 3 3 1
∫ ∫ 0
1
0
2
∇g i1 du dv =
The E(g i2), E(g i3), and E(g i4) have similar representations. Let f=
1 ( g i1 + g i2 + g i3 + g i 4 ) . 4
Then the energy of f defined on Qi can be formulated as Ei ( f ) =
(
1 2 2 2 4 ( f i − f j′ ) + 4 ( f j − f j+1 ) + 2( f i − f j ) 12 2
2
2
+2( f j+1 − f j′ ) + 2( f i − f j′ ) + 2( f i − f j+1 )
)
∑ 6 E ( f ) = ∑ (2 ( f − f ) i
i
j′
2
2
2
)
+ ( f i − f j ) + ( f j+1 − f j′ ) + ( f i − f j′ ) + ( f i − f j+1 ) .
g6xie.indd 55
(b)
Figure 4. Two kinds of considered edges in the harmonic energy calculation. (a) The edge [PiPj] shared by two adjacent quadrilaterals. (b) The diagonal edge [PiPj′] dividing a quadrilateral into two quadrilaterals Ta, Tb. Their weights for harmonic energy calculation are the summation of the cotangents of related angles.
According to the Tuette energy, we also could define a Laplacian, and we call it the Tuette-Laplacian, which is defined as the derivative of ET(f) with respect to f. ∆T f (Pi ) =
∑
4 ( fi − f j ) +
[ Pj Pj ]∈E
∑
4 ( f i − f j′ ) ,(4)
Pj Pj ′ ∈Ed
where Ed is the set of diagonal edges. With reference to the relationship between the Tuette energy and the harmonic energy defined in a triangular mesh, we define the harmonic energy of quadrilateral meshes as follows:
∑
2
∑
2
EH ( f ) = kPi ,Pj f i − f j + kPi ,Pj ′ f i − f j′ . Pi Pj ′ ∈Ed [ PiPj ]∈E (5) Suppose [P iP j] has two adjacent quadrilaterals Qa = {P iP jP j′P j+1} and Qb = {P iP jPk′Pk} (see Figure 4a). Let ∠PiPj Pj′ ∠Pi Pj′ Pj+1 + cot 4 4 ∠Pj′ Pj+1Pi ∠Pj+1PiPj + cot + cot . 4 4
k aPiPj = cot
Similarly, kPbi ,Pj is defined. Thus, kPi ,Pj = kPai ,Pj + kPbi ,Pj . In addition, a diagonal edge [P iP j′] divides a quadrilateral into two triangles T a, T b (see Figure 4b), then we define ∠PiPj Pj′ ∠Pi Pi+1Pj′ + cot . 4 4
Spherical Discrete Conformal Mapping
2
+ 2( f j − f j+1 ) 2
(a)
kPi ,Pj ′ = cot
2
2
Pj ′
Pj ′
.
Therefore, we define the Tuette energy of f on M as follows: ET ( f ) =
Pj
Pj
where u and v are any two vertices of M, [u, v] is an edge of M, and ku,v is the coefficient that is related to the two adjacent triangles of [u, v]. The energy is known as the Tuette energy, especially when ku,v = 1. Suppose f is a function defined on a quadrilateral mesh M, and f i is the function value of vertex P i. For given quads Qi = {P iP jP j′P j+1} (see Figure 3d), we use the following formulations to approximate f on the quads Qi:
E ( g i1 ) =
Tα Qα
ku ,v f (u) − f ( v ) ,
[ u ,v ]∈M
Tβ
Pi+1
Qβ
Pk ′
Pi+1
(3)
For a genus zero closed mesh M1 and a unit spherical mesh S2, there are many conformal mappings between them. That is, the conformal mapping between two meshes is not unique. The basis of IEEE Computer Graphics and Applications
55
10/22/15 1:38 PM
Human-Centered Data Visualization
Input: A quadrilateral mesh M1, time step dt, threshold e Output: A spherical Tuette mapping h 1. Compute the Guassian map n: M1 → S2 . Let h = n, compute the Tuette energy ET (h). 2. Update the map: h = h – dt × DTh. 3. Compute the Tuette energy ET′ (h) . 4. If δE = ET′ (h) − ET (h) < ε , return h. Otherwise let ET (h) = ET′ (h) , and repeat steps 2–4.
component DN(h) and the tangential component DT(h)—that is, D(h) = DN(h) + DT(h). To ensure the unique solution, an extra constraint is required. The constraint is that the mass center of S2 is coincided with the origin:
∫
S2
h dσ M1 = 0 ,(9)
Figure 5. The algorithm flow chart for building Spherical Tuette mapping. The key of the algorithm is to compute the homeomorphism h by minimizing the Tuette energy in Equation 3.
where σ M1 is the area element. Then we can construct the geometric partial differential equation
Input: A quadrilateral mesh M1, time step dt, threshold e Output: A spherical conformal mapping h 1. Compute the spherical Tuette map ht: M1 → S2 . Let h = ht, compute the harmonic energy Eh = EH (h). 2. Update the map: h = h – dt × DTh. 3. Compute a Möbius transformation T, such that T ° h satisfies the constraint in Equation 9 of the mass center. 4. Compute the harmonic energy Eh′ = EH (h) . 5. If δE = Eh′ − E h < ε , return h. Otherwise let , and repeat steps 2–5.
with the constraint of Equation 9. The steady state solution h of the geometric flow is the conformal mapping from M1 to S2. Equation 10 can be efficiently solved by the steepest descendent algorithm. The basic idea of computing the discrete spherical conformal mapping has two steps: First, we compute a homeomorphism h from M1 to S2, and then we diffuse h to be harmonic. To compute the homeomorphism h between M1 and S2, we construct the spherical Tuette mapping as in Figure 5, which minimizes the Tuette energy defined in Equation 3. Viewing the spherical Tuette mapping as the initial homeomorphism, we construct the spherical conformal mapping for genus zero surfaces as in Figure 6.
Figure 6. The algorithm flow chart for building Spherical harmonic mapping. The spherical Tuette mapping is used as the initialization.
computing conformal mappings is that degree one harmonic mappings are conformal for genus zero surfaces.9 For a map h : M1 → S2 where h = (h0, h1, h2), the Tuette energy Er(h) of h can be defined as ET ( f ) =
2
∑E
T
(hi ) ,(6)
i=0
where ET(hi) is the same as in Equation 3. Also, the harmonic energy EH(h) of h is given by EH (h) =
2
∑E
H
(hi ) ,(7)
i=0
where EH(hi) is defined as in Equation 5. The Laplacian of h is D(h) = (D(h0), D(h1), D(h2)),(8) where D is the Laplacian operator in Equations 2 or 4. In practice, we use Equation 4 for spherical Tuette mapping (see Figure 5), whereas we use Equation 2 for spherical harmonic mapping (see Figure 6). Because a map h : M1 → S2 is harmonic if and only if D(h) has vanishing tangential component, h could be written as the sum of the normal 56
g6xie.indd 56
∂h = −∆T h (10) ∂t
Experimental Results In this section, we demonstrate the proposed quadrilateral mesh parameterization methods on three 3D models—Nicolo, Venus, and Moai—by displaying the texture mapping results based on specific parameterization manners. We also report the quantification evaluation based on the angle distortion. Some related parameterization methods were also carried out for comparison. Figure 7 illustrates the parameterization results on the Nicolo model when using different disk parameterization methods based on triangular meshes: least squares conformal Maps (LSCM),10 intrinsic,1 and low-stretch methods.11 To execute the triangular mesh based parameterization algorithms, we converted the quadrilateral Nicolo model (left image in Figure 7a) into the triangle mesh model (left image in Figure 7c) by connecting diagonals. As shown, our method can preserve the angles better than other methods. Furthermore, because our method adopts the free boundary condition, it is in a fair way to attain a very small change in the boundary length. In the experiments, the maximum boundary length ratios between the parameterized meshes and the original meshes are below 3 percent.
November/December 2015
10/22/15 1:38 PM
(a)
(b)
(c) Figure 7. The disk parameterization on the Nicolo model. (a) The quadrilateral Nicolo model and its texture mapping result achieved by our parameterization method. (b) The corresponding disk parameterization domain of the Nicolo model. (c) The triangulated Nicolo model and the parameterization results generated by least squares conformal maps (LSCM),10 intrinsic parameterization,1 and low-stretch parameterization methods,11 respectively, from left to right.
The spherical parameterization was implemented on the Venus and Moai models. Figure 8 shows the texture mapping results on Venus model based on different parameterization methods. For the sake of comparison, the figure shows the result on the quadrangular mesh when using the Maya software and our method and the result on the triangulated mesh when using the conformal spherical parameterization algorithm. 2 Our method attains better angle preservation
(a)
(b)
than the Maya software (especially in the neck region), and our result are comparable with the Gu parameterization method.2 The quadrangular Venus and Moai models were also parameterized onto unit spheres using several methods. As shown in Figures 8 and 9, our method attains slightly better result than the Maya software, but its results are worse than the Gu triangular mesh based parameterization method. 2
(c)
Figure 8. The spherical parameterization on the Venus model. (a) The Venus model and the parameterization result by our method. (b) The texture mapping results when using Maya software (left) and our method (right). (c) The triangulated Venus model. (d) The texture mapping based on the Gu parameterization method2 on the triangular mesh.
g6xie.indd 57
IEEE Computer Graphics and Applications
57
10/22/15 1:38 PM
Human-Centered Data Visualization
(a)
(b)
(c)
Figure 9. The spherical parameterization on the Moai model. (a) The Moai model. (b) The spherical parameterization domain of the Moai model (left), the texture mapping results based on the spherical parameterization by using Maya software (middle), and our method (right). (c) The triangulated Moai model. (d) The texture mapping result based on the Gu parameterization method on the triangular mesh.2
Table 1. Angle distortion of different parameterization methods. Model
Planar parameterization methods
Spherical parameterization methods
Low stretch
Intrinsic
LSCM
Ours
Gu
Nicolo
0.9312
0.3230
0.2995
0.0879
–
Venus
–
–
–
–
0.0025
0.02175
Moai
–
–
–
–
0.074
0.24
To quantify the parameterization distortion, we use the following function to measure the angle distortion: TA =
n
k
∑∑ f =1 i=1
σp 2 ρ f io − 1 ,(11) σi
where n is the face number and k = 3(4) when the mesh is triangular (quadrilateral). The σip and σio are the angle after parameterization and the original angle, respectively. The weight rf is defined as rf = Af/Sf A f, where Af is the area of the triangle or quadrilateral f. Table 1 gives the angle distortions of different models with respect to various methods, including the low-stretch,11 intrinsic parameterization,1 LSCM,10 Gu conformal spherical parameterization,2 and our planar and spherical parameterization. The table shows that our spherical parameterization method gets larger distortion values than the Gu method,2 but our planar parameterization method attains the smallest distortion value among all the compared planar parameterization algorithms. These examples all have demonstrated the feasibility and practicability of our parameterization methods for quadrilateral meshes. For the planar case, unlike existing parameterization methods that 58
g6xie.indd 58
Ours –
simultaneously determine the angles and boundaries,1,10,11 our method specifies the angle system in advance so that it has better conformal capability. For the spherical case, because the harmonic mapping on a quadrilateral mesh is not well-defined in mathematics, our method in fact attains slightly worse performance compared with the triangle mesh based parameterization methods.2
W
ith the growing use of quadrilateral meshes in 3D modeling, direct parameterization methods for quadrilateral meshes are required. In the future, we will investigate the angle-preserving parameterization globally for high genus quadrilateral meshes. Furthermore, we would like to extend the proposed scheme to the parameterization of a point cloud.
Acknowledgments We thank Joel Daniels II from the New York University Polytechnic Institute for providing quad meshes (Moai and Venus quad models). This work is supported by the National Nature Science Foundation of China (grants 61202223 and 61373003) and the GuangZhou Program (grant 201508010032).
November/December 2015
10/22/15 1:38 PM
References 1. M. Desbrun, M. Meyer, and P. Alliez, “Intrinsic Parameterizations of Surface Meshes,” Computer Graphics Forum, vol. 21, no. 3, 2002, pp. 209–218. 2. X. Gu et al., “Genus Zero Surface Conformal Mapping and Its Application to Brain Surface Mapping,” IEEE Trans. Medical Imaging, vol. 23, no. 8, 2004, pp. 949–958. 3. A. Khodakovsky, N. Litke, and P. Schröder, “Globally Smooth Parameterizations with Low Distortion,” ACM Trans. Graphics, vol. 22, no. 3, 2003, pp. 350–357. 4. A. Sheffer and E. de Sturler, “Parameterization of Faceted Surfaces for Meshing Using Angle Based Flattening,” Engineering with Computers, vol. 17, no. 3, 2001, pp. 326–337. 5. L. Kharevych, B. Springborn, and P. Schröder, “Discrete Conformal Mappings via Circle Patterns,” ACM Trans. Graphics, vol. 25, no. 2, 2006, pp. 412–438. 6. L. Liu, C.M. Zhang, and F. Cheng, “Parametrization of Quadrilateral Meshes,” Proc. 7th Int’l Conf. Computational Science (ICCS), Part II, LNCS 4488, Springer, 2007, pp. 17–24. 7. M. Desbrun et al., “Implicit Fairing of Irregular Meshes Using Diffusion and Curvature Flow,” Proc. 26th Ann. Conf. Computer Graphics and Interactive Techniques (SIGGRAPH), 1999, pp. 317–324. 8. R. Schoen and S.T. Yau, Lectures on Harmonic Maps, Harvard Univ. Press, 1997. 9. X. Gu and S.T. Yau, “Computing Conformal Structures of Surfaces,” Comm. Information and Systems, vol. 2, no. 2, 2002, pp. 121–146. 10. B. Levy et al., “Least Squares Conformal Maps for Automatic Texture Atlas Generation,” ACM Trans. Graphics, vol. 21, no. 3, 2002, pp. 362–371.
11. S. Yoshizawa, A. Belyaev, and H.P. Seidel, “A Fast and Simple Stretch-Minimizing Mesh Parameterization,” Proc. Shape Modeling and Application, 2004, pp. 200–208. Wenyong Gong is a senior algorithm engineer at AEE in Shenzhen, China. His research interests include computer graphics, image reconstruction, and image recognition. Gong has a PhD from the Institute of Mathematics at Jilin University. Contact him at [email protected]. Xiaohua Xie (corresponding author) is a research professor at Sun Yat-Sen University. His research interests include image processing, computer vision, pattern recognition, and computer graphics, especially focusing on image understanding and object modeling. Xie has a PhD in applied mathematics from Sun Yat-Sen University. He received the Overseas High-Caliber Personnel Award in Shenzhen in 2012. Contact him at [email protected]. Rui Ma is a PhD candidate in the School of Computing Science at Simon Fraser University. His research interests include computer graphics, high-level geometry processing, and shape analysis. Ma has an MS in mathematics from Jilin University. Contact him at [email protected]. Tieru Wu is a professor in the Institute of Mathematics at Jilin University, China. He research interests include developing various techniques of computer graphics, geometry processing, machine learning in image processing, and image retrieval. His research is supported in part by the National Science Foundation of Chinese. Wu has a PhD in quantitative economics from Jilin University. Contact him at [email protected].
EXECUTIVE STAFF PURPOSE: The IEEE Computer Society is the world’s largest association of computing professionals and is the leading provider of technical information in the field. MEMBERSHIP: Members receive the monthly magazine Computer, discounts, and opportunities to serve (all activities are led by volunteer members). Membership is open to all IEEE members, affiliate society members, and others interested in the computer field. COMPUTER SOCIETY WEBSITE: www.computer.org
Next Board Meeting: 15–16 November 2015, New Brunswick, NJ, USA EXECUTIVE COMMITTEE President: Thomas M. Conte President-Elect: Roger U. Fujii; Past President: Dejan S. Milojicic; Secretary: Cecilia Metra; Treasurer, 2nd VP: David S. Ebert; 1st VP, Member & Geographic Activities: Elizabeth L. Burd; VP, Publications: Jean-Luc Gaudiot; VP, Professional & Educational Activities: Charlene (Chuck) Walrad; VP, Standards Activities: Don Wright; VP, Technical & Conference Activities: Phillip A. Laplante; 2015–2016 IEEE Director & Delegate Division VIII: John W. Walz; 2014–2015 IEEE Director & Delegate Division V: Susan K. (Kathy) Land; 2015 IEEE Director-Elect & Delegate Division V: Harold Javid
BOARD OF GOVERNORS Term Expiring 2015: Ann DeMarle, Cecilia Metra, Nita Patel, Diomidis Spinellis, Phillip A. Laplante, Jean-Luc Gaudiot, Stefano Zanero Term Expriring 2016: David A. Bader, Pierre Bourque, Dennis J. Frailey, Jill I. Gostin, Atsuhiro Goto, Rob Reilly, Christina M. Schober Term Expiring 2017: David Lomet, Ming C. Lin, Gregory T. Byrd, Alfredo Benso, Forrest Shull, Fabrizio Lombardi, Hausi A. Muller
g6xie.indd 59
Executive Director: Angela R. Burgess; Director, Governance & Associate Executive Director: Anne Marie Kelly; Director, Finance & Accounting: Sunny Hwang; Director, Information Technology Services: Ray Kahn; Director, Membership: Eric Berkowitz; Director, Products & Services: Evan M. Butterfield; Director, Sales & Marketing: Chris Jensen
COMPUTER SOCIETY OFFICES Washington, D.C.: 2001 L St., Ste. 700, Washington, D.C. 20036-4928 Phone: +1 202 371 0101 • Fax: +1 202 728 9614 • Email: [email protected] Los Alamitos: 10662 Los Vaqueros Circle, Los Alamitos, CA 90720 Phone: +1 714 821 8380 • Email: [email protected] MEMBERSHIP & PUBLICATION ORDERS Phone: +1 800 272 6657 • Fax: +1 714 821 4641 • Email: [email protected] Asia/Pacific: Watanabe Building, 1-4-2 Minami-Aoyama, Minato-ku, Tokyo 107-0062, Japan • Phone: +81 3 3408 3118 • Fax: +81 3 3408 3553 • Email: [email protected]
IEEE BOARD OF DIRECTORS President & CEO: Howard E. Michel; President-Elect: Barry L. Shoop; Past President: J. Roberto de Marca; Director & Secretary: Parviz Famouri; Director & Treasurer: Jerry Hudgins; Director & President, IEEE-USA: James A. Jefferies; Director & President, Standards Association: Bruce P. Kraemer; Director & VP, Educational Activities: Saurabh Sinha; Director & VP, Membership and Geographic Activities: Wai-Choong Wong; Director & VP, Publication Services and Products: Sheila Hemami; Director & VP, Technical Activities: Vincenzo Piuri; Director & Delegate Division V: Susan K. (Kathy) Land; Director & Delegate Division VIII: John W. Walz revised 5 June 2015
IEEE Computer Graphics and Applications
59
10/22/15 1:38 PM
