Non-uniform interpolatory subdivision schemes with improved smoothness

Subdivision schemes are used to generate smooth curves or surfaces by iteratively reﬁning an initial control polygon or mesh. We focus on univariate, linear, binary subdivision schemes, where the vertices of the reﬁned polygon are computed as linear combinations of the current neighbouring vertices. In the classical stationary setting, there are just two such subdivision rules, which are used throughout all subdivision steps to construct the new vertices with even and odd indices, respectively. These schemes are well understood and many tools have been developed for deriving their properties, including the smoothness of the limit curves. For non-stationary schemes, the subdivision rules are not ﬁxed and can be different in each subdivision step. Non-uniform schemes are even more general, as they allow the subdivision rules to be different for every new vertex that is generated by the scheme. The properties of non-stationary and non-uniform schemes are usually derived by relating the scheme to a corresponding stationary scheme and then exploiting the fact that the properties of the stationary scheme carry over under certain proximity conditions. In particular, this approach can be used to show that the limit curves of a non-stationary or non-uniform scheme are as smooth as those of a corresponding stationary scheme. In this paper we show that non-uniform subdivision schemes have the potential to generate limit curves that are smoother than those of stationary schemes with the same support size of the subdivision rule. For that, we derive interpolatory 2-point and 4-point schemes that generate C 1 and C 2 limit curves, respectively. These values of smoothness exceed the smoothness of classical interpolating schemes with the same support size by one.


Introduction
Given some initial data f 0 = (f 0 i ) i ∈ at level zero with f 0 i ∈ , a linear, uniform, stationary, binary subdivision scheme S generates refined data f k +1 = (f k +1 i ) i ∈ at level k +1 for any k ∈ 0 according to the two subdivision rules with certain weights a j , b j ∈ for j ∈ .These rules are often referred to as the even and the odd subdivision rule, because they determine the refined data with even and odd indices, respectively.In most applications, only a finite number of these weights are non-zero, and we shall adhere to this assumption.The scheme S is called convergent, if the sequence (F k ) k ∈ 0 of functions F k : → that are piecewise linear over /2 k and interpolate f k i at the parameter values t k i = i /2 k for i ∈ converges uniformly for any initial data f 0 and if the limit function is non-trivial for at least one choice of f 0 [15].The function S ∞ f 0 = lim k →∞ F k , which is necessarily continuous, is then called the limit function of S for f 0 .The scheme S can also be applied to initial data consisting of points rather than scalar values, by simply applying the subdivision rules to each coordinate, thus resulting in limit curves, and the piecewise linear interpolant F 0 is then called the control polygon of the limit curve.
Apart from early work by de Rham [10], the first subdivision schemes were designed to generate uniform B-spline curves of arbitrary degree [23].For these B-spline schemes, the weights a j and b j of the subdivision rules in (1) can be derived from the two-scale relations of the uniform B-spline basis functions, and the properties of the limit functions or curves are evident.Later research studied subdivision independently of B-splines, and simple tools are now available for checking the necessary and sufficient conditions for the convergence of S and the smoothness of the limit functions or curves [5,15], as well as for deriving the degrees of polynomial generation and reproduction [13], which are intimately related to the approximation order of the scheme S .
In order to generate more general classes of limit curves, several modifications of the rules in (1) have been considered.In the non-stationary setting, the weights of the rules are allowed to depend on the subdivision level k , and this flexibility can be used to reproduce circles [21] or exponential polynomials [17].Circles can as well be generated by non-linear subdivision schemes [28], where the subdivision rules may depend on the data f k i and hence vary not only for each subdivision level k , but also for each index i .Besides generating circles and Gaussian functions [31], non-linear schemes are capable of eliminating shape artefacts [12,24] and creating smooth curves on manifolds and curved surfaces [25].A third modification are non-uniform subdivision schemes [8,35], which can be used, for example, to generate non-uniform B-splines [4,19,30].For these schemes, the parameter values t k i corresponding to the data f k i are usually no longer assumed to be uniformly spaced, and consequently the subdivision rules depend on k and i , but they remain linear.
Several tools have been developed for proving the convergence and smoothness of a non-linear, nonstationary, or non-uniform scheme S , including divided differences [8], asymptotic equivalence [14,18], asymptotic similarity [7], and proximity conditions [32,34].In a nutshell, the common aim of these techniques is to establish a relation to a corresponding linear, uniform, stationary scheme S and to show that convergence and smoothness properties of S are inherited by S .For example, non-uniform corner cutting gives limit curves that are as smooth as Chaikin's quadratic B-spline scheme [6,27], namely C 1 , under certain constraints of the cut proportions [20], and several variants of the classical interpolatory 4-point scheme [11,16] are known to also give C 1 limit functions or curves [2,3,8,24].
In this paper, we show that it is possible to go beyond these results and to increase the smoothness of a linear scheme, if we allow the weights a j and b j of the subdivision rules in (1) to depend on k and i .That is, we consider non-uniform schemes, but in contrast to most other non-uniform schemes, we stick to the uniform, dyadic parameter values t k i = i /2 k .In particular, we focus on schemes that are interpolatory, that is, with weights a k i ,0 = 1 and a k i , j = 0 for j = 0, independently of k and i , as well as periodic, in the sense that the weights b k i , j repeat at a frequency of 2 k in i , so that the same refinement strategy is applied over each of the initial intervals Overall, this means that we consider subdivision schemes that generate the refined data ) i ∈ at level k + 1 for any k ∈ 0 according to the even, interpolatory subdivision rule and the 2 k odd subdivision rules where l = i mod 2 k , with certain finite sets J k i ⊂ and weights b k i , j ∈ for i ∈ {0, 1, . . ., 2 k − 1}, j ∈ J k i , and k ∈ 0 .We further call the subdivision scheme an n -point scheme, if the cardinality of all sets J k i is at most n .Within this setting, we first derive two novel interpolatory 2-point schemes with C 1 limit functions, as opposed to the stationary, interpolatory 2-point scheme that simply reproduces the piecewise linear interpolant F 0 of the initial data and hence is only C 0 (Section 2).The main idea behind our construction is inspired by the B-spline schemes, which were designed to converge to uniform B-spline curves, in that we also start with a certain function that we would like the non-uniform subdivision scheme to generate in the limit and then derive the appropriate weights b k i , j that achieve this.We then use this approach to derive two interpolatory 4-point schemes with C 2 limit functions (Section 3), and finally discuss how these schemes can be used for generating interpolating curves on surfaces (Section 4).

