PIA and PPIA for Interpolating Points and Derivatives at Endpoints by Bézier Curves

In this study, we are concerned with the interpolation problem that interpolates not only a given set of points but also derivatives at endpoints by using Bézier curves.*e progressive iterative approximation (PIA) is proposed to interpolate the given points and derivatives. To speed up the convergence rate of PIA, we exploit the preconditioned PIA (PPIA), in which the diagonally compensated reduction is used to construct the preconditioner. *e convergence of PIA and PPIA is also analyzed in this study. *e proposed PPIA is applied to approximate higher order or/and rational Bézier curves. Numerical examples are given to illustrate the effectiveness of the proposed methods.


Introduction
We begin with the description of the data interpolation problem with constrains at endpoints. Given a set of organized points p i n i�0 ∈ R 2 or R 3 , r th derivatives d r 0 (r � 1, 2, . . . , u) ∈ R 2 or R 3 at the endpoint p 0 and s th derivatives d s n (s � 1, 2, . . . , v) ∈ R 2 or R 3 at the endpoint p n . For i � 0, 1, . . . , n, the i th point p i assigned an ordered parameter t i , i.e., t 0 < t 1 < · · · < t n . We want to find a m degree Bézier curve, that interpolates these points as well as the derivatives at the endpoints, i.e., C t i � p i , i � 0, 1, . . . , n; d r C(t) dt r | t�t 0 � d r 0 , r � 1, 2, . . . , u; d s C(t) dt s | t�t n � d s n , s � 1, 2, . . . , v.
e terms B m j (t) in (1) are the Bernstein polynomials of degree m, i.e., B m j (t) � t j (1 − t) m− j . According to the existence and uniqueness of polynomial interpolant, since there are n + u + v + 1 different interpolation conditions in (2), we have m � n + u + v.
Data interpolation and interpolation with tangents, curvatures, and other derivatives appear frequently in the fields of science and engineering [1,2]. It is difficult to represent curves with complex shapes by a single curve, and a composite curve composed of several pieces of curves is necessary [3]. Hence, we need to control the smoothness of the composite curve when piecing curves together. It would be great if the derivatives at the endpoints are known. Very often, this knowledge is not given in practice, and we have to determine the derivatives at the endpoints. A careful choice of end conditions is important because it determines the shape of the interpolating curve near the endpoints. ere are several strategies, e.g., natural end condition, Bessel end condition, quadratic end condition [4,5]. In this study, we are especially interested in exploiting iterative methods for solving the interpolation curve (1) and do not consider the determination of constraints at the endpoints.
It is known to all that the interpolation Bézier curve (1) can be obtained by solving a system of linear equations directly when the coefficient matrix is well-conditioned. However, the condition number of the linear system will be very large for large m; therefore, it is required to introduce an efficient algorithm to solve the linear system. In recent years, a geometric iterative method, named PIA, plays an important role in data interpolation. Due to its clear geometric meaning, stable convergence, and simple iterative scheme, PIA has intrigued the researchers for decades. For more details about PIA, readers can refer to a recent survey study [6], in which the authors summarized PIAs and their successful applications.
In this study, we exploit the PIA format that interpolates not only a given set of points but also derivatives at endpoints by using Bézier curves. e study is organized as follows. In Section 2, we propose the PIA format and then analyze its convergence. In order to accelerate the convergence rate of PIA, we exploit the preconditioning technique for PIA in Section 3. In Section 4, the proposed method is applied to approximate higher order or rational Bézier curves with constraints. Numerical examples are given to illustrate the effectiveness of our proposed methods in Section 5, and some conclusion remarks are given in the last section.

Interpolating Derivatives Conditions at Endpoints.
We note that the k th derivative of (1) is given by where Δ k q j is defined according to the recurrence relative formula: By direct calculation, from (2) and (3), we have Hence, the control points q j (j � 0, 1, . . . , u, m− v, . . . , m) in (1) can be determined according to (5). As an example, we outline the control points with respect to the case of u � v � 3.
erefore, the Bézier curve (1) can be expressed by where the control points q j (j � 0, 1, . . . , u, m − v, . . . , m) can be obtained according to (5), and the control points q j (j � u + 1, u + 2, . . . , m − v − 1) need to be determined. We will give a computing method for solving these points in the following subsection.

