A general class of $C^1$ smooth rational splines: Application to construction of exact ellipses and ellipsoids

In this paper, we describe a general class of $C^1$ smooth rational splines that enables, in particular, exact descriptions of ellipses and ellipsoids - some of the most important primitives for CAD and CAE. The univariate rational splines are assembled by transforming multiple sets of NURBS basis functions via so-called design-through-analysis compatible extraction matrices; different sets of NURBS are allowed to have different polynomial degrees and weight functions. Tensor products of the univariate splines yield multivariate splines. In the bivariate setting, we describe how similar design-through-analysis compatible transformations of the tensor-product splines enable the construction of smooth surfaces containing one or two polar singularities. The material is self-contained, and is presented such that all tools can be easily implemented by CAD or CAE practitioners within existing software that support NURBS. To this end, we explicitly present the matrices (a) that describe our splines in terms of NURBS, and (b) that help refine the splines by performing (local) degree elevation and knot insertion. Finally, all $C^1$ spline constructions yield spline basis functions that are locally supported and form a convex partition of unity.


Introduction
Multivariate splines are used extensively for computeraided design (CAD) and, more recently, for computer-aided engineering (CAE). Smoothness of such splines is a particularly valuable trait. When the aim is to create a (freeform) geometric model for a smooth object, it helps if the splines used for the task are smooth themselves. For instance, this circumvents situations where small displacements to control points may produce 'non-smooth' features such as C 0 kinks, loss of curvature continuity, etc. Similarly, when the aim is to numerically approximate the solution to high-order partial differential equations (PDEs) using isogeometric analysis (IGA) -a generalization of classical finite elements [8] high smoothness of the approximating spaces can be beneficial. For instance, it can allow us to directly discretize the PDEs without any auxiliary variables, thus yielding simpler and more efficient implementations.
In this paper, we discuss a general class of C 1 smooth rational splines that allow for the construction of C 1 smooth curves and surfaces. These are an extension of classical C 1 non-uniform rational B-splines (NURBS) as they enjoy the flexibility of choosing locally unrelated weight functions as well as the option of local degree elevation -they can be roughly regarded as piecewise-NURBS. At the same time, they maintain intuitive control-point-based design. Moreover, they enable simple (low-degree) and smooth descriptions of some of the most important primitives for CAD and CAE (but also for computer vision, graphics and robotics): closed, real, non-degenerate quadrics -that is, ellipses in two dimensions and ellipsoids in three dimensions.
The ideas we present here build upon those from [25], in multiple directions, and their presentation is motivated by our primary objectives: self-contained, explicit, NURBScompatible descriptions that can be easily and efficiently implemented within existing CAD software. The most important novel contributions are the following.
• We describe the usage of classical univariate NURBS to assemble C 1 rational multi-degree spline basis functions using an extraction matrix. The general framework was explained in [25], but we provide here a simplified exposition of the construction and a formal proof of the properties; see Remark 2.3. We mainly stick to parametric smoothness, but a construction centred around the notion of geometric smoothness can be formulated as well; see Remark 2.5.
• We describe efficient refinement of the C 1 splines leveraging classical NURBS refinement. The novelty here relies in an explicit and simple construction of the refinement matrices.
• We describe how tensor-product bivariate C 1 rational splines can be used to build C 1 smooth geometries that may contain one or two polar singularities; the C 1 smooth splines describing the geometries are called polar splines. As above, the idea is based on building an extraction matrix.
• We describe efficient refinement of polar splines. In particular, we provide an explicit and simple construction of the refinement matrices.
• We provide explicit descriptions of ellipses and ellipsoids built using low-degree C 1 splines, and we detail their extraction in terms of NURBS so that they can be readily implemented and used in CAD or CAE software. Table 1 summarizes the descriptions included in this paper.

Extraction matrices
At the core of our approach is the notion of the socalled design-through-analysis (DTA) compatible extraction matrix. 1 Roughly speaking, such matrix helps us assemble 'simple splines' into 'more general splines.' Examples are the Bézier extraction matrix introduced to assemble Bernstein polynomials into B-/T-splines [3,19]; the multidegree extraction matrix for assembling elements of extended Tchebycheff spaces into generalized Tchebycheffian B-splines [22,26,7]; and the unstructured spline extraction matrices for assembling tensor-product splines into splines on unstructured quadrilateral meshes [25,27,24].
Here, we apply the concept of extraction in the following context. We start from multiple sets of (univariate or bivariate) NURBS basis functions defined on adjacent domains, and collect all of these functions in the set {b j : j = 1, . . . , m}. Then, we assemble them into more general C 1 rational (polar) splines using a matrix C (with entries C i j ), called the extraction matrix. Denote this new set of splines by {N i : i = 1, . . . , n}, where n < m. These are defined as follows, We are particularly interested in matrices C such that the functions N i • satisfy certain smoothness constraints that may or may not be satisfied by the b j , and • possess the properties of non-negativity, locality, linear independence and partition of unity that the b j already possess.
Such extraction matrices are called DTA-compatible. It is easy to see that the action of a DTA-compatible extraction matrix on a convex partition of unity, local basis gives rise to another local basis that also forms a convex partition of unity. Indeed, by summing over i in Equation (1), we have as the b j form a partition of unity. Since C has non-negative entries and is a full-rank matrix, non-negativity and linear independence of N i follow from the non-negativity and linear independence of b j .

