G2  B-spline interpolation to a closed mesh

Computer-Aided Design 43 (2011) 145–160

Contents lists available at ScienceDirect

Computer-Aided Design journal homepage: www.elsevier.com/locate/cad

G2 B-spline interpolation to a closed mesh Kan-Le Shi a,b,c,d,∗ , Sen Zhang a,b,c,d , Hui Zhang a,c,d , Jun-Hai Yong a,c,d , Jia-Guang Sun a,c,d , Jean-Claude Paul a,c,d,e a

School of Software, Tsinghua University, Beijing 100084, PR China

b

Department of Computer Science and Technology, Tsinghua University, Beijing 100084, PR China

c

Key Laboratory for Information System Security, Ministry of Education of China, Beijing 100084, PR China

d

Tsinghua National Laboratory for Information Science and Technology, Beijing 100084, PR China

e

INRIA, Nancy 54600, France

article

info

Article history: Received 2 February 2010 Accepted 20 October 2010 Keywords: Smooth Mesh interpolation Coons B-spline surface G2 -continuity

abstract This paper focuses on interpolating vertices and normal vectors of a closed quad-dominant mesh1 G2 continuously using regular Coons B-spline surfaces, which are popular in industrial CAD/CAM systems. We first decompose all non-quadrangular facets into quadrilaterals. The tangential and second-order derivative vectors are then estimated on each vertex of the quads. A least-square adjustment algorithm based on the homogeneous form of G2 continuity condition is applied to achieve curvature continuity. Afterwards, the boundary curves, the first- and the second-order cross-boundary derivative curves are constructed fulfilling G2 continuity and compatibility conditions. Coons B-spline patches are finally generated using these curves as boundary conditions. In this paper, the upper bound of the rank of G2 continuity condition matrices is also strictly proved to be 2n − 3, and the method of tangentvector estimation is improved to avoid petal-shaped patches in interpolating solids of revolution. Several examples demonstrate its feasibility. © 2010 Elsevier Ltd. All rights reserved.

1. Introduction There are many practical design situations in which smooth parametric surfaces are constructed by interpolating meshes that control the global geometry of surfaces. Since mesh has become a popular representation of 3D models in many fields, such as geometric modeling, reverse engineering, animation and mechanical analysis, the hybrid modeling of mesh and spline has attracted special attention from academia and industry in recent years. Converting a mesh into parametric surfaces is an essential operation. It has been widely discussed. The key problem is to construct geometric continuous surface patches interpolating boundary conditions derived from the original mesh. G2 continuity is usually required in many engineering or scientific problems because of some background in arts and physics, such as car-body and ship-hull designing. As B-spline is a standard surface form in most CAD/CAM or geometric modeling systems, interpolating mesh using B-spline surfaces has its distinct importance.

∗ Corresponding author at: School of Software, Tsinghua University, Beijing 100084, PR China. E-mail addresses: [email protected], [email protected] (K.-L. Shi). 1 A mesh is quad-dominant if it contains quad facets as its majority. The quantity of quad facets is much more than the quantity of triangular and other multi-sided facets. 0010-4485/$ – see front matter © 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.cad.2010.10.004

Some notable related work has been done in the past decades. For G1 and G2 continuity2 between two adjacent parametric surfaces, criterions and construction methods have been proposed [1–5]. Peters [6] analyzed the C k continuity problem of n patches around a common vertex. Jones [7] introduced the G2 continuity conditions on the common vertex as a linear system, and presented a conjecture that the rank of the condition matrix is 2n − 3. In order to smooth meshes of arbitrary topology, different types of subdivision schemes were proposed. The original mesh vertices are considered as the control points of these subdivision surfaces. Popular subdivision algorithms like Catmull–Clark [8] can achieve C 2 almost everywhere except for the so-called extraordinary vertices. A lot of efforts had been made to solve this problem. Finally the global G2 -continuity had been achieved by Loop, Peters, etc., [9–12]. Methods for interpolating vertices of a prescribed mesh using subdivision surfaces were also proposed [13]. Since these methods face the difficulties in data exchange, accurate evaluation and geometric

2 A surface or adjacent surfaces are Gn continuous on a specified point if there is a parameterization, under which the nth-order parametric continuity is satisfied. The geometric meaning of G1 continuity, or the tangential/normal continuity, is that the normal direction is unique on specified point. Also, for the G2 continuity, or the curvature continuity, it requires a single fixed directional curvature in an arbitrary direction.

146

K.-L. Shi et al. / Computer-Aided Design 43 (2011) 145–160

processing (such as computing intersection curves), they are rarely applied in industrial modeling systems. Peters [14], Bajaj and Ihm [15] proposed new splines and implicit algebraic surfaces to smooth meshes besides subdivision. These methods present good results but suffer from the compatibility problem in data exchange with classical surface modeling techniques. Since the mesh interpolation problem is similar to interpolating curve networks, which was studied in depth for ship-hull design, there are many worthy papers that can be used for reference. Peters [16,17] handled the vertex enclosure constraints and gave an interpolation method of a mesh of curves using G1 regular parametric surfaces with one polynomial piece per facet. A G1 Bézier solution over irregular curve networks was given by Cho et al. [18] several years later. Ye [19] achieved G2 continuity in constructing surfaces by interpolating given quadrilateral GC 2 quintic curve meshes. His method handles G2 compatibility problems, and generates (9, 9)-degree (meaning in both u- and v direction, the degrees are 9) Bézier surfaces with prescribed quintic curve network. But it only handles the rectangular curvature continuous quintic curve mesh, and does not discuss the continuity problem on a vertex shared by n (n ̸= 4) surfaces (expressed as the complicated-topology problem). Hermann [20] provided an interpolation method using G2 biquintic Gregory patches, which are not easily imported into popular modeling systems because of the surface representation. On local interpolation of point meshes, Li and Liu [21] introduced a theoretical analysis, and gave a minimum-degree declaration of this problem. In their article, G2 conditions and criterions of derivative vectors are derived but it is short of concrete construction process for parametric surfaces. For a G1 smooth join of tensor product parametric surfaces, both simple and complicated topological types were considered [22–24]. The most related papers mainly focus on triangular meshes. Hahmann and Bonneau [25,26] interpolated triangular meshes using G1 continuous triangular parametric surfaces. Recently, Tong and Kim [27] applied the technique of G1 smoothing triangular meshes to approximating implicit surfaces with highquality, and they achieved good results. But a notable shortcoming is that the triangular Bézier surface used by their method is not well supported by most industrial CAD/CAM or geometric modeling systems. To solve this problem, we focus on quad-dominant rough meshes and use regular B-spline form. In this paper, we smooth a closed mesh by G2 -continuous (9, 9)-degree quadrilateral Coons B-spline surfaces with one piece per facet, interpolating all vertices and corresponding normal vectors. The input mesh needs to be closed and quad-dominant, with arbitrary topology. The normal vectors of the vertices can be auto-computed or interactively prescribed. We first decompose each non-quadrangular facet into co-vertex quads by estimating a central vertex and splitting boundary edges. T-shaped edges are then manipulated to avoid singularity problem in the following steps. For each vertex, the first- and the second-order derivative vectors are estimated and then least-squarely adjusted to satisfy G2 continuity conditions. After that, the boundary curves, their first- and second-order crossboundary derivatives of the Coons patches are constructed using these vectors, fulfilling G2 compatibility conditions. Finally, we generate the B-spline surfaces with these curves as boundary conditions. The main contributions of our work can be summarized as follows:

• Based on the reformulation of geometric smoothness constraints, the upper bound of the rank of the condition matrix of G2 continuity and compatibility on vertices is strictly proved to be 2n − 3. The construction method of the derived homogeneous form of G2 continuity conditions is also proposed.

• A complete surface interpolation method to a closed quaddominant mesh is proposed. It achieves global G2 continuity. The constructed parametric surfaces are all in regular Coons B-spline form, which benefits data exchange and further manipulation. All the temporary curves are also in B-spline form. The method of tangent-vector estimation is optimized for solids of revolution so that better shape can be achieved. The rest of this paper is organized as follows. Section 2 introduces the problem definition, the terms, and the main framework of the construction method. Section 3 presents the compact expression of G2 continuity conditions on vertices. Section 4 focuses on how to compute the tangent vectors, the normal curvatures and the mixed partial derivatives on vertices, which are named structural derivatives, as they mainly decide the shape of the patches. To handle the non-quadrangular facets, Section 5 presents a simple method to decompose all k-sided facets (where k ̸= 4), and also gives a method to manipulate Tshaped edges. Then we propose a method to solve the high-order compatibility problems, and construct the spanning curves, their first- and second-order cross-boundary derivatives in Section 6. Several examples in Section 7 underline the curvature continuity of the generated patches. Finally, the paper concludes in Section 8. 2. Algorithm framework The input of this problem is a quad-dominant mesh, with nv vertices {vi |i = 0, 1, . . . , nv − 1} and nf facets {fj |j = 0, 1, . . . , nf − 1}. The normal vectors {ni |i = 0, 1, . . . , nv − 1} can be prescribed or automatically computed. After the quadrangulation step, each k-sided facet (k ̸= 4) is decomposed into k quads (see Section 5). Our algorithm generates G2 B-spline patches with one piece per quadrangular facet. The four corner points of a constructed surface are the vertices of the corresponding quadrilateral exactly. The normal vectors of these corners points also interpolate the prescribed or pre-computed ones of the original vertices. The interpolation algorithm has the following steps. 1. Quadrangulate all non-quadrangular facets Each non-quadrangular facet is decomposed into quadrangular facets by connecting an estimated central vertex to middle points of the spanning curves of edges. The central vertex estimating algorithm utilizes the mean tangential direction of each edge such that sub-facets approximately interpolate the boundary normal vectors, where better shape can be achieved. 2. Manipulate T-shaped edges T-shaped edges, which lead to singularity problems while computing inverse matrices in the following steps, may exist in the input mesh as well as being generated in the previous step. We split the corresponding facet into two sub-facets to handle it. Section 5 gives the detailed algorithm about these two steps. After quadrangulation and removing T-shaped edges, the mesh only contains quadrangular facets. So the following steps only consider quadrangular meshes. 3. Compute structural derivatives and construct spanning curves First, we estimate tangent vectors and second-order derivative vectors for each vertex. A least-square adjustment algorithm is then applied to them to ensure that they can fulfill G2 continuity conditions on vertices. Since these vectors determine the shape of the final surfaces, we call them the structural derivatives in this paper (see the orange vectors in Fig. 1(b)). Furthermore, the algorithm is optimized for solids of revolution to achieve better shape, as they are usually represented by quad-dominant meshes. Using the derivative vectors as boundary conditions, we can generate the spanning curves, the boundary of the final patches (see Fig. 1(b)).