Approximate Interpolation Algorithm.
First, the given points p i (i � 1, 2, . . . , n − 1) and q j (j � 0, 1, . . . , u, m− v, . . . , m) are interpreted as the control points of a Bézier curve, and therefore, we can generate an initial approximate interpolation curve: be the adjusting vector of the j th control point. en, we can update the j th point q 1 j � q 0 j + δ 0 j and generate the first approximate interpolation curve: 2 Mathematical Problems in Engineering We repeat these procedures; thus, we obtain a sequence of approximate interpolation curves: where us, we obtain a sequence of approximate interpolation curves C k (t), k � 0, 1, . . . ,. e initial Bézier curve is said to have the property of PIA with constraints at endpoints if lim k⟶∞ C k t i � p i , i � 0, 1, . . . , n; Remark 1. By simple calculation, it is easy to verify that these approximate interpolation curves C k (t) interpolate the endpoints p 0 and p n as well as the derivatives d r 0 (r � 1, 2, . . . , u) and d s n (s � 1, 2, . . . , v). Let en, the iterative process (11) can be written in the matrix form: where I is the (n − 1) × (n − 1) identity matrix and B is said to be the collocation matrix resulting from the Bernstein basis at parameters t i e iterative process (14) is equivalent to the Richardson iteration for solving the system 2.3. Convergence Analysis of PIA. Before analyzing the convergence, we review some definitions and conclusions.
Definition 1 (see [7,8]). A matrix is called totally positive if all its minors are positive and positive stable (semistable) if all its eigenvalues have positive (nonnegative) real parts. It is easy to verify that a positive stable matrix is nonsingular. We note that all the eigenvalues of a totally positive matrix are positive ( [9]); thus, a totally positive matrix is also nonsingular.
Definition 2 (see [7]). A basis b i (t) n i�0 is called normalized if b i (t) ≥ 0(i � 0, . . . , n) and n i�0 b i (t) � 1 and totally positive if its collocation matrix at any increasing sequence is a totally positive matrix. Lemma 1 (see [8]). Let A be a given matrix. en, A is positive stable if and only if for each nonzero vector x, there exists a positive definite matrix K, such that Re(x * KAx) > 0.
Lemma 2 (see [10]). Let A be an n × n matrix with entries a ij . For i ∈ 1, . . . , n { }, let Λ i � j≠i |a ij | be the sum of the absolute values of the nondiagonal entries in the i th row. Let D(a ii , Λ i )⊆C be a closed disc centered at a ii with radius Λ i . en, every eigenvalue of A lies within at least one of the discs D(a ii , Λ i ).

Theorem 1.
e Bézier curve has the property of PIA with constraints at endpoints.
Proof. As stated in Remark 1, the approximate interpolation curves C k (t) interpolate endpoints p 0 and p 1 and derivatives d r 0 (r � 1, 2, . . . , u) and d s n (s � 1, 2, . . . , v). On the other hand, we insert u + v different parameters into the ordered parameters set t j n j�0 ; hence, we obtain a sequence of increasing parameters: Note that the Bernstein basis is normalized and totally positive [7]. Let B be the collocation matrix resulting from the Bernstein basis at the parameter sequence (17). en, it follows from Definition 2 that B is totally positive. Since the matrix B defined in (15) is a submatrix of B, then B is also totally positive according to Definition 1. Hence, all the eigenvalues of B are positive.
Again, since the Bernstein basis is normalized, we have ‖B‖ ∞ � 1; thus, the infinite norm of the submatrix B is less is implies that the spectral radius of the iteration matrix I − B is less than 1, i.e., ρ(I − B) < 1. erefore, the iterative format (14) converges. e proof is thus complete. eorem 1 tells us that the iteration process (14) always converges and we can obtain a satisfying approximate interpolant without solving a system of linear equations. However, as stated in [11], the collocation matrix resulting from the Bernstein basis is usually ill-conditioned, and the convergence rate of PIA is very slow. erefore, it is necessary to introduce a preconditioner to speed up the convergence rate of PIA. In the following section, we propose a preconditioning technique for PIA by exploiting the properties of collocation matrix B.

