A shape-preserving variant of Lane-Riesenfeld algorithm

1.   Introduction

Subdivision schemes (SSs) are systematic algorithms for producing smooth and pleasant curves or surfaces from the set of initial control points or control array. They are easy to use, simple to investigate, and highly flexible. With these and other attractive mathematical properties of SSs, their popularity is increasing in applications, such as in computer graphics, signal and image processing, computer animation and commercial industry. Generally, SSs can be categorized as approximating and interpolating schemes. Approximating schemes generate more smooth curves and surfaces as compared to interpolating ones.

Lane and Riesenfeld ^[16] introduced an algorithm for creating uniform B-spline curve of order $k$ $(k\in\mathbb{N})$. The algorithm executes on the set of initial control points in two steps. In the first step, initial data is refined by simply doubling the initial control points. This step is followed by successive smoothing stages by taking their mid-point averaging. Due to its simple and efficient execution, the Lane-Riesenfeld algorithm (LRA) became a celebrated algorithm to generate smooth curves and surfaces. Many researchers proposed diverse variants of LRA to satisfy different requirements. For instance, Schaefer et al. ^[18] applied binary non-linear averaging operator instead of binary linear smoothing operator to get SSs for exponential functions. Hormann and Sabin ^[9] constructed a family of SSs with cubic precision having a construction similar to that of uniform B-spline schemes.

Cashman et al. ^[5] provided another generalization of LRA in which the same operator was used for both the refinement and the smoothing stage. Ashraf et al. ^[1] offered a local quintic interpolant variant of LRA to develop a family of SSs which reproduce quintic polynomials. Mustafa et al. ^[12] presented their generalization of LRA by using an interpolating SS ^[6] as a refine operator and an approximating SS ^[9] as smoothing operator. Romani ^[15] presented a variant of LRA in which a perturbed version of Chaikin's scheme was used in the refining stage and smoothing stages remained unchanged. Mustafa and Rabia ^[10] used quartic B-spline based SS as a refining operator and mid-point averaging as a smoothing operator. Shi et al. ^[17] presented a parameter-dependent variant of LRA.

Generally speaking, $C^2$ continuity is enough in industrial design, but the design of precision instruments in manufacturing industry requires higher continuous curves and surfaces, which motivated us to present a shape-preserving variant of LRA. The new variant is used in the refining stage and as well as in the smoothing stage and generates a family of shape-preserving SSs with high continuity. Shape-preserving SSs have remarkable significance in geometric shape design. They are widely used in the design of curves to predict and obtain their shape as per the shape of initial control points. Although there are many shape-preserving properties, monotonicity and convexity preservation are two foremost shape-preserving properties. For further insight into shape-preserving SSs, one may refer to ^[2,3,4,8].

The proposed family of SSs preserves both monotonicity and convexity, where different family members are obtained by setting one parameter in the expression of the refinement rule. Every family member has different properties compared to the others and so, by choosing the parameter, one can pick the best-suited scheme according to the context of the application. This paper is organized as follows: In Section 2, a new variant of LRA with the refine stage and the smoothing stage composed of a shape-preserving scheme is proposed. In Section 3, the analysis of some basic properties of the proposed schemes is presented. The shape-preserving properties of the variant are given in Section 4. In Section 5, the limit stencil and the artifact analysis of the proposed schemes are presented. In Section 6, we discuss some numerical examples. Conclusions are drawn in Section 7.

2.   Construction of the new variant

LRA executes by first applying a refine stage $R_L$ on the set of initial control points $Y$ = $\{y_{i}\}_{i\in \mathbb{Z}}$ and then by applying $p$ smoothing stages $Q_L$. Thus a refinement stage $T$ can be considered as $T = Q^{p}R$, where in refining stage $R_L$ every control point is doubled such that

and smoothing stage $Q_L$ is based on mid-point averaging rule, i.e.,

Now we propose a new variant of LRA which is based on a shape-preserving SS. ^[14] suggested a relaxed four-point SS preserving both monotonicity and convexity of initial data. This scheme offers a refine operator $R_s$ as

and smoothing operator $Q_s$ as

By applying a single time refine stage followed by $p$-time smoothing stages, we have a family of shape-preserving SSs, named as $S_p$-scheme. The symbol (Laurent polynomial) of the $S_p$-scheme is

where $Q(z)$ is the symbol of smoothing stage $Q_{s}$

and $R(z)$ is the symbol of refine stage $R_{s}$