K.-L. Shi et al. / Computer-Aided Design 43 (2011) 145–160

147

Fig. 1. Smoothing a cube: (a) is the original cube, which contains 6 quads; (b) shows the structural derivatives and the spanning curves; (c) shows the first-order crossboundary derivative curves of one spanning curve; (d) presents the control points of each B-spline surface; (e) shows the constructed B-spline patches and their isoparametric curves; (f) is the rendered image of these surfaces; (g) shows their reflection lines, which underlines their G2 continuity; (h) is the Gaussian curvature graph.

4. Handle high-order compatibility problems and construct crossboundary derivatives The original vertices are also the common corners of adjacent parametric surfaces. It is necessary to handle the compatibility problem on each vertex. In our algorithm, we assume that two adjacent parametric surfaces G2 contact a common imaginary surface SI both (see Section 6). The boundary conditions of them are constructed by applying two different G2 -preserving transforms to SI , such that their G2 continuity can be preserved. In this paper, we derived the tangent, the twist, the second-order tangent, the secondorder partial/mixed twist, and the second-order total twist compatibility constraints. The boundary conditions of the common curves, the common cross-boundary derivatives and the common second-order cross-boundary derivatives of each pair of adjacent patches are then computed by solving these compatibility conditions. Then, the first- and the second-order cross-boundary derivative curves of the final surfaces are constructed in B-spline form (see Fig. 1(c)). 5. Construct Coons B-spline surfaces The previous two steps prepare the 0th- (the spanning curves), the first- and the second-order cross-boundary derivatives of the four boundaries of each parametric patch. So the second-order Coons construction method can be directly applied to construct the final patches using the boundary condition curves (see Fig. 1(d)–(f)). In summary, the main idea is to construct boundary curves and cross-boundary derivatives for each Coons B-spline surface using geometric and topology information of the prescribed mesh, satisfying global G2 continuity conditions (see Fig. 1(f)–(h)). We prove the continuity and compatibility properties theoretically. The following sections present detailed algorithms, derivations and proofs on each step.

Fig. 2. The normal vertex N and tangent vectors Vi on a vertex.

3. G 2 continuity conditions on vertices This section discusses the G2 continuity conditions on a n-sided vertex, where the normal vector N and all tangent vectors Vi (i = 0, 1, . . . , n − 1) are fixed. Shown in Fig. 2, we define

  ∂ Si (u, v) = Vi , u=0 ∂u v=0   ∂ 2 Si (u, v) = Cii , 2 u=0 ∂u v=0   ∂ 2 Si (u, v) = Cij , u=0 ∂ u∂v

  ∂ Si (u, v) = Vj , u=0 ∂v v=0   ∂ 2 Si (u, v) = Cjj , 2 u=0 ∂v v=0

(1)

v=0

where Si (i = 0, 1, . . . , n − 1) be parametric surfaces incident to O, and j = (i + 1) mod n. Suppose two parametric surfaces Si and Sp share the common vertex O (see Fig. 3). We can expand the total differential form. That is, d dt d2 dt 2

∂ du ∂ dv + ∂ u dt ∂v dt  2 2  2 2 du ∂ du dv ∂ 2 dv ∂ = + 2 + dt ∂ u2 dt dt ∂ u∂v dt ∂v 2 d2 u ∂ d2 v ∂ + 2 . + 2 dt ∂ u dt ∂v

=

148

K.-L. Shi et al. / Computer-Aided Design 43 (2011) 145–160

For dSi /dt, on O, we denote α = du/dt and β = dv/dt. Analogously, for dSp /dt, γ = −du/dt and δ = −dv/dt. According to the definition of G1 continuity [2,3], combining with Eq. (1), for arbitrary α and β , there exist γ and δ satisfying dSi /dt = α Vi + β Vj = dSp /dt = −γ Vp − δ Vq .

(2)

It follows that Vi , Vj , Vp and Vq are coplanar. So the G continuity condition is Vi · N = 0 (i = 0, 1, . . . , n − 1). Then, in the expression of d2 Si /dt 2 , on O, we denote αˆ = d2 u/dt 2 and βˆ = d2 v/dt 2 . Analogously, for d2 Sp /dt 2 , we denote γˆ = d2 u/dt 2 and δˆ = d2 v/dt 2 . According to the definition of G2 continuity, for arbitrary αˆ and βˆ , there exist γˆ and δˆ satisfying 1

 2 d Si /dt 2 = α 2 Cii + 2αβ Cij + β 2 Cjj + αˆ Vi + βˆ Vj d2 S /dt 2 = γ 2 Cpp + 2γ δ Cpq + δ 2 Cqq + γˆ Vp + δˆ Vq  2 p 2 d Si /dt = d2 Sp /dt 2 .

3.1. Reflexivity

Note that Vw · N = 0 (w ∈ {i, j, p, q}). Applying ·N to both sides of the equations above, we have

α 2 Cˆ ii + 2αβ Cˆ ij + β 2 Cˆ jj = γ 2 Cˆ pp + 2γ δ Cˆ pq + δ 2 Cˆ qq ,

(3)

where Cˆ w = Cw · N (w ∈ {ii, ij, jj, pp, pq, qq}). The geometric meaning of this condition is that, in an arbitrary direction Vij +

⃗ (Vij = α Vi + β Vj , Vpq = −γ Vp − δ Vq , see Fig. 3), the two Vpq = 0 corresponding curves Si (u(t ), v(t )), Sp (u(t ), v(t )) share the same normal curvature on the common vertex O. According to Eq. (2), we have [Vi , Vj , N ][α, β, 0]T + [Vp , Vq , N ][γ , δ, 0]T = 0⃗. Thus,

[γ , δ, 0]T = −[Vp , Vq , N ]−1 [Vi , Vj , N ][α, β, 0]T . Vp , Vq and N are not coplanar. It follows that its inverse matrix always exists. We denote mα→γ mα→δ

mβ→γ mβ→δ

 −[Vp , Vq , N ] [Vi , Vj , N ] = −1

.

.

 . . , .

[

mα→γ mα→δ

The reflexivity property shows that a surface is G2 continuous with itself on the corner point. That is Theorem 1. Mij→ij = I. We only need to choose p = i and q = j. Then, Rij→ij = I. Substitute the values for mα→γ , mα→δ mβ→γ and mβ→δ in Mij→ij , Theorem 1 can be directly proved. 3.2. Symmetry The symmetry property of the compact form indicates that if surface Si is G2 continuous with surface Sp , Sp is also G2 continuous with Si . That is Theorem 2. [Cˆ ii , Cˆ ij , Cˆ jj ]T = Mij→pq [Cˆ pp , Cˆ pq , Cˆ qq ]T if and only if [Cˆ pp , Cˆ pq , Cˆ qq ]T = Mpq→ij [Cˆ ii , Cˆ ij , Cˆ jj ]T . It means that the two equations are equivalent. One is redundant if the other one is used in the set of G2 conditions.

where

[γ , δ]T = Rij→pq [α, β]T =

Fig. 3. G2 continuity of two corresponding curves on Si and Sp , with two negative directions Vij and Vpq .

Proof. The equivalent expression of Theorem 2 is

]

mβ→γ [α, β]T . mβ→δ

1 Mij−→ pq = Mpq→ij .

Rij→pq is the first-order transition matrix, which represents the coordinate stretching of the linear reparameterization. Combining with (3), we have

Assume

α 2 Kα2 + 2αβ Kαβ + β 2 Kβ 2 ≡ 0, for arbitrary α and β , where Kα 2 = m2α→γ Cˆ pp + 2mα→γ mα→δ Cˆ pq + m2α→δ Cˆ qq − Cˆ ii ,

It follows that

(4)

mγ →α Rpq→ij = mγ →β

