An Unified λ-subdivision Scheme for Quadrilateral Meshes with Optimal Curvature Performance in Extraordinary Regions

We propose an unified λ-subdivision scheme with a continuous family of tuned subdivisions for quadrilateral meshes. Main subdivision stencil parameters of the unified scheme are represented as spline functions of the subdominant eigenvalue λ of respective subdivision matrices and the λ value can be selected within a wide range to produce desired properties of refined meshes and limit surfaces with optimal curvature performance in extraordinary regions. Spline representations of stencil parameters are constructed based on discrete optimized stencil coefficients obtained by a general tuning framework that optimizes eigenvectors of subdivision matrices towards curvature continuity conditions. To further improve the quality of limit surfaces, a weighting function is devised to penalize sign changes of Gauss curvatures on respective second order characteristic maps. By selecting an appropriate λ, the resulting unified subdivision scheme produces anticipated properties towards different target applications, including nice properties of several other existing tuned subdivision schemes. Comparison results also validate the advantage of the proposed scheme with higher quality surfaces for subdivision at lower λ values, a challenging task for other related tuned subdivision schemes.

Fig. 1.The proposed unified -subdivision scheme covers a continuous family of tuned subdivisions in , the subdominant eigenvalue of the respective subdivision matrix.Main subdivision coefficients,  2 () and  3 (), of the proposed scheme are represented by spline functions of .Subdivision schemes with different properties can be conveniently extracted by proper selection of  values.The  2 () and  3 () functions for  = 6 are illustrated (with different vertical shift and in different vertical scale for illustration only), and several typical selections of  to accommodate different applications are also highlighted.
We propose an unified -subdivision scheme with a continuous family of tuned subdivisions for quadrilateral meshes.Main subdivision stencil parameters of the unified scheme are represented as spline functions of the subdominant eigenvalue  of respective subdivision matrices and the  value can be selected within a wide range to produce desired properties of refined meshes and limit surfaces with optimal curvature performance in extraordinary regions.Spline representations of stencil parameters are constructed based on discrete optimized stencil coefficients obtained by a general tuning framework that optimizes eigenvectors of subdivision matrices towards curvature continuity conditions.To further improve the quality of limit surfaces, a weighting function is devised to penalize sign changes of Gauss curvatures on respective second order characteristic maps.By selecting an appropriate , the resulting unified subdivision scheme produces anticipated properties towards different target applications, including nice properties of several other existing tuned subdivision schemes.Comparison results also validate the advantage of the proposed scheme with higher quality surfaces for subdivision at lower  values, a challenging task for other related tuned subdivision schemes.

INTRODUCTION
Subdivision schemes provide an elegant solution for the representation of models with arbitrary topology.With a given set of subdivision rules, an input mesh is subdivided to generate a series of refined meshes that converge to a smooth limit surface.Subdivision schemes relax the rigid grid structure of B-spline control meshes and allow the use of extraordinary vertices (EVs) on a quadrilateral control mesh with the number of incident edges meeting at an internal extraordinary vertex (EV) being three or more than four.The number of edges connected to a vertex is usually called its valence.A popular subdivision scheme for quadrilateral meshes is the Catmull-Clark subdivision [Catmull and Clark 1978] which has been integrated in various modeling software, such as Rhino, Blender and Maya, for model definition using unstructured meshes with arbitrary topology.
The properties of Catmull-Clark subdivision at extraordinary vertex positions, however, can be further improved for some practical applications.First, Catmull-Clark subdivision cannot produce bounded curvature near extraordinary points [Sabin et al. 2003], which is a necessary condition for producing high quality limit surfaces with  2 / 2 continuity conditions.Second, Catmull-Clark surfaces might exhibit kinked highlight lines indicating less competitive surface qualities.Third, Catmull-Clark subdivision cannot reproduce convex shapes [Karciauskas et al. 2004] for generic data when the valence is  ≥ 5. To remedy these problems, a simple idea is to directly modify subdivision rules so that the resulting subdivision scheme preserves desired properties.Such a method is usually referred to as subdivision tuning.Subdivision rules as linear operations of vertices can be written in a matrix form, and the operation of local mesh subdivision can be performed through a matrix multiplication operation between the subdivision matrix and the array of local control vertices.As subdivision surfaces can be considered as infinitely refined control meshes, surface properties can thus be determined by the eigenstructure of respective subdivision matrices.Tuning of subdivision schemes is thus to modify subdivision rules so that the eigenstructure of the subdivision matrix follows specific patterns that lead to desired properties of limit surfaces.
A key parameter in the eigenstructure analysis of subdivision matrices is the subdominant eigenvalue , i.e., the twofold second largest eigenvalue that has paramount importance in subdivision tuning.Tuned subdivision schemes usually produce different  values suitable for different target applications.Intuitively,  quantifies the speed of contraction for 1-ring vertices in the process of subdivision.As illustrated in Fig. 1, a lower  value leads to faster contraction of 1-ring vertices.Existing tuned subdivision schemes with high surface qualities, e.g.[Augsdörfer et al. 2006] and [Ma and Ma 2018], usually result in polar artifacts in the mesh structure which might cause inconveniences for rendering [Augsdörfer et al. 2009] due to slightly higher  values.At  = 0.5, polar artifacts could be eliminated with uniform refined meshes [Ma and Ma 2019b], see Fig. 1, at the expense of slightly reduced limit surface qualities.So in practice, subdivision schemes with intermediate  values might be necessary for different modeling and graphics applications to balance between different desired priorities.In recent years, subdivision schemes have also been widely applied in engineering analysis.The use of subdivision schemes with a lower  value usually produces analysis results with faster convergence of solution errors [Ma and Ma 2019a;Wei et al. 2021], and the use of subdivision schemes at different  values would also be preferred for producing different target analysis solutions.
These observations motivate us to develop a continuous family of -subdivision schemes that cover a wide range of  values to accommodate different applications.Inspired by the clear geometric meaning of the subdominant eigenvalue , we use  as an intuitive indexing parameter to extract respective subdivision schemes with specific subdominant eigenvalue .For ease of implementation, main subdivision coefficients are written as spline functions of .By selecting an appropriate  value, the respective tuned subdivision scheme can be extracted by simple spline evaluation and followup arithmetic calculations.As to the selection of appropriate  values for different target applications, results of several featured -subdivisions are provided as illustrations for guiding the selection.Main features and contributions of this work are as follows.
• This paper presents a novel unified -subdivision scheme with a continuous family of tuned subdivision schemes unified under a single parameter  that has explicit and intuitive meanings for various practical applications.• Main subdivision coefficients of the unified -subdivision scheme are represented as spline functions in .For a given  value, the full set of respective subdivision parameters can be easily extracted through spline evaluations and further simple arithmetic calculations.• The proposed unified -subdivision scheme includes tuned subdivisions with better or similar properties to that in [Ma and Ma 2018], [Ma and Ma 2019b] and [Ma and Ma 2019a].
In addition, one can extract tuned subdivision schemes in a continuous wide range of  for various practical applications in modeling, computer graphics and engineering analysis.• The proposed unified subdivision scheme generates tangent plane continuous limit surfaces at extraordinary points with optimal or near optimal curvature performances for any feasible  by optimization towards curvature continuity conditions.To accommodate wider applications, the feasible domain of  is expanded to lower  region by proper relaxation of bounded curvature constraints.• We also construct a weighted objective function to further improve surface qualities at lower  values, e.g., [Ma and Ma 2019a].The new objective function includes weighting functions that penalize undesired sign changes of Gauss curvatures of second order maps for better surface performances.
The rest of the paper is organized as follows.Some further preliminaries on subdivision-based modeling and other related works can be found in Section 2. To construct a continuous family of the proposed subdivision schemes, we first perform subdivision tuning at a set of discrete  values in Section 3. Spline representations of main subdivision coefficients in , the so-called spline stencils, are further constructed from the resulting discrete tuned subdivision schemes in Section 4. Subdivisions at some featured  values with method for the evaluation of the full set of stencil parameters at an arbitrary  value are also highlighted in Section 4. Examples and some further discussions with comparisons can be found in Section 5. Conclusions and future works are discussed in Section 6.

