Spline Approximation, Part 2: From Polynomials in the Monomial Basis to B-splines—A Derivation

: In a series of three articles, spline approximation is presented from a geodetic point of view. In part 1, an introduction to spline approximation of 2D curves was given and the basic methodology of spline approximation was demonstrated using splines constructed from ordinary polynomials. In this article (part 2), the notion of B-spline is explained by means of the transition from a representation of a polynomial in the monomial basis (ordinary polynomial) to the Lagrangian form, and from it to the Bernstein form, which ﬁnally yields the B-spline representation. Moreover, the direct relation between the B-spline parameters and the parameters of a polynomial in the monomial basis is derived. The numerical stability of the spline approximation approaches discussed in part 1 and in this paper, as well as the potential of splines in deformation detection, will be investigated on numerical examples in the forthcoming part 3.


Introduction
In engineering geodesy, the use of point clouds derived from areal measurement methods, such as terrestrial laser scanning or photogrammetry, results in the necessity to approximate them by a curve or surface that can be described using a continuous mathematical function, often by means of splines. In part 1 of a three-part series of articles, presented by Ezhov et al. [1], the basic methodology of spline approximation was shown by means of ordinary cubic polynomials concatenated under constraints for continuity, smoothness, and continuous curvature. The resulting linear adjustment problem can be solved within the Gauss-Markov model with constraints for the unknowns.
However, a starting point for advanced considerations in engineering geodesy are almost always the formulas for B-spline curves and B-spline surfaces given in the textbook by Piegl and Tiller [2] (pp. 81 and 100), where the functional values of the B-spline basis functions are recursively computed according to the formulas by de Boor [3] and Cox [4]. As these formulas have a very complex mathematical derivation, but are still very easy to use, they are mostly used like a given recipe. Attempts to explain B-splines in an illustrative way often only include explanations of the application of the de Boor's algorithm, as, for example, in the contribution by Lowther and Shene [5]. A derivation of the B-spline with the help of quadratic Bézier spline curves was presented by Berkhahn [6] (p. 73 ff.).
In this paper, we take the spline representation presented by Ezhov et al. [1] as a starting point to derive the B-spline form by means of basis changes. This is to show the user that the de Boor's algorithm can be interpreted as another, very elegant and numerically stable, representation of the simple and intuitive representation based on ordinary polynomials. is fulfilled. The choice of a suitable function f depends on the properties of the data, e.g., its smoothness or periodicity. Furthermore, the choice of f is influenced by computational aspects, i.e., the computational effort for the determination of the coefficients and numerical stability of the resulting equation system. Functions that are commonly used for interpolation are: Ezhov et al. [1] elaborated a spline approximation using piecewise ordinary cubic polynomials. In the following, we will take a closer look at the interpolation using polynomials in different representations to bridge the gap to the B-spline representation.
In general, the set of functions for interpolating given data points consists of coefficients a i and so-called basis functions φ 0 (x), . . . , φ n (x) and the interpolating function is defined as a linear combination of these basis functions. Considering the interpolation condition (1), we obtain The simplest and most common type of interpolation uses polynomials. These polynomials can be represented in different ways, as shown by Gander [31], but all of them will give the same result for the interpolation.
To avoid excessive formula derivations and to illustrate the geometric relationships in the formula derivations, quadratic splines that consist of piecewise parabola segments are considered in the following. The generalization to splines of a higher degree is relatively straightforward.
Using the example of a parabola, we start with the representation of a polynomial in the monomial basis (ordinary polynomial) and perform the following transitions: (i) From ordinary polynomial to Lagrangian form (Section 2); (ii) From Lagrangian form to Bernstein form (Section 3); (iii) From Bernstein form to B-spline representation (Section 4).
The reason for the representation of a polynomial in a different basis than the wellknown monomial basis applied by Ezhov et al. [1] often lies in the fact that for the computation of the polynomial coefficients equation systems result, whose solution requires less computational effort and is numerically more stable. However, it is important to understand that the change of the basis functions still results in the same interpolating polynomial for the given data, and only the representation of the polynomial changes. After the derivation of the B-spline representation, this form is used for spline approximation using least squares adjustment (Section 5). The interesting fact that the determined spline parameters can be used for a transition "backwards" from B-spline to ordinary polynomial is shown in Appendix A. The transition "forwards" from ordinary polynomial to B-spline is shown in Appendix B.

Transition from Ordinary Polynomial to Lagrangian Form
With the monomial basis functions φ i (x) = x i , i = 0, 1, . . . , n, we obtain for the parabola with n = 2 from (3) A polynomial in the monomial basis is often referred to as standard form or ordinary polynomial. Arranging the monomials in a vector m(x) = [1, x, x 2 ] T and the coefficients in a vector a = [a 0 , a 1 , a 2 ] T , we obtain A visualization of the monomial functions for the parabola is shown in Figure 1.
polynomial to B-spline is shown in Appendix B.