Kβ 2 = m2β→γ Cˆ pp + 2mβ→γ mβ→δ Cˆ pq + m2β→δ Cˆ qq − Cˆ jj . Consequently, the sufficient and necessary condition of Eq. (4) is Kα 2 = Kαβ = Kβ 2 = 0. That is,

[Cˆ ii , Cˆ ij , Cˆ jj ]T = Mij→pq [Cˆ pp , Cˆ pq , Cˆ qq ]T ,

(5)

where Mij→pq denotes the second-order transition matrix, having the following expression m2α→γ mα→γ mβ→γ m2β→γ

2mα→γ mα→δ mα→γ mβ→δ + mβ→γ mα→δ 2mβ→γ mβ→δ

mδ→α mδ→β

[α, β]T = Rpq→ij [γ , δ]T . ]

[

mα→γ = mα→δ

mβ→γ mβ→δ

]−1

.

1 Thus, both Mij−→ pq and Mpq→ij can be represented by mα→γ , mβ→γ , mα→δ and mβ→δ . Combining with Eq. (6), we find that the final expression is identical after simplification. 

+ mβ→γ mα→δ Cˆ pq + mα→δ mβ→δ Cˆ qq − Cˆ ij ,



[γ , δ]T = Rij→pq [α, β]T , [

Kαβ = mα→γ mβ→γ Cˆ pp + mα→γ mβ→δ Cˆ pq

(6)

m2α→δ mα→δ mβ→δ  . m2β→δ



The equation above is called the compact form of G2 continuity on vertex in this paper.

3.3. Transitivity The geometric meaning of the transitivity property is that, on the common vertex O, if Si and Sp are G2 continuous, and if Sp and Sr are G2 continuous, then Si and Sr are also G2 continuous. That is, Theorem 3. [Cˆ ii , Cˆ ij , Cˆ jj ]T = Mij→rs [Cˆ rr , Cˆ rs , Cˆ ss ]T if [Cˆ ii , Cˆ ij , Cˆ jj ]T = Mij→pq [Cˆ pp , Cˆ pq , Cˆ qq ]T and [Cˆ pp , Cˆ pq , Cˆ qq ]T = Mpq→rs [Cˆ rr , Cˆ rs , Cˆ ss ]T . Proof. The equivalent expression of Theorem 3 is Mij→rs = Mij→pq Mpq→rs .

(7)

K.-L. Shi et al. / Computer-Aided Design 43 (2011) 145–160

149

The overall compact form of G2 continuity conditions on an nsided vertex is

⃗. M [Cˆ 00 , Cˆ 01 , Cˆ 11 , . . . , Cˆ n−1,n−1, Cˆ n−1,0 ]T = 0 It is important to study the rank of M, which leads us to choose the minimum number of independent conditions to form a simplest M. Preserving connectivity of the closure graph, we choose edges (Si , S0 ), i = 1, 2, . . . , n − 1 (see the bold edges in Fig. 4). The condition matrix of edge (S0 , S0 ) is an identity. That is M01→01 = I3 (see the reflexivity property). For the conditions (Si , S0 ) and (Sj , S0 ) (j = (i + 1) mod n), their compact expressions are Fig. 4. Continuity closure of 5-sided corner vertex, where the bold edges show one of a minimum set of complete constraints.

Analogous to the definition of Vij and Vpq , we define Vrs = ξ Vr + ψ Vs . Since [γ , δ]T = Rij→pq [α, β]T and [ξ , ψ]T = Rpq→rs [γ , δ]T , we have [ξ , ψ]T = Rpq→rs Rij→pq [α, β]T . That is,

[

mα→ξ Rij→rs = Rpq→rs Rij→pq = mα→ψ

]

mβ→ξ . mβ→ψ

Thus, mα→ξ , mβ→ξ , mα→ψ and mβ→ψ can be represented by other coefficients without ξ and ψ . Substitute them in the left and the right parts of Eq. (7). After simplification, we have an identical equation. 

[Cˆ ii , Cˆ ij , Cˆ jj ]T = Mij→01 [Cˆ 00 , Cˆ 01 , Cˆ 11 ]T ,

(8)

[Cˆ jj , Cˆ jk , Cˆ kk ]T = Mjk→01 [Cˆ 00 , Cˆ 01 , Cˆ 11 ]T ,

(9)

where k = (j + 1) mod n. Note that the last scalar equation of Condition (8) and the first one of Condition (9) are the same. We can remove one of them for i = 0, 1, . . . , n − 1. The left 3(n − 1) − n = 2n − 3 equations are the final G2 continuity conditions. M is a (2n − 3) × 2n matrix. So, we have the following theorem.3 Theorem 5. The upper bound of the rank of M is 2n − 3, where M is a G2 condition matrix of an n-sided corner.

3.4. Nonsingularity

It follows that the remaining rows of M are independent for nonsingular occasions, and that the final M leads the simplest form of G2 constraints on a common vertex.

Mij→pq is never singular. First, we introduce a theorem about the determinant of Mij→pq . That is

3.6. Least-square construction of second-order derivative vectors using a homogeneous form of G2 conditions

Theorem 4. |Mij→pq | = |Rij→pq |3 .

This subsection focuses on the problem: how to construct the second-order derivative vectors to satisfy G2 continuity conditions by least-squarely moving the specified (pre-computed or estimated) corresponding vectors. We denote

This theorem can be proved by expanding the left and the right parts of the equation using variables mα→γ , mβ→γ , mα→δ and mβ→δ . Since matrix Rij→pq is a coordinate-transition matrix, it is singular if and only if Vi and Vj are linearly dependent, or Vp and Vq are linearly dependent. In application, all the neighboring tangent vectors point to different directions. That ensures the nonsingularity of Rij→pq . That is, |Rij→pq | ̸= 0. It follows that |Mij→pq | ̸= 0, which means that Mij→pq is nonsingular. 3.5. Closure of continuity constraints The compact form (see Eq. (5)) ensures G2 continuity of surfaces Si and Sp on the common vertex. Consequently, conditions from all surface pairs form a closure of continuity constraints. Shown in Fig. 4, each edge (Sa , Sb ) denotes the relation of G2 continuity of a surface pair and the corresponding compact form of continuity conditions. If the diagram is a complete graph, which means that arbitrary two surfaces are G2 continuous on the common vertex, it follows that all the surfaces are G2 continuous. The reflexivity, symmetry and transitivity properties indicate that it is necessary, in minimum, to choose (n − 1) edges to form a complete graph, as all the other edges can be derived using these properties. It follows that we only need to choose (n − 1) groups of conditions to sufficiently represent all the G2 continuity constraints. In fact, different cases require various edge set, especially for different parities. Here, combining with the construction method in the following sections, it is required that the adjacent surfaces are G2 continuous along the boundary. So adjacent surfaces also need to be G2 continuous on end points. In this application, our closure graph is complete, because the continuity conditions of (Sk , S(k+1) mod n ), k = 0, 1, . . . , n − 1, are all sufficient and necessary.

⃗x = [Cˆ 00 , Cˆ 01 , Cˆ 11 , . . . , Cˆ n−1,n−1, Cˆ n−1,0 ]T . ⃗ = [∆0 , ∆1 , . . . , ∆2n−1 ]T , fulfilling M · (⃗x + The target is to find ∆ ⃗ where ‖∆ ⃗ ) = 0, ⃗ ‖2 is the least. ∆ The simplest M derived in the previous subsection is not the best choice, because we removed redundant conditions asymmetrically, and this may cause the imbalance of movement distribution and the instability of numerical calculation. The most complete and well-proportioned method is to list all constraints of of two surfaces, and that leads to a  3 possible combinations  n(n − 1) × 2n condition matrix MF . But MF is not square and 2 there are too many redundant rows. Thus, we introduce another symmetrical MH , called the homogeneous form of condition matrix. In the closure graph, we choose the following n edges (instead of n − 1): (S0 , S1 ), (S1 , S2 ), . . . , (Sn − 1, S0 ). For each edge (Si , Sj ) (j = (i + 1) mod n), the condition is [Cˆ ii , Cˆ ij , Cˆ jj ]T = Mij→jk [Cˆ jj , Cˆ jk , Cˆ kk ]T , where k = (j + 1) mod n. Note that the last row of Mij→jk must be

[1, 0, 0], since Cˆ jj ≡ 1 × Cˆ jj + 0 × Cˆ jk + 0 × Cˆ kk . Therefore, these n conditions can be removed, with one row per group symmetrically. The left 2n rows form a new (2n × 2n) square condition matrix, and that is the homogeneous form MH . All surfaces play symmetric roles in this form as they contribute the same number of variables and equations. If we shift the subscripts of the surfaces, the result

3 Jones’ numerical experiments [7] also showed that the rank of G2 continuity conditions is 2n − 3 exactly for nonsingular cases.

150

K.-L. Shi et al. / Computer-Aided Design 43 (2011) 145–160

Fig. 6. Computing a tangent vector using a vector-projection method. Fig. 5. Structural vectors around facet f0 : v0−3 are the original vertices, N0−3 are their normal vectors, e0,0−3 are the mesh edges around f0 , Vi,j are the tangent vectors of vertices along their corresponding edge-directions.

⃗ ⃗ , fulfilling MH · (⃗x + ∆ ⃗ ) = 0, remains. So we solve a least-square ∆ where rankMH ≤ 2n − 3 (discussed in the previous subsection). ⃗ , satisfying The equivalent problem is to find a least-square ∆ ⃗ Compute first the general solution of ⃗ = −MH · ⃗x = b. MH · ∆ ⃗ by Gaussian elimination and Schimidt orthogonalization. That ∆ ⃗ = ∆ ⃗ 0 + [∆] · ⃗k, where ∆ ⃗ 0 is a particular solution, [∆] = is, ∆ ⃗ 1, ∆ ⃗ 2, . . . , ∆ ⃗ m ] are the coordinate axes of the solution space [∆ with dimension m, and ⃗ k = [k1 , k2 , . . . , km ]T is the coordinate vector. Then, it is equivalent to finding a proper ⃗ k, fulfilling ⃗ 0 , when that [∆] · ⃗ k is the least-square approximation to −∆ ⃗ The Moore–Penrose ⃗ is the least-square approximation to 0. ∆ pseudoinverse [28] just solves this problem. We first represent [∆] = U Σ V T by singular value decomposition. Then, set ⃗k = ⃗ 0 , where Σ is a diagonal matrix, and Σ + is constructed −V Σ + U T · ∆ by taking the reciprocal of each non-zero element of Σ . We can get ⃗ by substituting ⃗k in the general solution. the final result of ∆ ⃗ Since all derivatives are prescribed as vectors, we only apply ∆ along the common normal direction N to the original ones. That is, Cij ⇐ Cij + N · ∆i+j ,

i = 0, 1, . . . , n − 1; j = i, i + 1.

Thus, after the least-square movement, all adjusted derivative vectors satisfy G2 continuity conditions on the common vertex. In application, we use MH as the condition matrix and apply this method to compute G2 compatible second-order derivative vectors. 4. Structural derivatives Our algorithm generates one Coons B-spline surface per facet (see Fig. 5). Firstly, for each vertex vi , we compute its normal vector Ni , tangent vectors Vi,j (j = 0, 1, . . . , ne,i − 1, where ne,i is the number of in-directional half-edges, or the degree of vi ), and its second-order derivative vectors. These vectors are called the structural derivatives in this paper, since they control the shape of the patches to construct. This section presents the details on how to compute these derivative vectors. For a prescribed vertex v , assume its neighboring vertices are [v0 , v1 , . . . , vnn −1 ] anti-clockwise, where nn is the degree of v . The normal vector N can be estimated using the following classical formula [29].



n n −1

N =

−

Fig. 7. Compute a twist vector.



(vi − v) × (v(i+1) mod nn − v)

i=0

where operator [[⃗ x]] = ⃗ x/‖⃗ x‖2 . A tangent vector (see Figs. 5 and 6) can be computed by mapping the original vertex-to-vertex vector to the tangential plane π . That is, Vi = (vi − v) − ((vi − v) · N ) · N .

(10)

Fig. 8. Shape optimization for solids of revolutions.

Suppose that there exists a parabola on the normal plane spanned by N and Vi , which interpolates the vertices v and vi , and also the derivative vector Vi . We estimate its second-order derivative on v , and that is, C˜ ii = 2((vi − v) · N ) · N .

(11)

Twist (or the mixed second-order partial) derivative vectors Cij (j = (i + 1) mod nn ) can be estimated by computing the mean value of tangential movements in two directions. Shown in Fig. 7, C˜ ij =

1 2

((Vˆ j − Vj ) + (Vˆ i − Vi )),

(12)

where Vˆ i and Vˆ j are the corresponding tangent vectors of neighboring vertices. Since the estimated second-order derivatives C˜ ii and C˜ ij (in Eqs. (11) and (12)) may not satisfy G2 continuity conditions on vertex v , the least-square interpolation method presented in Section 3.6 needs to be applied. Finally, we can compute all second-order derivatives Cii and Cij for all vertices, fulfilling G2 continuity. This estimation algorithm works, but does not always lead to perfect shapes. As quad-dominant meshes are usually used for representing solids of revolution, we propose an enhanced algorithm to compute tangent vectors on four-degree vertices for revolutions. The spanning curves (bolded in Figs. 5 and 8), which

K.-L. Shi et al. / Computer-Aided Design 43 (2011) 145–160

151

a

b

Fig. 11. Decompose a T-shaped corner: (a) the facet to decompose contains even edges; (b) the other occasion, where the facet contains odd edges. Fig. 9. Compute new normal for revolutions.

central vertex to each middle point of spanning curves. vi denote the vertices of a multi-sided facet (see Fig. 10). For each vertex, we have a normal vector ni . Choose the arithmetical mean value of normal vectors of two end-points as the normal of a half-edge. That is, ni,j = 12 (ni + nj ), j = (n + 1) mod n. It is also the normal vector of point (vi + vj )/2 if the normal linearly varies along the edge. (0)

Fig. 10. Decompose a 5-sided facet: v0−4 are the original vertices, and m0−4 are the middle points of the spanning curves. After the edge-tangent extension, we have w0−4 . The final central vertex m is the arithmetical mean value of them.

are the boundaries of parametric patches, interpolate the tangent vectors and the second-order derivatives. For revolutions, shown in Fig. 8, the best choice of one curve is the dashed circle, or an approximate curve on its plane. The classical method incurs the petal-shaped curves, key problem of which is that we choose bad normal vectors. Taking one vertex v for example, we denote Ti = vi − v , where vi (i = 0, 1, 2, 3) are its neighboring vertices. Compute a new normal Nˆ = [[(T0 × T2 ) × (T1 × T3 )]]. Shown in Fig. 9, the classical normal N and the new computed one Nˆ span an approximately orthogonal plane to the dashed circle (in Fig. 9, and it is the plane spanned by T1 and T3 ). Denote the direction l = [[N × Nˆ ]]. If Ti is in this direction (|[[Ti ]] · l| > h), we must use the new normal Nˆ to compute the tangent vectors according to Eq. (10), where h is a prescribed threshold value, and we can choose h = 0.5 by default. Else, we choose the classical normal N. This enhancement ensures the approximation of the spanning curves to the circles of revolutions. 5. Quadrangulation and T-shaped corner manipulation The input mesh may contain multi-sided facets and T-shaped corners, while our method outputs a regular quad mesh without Tshaped corners. The main idea is to decompose all irregular facets.

(1)

Each edge has two normals ni,j and ni,j , with one per neighboring facet. We can choose their arithmetical mean value as the final edge (0) (1) normal n¯ i,j = 12 (ni,j + ni,j ). It follows that a cross-edge tangent vector is ti,j = (vi − vj ) × n¯ i,j , which is used to compute the central point of the multi-sided facet (see Fig. 10). The middle points mi can be set to the parametric middle points of the corresponding spanning curves (see the previous section). As the spanning curves interpolate the derivatives of their end-points, we can use a three-degree B-spline form as follows.

Υ (t ) =

5 −

N3,i (t )Υi ,

i=0

where N3,i (t ) are the B-spline basis functions, and their knot-vector is U = {0, 0, 0, 0, 1/3, 2/3, 1, 1, 1, 1}. v0 and v1 denote the two end-points of a spanning curve. Their tangential vectors along this curve are V0,p and V1,q , respectively. The corresponding secondorder derivative vectors are C0,pp and C1,qq (see Section 3). Let

 Υ (0) = v0 , Υ (1) = v1 , Υ ′ (0) = V0,p , Υ ′ (1) = −V1,q , Υ ′′ (0) = C0,pp , Υ ′′ (1) = C1,qq . We construct the B-spline curve (see Appendix E) and evaluate the middle point m0 = Υ (1/2) =

+

23 144

1 2

(v0 + v1 )

(V0,p + V1,q ) +

5 288

(C0,pp + C1,qq ).

All middle points mi can be computed analogously. The central vertex can be estimated by averaging all quasicentral-points wi . In our experience, wi can be chosen as follows. 1

wi = mi + (‖vj+ − vj ‖2 + ‖vi− − vi ‖2 ) · [[ti,j ]], 8

where i− = (i + n − 1) mod n and j+ = (j + 1) mod n. 5.2. T-shaped corners

5.1. Quadrangulate multi-sided facets

Section 3 requires that two adjacent tangent vectors Vi = ∂ Si (u, v)/∂ u| u=0 and Vj = ∂ Si (u, v)/∂v| u=0 (j = (i + 1) mod n)

A facet that contains n half-edges (n ̸= 4) is called a multisided facet. There are various methods to decompose this kind of facet. We use edge-tangential vectors to estimate a central vertex, and then split the original facet into n quads by connecting the

on a common n-sided vertex are linear independent. A T-shaped corner is a vertex that contains such linearly dependent tangent vectors (see Fig. 11). There are two occasions. For even-edge facet, we can connect the T-shaped corner to its opposite vertex vm (see

v=0

v=0

152

K.-L. Shi et al. / Computer-Aided Design 43 (2011) 145–160

have the following form.



L(u) = p(u)I ′ (u) + q(u)T (u) R(u) = r (u)I ′ (u) + s(u)T (u),

(13)

where p(u), q(u), r (u) and s(u) are scalar functions. Analogously, we denote the second-order cross-boundary derivative curves as follows.

  F (u) = p2 (u)I ′′ (u) + 2p(u)q(u)T ′ (u)   + q2 (u)E (u) + a(u)I ′ (u) + b(u)T (u) 2 G(u) = r (u)I ′′ (u) + 2r (u)s(u)T ′ (u)    + s2 (u)E (u) + c (u)I ′ (u) + d(u)T (u) Fig. 12. Spanning curves and corresponding cross-boundary derivatives.

Fig. 13. Cross-boundary derivatives along a common boundary.

Fig. 11(a)), and decompose the facet f into two sub-facets. If the facet contains odd edges, we need first to create a middle point vm on the opposite edge (see Fig. 11(b)) using a similar algorithm introduced in the previous subsection. Then, the facet turns out to be an even-edge one. Note that for the odd case, the generated vm is also a T-shaped corner of the neighboring facet (see Fig. 11(b)). So we need to continue applying the algorithm until there exists no Tshaped corners. In application, the T-shaped corner can be detected with a prescribed angular threshold ϵ , satisfying [[Vi ]] · [[Vj ]] < − cos ϵ . As mentioned before, T-shaped corners can be generated by handling multi-sided facets, and also by decomposing odd-edge facets while manipulating T-shaped corners. We usually first handle the input T-shaped corners and then decompose all multisided facets. Finally, we handle all T-shaped corners generated by the previous step. Note that both the two operations do not generate middle points for the same original edge twice. That means, the T-shaped corner manipulation algorithm is convergent. The original k-sided facet may be decomposed into k sub-facets at most. 6. Spanning curves and cross-boundary derivatives Coons B-spline surface construction requires compatible boundary curves and cross-boundary derivatives. The boundary curves, which are the so-called spanning curves (see the bolded curves in Fig. 12), have already been discussed in Section 5.1. In this section, we focus on how to generate first- and second-order cross-boundary derivative curves, fulfilling (a) G2 continuity on a common boundary; and (b) G2 compatibility conditions on common vertices. Shown in Fig. 13, we assume that two neighboring patches share a common boundary I (u). Suppose an imaginary normal cross-boundary derivative curve T (u). The cross-boundary derivative curves R(u) and L(u) of the upper and the lower patches can

(14)

where a(u), b(u), c (u) and d(u) are scalar functions, and E (u) is called the second-order normal cross-boundary derivative curve. T (u) and E (u) are the first- and the second-order cross-boundary derivative curves of the imaginary surface that G2 contacts the two adjacent surfaces on the common boundary, the corresponding along- and cross-boundary derivative vectors of which are orthogonal on two end-points, respectively. This representation ensures G2 continuity of the two adjacent patches on the common boundary because of the transitivity of geometric continuity of surfaces sharing a common boundary, which is strictly proved in Appendix A. Then, we only need to discuss the compatibility conditions on a common n-degree vertex (such as the right vertex in Fig. 13). The constructed Coons surface satisfies all boundary conditions exactly only if these compatibility conditions are fulfilled. Adding indices to the representations of I (u), L(u), R(u), F (u) and G(u), we derive the boundary conditions of these curves on a vertex, where the parametric value is set to 1.0. 6.1. Positional compatibility First we assume some constants on a prescribed vertex v . O denotes its positional vector, while N denotes its normal vector. Vi (i = 0, 1, . . . , n − 1) are the out-direction tangent vectors computed in Section 4. Cii and Cij (j = (i + 1) mod n) are the corresponding second-order derivative vectors. The positional compatibility claims that Ii (1) = O, for arbitrary i = 0, 1, . . . , n − 1. 6.2. Tangent compatibility The tangent compatibility claims that



Li (1) = Vj Ri (1) = Vk Ii′ (1) = −Vi ,

where k = (i+n−1) mod n and j = (i+1) mod n. Assume that the value of the normal cross-boundary derivative curve on this vertex is on the tangential plane (the plane spanned by all tangent vectors Vi ), and is orthogonal to the along-boundary derivative vector.4 That is, Ti (1) = [[N × Vi ]]. Combining the formula above with Eq. (13), we have



pi (1)Ii′ (1) + qi (1)Ti (1) = Vk ri (1)Ii′ (1) + si (1)Ti (1) = Vj .

4 In this paper, the along-boundary derivative means the derivative vector or derivative curve of a specified curve, comparing with the concept of cross-boundary derivative.

K.-L. Shi et al. / Computer-Aided Design 43 (2011) 145–160

153

Since Ii′ (1) and Ti (1) are already known, it follows that

6.5. Second-order partial twist compatibility

 [pi (1), qi (1)]T = σ˜ Ii′ (1),Ti (1) (Vk ) [ri (1), si (1)]T = σ˜ Ii′ (1),Ti (1) (Vj ),

The second-order partial twist compatibility comes from a common property of surfaces S: ∂ 3 S (u, v)/∂ u2 ∂v ≡ ∂ 3 Si (u, v)/∂v∂ u2 . This compatibility condition claims that

where σ˜ ⃗x,⃗y (⃗ v ) is a linear mapping operator that returns the coordinates of v ⃗ in the space spanned by ⃗x and y⃗. As Vi , Vj , Vk and Ti are coplanar, the operator always works.



Combining the formula above with Eq. (14), it follows that

 6.3. Twist compatibility The twist compatibility claims that R′i (1) = −Cij L′i (1) = −Cki .



Combining with Eq. (13), we have p′i Ii′ + q′i Ti + qi Ti′ = −Cki − pi Ii′′



ri′ Ii′ + s′i Ti + si Ti′ = −Cij − ri Ii′′ , where ‘(1)’ following each function name is omitted for short (e.g., pi stands for pi (1)), and Ii′′ (1) are prescribed in the following subsection. Then, we have



Thus, we can derive two equations concerning Ti′ (1).



Ai = Ai,N + Ai,T Bi = Bi,N + Bi,T ,

where

Ti′ = −(Cki + pi Ii′′ + p′i Ii′ + q′i Ti )/qi Ti′ = −(Cij + ri Ii′′ + ri′ Ii′ + s′i Ti )/si .



(15)

They lead to the same Ti′ (1) if all second-order derivative vectors satisfy G2 continuity conditions (see Section 3) on the common vertex. This property is proved in Appendix B.

The second-order tangent compatibility claims that

Ai,N = N · (Ai · N ), Bi,N = N · (Bi · N ),

Ai,T = Ai − Ai,N Bi,T = Bi − Bi,N .

Note that Ai,T and Bi,T are on the plane spanned by Ti (1) and Ii′ (1). For an arbitrary i, these planes are the same one, as they share the same normal vector N. Thus, the original equations can be decomposed into following two linear systems.



a′j Ij′ + b′j Tj = Ai,T ci′ Ii′ + d′i Ti = Bi,T ,



6.4. Second-order tangent compatibility

2pj qj Tj′′ + si Ti′′ + q2j Ej′ = Ai,N 2ri si Ti′′ + qj Tj′′ + s2i Ei′ = Bi,N .

For the first linear system, we have

 ′ ′T [aj , bj ] = σ˜ Ij′ ,Tj (Ai,T ) [cj′ , d′j ]T = σ˜ Ij′ ,Tj (Bi,T ).

Fi (1) = Ckk Gi (1) = Cjj Ii′′ (1) = Cii .



Combining with Eq. (14), we have



q2i Ei + ai Ii′ + bi T = Ckk − p2i Ii′′ − 2pi qi Ti′ s2i Ei + ci Ii′ + di T = Cjj − ri2 Ii′′ − 2ri si Ti′ .

The second-order total twist compatibility is according to the fact: ∂ 4 S (u, v)/∂ u2 ∂v 2 ≡ ∂ 4 Si (u, v)/∂v 2 ∂ u2 . It claims that

 [ai , bi ]T = σ˜ Ii′ ,T (Ckk − p2i Ii′′ − 2pi qi Ti′ ) [ci , di ]T = σ˜ Ii′ ,T (Cjj − ri2 Ii′′ − 2ri si Ti′ ).

G′′i (1) = Fj′′ (1).

Analogously, we have two expressions concerning Ei (1).

− 2pi qi Ti − ai Ii − bi T )/ − 2ri si Ti′ − ci Ii′ − di T )/ .