RELATED WORK 2.1 Subdivision prerequisites
Subdivision schemes work by recursive refinement of input meshes producing a sequence of refined meshes converging to smooth limit surfaces.Each of refined vertices is computed as linear combinations of vertices in the old mesh, and coefficients that describe the influence from old vertices to a new vertex is referred to as a subdivision stencil [Sabin et al. 2007], as illustrated in Fig. 2. Subdivision of a given mesh can be performed by matrix multiplication, and the limit properties of subdivision surfaces can be analyzed through the spectrum of subdivision matrices.With proper organization of local control vertices, the respective subdivision matrix S is usually block-circulant that can be transformed to a block-diagonal form Ŝ through the help of Discrete Fourier transform [Ball and Storry 1988].
Eigenvalues of subdivision matrix S can thus be computed by solving eigenvalues of the diagonal submatrices in Ŝ which is much easier, and if the eigenvalue   is an eigenvalue of the -th diagonal block in Ŝ, then   has a Fourier index of , denoted by F (  ) = { } [Peters and Reif 2008].For a general subdivision scheme with smooth limit surfaces, the first three eigenvalues in descending order should be 1, , , and  < 1 is called the subdominant eigenvalue.Similarly, the one after  is called the subsubdominant eigenvalue.There are two right eigenvectors ì  1 and ì  2 for the subdominant eigenvalue , which can be considered as lists of two coordinates of vertices in R 2 .The set of vertices form the natural configuration of the subdivision scheme, see Fig. 1, and the spline rings defined by the natural configuration is referred to as the characteristic map [Reif 1995].Similarly, higher order maps can be defined if we consider the eigenvector for the subsubdominant eigenvalue as the third coordinate, see Fig. 4(b).Limit properties, e.g.,  1 continuity, can be analyzed through analysis of characteristic maps, while higher order properties are related to higher order maps.For a comprehensive tutorial on subdivision surfaces, please refer to [Zorin 2000] [Peters andReif 2008].

Unified subdivision schemes
With the introduction of Catmull-Clark subdivision [Catmull and Clark 1978] and Doo-Sabin subdivision [Doo and Sabin 1978], a wide range of subdivision schemes have been further proposed, including Loop subdivision [Loop 1987], √ 2-subdivision [Li et al. 2004], √ 3subdivision [Kobbelt 2000], Butterfly subdivision [Dyn et al. 1990], Kobbelt subdivision [Kobbelt 1996], and so on.For further introduction on subdivision schemes, please refer to [DeRose et al. 1998;Reif and Sabin 2019;Sabin 2005] and references therein for information.These subdivision schemes are mostly further generalization of certain class of splines or other basis functions.In the literature, there are also families of unified subdivision schemes that are proposed to unify and further generalize various individual classes of subdivision schemes, aiming at (1) subdivision with arbitrary continuity in regular regions [Stam 2001] [Zorin and Schröder 2001] [Deng and Ma 2013]; (2) generalizations of interpolatory and approximating subdivisions [Maillot and Stam 2001][Shi et al. 2008][Zhang et al. 2019][Novara and Romani 2016]; (3) subdivision for mixed triangle/quad meshes [Stam and Loop 2003][Lin et al. 2013] [Peters and Shiue 2004]; and (4) unifications for existing non-stationary schemes [Zheng and Zhang 2017].In this work, we present a unified -subdivision scheme that can be considered as a generalized unification of several existing tuned subdivision schemes, such as [Ma and Ma 2018], [Ma and Ma 2019b], and [Ma and Ma 2019a] for different target applications.