Non-uniform C 1 interpolatory 2-point schemes
Let us start by considering interpolatory 2-point schemes, where J k i = {−1, 0} for any k ∈ 0 and i ∈ {0, 1, . . ., 2 k − 1} in (3), so that each new data value f k +1 2i +1 is some linear combination of the neighbouring data values f k i and f k i +1 .It follows by induction over k ∈ 0 that f k i for any i ∈ {2 k m, . . ., 2 k (m + 1)} is some linear combination of f 0 m and f 0 m+1 , which implies that the whole limit function over the interval [m, m +1] depends only on these two values and is independent of the limit function over the neighbouring intervals [m − 1, m ] and [m + 1, m + 2].Therefore, the only chance to generate a limit function that is C 1 at m ∈ is to force the limit function to match a prescribed derivative at m over both adjacent intervals.The simplest unbiased choice is to set the derivative of the limit function to zero at all m ∈ .Hence, we would like to construct a subdivision scheme that generates, over each interval [m, m + 1], a limit function that interpolates the initial data and has vanishing first derivative at the endpoints of this interval.
Without loss of generality, let us focus on the unit interval [0, 1].Our goal now is to determine the weights of the odd subdivision rules in (3), such that the subdivision scheme generates a certain C 1 function p as the limit function over [0, 1].
Theorem 1.Let S be the non-uniform, periodic, interpolatory 2-point scheme with the odd subdivision rules in (3) and weights Then the scheme S generates data is convergent, and the limit functions f = S ∞ f 0 are C 1 functions with f (m) = f 0 m and f (m) = 0 for m ∈ .
Proof.We first prove (7) by induction over k and note that the base case (k = 0) follows from (4).For the inductive step, we distinguish the even and the odd case.On the one hand, by (2), the induction hypothesis, and the fact that On the other hand, by (3) with the weights in (5), the induction hypothesis, and (6), we have It then follows from (7) that the functions F k are piecewise linear interpolants of the function p at the dyadic points t k i and as such converge uniformly to p over [0, 1] as k → ∞ [9].By the periodicity of the scheme, the same applies to all intervals [m , m + 1] for m ∈ , thus proving the rest of the statement.It remains to understand under which conditions the weights in (5) can always be chosen such that they satisfy condition (6).To this end, let us express p in terms of the initial data f 0 0 and f 0 1 as Comparing the coefficients of f 0 0 and f 0 1 on both sides of ( 6) then yields the linear system Lemma 2. If r 0 and r 1 form a partition of unity and r 1 is strictly increasing, then the linear system in (8) has a unique solution for any k ∈ 0 and i ∈ {0, 1, . . ., 2 k − 1}, which is given by These weights α k i and β k i are both positive and sum to one.
Proof.Using the partition of unity property, we find that det( . This implies the invertibility of A k i for any k ∈ 0 and i ∈ {0, 1, . . ., 2 k − 1}.The formulas in ( 9) can be derived from (8) after inverting A k i explicitly, and α k i and β k i are positive, because r 1 is strictly increasing and They obviously sum to one.
Under the assumptions of Lemma 2, it follows that the odd subdivision rules (3) of the corresponding 2point scheme generate the new data f k +1 2i +1 as a convex combination of the previous data and any m ∈ .We shall now discuss two particular choices of basis functions r 0 and r 1 that lead to rather simple expressions of the weights α k i and β k i in (9).
Example 1.One possibility to satisfy the conditions in ( 4) is to choose as p the cubic Hermite interpolant, which is given by the basis functions These functions form a partition of unity and r 1 is strictly increasing, because r 1 (x ) = 6x (1 − x ) > 0 for x ∈ (0, 1).A straightforward calculation reveals that the formulas in ( 9) can be simplified to and they converge to lim next to the centre of the interval.
Because of this limit behaviour, the known tools for analysing the smoothness of non-uniform schemes are not applicable in this case, but it is evident from the construction that the limit functions of this nonuniform, periodic, interpolating 2-point scheme are the cubic C 1 B-splines with duplicated control points . .) and uniform double knots (. . ., −1, −1, 0, 0, 1, 1, . ..).