Construction of Preconditioner.
To derive a preconditioner for PIA, we first split B as and B q � B − R is a band matrix with bandwidth 2q + 1 (0 ≤ q ≤ n − 2). en, we define a diagonal matrix such that De � Re with e � (1, 1, . . . , 1) T . Finally, let be the diagonally compensated reduction matrix of B. en, the matrix R is called the reduced matrix and D is called the compensation matrix for R [10,12]. If M q is invertible, it can serve as a preconditioner for (16) and yield the preconditioned system 3.2. Preconditioned PIA. By applying Richardson method to the preconditioned system (21), we obtain We refer the iterative process (22) as PPIA.

Convergence Analysis of PPIA
Proof. In the proof of eorem 1, we show that B is positive stable. By Lemma 1, for any nontrivial y, there exists a positive definite matrix K, such that Re[y * KBy] > 0. According to the definitions of D and R, both the diagonal entries of D − R and the sum of the absolute values of the nondiagonal entries in the j th row of is means that all the eigenvalues of D − R have nonnegative real parts and D − R is positive semistable. Hence, for the positive definite matrix K and the nontrivial x, we have Re It follows from Lemma 1 that M q is positive stable and there is no eigenvalue equal to zero. is completes the proof.
As it is stated in [13], the error matrix E � M q − B or the residual matrix R � I − M − 1 q B is often used to measure the efficiency of preconditioner M q .
From (18) and (19), we have As is stated in [12], us, the sum of the entries of each column outside of the q diagonals decays to 0. is means that ‖E‖ ∞ approaches to 0 as q increases. Combined with (23), we conclude that the matrix M q is a good approximation to B for large q. us, the spectrum of M − 1 q B is clustered around 1 for large q. □ Remark 2. As it is known to all, the iterative process (22) is convergent if the spectral radius of the iteration matrix is less than 1, that is, In particular, when the ei- . erefore, the convergence of (22)  Recalling Note that when According to (25) and the conjecture in Remark 2, we have us, PPIA converges. also tends to 0. In particular, ρ(I − M − 1 q B) equals to 0 when q � n − 2. Hence, the spectral radius ρ(I − M − 1 q B) may decrease as q increases, i.e., the bigger the value of q is, the faster the PPIA format converges.
On the other hand, it should be pointed out that at the k(k � 0, 1, . . . , ) th iteration of PPIA, we have to calculate M − 1 q Δ (k) or solve the equation M q X (k) � Δ (k) equivalently. It is costly, especially for large q. erefore, we need to seek a tradeoff between the convergence rate and computational complexity. Experimentally, we found that q ≈ (n/2) is a suitable choice.
As a sample-based polynomial approximate method, Lu employed the weighted PIA (WPIA) to iteratively approximate the points sampled from the rational Bézier curves and achieved good approximations [17]. As mentioned in [19], Lu's method is slow due to the slow convergence of WPIA and cannot deal with polynomial approximation with higher order continuity conditions at the endpoints. In fact, polynomial approximations with constraints are more practical than those without constraints. Due to the effectiveness in interpolation, these problems resulted from Lu's method can be solved efficiently by employing the presented PPIA with constraints.

Polynomial Iterative Approximation for Higher Order and
Rational Bézier Curves with Constraints. Given a rational Bézier curve of degree N where v j and ω j ∈ R + are the control points and the associated weights, respectively. First, we sample n + 1 points R(t j ) n j�0 at the parameter values t j (j � 0, 1, . . . , n) and the derivatives R (r) (t 0 ) u r�1 and R (s) (t 1 ) v s�1 at endpoints, where 0 � t 0 < t 1 < · · · < t n � 1. en, PPIA for Bézier curves with constraints is employed to interpolate the points R(t j ) n j�0 and the derivatives R (r) (t 0 ) u r�1 and R (s) (t 1 ) v s�1 . erefore, we can generate a sequence of m degree Bézier curves C k (t), which are the polynomial approximations of R(t).