Limit surface properties of subdivision schemes
Properties of subdivision surfaces can be analyzed from the eigenstructures of the respective subdivision matrices [Doo and Sabin 1978].A sufficient condition for  1 continuity at extraordinary positions is given in [Reif 1995] that requires two identical subdominant eigenvalues, denoted by , and regular and injective characteristic maps.It has been verified that Catmull-Clark, Doo-Sabin, and Loop schemes, satisfy the  1 condition and thus are  1 continuous at extraordinary positions [Peters and Reif 1998] [ Umlauf 2000].For surfaces to be   continuous at extraordinary positions, a lower bound for the degree of respective polynomial patches is given in [Prautzsch and Reif 1999], and a stationary subdivision scheme requires a minimum degree of 6 to be curvature continuous at extraordinary positions.The Catmull-Clark scheme is thus only  1 continuous at extraordinary positions, and the curvatures are unbounded [Karciauskas et al. 2004].Criteria for bounded curvature are presented in [Doo and Sabin 1978] and subdivision rules can be tuned to satisfy such properties.Although curvatures may not be bounded for Catmull-Clark subdivision near extraordinary positions, the principle curvatures are square-integrable [Reif and Schröder 2001], indicating that the scheme can be used for engineering analysis.There are also other constructions for  2 continuity at extraordinary positions, e.g., guided subdivision [Karčiauskas and Peters 2007][ Karčiauskas and Peters 2018] and patchwork methods [Loop and Schaefer 2008] [Ma and Ma 2020] [Yang et al. 2023].These schemes may produce higher-quality limit surfaces at the cost of higher complexity for implementation.Patchwork methods may not always produce nested limit surfaces.

Tuning of subdivision schemes
Subdivision coefficients can be tuned to have desired properties.The Catmull-Clark scheme and the Loop scheme are tuned for bounded curvature in [Sabin 1991] and [Loop 2002], respectively, and bounded curvature is also achieved for non-uniform subdivisions [Cashman et al. 2009b].To incorporate a wider range of possible desired properties, a general tuning framework through optimization of eigenvectors is proposed in [Barthe and Kobbelt 2004] that produces surfaces with bounded curvature and alleviated polar artifact.Performances of limit surfaces at extraordinary positions can be further improved by optimizing respective eigenvectors towards curvature continuity conditions while maintaining minimum curvature variations [Ma and Ma 2018] and least polar artifact [Ma and Ma 2019b].In [Augsdörfer et al. 2006], the Gauss curvature near extraordinary positions is minimized through optimization of Gauss curvatures on a set of representative central surfaces.

Tuning for applications in Isogeometric Analysis (IGA)
With the popularity of subdivision in engineering analysis, tuning of subdivision schemes for analysis purposes have attracted much attention.In [Zhang et al. 2018], the Catmull-Clark scheme is tuned through optimization of second order characteristic maps using a thin-plate energy, and the tuned scheme yields a significant reduction of solution errors in  2 -norms.For higher convergence rates in isogeometric analysis (IGA), a lower subdominant eigenvalue of  = 0.39 is used in [Ma and Ma 2019a] to increase local mesh densities near extraordinary vertices, and optimal convergence rates are observed for  2 -approximation problems.Similar ideas are also used for tuning of non-uniform subdivisions in [Li et al. 2019] [Wei et al. 2021] and Loop subdivision in [Kang et al. 2022].In general, these schemes require  values lower than 0.5 which might be undesirable for high quality limit surfaces.Further discussions on subdivision for IGA applications can also be found in [Dietz et al. 2023].

Other tuning and construction schemes
To balance between surface qualities and engineering analysis, a tuned scheme [Wang and Ma 2023] further relaxes Catmull-Clark rules for 2-ring refined vertices so that additional degrees of freedom are introduced in optimization for improved surface qualities at the cost of increased complexity in implementation for both tuning and later applications.In the literature, one can also find special rules for 2-ring refined vertices in [Karčiauskas and Peters 2023a][Karčiauskas and Peters 2023b][Karčiauskas and Peters 2022] that are constructed based on the idea of guided subdivision [Karčiauskas and Peters 2007] with refined vertices that are dependent on all 3-ring vertices, leading to much more degrees of freedom for improving desired limit surface properties and with well behaved limit surface quality benefiting both modeling and relevant applications.

CONSTRUCTION OF SUBDIVISION STENCILS AT SPECIFIED 𝜆 VALUES
In this section, we present a general tuning framework to compute stencil coefficients through optimization towards curvature continuity conditions at extraordinary positions.The tuning framework is similar to that in [Ma and Ma 2018], but a different weighted objective function with integrated penalty for reducing sign changes in Gauss curvature of the second order characteristic maps is used for stencil optimization in this work, which leads to improved properties of limit surfaces, especially for subdivision at lower  values.

Subdivision stencils and subdivision matrices
The proposed family of subdivision schemes generalizes Catmull-Clark subdivision [Catmull and Clark 1978] by modifying subdivision rules for 1-ring refined vertices.The respective stencil coefficients are written as functions of , the subdominant eigenvalue of respective subdivision matrices.The symbolized subdivision stencils are illustrated in Fig. 2. For simplicity, we omit the variable  of subdivision coefficients sometimes, and use  2 ,  3 , etc., instead.
The subdivision stencils in Fig. 2 are for quadrilateral meshes of arbitrary topology with separated extraordinary vertices.In cases of input meshes with general topological structure or quadrilateral meshes with connected extraordinary vertices, we can subdivide the mesh once using Catmull-Clark rules [Catmull and Clark 1978] to separate extraordinary vertices, and the proposed unified subdivision scheme can thus be applied.
Based on the irregular rules in Fig. 2 and regular rules for Catmull-Clark subdivision, local subdivision matrices, denoted by S, can be constructed to perform local mesh subdivision, written as where ì   is the vector of control vertices at level  = 0, 1, 2, . ... The respective eigenvalues and eigenvectors of the subdivision matrix S can be explicitly written as functions of subdivision stencils in Fig. 2, which will be further used for subsequent tuning of desired subdivision stencils.A detailed formulation of subdivision matrices with eigenstructure analysis can be found in the attached supplementary materials (Section A).
e 1 (6) e 2 ( 6) Hollow circles are regular vertices from the original mesh and red pentagram represents the extraordinary vertex in the original mesh.Red, blue and green solid circles are the refined face, edge, and vertex vertices, respectively.Note that all subdivision coefficients can be written as functions of subdominant eigenvalues , respectively, and details will be given in Section 4.