p2i Ii′′ ri2 Ii′′

′

′

For the second one, we note that there are 2n vector variables (Ti′′ and Ei′ ) and also 2n equations. It can be guaranteed to be a full-rank linear system since si , qi ̸= 0 for nonsingular cases. 6.6. Second-order total twist compatibility

Thus,

Ei = (Ckk − Ei = (Cjj −

2pj qj Tj′′ + si Ti′′ + q2j Ej′ + a′j Ij′ + b′j Tj = Ai 2ri si Ti′′ + qj Tj′′ + s2i Ei′ + ci′ Ii′ + d′i Ti = Bi ,

where Ai and Bi are expressions declared in Appendix C. In their definition, p′′i (1), q′′i (1), ri′′ (1), s′′i (1) and Ii′′′ (1) can be specified by evaluating curves that interpolate the boundary conditions derived before. For example, pi (1) and p′i (1) are computed in Sections 6.2 and 6.3. On the opposite end-point, pi (0) and p′i (0) are also computed similarly for the neighboring vertex. Thus, a threedegree Bézier or B-spline curve can be constructed interpolating these boundary conditions (see Section 6.7 and Appendix E). So the second-order derivative vector p′′i (1) can be evaluated. For the spanning curve I (u), we use three-degree-six-point uniform B-spline to interpolate the computed boundary conditions. Thus, Ii′′′ (1) also can be evaluated from this curve. Consequently, Ai and Bi can be evaluated. Ai and Bi are decomposed into two parts, respectively. That is,

 [p′i , q′i ]T = σ˜ Ii′ ,Ti (−Cki − pi Ii′′ ) [ri′ , s′i ]T = σ˜ Ii′ ,Ti (−Cij − ri Ii′′ ).



Fj′ (1) + R′′i (1) = 0 G′i (1) + L′′j (1) = 0.

q2i 2 si

They are compatible if G2 continuity conditions are fulfilled on the common vertex, and this also can be proved similarly to the proof of compatibility in the previous subsection.

Combine it with Eq. (14) and then expand the second-order derivatives. We have s2i Ei′′ − q2j E ′′ + ci′′ Ii′ + d′′i Ti − a′′j Ij′ − b′′j Tj = Ki , where Ki can be evaluated according to the expression shown in Appendix D. Note that a′′i (1), b′′i (1), ci′′ (1) and d′′i (1) also can be evaluated using the method presented in the previous subsection. Thus, the linear system has only n vector variables Ei′′ , as well as n equations. Analogously, it also can be guaranteed to be a full-rank linear system since si , qi ̸= 0 for nonsingular cases.

154

K.-L. Shi et al. / Computer-Aided Design 43 (2011) 145–160

Fig. 14. Example 1: (a) is the original mesh; (b) shows the derived quad mesh after quadrangulation; (c) is the rendered image of B-spline surfaces.

6.7. Construct Coons surfaces Previous subsections derive the following boundary conditions. – – – – –

Ii , I ′ , I ′′ Ti , T ′ , T ′′ Ei , E ′ , E ′′ pi , p′i , qi , q′i , ri , ri′ , si , s′i , ai , a′i , bi , b′i , ci , ci′ , di , d′i .

For each vertex, these boundary conditions can be computed. Since all the derivatives are considered as in-direction, while constructing curves, some value should be reversed. v0 and v1 denote the two end-points of a prescribed edge e. The corresponding indices of e are h and l, respectively. We construct the curves as follows. Ie = BSpline(v0 .Ih , −v0 .Ih′ , v0 .Ih′′ , v1 .Il , v1 .Il′ , v1 .Il′′ ) Te = BSpline(v0 .Th , −v0 .Th′ , v0 .Th′′ , −v1 .Tl , −v1 .Tl′ , −v1 .Tl′′ ) Ee = BSpline(v0 .Eh , −v0 .Eh′ , v0 .Eh′′ , v1 .El , v1 .El′ , v1 .El′′ ) pe = Bezier(−v0 .ph , v0 .ph , v1 .rl , v1 .rl ) ′

′

qe = Bezier(v0 .qh , −v0 .q′h , −v1 .sl , −v1 .s′l ) re = Bezier(−v0 .rh , v0 .rh′ , v1 .pl , v1 .p′l ) se = Bezier(v0 .sh , −v0 .s′h , −v1 .ql , −v1 .q′l ) ae = Bezier(−v0 .ah , v0 .a′h , v1 .cl , v1 .cl′ ) be = Bezier(v0 .bh , −v0 .b′h , −v1 .dl , −v1 .d′l ) ce = Bezier(−v0 .ch , v0 .ch′ , v1 .al , v1 .a′l ) de = Bezier(v0 .dh , −v0 .d′h , −v1 .bl , −v1 .b′l ). The constructed Ie , Te and Ee are three-degree uniform Bspline curves, and the others are in three-degree Bézier form (see Appendix E for definitions of the construction operators). According to Eqs. (13) and (14), we can write the expressions of the first-order and the second-order cross-boundary derivatives for the upper patch (see Fig. 13). Re (u) = re (u)Ie′ (u) + se (u)Te (u) Ge (u) = re2 (u)Ie′′ (u) + 2re (u)se (u)Te′ (u)

+ s2e (u)Ee (u) + ce (u)Ie′ (u) + de (u)Te (u).

Since all curves are in B-spline form, symbolic operators for Bspline curves [30] are applied to construct these curves. For a prescribed quadrangular facet, we compute the spanning curve Ie (u), the cross-boundary derivative Re (u), and the secondorder cross-boundary derivative Ge (u) for each edge. The final

Fig. 15. Local details of Example 1: (a) the quads after quadrangulation; (b) shows the rendered B-spline surfaces; continuous reflection lines in (c) shows their G2 continuity; (d)–(f) are their Gaussian, maximum and mean curvature charts. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

G2 Coons surface can be constructed using these boundary conditions [31]. Because Ge has the maximum degree nine, the Coons surface is in (9, 9)-degree B-spline form. Of course, we may only use Ie and Re to construct (6, 6)-degree Coons surfaces, which only preserve G1 continuity. In our experience, this also leads to better results than other G1 methods, as G2 conditions are considered while constructing these curves. 7. Examples In this section, several examples are presented to show feasibility of our method. Fig. 1 introduces a simplest example, where a cube is smoothed by B-spline surfaces, satisfying G2 continuity. Example 1 gives a more complicated model, containing 456 facets (see Fig. 14(a)). After quadrangulation, we have a 936facet quad mesh (see Fig. 14(b)). Note that the quadrangulation method preserves the feature of revolution, instead of just decomposing all non-quads. Smooth and G2 connected B-spline surfaces are shown in Fig. 14(c). Fig. 15 is the local view of this model. The continuous reflection lines (in Fig. 15(c)) and three

K.-L. Shi et al. / Computer-Aided Design 43 (2011) 145–160

155

Fig. 16. Example 2: Venus (126 facets). (a) is the original mesh; (b) and (c) show the continuous patches and their boundaries; and (d) shows the reflection lines of the model.

Fig. 17. Example 2: Venus (350 facets). (a) is the original mesh; (b) and (c) show the continuous patches and their boundaries; and (d) shows the reflection lines of the model.

Fig. 18. Example 3: smoothing non-coplanar quads. (a) is the original quads (note that some quads are non-coplanar); (b) shows the constructed B-spline patches and their isoparametric curves; (c) is the rendered image of these surfaces and (d) is the zebra stripes in polar form, which show their G2 continuity.

kinds of curvature graphs5 (in Fig. 15(d)–(f)) show global G2 continuity. Example 2 is the body of the classical model of Venus. Fig. 16(a) contains only 126 quad facets, which is extremely rough and discontinuous. Our algorithm is used to smooth these kinds of

5 These curvature graphs are generated by Pro ENGINEER 2002. The colortransition from red to blue shows the corresponding curvature values at each point on the surface, from their upper bound to the lower bound. The continuity of the surface color presents the curvature continuity of the surface.

rough meshes, in order to convert them into fair, continuous parametric surface patches. After smoothing, G2 continuous patches are generated. The continuous reflection lines in Fig. 16(c) underline their curvature continuity. Fig. 17 shows the result of smoothing a 360-facet quad mesh of the same model. More details are presented in this figure. Example 3 is a simple irregular model, containing 6 noncoplanar quads. The results in Fig. 18 show that our algorithm can handle this kind of irregular case, and the shape remains continuous and smooth. Example 4 shows a complete pipeline of our method. The prescribed mesh model is quad-dominant, with 4 non-quadrangular

156

K.-L. Shi et al. / Computer-Aided Design 43 (2011) 145–160

Fig. 19. Example 4: smoothing a cube with a trimmed corner. (a) is the input mesh (there are three pentagonal facets and one triangular facet); (b) shows the result of quadrangulation; (c) presents the spanning curves; (d) shows the colored patches; (e) shows the rendered B-spline surfaces and (f) presents the continuous reflection lines, which underline the G2 continuity of the patches.

Fig. 20. Example 5: smoothing a 24-facet torus-shaped quads. (a) is the original quads; (b) shows the constructed B-spline patches and their isoparametric curves; (c) is the rendered image of these surfaces; the reflection lines in (d) show their G2 continuity; (e) and (f) are Gaussian and mean curvature graphs.

facets. The quadrangulation step decomposes non-quadrangular facets and generates new vertices, interpolating the positional and normal vectors of neighboring vertices. After computing the structural vectors, spanning curves are generated. Then, we can construct the B-spline surfaces with computed Coons boundary conditions. Fig. 19 shows the results. Example 5 (see Fig. 20) shows a 24-facet torus-shaped model. The reflection lines and curvature graphs also emphasize their G2 continuity. In this example, we compute the tangent vector on each vertex using our optimized method for solids of revolution. Fig. 21 gives a comparison between this method and the classical

vector-projection method (which only projects the edges to the tangent plane for generating tangent vectors). The classical method leads to the petal-shaped unsmooth surfaces (see Fig. 21(a) and (b)), since the tangent vectors are constructed without considering the feature of revolution. Instead, the spanning curves generated by our method are much more like the natural circles on the expected torus (see Fig. 21(c)). Thus, we achieve a better shape (see Fig. 21(d)). Example 6 is a 80-facet mesh, approximating a torus. Our method can generate B-spline patches G2 smoothly. The final shape is close to the expected torus (Fig. 22).

K.-L. Shi et al. / Computer-Aided Design 43 (2011) 145–160

157

Fig. 21. A comparison between the classical vector-projection method, and our improved method in computing tangent vectors: (a) and (b) show the results by the classical method (note that the spanning curves are not always arc-shaped, and the generated surfaces are a bit twisted, although they are continuous exactly); (c) and (d) are our results.

Fig. 22. Example 6: a 80-facet torus approximation. (a) is the original quads; (b) shows the spanning curves; (c) shows the constructed B-spline patches and their isoparametric curves; (d) is the rendered image; the polar zebra stripes in (e) shows their G2 continuity; and (f) is their Gaussian curvature graph.

Table 1 Time cost (in seconds). The column ‘‘facets’’ shows the count of the facets of the model; the first and the third rows show the comparison time-costs of constructing G1 continuous patches using the same algorithm; the rows named ‘‘quad’’, ‘‘vectors’’, ‘‘curves’’ and ‘‘surfaces’’ show the time-costs of the steps of quadrangulation, constructing structural vectors, constructing derivative curves and constructing the final B-spline surfaces, respectively. Model 1

Example 1 (G ) Example 1 (G2 ) Example 2 (G1 ) Example 2 (G2 )

Facets

Quad.

Vectors

Curves

Surfaces

Total

456 456 350 350

0.58 0.61 1.07 1.14

1.14 1.16 1.83 1.88

0.52 4.56 4.01 7.54

2.48 7.96 0.53 20.68

4.72 14.29 7.44 31.24

Table 1 is a comparison of time-cost.6 It shows that the bottleneck is the construction of derivative curves and the final surface patches using symbolic operators of B-spline. So the algorithm spends more time achieving G2 continuity than only to create G1 continuous surfaces. The topology, not only the count of the facets, also determines the overall time-cost, since facets may be decomposed in the steps of quadrangulation and T-shaped-edge manipulation.

6 The algorithm was implemented in Microsoft C# 2.0, Intel T2300 1.6 GHz, 1G Memory, Windows XP. The values only present the comparison of steps and models, since the speed will be much faster if we re-implement the algorithm using C++.

8. Conclusions In this paper, we first present the compact form of G2 continuity conditions on vertices and prove that the maximum rank of the condition matrix is 2n − 3 by introducing its reflexivity, symmetry, transitivity and nonsingularity properties. Then, the least-square interpolating method using the homogeneous form of G2 continuity conditions is proposed to generate the second-order derivative vectors fulfilling G2 continuity conditions on vertices by adjusting prescribed derivative vectors. Spanning curves and their first- and second-order cross-boundary derivatives are constructed by computing all compatibility conditions on vertices. The final (9, 9)-degree Coons B-spline surfaces achieve global G2

158

K.-L. Shi et al. / Computer-Aided Design 43 (2011) 145–160

continuity exactly. The method is optimized for solids of revolution to get better shapes. The input mesh is only required to be closed, manifold and quad-dominant. We provide several examples to show its feasibility. This algorithm can be used for many occasions, such as smoothing a rough mesh, high-quality mesh rendering, converting a mesh into G2 continuous parametric surfaces in hybrid designing, etc. Acknowledgements The research was supported by Chinese 973 Program (2010CB328001), the National Science Foundation of China (60625202, 90715043), and Chinese 863 Program (2007AA040401). The fourth author was supported by the Fok Ying Tung Education Foundation (111070). The last author was supported by ANR-NSFC (60911130368). We wish to thank Peng Liu, Jia-Ting Chen and YaMei Wen of Tsinghua University for their valuable technical suggestions. Appendix A. Proof on G 1 and G 2 continuity on common boundary Assume that two parametric surfaces B and C share a common boundary I (u), and the cross-boundary derivative curves are L(u) and R(u), respectively. In application, the sign of the values of q(u) and s(u) will not change and there is no zero-value point for arbitrary u. Supposing the boundaries are not singular, we have the following theorem. Theorem 6. B and C are G1 continuously connected on the common boundary I (u) if there exist a curve T (u) and four scalar functions p(u), q(u), r (u) and s(u), satisfying the definition Eq. (13) and q(u), s(u) ̸= 0, ∀u. Proof. The G1 continuity condition on the common boundary [2– 5] is: there exist three scalar functions α(u), β(u) and γ (u), fulfilling

α(u)I (u) + β(u)L(u) + γ (u)R(u) ≡ 0. ′

And then,



L= p−

q rq  ′ I + R = α I ′ + β R. s s

Combining with Eq. (A.1), we choose



α = p − rq/s β = q/s.

(A.3)

Put the definitions of F and G (see Eq. (14)) into Eq. (A.2), and unite Eq. (A.3). After simplification, we have 2qr (qr − ps)I ′′ + 2(sp − rq)qR′

+ 2(qr − ps)qsT ′ + s2 γ I ′ + s2 δ T = 0. Since R′ = r ′ I ′ + rI ′′ + s′ T + sT ′ , we get

(2(sp − rq)qs′ + s2 δ)T + (s2 γ + 2(sp − rq)qr ′ )I ′ = 0. Thus, we can choose



δ = 2(qr − ps)qs′ /s2 γ = 2(ps − qr )qr ′ /s2 ,

so that Eq. (A.2) is identical. Accordingly, α(u), β(u), γ (u) and δ(u) always exist, and the G2 continuity conditions are fulfilled.  Appendix B. Proof on the compatibility of Ti′ First, we propose a theorem about the compatibility of Ti′ in Eq. (15). Assume that i is fixed, j = (i + 1) mod n and k = (i + n − 1) mod n, where n is the degree of the vertex. According to the G1 continuity condition, there exist three scalars α, β and γ , fulfilling

α Vi + β Vk + γ Vj = 0.

(B.1)

Theorem 8. Ti′ in Eq. (15) is compatible if

α Cˆ ii + β Cˆ ki + γ Cˆ ij = 0.

Here, we only need to choose

 α(u) = q(u)r (u) − p(u)s(u) β(u) = s(u) γ (u) = −q(u).

Proof. Section 6.2 declares that Li (1) = Vk and Ri (1) = Vj . Combining with Eq. (13), we have

Combining the formula above with Eq. (13), we have an identical equation. That means, the three scalar functions always exist. The G1 continuity condition is also fulfilled.  Suppose that the second-order cross-boundary derivative curves of B and C are F (u) and G(u), respectively. We have the following theorem.

si (1)Vk − qi (1)Vj = (qi (1)ri (1) − pi (1)si (1))Vi . Thus, we choose



α = pi (1)si (1) − qi (1)ri (1) β = si (1) γ = −qi (1).

(B.2)

Combining the two equations in (15), we have Theorem 7. B and C are G2 continuously connected on the common boundary I (u) if there exist two curves T (u) and E (u), and eight scalar functions p(u), q(u), r (u), s(u), a(u), b(u), c (u) and d(u), satisfying the definition Eqs. (13) and (14), and q(u), s(u) ̸= 0, ∀u. Proof. G2 continuity condition on the common boundary [2,3] is that there exist four scalar functions α(u), β(u), γ (u) and δ(u), fulfilling L = αI ′ + β R

(A.1)

F = α 2 I ′′ + 2αβ R′ + β 2 G + γ I ′ + δ R

(A.2)

where ‘(u)’ following each function name is omitted for short. According to Eq. (13), we have r 1 T = − I ′ + R. s s

si p′i Ii′ + si q′i Ti + si Cki + si pi Cii

= qi ri′ Ii′ + qi s′i Ti + qi Cij + qi ri Cii

(B.3)

where ‘(1)’ following each function name is omitted for short. Since the second-order derivatives Cij can be represented as Cij = N · Cˆ ij + Cij,T , where Cij,T lies on the common tangent plane spanned by Ii′ (1) and Ti (1). The G2 continuity conditions only constrain Cˆ ij . Cij,T can be chosen freely in order to satisfy Eq. (B.3). Thus, applying ·N to both sides of Eq. (B.3), we have

(si pi − qi ri )Cˆ ii + si Cˆ ki − qi Cˆ ij = 0. Combining it with Eq. (B.2), we have the final result

α Cˆ ii + β Cˆ ki + γ Cˆ ij = 0. 

K.-L. Shi et al. / Computer-Aided Design 43 (2011) 145–160

Theorem 9. α Cˆ ii +β Cˆ ki +γ Cˆ ij = 0, if the patches are G2 continuously connected on the common vertex. Proof. Since the patches are G2 continuous, it is equivalent to one C 2 single surface after proper reparameterizations. We suppose this surface is S (u(t1 , t2 , τ ), v(t1 , t2 , τ )), where u(t1 , t2 , τ ) = u0 + a1 t1 + a2 t2 + aτ v(t1 , t2 , τ ) = v0 + b1 t1 + b2 t2 + bτ ,



Appendix E. Definitions of operators BSpline and Bezier The operator BSpline(v0 , d0 , e0 , v1 , d1 , e1 ) constructs a threedegree-six-control-point uniform B-spline curve c, satisfying c (0) = v0 ,

c ′ (0) = d0 ,

c (1) = v1 ,

c ′ (1) = −d1 ,

c ′′ (0) = e0 , c ′′ (1) = d1 .

The parameters of c are listed as follows [31].

fulfilling Vk = ∂ S (u0 , v0 )/∂ t1 = a1 ∂ S (u0 , v0 )/∂ u + b1 ∂ S (u0 , v0 )/∂v Vj = ∂ S (u0 , v0 )/∂ t2 = a2 ∂ S (u0 , v0 )/∂ u + b2 ∂ S (u0 , v0 )/∂v Vi = ∂ S (u0 , v0 )/∂τ = a∂ S (u0 , v0 )/∂ u + b∂ S (u0 , v0 )/∂v,



and S (u0 , v0 ) is the original common vertex. Choose

 α = a1 b 2 − a2 b 1 β = a2 b − ab2 γ = ab1 − a1 b, so that Eq. (B.1) is always fulfilled. According to the property of C 2 continuity, there exist second-order total and mixed partial derivatives on S (u0 , v0 ). Then, we compute the mixed partial derivatives of S on (u0 , v0 ). That is,

 2 ∂ S (u0 , v0 )/∂τ 2 = a2 Duu + 2abDuv + b2 Dvv ∂ 2 S (u , v )/∂ t1 ∂τ = a1 aDuu + (a1 b + ab1 )Duv + b1 bDvv  2 0 0 ∂ S (u0 , v0 )/∂ t2 ∂τ = a2 aDuu + (a2 b + ab2 )Duv + b2 bDvv , where

 c .degree = 3,   c .knots = {0, 0, 0, 0, 1/3, 2/3, 1, 1, 1, 1}, c [0] = v0 , c [1] = v0 + d0 /9, c [2] = v0 + d0 /3 + e0 /27,   c [5] = v1 , c [4] = v1 + d1 /9, c [3] = v1 + d1 /3 + e1 /27. Square brackets here stand for the indices of the control points. The Bézier operator Bezier(v0 , d0 , v1 , d1 ) constructs a three-degreefour-control-point Bézier curve b, satisfying b(0) = v0 ,

b′ (0) = d0 ,

b(1) = v1 ,

b′ (1) = −d1 .

The parameters of b are listed as follows.

 b.degree = 3,   b.knots = {0, 0, 0, 0, 1, 1, 1, 1}, b[0] = v0 , b[1] = v0 + d0 /3,   b[3] = v1 , b[2] = v1 + d1 /3. References

 Duu = ∂ 2 S (u0 , v0 )/∂ u2 D = ∂ 2 S (u0 , v0 )/∂ u∂v  uv Dvv = ∂ 2 S (u0 , v0 )/∂v 2 . According to the definition, Cki is the derivative of Vk in the direction of Vi , Cij is the derivative of Vj in direction Vi , and Cii is the second-order derivative vector in Vi direction. So we have

 Cii = ∂ 2 S (u0 , v0 )/∂τ 2 C = ∂ 2 S (u0 , v0 )/∂ t1 ∂τ  ki Cij = ∂ 2 S (u0 , v0 )/∂ t2 ∂τ . Combining all equations above, we have

α Cˆ ii + β Cˆ ki + γ Cˆ ij ≡ 0.  Since we have ensured that all common vertices are G2 continuous, according to Theorems 8 and 9, Ti′ (1) is always compatible. Appendix C. The definitions of Ai and Bi Ai = −(2pj p′j + aj )Ij′′ − 2ri′ Ii′′ − (2p′j qj + 2pj q′j + bj )Tj′

− 2s′i Ti′ − 2qj q′j Ej − p2j Ij′′′ − ri Ii′′′ − ri′′ Ii′ − s′′i Ti Bi = −(2pj p′j + aj )Ij′′ − 2ri′ Ii′′ − (2p′j qj + 2pj q′j + bj )Tj′

− 2s′i Ti′ − 2qj q′j Ej − p2j Ij′′′ − ri Ii′′′ − ri′′ Ii′ − s′′i Ti . Appendix D. The definition of Ki 2

159

(4)

Ki = 2(p′ j + pj p′′j + a′j )Ij′′ + (4pj p′j + aj )Ij′′′ + p2j Ij

+ 2(p′′j qj + 2p′j q′j + pj q′′j + b′j )Tj′ + (4p′j qj Tj′′ + bj + 4pj q′j )Tj′′ + 2pj qj Tj′′′ + 4qj q′j Ej′ + 2(q′j q′j + qj q′′j )Ej − 2(ri′2 + ri ri′′ + ci′ )Ii′′ − (4ri ri′ + ci )Ii′′′ − ri2 Ii(4) − 2(ri′′ si + 2ri′ s′i + ri s′′i + d′i )Ti′ − (4ri′ si + 4ri s′i + di )Ti′′ − 2ri si Ti′′′ − 4si s′i Ei′ − 2(s′i s′i − si s′′i )Ei .

[1] Zheng JM, Wang GZ, Liang YD. GC n continuity conditions for adjacent rational parametric surfaces. Computer Aided Geometric Design 1995;12:111–29. [2] Ye XZ. The Gaussian and mean curvature criteria for curvature continuity between surfaces. Computer Aided Geometric Design 1996;13:549–67. [3] Ye XZ, Liang YD, Nowacki H. Geometric continuity between adjacent Bézier patches and their constructions. Computer Aided Geometric Design 1996;13: 512–48. [4] Liang XZ, Che XJ, Li Q. G2 continuity conditions for two adjacent B-spline surfaces. In: Proceedings of the geometric modeling and processing. 2004. [5] Che XJ, Liang XZ, Li Q. G1 continuity conditions of adjacent NURBS surfaces. Computer Aided Geometric Design 2005;22:285–98. [6] Peters J. Joining smooth patches around a vertex to form a C k surface. Computer Aided Geometric Design 1992;9:387–411. [7] Jones AK. Nonrectangular surface patches with curvature continuity. Computer-Aided Design 1988;20(6):325–35. [8] Catmull E, Clark J. Recursively generated B-spline surfaces on arbitrary topological meshes. Computer-Aided Design 1978;10:350–5. [9] Loop C. Second order smoothness over extraordinary vertices. In: Symposium on geometry processing. 2004. p. 169–78. [10] Karčiauskas K, Myles A, Peters J. A C 2 polar jet subdivision. In: Proceedings of the fourth eurographics symposium on geometry processing. 2006. p. 173–80. [11] Karčiauskas K, Peters J. Bicubic polar subdivision. ACM Transactions on Graphics 2007;26: 14:1–6. [12] Karčiauskas K, Peters J. Bi-3 C 2 polar subdivision. ACM Transactions on Graphics 2009;28:48:1–48:12. [13] Zorin D, Schroder P, Sweldens W. Interpolating subdivision for meshes with arbitrary topology. In: Proceedings of the 23rd annual conference on computer graphics and interactive techniques. 1996. p. 189–92. [14] Peters J. Curvature continuous spline surfaces over irregular meshes. Computer Aided Geometric Design 1996;13:101–31. [15] Bajaj CL, Ihm I. Smoothing polyhedra using implicit algebraic splines. In: Proceedings of the 19th annual conference on computer graphics and interactive techniques. 1992. p. 79–88. [16] Peters J. Smooth mesh interpolation with cubic patches. Computer-Aided Design 1990;22:109–20. [17] Peters J. Smooth interpolation of a mesh of curves. Constructive Approximation 1991;7:221–46. [18] Cho DY, Lee KY, Kim TW. Interpolating G1 Bézier surfaces over irregular curve networks for ship hull design. Computer-Aided Design 2006;38:641–60. [19] Ye XZ. Curvature continuous interpolation of curve meshes. Computer Aided Geometric Design 1997;14:169–90. [20] Hermann T. G2 interpolation of free form curve networks by biquintic Gregory patches. Computer Aided Geometric Design 1996;13:873–93. [21] Li H, Liu SQ. Local interpolation of curvature-continuous surfaces. ComputerAided Design 1992;24:491–503. [22] Du WH, Schmitt FJM. On the G1 continuity of piecewise Bézier surfaces: a review with new results. Computer-Aided Design 1990;22:556–73. [23] Peters J. Local smooth surface interpolation: a classification. Computer Aided Geometric Design 1990;7:191–5.

160

K.-L. Shi et al. / Computer-Aided Design 43 (2011) 145–160

[24] Shi XQ, Wang TJ, Yu PQ. A practical construction of G1 smooth biquintic B-spline surfaces over arbitrary topology. Computer-Aided Design 2004;36: 413–24. [25] Hahmann S, Bonneau GP. Triangular G1 interpolation by 4-splitting domain triangles. Computer Aided Geometric Design 2000;17:731–57. [26] Hahmann S, Bonneau GP. Polynomial surfaces interpolating arbitrary triangulations. IEEE Transactions on Visualization and Computer Graphics 2003;9: 99–109.

[27] Tong WH, Kim TW. High-order approximation of implicit surfaces by G1 triangular spline surfaces. Computer-Aided Design 2009;41:441–55. [28] Greville BIAT. Generalized inverses. Springer; 2003. [29] Piegl LA, Tiller W. Filling n-sided regions with NURBS patches. The Visual Computer 1999;15:77–89. [30] Piegl LA, Tiller W. Symbolic operators for NURBS. Computer-Aided Design 1997;29(5):361–8. [31] Piegl LA, Tiller W. The NURBS book. 2nd ed. Springer; 1997.