mathematics
Article
Spline Approximation, Part 2: From Polynomials
in the Monomial Basis to B-splines—A Derivation
Nikolaj Ezhov 1,*, Frank Neitzel 1and Svetozar Petrovic 1,2
Citation: Ezhov, N.; Neitzel, F.;
Petrovic, S. Spline Approximation,
Part 2: From Polynomials in the
Monomial Basis to B-splines—A
Derivation. Mathematics 2021,9, 2198.
https://doi.org/10.3390/math
9182198
Academic Editor: Clemente Cesarano
Received: 13 July 2021
Accepted: 25 August 2021
Published: 8 September 2021
Publisher’s Note: MDPI stays neutral
with regard to jurisdictional claims in
published maps and institutional affil-
iations.
Copyright: © 2021 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).
1Institute of Geodesy and Geoinformation Science, Technische Universität Berlin, 10623 Berlin, Germany;
2Section 1.2: Global Geomonitoring and Gravity Field, GFZ German Research Centre for Geosciences,
14473 Potsdam, Germany
*Correspondence: [email protected]
Abstract:
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 finally 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.
Keywords:
spline; B-spline; polynomial; monomial; basis change; Lagrange; Bernstein; interpolation;
approximation; least squares adjustment
1. 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.
Mathematics 2021,9, 2198. https://doi.org/10.3390/math9182198 https://www.mdpi.com/journal/mathematics
Mathematics 2021,9, 2198 2 of 24
A brief overview of the historical development of the B-spline can be found, for
example, in the textbook by Farin [
7
] (p. 119). He explains that the forerunner of the
B-spline first appeared in a publication on approximation theory by the mathematician
Lobatschewsky [
8
]. It was constructed as a convolution of certain probability distributions.
However, at that time the term “spline” has not been introduced yet. Subsequently, in
1946, Schoenberg used B-splines in his seminal publications [
9
,
10
] for data approximation.
It is commonly accepted that, with these articles, the modern representation of spline
approximation began and the term “spline” was used for the first time, see the historical
notes in the textbook by Schumaker [
11
] (p. 10). Further on, many authors dealt with the
utilization and the development of the B-splines. Farin [
7
] (p. 119) points out that one of
the most important developments in this theory was the introduction of the recurrence
relations, discovered by de Boor [
3
], Cox [
4
], and Mansfield. The contribution of Mansfield
to this development is mentioned by de Boor in [
3
]. However, the recurrence relation had
already appeared in contributions by Popoviciu and Chakalov in the 1930s, see the article
by de Boor and Pinkus [12] and the literature cited therein.
B-splines were later fully exploited in computer graphics. Consequently, most papers
about B-splines were published in this field, see for example, the contributions by the
authors Piegl and Tiller [
2
], Szilvasi-Nagy [
13
], Farin [
7
], and Zheng et al. [
14
]. However,
there are several publications where B-splines were used for statistical data analysis, e.g.,
the contributions by Wold [15] or Anderson and Turner [16].
Applications of (B)-splines in geodesy can be traced back to 1975, when Sünkel began
using them for various problems. In [
17
], bicubic spline functions were used for the
reconstruction of functions from discrete data. In [
18
], the representation of geodetic
integral formulas by bicubic spline functions was shown using splines constructed from
ordinary cubic polynomials, see [
18
] (p. 10 ff.), similar to the representation used later
by Ezhov et al. [
1
]. In [
19
], the local interpolation of gravity data by spline functions was
derived, and in [
20
] (p. 18 ff.) B-splines were used for smooth surface representation. As
references for the B-spline functions, the publications of Schoenberg [
21
], Schumaker [
22
],
and Späth [23] are listed by Sünkel [20].
On the basis of B-splines on the line, Jekeli [
24
] (p. 7 ff.) showed the use of spherical
B-splines for representations of functions on a sphere for geopotential modeling. Mautz
et al. [
25
] used B-splines to develop a representation of global ionosphere maps based on
B-spline wavelets. The local improvement of the International Reference Ionosphere (IRI)
by means of an N-dimensional B-spline surface was developed by Koch and Schmidt [
26
].
The application of (B)-splines in engineering geodesy can be traced back to 1985, when
Dzapo et al. [
27
] used cubic splines for the determination of the lengths of railroad tracks. A
recent application of B-spline curves and B-spline surfaces is the approximation of 3D point
cloud data, as shown e.g., by Bureick et al. [
28
]. Kermarrec et al. [
29
] applied hierarchical
B-splines for the approximation of 3D point cloud data obtained from measurements
with a terrestrial laser scanner. Furthermore, B-spline models play an important role in
spatio-temporal deformation analysis as shown by Harmening [30].
The main goal of this paper is to present an alternative derivation of the B-spline
function to the purely mathematical derivation presented by Schoenberg [
9
] and the
recursion formulas for the B-spline function developed by de Boor [
3
]. This derivation
is based on the transition from the representation of a polynomial in the monomial basis
(ordinary polynomial) to the Lagrangian form and, from there, to the Bernstein form,
which finally results in the B-spline representation. In all investigations, univariate splines
are used in form of a spline function
y=f(x)
. For the derivation of the formulas, we
first consider the case of interpolation. The developed formulas are then used for spline
approximation using least squares adjustment.
The general interpolation problem is, to find for
n+
1 given data points
(x0
,
y0)
,
(x1
,
y1)
,
. . .
,
(xn
,
yn)
, with
x0<x1<. . . <xn
a function
f:R→R
in such a way that the
interpolation condition
f(xi) = yi,i=0, 1, . . . , n(1)
Mathematics 2021,9, 2198 3 of 24
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 computa-
tional 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:
- Polynomials;
- Trigonometric functions;
- Exponential functions;
- Logarithmic functions;
- Rational functions.
Ezhov et al. [
1
] elaborated a spline approximation using piecewise ordinary cubic poly-
nomials. 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 coeffi-
cients
ai
and so-called basis functions
φ0(x)
,
. . .
,
φn(x)
and the interpolating function is
defined as a linear combination
f(x) =
n
∑
i=0
aiφi(x)(2)
of these basis functions. Considering the interpolation condition (1), we obtain
f(xi) =
n
∑
i=0
aiφi(xi) = yi,i=0, 1, . . . , n. (3)
The simplest and most common type of interpolation uses polynomials. These poly-
nomials 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 well-
known monomial basis applied by Ezhov et al. [
1
] often lies in the fact that for the com-
putation 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.
2. Transition from Ordinary Polynomial to Lagrangian Form
With the monomial basis functions
φi(x) = xi,i=0, 1, . . . , n, (4)
Mathematics 2021,9, 2198 4 of 24
we obtain for the parabola with n=2 from (3)
P2(x) = a0+a1x+a2x2. (5)
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,x2]T
and the coefficients in
a vector a= [a0,a1,a2]T, we obtain
P2(x) = aTm(x). (6)
A visualization of the monomial functions for the parabola is shown in Figure 1.
Mathematics 2021, 9, x FOR PEER REVIEW 4 of 27
dinary polynomial is shown in Appendix A. The transition “forwards” from ordinary
polynomial to B-spline is shown in Appendix B.
2. Transition from Ordinary Polynomial to Lagrangian Form
With the monomial basis functions
() , 0,1, ,
i
i
x
xi n
φ
==, (4)
We obtain for the parabola with 2n= from Error! Reference source not found.
2
2012
()
P
xaaxax=+ + . (5)
A polynomial in the monomial basis is often referred to as standard form or ordi-
nary polynomial. Arranging the monomials in a vector 2T
() [1,, ]
x
xx=m and the coeffi-
cients in a vector T
012
[,,]aaa=a, we obtain
T
2() ()
P
xx=am . (6)
A visualization of the monomial functions for the parabola is shown in Error! Ref-
erence source not found..
Figure 1. Monomial basis functions for a parabola in the interval [0, 1].
Considering the interpolation condition Error! Reference source not found., the
coefficients can be determined from the solution of the linear equation system
=
y
Va , (7)
with T
012
[,,]yyy=y and the Vandermonde matrix
2
00
2
11
2
22
1
1
1
x
x
x
x
x
x
=
V, (8)
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 Error! Ref-
erence source not found., the determination of the coefficients will be demonstrated.
Table 1. Numerical example for interpolating a parabola.
i xi y
i
0 0 1
1 1 3
2 3 2
Figure 1. Monomial basis functions for a parabola in the interval [0, 1].
Considering the interpolation condition (1), the coefficients can be determined from
the solution of the linear equation system
y=Va, (7)
with y= [y0,y1,y2]Tand the Vandermonde matrix
V=
1x0x2
0
1x1x2
1
1x2x2
2
, (8)
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.
Table 1. Numerical example for interpolating a parabola.
i xiyi
0 0 1
1 1 3
2 3 2
From Equation (7), we obtain for this numerical example
1
3
2
=
100
111
139
a0
a1
a2
. (9)
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
a0=
1,
a1=17
6
,
a2=−5
6
. Thus, the equation of the
parabola reads
P2(x) = 1+17
6x−5
6x2. (10)
It can be stated that, when using the monomial basis, the numerical effort to determine
the coefficients requires roughly
n3
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.
2.1. Transition by Means of Basis Transformation
We want to solve the problem of basis transformation between two different sets of
coefficients aiand bifulfilling
P2(x) =
2
∑
i=0
aiφi(x) =
2
∑
i=0
biψi(x), (11)
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
P2(x) =
2
∑
i=0
aiφi(x) = a0+a1x+a2x2=a0a1a2
1
x
x2
, (12)
as already shown in (5) and (6). Multiplying out the second summation formula in
(11) yields
P2(x) =
2
∑
i=0
biψi(x) = b0ψ0(x) + b1ψ1(x) + b2ψ2(x) = b0b1b2
ψ0(x)
ψ1(x)
ψ2(x)
(13)
and, according to (11), we can write
a0a1a2
1
x
x2
=b0b1b2
ψ0(x)
ψ1(x)
ψ2(x)
. (14)
In the case of monomial basis, we selected basis functions and obtained the coefficients,
see (7), while here we do the opposite. We select
b0=y1
,
b1=y2
,
b2=y3
for the coefficients
to obtain corresponding basis functions for which we introduce the new notation
li(x)
with
i=
0, 1, 2. Furthermore, for
a0
,
a1
,
a2
we introduce the solution according to (7) and
we obtain
y0y1y2
111
x0x1x2
x2
0x2
1x2
2
−1
1
x
x2
=y0y1y2
l0(x)
l1(x)
l2(x)
. (15)
Mathematics 2021,9, 2198 6 of 24
By comparing the coefficients of the vectors
[y0y1y2]
, we see that the basis
functions li(x)can be computed from
l0(x)
l1(x)
l2(x)
=
111
x0x1x2
x2
0x2
1x2
2
−1
1
x
x2
. (16)
Introducing the vector l(x) = l0(x)l1(x)l2(x)Tyields
l(x) = (VT)−1m(x). (17)
The same formula can also be obtained taking the basis transformation from La-
grange to monomials
VTl(x) = m(x)
, as presented by Gander [
31
], and solving it for
l(x)
.
From (17), we obtain the result
l(x) = h(x−x1)(x−x2)
(x0−x1)(x0−x2)
(x−x0)(x−x2)
(x1−x0)(x1−x2)
(x−x0)(x−x1)
(x2−x0)(x2−x1)iT, (18)
which is called the Lagrange basis. A visualization of the Lagrange basis functions for the
parabola is shown in Figure 2.
Mathematics 2021, 9, x FOR PEER REVIEW 6 of 27
In the case of monomial basis, we selected basis functions and obtained the coeffi-
cients, see (7), while here we do the opposite. We select 01
by=, 12
by=, 23
by= for the
coefficients to obtain corresponding basis functions for which we introduce the new no-
tation ()
i
lx with 0,1, 2i=. Furthermore, for 012
,,aaa
we introduce the solution ac-
cording to Error! Reference source not found. and we obtain
[] []
1
0
012 01 2 0121
222 2
01 2 2
111 1 ()
()
()
lx
yyyxxx x yyylx
x
xx x lx
−
=
. (15)
By comparing the coefficients of the vectors 012
[]yyy
, we see that the basis
functions ()
i
lx
can be computed from
1
0
1012
222 2
2012
() 1 1 1 1
()
()
lx
lx x x x x
lx x x x x
−
=
. (16)
Introducing the vector
[]
T
012
() () () ()
x
lx lx lx=l yields
T1
() ( ) ()
x
x
−
=lVm
. (17)
The same formula can also be obtained taking the basis transformation from La-
grange to monomials T() ()
x
x=Vl m , as presented by Gander [31. ], and solving it for
()
x
l. From Error! Reference source not found., we obtain the result
T
02 01
12
0102 1012 2021
()()()()
()()
() ()()()()()()
xx xx xx xx
xxxx
xxxxx xxxx xxxx
−− −−
−−
=
−− −− −−
l, (18)
which is called the Lagrange basis. A visualization of the Lagrange basis functions for the
parabola is shown in Error! Reference source not found..
Figure 2. Lagrange basis functions for a parabola in the interval [0, 3] for the example 00x=,
11x=and 23x=.
Finally, the resulting polynomial reads
02 01
12
2012
0102 1012 2021
()() ()()
()()
() ()()()()()()
xx xx xx xx
xxxx
Px y y y
xxxx xxxx xxxx
−− −−
−−
=++
−− −− −−. (19)
An alternative derivation of the Lagrange representation of a polynomial is pre-
sented in the following subsection.
Figure 2.
Lagrange basis functions for a parabola in the interval [0, 3] for the example
x0=
0,
x1=
1
and x2=3.
Finally, the resulting polynomial reads
P2(x) = (x−x1)(x−x2)
(x0−x1)(x0−x2)y0+(x−x0)(x−x2)
(x1−x0)(x1−x2)y1+(x−x0)(x−x1)
(x2−x0)(x2−x1)y2. (19)
An alternative derivation of the Lagrange representation of a polynomial is presented
in the following subsection.
2.2. 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) = a0+ (a1+a2x)x
, which represents a combination (product) of two
lines plus some constant value
a0
. 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 yicoordinates.
We consider the points (x0,y0),(x1,y1),(x2,y2), where x0<x1<x2. Now two lines
f0,1(x) = a0+a1x(20)
Mathematics 2021,9, 2198 7 of 24
within the interval [x0,x1)and
f1,1(x) = b0+b1x(21)
within the interval
[x1
,
x2]
are defined, see Figure 3. The subscripts
i
and
d
in
fi,d
represent
the number of the polynomial ifor a given degree d.
Mathematics 2021, 9, x FOR PEER REVIEW 7 of 27
2.2. Transition by Linear Combinations
In this subsection, we present an alternative derivation of the polynomial in the La-
grange representation. We use the fact that each parabola can be represented as a com-
bination of two lines. In fact, if we take a look at Error! Reference source not found., we
find that it can be rearranged into the form
012
() ( )
f
xa aaxx=+ +
, which represents a
combination (product) of two lines plus some constant value
0
a
. By using this idea, the
Lagrange representation of a quadratic polynomial can be derived, which results in a
linear combination of three quadratic functions
()px
and the respective i
y
coordi-
nates.
We consider the points
00 11 22
( , ),( , ),( , )
x
yxyxy
, where
012
x
xx
<<
. Now two lines
0,1 0 1
()fxaax=+ (20)
within the interval
01
[,)
x
x
and
1,1 0 1
()fx bbx=+ (21)
within the interval
12
[, ]
x
x
are defined, see Error! Reference source not found.. The
subscripts
i
and
d
in
,id
f
represent the number of the polynomial
i
for a given de-
gree
d
.
Figure 3. Line interpolation.
From two given points, the slope of the line Error! Reference source not found. can
be written as
10
1
10
yy
axx
−
=−
. (22)
The y-intercept can be computed from
10
000
10
yy
ayx
xx
−
=− −
or 10
011
10
yy
ayx
xx
−
=− −
. (23)
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
0
1
0,1 0 1
10 10
() xx
xx
fx y y
xx xx
−
−
=+
−−
. (24)
Analogously, the equation for the second line can be written as
Figure 3. Line interpolation.
From two given points, the slope of the line (20) can be written as
a1=y1−y0
x1−x0
. (22)
The y-intercept can be computed from
a0=y0−x0
y1−y0
x1−x0
or a0=y1−x1
y1−y0
x1−x0
. (23)
Inserting (22) and (23) into (20) yields, after some rearrangement, for the first line
f0,1(x) = x1−x
x1−x0
y0+x−x0
x1−x0
y1. (24)
Analogously, the equation for the second line can be written as
f1,1(x) = x2−x
x2−x1
y1+x−x1
x2−x1
y2. (25)
By taking a combination of these two line equations, where
x
varies within the interval
[x0,x2], the expression
f0,2(x) = x2−x
x2−x0
f0,1(x) + x−x0
x2−x0
f1,1(x)(26)
for the parabola is derived, which represents a linear combination of three quadratic
polynomial functions and the coordinates y0,y1and y2.
To explain how the factors in front of the terms
f0,1(x)
and
f1,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
(x0
,
y0)
,
(x1
,
y1)
,
. . .
,
(xn
,
yn)
, 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
xi
,
the following polynomials
px0(x) = (x−x1) (x−x2)
(x0−x1) (x0−x2), (27)
Mathematics 2021,9, 2198 8 of 24
px1(x) = (x−x0) (x−x2)
(x1−x0) (x1−x2), (28)
px2(x) = (x−x0) (x−x1)
(x2−x0) (x2−x1)(29)
can be created.
As previously mentioned, these functions
pxi(x)
are equal to 1 when the variable
x
is
equal to the corresponding coordinate
xi
of the given point and 0 at all other given points.
Therefore, by multiplying each of these functions
pxi(x)
with their corresponding values
yi
ensures that, at the particular given point
xi
, the result will be
yi
and 0 at all other given
points. The summation of all these functions
pxi(x)
, multiplied by their corresponding
values yi, results in a Lagrangian interpolation polynomial
f(x) = (x−x1)(x−x2)
(x0−x1)(x0−x2)y0+(x−x0)(x−x2)
(x1−x0)(x1−x2)y1+(x−x0)(x−x1)
(x2−x0)(x2−x1)y2. (30)
The quadratic polynomial (30) represents a linear combination of the Lagrangian basis
polynomials with the constants
y0
,
y1
and
y2
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
f0,1(x) = x−x1
x0−x1
y0+x−x0
x1−x0
y1, (31)
f1,1(x) = x−x2
x1−x2
y1+x−x1
x2−x1
y2(32)
are derived. If we follow the already described Lagrange’s idea to obtain a quadratic
polynomial, the following multiplications
f0,2(x) = x−x2
x0−x2
f0,1(x) + x−x0
x2−x0
f1,1(x)(33)
must be performed.
Consequently, Equation (33) results in a Lagrangian polynomial of second order.
Looking at (33), the factors
(x−x2)/(x0−x2)
and
(x−x0)/(x2−x0)
are those that appear
in front of the terms
f0,1(x)
and
f1,1(x)
in (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.
(x−x2)/(x0−x2)
. In (24), (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
f0,1(x)
and
f1,1(x)
in (26) is explained, inserting
(24) and (25) into (26) finally yields
f0,2(x) = (x1−x)(x2−x)
(x1−x0)(x2−x0)y0+(x−x0)(x2−x)
(x1−x0)(x2−x1)y1+(x−x0)(x−x1)
(x2−x0)(x2−x1)y2. (34)
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.
Mathematics 2021,9, 2198 9 of 24
Mathematics 2021, 9, x FOR PEER REVIEW 9 of 27
21
1,1 1 2
12 21
() xx xx
fx y y
xx xx
−−
=+
−−
(32)
are derived. If we follow the already described Lagrange’s idea to obtain a quadratic
polynomial, the following multiplications
0
2
0,2 0,1 1,1
02 20
() () ()
xx
xx
fx fx fx
xx xx
−
−
=+
−−
(33)
must be performed.
Consequently, Equation Error! Reference source not found. results in a Lagrangian
polynomial of second order. Looking at Error! Reference source not found., the factors
202
()/( )
x
xxx
−−
and
020
()/( )
x
xxx
−−
are those that appear in front of the terms
0,1
()
f
x and
1,1
()
f
x in Error! Reference source not found.. 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.
202
()/( )
x
xxx
−−
. In (24), (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
0,1
()
f
x and
1,1
()
f
x in
Error! Reference source not found. is explained, inserting
Error! Reference source not found. and Error! Reference source not found. into
Error! Reference source not found. finally yields
02 0 1
12
0,2 0 1 2
1020 1021 2021
()() ()()
()()
() ()()()()()()
xx x x xx xx
xxxx
fx y y y
xxxx xxxx xxxx
−− −−
−−
=++
−− −− −−
. (34)
This form coincides with Error! Reference source not found. and represents a
second degree Lagrange interpolating polynomial through three points, see the hand-
book by Bronshtein et al. [32. ] (p. 918). The plot of the parabola is depicted in Error!
Reference source not found..
Figure 4. Lagrange interpolation by lines (black) and by a parabola (red).
Inserting the values of the numerical example from Error! Reference source not
found. into Error! Reference source not found. yields the Lagrangian form
0,2
(1)(3) (0)(3) (0)(1)
() 1 3 2
(1 0)(3 0) (1 0)(3 1) (3 0)(3 1)
xx x x xx
fx −− −− −−
=⋅+⋅+⋅
−− −− −−
, (35)
which can be “simplified” to the monomial basis form
2
0,2
17 5
() 1 66
fx x x=+ −
. (36)
Figure 4. Lagrange interpolation by lines (black) and by a parabola (red).
Inserting the values of the numerical example from Table 1into (34) yields the La-
grangian form
f0,2(x) = (1−x) (3−x)
(1−0) (3−0)·1+(x−0) (3−x)
(1−0) (3−1)·3+(x−0)(x−1)
(3−0)(3−1)·2, (35)
which can be “simplified” to the monomial basis form
f0,2(x) = 1+17
6x−5
6x2. (36)
The result is the same as in (10), but the values of the parameters
a0
,
a1
,
a2
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.).
3. 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
[x0
,
x2]
in which the polynomial is defined. At the end of
the extended line segments, two additional y-coordinates for the values
x0
and
x2
are
computed, as seen in Figure 5. To distinguish these two y-coordinates from the given
values, a different notation, y0
0and y0
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
x1
is not necessarily in the middle of the interval
[x0
,
x2]
. Therefore, in general
x16= (x0+x2)/2.
Mathematics 2021,9, 2198 10 of 24
Mathematics 2021, 9, x FOR PEER REVIEW 10 of 27
The result is the same as in Error! Reference source not found., but the values of the
parameters
0
a
,
1
a
,
2
a
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 Error! Reference source not found., resp.
Error! Reference source not found., namely that each time a value x changes, the La-
grange basis polynomials must be recalculated. Further general limits of Lagrange in-
terpolation, e.g., that polynomial interpolants may oscillate, are explained by Farin [7. ]
(p. 101 ff.).
3. 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 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 Error! Reference source not found., the interpolating parts of the line
segments can be extended to both ends of the interval 02
[,]
x
x
in which the polynomial is
defined. At the end of the extended line segments, two additional y-coordinates for the
values 0
x
and 2
x
are computed, as seen in Error! Reference source not found.. To dis-
tinguish 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 02
[,]
x
x
. Therefore, in general 102
()/2
xxx
≠+
.
Figure 5. Extended line segments from Lagrangian form.
From Error! Reference source not found., as explained previously, the line equa-
tions
0
2
0,1 0 2
20 20
() xx
xx
fx y y
xx xx
−
−′
=+
−−
, 0
2
1,1 0 2
20 20
() xx
xx
fx y y
xx xx
−
−′
=+
−−
(37)
are derived. By taking combinations of these two lines, it follows
2 2
00 0
22 2
0,2 0 2 0 2
20 2020 2020 20
() xx xx xx
xx xx xx
fx y y y y
xx xxxx xxxx xx
−− −
−− −
′′
=+ + +
−−−−−−
, (38)
which can be brought into the form
[0,1]
Figure 5. Extended line segments from Lagrangian form.
From Figure 5, as explained previously, the line equations
f0,1(x) = x2−x
x2−x0
y0+x−x0
x2−x0
y0
2,f1,1(x) = x2−x
x2−x0
y0
0+x−x0
x2−x0
y2(37)
are derived. By taking combinations of these two lines, it follows
f0,2(x) = x2−x
x2−x02
y0+x2−x
x2−x0
x−x0
x2−x0
y0
2+x−x0
x2−x0
x2−x
x2−x0
y0
0+x−x0
x2−x02
y2, (38)
which can be brought into the form
f0,2(x) = x2−x
x2−x02
y0+2x2−x
x2−x0
x−x0
x2−x0y0
2+y0
0
2+x−x0
x2−x02
y2. (39)
A geometrical interpretation of the term (y0
2+y0
0)/2 is depicted in Figure 6.
Mathematics 2021, 9, x FOR PEER REVIEW 11 of 27
22
020 0
22
0,2 0 2
20 2020 20
() 2 2
xx y y xx
xx xx
fx y y
xx xxxx xx
′′
−+ −
−−
=+ +
−−− −
. (39)
A geometrical interpretation of the term 20
()/2
yy
′′
+
is depicted in Error! Reference
source not found..
Figure 6. Geometrical interpretation of the term
20
()/2yy
′′
+
.
Equation Error! Reference source not found. is a more general representation of the
Bernstein form, independent of the limits of the interval 02
[,]
x
x
, which means that x
can vary between any two real values. For proving the Weierstrass approximation theo-
rem, Bernstein used a special interval of length 1, i.e., 02
01
[,][,]
xx =
. Using this interval
simplifies the equations, which were later used by Bézier and de Casteljau. Therefore,
when using 02
01
[,][,]
xx =
, Equation Error! Reference source not found. results in
0,1 0 2
() (1 )fx xyxy
′
=− + ,
1,1 0 2
() (1 )fx xy xy
′
=− + (40)
and Error! Reference source not found. obtains the form
22
20
0,2 0 2
( ) (1 ) 2 (1 ) 2
yy
fx xy x x xy
′′
+
=− + − +
, (41)
in which
2
(1 )
x
−,
2(1 )xx−
and
2
x
are the Bernstein basis polynomials and 0
y
,
20
()/2yy
′′
+
and 2
y
are the Bernstein coefficients. In general, a polynomial in Bernstein
form of degree
n
can be written as
,
0
() ()
n
niin
i
Px B x
β
=
=
, (42)
where
i
β
are the Bernstein coefficients and
,
() (1 ) , 0,1, ,
ini
in
n
Bx x x i n
i
−
=− =
(43)
are the Bernstein basis polynomials; see e.g., the handbook by Bronshtein et al. [32. ] (p.
935). For a parabola, we obtain 2
0,2
(1 )Bx=−
,
1,2
2(1 )
B
xx=−
and
2
2,2
Bx=
, cf.
Error! Reference source not found.. These basis polynomials for a parabola are depicted
in Error! Reference source not found..
Figure 6. Geometrical interpretation of the term (y0
2+y0
0)/2.
Equation (39) is a more general representation of the Bernstein form, independent of
the limits of the interval
[x0
,
x2]
, 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.,
[x0
,
x2] = [
0, 1
]
. Using this interval simplifies the equations, which were later
used by Bézier and de Casteljau. Therefore, when using
[x0
,
x2]=[
0, 1
]
, Equation (37)
results in
f0,1(x) = (1−x)y0+xy0
2,f1,1(x) = (1−x)y0
0+xy2(40)
and (39) obtains the form
f0,2(x) = (1−x)2y0+2x(1−x)y0
2+y0
0
2+x2y2, (41)
Mathematics 2021,9, 2198 11 of 24
in which
(1−x)2
, 2
x(
1
−x)
and
x2
are the Bernstein basis polynomials and
y0
,
(y0
2+y0
0)/
2
and
y2
are the Bernstein coefficients. In general, a polynomial in Bernstein form of degree
ncan be written as
Pn(x) =
n
∑
i=0
βiBi,n(x), (42)
where βiare the Bernstein coefficients and
Bi,n(x) = n
ixi(1−x)n−i,i=0, 1, . . . , n(43)
are the Bernstein basis polynomials; see e.g., the handbook by Bronshtein et al. [
32
] (p. 935).
For a parabola, we obtain
B0,2 = (1−x)2
,
B1,2 =
2
x(
1
−x)
and
B2,2 =x2
, cf. (41). These
basis polynomials for a parabola are depicted in Figure 7.
Mathematics 2021, 9, x FOR PEER REVIEW 12 of 27
Figure 7. Bernstein basis polynomials for a parabola in the interval [0, 1].
In Error! Reference source not found., the term 20
()/2yy
′′
+, 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 in-
tersect 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 00
(, )
x
y and 22
(, )
x
y, replacing the point 11
(, )
x
y by
any other point on the same parabola, we will again obtain the equation of the form
Error! Reference source not found.. The equation of the parabola remains unchanged
although 0
y′ and 2
y′ have changed, since the only term depending on the new point is
02
()/2yy
′′
+. A parabola is described by a differentiable function. By computing its de-
rivative, equating it with zero and solving the resulting equation for
x
, we obtain the
coordinate of the extremum e
x
. Consequently, the resulting equation for the coordinate
of the extremum e
y can be solved. By substituting these two values for e
x
and e
y in a
line equation constructed for 0
y′ and 2
y′, where 0
x
x= and 2
x
x= respectively, it can
be shown that 02 y
yyC
′′
==. Therefore, 02
yy
′′
≠ implies 02
()/2
y
Cyy
′′
=+ . This proves
the above statement, which is illustrated in Error! Reference source not found..
M
and
N represent the fixed points at the beginning and at the end of the parabola, with
16
,,
mm and 16
,,
nn as end points at the opposite sides of the interval. Every pair of
line segments is defined as
()
0,
M
nNM,
()
11
,
M
nNm ,
()
22
,
M
nNm , …,
()
6
,
M
NNm
and depicted in light blue color. Error! Reference source not found. illustrates that all
these pairs of line segments intersect on the parabola.
Figure 7. Bernstein basis polynomials for a parabola in the interval [0, 1].
In (39), the term
(y0
2+y0
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
(x0
,
y0)
and
(x2
,
y2)
, replacing the point
(x1
,
y1)
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
y0
0
and
y0
2
have changed, since
the only term depending on the new point is
(y0
0+y0
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
xe
. Consequently, the
resulting equation for the coordinate of the extremum
ye
can be solved. By substituting
these two values for
xe
and
ye
in a line equation constructed for
y0
0
and
y0
2
, where
x=x0
and
x=x2
respectively, it can be shown that
y0
0=y0
2=Cy
. Therefore,
y0
06=y0
2
implies
Cy= (y0
0+y0
2)/
2. This proves the above statement, which is illustrated in Figure 8.
M
and
N
represent the fixed points at the beginning and at the end of the parabola, with
m1
,
. . .
,
m6
and
n1
,
. . .
,
n6
as end points at the opposite sides of the interval. Every pair of
line segments is defined as
Mn0,NM
,
Mn1,Nm1
,
Mn2,Nm2
,
. . .
,
MN,Nm6
and
depicted in light blue color. Figure 8illustrates that all these pairs of line segments intersect
on the parabola.
Mathematics 2021,9, 2198 12 of 24
Mathematics 2021, 9, x FOR PEER REVIEW 13 of 27
Figure 8. Construction of a Bernstein polynomial. The extremum of the parabola is marked with a
red dot and the control point with a blue dot. Green dots indicate points on the parabola, resp.
points resulting from extension of straight lines between these points to the boundaries of the in-
terval
02
xx
[,]
. The convex hull formed by the control points
00
(, )xy
,
(,)
xy
CC
and
22
(, )xy
is
highlighted in grey.
Considering the above, there are an infinite number of combinations defining the
term 20
()/2yy
′′
+
depending on the position of the point 11
(, )
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
2
mn
y
yy
C′′
+
=
, (44)
where m
y′
and n
y′
represent any pair of
y′
-coordinates from the mentioned line seg-
ments that intersect on the parabola. Correspondingly, the coordinate
x
C
is defined as
1
2
ii
x
xx
C
+
+
=, (45)
where
i
x
and 1i
x
+ define the ends of the interval in which the parabola varies. Hence,
x
C
is always in the middle of the interval.
There is an additional property of the line segments 0
M
n
and 6
Nm
. They inter-
sect at the point
(, )
xy
CC
and are tangential to the curve. This can be easily proven. We
can compute the first derivative and solve for m
y′
and n
y′
from
Error! Reference source not found., then construct two line equations for 00 2
(, )(, )
n
x
yxy
and 022
(, )(, )
m
x
yxy
and solve for their intersection.
In the context of geometric modelling, the point defined by
(, )
xy
CC
is called control
point and, concerning basis splines (B-splines), it is named de-Boor-point. In Error! Refer-
ence 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
Figure 8.
Construction of a Bernstein polynomial. The extremum of the parabola is marked with
a red dot and the control point with a blue dot. Green dots indicate points on the parabola, resp.
points resulting from extension of straight lines between these points to the boundaries of the interval
[x0
,
x2]
. The convex hull formed by the control points
(x0
,
y0)
,
(Cx
,
Cy)
and
(x2
,
y2)
is highlighted
in grey.
Considering the above, there are an infinite number of combinations defining the term
(y0
2+y0
0)/
2 depending on the position of the point
(x1
,
y1)
, 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
Cy=y0
m+y0
n
2, (44)
where
y0
m
and
y0
n
represent any pair of
y0
-coordinates from the mentioned line segments
that intersect on the parabola. Correspondingly, the coordinate Cxis defined as
Cx=xi+xi+1
2, (45)
where
xi
and
xi+1
define the ends of the interval in which the parabola varies. Hence,
Cx
is
always in the middle of the interval.
There is an additional property of the line segments
Mn0
and
Nm6
. They intersect
at the point
(Cx
,
Cy)
and are tangential to the curve. This can be easily proven. We can
compute the first derivative and solve for
y0
m
and
y0
n
from (44), then construct two line
equations for (x0,y0)(x2,yn)and (x0,ym)(x2,y2)and solve for their intersection.
In the context of geometric modelling, the point defined by
(Cx
,
Cy)
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
f0,2(x) = x2−x
x2−x02
y0+2x2−x
x2−x0
x−x0
x2−x0
Cy+x−x0
x2−x02
y2. (46)
As a more general representation of the Bernstein polynomial, (46) is a convex com-
bination as well. The proof is obvious. Since
x0<x2
, all functions in front of
y0
,
Cy
and
y2
are non-negative and have a sum of 1. The control points
(x0
,
y0)
,
(Cx
,
Cy)
and
(x2
,
y2)
form a convex hull, highlighted in grey in Figure 8, that encloses all realizations of
xin (46)
,
Mathematics 2021,9, 2198 13 of 24
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 [x0,x2]=[0, 1].
If we now consider the example from Table 1once again, we get, with
y0
0=
3.5 and
y0
2=7.0, from (44) the result Cy=5.25 and therefore
f0,2(x) = 3−x
3−02
·1+2·3−x
3−0·x−0
3−0·5.25 +x−0
3−02
·2, (47)
which can be “simplified” to the monomial basis form
f0,2(x) = 1+17
6x−5
6x2. (48)
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.
4. 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
Cyi
. 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
(x0
,
y0)
,
(x1
,
y1)
and
(x2
,
y2)
. 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
x0
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 (x0,y0)and (x1,y1); and
- One knot (x2,y2)is given that is to be interpolated.
Mathematics 2021, 9, x FOR PEER REVIEW 15 of 27
Figure 9. Initial situation before the construction of the second part of the spline. 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.
The task is now to construct a second parabola between the knots 11
(, )
x
y
and
22
(, )
x
y
which is tangential to the first one.
The second parabola is constructed as follows:
− The line segment
1
011
(, )(,)
y
xC xy is extended from the end point 11
(, )xy
of the pa-
rabola to the end of the interval 12
[,]
xx
, yielding the additional coordinates
2
2
(, )
y
xC
and therefore the new line segment
2
11 2
(, )(, )
y
x
yxC , see Error! Reference source not
found.;
− Between points
2
1
(, )
y
x
C
and 22
(, )xy
, the line segment
2
122
(, )(, )
y
x
Cxy
can be con-
structed, see Error! Reference source not found.;
− By taking combinations of these two line segments, within the interval 12
[,]
x
x
, an-
other parabola will be constructed, smoothly connected to the previous one. The
whole function constructed from these two smoothly connected parabolas repre-
sents a spline. Both parabolas end tangentially at the point 11
(, )
x
y
w.r.t the line
segment, defined as
11 2 2
(, )(, )
xy xy
CC CC , see Error! Reference source not found..
Figure 10. Line segment 2
11 2
(, )(, )
y
x
yxC . 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.
Figure 9.
Initial situation before the construction of the second part of the spline. 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.
The task is now to construct a second parabola between the knots
(x1
,
y1)
and
(x2
,
y2)
which is tangential to the first one.
Mathematics 2021,9, 2198 14 of 24
The second parabola is constructed as follows:
-
The line segment
(x0,Cy1)(x1,y1)
is extended from the end point
(x1
,
y1)
of the
parabola to the end of the interval
[x1
,
x2]
, yielding the additional coordinates
(x2
,
Cy2)
and therefore the new line segment (x1,y1)(x2,Cy2), see Figure 12;
-
Between points
(x1
,
Cy2)
and
(x2
,
y2)
, the line segment
(x1,Cy2)(x2,y2)
can be con-
structed, see Figure 10;
-
By taking combinations of these two line segments, within the interval
[x1
,
x2]
, 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
(x1
,
y1)
w.r.t the line segment, defined as
(Cx1,Cy1)(Cx2,Cy2), see Figure 11.
Mathematics 2021, 9, x FOR PEER REVIEW 16 of 27
Figure 11. Line segment 2
122
(, )(, )
y
x
Cxy
. 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.
Figure 12. Second parabola within the interval
12
[,]
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
1
0
1
0,1 0
10 10
() y
xx
xx
fx y C
xx xx
−
−
=+
−−
,
12
0
2
1,1
20 20
() yy
xx
xx
fx C C
xx xx
−
−
=+
−−
. (49)
By taking a combination of these two lines within the interval 01
[,]
x
x
, it follows
12
22
00 0
11 2
0,2 0
10 1010 1020 2010
()
() ()()
yy
xx xx xx
xx xx xx
f
xy C C
xx xxxx xxxx xxxx
−− −
−− −
=+ + +
−−−−−−−
.(50)
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
12
0
2
2,1
20 20
() yy
xx
xx
fx C C
xx xx
−
−
=+
−−
,
2
21
3,1 2
21 21
() y
xx xx
fx C y
xx xx
−−
=+
−−
. (51)
Analogously to Error! Reference source not found., but in this case within the in-
terval 12
[,]
x
x
, we obtain
Figure 10.
Line segment
(x1,Cy2)(x2,y2)
. 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.
Mathematics 2021, 9, x FOR PEER REVIEW 16 of 27
Figure 11. Line segment 2
122
(, )(, )
y
x
Cxy
. 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.
Figure 12. Second parabola within the interval
12
[,]
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
1
0
1
0,1 0
10 10
() y
xx
xx
fx y C
xx xx
−
−
=+
−−
,
12
0
2
1,1
20 20
() yy
xx
xx
fx C C
xx xx
−
−
=+
−−
. (49)
By taking a combination of these two lines within the interval 01
[,]
x
x
, it follows
12
22
00 0
11 2
0,2 0
10 1010 1020 2010
()
() ()()
yy
xx xx xx
xx xx xx
f
xy C C
xx xxxx xxxx xxxx
−− −
−− −
=+ + +
−−−−−−−
.(50)
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
12
0
2
2,1
20 20
() yy
xx
xx
fx C C
xx xx
−
−
=+
−−
,
2
21
3,1 2
21 21
() y
xx xx
fx C y
xx xx
−−
=+
−−
. (51)
Analogously to Error! Reference source not found., but in this case within the in-
terval 12
[,]
x
x
, we obtain
Figure 11.
Second parabola within the interval
[x1
,
x2]
, 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.
Mathematics 2021,9, 2198 15 of 24
Mathematics 2021, 9, x FOR PEER REVIEW 15 of 27
Figure 9. Initial situation before the construction of the second part of the spline. 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.
The task is now to construct a second parabola between the knots 11
(, )
x
y
and
22
(, )
x
y
which is tangential to the first one.
The second parabola is constructed as follows:
− The line segment
1
011
(, )(,)
y
xC xy is extended from the end point 11
(, )xy
of the pa-
rabola to the end of the interval 12
[,]
xx
, yielding the additional coordinates
2
2
(, )
y
xC
and therefore the new line segment
2
11 2
(, )(, )
y
x
yxC , see Error! Reference source not
found.;
− Between points
2
1
(, )
y
x
C
and 22
(, )xy
, the line segment
2
122
(, )(, )
y
x
Cxy
can be con-
structed, see Error! Reference source not found.;
− By taking combinations of these two line segments, within the interval 12
[,]
x
x
, an-
other parabola will be constructed, smoothly connected to the previous one. The
whole function constructed from these two smoothly connected parabolas repre-
sents a spline. Both parabolas end tangentially at the point 11
(, )
x
y
w.r.t the line
segment, defined as
11 2 2
(, )(, )
xy xy
CC CC , see Error! Reference source not found..
Figure 10. Line segment 2
11 2
(, )(, )
y
x
yxC . 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.
Figure 12.
Line segment
(x1,y1)(x2,Cy2)
. 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
f0,1(x) = x1−x
x1−x0
y0+x−x0
x1−x0
Cy1,f1,1(x) = x2−x
x2−x0
Cy1+x−x0
x2−x0
Cy2. (49)
By taking a combination of these two lines within the interval [x0,x1], it follows
f0,2(x) = x1−x
x1−x02
y0+x1−x
x1−x0
x−x0
x1−x0
+x−x0
x1−x0
x2−x
x2−x0Cy1+(x−x0)2
(x2−x0) (x1−x0)Cy2(50)
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
f2,1(x) = x2−x
x2−x0
Cy1+x−x0
x2−x0
Cy2,f3,1(x) = x2−x
x2−x1
Cy2+x−x1
x2−x1
y2. (51)
Analogously to (50), but in this case within the interval [x1,x2], we obtain
f1,2(x) = x2−x
x2−x1
x2−x
x2−x0
Cy1+x2−x
x2−x1
x−x0
x2−x0
+x−x1
x2−x1
x2−x
x2−x1Cy2+x−x1
x2−x12
y2. (52)
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
Cy1
and
Cy2
, 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
a0+a1x+a2x2
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
Mathematics 2021,9, 2198 16 of 24
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)
y0
to
Cy0
and
y2
to
Cy3
. 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
x−2,x−1andx3,x4, (53)
which is called knot multiplication. Additional knots with e.g.,
x−2=x−1=x0and x2=x3=x4(54)
are especially used to gain endpoint interpolation. The additional knots do not necessarily
have to be equal to
x0
and
x2
, they are just required to follow the order
x−2≤x−1≤x0
and
x2≤x3≤x4
. 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
Bi,0(x) = (1, xi≤x<xi+1,
0, otherwise, (55)
with a recursive step to compute the basis functions
Bi,d(x) = x−xi
xi+d−xi
Bi,d−1(x) + xi+d+1−x
xi+d+1−xi+1
Bi+1,d−1(x), (56)
which is summed over the k-th interval of the spline fk,d(x), such as
fk,d(x) =
k
∑
i=k−d
Cyi+dBi,d(x). (57)
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
[xk
,
xk+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
Cyi
in (57) correspond to
each basis function terminated by
Bi,0(x)
in (55) at the end of the recursion. An example
of quadratic basis functions using the knots
xk0=
0,
xk1=
1 and
xk2=
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.
Mathematics 2021,9, 2198 17 of 24
Mathematics 2021, 9, x FOR PEER REVIEW 18 of 27
with a recursive step to compute the basis functions
1
,,1 1,1
11
() () ()
iid
id id i d
id i id i
xx x x
Bx B x B x
xx x x
++
−+−
++++
−−
=+
−−
, (56)
which is summed over the k-th interval of the spline
,
()
kd
f
x, such as
,,
() ()
id
k
kd y id
ikd
fCxBx
+
=−
=
. (57)
At the beginning of the recursion Error! Reference source not found., index
d
denotes the degree of the spline and i is a computed index
,,ikd k=−
. In each re-
cursive step,
d
is lowered down for 1 until it reaches 0. At
0d=
, the recursion formula
terminates and the base case (55) is evaluated.
The index
k
is the index of the interval 1
[,]
kk
x
x
+ with
0,1, , 1kl=−
, where
l
is
the number of knots. Therefore, for
0k=
, the recursion formula
Error! Reference source not found. is computed for every i of the sum
Error! Reference source not found. and evaluated in Error! Reference source not found.
. Afterwards, for
1k=
, the same computation is performed, and so on, until the last in-
terval. The parameters
i
y
C in Error! Reference source not found. correspond to each
basis function terminated by
,0
()
i
B
x in Error! Reference source not found. at the end of
the recursion. An example of quadratic basis functions using the knots
0
0
k
x=
,
1
1
k
x=
and
2
3
k
x=
, resulting in two intervals
0k=
and
1k=
, and the knot sequence
[0,0,0,1,3,3,3]
, is shown in Error! Reference source not found..
Figure 13. Quadratic B-Spline basis functions using the knots 00
k
x
=
, 11
k
x
=
and 23
k
x
=
, result-
ing in two intervals 0k= and 1k=, and the knot sequence
[0,0,0,1,3,3,3]
. The knots are marked
with black triangles on the x-axis.
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
2d=
in (56), then we follow all recursive
steps until
0d=
for every
,
()
id
B
x, which yields the explicit expression of a generic
quadratic B-spline formula
Figure 13.
Quadratic B-Spline basis functions using the knots
xk0=
0,
xk1=
1 and
xk2=
3, resulting
in two intervals
k=
0 and
k=
1, and the knot sequence
[
0, 0, 0, 1, 3, 3, 3
]
. The knots are marked with
black triangles on the x-axis.
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
Bi,d(x)
, which yields the explicit expression of a
generic quadratic B-spline formula
Bi,2(x) = (x−xi)2
(xi+2−xi) (xi+1−xi)Bi,0(x)
+(xi+3−x)2
(xi+3−xi+1) (xi+3−xi+2)Bi+2,0(x)
+(x−xi) (xi+2−x)
(xi+2−xi) (xi+2−xi+1)+(xi+3−x) (x−xi+1)
(xi+3−xi+1) (xi+2−xi+1)Bi+1,0(x).
(58)
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
Cy0B−2,2(x) = Cy0(x−x−2)2
(x0−x−2) (x−1−x−2)B−2,0(x)
+(x1−x)2
(x1−x−1) (x1−x0)B0,0(x)
+(x−x−2) (x0−x)
(x0−x−2) (x0−x−1)+(x1−x) (x−x−1)
(x1−x−1) (x0−x−1)B−1,0(x)i.
(59)
With i=−1, the second element of the sum is
Cy1B−1,2(x) = Cy1(x−x−1)2
(x1−x−1) (x0−x−1)B−1,0(x)
+(x2−x)2
(x2−x0) (x2−x1)B1,0(x)
+(x−x−1) (x1−x)
(x1−x−1) (x1−x0)+(x2−x) (x−x0)
(x2−x0) (x1−x0)B0,0(x)i
(60)
and, with i=0, the last element of the sum is
Cy2B0,2(x) = Cy2(x−x0)2
(x2−x0) (x1−x0)B0,0(x)
+(x3−x)2
(x3−x1) (x3−x2)B2,0(x)
+(x−x0) (x2−x)
(x2−x0) (x2−x1)+(x3−x) (x−x1)
(x3−x1) (x2−x1)B1,0(x)i.
(61)
Mathematics 2021,9, 2198 18 of 24
Taking the knot multiplications (53) and (54) into account for all elements of the sums
(59)–(61) yields
Cy0B−2,2(x) = Cy0(x−x0)2
(x0−x0) (x0−x0)B−2,0(x)
+(x1−x)2
(x1−x0) (x1−x0)B0,0(x)
+(x−x0) (x0−x)
(x0−x0) (x0−x0)+(x1−x) (x−x0)
(x1−x0) (x0−x0)B−1,0(x)i,
(62)
Cy1B−1,2(x) = Cy1(x−x0)2
(x1−x0) (x0−x0)B−1,0(x)
+(x2−x)2
(x2−x0) (x2−x1)B1,0(x)
+(x−x0) (x1−x)
(x1−x0) (x1−x0)+(x2−x) (x−x0)
(x2−x0) (x1−x0)B0,0(x)i,
(63)
Cy2B0,2(x) = Cy2(x−x0)2
(x2−x0) (x1−x0)B0,0(x)
+(x2−x)2
(x2−x1) (x2−x2)B2,0(x)
+(x−x0) (x2−x)
(x2−x0) (x2−x1)+(x2−x) (x−x1)
(x2−x1) (x2−x1)B1,0(x)i.
(64)
The base case (55) evaluates only the functions in front of
B0,0(x)
, all other base
functions are not taken into consideration nor are the fractions where division by zero
occurred. Finally, by summing up all results (62), (63) and (64) for the first interval, as
stated in (57), and evaluating by (55), we obtain
f0,2(x) = 0
∑
i=−2
Cyi+2Bi,2(x)
=(x1−x)2
(x1−x0) (x1−x0)Cy0
+(x−x0) (x1−x)
(x1−x0) (x1−x0)+(x2−x) (x−x0)
(x2−x0) (x1−x0)Cy1
+(x−x0)2
(x2−x0) (x1−x0)Cy2.
(65)
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
y0
is renamed to
Cy0
. If the same
procedure is applied to
k=
1, the result is identical to (52), considering that
y2
is renamed
to Cy3.
Since spline functions
y=f(x)
are considered, only the component
Cyi
of a control
point
(Cxi
,
Cyi)
appears as a factor in front of the basis functions
Bi,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.
5. 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.
Mathematics 2021,9, 2198 19 of 24
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.
5.1. Definition of the Problem
Let:
-y0,y1, . . . , yj, . . . , ymbe a sequence of observed values;
-ΣLL =σ2
0QLL be the variance-covariance matrix of the observations yj;
-x0<x1<. . . <xj<. . . <xm
be a sequence of non-decreasing error-free values
referring to the observations;
-xk0<. . . <xkn−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:
-Cy0
,
Cy1
,
. . .
,
Cyn
, 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 overde-
termined configuration using the Gauss-Markov model.
5.2. 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
xj
belongs. One has to consider that the
knot sequence must be given in advance.
Considering
xj
as fixed values, yields also
Bi,d(x)
, computed from (55), (56), and
(57)
,
as fixed values. Here, we assume that
Bi,2(x)
corresponds to the given parameter
Cyi
.
Explaining the observation equations explicitly by means of the Cox-de Boor’s formula is
too cumbersome. Hence, with yjas observations, the observation equations
y0+v0=n
∑
i=1
ˆ
CyiBi,2(x0),
y1+v1=n
∑
i=1
ˆ
CyiBi,2(x1),
.
.
.
ym+vm=n
∑
i=1
ˆ
CyiBi,2(xm)
(66)
can be set up. The observation vector
L=y0y1y2· · · ymT(67)
contains all observations, while the vector of unknown parameters can be written as
ˆ
X=ˆ
Cy0ˆ
Cy1· · · ˆ
CynT. (68)
From the stochastic model
ΣLL =σ2
0QLL, (69)
with σ2
0as theoretical variance factor, the corresponding weight matrix is obtained from
P=Q−1
LL , (70)
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
L+v=Aˆ
X, (71)
Mathematics 2021,9, 2198 20 of 24
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
A=
B0,2(x0)B1,2(x0)· · · Bn,2(x0)
B0,2(x1)B1,2(x1)· · · Bn,2(x1)
.
.
..
.
.....
.
.
B0,2(xm)B1,2(xm)· · · Bn,2(xm)
. (72)
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
xj
that belong to a particular interval and everything else is
equal to zero.
5.3. Least Squares Adjustment
One feature of the B-splines is that constraints for continuity, smoothness and con-
tinuity 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 un-
knowns can be obtained from
ˆ
X= (ATPA)−1ATPL. (73)
With the residuals
v=Aˆ
X−L(74)
from (71), the a posteriori estimate of the reference standard deviation
s0=rvTPv
r(75)
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.
6. 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
Mathematics 2021,9, 2198 21 of 24
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
yi
as observations and
xi
as error-free values, the resulting
linear adjustment problem could be solved within the Gauss-Markov model. Finally,
in Appendix Ait 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
yi
and
xi
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.
Acknowledgments: The authors acknowledge support from the German Research Foundation and
the Open Access Publication Fund of TU Berlin.
Conflicts of Interest: The authors declare no conflict of interest.
Appendix A. Transition from B-spline to Ordinary Polynomial
Through Sections 2–4, it was explained how a constituent part of the B-spline can be
derived from an ordinary second-degree polynomial. In Appendix A, the parameters of an
ordinary second-degree spline are to be expressed as functions of the B-spline parameters.
We consider the following case, cf. Figure 11:
-
Three given points
(x0
,
y0)
,
(x1
,
y1)
,
(x2
,
y2)
are to be interpolated by a quadratic spline;
-
using the given B-spline parameters
Cy1
,
Cy2
, i.e., the y-component of the control points.
Since the parameters of each polynomial of an ordinary second-degree spline are
derived in the same way, in the text that follows, only the derivation of the first polynomial
that varies within the interval [x0,x1]is described.
From the line segments
(x0,y0)(x1,Cy1)
and
(x0,Cy1)(x2,Cy2)
, see Figure 11, two line
equations of the form
f0,1(x) = a0+a1x
resp.
f1,1(x) = b0+b1x
can be derived. From (49),
the slope for the first line equation is
a1=Cy1−y0
x1−x0
(A1)
and the y-intercept is
a0=y0−x0
Cy1−y0
x1−x0
. (A2)
Analogously, the parameters for the second line equation are
b1=Cy2−Cy1
x2−x0
(A3)
and
b0=Cy1−x0
Cy2−Cy1
x2−x0
. (A4)
Mathematics 2021,9, 2198 22 of 24
The spline segment of the form of a second degree polynomial is derived by taking
a combination
f0,2(x) = x1−x
x1−x0
(a0+a1x) + x−x0
x1−x0
(b0+b1x)(A5)
of these two lines in the same way as in (50). Since
x
is the variable, after some rearrange-
ment, the coefficients in front of the variable can be expressed as
f0,2(x) = a0x1−b0x0
x1−x0
+b0−a0+a1x1−b1x0
x1−x0x+b1−a1
x1−x0x2. (A6)
The derived coefficients of the polynomial that varies within the interval
[x0
,
x1]
read
A0=a0x1−b0x0
x1−x0
,A1=b0−a0+a1x1−b1x0
x1−x0
,A2=b1−a1
x1−x0
. (A7)
By inserting the terms from (A7) into (A6), we obtain the ordinary second-
degree polynomial
f0,2(x) = A0+A1x+A2x2. (A8)
Only the interval
[x0
,
x1]
is considered as the domain of the definition of the spline.
The parameters of the ordinary second-degree polynomial, which varies within the interval
[x1,x2], can be derived analogously to the first one and would take the form
f1,2(x) = B0+B1x+B2x2. (A9)
The spline function expressed by these two polynomials is identical with the one
expressed by a B-spline.
A more general derivation of the B-spline representation of polynomials was devel-
oped by Lyche et al. [
37
] (pp. 15–18), and it was shown that polynomials can be represented
in terms of B-splines of at least the same degree.
Appendix B. Transition from Ordinary Polynomial to B-spline
In Appendix A, we showed how second-degree B-spline parameters can be used to
derive the parameters for an equivalent representation in the form of ordinary second-
degree polynomials. In Appendix B, we explain the reverse procedure: how the ordinary
quadratic spline parameters
A0
,
A1
,
A2
, respectively
B0
,
B1
,
B2
, can be used for derivation
of the B-spline parameters y0,C1,C2and y2.
With Equation (A8), we can express each parameter of the ordinary cubic polynomial
as follows:
A0(x1−x0) = y0 x2
1
x1−x0!+C1−x0x1
x1−x0
−x0x2
x2−x0+C2 x2
0
x2−x0!, (A10)
A1(x1−x0) = y0−2x1
x1−x0+C12x0
x2−x0
+2x1
x1−x0+C2−2x0
x2−x0, (A11)
A2(x1−x0) = y01
x1−x0+C1−1
x2−x0
−1
x1−x0+C21
x2−x0. (A12)
According to this approach, the parameters in (A9) can be expressed as:
B0(x2−x1) = C1 x2
2
x2−x0!+C2x0x2
x2−x0
−x1x2
x2−x1+y2 x2
1
x2−x1!, (A13)
B1(x2−x1) = C1−2x2
x2−x0+C22x1
x2−x1
+2x2
x2−x0+y2−2x1
x2−x1, (A14)
Mathematics 2021,9, 2198 23 of 24
B2(x2−x1) = C11
x2−x0+C2−1
x2−x1
−1
x2−x0+y21
x2−x1. (A15)
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
y0=A0+A1x0+A2x2
0, (A16)
C1=A0+A1
x0+x1
2+A2x0x1, (A17)
C2=A0+A1
x1+x2
2+A2x1x2(A18)
and for the system (A13)–(A15), the results are
C1=B0+B1
x0+x1
2+B2x0x1, (A19)
C2=B0+B1
x1+x2
2+B2x1x2, (A20)
y2=B0+B1x2+B2x2
2. (A21)
In the terms of B-spline,
y0
and
y2
are parameters and, in this derivation, they are
treated as such. If
y0
and
y2
are considered as known values, we can also solve only for
C1
and
C2
as a system of two equations. However, in such case, the solutions for
C1
and
C2
are not as elegant as those from the equations above.
References
1.
Ezhov, N.; Neitzel, F.; Petrovic, S. Spline Approximation—Part 1: Basic Methodology. J. Appl. Geod.
2018
,12, 139–155. [CrossRef]
2. Piegl, L.; Tiller, W. The NURBS Book, 2nd ed.; Springer: Berlin/Heidelberg, Germany, 1997.
3. de Boor, C. On calculating with B-splines. J. Approx. Theory 1972,6, 50–62. [CrossRef]
4. Cox, M. The numerical evaluation of B-splines. IMA J. Appl. Math. 1972,10, 134–149. [CrossRef]
5. Lowther, J.; Shene, C.-K. Teaching B-splines is not difficult! ACM SIGCSE Bull. 2003,35, 381–385. [CrossRef]
6.
Berkhahn, V. Geometric Modelling in Civil Engineering Informatics. Habilitation Thesis, University of Hannover, Hanover,
Germany, Shaker Verlag, Aachen, Germany, 2005. (In German).
7. Farin, G. Curves and Surfaces for CAGD: A Practical Guide, 5th ed.; Morgan Kaufmann Publishers: San Francisco, CA, USA, 2002.
8. Lobatschewsky, N. Probabilitédes résultats moyens tirés d’observations répetées. J. Die Reine Angew. Math. 1842,24, 164–170.
9.
Schoenberg, I.J. Contributions to the problem of approximation of equidistant data by analytic functions, Part A: On the problem
of smoothing of graduation, a first class of analytic approximation formulae. Quart. Appl. Math. 1946,4, 45–99. [CrossRef]
10.
Schoenberg, I.J. Contributions to the problem of approximation of equidistant data by analytic functions, Part B: On the problem
of osculatory interpolation, a second class of analytic approximation formulae. Quart. Appl. Math. 1946,4, 112–141. [CrossRef]
11. Schumaker, L. Spline Functions: Basic Theory, 3rd ed.; Cambridge University Press: Cambridge, UK, 2007.
12.
de Boor, C.; Pinkus, A. The B-spline recurrence relations of Chakalov and of Popoviciu. J. Approx. Theory
2003
,124, 115–123.
[CrossRef]
13. Szilvasi-Nagy, M. Almost curvature continuous fitting of B-spline surfaces. J. Geom. Graph. 1998,2, 33–43.
14.
Zheng, W.; Bo, P.; Liu, Y.; Wang, W. Fast B-spline curve fitting by L-BFGS. Comput. Aided Geom. Des.
2012
,29, 448–462. [CrossRef]
15. Wold, S. Spline functions in data analysis. Technometrics 1974,16, 1–11. [CrossRef]
16.
Anderson, I.; Turner, D. Discrete B-spline approximation in a variety of norms. In Advanced Mathematical and Computational Tools
in Metrology V; Ciarlini, P., Cox, M.G., Filipe, E., Pavese, F., Richter, D., Eds.; World Scientific Publishing Company: Singapore,
2001; pp. 1–15.
17.
Sünkel, H. Reconstruction of functions from discrete mean values using bicubic spline functions. In Mitteilungen der Geodätischen
Institute der Technischen Universität Graz; Contributions of the Graz Group to the XVI General Assembly of IUGG/IAG in Grenoble
1975, Folge 20; Meissl, P., Moritz, H., Rinner, K., Eds.; Graz University of Technology: Graz, Austria, 1975; pp. 267–304.
18.
Sünkel, H. The representation of geodetic integral formulae by bicubic spline functions. In Mitteilungen der Geodätischen Institute
der Technischen Universität Graz; Folge 28; Graz University of Technology: Graz, Austria, 1977. (In German)
19.
Sünkel, H. Local interpolation of gravity data using spline functions. In Tagungsbericht über das 1. Alpengravimetrie-Kolloquium;
Steinhauser, P., Ed.; Zentralanstalt für Meteorologie und Geodynamik: Vienna, Austria, 1977; pp. 83–101. (In German)
Mathematics 2021,9, 2198 24 of 24
20.
Sünkel, H. A General Surface Representation Module Designed for Geodesy (GSPP); Report No. 292; Department of Geodetic Science,
The Ohio State University: Columbus, OH, USA, 1980.
21.
Schoenberg, I.J. Cardinal Spline Interpolation; Regional Conference Series in Applied Mathematics, No. 12; SIAM: Philadelphia,
PA, USA, 1973.
22.
Schumaker, L.L. Fitting surfaces to scattered data. In Approximation Theory II; Lorentz, G.G., Chui, C.K., Schumaker, L.L., Eds.;
Academic Press: New York, NY, USA, 1976; pp. 203–268.
23.
Späth, H. Spline Algorithms for Construction of Smooth Curves and Surfaces; Oldenbourg: München-Vienna, Austria, 1973.
(In German)
24.
Jekeli, C. Spline Representations of Functions on a Sphere for Geopotential Modeling; Report No. 475; Department of Civil and
Environmental Engineering and Geodetic Science, The Ohio State University: Columbus, OH, USA, 2005.
25.
Mautz, R.; Ping, J.; Heki, K.; Schaffrin, B.; Shum, C.; Potts, L. Efficient spatial and temporal representations of global ionosphere
maps over Japan using B-spline wavelets. J. Geod. 2005,78, 660–667. [CrossRef]
26.
Koch, K.R.; Schmidt, M. N-dimensional B-spline surface estimated by lofting for locally improving IRI. J. Geod. Sci.
2011
,1, 41–51.
[CrossRef]
27.
Dzapo, M.; Junasevic, M.; Lapaine, M.; Petrovic, S. Use of splines for determination of lengths from starting point along railway
lines. Geod. List. 1985,39, 99–105. (In Croatian)
28.
Bureick, J.; Neuner, H.; Harmening, C.; Neumann, I. Curve and surface approximation of 3D point clouds. Avn Allg. Vermess.
Nachr. 2016,123, 315–327.
29.
Kermarrec, G.; Hartmann, J.; Faust, H.; Hartmann, K.; Besharat, R.; Samuel, G.; Lixian, C.; Alkhatib, H. Understanding hierarchical
B-splines with a case study: Approximation of point clouds from TLS observations. Z. Geodäsie Geoinf. Landmanagement
2020
,145,
224–235.
30.
Harmening, C. Spatio-Temporal Deformation Analysis Using Enhanced B-Spline Models of Laser Scanning Point Clouds.
Ph.D. Thesis, Technische Universität Wien, Faculty of Mathematics and Geoinformation, Vienna, Austria, 2020.
31. Gander, W. Change of basis in polynomial interpolation. Numer. Linear Algebra Appl. 2005,12, 769–778. [CrossRef]
32.
Bronshtein, I.N.; Semendyayev, K.A.; Musiol, G.; Muehlig, H. Handbook of Mathematics, 5th ed.; Springer: Berlin/Heidelberg,
Germany, 2007.
33. Rockafellar, R.T. Convex Analysis; Princeton Mathematical Series 28; Princeton University Press: Princeton, NJ, USA, 1970.
34.
Bureick, J.; Alkhatib, H.; Neumann, I. Robust spatial approximation of laser scanner point clouds by means of free-form curve
approaches in deformation analysis. J. Appl. Geod. 2016,10, 27–35. [CrossRef]
35. Niemeier, W. Adjustment Computation, Statistical Evaluation Methods, 2nd ed.; De Gruyter: Boston, MA, USA, 2008. (In German)
36. Neitzel, F.; Ezhov, N.; Petrovic, S. Total least squares spline approximation. Mathematics 2019,7, 462. [CrossRef]
37.
Lyche, T.; Manni, C.; Speleers, H. B-Splines and Spline Approximation. 2017. Available online: https://www.mat.uniroma2.it/
~{}speleers/cime2017/material/notes_lyche.pdf (accessed on 9 July 2021).