Expected properties for subdivision schemes
For a given subdivision scheme, properties of limit surfaces are closely related to eigenstructures of the respective subdivision matrices.In the following, we highlight relevant conditions and requirements of respective eigenstructures for meeting expected properties of desired subdivision schemes, which will be used for later stencil tuning and optimization.For clarity, eigenvalues of the local subdivision matrix S are arranged with descending absolute values as 3.2.1 Convergence.The subdivision scheme converges and is affine invariant [Doo and Sabin 1978] if and only if 3.2.2 1 continuity.The subdominant eigenvalue  is positive and has algebraic and geometric multiplicity two, with Fourier indices 1,  − 1, and the characteristic map defined by the corresponding eigenvectors is regular and injective [Peters and Reif 1998].

Convex hull property.
All subdivision stencils in Fig. 2 should be nonnegative, i.e., Note that the variable  in these stencils is omitted for brevity.
3.2.5  2 continuity.The necessary and sufficient conditions [Reif 2007][Peters andReif 2008] for a subdivision surface to be  2 at an extraordinary position for all generic initial data is given as The second condition can be written as where   is the eigenfunctions corresponding to   , and  ≤ 5 due to the linear independence.At least triple subsubdominant eigenvalues with Fourier indices {0, 2,  − 2} are required to produce all possible basic shapes, which are always satisfied by our scheme.

Stencil optimization with desired properties
The subdivision tuning can be performed through an optimization framework with objective functions being the main expected properties subjecting to a group of necessary constraints for stencil coefficients to meet other desired properties of subdivision schemes.
3.3.1 Selection of free variables for optimization.Some of the aforementioned properties can be strictly satisfied, such as convergence of refined meshes,  1 continuity, the convex hull property, and bounded curvature.These properties usually lead to hard constraints in the optimization of stencil coefficients.
As illustrated in Fig. 2, there are nine subdivision coefficients, i.e.,  1 ,  2 ,  3 ,  1 ,  2 ,  1 ,  2 , , .If all the hard constraints are satisfied, seven degrees of freedom will be eliminated, leaving only two independent parameters for optimization.In this work, we choose  2 and  3 , the stencil coefficients for face vertices, as independent variables, and other coefficients for prescribed  and  can be written as functions of  2 and  3 in terms of Detailed expression for Eq. ( 7) can be found in Eqs. ( 25)-( 28) and in Section 4.5, and hard constraints can be transformed into parameter bounds of the following form where details are given in supplementary materials (Section B).

Constraint relaxation for bounded curvature at lower 𝜆 values.
The lower bounds   for  of valences  = 3, . . ., 20 to produce schemes with bounded curvature are illustrated in Fig. 3.If bounded curvature is strictly required, one can only select  in a narrow range, especially for higher valences.To further extend feasible domains for , the bounded curvature constraint is relaxed as where  ≥ 1 is a relaxation coefficient for the bounded curvature criterion of  =  2 .In this work, the value of  is computed as where   is the lower bound for  to satisfy other desired properties and Δ is a small tolerance necessary in eigenstructure analysis.The lower bounds   for all valences have consistent forms with explicit formula in Eq. ( 47) given in the supplementary materials (Section B).A typical  −  curve is illustrated in Fig. 6(e) for valence  = 6.We empirically select Δ = Δ 0 = 10 −3 at  ≤   and Δ = 0 at  ≥   + Δ 0 in this work, where Δ 0 = 10 −2 .Further investigations regarding the selection of  for producing well behaved subdivisions will be conducted in our future work.

3.3.3
The objective function for optimization.The requirements for  2 continuity are adapted to construct the objective function for stencil optimization.The  2 continuity condition in Eq. ( 6) requires computations of eigenfunctions   for eigenvalues  and .The eigenfunction   has intuitive geometric meanings.If we consider eigenvectors of  as arrays of coordinates for control vertices, the resulting subdivision surface is just the span of eigenfunctions ( 1 , 2 ), also written as (, ) in Fig. 4(a).Similarly, if we further consider one of the eigenvectors of  as the third coordinate for respective control vertices over (, ) in Fig. 4(a), the resulting subdivision surface will be the span of eigenfunctions ( 1 , 2 ,  ),  = 3, 4, 5, also written as (, , ) in Fig. 4(b), corresponding to one of the three independent eigenvectors of , respectively.We can write   ,  = 3, 4, 5, as   ,  = 1, 2, 3, but for simplicity as   ,  = 1, 2, 3 in the rest of the paper, see Fig. 4(b) for the cup case of valence  = 7.
The quadratic precision error.As  2 continuity at extraordinary positions requires eigenfunctions corresponding to  to be standard quadratics of eigenfunctions for , the condition in Eq. ( 6) can be rewritten into the following form that can be used for optimization, where is the column vector of basis functions to parameterize the integral domain Ω in Fig. 4 following the method in [Stam 1998], ì  , ì  are eigenvectors corresponding to , ì   ,  = 1, 2, 3 are eigenvectors for , and   ,   ,   ∈ R for  = 1, 2, 3 are coefficients of the -th standard quadratics that can be obtained by minimizing respective integrals of   (, ) for  = 1, 2, 3 as, min If we consider (, ,   ) as a subdivision surface, the minimization in Eq. ( 12) gives standard quadratics of  and  that best fit the subdivision surface (, ,   ).The partial derivatives in Eq. ( 11) can be computed following the chain rule with the help of basis functions in ì  and invertible Jacobian matrices guaranteed by  1 conditions.
As the eigenstructure analysis is based on a 5-ring configuration, the integral in Eq. ( 12) should be performed over all 4-ring patches Ω, i.e., the shaded region in Fig. 4(a) with the hole enclosed.In our implementation, we can perform the integral only on the highlighted region Ω 0 in Fig. 4(a) which contains only regular patches for easy evaluation, and the remaining integral on the region enclosed by Ω 0 can be easily computed as the sum geometric sequences which are convergent due to the parameter bounds in Eq. ( 8).12): (a) the shaded base domain Ω 0 used for numerical quadrature and (b) an eigenfunction for the subsubdominant eigenvalue .Due to the scaling nature of the refinement process for eigenvectors, the integral over Ω, i.e., the union of Ω 0 and the hole region enclosed, can often be computed from the integral over Ω 0 as the sum of geometric sequences.
Weighting function to further improve second order characteristic maps.Curvature performances near extraordinary positions are heavily influenced by second order characteristic maps.If second order characteristic maps have sign-changing Gauss curvatures, the resulting limit surface might have hybrid curvature behavior that leads to artifacts [Peters and Reif 2004].
We use a weighting function to penalize sign changes for Gauss curvature on second order characteristic maps.Following the same notations in Eq. ( 11), Gauss curvature for three second order characteristic maps is computed as [Pressley 2010] As we are mostly interested in the sign of  *  (, ), we only use the numerator of  *  (, ) for subsequent computations, written as Ideally, if  2 continuity conditions are satisfied,   (, ) should be a constant equal to     −  2  .For the cup case of the second order characteristic map,     −  2  > 0, while for the saddle cases,     −  2  < 0. Now we redefine   (, ) as Through multiplication by the sign of the Gauss curvature for the ideal case, i.e., sign(    −  2  ), the new form of   (, ) in Eq. ( 15) will only be negative if Gauss curvature takes signs other than sign(    − 2  ).We thus only need to penalize the case for   (, ) < 0. A regularization operator for   (, ) can be defined as where   is the regularization parameter.With lower   values,  (  ,   ) approximates   (, ) better for   (, ) > 0, while for   (, ) < 0,  (  ,   ) becomes a small positive number that can be used in the denominator as a penalty for   (, ) < 0, see Fig. 5(a).
Based on the properties of  (  ,   ), a dimensionless weighting function is thus constructed of the following form where  is a parameter to control the maximum weight applied, and   describes the overall slope for   near   (, ) = 0.The maximum weight is  2 and lower   values produce steeper   curves, as illustrated in Fig. 5(b).A penalty will be applied if   (, ) < 0 when sign changes appear in Gauss curvatures, and the influence of the weighting function would be negligible for larger   (, ) values, which is desirable for stencil tuning.The weighted objective function for stencil optimization.In this work, we use the following integral of weighted quadratic precision error to quantify the quality of each map, where   (, ) and   (  ,   , ) are quadratic precision errors and the weighting function in Eq. ( 11) and Eq. ( 17), respectively, and the denominator is the area of the integral domain used for normalization.
The optimized stencil coefficients ( 2 ,  3 ) for specified  of valence  can be obtained through the following minimization problem, min where G( , ,  2 ,  3 ) contains parameter bounds in Eq. ( 8).
The minimization problem in Eq. ( 19) is solved by the 'fmincon' function in Matlab where the 'bfgs' algorithm is used for related Hessian matrix approximation.The objective function in Eq. ( 19) is a weighted version of that in [Ma and Ma 2019b] whose global minima have been verified.The introduction of the weighting function Eq. ( 17) slightly modifies the behavior of the objective function at lower  values, creating a local 'flat' region near the global minima that might lead to premature convergence.Stencils in this local 'flat' region exhibit almost identical performances according to our numerical tests.To guarantee a reliable optimal solution, we also perform optimizations for a number of initial solutions uniformly sampled from the feasible domain of ( 2 ,  3 ) (with a number of 10 × 10 samples) and use the locally optimized result that best fits the overall curves in Fig. 7.