Transition from Ordinary Polynomial to Lagrangian Form
With the monomial basis functions ( ) , 0,1, , We obtain for the parabola with 2 n = from Error! Reference source not found.
A polynomial in the monomial basis is often referred to as standard form or ordinary polynomial. Arranging the monomials in a vector 2 T ( ) [1, , ] x x x = m and the coefficients in a vector A visualization of the monomial functions for the parabola is shown in Error! Reference source not found.. Considering the interpolation condition Error! Reference source not found., the coefficients can be determined from the solution of the linear equation system with T 0 1 2 [ , , ] y y y = y and the Vandermonde matrix  Considering the interpolation condition (1), the coefficients can be determined from the solution of the linear equation system with y = [y 0 , y 1 , y 2 ] T and the Vandermonde matrix cf. the textbook by Farin [7] (p. 100). Using a numerical example with three points that define a parabola from the handbook by Bronshtein et al. [32] (p. 918), see Table 1, the determination of the coefficients will be demonstrated. From Equation (7), we obtain for this numerical example Mathematics 2021, 9, 2198 5 of 24 This system of equations can be solved with, for example, the Gaussian elimination method and the numerical result is a 0 = 1, a 1 = 17 6 , a 2 = − 5 6 . Thus, the equation of the parabola reads It can be stated that, when using the monomial basis, the numerical effort to determine the coefficients requires roughly n 3 operations using direct methods for solving systems of linear equations, as stated e.g., by Farin [7] (p. 102), and is therefore comparably high. It can be reduced by employing iterative methods. In addition, the matrix of the equation system is increasingly ill-conditioned as the degree of the polynomial increases. Both the conditioning of the resulting linear equation system and the computational effort for determining the coefficients can be improved by using other basis functions. The result is the same interpolating polynomial, but in a different representation.
In the following subsections we will illustrate two ways how to perform a transition from the monomial basis to a representation with Lagrange polynomials.