Related literature
As mentioned in the previous section, the construction of smooth univariate splines by joining simpler pieces has been recently explored in [25,22,26] for polynomial multidegree splines, and in [7] for generalized Tchebycheffian splines. These approaches have conceptual similarities with the notion of beta-splines [2]. The main differences are that the former approaches do not rely on symbolic computations while the latter does, and the former approaches consider parametric continuity while the latter studies geometric continuity. The use of smooth univariate rational splines for construction of circles has been previously explored in [1,12,13]. It is known that a circle cannot be represented by a single (symmetric) periodic C 1 quadratic NURBS curve [17, Section 7.5] nor a C 2 cubic NURBS curve [5,Section 13.7]. However, it is possible to find C p smooth descriptions using NURBS of degree 2(p+1), which is shown to be the minimal degree in [1]. On the other hand, [12,13] presented a C 1 piecewise quadratic NURBS description of the circle and used it for IGA. Our rational multi-degree splines form a flexible extension of the latter framework, and allow for a variety of exact descriptions of circles using low (multi-)degrees, as indicated in Table 1.
In two dimensions, closed quadrics or, more generally, smooth closed surfaces of genus zero can be built using tensor-product splines by introducing polar singularities. For such polar surfaces, subdivision schemes producing C 1 surfaces [10,14] and C 2 surfaces [9,15] have been previously worked out. The corresponding limit surfaces consist of an infinite sequence of surface rings where the faces shrink to a point in the limit. A more CAD-friendly finite construction was developed in [16]; this approach constructs 'shape' basis functions for C 2 polar splines with bi-degree (6, 3). These basis functions correspond to unique Fourier frequencies in the polar expansion of a quadratic surface. The 'shape' basis does not enjoy non-negativity and does not form a partition of unity, and extensions of it to higher smoothness leads to degrees of freedom that control non-intuitive shape parameters. Similar recent constructions for obtaining C 1 polar spline caps can be found in [11]. Curvature continuous polar NURBS surfaces were discussed in [21], and [20] presented a construction of polar caps using periodic B-spline surfaces with G n continuity for arbitrary n. On the CAE side, a standard circular serendipity-type element for IGA was proposed in [12], and C k smooth basis functions over singular parametrizations of triangular domains were constructed in [23]. A design-through-analysis friendly construction of C k smooth polar surfaces was recently proposed in [25], and the current work builds further upon this construction.
A completely different approach for dealing with curves and surfaces is the use of implicit representations [6]. Such representations enjoy nice geometric properties (especially for simple shapes). For instance, they allow for a straightforward point membership classification. On the other hand, explicit smooth B-spline representations are more convenient for direct geometric modeling and (local) modification.

Outline
In Section 2, we present the construction and refinement of C 1 smooth rational multi-degree spline curves via explicitly defined extraction and refinement matrices. The construction and refinement of C 1 smooth polar surfaces using

Piecewise-rational curves
In this section, we focus on a multi-degree extension of univariate NURBS splines. The multi-degree spline space is defined as a collection of classical NURBS spaces (with possibly different polynomial degrees and weight functions) glued together C 1 smoothly. For such space we present a construction of a set of basis functions, with similar properties to classical NURBS. After discussing some preliminary material on NURBS in Section 2.1, we elaborate how these basis functions can be computed through a DTA-compatible extraction matrix in Section 2.2. A more general but also more complex algorithmic construction has been detailed in [25,Section 2] and further explored in [22,26] for polynomial multi-degree splines. Then, in Section 2.3, we give an explicit procedure how to compute a refined representation of a given curve. Finally, in Section 2.4, we illustrate how this tool can be used to describe arbitrary ellipses in a C 1 smooth fashion using low-degree piecewise-rational curve representations suited for integrated design and analysis.