SPLINE REPRESENTATION OF MAIN STENCIL COEFFICIENTS
Based on the tuning framework in the previous section, we perform subdivision tuning at a series of discrete  values, and further represent the optimized  2 and  3 coefficients as continuous B-spline functions of  through spline approximation.A continuous family of subdivision stencils can then be retrieved at any feasible  value through spline evaluation, which comes to the proposed unified -subdivision scheme.

Parameter setting for tuning at discrete 𝜆 values
4.1.1Domain for  sampling.As indicated in Fig. 3, 1/8 ≤  ≤ 1 are required for  1 continuity at extraordinary positions.We select samples for tuning subdivision stencils from  ∈ [0.15, 0.9], which would cover most scenarios for practical applications.Note that for consistency, we denote by [  ,   ] the domain of [0.15, 0.9] for  sampling in subsequent discussions, see Table 1.
4.1.2Strategies for  sampling.As bounded curvature is an important property for subdivision schemes, we divide the interval of  ∈ [  ,   ] at  =   , the lower bound for bounded curvature, and select 50 uniformly distributed sample  values in each subinterval, collectively denoted by , which should be sufficient for representing the entire family of unified subdivision scheme.
4.1.3Parameters   and  for the weighting function.There are two parameters for the weighting function in Eq. ( 17), namely   and .To ensure that the regularization have similar influences for all three second order characteristic maps, we use where  = 0.2, and   ,   ,   are coefficients in Eq. ( 15).For the parameter  in Eq. ( 17), we use  = 5 in this work.
4.1.4Simplified evaluation of weighted objective functions.To simplify the computation of Eq. ( 18) by making the most of the scaling property of eigenvectors, the weighting functions defined on inner spline rings are considered as scaled copies of those defined on Ω 0 in Fig. 4(a).The integrals in Eq. ( 18) can then be computed as sums of geometric sequences similar to that in [Ma and Ma 2018].Numerical verifications show that such a treatment produces reliable and consistent results for the optimization in Eq. ( 19).All integrals for computing the objective function in Eq. ( 19) are finally produced with 5 × 5 Gaussian points for each patch, and detailed procedures for the optimization are illustrated in Algorithm 1.

Calculation of global optimal stencil coefficients
Apart from stencil coefficients at discrete sample  values, we also compute global optimal subdivision stencils ( *  ,  2 ,  3 ) over the entire feasible domain for  for each valence  , which will be recovered by spline stencils in this work.The optimal  *  value for a specified  , along with corresponding optimal stencil  2 and  3 values, can be obtained by another layer of one-dimensional optimization with  being the variable for minimization by Eq. ( 19).The optimization framework in this case is slightly different from Algorithm 1, but similar to that in [Ma and Ma 2018] with a different objective function in this work.As expected for valence  = 4, the above optimization in  produces stencils for bi-cubic B-spline subdivision with  2 =  3 = 1/4 at  *  = 1/2, i.e., the regular case.