The parameter $p$ enumerates the number of applied smoothing stages. It also classifies family members. The family member with smoothing stage $p = 0$ is the primary four-point scheme ^[14]. Table 1 presents the mask (set of non-zero coefficients of the scheme) of first four family members by substituting $p = 0, 1, 2$ and $3$ in (2.5)–(2.7).

3.   Basic properties of the $S_p$-scheme

In this section, we present smoothness analysis, Hölder continuity analysis, and the support of the basic limit function of the $S_p$-scheme.

3.1.   Smoothness analysis

Smoothness analysis is a fundamental question for subdivision schemes. Dyn and Levin ^[7] presented an effective method to verify convergence by using the symbol $S_p (z)$. Later this method became known as Laurent polynomial method. We adopt the main idea of this classical approach and further analyze the smoothness of the $S_p$-scheme.

Theorem 3.1. The $S_p$-scheme has $C^{p+3}$ continuity for $p = 0, 1, 2, \ldots$.

Proof. The symbol of the $S_{p}$-scheme can be written as

with

Let $S_{\xi_{p}}$ be the subdivision scheme corresponding to the symbol $\xi_{p}(z)$, defining

where $\xi_{i}^{[l]}$ denote the coefficients of the symbol $\xi_{p}^{[l]}(z)$, and $\xi_{p}^{[l]}(z) = \xi_{p}(z)\xi_{p}(z^2)\ldots\xi_{p}(z^{2^{l-1}})$.

For $p = 0$, we have

where

Let $S_{\xi_{0}}$ be the scheme corresponding to the symbol $\xi_{0}(z)$ and $\left\|S_{\xi_{0}}^2\right\|_{\infty} = \frac{3}{4} < 1$. Since $S_{\xi_{0}}$ is contractive, so by Dyn and Levin (^[7], Corollary 4.14), $S_0$-scheme has $C^3$ continuity.

Similarly for other values of $p$, we can easily have if the subdivision scheme $S_{\xi_{p}}$ is contractive, i.e. $\left\|S_{\xi_{p}}^{l}\right\|_{\infty} = \max\left\{\sum\limits_{j\in \mathbb{Z}}|\xi_{i-2^lj}^{[l]}|: 0\leq i < 2^l \right\} < 1,$ then $S_p$-scheme has $C^{p+3}$ continuity.

3.2.   Hölder continuity analysis

Hölder continuity can be considered as extended continuity. There are upper and lower bounds on it. First, we discuss the upper bound on the Hölder continuity of the $S_p$-scheme which is based on the idea presented in Rioul ^[13] and Dyn and Levin ^[7].

Theorem 3.2. The upper bound on Hölder continuity of the $S_p$-scheme is

where $p = 0, 1, 2, \ldots,$ and $\zeta_{p}$ is the joint spectral radius of the matrices $Q_0$ and $Q_1$ where these matrices can be obtained by using symbol of the $S_p$-scheme.

Proof. By (2.4), symbol of the $S_p$-scheme is given by

where $Q_{p}(z) = \left(\frac{-z^{-1}+18-z }{16}\right)^{p}(-z^{-1}+4-z)$ and $p = 0, 1, 2, \ldots$. Let $q_0, q_1, \ldots, q_d$ be the non-zero coefficients of $Q_{p}(z)$. Also let $Q_0$ and $Q_1$ be the matrices of order $d\times d$. The entries of these matrices are given by

Let $\zeta_{p}$ be the joint spectral radius of the matrices $Q_0$ and $Q_1$, then by Rioul ^[13] and Dyn and Levin ^[7] upper bound on Hölder continuity of the $S_{p}$-scheme is $p+6-\log_2(\zeta_{p})$.

The following theorem estimates the lower bound on the Hölder continuity of the $S_p$-scheme.

Theorem 3.3. The lower bound on the Hölder continuity of the $S_p$-scheme is

where $p = 0, 1, 2, \ldots.$

Proof. The symbol of the $S_p$-scheme is given by

where $U(z) = (g(z))^{p}h(z)$, $g(z) = \frac{-z^{-1}+18-z }{16}$ and $h(z) = -z^{-1}+4-z.$ So Hölder continuity of $S_p$-scheme is bounded from below by

Since $g(z)$ and $h(z)$ both are alternating symbols, so their product $U(z)$ is also alternating and

where $U^{\diamond}$ and $U_{\diamond}$ are sum of even and odd coefficients of $U(z)$ respectively. In matrix-vector notation they can be written as