Preliminaries on NURBS
We start by defining notation for NURBS basis functions, and introduce some classical relations that can be found, e.g., in [18,17].
Given a basic interval I := [x 1 , x 2 ] ⊂ , let us denote with ξ an open knot vector of degree p ∈ and length n+p+1 ∈ , i.e., ξ := [ξ 1 , ξ 2 , . . . , ξ n+p+1 ], ξ i+1 ≥ ξ i , The number of times a knot value ξ i is duplicated in the knot vector is called the knot's multiplicity. The multiplicity of ξ i is denoted with m i , and we assume that 1 ≤ m i ≤ p − 1.
The corresponding set of B-splines {b j,p : j = 1, . . . , n} are defined using the recursive relation, and under the convention that fractions with zero denominator have value zero. With the above definition, all the B-splines take the value zero at the end point x 2 . Therefore, in order to avoid asymmetry over the interval I, it is common to assume the B-splines to be left continuous at x 2 . We will follow suit. Let us denote with w a weight vector of length n, i.e., w := [w 1 , w 2 , . . . , w n ], w i > 0.
The corresponding set of NURBS {b w j,p : j = 1, . . . , n} are defined by .
Each b w j,p is non-negative on I and is locally supported on [ξ j , ξ j+p+1 ]. Moreover, the functions b w j,p are linearly independent and form a partition of unity. They satisfy the following end-point conditions: The NURBS space corresponding to ξ and w is denoted with [ξ, w ] and is defined as the span of {b w j,p : j = 1, . . . , n}. This is a space of piecewise-rational functions of degree p with smoothness C p−m i at knot ξ i and its dimension is n. The assumption on the multiplicity will ensure us global C 1 smoothness. Note that when w 1 = · · · = w n , the members of this space are piecewise-polynomial.
Remark 2.1. The structure of ξ in Equation (2) is such that p, m i , n and I are embedded in it. Therefore, we will assume that once a knot vector ξ is known, so are the degree, smoothness, and dimension of the corresponding NURBS space [ξ, w ].
We identify a function f ∈ [ξ, w ] with the vector of its coefficients [ f 1 , . . . , f n ], Only the first (last) k + 1 basis functions contribute towards the k-th order derivative at the left (right) end point of I. In particular, we have A NURBS curve embedded in d , d ≥ 2, can be constructed as where f j ∈ d are the control points assigned to each basis function. All coordinate functions of this curve belong to [ξ, w ] and therefore all the above relations hold for them.