Spline representation of main stencil coefficients
With optimized stencil coefficients at sample  values {, f 2 , f 3 }, as shown in Fig. 6(a) for discrete points, including the global optimal stencil at  *  in Section 4.2, we further construct spline representations of  2 () and  3 () in  through B-spline approximation [Piegl and Tiller 1997]  Observed in Fig. 7 that there are feature  values at which the  2 and  3 functions are  0 continuous, mostly caused by parameter bounds.We place -multiple knots at such feature  values.
It should be noted that for consistency of knot structures, three -multiple knots, denoted by   ,   =   , and   , are assigned for all knot vectors of  = 3, . . ., 20, even though the functions for  2 and  3 are smooth at such  values for some valences  .In such cases, the required -multiple knots are inserted by knot insertion after obtaining the desired  2 or  3 functions.Between each pair of multiple knots, an intermediate knot is also inserted to enrich the knot vector for better representation of  2 and  3 functions.The knot vectors for  2 and  3 functions are structured in the form of where   = 17 is the number of control coefficients for  2 and  3 following the proposed knots structure.Details for each knot value can be found in Table 1.Knots with superscripts + and − are intermediate knots inserted at the center of neighboring multiple knots.
4.3.2Calculation of control coefficients for  2 and  3 functions.When constructing the spline stencils, we require resulting  2 () and  3 () functions to exactly recover the global optimal stencil at  *  for all valence  , including the regular stencil for  = 4 at  *  = 1/2.Given knots  with degree  defined in Eq. ( 21) and    optimized discrete stencil parameters of  2 and  3 with respective  values obtained for  ∈ [  ,   ] from the previous subsection, collectively denoted here as M, we thus further construct B-spline functions  2 () and  3 () approximating M with  through constrained B-spline approximation while interpolating the globally optimal stencil at  *  for all valences  .During the process of approximation, we also apply inequality constraints at all discrete sample stencils in M such that the resulting spline stencils of  2 () and  3 () would fall within respective upper and lower bounds in Eqs. ( 53)-( 62), see supplementary materials (Section B).During the approximation process in the implementation, we only apply multiple knots at visually  0 positions, if exist for certain valences  , and at the two end knots.We later insert respective multiple knots in the knot vector to recover the knot structure in  afterwards.Also, as indicated in Eqs. ( 61)-( 62) in supplementary materials (Section B), parameter bounds for  2 are influenced by  3 , so for each valence  , we produce the spline approximation of  3 first followed by the computation of  2 spline function.
The results of spline approximations are control coefficients of   (),  = 2, 3 for each of the respective valence  ≤ 20, that can be organized as a vector ì . Combined with knots  with degree , the B-spline representations of   (),  = 2, 3 are thus fully defined.Full lists of control coefficients ì     for   (),  = 2, 3 can be found in Table 5 and Table 6 in supplementary materials (Section C).See also Section 4.5 for further information on the evaluation of   (),  = 2, 3 and further construction of desired full subdivision stencils.

4.3.3
Illustration of spline functions of  2 () and  3 ().The obtained spline functions of  2 () and  3 () for  = 6 are illustrated in Fig. 6.From Fig. 6(b) and (d), it can been seen that spline functions  2 () and  3 () approximate discrete optimized  2 and  3 values well over the entire feasible domain for  ∈ [  ,   ].The relaxation parameter  in Eq. ( 10) follows the same pattern of Fig. 6(e) for all valences  .The angle  is the half-angle of the outermost corner of 1-ring quadrilaterals in the respective natural configurations, and  should be  ≤  2 to avoid concave corners and to guarantee regularity and injectivity of characteristic maps for  1 continuity.In this work, we require that  ≤  2 −  2 , leading to parameter bounds of Eq. ( 62) and Eq. ( 60) in supplementary materials (Section B).The application of  constraints in stencil optimization and for regularity check are similar to that in [Ma and Ma 2019a].Fig. 6(f) illustrates the resulting spline stencil parameters  2 (),  3 (), and respective  ( , ,  2 ,  3 ), with marked parameters at some featured  for valence  = 6.Fig. 7 provides further illustrations of  2 () and  3 () similar to that in Fig. 6(f) for some other selected valences.

Highlights of featured 𝜆 values and regions
Properties of the proposed unified -subdivision scheme depend upon the selection of the  value.Table 1 highlights several featured  values that are useful for practical applications.
Table 1.Summary of important and featured  values, see also Fig. 3. Multiple knots used in Eq. ( 21), i.e.,   ,   =   , and   , are shaded with light gray, while   = 0.15 and   = 0.9 are not listed in the •   is the lower bound for  to have bounded curvature.Values of   can be found in Table 1, and explicit computation for   can also be found in Eq. ( 23).

• 𝜆 *
is the value of  with global optimal bounded curvature for any valence  over the entire feasible domain for .By using  =  *  , the  2 and  3 functions will return subdivision stencils that produce the best possible surface qualities.
•  *  = 0.39 is the value of  that produces IGA solutions with optimal convergence rates in  2 -norms [Ma and Ma 2019a].If errors in other norms are considered, such as  1 -and  ∞norms, even lower  values should be used.
•  *  = 1/2 is the value of  that produces uniform refined meshes with the least polar artifacts.In summary, while the feasible region of the resulting spline stencil parameters of  2 () and  3 () is  ∈ [  ,   ], the recommended region for practical applications is  ∈ [  1 ,   2 ] for producing quality subdivisions.While the proposed unified  subdivision scheme produces bounded curvature subdivisions in theory for  ∈ [  ,   ] as shown in Fig. 3, it produces well behaved bounded curvature subdivisions for  ∈ [  ,   2 ].
Also note that at  *  ,  *  , and  *  , the proposed unified  subdivision scheme produces optimized subdivision stencils having similar properties to that in [Augsdörfer et al. 2006;Ma and Ma 2018] with the best curvature performance, [Ma and Ma 2019a;Wei et al. 2021] at lower  values with improved performance in IGA, and [Ma and Ma 2019b;Reif and Sabin 2019] at  close to 1/2 with uniform refined meshes, respectively.