Example 2.
Even simpler weights can be found if we take as p the piecewise quadratic C 1 function given by the basis functions for which the weights in (9) turn out to be (α 0 0 , β 0 0 ) = 1 2 , 1 2 and 1  4 not only in the limit, but for any k > 0, and likewise (α 3  4 .As in Example 1, the known tools are not able to prove the smoothness of the limit functions generated by this non-uniform, periodic, interpolating 2-point scheme, but it follows from the construction that they are the quadratic C 1 B-splines with duplicated control points (. . ., In the geometric setting, the fact that each new point f k +1 2i +1 is generated as a convex combination of the old points f k i and f k i +1 implies that the limit curve has the same shape as the initial control polygon, with apparent C 1 discontinuities at the initial control points.The explanation for this behaviour is that the new points generated by the non-uniform scheme are non-uniformly distributed along the edges of the initial control polygon and clustered towards the endpoints, which differs from the uniform distribution of points given by the classical interpolatory 2-point scheme with the stationary weights 1  2 , 1 2 (see Figure 1).Hence,  the limit curve of the non-uniform scheme has a continuous rather than a piecewise constant derivative with respect to the dyadic parameterization, with vanishing first derivative at m ∈ .This fact is more evident in the functional setting, where it leads to flat spots at the integers (see Figure 2).We should stress that we do not believe these non-uniform 2-point schemes to be of any practical use, but rather consider them as "academic examples", which showcase that non-uniform schemes have the potential to generate limit functions and curves with higher smoothness than those generated by the corresponding stationary scheme.And it turns out that this potential is not restricted to 2-point schemes.

Non-uniform C 2 interpolatory 4-point schemes
For the non-uniform interpolatory 2-point schemes above, we can guarantee the limit functions to be C 1 only by forcing the first derivative to vanish at the integers.The reason for this arguably arbitrary choice is the locality of the subdivision rule, but we can overcome this problem if we enlarge the support of the rules.
Hence, let us now consider the interpolatory 4-point schemes with J k i = {−2, −1, 0, 1} and weights for the odd subdivision rules in (3), so that the new data value , and f k i +2 , just as for the classical interpolatory 4-point scheme [11], which uses the stationary weights − 1 16 , 9  16 , 9  16 , − 1 16 .As before, we first focus on the unit interval [0, 1], and this time we aim at reproducing in the limit a certain C 2 function p : [0, 1] → that satisfies Again, the key trick is to choose the weights in (10) such that Expressing p in terms of the initial data f 0 −1 , f 0 0 , f 0 1 , and f 0 2 as and comparing the coefficients of f 0 −1 , . . ., f 0 2 in (12) leads to the linear system Lemma 3. If r −1 , r 0 , r 1 , r 2 reproduce quadratic polynomials in the sense that and the third-order divided difference of r 2 satisfies for some k ≥ 2 and i ∈ {1, 2, . . ., 2 k − 2}, then the corresponding linear system in (13) has a unique solution.
Proof.Using ( 14), we first transform We then exploit the fact that the parameter values t k i −1 , . . ., t k i +2 are uniformly spaced with distance 2 −k , so that replacing the second, third, and fourth column of this matrix with the first, second, and third order differences of its columns gives the lower triangular matrix hence B k i is non-singular, because of (15).
As in the proof of Theorem 1, we would now like to conclude by induction that the scheme generates data which in turn implies that p is the limit curve over [0, 1], and for the even case the arguments are the same.Alas, the inductive step for showing that ) relies on the identity f k −1 = p (t k −1 ), which is not supported by the induction hypothesis for k > 0, and likewise for f k +1 2 k +1 −1 .To fix this problem, let us take a look at the first subdivision step.It is not hard to verify that for k = 0 and i = 0, the correct weights, which guarantee generating the data Therefore, computing f 2 1 in the standard fashion as results in a value that depends on the initial data f 0 −2 , . . ., f 0 2 , but in order to have f 2 1 = p (t 2 1 ), it must not depend on f 0 −2 .The trick now is to replace the index set J 1 0 = {−2, −1, 0, 1} by the set J 1 0 = {−2, −1, 0, 2}, so that the subdivision rule (3) for f 2  1 becomes Comparing the coefficients of f 0 −1 , . . ., f 0 2 with those on the right hand side of then reveals that for k = 1 and i = 0 the correct choice of weights To compute f 2 3 , we proceed analogously, using the index set J 1 1 = {−3, −1, 0, 1}, and it turns out that the correct weights b are So far, this ensures the interpolation property in ( 16) to hold for k = 0, 1, 2, but for k > 2, we need yet another approach.If we replace again the index set J k 0 for k ≥ 2 with {−2, −1, 0, 2}, then we run into the same problem as before, because To resolve this, we replace J k 0 with {−3, −2, −1, 0} instead, so that the subdivision rule for f k +1 1 becomes Note that this choice would not work for k = 1, because f 1 3 depends on f 0 3 .Using the induction hypothesis that f k j = p (t k j ) for j = 0, 1, 2, 3, the weights ) are then given by the linear system in ( 13) for i = 0, with B k 0 replaced by B k 1 , which still leads to a unique solution under the assumptions of Lemma 3. Similarly, we replace the index set J k 2 k −1 with {−1, 0, 1, 2} for k ≥ 2, and the weights , are given by (13 Based on the discussion above, we are now ready to summarize our construction.Theorem 4. Let r −1 , r 0 , r 1 , r 2 be four basis functions that satisfy as well as the assumptions of Lemma 3 for any k ≥ 2 and i ∈ {1, 2, . . ., 2 k − 2}.Further let S be the non-uniform, periodic, interpolatory 4-point scheme with the odd subdivision rules in (3) and the weights in (17), (18), and (19) for k = 0 and k = 1, as well as Then the scheme S generates the data in (16), is convergent, and the limit functions are C 2 functions with f As in the previous section, we shall now discuss two particular choices of basis functions r −1 , r 0 , r 1 , r 2 , for which rather simple expressions of all weights α k i , β k i , γ k i , δ k i can be derived.
Example 3. One possibility to satisfy the conditions in (11) is to choose as p the quintic Hermite interpolant, which is given by the basis functions It is not hard to verify that these functions satisfy the conditions in ( 14) and ( 20) and that For k = 0, these values are imaginary, and if k > 0, then 2 k − 1 is odd and the root would have to be odd, too, in order for i to be an integer, that is, for some m ∈ , but this is impossible, because 4 k does not divide by 5.This asserts that these basis functions satisfy the assumptions of Theorem 4, and the weights of the resulting subdivision scheme turn out to be for k = 0 and k = 1, as well as 15  16 , − 5 16 , 1   16   and lim at the endpoints and to the weights of the classical interpolatory 4-point scheme, 9  16 , 9  16 , − 1  16 , at all interior points with i ∈ {1, 2, . . ., 2 k − 2}.
As the limit weights in (21) are the non-symmetric variants of the stationary 4-point weights, in the sense that they can also be derived from locally fitting and evaluating a cubic polynomial, it might be possible to extend one of the known tools to show that the limit functions of this non-uniform, periodic, interpolating 4-point scheme are as smooth as those of the stationary 4-point scheme, namely C 1 , but not that they are C 2 .In fact, it is not hard to see that the limit functions are the quintic C 2 B-splines with control points . . ., A simple calculation reveals that these functions satisfy the conditions in ( 14) and (20) and that for the new points at m for the new points at m + 1 4 − 2 −k −1 and m + 3 4 + 2 −k −1 .All other b k i are just the classical stationary weights − 1 16 , 9  16 , 9  16 , − 1  16 , because they are used to compute a new point from four old points that lie on a common cubic piece of the limit function.
Because of the special weights in (22), the usual tools for analysing the smoothness of non-uniform schemes cannot be applied to this scheme, but it is clear by construction that the limit functions of this non-uniform, periodic, interpolating 4-point scheme are the cubic C 2 B-splines with control points . . ., , . . .and uniform knots . . ., −1, − 3 4 , − 1 2 , − 1 4 , 0, 1 4 , 1 2 , 3 4 , 1, . . .= /4.Figures 3 and 4 show a comparison between the classical interpolatory 4-point scheme and the two nonuniform variants from Examples 3 and 4 in the geometric and the functional setting, respectively.While the shapes of the limit functions for the non-uniform schemes are slightly negatively affected by the prescribed first derivative at the integers, the limit curves look rather pleasing and smoother than the limit curve of the classical scheme, as expected due to the improved smoothness.Note that the limit functions and curves generated by the two non-uniform schemes are visually almost indistinguishable, but the differences are evident for the second derivatives of the limit functions, which are piecewise cubic and piecewise linear, respectively (see Figure 4).

Curves on surfaces
One may argue that the schemes derived above are questionable, because the exact limit curve is known a priori and can thus be evaluated directly.However, the same argument would apply to the uniform B-spline schemes [23], and the situation changes once we leave the non-Euclidean setting, for example, if we consider subdivision on a manifold surface M. In this case, the usual definition of the affine average of n points p 1 , p 2 , . . ., p n for certain weights w i that sum to one,  is not well-defined anymore, because the points on M do not form a vector space.Instead, one can define the affine average using the weighted Fréchet mean, where d : M×M → is the distance function, which measures the length of the geodesic, that is, the shortest path between two surface points.If d is the Euclidean distance, then both definitions are equivalent, but in general, the minimizer p in (23) may not be unique and may not depend continuously on the points p i and the weights w i [22,1].Therefore, even though the limit curves of the non-uniform schemes in Examples 1-4 are B-spline curves and every point on these curves is an affine combination of the initial points f 0 i , evaluating them directly via the Fréchet mean can result in curves that are not even C 0 .If we generate the limit curve instead using a subdivision process, then we are guaranteed to have a sequence of continuous, piecewise geodesic curves, which usually converge to a limit curve that is as smooth as the one in the Euclidean setting [34,33].Moreover, in order to properly define the subdivision rules for points on surfaces, it suffices to have a generalization of the affine average of two points.In the Euclidean setting, the affine average A(p , q , t ) = (1 − t )p + t q (24) of p and q with weights (1 − t ) and t is simply the unique point r = A(p , q , t ) on the straight line through p Since the equivalent of a straight line on a surface M is a geodesic, it is then clear how to generalize the average operator in (24) to two points p , q ∈ M. Note that if t ∈ [0, 1], then r = A(p , q , t ) is a point on the geodesic between p and q , while this geodesic needs to be extrapolated beyond p for t < 0 and beyond q for t > 1.A formal description of this procedure is given by Wallner [33] for smooth surfaces and by Polthier and Schmies [26] for surface meshes.The ambiguity of this definition, which may arise whenever the geodesic between p and q is not unique, can be resolved by making an arbitrary but deterministic choice.With the average operator A in (24) defined on M, we now follow Schaefer and Goldman [29] and express the odd subdivision rules (3) of our non-uniform 4-point schemes with J k l = { j 2 , j 1 , j 0 , j −1 } in terms of three 2-point-averages as For this order of averages, the first two averages are usually geodesic interpolations with µ 1 , µ 2 ∈ (0, 1) and only the last average requires geodesic extrapolation with a rather small extrapolation factor −µ 3 > 0. For example, this is the case (with −µ 3 = 1 8 ) for almost all subdivision rules of the non-uniform schemes in Examples 3 and 4.Only at the endpoints of each interval, for k > 1 and i ∈ {0, 2 k − 1}, this order applies extrapolation for the first or the second (or both) averages and interpolation for the third average.
When implementing the scheme from Example 4, we observed that the rather large extrapolations that may occur during the first three subdivision steps (e.g., with factor −µ 2 = 41 10 for k = 1 and i = 1 and −µ 1 = 163 230 for k = 3 and i = 0) can lead to some unexpected artefacts of the limit curve.We therefore propose to circumvent these subdivision steps and to compute instead the data f 4 = (f 4 i ) i ∈ directly as affine combinations of the initial data, using the fact that for any i ∈ .This avoids extrapolation with factors greater than 1  8 and generally results in artefact-free and presumably smooth curves (see Figure 5).Proving that these limit curves are C 2 is beyond the scope of this paper and left to future work.

Figure 1 :
Figure 1: Initial points, points generated by the first four subdivision steps, and limit curve (from left to right) for the stationary (top) and the non-uniform interpolatory 2-point schemes from Example 1 (middle) and Example 2 (bottom).

Figure 2 :
Figure 2: Initial data, data generated by the first three subdivision steps, limit function, and first derivative of the limit function (from left to right) for the stationary (top) and the non-uniform interpolatory 2-point schemes from Example 1 (middle) and Example 2 (bottom).

Figure 3 :
Figure 3: Initial points, points generated by the first four subdivision steps, and limit curve (from left to right) for the stationary (top) and the non-uniform interpolatory 4-point schemes from Example 3 (middle) and Example 4 (bottom).

Figure 4 :
Figure 4: Initial data, data generated by the first two subdivision steps, limit function, first and second derivative of the limit function (from left to right) for the stationary (top) and the two non-uniform interpolatory 4-point schemes from Example 3 (middle) and Example 4 (bottom).

Figure 5 :
Figure 5: Initial data (black dots) and interpolating curves on a manifold surface after 8 subdivision steps with the stationary (blue) and the non-uniform (red) 4-point scheme from Example 4.