Rational multi-degree B-splines
Consider m open knot vectors ξ (i) of degree p (i) , i = 1, . . . , m, defined as in Equation (2). We denote the left and right end points of the interval I (i) associated to ξ (i) with x (i) 1 and x (i) 2 , respectively. The collection Ξ := (ξ (1) , . . . , ξ (m) ) is called an m-segment knot vector configuration. The multidegree spline spaces will be constructed by considering spline spaces over the knot vectors ξ (i) , which are glued together with certain smoothness requirements at the end points x is called the i-th segment join. We define the mapping φ (i) for each segment i = 1, . . . , m, for an arbitrarily chosen origin τ ) ⊂ , and we construct the composed interval (1) , . . . , w (m) ) be a sequence of weight vectors defined as in Equation (3). We refer the reader to Figure 1 for a visual illustration of the notation of the above concepts, in case m = 2 and p (1) = 2, p (2) = 3.
The space of rational multi-degree splines is defined as and the periodic space of rational multi-degree splines as Figure 1: A visual illustration of the notation and the construction of rational multi-degree B-splines as described in Section 2.2. Here, the quadratic (blue) and cubic (red) NURBS shown at the bottom are used to build the C 1 multi-degree B-splines shown at the top.
The elements of [Ξ, W ] and per [Ξ, W ] are piecewise-NURBS functions such that the pieces meet with C 1 continuity at each segment join. It is clear that classical NURBS spaces are a special case of the rational multi-degree spline spaces.
In the following, we build a suitable basis for the spaces In the first step, we map these basis functions from I (i) to Ω (i) using φ (i) in Equation (5), and extend them on the entire interval Ω by defining them to be zero outside Ω (i) . More precisely, specifying the cumulative local dimensions µ i for i = 0, . . . , m, we define for i = 1, . . . , m and j = 1, . . . , n (i) , For the sake of simplicity, we dropped the reference to the (local) degree and weight in the notation. From the properties of NURBS, it is clear that the functions b 1 , . . . , b µ m are linearly independent and form a non-negative partition of unity on Ω. We arrange these basis functions in a column vector b of length µ m . We refer the reader again to Figure 1 for a visual illustration of the notation of the above concepts. Now, we construct extraction matrices H and H per such that the functions in {B i : i = 1, . . . , n} and {B per i span [Ξ, W ] and per [Ξ, W ], respectively. The key here, and the reason our approach can be efficiently implemented by design, is that these extraction matrices can be explicitly specified. To this end, we define counters η i for i = 0, . . . , m, and parameters α (i) and β (i) for i = 1, . . . , m − 1, In the periodic setting, α (m) and β (m) are computed using the above equations by identifying the index i + 1 with 1. Recall Equation (4) to see the motivation behind the definition of the above parameters. Then, we define a common sparse matrix H c of size η m × (µ m − 2), whose non-zero entries H c i j are identified as follows: for i = 1, . . . , m and j = 1, . . . , and for i = 1, . . . , m − 1, The desired extraction matrices in Equation (6) are then specified as follows: The number of rows in the two matrices are denoted with n := η m + 2 and n per := η m , respectively. The sparse and simple structure of both matrices means that it is easy to verify that both have full rank. Indeed, this conclusion can be directly deduced from the full rank of H c , which in turn is implied by Equation (7). Moreover, their entries are nonnegative, and the column sum is equal to one. Hence, we conclude that these matrices are DTA-compatible.
How these matrices help us build C 1 splines can be understood by taking into account Equation (4) More precisely, only these four functions have non-vanishing values and first derivatives here. In view of (4), a spline f given by We can verify that the entries of H satisfy exactly such relations. Indeed, for some j, the matrix H defines two new functions B j and B j+1 such that When setting [ f 1 , f 2 , f 3 , f 4 ] equal to the first or the second row ofH (i) , we see that Equation (10) is satisfied. The following result follows from the above discussion. Remark 2.3. In [25, Section 2.3.5] it was observed that the C 1 smooth piecewise-rational basis functions enjoy the properties described in Theorem 2.2. However, a formal proof was missing. It was also pointed out that the property of nonnegativity is in general not present in case of C 2 or higher smoothness. On the other hand, this is possible when restricting to polynomial pieces [26].
Remark 2.4. With the aim of designing quadric curves, also called conics, it is natural to choose local NURBS spaces of the same degree p and defined on the same uniform knot vector ξ. Moreover, it is common to set w Under these circumstances, the ratios in Equation (8) read as Finally, if there is additional symmetry in the choice of weights, so w (i) , we simply get Once we have computed a DTA-compatible extraction matrix H (or H per ), given n control points f i ∈ d , d ≥ 2, we can construct a piecewise-rational curve f embedded in d , For a fixed curve, the transpose of H (or H per ) defines the relationship between control points of the b j (discontinuous at the segment joins) and control points of the smooth B i . More precisely, if Remark 2.5. When dealing with curves, the proposed piecewise-NURBS framework can also be formulated in the context of geometric continuity [4]. In such case, the C 1 smoothness condition at the segment join in (10) is replaced by the G 1 smoothness condition for a given geometric shape parameter γ (i) > 0, resulting in the matrixH It is clear that this matrix is still DTA-compatible.

Refinement of piecewise-rational curves
The rational spline spaces defined in the previous section can be refined in a multitude of ways. We could reduce the smoothness at segment joins, raise the polynomial degrees of local NURBS spaces, and/or insert new knots in local NURBS spaces [25,Section 2.4.3]. A combination of these possibilities could be judiciously employed to achieve spline spaces that provide higher resolution or approximation power exactly where needed. In this section, we present an explicit construction of refined representations of a given piecewiserational curve.
Before delving into the details of the refinement procedure, we first define two matrices G and G per that can be regarded as right inverses of the extraction matrices H and H per , respectively. Looking at the structure of the matrix H c specified in Equations (7)-(8), we can define a sparse matrix G c of size (µ m − 2) × η m , whose non-zero entries G c i j are identified as follows: for i = 1, . . . , m and j = 1, . . . , n (i) − 2, From its construction it is clear that the product H c G c is equal to the identity matrix. Similarly, keeping in mind Equation (9), the matrices give rise to products HG and H per G per that are equal to identity matrices. Now, let be a given spline space and let us denote the target refined space with˜ . For simplicity of notation, we drop the superscript per in case of periodicity. Then, we consider the two unique representations of a curve f with coordinate functions in ⊂˜ , . We now seek the refinement matrix R of size n ×ñ that helps us computeF from F , i.e.,F = F R.
Assume that H andH are the extraction matrices corresponding to the spaces and˜ , respectively. Incorporating these matrices in the representations in Equation (13)  This implies that we can compute R by solving the following (overdetermined) linear system with a unique solution,

RH = HS.
After multiplication of both sides of this system with the ma-trixG (corresponding toH) as defined in Equation (12), we arrive at R = HSG. (14) Note that the application ofG in Equation (14) means that a subset ofñ columns of HS are selected to form R.
Remark 2.6. The definition ofG is done for the sake of simplicity of computation of R in Equation (14), but is not unique. Any matrix that is a right inverse ofH would be a valid choice as well, such as the standard Moore-Penrose right inverseH T (HH T ) −1 .

Circles and ellipses
We now present the general construction of ellipses (and as a special case also circles) using the C 1 rational splines introduced thus far. We present three approaches for doing so using splines of low(est) degree, i.e., C 1 splines of quadratic degree, cubic degree and mixed quadratic/cubic multidegree. All approaches will construct four C 1 piecewise-NURBS functions B i and associated control points f i , i = 1, . . . , 4, such that the curve f , describes the exact ellipse centred at (0, 0) and with axis lengths (a x , a y ), Since the splines B i form a partition of unity, these ellipses can be affinely transformed by directly applying the transformation to the control points f i . Subdivided or higher-degree representations can be easily obtained by refining the representations provided here (see Section 2.3).
To visually illustrate the smoothness of f , we will also show the curvef obtained by perturbing one of the control points. Since all B i are smooth, the perturbed curve will also be smooth. For uniformity throughout the examples, we will choose the control points of the perturbed curve as

C 1 description of degree 2
Here we present a C 1 quadratic description of the ellipse in Equation (15) Figure 2 (a, top row). Finally, we can build an ellipse centred at (0, 0) and with axis lengths (a x , a y ) by combining the splines B i with the control points f i defined as Remark 2.7. To verify that the curve f satisfies Equation (15), we can proceed as follows. The simplest approach is to numerically evaluate f (t) at all t and plug the result in that equation. Alternatively, this verification can also be performed analytically by looking at the explicit expressions of the rational pieces that form f . For instance, consider the first quadratic rational piece, g (1) , that is a part of f . As discussed in Equation (11), we can get the control points of this piece, denoted with g (1) j , j ∈ {1, 2, 3}, by applying the transpose of (a submatrix of) H per from Equation (17) to a vector containing the points f i , i.e., g (1) This yields the control points 3 = (a x , 0). Combining the above control points with the NURBS basis defined on the first segment, where w(t) := (1− t) 2 + 2t(1− t)+ t 2 , some simple algebra shows that g (1) indeed satisfies Equation (15). Verifications for the other pieces of f can be similarly done.
Example 2.8. Choosing a x = a y = 1, we obtain a circle of radius 1, as shown in Figure 2 (a, middle row). This C 1 quadratic description is equivalent to the one used in [12].
The choice a x = 2a y = 1 yields an ellipse with axis lengths (1, 1 2 ), as shown in Figure 2 (a, bottom row). The perturbed versions of these conics, with the control points chosen as in Equation (16), are shown as well and they remain clearly smooth.

C 1 description of degree 3
Here we present a C 1 cubic representation of the ellipse in Equation (15) These basis functions are shown in Figure 2 (b, top row). Choosing the associated control points f i as we get a C 1 cubic description of an ellipse centred at (0, 0) with axis lengths (a x , a y ). This can be verified in the vein of Remark 2.7.
Example 2.9. Choosing a x = a y = 1, we obtain a circle of radius 1, as shown in Figure 2 (b, middle row). This C 1 cubic description is equivalent to the one used in [25]. The choice a x = 2a y = 1 yields an ellipse with axis lengths (1, 1 2 ), as shown in Figure 2 (b, bottom row). The perturbed versions of these conics, with the control points chosen as in Equation (16), are shown as well and they remain clearly smooth. It can be observed that, compared to the description from Section 2.4.1, the control points here are at a greater distance from the curve. This is completely analogous to the behavior of classical NURBS. Then, we can define four C 1 multi-degree piecewise-NURBS functions B i on Ω using the extraction matrix Choosing the associated control points f i as we get a C 1 multi-degree description of an ellipse centred at (0, 0) with axis lengths (a x , a y ). This can be verified in the vein of Remark 2.7.
Example 2.11. Choosing a x = a y = 1, we obtain a circle of radius 1, as shown in Figure 2 (c, middle row). The choice a x = 2a y = 1 yields an ellipse with axis lengths (1, 1 2 ), as shown in Figure 2 (c, bottom row). The perturbed versions of these conics, with the control points chosen as in Equation (16), are shown as well and they remain clearly smooth. Once again, compared to the description from Section 2.4.1, the control points here lie at a greater distance from the cubic portion of the curve.

Piecewise-rational polar surfaces
In this section, we describe how to construct C 1 smooth representations for polar surfaces containing single or double polar singularities (e.g., hemispheres and spheres, respectively). Such surfaces can be obtained by starting from a bivariate tensor-product (piecewise-NURBS) spline patch and collapsing one or two of its edges, respectively, as illustrated in Figure 3. Each of such edge collapses creates a polar point and can be achieved by coalescing the control points related to basis functions with non-zero values on the edge. In general, however, this control-point coalescing will introduce kinks at the poles and the surface representation will not be smooth. To achieve overall smoothness, additional conditions need to be satisfied by the control points [25, Section 3]. In Section 3.1, we derive C 1 smoothness conditions at a polar point, and they enable us to build smooth polar splines as linear combinations of bivariate tensor-product splines in Section 3.2. Then, in Section 3.3, we give an explicit procedure how to compute a refined representation of a given polar surface. Finally, in Section 3.4, we present explicit descriptions of arbitrary ellipsoids using C 1 smooth low-degree polar spline representations suited for integrated design and analysis.

Smoothness conditions at the polar points
A polar surface will be smooth at a polar point if it can be locally (re)parameterized in a smooth way. Such parameterizations can be specified in a constructive manner and we elaborate upon it in this section. The resulting conditions will help us build smooth polar B-splines in the next section.
As shown in Figure 3, we first describe the initial setupa tensor-product spline space on a rectangular domain. We start from two univariate C 1 rational spline spaces s , t defined on the univariate domains Ω s := [s 1 , s 2 ] and Ω t := [t 1 , t 2 ], respectively; the superscripts of s and t are meant to indicate the symbols used for the respective coordinates. Using a Cartesian product, we build the rectangular domain Ω := Ω s × Ω t , and on Ω we define the tensor-product spline space := s ⊗ t . Without loss of generality, we assume that s 1 = t 1 = 0. This tensor-product spline space is spanned by tensor-product B-spline basis functions B i j , i = 1, . . . , n s ; j = 1, . . . , n t . Here, n s and n t denote the respective dimensions of the chosen univariate spline spaces; the basis functions spanning these spaces are denoted with B s i and B t j . Then, the tensor-product basis function B i j is simply the product B s i B t j . The functions B i j are assumed to be periodic in s and non-periodic in t. Now, let us use the functions B i j to map the domain Ω to a polar surface using edge-collapse. Then, the smoothness conditions at a collapsed edge will only involve those B i j that have non-zero first derivatives there. Observe that, if n t ≥ 4, then any B i j with non-zero first derivatives at the bottom edge of Ω will have zero first derivatives at its top edge, and vice versa. The upshot is that, when we are collapsing both the bottom and top edges of Ω into two polar points, as in Figure 3 (b), the smoothness conditions at those points are independent of each other and can be resolved separately for n t ≥ 4.

Single polar point
In light of the above discussion, in the following we first focus on the case of a single collapsed edge, i.e., the one shown in Figure 3 (a). We derive smoothness conditions that will help us build smooth polar spline functions (and, using them, smooth polar surfaces). This is done by explicitly specifying the parameterization with respect to which the spline functions are deemed smooth. First, we construct a planar disklike domain Ω pol , called the polar parametric domain, via a suitable polar map F ; see Figure 3 (a). Next, for an arbitrary In general, f pol will be multivalued at the pole. Finally, we derive the required smoothness conditions by asking for f pol to be C 1 smooth at the polar point.
We start by building F . Assign the control point F i j := (ρ j cos(θ i ), ρ j sin(θ i )) ∈ 2 to the basis function B i j , where and The above choice of control-point values has been made in the interest of standardization and is not unique. Using these control points, we can construct the disk-like domain Ω pol with the aid of the map F from Ω to Ω pol , t). (21) Note that the above construction will not necessarily yield an exactly circular domain Ω pol ; its shape will depend on the choice of . This domain will serve as the reference element for the smoothness of polar configurations, i.e., we will define polar splines such that they are C 1 smooth functions over Ω pol . It is clear that for all s ∈ Ω s , where (0, 0) ∈ Ω pol is the polar point. Note that this implies Let B pol i j be the image of B i j under the polar map F : Ω → Ω pol in Equation (21) so that Then, for given coefficients f i j , a polar spline function f pol over Ω pol can be constructed as We can pull f pol back to Ω as follows, Moreover, by using the chain rule we can also relate the partial derivatives of f and f pol : For f pol to be C 1 smooth at the polar point, there must exist real values α, β, γ such that In view of (22) and (23) this means for all s ∈ Ω s , In particular, since only B i j , j ≤ 2, have non-zero values and derivatives when t = 0, the above condition translates to the following requirement for all s ∈ Ω s ,