Stencil evaluation at given 𝜆 values
As stencil coefficients  2 () and  3 () are represented as B-spline functions, they can be efficiently evaluated [Piegl and Tiller 1997].If the column vector ì  () contains B-spline basis functions corresponding to  at , then the  2 and  3 values at  for valence  , denoted by   ( , ),  = 2, 3, can be written as where ì     is the vector of control coefficients for   of valence  ,  = 2, 3, and ì     is available in supplementary materials (Section C).In Table 2 and Table 3 of the attached supplementary materials (Section C), we also include selected evaluations of  2 () and  3 () at some featured  values that can be readily used for further computing relevant full subdivision stencils.
With  , ,  2 , and  3 given, other stencil parameters or subdivision coefficients can be further computed by simple arithmetic formulae as follows.First we compute the lower bound   for bounded curvature.For a specific valence  ,   is given as where with  ,1 = cos(2/ ) and  ,2 = cos(4/ ).Note that values of   are also recorded in Table 1, with four decimal places.
Then we compute  following Eq.( 10), where   is used.The remaining stencil coefficients can thus be computed as where The intermediate variables,  1 ,  2 , and  3 , are generalized Catmull-Clark coefficients [Ma and Ma 2018] for refined edge vertices where  0 is the original extraordinary vertex,  0 and  1 are end vertices of the corresponding edge in the original mesh, and   and  ′  are two neighboring face vertices in the refined mesh.

EXAMPLES AND FURTHER DISCUSSIONS 5.1 Performance at various 𝜆 values
Fig. 8 provides illustrations of limit surfaces of two general models with subdivision in extraordinary regions at a range of  values using our method.One can observe that, the -subdivision scheme produces the most favorable models when  falls in the middle region, while minor visual surface twisting might appear when  is selected either too small or too large on the two sides.Fig. 9 provides some further illustrations of the proposed method on two other general models at some important featured  values,  *  ,  *  and  *  in Table 1 of Section 4.4, with also highlight lines display.The limit surfaces appear to be satisfactory for all three candidate  values while the highlight lines for  *  and  *  usually behave better than that for  *  .The result is consistent with observations in Fig. 8.In terms of refined mesh structure, the value of  directly influences the contraction of 1-ring vertices, which can be seen both in Fig. 10 with selected natural configurations and in Fig. 11(a-e) with refined meshes.The proposed -subdivision at  =  *  for each valence  in Table 1 produces the smoothest limit surface as shown in Fig. 9   .Two examples of the proposed method at various  values: (a-f) results for the mesh "tripleTorus" and (g-l) results for the mesh "brick".Both meshes in (a) and (g) are from [Cashman et al. 2009a].The mesh "brick" contains extraordinary vertices of valences  = 3, 6, 8, and for the mesh "tripleTorus", all the extraordinary vertices are of valence  = 6.
For detailed illustrations of performances of our method, we use the  = 7 mesh in Fig. 12 and plot refined meshes, surface renderings, Gauss curvature distributions, and highlight lines, in Fig. 11   ,  *  ] produces good quality limit surfaces and is recommended for surface representations, while small values of  <  *  would be preferred for IGA applications, such as at  *  mentioned in Section 4.4.

Comparisons with existing subdivision schemes
The proposed -subdivision scheme is an improved generalization of existing schemes in [Ma and Ma 2019a], [Ma and Ma 2018], and [Ma and Ma 2019b], with corresponding  denoted by   ,   , and   .Here we compare the proposed method with these schemes and the classical Catmull-Clark subdivision [Catmull and Clark 1978].
Comparison of second order characteristic maps: Compared with previous schemes [Ma and Ma 2019a], [Ma and Ma 2018], and [Ma and Ma 2019b], a major difference is that we use a weighted objective function in subdivision tuning to penalize the undesired sign changes of second order characteristic maps for producing better surfaces.We use   defined in Eq. ( 33) to measure the extent of overall sign change for three second order characteristic maps, 1 is the minimum Gauss curvature for the cup case which ideally should be positive, and   2 and   3 are the maximum Gauss curvature for two saddle cases which ideally should be negative.The smaller the   value, the better the respective subdivision with less Gauss curvature sign changes.
The comparison of   between the proposed scheme and [Ma and Ma 2019a] is illustrated in Fig. 13, which shows that the proposed scheme produces better maps with less Gauss curvature sign changes than that in [Ma and Ma 2019a].The results for comparison among the proposed scheme and three previous schemes [Catmull and Clark 1978], [Ma and Ma 2019b] and [Ma and Ma 2018] are omitted since   vanishes for all of them, which indicate that all these schemes, including ours at higher  values, produce second order maps with no Gauss curvature sign changes.
Further comparisons on meshes with a single EV: Our method is also compared with [Ma and Ma 2019a] on challenging single-EV meshes in Fig. 12 for further validation of performances at lower  values.Comparisons with Catmull-Clark subdivision are also included as a baseline, see Fig. 14 and Fig. 15.
At lower  values, e.g.,  =   , the proposed method produces tuned subdivisions with better second order maps and better limit surfaces than [Ma and Ma 2019a], as illustrated in Fig. 14(a-f) and Fig. 15(a-f).The ridges in surface renderings are alleviated, and the highlight lines are smoother for our method.Gauss curvatures by our method have tighter bounds and are more uniformly distributed.By using the same  value of Catmull-Clark subdivision, our scheme generates surfaces with smoother highlight lines and bounded Gauss curvature, see comparisons in Fig. 14(g-l) and Fig. 15(g-l).
For tuned subdivisions at  values in the middle region, as the influence of weighting function in Eq. ( 17) is limited, the proposed method produces almost identical results at  *  and  *  compared with tuned subdivisions in [Ma and Ma 2019b] and [Ma and Ma 2018], respectively, and relevant comparisons are thus omitted.

