A closed-form solution to the inverse problem in interpolation by a Bézier-spline curve

A geometric construction of a Bézier curve is presented by a unifiable way from the mentioned literature with some modification. A closed-form solution to the inverse problem in cubic Bézier-spline interpolation will be obtained. Calculations in the given examples are performed by a Maple procedure using this solution.


Preliminaries
Spline and B-spline are concepts that have been studied deeply and comprehensively, since spline functions were realized to be a mathematical tool to model the physical process of drawing a smooth curve in the early 1960s, although the terminology "spline function" was first introduced by Schoenberg in 1947. There is a huge literature of the subject and it would be impossible to have a complete list of research publications even on main developments in the theory of spline functions and their applications. Basic knowledge and important applications of spline functions in approximation and interpolation theory can be found in the two great books, [1] and [6]. The main theoretical results obtained during the development process of this active research area are also provided by the historical notes in these books. Besides, the work [1] introduces algorithms and packages for computation, as well as Fortran codes to implement these algorithms.
Of numerous related topics on spline functions, we are especially interested in the one on interpolation using natural cubic spline functions called "natural spline interpolation", as shortly mentioned in [1,Chapter IV]. We refer the reader to this monograph and many references therein for further reading on this topic. One of the basic reasons to write this paper is that the graph of a natural spline interpolant has the same properties as those of a piecewise Bézier 1 curve that can be constructed easily and naturally in a geometric way.
We recall that in Of modern texts on Bézier curves and related topics, we follow [4, Chapter 15] for a "surprisingly simple geometric construction of a Bézier curve" and just recall here the general results and concepts with some modification. To define a Bézier curve of degree n ≥ 2, we begin with the case of n = 2. Consider distinct points A 0 , A 1 , and A 2 , and their position vectors We will find a curve, called a Bézier curve of degree 2, that begins at A 0 tangent to segment A 0 A 1 at A 0 and ends at A 2 tangent to segment A 1 A 2 at A 2 . Taking an arbitrary value t i of a variable t ∈ [0, 1], we determine points A, B respectively on segments A 0 A 1 , A 1 A 2 , such that: To construct a simple rule, we rewrite (1) in the form: where a and b are the position vectors of A and B, respectively. Then, we define a point C of the curve on segment AB whose position vector c satisfies c = a + t i (b − a). Therefore, we have found the point C, corresponding to the chosen value t i ∈ [0, 1]. Note that if t i = 0, then C ≡ A 0 and if t i = 1, then C ≡ A 2 . By repeating the above steps for other values of the variable t, we have a series of corresponding points Cs of the curve. Thus, the curve can now be seen as the locus of all Cs, replacing t i s with t, and has, noticing (2), its vector function r(t) as follows: for t ∈ [0, 1]. It is obvious that we can write r(t) in the form: We can directly examine if the vector function r(t) satisfies the required properties of a Bézier curve. Indeed, since r(0) = r 0 , r(1) = r 2 , r (0) = 2(r 1 − r 0 ) and r (1) = 2(r 2 − r 1 ), the curve is, respectively, tangent to A 0 A 1 at A 0 and A 1 A 2 at A 2 . The points A 0 , A 1 , A 2 are called control points of the resulting Bézier curve. If r 1 = r 0 or r 2 = r 1 (that is A 1 ≡ A 0 or A 2 ≡ A 1 ), the formula (4) may be still valid, because the zero vector is considered parallel to all vectors by convention.
Next, we consider points A 0 , A 1 , A 2 , and A 3 , corresponding to position vectors To recognize some rule when determining a Bézier curve of degree 3, we denote by a, b, and c the position vectors corresponding to the oriented couples r 0 → r 1 , r 1 → r 2 and r 2 → r 3 . Then, we denote by d and e the position vectors corresponding to the couples a → b and b → c. Finally, we denote by f the position vector of a point on the curve, corresponding to the couple d → e. According to the above steps, a position vector corresponding to an oriented couple is written, for instance, a corresponding to r 0 → r 1 , as a = r 0 + t (r 1 − r 0 ). Similarly, we have the relations: that can be described formally in Table 1.
In view of (6) finally becomes Now, we may have a general definition of Bézier curves. Given n + 1 points A 0 , A 1 , …, A n , a Bézier curve of degree n, taking these points as its control points, is the curve C that can be represented by the vector function defined on [0, 1]: where r i is the corresponding position vector of A i , i = 0, 1, . . . , n. It is easy to see that r(0) = r 0 , r(1) = r n , r (0) = n(r 1 − r 0 ), and r (1) = n(r n − r n−1 ). By these properties, in the increasing direction of t, C begins at point A 0 tangent to segment A 0 A 1 at A 0 , and ends at A n tangent to segment A n−1 A n at A n . Moreover, from its definition, since 0 ≤ t ≤ 1 and C is always contained within the convex hull of all the control points. In the case of n = 3, C is called a cubic Bézier curve and has, from the expanded form of (7), its vector function