Double polar point
Equation (26) shows the required smoothness conditions when the bottom edge of Ω is being collapsed. Next, if we also want to collapse the top edge of Ω, we can repeat the previous argument with minor changes. We would, of course, need to choose a mapF that collapses the edge Ω s × {t 2 } instead. One way of achieving this could be by choosing the control pointsF i j := (ρ j cos(θ i ),ρ j sin(θ i )) ∈ 2 , whereρ j := 1 − ρ j ,θ i := 2π − θ i . Then, we can follow the same argument as in Section 3.1.1.
Asking for C 1 smoothness off pol is equivalent to asking that there exist real valuesα,β,γ such that for all s ∈ Ω s , Note once again that the smoothness at the polar point corresponding to t = 0 is imposed with respect to the parameterization F (Ω), while that at the polar point corresponding to t = t 2 is imposed with respect to the parameterizationF (Ω). The corresponding smoothness conditions in Equations (26) and (27) involve different coefficients f i j for n t ≥ 4, and so can be resolved separately. Remark 3.1. The choices of θ i , ρ j ,θ i ,ρ j are such that the maps F andF preserve the orientation of the parametric domain Ω.

Rational polar B-splines at the polar points
We now elaborate how the derived C 1 smoothness constraints at a polar point will enable the computation of a DTA-compatible extraction matrix. This matrix represents a linear map to a set of polar spline basis functions that are C 1 smooth on the polar parametric domain.