Other observations and further discussions
The proposed -subdivision scheme produces tuned subdivisions and further generalization of Catmull-Clark subdivision.All existing schemes for sharp features for Catmull-Clark subdivision can be directly integrated into our scheme.From Fig. 8(b)(h), we also observe that soft crease features appear on limit surfaces in extraordinary regions with lower  values.This property could be further addressed to develop an alternative scheme for defining soft or semi-sharp features over -subdivision surfaces.Properly organized subdivisions with a lower  value can be applied at tagged vertices/edges with which newly inserted vertices can be positioned close to the respective crease edge forming a soft crease.Sharpness of the soft crease can be controlled by . .Natural configurations at selected  values for  = 5, 8, 11 in columns 1 to 3, respectively.Note that respective maps are organized following an ascending order from the first row to the fifth row.Featured  values can be found in Table 1.
The proposed -subdivision scheme could also be used for mesh reparameterization with desired distribution of control vertices for producing field-aligned solutions in isogeometric analysis.For fluid flow simulation or crack animation, for examples, one can perform mesh refinement or develop mesh reparameterization algorithms using -subdivision rules to increase mesh density near crack or at other desired positions for producing solutions with reduced local errors and balanced solution error distribution.Fig. 13.Comparisons of   in Eq. ( 33) for measuring overall Gauss curvature sign changes of the three second order characteristic maps between the proposed scheme and MM2019a ( [Ma and Ma 2019a]).For comparison, we use the same  value as [Ma and Ma 2019a] for the proposed scheme.
potential users, and the evaluation of subdivision stencils at given  values can be performed by B-spline evaluations followed by simple arithmetic calculations.A tuning framework has been proposed for subdivision stencil optimization towards curvature continuity conditions.The bounded curvature constraint has been relaxed for extended feasible domains of  with minimum possible curvature The proposed unified -subdivision scheme can be applied in both graphics and engineering analysis.As future works, we will apply the proposed method for the creation of semi-sharp features, and use the method to solve relevant engineering problems.

Fig. 2 .
Fig. 2. Symbolized subdivision stencils near an extraordinary vertex: (a)-(c) stencil coefficients for refined face, edge, and vertex vertices, respectively.Hollow circles are regular vertices from the original mesh and red pentagram represents the extraordinary vertex in the original mesh.Red, blue and green solid circles are the refined face, edge, and vertex vertices, respectively.Note that all subdivision coefficients can be written as functions of subdominant eigenvalues , respectively, and details will be given in Section 4.

Fig. 4 .
Fig. 4. Illustration of the integral domain in Eq. (12): (a) the shaded base domain Ω 0 used for numerical quadrature and (b) an eigenfunction for the subsubdominant eigenvalue .Due to the scaling nature of the refinement process for eigenvectors, the integral over Ω, i.e., the union of Ω 0 and the hole region enclosed, can often be computed from the integral over Ω 0 as the sum of geometric sequences.

Fig. 5 .
Fig. 5. Illustrations for the weighting functions in Eq. (17): (a) the regularization function in Eq. (16) with different regularization parameters   ; and (b) the weighting functions in Eq. (17) with different  and   values.
4.3.1 Determination of knot vectors.We propose to use B-spline functions with degree  = 3 and with open knots for representing  2 and  3 against  for all valences  = 3, . . ., 20.

Fig. 7 .
Fig. 7. Stencil coefficients  2 () and  3 () for selected valences  ≤ 20.Regions shaded with orange and blue are the feasible domains for  2 and  3 , respectively.The recommended regions for  ∈ [  1 ,   2 ] in Section 4.4 are highlighted by a piece of thicker lines in the plots for  2 and  3 , and their boundaries are highlighted by black dashed lines.Some other featured  values in Section 4.4 are also marked for reference.
Fig.8.Two examples of the proposed method at various  values: (a-f) results for the mesh "tripleTorus" and (g-l) results for the mesh "brick".Both meshes in (a) and (g) are from[Cashman et al. 2009a].The mesh "brick" contains extraordinary vertices of valences  = 3, 6, 8, and for the mesh "tripleTorus", all the extraordinary vertices are of valence  = 6.
again at featured  values  *  ,  *  and  *  in Section 4.4.Surface renderings have consistent performances to those in Fig. 8, and the best highlight lines and Gauss curvature performances are also observed at  =  *  .Deviations from  =  *  might lead to possible degradation of surface qualities, in both highlight lines and curvature distributions.So when selecting  values, there is usually a balance between different priorities.Empirically,  ∈ [ * Fig. 9. Results from the proposed scheme on general meshes at selected : (a-b) initial control meshes of "Kitten" (from [Jakob et al. 2015] with EVs of  = 3, 5, 6) and "Propeller" (from [Scott et al. 2013] with EVs of  = 3, 5); (c-h) results for "Kitten"; and (i-n) results for "Propeller".Results for  =  *  ,  *  and  *  are illustrated in columns 1 to 3, respectively.The highlighted regions illustrate highlight lines near EVs of valence  = 5.
Fig.10.Natural configurations at selected  values for  = 5, 8, 11 in columns 1 to 3, respectively.Note that respective maps are organized following an ascending order from the first row to the fifth row.Featured  values can be found in Table1.
Fig. 11.Illustration of results for the  = 7 mesh at featured  values: (a-e) meshes after two steps of subdivision; (f-j) rendering of limit surfaces; (k-o) Gauss curvature distributions; and (p-t) highlight lines.Results for  = 0.35,  *  ,  *  ,  *  , and  *  + Δ with Δ = (  2 −  *  )/2, are illustrated in columns 1 to 5, respectively.Gauss curvature is denoted by   , and the same color scale is used for all Gauss curvature plottings.
Fig. 14.Comparison using the  = 5 mesh in Fig. 12: (a-f) comparison with MM2019a [Ma and Ma 2019a], and (g-l) comparison with CC1978 [Catmull and Clark 1978].The values   and   are the subdominant eigenvalue for MM2019a and CC1978, respectively.Gauss curvature is denoted by   , and the same color scale is used for all Gauss curvature plottings.
Fig. 3. Illustration of lower bounds   of  for different valence  for producing subdivisions with bounded curvature in  ≥   together with other featured  values inTable 1 (see later Section 4.4).In this figure,   1 is the lower bound of  to produce  1 limit surfaces.
In this figure, the superscript "#" denotes results from direct tuning through Eq. (19) at discrete sample  values, some featured  values are marked with vertical dashed lines, and recommended regions in Section 4.4 are shaded with light green in (a)-(e).
table.The column   contains  values for Catmull-Clark subdivision for reference.
•   1 is an empirical lower bound for the recommended region of .For  <   1 , the subdivision schemes are less stable.A slight decrease in  might lead to rapid increase of  ( , ,  2 ,  3 ) in Eq. (19), especially for higher valences, see Fig.7(g) and (h).•   2 is an empirical upper bound for the recommended region of .For  >   2 , the subdivision schemes are also less stable with possible overlapping or concave 2-ring quadrilaterals in respective natural configurations.