Introduction
Another kind of Bézier curve called a cubic and uniform Bézier-spline curve with control points B 0 , B 1 , …, B n (n ≥ 3), beginning at B 0 and ending at B n , can be described as follows: Fig. 1 A part of the curve C with its control points S k−1 , P k−1 , Q k , and S k -For 0 < k < n, divide each segment B k−1 B k into three equal parts with subdivision points P k−1 , Q k , such that, in the direction from B 0 to B n , each B k has Q k and P k as its immediate neighbor to the left and to the right, respectively; denote by S k the midpoint of segment Q k P k and put S 0 = B 0 , S n = B n . -For 0 < k ≤ n, take a cubic Bézier curve C k with control points S k−1 , P k−1 , Q k , and S k represented by a formula like the one in (8). -All C k , k = 1, . . . , n are then joined by the way of constructing them to form the desired curve C that we just call here the Bézier-spline curve, for brevity.
A sample of this curve is illustrated in Fig. 1.
We, respectively, denote by b 0 , b 1 , …, b n the position vectors corresponding to the points B 0 , B 1 , …, B n , and by p k , q k , s k the position vectors corresponding to the points P k , Q k , S k . Then, we have the following: By (9), we get that For 2 ≤ k ≤ n − 1, we find the vector function f k (t) of C k defined on [k − 1, k] by the use of (8) with control points S k−1 , P k−1 , Q k , S k instead, replacing t by t + 1 − k. Hence, from (9) and (10), we can write the following: Similarly, we have found the vector function f 1 (t) of C 1 and f n (t) of C n , respectively defined on [0, 1] and [n − 1, n]. Here, these are functions: Now, the vector function f(t) of C defined on [0, n] can be written as follows: and it is easy to check that Hence, the curvatures of C at t = 0 and t = n are both zero. Actually, the Bézier-spline curve C here can be seen as a piecewise cubic Bézier curve whose vector function has the same meaning of "a spline of degree 3" as given in [5, Definition 6.1.1], or "a natural cubic spline function" as given in [1,Chapter IV]. However, the notions of "spline function", "B-spline function", and their representation in these cited references are deeply investigated by the use of truncated power functions that we will not mention here about. An illustration of C in the case of n = 6 is given in Fig. 2 and depicted by PSTricks 2 macros in the guiding document [7]. Given distinct points S 0 , S 1 , …, S n , we need to find the Bézier-spline curve C with control points B 0 , B 1 , …, B n , such that C interpolates S k , k = 0, . . . , n, in the same meaning as presented above. This is an inverse problem and we reduce it, from (10), to find position vectors b k corresponding to B k , k = 0, . . . , n, as the solution to the system of equations: where s k is position vector of S k , k = 0, . . . , n. Because b 0 = s 0 and b n = s n . we need to find b 1 , …, b n−1 only, and the above system may be rewritten in the form: Since the coefficient matrix A of (12) is both strictly dominant and tridiagonal, we have used Crout factorization (see [2, Algorithm 6.7]), as a traditional way so far, to solve (12) for b k , k = 1, . . . , n − 1. By this way, the control points are obtained iteratively. Therefore, the dependence of each b k on the whole set of s k is just a numeric dependence, although that is a linear and functional one. If the system contains irrational coefficients, it could be very difficult to control the evaluation errors after each step for such an iterative method. This could lead to large errors. In addition, there is no computer algebra system that can provide an explicit expression for the solution depending on the given data set when n is a large number. To our best knowledge, at this moment, there has not been any publication that provides a closed-form solution in general to such a system of linear equations. Actually, we can have global control on both of approximation and evaluation errors with a closed-form solution of the system. Especially, if we used such a solution for our Bézier-spline curve to interpolate data points, for example, in statistics, we would easily give estimates for errors from wrong measurement data. Thus, our main task here is how to solve the system (12) for a closed-form solution and find the most convenient representation of that solution.