Single polar point
As before, we start by considering the case of a single collapsed edge, i.e., the one shown in Figure 3 are C 1 at the polar point. For fixed j, the set {B pol i j : i = 1, . . . , n s } is called the ( j − 1)-th polar ring of basis functions. When j > 2, all basis functions in the ( j − 1)-th ring already satisfy the C 1 continuity conditions at the polar point (their derivatives are identically zero there), so they can be included without modifications in the set of polar basis functions being created. The others will be substituted by three smooth polar basis functions. This dictates that E will be a matrix, with n := n s (n t − 2) + 3 rows and n s n t columns, taking the following sparse block-diagonal form: where I is the identity matrix of size n s (n t −2)×n s (n t −2) and E is a matrix of size 3 × 2n s . The entry ofĒ corresponding to its l-th row and (i + ( j − 1)n s )-th column is denoted with E l,(i j) . We can then rewrite Equation (28) as follows for l = 1, 2, 3, We can pull these back to Ω using Equation (24) to obtain the equivalent representation for l = 1, 2, 3, We will enforce C 1 continuity at the polar point by requiring the basis functions N pol l to satisfy a linearly independent Hermite data set at the polar point, in the spirit of Equation (25). To this end, we will use three source basis functions {T l : l = 1, 2, 3}, that provide us with the appropriate Hermite data. Given a non-degenerate triangle with vertices v 1 , v 2 and v 3 , let (λ 1 , λ 2 , λ 3 ) be the unique barycentric coordinates of point (u, v) with respect to such that Then, we define T l (u, v) := λ l , l = 1, 2, 3.
These functions can be interpreted as triangular Bernstein polynomials of degree 1. They are non-negative on the domain triangle . Moreover, they are linearly independent, form a partition of unity, and span the space of bivariate polynomials of total degree less than or equal to 1. Then, we require that N l in Equation (30) is a spline function f such that it satisfies the continuity constraints in Equation (26), with for l = 1, 2, 3. In the interest of standardization, we choose the triangle as equilateral with vertices recall the definition of ρ 2 from Equation (19). After some calculations, we deduce that This relation says that Ē 1,(i2) ,Ē 2,(i2) ,Ē 3,(i2) are simply the barycentric coordinates of the control point F i2 := (ρ 2 cos(θ i ), ρ 2 sin(θ i )) with respect to . It is easily checked that encloses the circle centred at (0, 0) with a radius of ρ 2 , and hence Ē 1,(i2) ,Ē 2,(i2) ,Ē 3,(i2) are guaranteed to be nonnegative. In summary,Ē is specified as · · ·Ē 1,(i2) · · ·Ē 1,(n s 2) 1 3 · · · 1 3Ē 2, (12) · · ·Ē 2,(i2) · · ·Ē 2,(n s 2) (31) This matrix has full rank and the column sum is equal to one, thus confirming that E is DTA-compatible. The following result follows from the above discussion. : l = 1, . . . , n} are linearly independent, locally supported, and form a convex partition of unity on Ω pol . Remark 3.3. As long as is chosen to be a triangle enclosing the first polar ring of control points F i j for a given configuration, we are guaranteed non-negative extraction coefficients. It is only in the interest of standardization that we have chosen to fix as an equilateral triangle with a fixed pattern of vertices.
Given n control points f l ∈ d , d ≥ 3, we can construct a C 1 polar surface f embedded in d , or, equivalently, after pulling back to Ω, The behavior of f at the polar point is going to be fully specified by the first three control points f 1