thus we have

By eigenvalue decomposition, we have

which implies that

thus we have

Therefore, lower bound on the Hölder continuity of the $S_p$-scheme is

where $p = 0, 1, 2, \ldots.$

Upper and lower bounds on the Hölder continuity of the $S_{p}$-scheme for different values of $p$ can be easily calculated by using Theorem 3.2 and Theorem 3.3 respectively. Table 2 summarizes the continuity analysis of the $S_{p}$-scheme. It is clear from Table 2 that as we increase value of $p$, levels of smoothness and Hölder continuity of the $S_{p}$-scheme go up steadily with $p$.

3.3.   Support of basic limit function

Proposition 3.4. Support width of basic limit function of the $S_p$-scheme is $3p+8.$

Proof. The basic limit function of a convergent subdivision scheme $S_p$ is defined as the limit function $S^{\infty}_p\delta$, where $\delta$ be the Kronecker delta sequence. Since the number of non-zero entries in the masks of smoothing stage $Q_{s}$ and refine stage $R_{s}$ are four and nine respectively, therefore by ^[11] support widths of basic limit functions for $Q_{s}$ and $R_{s}$ are three and eight respectively. Since the mask of $S_p$-scheme is obtained by applying refine stage $R_{s}$ on the initial data followed by $p$-smoothing stages $Q_{s}$. Hence the support width of the basic limit function for the $S_p$-scheme is $3p+8$.

4.   Shape-preserving properties of the $S_p$-scheme

Schemes holding shape-preserving properties are very useful in curve designing. Monotonicity and convexity preserving properties are two of the most important shape-preserving properties. The shape-preserving variant of LRA generates a family of shape-preserving SSs. For this, we show that $S_1$-scheme produces monotonic and convex curves under certain conditions imposed on the initial data. Considering $p = 1$ in (2.5), the $S_1$-scheme is as follows:

4.1.   Monotonicity

Monotonicity preserving property of a SS means that if we choose initial data monotone then the limit curve generated by using this data should also be monotone. To show that the $S_1$-scheme holds monotonicity preserving property we consider its first-order divided difference (DD), i. e., $n_{i}^k = y_{i+1}^k-y_{i}^{k}$. So the first order DD of the $S_1$-scheme is given by:

and

In the following theorem, we discuss the monotonicity preserving property of the $S_1$-scheme.

Theorem 4.1. Let the initial control points $\{y_i^0\}_{i\in \mathbb{Z}}$ be strictly monotone increasing such that

and $r_i^k = \frac{n_{i+1}^k}{n_i^k}, \ R^k = \max\{r_i^k, \frac{1}{r_i^k}\}, \ \forall \ k\in\mathbb{N}_0, i\in \mathbb{Z}$. Also, let $\lambda$ be such that $\lambda\in(-\infty, \frac{233-\sqrt{127417}}{66})\cup(1, \frac{357+\sqrt{146713}}{86})$. If $\frac{1}{\lambda}\leq R^0\leq \lambda,$ and $\{y_i^k\}_{i\in \mathbb{Z}}$ is defined by the scheme (4.1) then,

i. e., the scheme (4.1) generates strictly monotone increasing limit curves.

Proof. We establish the result by using induction method. As we know that the initial control points $\{y_i^0\}_{i\in \mathbb{Z}}$ are strictly monotone increasing so it is obvious that (4.2) is true for $k = 0$.

Now let us suppose that (4.2) is true for $k\geq 1$ and show that it is also true for $k+1$.

Consider

and

So, we have $n_{i}^{k+1} > 0, \ \forall \ i \in \mathbb{Z},$ which implies $n_{i}^{k} > 0, \ \forall \ k \in \mathbb{N}_{0}, i \in \mathbb{Z}$. Now, we show that $\frac{1}{\lambda}\leq R^{k+1}\leq\lambda \ \forall \ k \in \mathbb{N}_{0}$.

Let

we have

where

and

It follows from (4.3), the denominator $D_1$ is $\frac{2048n_{2i}^{k+1}}{n_{i-1}^k} > 0$, and the numerator $C_1$ fulfils

which implies $r_{2i}^{k+1}\leq \lambda, \ \forall \ i \in \mathbb{Z}$.

Since

we have

where

and

It follows from (4.4), the denominator $D_2$ is $\frac{2048n_{2i+1}^{k+1}}{n_{i-1}^k} > 0$, and the numerator $C_2$ fulfils