Remark 4.
We remark here that if all the weights ω j in (18) equal to 1, the rational Bézier curve (28) will degenerate into the Bézier curve of degree N. At this time, if m < N, it is the problem of approximating a given curve by a lower degree Mathematical Problems in Engineering Bézier curve, which is an important approximation problem in CAGD, i.e., degree reduction of Bézier curves.
As stated in [17], the iterative method guarantees the convergence of data interpolation but does not necessary converge to the given curve R(t). erefore, it is not necessary for PPIA to iterate infinitely. To ensure the computing efficiency, the iteration can be terminated if is satisfied, where L (k) p is the approximation error at the k th iteration given by Since there may exist rational functions in the integrand of (30), we need an effective and stable numerical method to calculate the integral (30). We note in [29] that the Gaussian rule or/and the adaptive quadrature method would be a better choice for rational cases. erefore, we employ the τ-order Gauss-Legendre rule to calculate the integral (30).
By solving the roots of the Legendre polynomial of degree τ, we can obtain Gauss points x i , i � 1, 2, . . . , and then calculate τ positive weights ω i . Finally, the integral in (30) can be evaluated approximately by We refer the reader to read [30] for more details about the τ-order Gauss-Legendre quadrature. In our tests, we used 15-order Gauss-Legendre quadrature to estimate the approximation error L (k) 2 , i.e., τ � 15 and p � 2. Finally, we summarize the approximation algorithm in the following Algorithm 1.

Numerical Examples
In this section, several numerical experiments are given to illustrate the effectiveness of the proposed methods. All the numerical experiments were carried out on a computer with Intel(R) Core(TM) i5-5200U CPU @2.20 GHz by Matlab R2012b.

Tests of PIA and PPIA with Constraints.
Since we have shown in eorem 1 that the approximate interpolation curves C k (t) interpolate the endpoints and the derivatives at the endpoints exactly, we need to measure the interpolation error at the points p i , i � 1, 2, . . . , n − 1. erefore, the interpolation error of the k th approximate interpolation curve C k (t) can be expressed by where the norm is the Euclidean norm.

(33)
Example 2. We sample 19 points p i (i � 0, 1, . . . , 18), r(r � 1, 2, 3) th derivatives d r 0 , and s(s � 1, 2, 3) th derivatives d s n from the helix of radius 5 in the following way: PIA and PPIA formats for constraints at the endpoints are employed to test Example 1 and Example 2. When we test PPIA, we take the half bandwidth q � 4 in Example 1 and q � 10 in Example 2. In Table 1, we list the spectral radii of the iteration matrices of PIA and PPIA with different constraints at the endpoints. For one thing, the spectral radii of the iteration matrices of PIA and PPIA are less than 1; this means that both PIA and PPIA formats converge. For another, the spectral radii of PPIA are much less that those of PIA; hence, we can expect that PPIA would have a better convergence behavior than PIA.
We list in Table 2 the interpolation errors of PIA and PPIA when we test Example 1, and we list in Table 3 the Input: original curve R(t), largest admissible number k max , and outer stopping tolerance θ in (19) Output: control points of the approximate curve, number of iterations k, and approximation error L (k+1) p (1) Sample the points R(t j ) n j�0 and the derivatives R (r) (2) Compute the initial error L (0) 2 according to (31) (3) For k � 0, 1, . . . , k max (i) Employ PPIA for Bézier curves with constraints to generate the (k + 1) th approximation curve 2 , break for End for ALGORITHM 1: Polynomial approximation for higher order and rational Bézier curves with constraints. interpolation errors of PIA and PPIA when we test Example 2. As we can see from Tables 2 and 3, the interpolation errors of PIA and WPIA decrease as the number of iterations increases, which indicates that both PIA and PPIA converge. With the same iterations, the interpolate errors of PPIA are much less than those of PIA.
is means that the preconditioning technique can accelerate the convergence rate of PIA significantly.
In Figures 1 and 2, we display the interpolation Bézier curves (dashed lines) when approximating the points with different derivatives given in Example 1. We also show in Figures 1 and 2 the semicircle (solid lines) for comparison. It should be pointed out that the approximate interpolation Bézier curves generated by PPIA almost coincide with the semicircle, no matter which endpoint conditions are used. In Figures 3 and 4, we display the interpolation Bézier curves when approximating the points combined with different derivatives given in Example 2. All the approximate interpolation Bézier curves in Figures 1-4 are obtained by preforming 5 PIA or PPIA iterations. Clearly, the IPPIA can yield better approximations than PIA.