A closed-form solution to the inverse problem
According to the factorization algorithm ([2], Algorithm 6.7) as indicated in Sect. 2, the matrix A of (12) is factored into the form LU , where L is lower triangular and U is upper triangular with all 1 on the diagonal as follows: From that algorithm, we easily find Now, we will solve the system (12) by performing elementary row operations to its augmented matrix, which has the last column with the entries 6s 1 − s 0 , 6s 2 , …, 6s n−2 and 6s n−1 − s n . Multiplying row 1 by − −1 1 and adding its result to row 2, then multiplying row 2 of the resulting matrix by − −1 2 , and adding its result to row 3, and so on, we finally obtain the following: ⎡ We will get a formula to evaluate i s by terms of a sequence β m , m = 0, 1, . . . , n − 2. We set β 0 = 1, β 1 = 4, and for k ≥ 1: Then, by induction, we easily check for k ≥ 1 that Moreover, we will need to take values of β k s directly for later use. Noting that (14) is a linear second-order difference equation with constant coefficients, we get its general solution having the form of C 1 r k 1 + C 2 r k 2 , where r 1 and r 2 are the distinct roots of the characteristic equation r 2 − 4r + 1 = 0, as given in [3, Section 2.3]. Since β 0 = 1 and β 1 = 4, the unique particular solution of (14) is as follows: We will use (16) for all k ≥ −1, noticing β −1 = 0. From (15), we also derive the following: Now, for brevity, we set the matrix (13) in the form: so that, we can have, by a backward substitution, its solution: − · · · + (−1) n−2 z n−1 1 2 · · · n−1 or, using (17): We will express b k s in the form of a linear combination of s j s by carefully analyzing z j s. Indeed, from (13) and (17), we may have for 1 ≤ k ≤ n − 2, and Putting (19) and (20) into (18), we have that On the right-hand side, we may add up the two terms in the second round brackets to the first summation sign to write: The term-by-term extraction of the second sum on the right-hand side of (21) can help us find the rule to rewrite this sum in a more convenient form. It is implemented for each index j: By adding up coefficients of terms of the form (−1) i β i−1 s i in "vertical lines", we get the result: In the case of k = 1, we have the following: that is, the first sum on the right-hand side of (22) may be referred to as zero when k = 1. Apparently, if we put (22) into (21), then we need a special summation formula to abbreviate the resulting expression of b k . Here is such a formula that we can check its correctness by induction: Indeed, to prove (23), we consider the special case when s = 1 and examine the correctness of thanks to β 0 = 1. Since is true when m = 1. Assume that (24) is true when m = k, for some k ≥ 1. Then, we have the following: From (16), we derive and it follows from (25) that (24) is true when m = k + 1. Thus, we have proved the correctness of (24) for all m ≥ 1. Now, by (24), we can write the following: noticing that the second sum on the right-hand side of (26) may be referred to as zero when s = 1. By the formula (16), we easily obtain the following: so we have proved (23), taking into account (27). After putting (22) into (21), we use (23) to simplify the new expression of b k , so we have finished the proof of the following theorem. Theorem 3.1 Let C be the Bézier-spline curve that interpolates given distinct points S k in the same meaning as introduced in Sect. 2 with control points B k having the position vectors s k and b k , respectively, k = 0, 1, . . . , n, with b 0 = s 0 and b n = s n . Then, we have that b k , k = 1, . . . , n − 1, is the unique solution of (9) and can be represented by the formula: where β k s is evaluated by (16) that can be now rewritten easily as follows: then the formula for b k in Theorem 3.1 can be rewritten in the components of (w 1 , . . . , w n−1 ) T as follows: If S k = (x k , y k ), k = 0, . . . , n, where y k = f (x k ) is the value of a function f (x) at the value x k of x, then the curve C in Theorem 3.1 gives us an approximation to the graph of f . If n is a large number, C is obviously a better choice beside the graph of the so-called Lagrange form of the interpolating polynomial between S k s.
Finally, it is also hoped that the obtained result might get attention from the need of seeking a simple and effective method for finding curve fitting in statistics or simulating motion orbits in mechanics, as well as for providing a similar treatment to numerically solve differential equations using the finite difference schemata.