Therefore $r_{2i+1}^{k+1}\leq \lambda, \ \forall \ i \in \mathbb{Z}$.

In the same way, we can have $\frac{1}{r_{2i}^{k+1}}\leq \lambda$ and $\frac{1}{r_{2i+1}^{k+1}}\leq \lambda$. By induction, we have $R^{k}\leq \lambda, \ \forall \ k \in \mathbb{N}_{0}.$ Considering $R^{k} = \max\{r_{i}^{k}, \frac{1}{r_{i}^{k}}\},$ which leads to $R^{k}\geq\frac{1} {\lambda}, \ \forall \ k \in \mathbb{N}_{0}.$ This completes the proof.

4.2.   Convexity

Convexity preserving property of a SS means that if we choose initial data convex then the limit curve generated by using this data should also be convex. To show that the $S_1$-scheme holds convexity preserving property we consider the second-order DD, i. e., $h_i^{k+1} = 2^{2k+1}(y_{i-1}^{k+1}-2y_{i}^{k+1}+y_{i+1}^{k+1})$. So the second-order DD of the $S_1$-scheme is given by:

and

In the following theorem, we discuss convexity preserving property of the $S_1$-scheme.

Theorem 4.2. Let the initial control points $\{y_i^0\}_{i\in \mathbb{Z}}$ be convex, i.e., $h_i^0 > 0,$ and $m_i^k = \frac{h_{i+1}^k}{h_i^k}, \ M^k = \max\{m_i^k, \frac{1}{m_i^k}\}, \ \forall \ k\in\mathbb{N}_0, i\in \mathbb{Z}$. Also, let $\mu$ be such that $\mu\in(1, \frac{23}{3})$. If $\frac{1}{\mu}\leq M^0\leq \mu,$ and $\{y_i^k\}_{i\in \mathbb{Z}}$ is given by the scheme (4.1) then,

i. e., the scheme (4.1) generates convex limit curves.

Proof. We establish the result by using induction method. As we know that the initial control points $\{y_i^0\}_{i\in \mathbb{Z}}$ are convex so it is obvious that (4.5) is true for $k = 0$.

Now let us suppose that (4.5) is true for $k\geq 1$ and show that it is also true for $k+1$.

Consider

and

So, we have $h_{i}^{k+1} > 0, \ \forall \ i \in \mathbb{Z},$ which implies $h_{i}^{k} > 0, \ \forall \ k \in \mathbb{N}_{0}, i \in \mathbb{Z}$. Now, we show that $\frac{1}{\mu}\leq M^{k+1}\leq\mu \ \forall \ k \in \mathbb{N}_{0}$.

Let

we have

where

and

It follows from (4.6), the denominator $D_3$ is $\frac{512h_{2i}^{k+1}}{h_{i-1}^k} > 0$, and the numerator $C_3$ satisfies,

Thus $m_{2i}^{k+1}\leq \mu, \ \forall \ i \in \mathbb{Z}$.

Since

we have

where

and

It follows from (4.7), the denominator $D_4$ is $\frac{512h_{2i+1}^{k+1}}{h_{i-1}^k} > 0$, and the numerator $C_4$ satisfies,

Therefore $m_{2i+1}^{k+1}\leq \mu, \ \forall \ i \in \mathbb{Z}$.

In the same way, we can have $\frac{1}{m_{2i}^{k+1}}\leq \mu$ and $\frac{1}{m_{2i+1}^{k+1}}\leq \mu$. By induction, we have $M^{k}\leq \mu, \ \forall \ k \in \mathbb{N}_{0}.$ Considering $M^{k} = \max\{m_{i}^{k}, \frac{1}{m_{i}^{k}}\},$ which leads to $M^{k}\geq\frac{1} {\mu}, \ \forall \ k \in \mathbb{N}_{0}.$ This completes the proof.

In the same way, we can show that rest of the members of the family also hold monotonicity preserving and convexity preserving properties.

5.   Limit stencil and artifact analysis of the $S_p$-scheme

In this section limit stencil and artifact analysis of the $S_p$-scheme are presented.

5.1.   Limit stencil analysis

Limit stencil is used to achieve a control point on the limit curve in the form of initial control points. By using limit stencil, it is very convenient to obtain point on limit curve with comparatively less number of computation. We obtain limit stencil by using

and

where $K$ is diagonal matrix of eigenvalues of subdivision matrix of the $S_p$-scheme and eigenvectors corresponding to the eigenvalues are given in the form of matrix $L$. By using (5.1), we calculate limit stencil of the $S_p$-scheme for different values of $p$, which are given in Table 3.