Tests of Polynomial Approximation for Higher Order and
Rational Bézier Curves. First, we employ Algorithm 1 to test the degree reduction of higher order Bézier curves in Example 3 and then test the polynomial approximation for rational Bézier curve in Example 4 and Example 5. In our tests, we take the endpoint interpolation conditions (u � v � 1) and take the half bandwidth q � ⌈n/2⌉ when we construct the preconditioner for PIA, where ⌈ · ⌉ is the round-up operation. When we test Example 3 and Example 4, we take θ � 0.9. Again, when we test Example 5, we take θ � 0.98.      By employing Algorithm 1, we use the Bézier curves of degree m(m � 6, 7, . . . , 14) to approximate the Bézier curve given in Example 3. In Table 4, we list the iteration number k, the approximation errors L (k) 2 , and the elapsed CPU time T (in seconds) of Algorithm 1; the numerical results run by Lu's method ( [17]) are also listed for comparison.
e results for degree reduction with constraints, reported in Table 4, indicate that Algorithm 1 outperforms Lu's method, in terms of the approximation error and the computing time.    In Figure 5(a), we show the approximations of degree 8 obtained by Algorithm 1 (blue dashed line) and Lu's method (red dashdot line). In Figure 5(b), we display the corresponding error distance curves of the approximations shown in Figure 5(a). It is evident from Figure 5 that Algorithm 1 can produce better approximations than Lu's method. What is more, the approximations obtained by Algorithm 1 are C 1 -continuity at two endpoints, whereas the approximations obtained by Lu's method only satisfy the simplest C 0 -continuity.        (11,8), and the associated weights 1, 2, 3, 6, 4, 5, 3, 4, 2, 1.
We employ Algorithm 1 and Lu's method to generate polynomial approximations for the rational Bézier curves given in Example 4 and Example 5. In Table 5, we list the iteration number k, the approximation error L (k) 2 , and the elapsed CPU time T (in seconds) of Algorithm 1 and Lu's method. e results indicate that Algorithm 1 also performs better than Lu's method for polynomial approximation for rational Bézier curves with constraints. e m(m � 9, 12, 15) degree approximations for the rational Bézier curve in Example 4 and their corresponding error distance curves are shown in Figures 6 and 7, respectively. Again, m(m � 10, 12, 15) degree approximations for the rational Bézier curve in Example 5 and their corresponding error distance curves are shown in Figures 8  and 9, respectively. It is already clear from these figures the reliability of the proposed approximation algorithm.

Conclusion
e presented PIA format for Bézier curves iteratively interpolates not only a set of points but also derivatives at endpoints. It can construct higher order continuous piecewise interpolation Bézier curves without solving a linear system directly. To accelerate the convergence rate of PIA, the preconditioning technique is introduced. Numerical examples demonstrate that the given methods are effective, and the convergence rate of PIA is accelerated significantly by the presented preconditioning technique.
Furthermore, due to the efficient performance in data interpolation, the presented PPIA for Bézier curves with  constraints is applied to approximate higher order and rational Bézier curves. Our numerical experiments show that the proposed method achieves higher efficiency and results in a better approximation than the similar sample-based polynomial approximation method. More importantly, they can satisfy C k (k � 1, 2, . . . , ) continuity at the endpoints.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.