Double polar point
When dealing with double polar surfaces, the spline construction can be obtained by collapsing a pair of two opposite edges as illustrated in Figure 3 (b). As explained in Section 3.1, the smoothness treatment of the two poles can be done separately for n t ≥ 4. In this case, each pole leads to a local extraction matrix by applying the same procedure as in Section 3.2.1 and the combined global extraction matrix takes the following sparse block-diagonal form: where I is the identity matrix of size n s (n t − 4) × n s (n t − 4) andĒ (i) , i = 1, 2, are matrices of size 3 × 2n s . By choosing the two polar parameterizations F (Ω) andF (Ω) specified in Section 3.1, it is easily verified that one can set whereĒ is the matrix defined in Equation (31) and J k is the exchange matrix of size k × k, i.e., an anti-diagonal matrix of the form The extraction matrix E can then be used to compute the set of spline functions {N l : l = 1, . . . , n} in terms of the tensorproduct functions {B i j : i = 1, . . . , n s ; j = 1, . . . , n t }. Similar to the single-pole result in Corollary 3.4, these spline functions have the following properties.

Refinement of piecewise-rational polar surfaces
A polar spline surface can be refined in a manner similar to the one discussed in Section 2.3. We begin with the following observation. Consider the 3 × 3 matrix for some angles θ ι and θ κ ∈ {θ ι , θ ι + π} selected from the set in Equation (20). This a submatrix ofĒ, defined in Equation (31), consisting of three linearly independent columns. Its inverse is given bȳ Then, we define a sparse matrixD of size 2n s × 3, whose non-zero entriesD i j are identified as follows: for j = 1, 2, 3, From its construction it is clear that the productĒD is equal to the identity matrix. Note that the product (J 3Ē J 2n s )(J 2n sD J 3 ) is also equal to the identity matrix. Hence, keeping in mind the definition of E in Equations (34)-(35), the matrix gives rise to a product ED that is equal to the identity matrix. A similar matrix D can be found (with only two diagonal blocks) for the matrix E defined in Equation (29).
Any refinement matrix of polar splines can then be built using the following procedure. First, we compute the control points of the tensor-product basis functions B i j as in Equation (33). Then, we refine the tensor-product control points using a tensor product of univariate refinement matrices (see Section 2.3). Denote this matrix with S, and denote the polar spline extraction matrices before and after refinement with E andẼ, respectively. Then, control points of the refined polar spline basis functions can be obtained by applying a matrix R to the original set of polar control points, where R is computed by solving the following (overdetermined) linear system with a unique solution, Example 3.8. The box at the top in Figure 4 (a) shows a bidegree (2, 2) unit sphere built by choosing a x = a y = a z = 1, while the box at the bottom shows a bi-degree (2, 2) ellipse with axis lengths (1, 1 2 , 1 3 ) built by choosing a x = 2a y = 3a z = 1. These descriptions use only 8 rational pieces. In each box, the figure at the top shows the exact quadric, while the figure at the bottom shows the deformed quadric obtained by perturbing the control points as per Equation (39). The exact and deformed surfaces are all C 1 smooth at the poles.
3.4.2. C 1 description of degree (2,3) For the second approach, we choose Ω s = [0, 4] and Ω t = [0, 1], and build the univariate rational spline spaces s (periodic) and t on them using the following sets of parameters:  (17) and H t = I 4 , respectively. The full tensor-product extraction matrix H is obtained as in Equation (40). Moreover, the polar extraction operator E is equal to the matrix in Equation (41). The latter matrix maps a total of 16 tensor-product C −1 piecewise-NURBS B j of degree (2, 3) to a total of 6 C 1 polar splines N l . Equivalently, EH maps a total of 48 tensor-product NURBS b j of degree (2, 3) to a total of 6 C 1 polar splines N l . These relations are encapsulated in the following equation, Note that, since we are using a higher-degree representation compared to Section 3.4.1, the control points move farther away from the spline surface, mimicking the behavior of classical NURBS.
Example 3.9. The box at the top in Figure 4 (b) shows a bidegree (2, 3) unit sphere built by choosing a x = a y = a z = 1, while the box at the bottom shows a bi-degree (2, 3) ellipse with axis lengths (1, 1 2 , 1 3 ) built by choosing a x = 2a y = 3a z = 1. These descriptions use 4 rational pieces. In each box, the figure at the top shows the exact quadric, while the figure at the bottom shows the deformed quadric obtained by perturbing the control points as per Equation (39). The exact and deformed surfaces are all C 1 smooth at the poles. (3,3) Finally, for the third approach, we choose Ω s = [0, 2] and Ω t = [0, 1], and build the univariate rational spline spaces s (periodic) and t on them using the following sets of parameters: The corresponding piecewise-NURBS extraction operators are H s = H per defined in Equation (18) and H t = I 4 , respectively. The full tensor-product extraction matrix H is obtained as in Equation (40). Moreover, the polar extraction operator E is equal to the matrix in Equation (41). The latter matrix maps a total of 16 tensor-product C −1 piecewise-NURBS B j of degree (3, 3) to a total of 6 C 1 polar splines N l . Equivalently, EH maps a total of 32 tensor-product NURBS b j of degree (3, 3) to a total of 6 C 1 polar splines N l . These relations are encapsulated in the following equation, Observe again that, since we are using a higher-degree representation compared to Sections 3.4.1 and 3.4.2, the control points move even farther away from the spline surface, mimicking the behavior of classical NURBS. Remark 3.11. The examples presented here have focused on the simplest possible C 1 descriptions of quadrics, namely descriptions that either use lowest-degree splines -bi-degree (2, 2) -or the smallest number of polynomial pieces -two. Unsurprisingly, these simplest descriptions can lead to large control triangles since each control point influences a large portion of the spline surface. Nevertheless, localized control is easily attained upon refinement (see Section 3.3) and, in particular, refinement also leads to much smaller control triangles that offer much finer geometric control. The surface shown in Figure 5 illustrates this point. This surface has been obtained by refining and modifying the control points of the bi-degree (2, 2) sphere from Figure 4

Conclusions
We have presented a general class of C 1 smooth rational splines that allow for the construction and refinement of C 1 smooth curves and (polar) surfaces. They are built by gluing together multiple sets of NURBS basis functions with C 1 smoothness using a DTA-compatible extraction matrix. The main features of the splines we have built are the following: • all standard properties of NURBS, including support for intuitive control-point-based design, • (local) degree elevation and knot insertion based on classical NURBS refinement, • low-degree C 1 descriptions of exact ellipses and ellipsoids, and • compatibility with CAD or CAE software through the explicit representation in terms of NURBS.
In particular, with regard to the last two bullets above, we believe that the explicit, NURBS-compatible C 1 descriptions of ellipses and ellipsoids provided herein will be of use to geometric modellers [11] and computational scientists [24] alike. For instance, the exact C 1 (re)parameterizations at polar points may make the design of algorithms more stable and efficient; it may also avoid the need for special treatment of polar points.