Transition by Means of Basis Transformation
We want to solve the problem of basis transformation between two different sets of coefficients a i and b i fulfilling where φ i (x) and ψ i (x) are two different sets of basis functions. Using the monomial basis function (4) in the first summation formula, we obtain as already shown in (5) and (6). Multiplying out the second summation formula in (11) yields and, according to (11), we can write In the case of monomial basis, we selected basis functions and obtained the coefficients, see (7), while here we do the opposite. We select b 0 = y 1 , b 1 = y 2 , b 2 = y 3 for the coefficients to obtain corresponding basis functions for which we introduce the new notation l i (x) with i = 0, 1, 2. Furthermore, for a 0 , a 1 , a 2 we introduce the solution according to (7) and we obtain By comparing the coefficients of the vectors [ y 0 y 1 y 2 ], we see that the basis functions l i (x) can be computed from Introducing the vector l( The same formula can also be obtained taking the basis transformation from Lagrange to monomials V T l(x) = m(x), as presented by Gander [31], and solving it for l(x).
From (17), we obtain the result which is called the Lagrange basis. A visualization of the Lagrange basis functions for the parabola is shown in Figure 2.
By comparing the coefficients of the vectors 0 Introducing the vector The same formula can also be obtained taking the basis transformation from Lagrange to monomials m , as presented by Gander [31. ], and solving it for ( ) x l . From Error! Reference source not found., we obtain the result T 0 2 0 1 1 2 which is called the Lagrange basis. A visualization of the Lagrange basis functions for the parabola is shown in Error! Reference source not found.. Finally, the resulting polynomial reads An alternative derivation of the Lagrange representation of a polynomial is presented in the following subsection. Finally, the resulting polynomial reads An alternative derivation of the Lagrange representation of a polynomial is presented in the following subsection.

Transition by Linear Combinations
In this subsection, we present an alternative derivation of the polynomial in the Lagrange representation. We use the fact that each parabola can be represented as a combination of two lines. In fact, if we take a look at (5), we find that it can be rearranged into the form f (x) = a 0 + (a 1 + a 2 x)x, which represents a combination (product) of two lines plus some constant value a 0 . By using this idea, the Lagrange representation of a quadratic polynomial can be derived, which results in a linear combination of three quadratic functions p(x) and the respective y i coordinates.
We consider the points (x 0 , y 0 ), (x 1 , y 1 ), (x 2 , y 2 ), where x 0 < x 1 < x 2 . Now two lines within the interval [x 0 , x 1 ) and within the interval [x 1 , x 2 ] are defined, see Figure 3. The subscripts i and d in f i,d represent the number of the polynomial i for a given degree d.

nates.
We consider the points 0 0  From two given points, the slope of the line Error! Reference source not found. can be written as The y-intercept can be computed from Inserting Error! Reference source not found. and Error! Reference source not found. into Error! Reference source not found. yields, after some rearrangement, for the first line Analogously, the equation for the second line can be written as From two given points, the slope of the line (20) can be written as The y-intercept can be computed from Inserting (22) and (23) into (20) yields, after some rearrangement, for the first line Analogously, the equation for the second line can be written as By taking a combination of these two line equations, where x varies within the interval [x 0 , x 2 ], the expression for the parabola is derived, which represents a linear combination of three quadratic polynomial functions and the coordinates y 0 , y 1 and y 2 .
To explain how the factors in front of the terms f 0,1 (x) and f 1,1 (x) in (26) are derived, the idea behind the Lagrange polynomials has to be explained. Lagrange was looking for an interpolation polynomial that could be constructed without solving a system of equations, e.g., with Vandermonde matrix as design matrix, see (7). There is a widely accepted assumption that his idea for the solution was to find a function that, at each given data point (x 0 , y 0 ), (x 1 , y 1 ), . . . , (x n , y n ), gives 1 and, at the rest of the other given points gives 0. For a mathematician, this type of polynomials was, presumably, relatively easy to find. We will take the case of a parabola into consideration. Therefore, for each given x i , the following polynomials can be created. As previously mentioned, these functions p x i (x) are equal to 1 when the variable x is equal to the corresponding coordinate x i of the given point and 0 at all other given points. Therefore, by multiplying each of these functions p x i (x) with their corresponding values y i ensures that, at the particular given point x i , the result will be y i and 0 at all other given points. The summation of all these functions p x i (x), multiplied by their corresponding values y i , results in a Lagrangian interpolation polynomial The quadratic polynomial (30) represents a linear combination of the Lagrangian basis polynomials with the constants y 0 , y 1 and y 2 as coefficients. A general formula for the definition of the Lagrange interpolation polynomials can be found e.g., in the handbook by Bronshtein et al. [32] (p. 918).
As previously explained, a quadratic polynomial can be represented as a combination of two lines. Using the Lagrangian form for a polynomial of first order, the following two line equations are derived. If we follow the already described Lagrange's idea to obtain a quadratic polynomial, the following multiplications must be performed. Consequently, Equation (33) results in a Lagrangian polynomial of second order. Looking at (33), the factors ( (26). However, if we compare the Lagrange polynomials in (31), (32), (33) with the expressions in (24), (25), (26), the only difference is that (31), (32) and (33) have negative denominators and numerators at some points in the Lagrange polynomials, e.g. (25) and (26) there are no negative denominators or numerators, although the equations are the same as in (31), (32) and (33). This is for reasons of better analogy with the B-spline form.
After the origin of the factors in front of f 0,1 (x) and f 1, (26) is explained, inserting (24) and (25) into (26) finally yields This form coincides with (19) and represents a second degree Lagrange interpolating polynomial through three points, see the handbook by Bronshtein et al. [32] (p. 918). The plot of the parabola is depicted in Figure 4.
This form coincides with Error! Reference source not found. and represents a second degree Lagrange interpolating polynomial through three points, see the handbook by Bronshtein et al. [32. ] (p. 918). The plot of the parabola is depicted in Error! Reference source not found.. Inserting the values of the numerical example from Error! Reference source not found. into Error! Reference source not found. yields the Lagrangian form which can be "simplified" to the monomial basis form Inserting the values of the numerical example from Table 1 into (34) yields the Lagrangian form which can be "simplified" to the monomial basis form The result is the same as in (10), but the values of the parameters a 0 , a 1 , a 2 are now derived with the advantage of not explicitly solving a linear equation system. A disadvantage of Lagrange interpolation in practical application is directly visible from (19), resp. (34), namely that each time a value x changes, the Lagrange basis polynomials must be recalculated. Further general limits of Lagrange interpolation, e.g., that polynomial interpolants may oscillate, are explained by Farin [7] (p. 101 ff.).

Transition from Lagrangian Form to Bernstein Form
The Bernstein form, named after S. N. Bernstein, was used to prove the Weierstrass approximation theorem; see the explanation by Farin [7] (p. 90 ff.). The coefficients of a Bernstein polynomial over the interval [0, 1] are all non-negative and sum up to 1, as Farin [7] (p. 57 ff.) showed, which is called a convex combination; see the textbook by Rockafellar [33] (p. 11). Therefore, by using the Bernstein form, numerical instabilities are avoided.
Looking at Figure 3, the interpolating parts of the line segments can be extended to both ends of the interval [x 0 , x 2 ] in which the polynomial is defined. At the end of the extended line segments, two additional y-coordinates for the values x 0 and x 2 are computed, as seen in Figure 5. To distinguish these two y-coordinates from the given values, a different notation, y 0 and y 2 , is used (not to be confused with the derivation of a function). In a general case, all five points have different coordinates which implies that the coordinate x 1 is not necessarily in the middle of the interval [x 0 , x 2 ]. Therefore, in general values 0 x and 2 x are computed, as seen in Error! Reference source not found.. To distinguish these two y-coordinates from the given values, a different notation, 0 y′ and 2 y′ , is used (not to be confused with the derivation of a function). In a general case, all five points have different coordinates which implies that the coordinate 1 x is not necessarily in the middle of the interval 0 2 [ , ] x x . Therefore, in general 1 From Error! Reference source not found., as explained previously, the line equations are derived. By taking combinations of these two lines, it follows which can be brought into the form From Figure 5, as explained previously, the line equations are derived. By taking combinations of these two lines, it follows which can be brought into the form A geometrical interpretation of the term (y 2 + y 0 )/2 is depicted in Figure 6.
is depicted in Error! Reference source not found..
and Error! Reference source not found. obtains the form in which and 2 x are the Bernstein basis polynomials and 0 y , and 2 y are the Bernstein coefficients. In general, a polynomial in Bernstein form of degree n can be written as Equation (39) is a more general representation of the Bernstein form, independent of the limits of the interval [x 0 , x 2 ], which means that x can vary between any two real values. For proving the Weierstrass approximation theorem, Bernstein used a special interval of length 1, i.e., [x 0 , x 2 ] = [0, 1]. Using this interval simplifies the equations, which were later used by Bézier and de Casteljau. Therefore, when using [x 0 , and (39) obtains the form in which (1 − x) 2 , 2x(1 − x) and x 2 are the Bernstein basis polynomials and y 0 , (y 2 + y 0 )/2 and y 2 are the Bernstein coefficients. In general, a polynomial in Bernstein form of degree n can be written as where β i are the Bernstein coefficients and are the Bernstein basis polynomials; see e.g., the handbook by Bronshtein et al. [32] (p. 935).
For a parabola, we obtain B 0,2 = (1 − x) 2 , B 1,2 = 2x(1 − x) and B 2,2 = x 2 , cf. (41). These basis polynomials for a parabola are depicted in Figure 7. In Error! Reference source not found., the term 2 0 ( )/2 y y ′ ′ + , see Error! Reference source not found., reveals one fundamental property of the quadratic polynomial that, to the best of our knowledge, is not described in the literature so far. It can be stated as: All pairs of line segments that have one end fixed at the beginning or at the end of the parabola and the other end at the opposite sides of the interval, and intersect themselves on the parabola, have an equal average of the coordinates at the ends that are not fixed. Moreover, the extremum of the parabola lies on the intersection of those line segments that have equal values at the opposite ends of the interval.  In (39), the term (y 2 + y 0 )/2, see Figure 6, reveals one fundamental property of the quadratic polynomial that, to the best of our knowledge, is not described in the literature so far. It can be stated as: All pairs of line segments that have one end fixed at the beginning or at the end of the parabola and the other end at the opposite sides of the interval, and intersect themselves on the parabola, have an equal average of the coordinates at the ends that are not fixed. Moreover, the extremum of the parabola lies on the intersection of those line segments that have equal values at the opposite ends of the interval.
The proof of this statement is straightforward. A parabola is uniquely defined by three points. If we retain end points (x 0 , y 0 ) and (x 2 , y 2 ), replacing the point (x 1 , y 1 ) by any other point on the same parabola, we will again obtain the equation of the form (39). The equation of the parabola remains unchanged although y 0 and y 2 have changed, since the only term depending on the new point is (y 0 + y 2 )/2. A parabola is described by a differentiable function. By computing its derivative, equating it with zero and solving the resulting equation for x, we obtain the coordinate of the extremum x e . Consequently, the resulting equation for the coordinate of the extremum y e can be solved. By substituting these two values for x e and y e in a line equation constructed for y 0 and y 2 , where x = x 0 and x = x 2 respectively, it can be shown that y 0 = y 2 = C y . Therefore, y 0 = y 2 implies C y = (y 0 + y 2 )/2. This proves the above statement, which is illustrated in Figure 8. . The convex hull formed by the control points 0 0 ( , ) x y , ( , ) x y C C and x y , between M and N , that defines the parabola, cf. Error! Reference source not found.. All of these pairs of line segments, as previously explained, have an equal average of the y-coordinates at the end points. Therefore, the term that uniquely represents all these averages is written more generally in the form where m y′ and n y′ represent any pair of y′ -coordinates from the mentioned line segments that intersect on the parabola. Correspondingly, the coordinate x C is defined as where i x and In the context of geometric modelling, the point defined by ( , ) x y C C is called control point and, concerning basis splines (B-splines), it is named de-Boor-point. In Error! Reference source not found., this point is marked with a blue dot. Moreover, the first and the last point of the parabola are also referred to as control points, see e.g., the textbook by Piegl and Tiller [2. ] (pp. 389-390).
With Error! Reference source not found., the general form of Error! Reference source not found. can be written as Considering the above, there are an infinite number of combinations defining the term (y 2 + y 0 )/2 depending on the position of the point (x 1 , y 1 ), between M and N, that defines the parabola, cf. Figure 6. All of these pairs of line segments, as previously explained, have an equal average of the y-coordinates at the end points. Therefore, the term that uniquely represents all these averages is written more generally in the form where y m and y n represent any pair of y -coordinates from the mentioned line segments that intersect on the parabola. Correspondingly, the coordinate C x is defined as where x i and x i+1 define the ends of the interval in which the parabola varies. Hence, C x is always in the middle of the interval. There is an additional property of the line segments Mn 0 and Nm 6 . They intersect at the point (C x , C y ) and are tangential to the curve. This can be easily proven. We can compute the first derivative and solve for y m and y n from (44), then construct two line equations for (x 0 , y 0 )(x 2 , y n ) and (x 0 , y m )(x 2 , y 2 ) and solve for their intersection.
In the context of geometric modelling, the point defined by (C x , C y ) is called control point and, concerning basis splines (B-splines), it is named de-Boor-point. In Figure 8, this point is marked with a blue dot. Moreover, the first and the last point of the parabola are also referred to as control points, see e.g., the textbook by Piegl and Tiller [2] (pp. 389-390).
With (44), the general form of (39) can be written as As a more general representation of the Bernstein polynomial, (46) is a convex combination as well. The proof is obvious. Since x 0 < x 2 , all functions in front of y 0 , C y and y 2 are non-negative and have a sum of 1. The control points (x 0 , y 0 ), (C x , C y ) and (x 2 , y 2 ) form a convex hull, highlighted in grey in Figure 8, that encloses all realizations of x in (46), which means that all points of the parabola lie within the convex hull. A convex hull is the smallest convex set that contains a given set, with a convex set being a set wherein the straight line segment, connecting any two points of the set is completely contained within the set; see e.g., the textbook by Farin [7] (p. 439). The advantage of the general representation of the Bernstein polynomial is that it can be applied directly without scaling the data set to the interval [x 0 , If we now consider the example from Table 1 once again, we get, with y 0 = 3.5 and y 2 = 7.0, from (44) the result C y = 5.25 and therefore which can be "simplified" to the monomial basis form This is of course the same as (10) and (36). For the case of the interpolation by a parabola, it can be concluded that both the Bernstein and the Lagrangian form yield the same result. However, the Bernstein form has the advantage that all points of the parabola lie in the convex hull defined by its parameters.

Transition from Bernstein Form to B-spline
In the previous section, the Bernstein form was derived with a more general approach, instead of the standard representation within the interval [0, 1]. This form can be easily used for derivation of the B-spline.
In this section, Equation (46) is used for further derivations. It is assumed that the line segments that define the quadratic polynomial intersect at the extremum of the parabola and have an identical coordinate C y i . Since the notion of spline is to be explained, only the endpoints of two smoothly connected parabolas are considered, which are better known as knots (x 0 , y 0 ), (x 1 , y 1 ) and (x 2 , y 2 ). The problem of the constraints and how they are imposed on the spline, for spline interpolation, is not taken into consideration. It is presumed that the first parabola of the spline is already defined, e.g., by imposing a constraint for the direction of the tangent in point x 0 as boundary condition, and in this case, is identical with the one from the previous section, see Figure 8.
The initial situation is depicted in Figure 9, where: -One parabola already exists between the knots (x 0 , y 0 ) and (x 1 , y 1 ); and -One knot (x 2 , y 2 ) is given that is to be interpolated. The task is now to construct a second parabola between the knots 1 1 ( , ) x y and 2 2 ( , ) x y which is tangential to the first one.
The second parabola is constructed as follows: − The line segment The task is now to construct a second parabola between the knots (x 1 , y 1 ) and (x 2 , y 2 ) which is tangential to the first one. The second parabola is constructed as follows: - The line segment (x 0 , C y 1 )(x 1 , y 1 ) is extended from the end point (x 1 , y 1 ) of the parabola to the end of the interval [x 1 , x 2 ], yielding the additional coordinates (x 2 , C y 2 ) and therefore the new line segment (x 1 , y 1 )(x 2 , C y 2 ), see Figure 12; -Between points (x 1 , C y 2 ) and (x 2 , y 2 ), the line segment (x 1 , C y 2 )(x 2 , y 2 ) can be constructed, see Figure 10; -By taking combinations of these two line segments, within the interval [x 1 , x 2 ], another parabola will be constructed, smoothly connected to the previous one. The whole function constructed from these two smoothly connected parabolas represents a spline. Both parabolas end tangentially at the point (x 1 , y 1 ) w.r.t the line segment, defined as (C x 1 , C y 1 )(C x 2 , C y 2 ), see Figure 11.  [ , ] x x , smoothly connected to the previous one.
The extrema of the parabolas are marked with a red dot, the control points with a blue dot, and the green dots indicate points to be interpolated.
Finally, this is a type of spline that is constructed by using convex combinations and it is exactly the same as the B-spline. This is shown by the equations derived further in the text. Algebraically, this derivation is similar to the one in the previous section. The first pair of line equations is written as By taking a combination of these two lines within the interval 0 1 [ , ] x x , it follows For the second parabola the same approach is applied. The second pair of line equations is defined in the same way as the previous one, yielding  [ , ] x x , smoothly connected to the previous one.
The extrema of the parabolas are marked with a red dot, the control points with a blue dot, and the green dots indicate points to be interpolated.
Finally, this is a type of spline that is constructed by using convex combinations and it is exactly the same as the B-spline. This is shown by the equations derived further in the text. Algebraically, this derivation is similar to the one in the previous section. The first pair of line equations is written as By taking a combination of these two lines within the interval 0 1 [ , ] x x , it follows 1 2 For the second parabola the same approach is applied. The second pair of line equations is defined in the same way as the previous one, yielding Figure 11. Second parabola within the interval [x 1 , x 2 ], smoothly connected to the previous one. The extrema of the parabolas are marked with a red dot, the control points with a blue dot, and the green dots indicate points to be interpolated.
− By taking combinations of these two line segments, within the interval 1 2 [ , ] x x , another parabola will be constructed, smoothly connected to the previous one. The whole function constructed from these two smoothly connected parabolas represents a spline. Both parabolas end tangentially at the point 1 1 ( , ) x y w.r.t the line segment, defined as ( , )( , ) x y x y C C C C , see Error! Reference source not found.. x y x C . The extremum of the first parabola is marked with a red dot, the control point with a blue dot, and the green dots indicate points to be interpolated. Finally, this is a type of spline that is constructed by using convex combinations and it is exactly the same as the B-spline. This is shown by the equations derived further in the text. Algebraically, this derivation is similar to the one in the previous section. The first pair of line equations is written as By taking a combination of these two lines within the interval [x 0 , For the second parabola the same approach is applied. The second pair of line equations is defined in the same way as the previous one, yielding Analogously to (50), but in this case within the interval [x 1 , x 2 ], we obtain If the above derived expressions are compared to an explicit expression of a quadratic B-spline it becomes clear that they are the same. However, for the sake of clarity we have used a particular notation that does not correspond to the usual B-spline notation.
Looking at (50) and (52), the only unknowns are C y 1 and C y 2 , all other values are given. It appears that a quadratic spline, constructed from two linear polynomials, can be estimated with only two equations. This type of reasoning would lead to a paradox. Therefore, for this type of a spline, since all constraints (smoothness and continuity) are hidden inside the equations, at least four observations are necessary for a solution. For comparison, a spline constructed from two polynomials of the form a 0 + a 1 x + a 2 x 2 requires at least three observations and three constraints.
The idea presented in this section can be used for constructing splines of higher degrees as well. In those cases, depending on the degree of the spline, one would need more knots and more control points for creating a spline segment, while the number of combinations increases. Moreover, the expressions became more complicated, which makes this approach unsuitable for practical applications. However, de Boor [3] developed an algorithm for computing the basis functions that, in our case, constructs the expressions in front of the coordinates of the control points in (50) and (52).
De Boor's algorithm is based on the Cox-de Boor recursion formula that we will use for explaining the computation of a B-Spline. Recursion itself is non-intuitive, which makes this formula difficult to be explained analytically. The higher the degree of the spline, the more cumbersome and harder to follow the equations become. Therefore, in order to make a comparison with (50) and (52), we will use an explicit expression of a quadratic B-spline that is derived from the Cox-de Boor recursion formula. For this purpose, we use the knots and the naming convention from the previous example (49)-(52). For this formula to work, we have to rename the parameters (control points) y 0 to C y 0 and y 2 to C y 3 . The Cox-de Boor recursion requires additional knots at the beginning and at the end of the spline based on the degree of the spline. If the spline degree is d = 2, two additional knots are presumed to be both at the beginning and at the end of the spline. In our case, they are which is called knot multiplication. Additional knots with e.g., are especially used to gain endpoint interpolation. The additional knots do not necessarily have to be equal to x 0 and x 2 , they are just required to follow the order x −2 ≤ x −1 ≤ x 0 and x 2 ≤ x 3 ≤ x 4 . However, in our example we follow the convention of the knot multiplication, and because of the knot multiplication, formula (56) has zero divisions, and the solution for this problem is to declare that "anything divided by zero is equal to zero". The base case of the recursion formula is with a recursive step to compute the basis functions which is summed over the k-th interval of the spline f k,d (x), such as At the beginning of the recursion (56), index d denotes the degree of the spline and i is a computed index i = k − d, . . . , k. In each recursive step, d is lowered down for 1 until it reaches 0. At d = 0, the recursion formula terminates and the base case (55) is evaluated.
The index k is the index of the interval [x k , x k+1 ] with k = 0, 1, . . . , l − 1, where l is the number of knots. Therefore, for k = 0, the recursion formula (56) is computed for every i of the sum (57) and evaluated in (55). Afterwards, for k = 1, the same computation is performed, and so on, until the last interval. The parameters C y i in (57) correspond to each basis function terminated by B i,0 (x) in (55) at the end of the recursion. An example of quadratic basis functions using the knots x k 0 = 0, x k 1 = 1 and x k 2 = 3, resulting in two intervals k = 0 and k = 1, and the knot sequence [0, 0, 0, 1, 3, 3, 3], is shown in Figure 13.
. Afterwards, for 1 k = , the same computation is performed, and so on, until the last interval. The parameters   To compare this approach with Error! Reference source not found. and Error! Reference source not found., we use an explicit expression of a generic quadratic B-spline formula. For deriving it, we consider  To compare this approach with (50) and (52), we use an explicit expression of a generic quadratic B-spline formula. For deriving it, we consider d = 2 in (56), then we follow all recursive steps until d = 0 for every B i,d (x), which yields the explicit expression of a generic quadratic B-spline formula For clearer description, we treat each part of (57) separately and, for the computation, we use the knot sequence (53). At the interval k = 0, with index i = −2, the first element of the sum is With i = −1, the second element of the sum is and, with i = 0, the last element of the sum is Taking the knot multiplications (53) and (54) into account for all elements of the sums (59)-(61) yields The base case (55) (55), we obtain If Equation (65) is compared to (50), it can be concluded that they are identical. It must just be taken into consideration, as stated previously, that y 0 is renamed to C y 0 . If the same procedure is applied to k = 1, the result is identical to (52), considering that y 2 is renamed to C y 3 .
Since spline functions y = f (x) are considered, only the component C y i of a control point (C x i , C y i ) appears as a factor in front of the basis functions B i,d (x). Using a parametric spline curve, both components of a control point appear in front of the basis functions as e.g., applied by Bureick et al. [34].
The Cox-de Boor recursion Formulas (55) and (56) are the basis of the de Boor's algorithm. The presented functions in front of the parameters are called basis functions and the curve is called B-spline curve.

B-spline Approximation
As pointed out by Ezhov et al. [1], in engineering geodesy it is not appropriate to apply a spline interpolation due to the high point density and the fact that measurements are affected by random measurement errors. Observation errors and other abrupt changes in the data points would be modeled, resulting in a strongly oscillating spline. The solution to this is to divide the sequence of points into not so many intervals, determined by predefined knots. The distribution of the knots is a fundamental problem in spline approximation and different knot placement strategies can be applied to solve it, as explained by Ezhov et al. [1]. Within the resulting intervals, we consider an overdetermined configuration and, hence, the parameters of a B-spline can be computed by least squares adjustment.
Since the Cox-de Boor Formulas (55)-(57) manage the whole process of computing the coefficients of the unknowns, the approximation can be performed with a simple linear functional model. Here, we consider quadratic splines that consist of piecewise parabolic segments. A generalization to splines of higher degree is straightforward.

Definition of the Problem
Let: y 0 , y 1 , . . . , y j , . . . , y m be a sequence of observed values; -Σ LL = σ 2 0 Q LL be the variance-covariance matrix of the observations y j ; - x 0 < x 1 < . . . < x j < . . . < x m be a sequence of non-decreasing error-free values referring to the observations; - x k 0 < . . . < x k n−1 be a non-decreasing sequence of user-defined values for the knots on the x-axis which are regarded as error-free. To determine: -C y 0 , C y 1 , . . . , C y n , which represent the y-component of the control "points" of the quadratic B-spline.
After defining these quantities, we can create a spline approximation from an overdetermined configuration using the Gauss-Markov model.

Formulation of the Adjustment Problem
The observation equations are of the same form for each quadratic polynomial within a spline. This is a consequence of the Cox-de Boor recursion formula. This approach (55), (56) and (57) determines to which interval a value x j belongs. One has to consider that the knot sequence must be given in advance.
Considering x j as fixed values, yields also B i,d (x), computed from (55), (56), and (57), as fixed values. Here, we assume that B i,2 (x) corresponds to the given parameter C y i . Explaining the observation equations explicitly by means of the Cox-de Boor's formula is too cumbersome. Hence, with y j as observations, the observation equations can be set up. The observation vector L = y 0 y 1 y 2 · · · y m T (67) contains all observations, while the vector of unknown parameters can be written aŝ From the stochastic model with σ 2 0 as theoretical variance factor, the corresponding weight matrix is obtained from supposing the cofactor matrix to be non-singular. By looking at the observation Equation (66), it is obvious that this spline approximation problem is linear and, hence, can easily be written in matrix notation where A is the design matrix that contains the coefficients of the unknowns. Equation (71), together with (69), is denoted as Gauss-Markov model; see e.g., the textbook by Niemeier [35] (p. 137). Considering the sequence of the unknowns in (68), the design matrix reads Even though the construction of this matrix is quite trivial, the presence of the basis functions makes it rather unusual. However, with the derivations in Section 4, it is at least clear how these basis functions can be obtained. Looking at the notation in (72), one might be misled that A is a full matrix, which is not the case. To put it simply, the Cox-de Boor formula evaluates only those x j that belong to a particular interval and everything else is equal to zero.

Least Squares Adjustment
One feature of the B-splines is that constraints for continuity, smoothness and continuity of curvature are accounted for implicitly in the functional model. Therefore, no additional constraints have to be introduced and the least squares solution for the unknowns can be obtained fromX With the residuals v = AX − L (74) from (71), the a posteriori estimate of the reference standard deviation can be computed, where r is the redundancy of the adjustment problem. Further details on the a posteriori analysis of adjustment results and on knot placement strategies can be found in the article by Ezhov et al. [1]. The fact that the solution for the unknowns (68) can be used for a transition "backwards" from B-spline to ordinary polynomial is shown in Appendix A.

Conclusions and Outlook
In engineering geodesy, point clouds derived from areal measurement methods, such as terrestrial laser scanning or photogrammetry, are often approximated by a continuous mathematical function for further analysis, such as deformation monitoring. In many cases, the formulas for B-spline curves and B-spline surfaces, given in the textbook by Piegl and Tiller [2] (pp. 81 and 100), are applied, where the functional values of the B-spline basis functions are recursively computed according to the formulas by de Boor [3] and Cox [4]. This approach is very easy to handle and results in a numerically stable solution for the unknowns to be determined. As these formulas have a very complex mathematical derivation, but are still very easy to use, they are mostly used like a given recipe without a deeper understanding of their derivation.
In part 1 of a series of three articles, Ezhov et al. [1] explained the basic methodology of spline approximation using splines constructed from ordinary polynomials. In this paper (part 2) the goal was to develop an alternative derivation of the B-spline. To avoid excessive formula derivations and to illustrate the geometric relationships, quadratic splines were considered that consist of piecewise parabolic segments. Starting with the representation of a parabola in the monomial basis (ordinary polynomial) a transition to the Langrangian form was performed and, from there, to the Bernstein form, which finally resulted in the B-spline representation. In all investigations univariate splines were used, in the form of a spline function y = f (x).
For the derivation of the formulas, we first considered the case of interpolation. The developed formulas were then used for spline approximation by means of least squares adjustment. With the values y i as observations and x i as error-free values, the resulting linear adjustment problem could be solved within the Gauss-Markov model. Finally, in Appendix A it was shown that the determined spline parameters can be used for a transition "backwards" from B-spline to ordinary polynomial. The transition "forwards" from ordinary polynomial to B-spline was shown in Appendix B.
The more general case, where both values y i and x i are introduced as observations into an approximation with a spline function, was already elaborated by Neitzel et al. [36]. Investigations of the numerical stability of the spline approximation approaches, based on ordinary polynomials and on truncated polynomials, discussed in part 1 and the B-splines explained in this article, as well as the potential of splines in deformation detection, will be presented in a forthcoming part 3 of this series of three articles.
Author Contributions: Conceptualization, F.N. and S.P.; methodology, N.E.; writing-original draft preparation, F.N. and N.E.; writing-review and editing, F.N. and S.P. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.
From Equations (A10)-(A15), two separate equation systems can be formed to solve for the B-spline parameters. Since the equations are linear, they can be solved by applying some of the methods from linear algebra, such as Gaussian elimination.
Once the systems are solved, the results for the parameters of the system (A10)-(A12) are and for the system (A13)-(A15), the results are In the terms of B-spline, y 0 and y 2 are parameters and, in this derivation, they are treated as such. If y 0 and y 2 are considered as known values, we can also solve only for C 1 and C 2 as a system of two equations. However, in such case, the solutions for C 1 and C 2 are not as elegant as those from the equations above.