5.2.   Artifact analysis

We consider different data patterns to evaluate the diverse behavior of the subdivision schemes by analyzing how do they respond to these data patterns. For example, we use monotone and convex data for the analysis of shape-preserving properties. Polynomial data led us to determine the degree of polynomial generation and reproduction. Cardinal data, where all the control points have zero value except one control point which has unit value, is used to determine support and shape of basis function. Similarly, when we extract data from a sinusoidal function, we observe how the curvature of the limit curve changes along its arc length. The ripples in the curvature plot are called artifacts of the scheme. Magnitude of artifact $G$ can be calculated from the stencil of the scheme and is well explained in ^[19]. We can calculate $G$ as follows: Let $N_p(z)$ be the Laurent polynomial of the limit stencil and $S_p(z)$ be the Laurent polynomial of the mask of the $S_p$-scheme. Considering $g$ be the highest power of $z$ in the product $N_p(z)S_p(z)$ and let $\hat{N_p(z)} = N_p(z)S_p(z)z^{-g/2}$. We express $\hat{N_p(z)}$ as polynomial $\hat{Q_p(\delta)}$ in $\delta = (1+z)/2\sqrt{z}$. Then the artifact magnitude can be obtained by considering $\sin(\pi \kappa/2)$ for $\delta$ in that polynomial, where $\kappa$ is the spatial frequency leading to

The artifact magnitude in the limit curve produced by the $S_{p}$-scheme for different values of $p$ is shown in Figure 1. Here we observe that as we increase the value of parameter $p$ the number of artifacts decreases. It also decreases with the increase in the number of initial control points.

6.   Numerical examples

In this section, we observe the performance of the proposed $S_{p}$-scheme when we apply it on a different set of control points. First of all, we discuss the performance of the $S_{0}$-scheme. In Figures 2 and 3, we consider two different initial control polygons. These control polygons are represented by dotted lines. We have applied the $S_{0}$-scheme on these initial polygons up to three subdivision levels. Both figures draw a comparison of the curves, which are obtained by the $S_{0}$-scheme, with the initial polygon. Sub-figures 2($a$) and 3($a$) show the curves which are generated by applying the $S_{0}$-scheme one time on these initial polygon. Sub-figures 2($b$) and 3($b$) show the curves which are generated by applying the $S_{0}$-scheme two times on these initial polygons. Finally, we have the limit curves presented in Sub-figures 2($c$) and 3($c$), which are obtained by applying the $S_{0}$-scheme three times on these initial polygons.

Now we discuss the performance of the $S_{1}$-scheme. In Figures 4 and 5, we consider two different initial control polygons. These control polygons are represented by dotted lines. We have applied the $S_{1}$-scheme on these initial polygons up to three subdivision levels. Both figures draw a comparison of the curves, which are obtained by the $S_{1}$-scheme, with the initial polygon. Sub-figures 4($a$) and 5($a$) show the curves which are generated by applying the $S_{1}$-scheme one time on these initial polygon. Sub-figures 4($b$) and 5($b$) show the curves which are generated by applying the $S_{1}$-scheme two times on these initial polygons. Finally, we have the limit curves presented in Sub-figures 4($c$) and 5($c$), which are obtained by applying the $S_{1}$-scheme three times on these initial polygons. Similarly, Figures 6 and 7 show performance of the $S_{2}$-scheme when applied on different polygons.

7.   Conclusions

Subdivision is a primary tool to describe smooth curves and surfaces in the field of computer aided geometric design. There are different variants of LRA in literature and we have presented a new variant of LRA which is based on a shape-preserving SS. We have shown that this shape-preserving variant of LRA generates a family of shape-preserving subdivision schemes. An integer parameter ($p$) is involved in the symbol of the proposed family and by changing this parameter we can get different SSs as per our requirement. We have also presented limit stencil and artifact analysis of the proposed schemes. We have proved that by increasing the value of the parameter $p$, SSs with higher smoothness, higher Hölder continuity, and less artifact magnitude can be obtained. We have also provided several examples to illustrate that the proposed schemes give great flexibility to geometric designers for the generation of smooth geometric models. In conclusion, we have a family of SSs generating high continuity limit curves and satisfying shape-preserving properties as a remarkable addition to the zoo of existing SSs.

Conflict of interest

The authors declare that they have no competing interests.