B-Spline Interpolation and Approximation

1. Parameter Selection and Knot Vector Generation

i. Overview

ii. The Uniformly Spaced Method

iii. The Chord Length Method

iv. The Centripetal Method

v. Knot Vector Generation

vi. The Universal Method

vii. Parameters and Knot Vectors for Surfaces

2. Solving Systems of Linear Equations

3. Curve Interpolation: Global Interpolation
4. Curve Approximation: Global Approximation
5. Surface Interpolation: Global Interpolation
6. Surface Approximation: Global Approximation

1
Parameter Selection Overview

The input to a B-spline interpolation/approximation algorithm usually consists of a set of data points. Thus,
the first step is to find a set of parameters that can "fixed" these points at certain values. More precisely, if
the data points are D0, ..., Dn, then n+1 parameters t0, ..., tn in the domain of the curve must be found so that
data point Dk corresponds to parameter tk for k between 0 and n. This means that if C(u) is a curve that
passes through all data points in the given order, then we have Dk = C(tk) for all 0 £ k £ n. In the figure
below, we have seven data points (i.e., n = 6) and seven parameters must be found to setup the desired
correspondence.
There are infinite number of possibilities for selecting these
parameters. For example, we can evenly divide the domain,
or randomly pick n+1 values from the domain. However,
poorly chosen parameters cause unpredictable results. The
following figure shows four data points and three
interpolating curves, each of which is obtained using a
different set of parameters. One of them bends outward too
much and creates an unnecessary bulge. The black curve has
a peak and a small bulge. Only the one between these two
follows the trend of data points closely. Thus, the choice of
parameters affects the shape of the curve, and, consequently,
affects the parameterization of the curve.

On the next few pages, we shall discuss a number of parameter

selection methods, which include the Uniformly Spaced Method,
Chord Length Method and Centripetal Method. Once a set of
parameters is obtained, we need to compute a knot vector. See
Knot Vector Generation for the details. We will also discuss the
Universal Method in which a uniformly spaced knot vector is
chosen and then used for computing the parameters. Please keep
in mind that there are other methods for selecting
parameters.

2
The Uniformly Spaced Method

The simplest parameter selection method is the uniformly spaced method. Suppose the domain, as usual,
is [0,1] and n+1 uniformly spaced parameters are required. The first and the last parameters must be 0 and 1
because we want the curve to pass through the first and the last data points. Therefore, we have t0 = 0 and
tn=1.
Since n+1 points divide the interval [0,1] into n subintervals evenly, each of which must be of length 1/n,
the division points are 0, 1/n, 2/n, 3/n, ..., (n-1)/n and 1. Thus, we have

ì t0 = 0
ï
ï i
í ti = for 1 £ i £ n - 1
ï n
ïîtn = 1

For example, if we need 5 parameters, n = 4, and the uniformly spaced parameters are 0, 1/4, 1/2, 3/4 and 1.
If we need 8 parameters, n = 7, and the uniformly spaced parameters are 0, 1/7, 2/7, 3/7, 4/7, 5/7, 6/7 and 1.
What if the domain is [a,b] rather than [0,1]? In this case, [a,b] is divided into n intervals by n+1 division
point, a and b included. Since the length of this interval is b - a, the length of each subinterval is (b-a)/n.
Hence, the division points (i.e., parameters) are

ì t0 = a
ï
ï b-a
í ti = a + i for 1 £ i £ n - 1
ï n
ïîtn = b

While the uniformly spaced method is simple, it is known to generate some unpleasant results. For example,
when data points are not uniformly spaced, using uniformly spaced parameters could generate erratic
shapes such as big bulges, sharp peaks and loops. In Figure (a) below, there is a loop at data point 3. In
Figure (b), the curve wiggles its way through data points 1, 2 and 3. While we cannot say that these
problems are unique to the uniformly spaced method, it does occur more frequently than with other
methods.

(a) (b)

3
The Chord Length Method

If an interpolating curve follows very closely to the data polygon, the length of the
curve segment between two adjacent data points would be very close to the length
of the chord of these two data points, and the length of the interpolating curve
would also be very close to the total length of the data polygon. In the figure below,
each curve segment of an interpolating polynomial is very close to the length of its
supporting chord, and the length of the curve is close to the length of the data
polygon. Therefore, if the domain is subdivided according to the distribution of the
chord lengths, the parameters will be an approximation of the arc-length
parameterization. This is the merit of the chord length or chordal method.
Suppose the data points are D0, D1, ..., Dn. The length between Di-1 and Di is |Di-Di-1|,
and the length of the data polygon is the sum of the lengths of these chords:
n
L = å Di - Di -1
i =1

Therefore, the ratio of the chord length from data point D0 to data point Dk, denoted as Lk, over the length of
the data polygon is
k
Lk = å Di - Di -1
i =1

If we prefer to have an arc-length parameterization of the interpolating curve, the domain has to be divided
according to the ratio Lk. More precisely, if the domain is [0,1], then parameter tk should be located at the
value of Lk:

t0 = 0
k

åD i - Di -1
tk = i =1

L
tn = 1

where L is the length of the data polygon. In this way, the parameters divide the domain into the ratio of the
chord lengths.
Let us look at an example. Suppose we have four data points (n = 3): D0 = < 0,0 >, D1 = < 1,2 >, D2 = < 3,4
> and D3 = < 4,0 >. The length of each chord is

D1 - D0 = 5 = 2.236
D2 - D` = 2 2 = 2.828
D3 - D2 = 17 = 4.123

and the total length is

L = 5 + 2 2 + 17 = 9.8176
Finally, we have the corresponding parameters:

4
 t0 = 0
D1 - D0
t1 = = 0.2434
L
D1 - D0 + D2 - D1
t2 = = 0.5512
L
t3 = 1

The following figures show the data points and the parameter distributions of the uniformly spaced method
and the chord length method.

What if the domain is [a,b] rather than [0,1]? Note that Lk is a ratio between 0 and 1. Since the length of [a,
b] is b-a, Lk(b-a), 0£k£ n, divide [0,b-a] into the same ratio as we did for [0,1]. Therefore, the following
parameters divide [a,b] according to the chord lengths:

t0 = a
ti = a + Lk (b - a )
tn = b

The chord length method is widely used and usually performs well. Since it is known (proved by R. Farouki)
that polynomial curves cannot be parameterized to have unit speed (i.e., arc-length parameterization), the
chord length can only be an approximation. Sometimes, a
longer chord may cause its curve segment to have a bulge
bigger than necessary. In the figure below, the black and
blue curves both interpolate 7 data points. As you can see,
both curves have very similar shape, except for the last
segment, and the one computed with chord length method
wiggling a little. The last curve segments are very different
and the curve using the chord length method has a large
bulge and twists away from the red curve produced by the
uniformly spaced method. This is a commonly seen
problem with the chord length method.

5
The Centripetal Method

E. T. Y. Lee proposed this centripetal method. Suppose we are driving a car through a slalom course. We
have to be very careful at sharp turns so that the normal acceleration (i.e., centripetal force) should not be
too large. Otherwise, our car may lose control. To drive safely, Lee suggested that the normal force along
the path should be proportional to the change in angle. The centripetal method is an approximation to this
model. In fact, we can consider the centripetal method as an extension to the chord length method.
Suppose the data points are D0, D1, ..., Dn. First of all, we should select a positive "power" value a. Usually,
it is a = 1/2 for square root. Then, the distance between two adjacent data points is measured by |Dk - Dk-1|a
rather than the conventional |Dk - Dk-1|. The length of the data polygon under this new measure is
n
L = å Di - Di -1
a

i =1

The ratio of the distance from D0 to Dk on the data polygon over the total length is
k

åD
a
i - Di -1
Lk = i =1

Therefore, L0 = 0, L1, ..., Ln = 1 divide [0,1] according to the length of data polygon under the new distance
measure. Hence, the parameters are

t0 = 0
k

åD
a
i - Di -1
tk = i =1

L
tn = 1

If a = 1, the centripetal method reduces to the chord length method, and, hence, the former can be
considered as an extension to the latter. If a < 1, say a = 1/2 (i.e., square root), |Dk - Dk-1|a is less than |Dk -
Dk-1|. Consequently, the impact of a longer chord (i.e., length > 1) on the length of the data polygon is
reduced, and the impact of a shorter chord (i.e., length < 1) on the length of the data polygon is increased.
Because of this characteristic, Lee claimed that centripetal method handle sharp turns better than the chord
length method.
Let us redo the example discussed in the chord length method. We have four data points (n = 3): D0 = <
0,0 >, D1 = < 1,2 >, D2 = < 3,4 > and D3 = < 4,0 >. Choose a = 1/2 and the length of each chord under this
new measure is
1/ 2
D1 - D0 = 5 = 1.495
1/ 2
D 2 - D` = 2 2 = 1.682
1/ 2
D3 - D 2 = 17 = 2.031

The total length is

L= 5+ 2 2+ 17 = 5.208

6
Therefore, the parameters are
t0 = 0
1/ 2
D1 - D0
t1 = = 0.2871
L
1/ 2 1/ 2
D1 - D0 + D2 - D1
t2 = = 0.6101
L
t3 = 1

The following gives the distributions of three sets of parameters computed using the uniformly spaced
method, the chord length method, and the centripetal method.

Let us take a look at an extreme example. The following

shows four data points interpolated by three B-spline curves
using the uniformly spaced method (in black), chord length
method (in blue) and centripetal method (in red). As you
can see, the uniformly spaced method has a peak, the chord
length method have two big bulges, and the centripetal
method interpolates the two very close adjacent points
nicely.

Can we say the centripetal method is superior to the other

two? The following example tells a different story. We
have 7 data points and the black, blue and red curves are
generated using the uniformly spaced method, chord
length method and centripetal method. Now, the uniformly
spaced method provides a very tight interpolation. The
centripetal method is slightly off the the tight result using
the uniformly spaced method. Finally, the chord length
method wiggles through the two longest chords too much!

7
Knot Vector Generation

Once a set of parameters is obtained, we can generate a knot vector. Suppose we have n+1 parameters t0,
t1, ..., tn and degree p. For a B-spline curve of degree p, we need m+1 knots, where m=n+p+1. If the curve is
clamped, these knots are u0 = u1 = .... = up = 0, up+1, ..., um-p-1, um-p = um-p+1 = ... = um = 1. More precisely, the
first p+1 and last p+1 knots are 0's and 1's, respectively. For the remaining n-p knots, they can be uniformly
spaced or chosen properly to achieve some desired conditions.
Suppose the remaining n-p internal knots are uniformly spaced. Then, up = 0, up+1, ..., um-p-1, um-p = 1 divide
[0,1] into n-p+1 subintervals. Therefore, the knots are
u0 = u1 = L = u p = 0
j
u j+ p = for j = 1, 2,L , n - p
n - p +1
um - p = um - p +1 = L = um = 1

For example, if we have 6 (n = 5) parameters and the degree is p = 3, we should find (n+p+1)+1 =
(5+3+1)+1 = 10 knots (i.e., m=9). Since we use clamped curves, the first and the last four (i.e., p+1) knots
should be equal to 0 and 1, respectively. More precisely, we have 0, 0, 0, 0, u4, u5, 1, 1, 1 and 1, and the two
internal knots divide [0,1] into three subintervals, each of which has length 1/3. Hence, the knot vector is
{ 0, 0, 0, 0, 1/3, 2/3, 1, 1, 1, 1}.
A uniformly spaced knot vector does not require the knowledge of the parameters, and is very simple to
generate. However, this method is not recommended because, if it is used with the chord length method for
global interpolation, the system of linear equations would be singular.
Another popular method for generating knot vector, suggested by de Boor, is to "average" the parameters.
Here is the computation formula:

u0 = u1 = L = u p = 0
j + p -1
1
u j+ p =
p
åt
i= j
i for j = 1, 2,L , n - p

um - p = um - p +1 = L = um = 1

Thus, the first internal knot is the average of p parameters t1, t2, ..., tp; the second internal knot is the average
of the next p parameters, t2, t3, ..., tp+1. Suppose the 6 (n = 5) parameters are

t0 t1 t2 t3 t4 t5

0 1/4 1/3 2/3 3/4 1

Suppose we need to generate a vector for degree p = 3. As mentioned earlier, 10 knots (m = 9) are required.
Of these 10, the first four and last four knots are 0's and 1's, respectively. Thus, the first internal knot is the
average of parameters 1/4, 1/3 and 2/3, and the second internal knot is the average of the next three
parameters 1/3, 2/3 and 3/4. The computed knot vector is

8
 u0 = u1 = u2 = u3 u4 u5 u6 = u7 = u8 = u9

0 (1/4+ 1/3 + 2/3)/3 = 5/12 (1/3 + 2/3 + 3/4)/3 = 7/12 1

The following diagram illustrates the positions of parameters and the generated knots. Note that 0(4) and
1(4) means 0 and 1 are quadruple knots (i.e., multiplicity = 4). As you can see, the first non-zero knot span
[0,5/12) contains two parameters, the second non-zero knot span [5/12,7/12) contains no parameter, and the
third non-zero knot span [7/12,1) contains two parameters.

9
The Universal Method

In 1999 Choong-Gyoo Lim proposed an interesting method for selecting parameters. In previously
discussed methods, we determine the parameters and then generate a knot vector. Lim's method, known as
the universal method, works the other way around by using uniformly spaced knots for computing the
parameters.
Suppose we need n+1 parameters, one for each data points, and the degree of the desired interpolating
B-spline curve is p. Then, the number of knots is m+1, where m = n + p + 1. Lim suggested that these knots
should be uniformly spaced. More precisely, the first p+1 knots are set to 0, the last p+1 knots are set to 1,
and the remaining n - p knots evenly subdivide the domain [0,1]. Therefore, the knots are:
u0 = u1 = L = u p = 0
i
u p+i = for i = 1, 2,L , n - p
n - p +1
um - p = um - p +1 = L = um = 1

This set of clamped type knots defines n+1 B-spline basis functions. Then, the parameters are chosen to be
the values at which their corresponding basis functions reach a maximum. Take a look at the following
figure, where n = 6 (i.e., 7 data points), p = 4, and m = 11 (i.e., 12 knots). Because we use clamped knots, 0
and 1 are two knots of multiplicity 5 (i.e., p+1) and there are only two internal knots at 1/3 and 2/3. Note
that there are n+1 B-spline basis functions, and the maximums of the first and the last are 0 and 1,
respectively. The other maximums are marked with blue bars, and the parameters are marked with yellow
dots.

Let us do a hand calculation example. Suppose we have 4 data points (i.e., n=3) and degree p=2. Thus, the
number of knots is 7 (i.e., m = n + p + 1 = 6). Since the knots are uniformly spaced, they are

u0 = u1 = u2 u3 u4 = u5 = u6

0 0.5 1

Then, we can compute the desired B-spline Basis functions. Let us start with degree 0 (i.e., p = 0). Click
here to review the definition of B-spline basis functions, and here for computation examples. The
computed results are:

Knot span Basis Function

[u0,u1) = [0,0) N0,0(u) = 0

10
 [u1,u2) = [0,0) N1,0(u) = 0

[u2,u3) = [0,0.5) N2,0(u) = 1

[u3,u4) = [0.5,1) N3,0(u) = 1

[u4,u5) = [1,1) N4,0(u) = 0

[u5,u6) = [1,1) N5,0(u) = 0

Next, we compute basis functions of degree 1. Since N0,0(u) and N1,0(u) are both 0, N0,1(u) is zero
everywhere. Similarly, N4,1(u) is also zero everywhere. Therefore, for basis functions of degree 1, we only
need to compute N1,1(u), N2,1(u) and N3,1(u) as shown below:

Basis Function Equation Non-zero Range

N0,1(u) 0 everywhere

N1,1(u) 1-2u [0,0.5)

2u [0,0.5)
N2,1(u)
2(1-u) [0.5,1)

N3,1(u) 2u-1 [0.5,1)

N4,1(u) 0 everywhere

The basis functions of degree 2 are computed as follows:

Basis Function Equation Non-zero Range

N0,2(u) (1-2u)2 [0,0.5)

2u(2-3u) [0,0.5)
N1,2(u)
2(1-u)2 [0.5,1)

2u2 [0,0.5)
N2,2(u)
-2(1-4u+3u2) [0.5,1)

N3,2(u) (2u-1)2 [0.5,1)

The following figure shows all four B-spline basis functions of degree 2.

11
It is not difficult to verify that the maximums of N0,2(u), N1,2(u), N2,2(u) and N3,2(u) are 1 at u = 0, 2/3 at u =
1/3, 2/3 at u = 2/3 and 1 at u = 1. Therefore, with the universal method, the knot vector is { 0, 0, 0, 0.5, 1, 1,
1} and the parameter vector is { 0, 1/3, 2/3, 1}. We have a uniformly spaced parameter in this case.
However, the parameters are not necessarily uniformly spaced in general, although the distribution of
parameters is usually very even. Because of this fact, the universal method performs like the uniformly
spaced method.
Because the knots are uniformly spaced, there is no multiple knots and all basis functions are of Cp-1
continuous. Click here to learn more about continuity of B-spline basis functions. Therefore, in practice,
the maximum of a B-spline basis function does not have to be computed precisely. Our experience shows
that sampling some values in the non-zero domain and choosing the one with maximum function value
usually provides rather accurate result. If accurate maximums are desirable, one-dimensional search
techniques such as the Golden Section Search can be used. The search domain is, of course, the non-zero
domain of a B-spline basis function.
The universal method has a very interesting property. That is, it is affine invariant. This means, the
transformed interpolating B-spline curve can be obtained by transforming the data points. This is similar to
the affine invariance property of B-spline curves. In fact, it can easily be shown that if the same set of knots
and parameters are used in the original and transformed interpolating B-spline curves, then transforming the
curve can be achieved by transforming the data points. From this result, we immediately know that the
uniformly spaced method is affine invariant, because the knot vector is computed from a set of uniformly
spaced parameters which are not changed before and after a transformation. However, the chord length
method and centripetal method are not affine invariant, because in the transformed curve, the distribution of
chord lengths may not be the same as the original. As a result, after transforming the data points, we have a
new set of chord lengths from which a new set of parameters must be computed.

12
Parameters and Knot Vectors for
Surfaces

Since a B-spline surface of degree (p,q) defined by e+1 rows and f+1 columns of control points has the
following equation
e f
S(u , v) = åå Pij N i , p (u ) N j , q (v)
i =0 j = 0

it requires two sets of parameters for surface interpolation and approximation. Suppose we have m+1 rows
and n+1 columns of data points Dij, where 0£ i£ m and 0£ j£ n. Thus, m+1 parameters s0, ..., sm in the
u-direction (i.e., one for each row of data points) and n+1 parameters t0, ..., tn in the v-direction (i.e., one for
each column of data points) are required so that point (sc,td) in the domain of the surface corresponds to
point S(sc,td) on the surface, which, in turn, corresponds to the data point Dcd. Putting this in an equation
form, the point on the surface, evaluated at (sc,td) as shown below,
e f
S( sc , td ) = åå Pij N i , p ( sc ) N j , q (td )
i =0 j = 0

corresponds to the data point Dcd, where sc and td are parameters in the u- and v-direction, respectively.
In the equation of the desired B-spline surface, the u-direction corresponds to the index i in Ni,p(u) and Pij.
Since i runs from 0 to m, N0,p(u), N1,p(u), ..., Nm,p(u) are coefficients of the columns of the control points.
Therefore, in the u-direction, we need m+1 parameters, and using a parameter computation method that
takes degree p and the data points on column j, we can compute m+1 parameters u0,j, u1,j, ..., um,j. This is
shown in the figure below. The desired parameters s0, s1, ..., sm are simply the average of each row. More
precisely, parameter si is the average of parameters on row i, si = (ui,0 + ui,1 + ... + ui,n)/(n+1).

The computation of parameters for the v-direction is similar. Each row of data points has n+1 points and
hence requires n+1 parameters. Thus, for data points on row i, we can compute n+1 parameter values vi,0,
vi,1, ..., vi,n. Since we have m+1 rows, these values can be organized into a (m+1)×(n+1) matrix, and
parameter tj is simply the average of the parameters on column j, tj = (v0,j + v1,j + ... + vm,j)/(m+1). In this
way, we obtain n+1 parameters for the v-direction.
This procedure is summarized as follows:

13
 Input: (m+1)×(n+1) data points Di,j;
Output: Parameters in the u-direction s0, ..., sm
and parameters in the v-direction t0, ..., tn;
Algorithm:
/* calculate s0, ..., sm */
for j := 0 to n do /* column j */
Compute a set of m+1 parameters u0,j, u1,j, ..., um,j;
for i := 0 to m do
si = (ui,0 + ui,1 + ... + ui,n)/(n+1);
/* parameters for the u-direction are available */

/* calculate t0, ..., tn */

for i := 0 to m do /* row i */
Compute a set of n+1 parameters vi,0, vi,1, ..., vi,n;
for j := 0 to n do
tj = (v0,j + v1,j + ... + vn,j)/(m+1);
/* parameters for the v-direction are available */

Then, parameter values s0, s1, ..., sm and degree p jointly define a knot vector U in the u-direction. Similarly,
parameters t0, t1, ..., tn and degree q jointly define a knot vector V in the v-direction. Click here to review
the way of computing a knot vector from a set of parameters.
Note that the above is only a conceptual algorithm and is inefficient. Once you know what is going on, it is
quite easy to convert it to a very efficient algorithm. Note also that this algorithm only works for the
uniformly spaced, chord length and centripetal methods. For the universal method, because there is no data
points involved, we can apply uniform knots to one row and one column of data points for computing the
parameters.

14
Solving Systems of Linear Equations

Since solving a system of linear equations is a basic skill that will be used for interpolation and
approximation, we will briefly discuss a commonly used technique here. In fact, what we will be using is a
slightly more general form. Suppose we have a n×n coefficient matrix A, a n×h "constant" term matrix B,
and a n×h unknown matrix X defined as follows:
é b11 b12 L b1h ù é x11 x12 L x1h ù é a11 a12 L a1n ù
êb b L b2 h úú êx x22 L x2 h úú êa a22 L a2 n úú
B = ê 21 22 X = ê 21 A = ê 21
êM M O M ú ê M M O M ú ê M M O M ú
ê ú ê ú ê ú
ëbn1 bn 2 L bnh û n´ h ë xn1 xn 2 L xnh û n´ h ë an1 an 2 L ann û n´ n

Suppose further that they satisfy the following relation:

B=A×X
If A and B are known, we need a fast method to solve for X. One may suggest the following: compute the
inverse matrix A-1 of A and the solution is simply X = A-1B. While this is a correct way to solve the
problem, it is a little overkill. One might also observe the following fact: column j of matrix B is the
product of matrix A and column j of matrix X as shown below:

é b1 j ù é a11 a12 L a1n ù é x1 j ù

êb ú ê êx ú
ê 2 j ú = ê a21 a22 L a2 n úú ê 2jú
ê M ú ê M M O M ú ê M ú
ê ú ê ú ê ú
êë bnj úû ë an1 an 2 L ann û êë xnj úû

With this in mind, we can solve column j of X using A and column j of B. In this way, the given problem is
equivalent to solving h systems of linear equations. Actually, this is not a good strategy either because in
solving for column 1 of X matrix A will be destroyed, and, as a result, for each column we need to make a
copy of matrix A before running the linear system solver.
So, we need a method that does not use matrix inversion and does not have to copy matrix A over and over.
A possible way is the use of the LU decomposition technique.

LU Decomposition

An efficient procedure for solving B = A.X is the LU-decomposition. While other methods such as
Gaussian elimination method and Cholesky method can do the job well, this LU-decomposition method can
help accelerate the computation.
The LU-decomposition method first "decomposes" matrix A into A = L.U, where L and U are lower
triangular and upper triangular matrices, respectively. More precisely, if A is a n×n matrix, L and U are also
n×n matrices with forms like the following:

é l11 0 L 0ù éu11 u12 L u1n ù

êl l22 L 0 úú ê0 u22 L u2 n úú
L = ê 21 U=ê
êM M O Mú êM M O M ú
ê ú ê ú
ëln1 ln 2 L lnn û ë0 0 L unn û

15
The lower triangular matrix L has zeros in all entries above its diagonal and the upper triangular matrix U
has zeros in all entries below its diagonal. If the LU-decomposition of A = L.U is found, the original
equation becomes B = (L.U).X. This equation can be rewritten as B = L.(U.X). Since L and B are known,
solving for B = L.Y gives Y = U.X. Then, since U and Y are known, solving for X from Y = U.X yields the
desired result. In this way, the original problem of solving for X from B = A.X is decomposed into two
steps:
1. Solving for Y from B = L.Y
2. Solving for X from Y = U.X

Forward Substitution

How easy are these two steps? It turns out to be very easy. Consider the first step. Expanding B = L.Y gives

é b11 b12 L b1h ù é l11 0 L 0 ù é y11 y12 L y1h ù

êb L b2 h úú êêl21 l22 L 0 úú êê y21 L y2 h úú
ê 21 b22 = ·
y22
êM M O M ú êM M O Mú ê M M O M ú
ê ú ê ú ê ú
ëbn1 bn 2 L bnh û ëln1 ln 2 L lnn û ë yn1 yn 2 L ynh û

It is not difficult to verify that column j of matrix B is the product of matrix A and column j of matrix Y.
Therefore, we can solve one column of Y at a time. This is shown below:

é b1 ù é l11 0 L 0 ù é y1 ù
êb ú êl L 0 úú êê y2 úú
ê 2 ú = ê 21 l22 ·
êMú êM M O Mú êMú
ê ú ê ú ê ú
ëbn û ëln1 ln 2 L lnn û ë yn û

This equation is equivalent to the following:

b1 = l11 y1
b2 = l21 y1 + l22 y2
M M O
bn = ln1 y1 + ln 2 y2 + ln3 y3 + L + lnn yn

From the above equations, we see that y1 = b1/l11. Once we have y1 available, the second equation yields y2
= (b2-l21y1)/l22. Now we have y1 and y2, from equation 3, we have y3 = (b3 - (l31y1 +l32y2)/l33. Thus, we
compute y1 from the first equation and substitute it into the second to compute y2. Once y1 and y2 are
available, they are substituted into the third equation to solve for y3. Repeating this process, when we reach
equation i, we will have y1, y2, ..., yi-1 available. Then, they are substituted into equation i to solve for yi
using the following formula:

1 é i -1
ù
yi =
lii êbi - å li ,k yk ú
ë k =1 û
Because the values of the yi's are substituted to solve for the next value of y, this process is referred to as
forward substitution. We can repeat the forward substitution process for each column of Y and its
corresponding column of B. The result is the solution Y. The following is an algorithm:
Input: Matrix Bn×h and a lower triangular matrix Ln×h
Output: Matrix Yn×h such that B = L.Y holds.
Algorithm:
/* there are h columns */

16
 for j := 1 to h do
/* do the following for each column */
begin
/* compute y1 of the current column */
y1,j = b1,j / l1,1;
for i := 2 to n do /* process elements on that column */
begin
sum := 0; /* solving for yi of the current column */
for k := 1 to i-1 do
sum := sum + li,k× yk,j;
yi,j = (bi,j - sum)/li,i;
end
end

Backward Substitution

After Y becomes available, we can solve for X from Y = U.X. Expanding this equation and only
considering a particular column of Y and the corresponding column of X yields the following:

é y1 ù éu11 u12 L u1n ù é x1 ù

êy ú ê 0 u L u2 n úú êê x2 úú
ê 2ú = ê 22
·
êMú êM M O M ú êMú
ê ú ê ú ê ú
ë yn û ë 0 0 L unn û ë xn û

This equation is equivalent to the following:

y1 = u11 x1 + u12 x2 + u13 x3 + L + u1n xn
y2 = L u22 x2 + u23 x3 + L + u2 n xn
M O
yn = unn xn

Now, xn is immediately available from equation n, because xn = yn/un,n. Once xn is available, plugging it into
equation n-1

yn-1 = un -1,n -1 xn -1 + un -1,n xn

and solving for xn-1 yields xn-1 = (yn-1- un-1,nxn)/ un-1,n-1. Now, we have xn and xn-1. Plugging them into
equation n-2

yn -2 = un - 2,n- 2 xn- 2 + un -2,n -1 xn -1 + un - 2,n xn

and solving for xn-2 yields xn-2 = [yn-2- (un-2,n-1xn-1 + un-2,nxn-)]/ un-2,n-2.
From xn, xn-1 and xn-2, we can solve for xn-3 from equation n-3. In general, after xn, xn-1, ..., xi+1 become
available, we can solve for xi from equation i using the following relation:

1 é n
ù
xi =
uii ê
ë
y i - å
k = i +1
ui ,k xk ú
û

17
Repeat this process until x1 is computed. Then, all unknown x's are available and the system of linear
equations is solved. The following algorithm summarizes this process:
Input: Matrix Yn×h and a upper triangular matrix Un×h
Output: Matrix Xn×h such that Y = U.X holds.
Algorithm:
/* there are h columns */
for j := 1 to h do
/* do the following for each column */
begin
/* compute xn of the current column */
xn,j = yn,j / un,n;
for i := n-1 downto 1 do /* process elements of that column */
begin
sum := 0; /* solving for xi on the current column */
for k := i+1 to n do
sum := sum + ui,k× xk,j;
xi,j = (yi,j - sum)/ui,i;
end
end
This time we work backward, from xn backward to x1, and, hence, this process is referred to as backward
substitution.
LU-decomposition and forward and backward substitutions, including their subroutines/functions, should
be available in many numerical method textbooks and mathematical libraries. Check your numerical
methods textbook for the details.

18
Curve Global Interpolation

The simplest method of fitting a set of data points with a B-spline curve is the global interpolation method.
The meaning of global will be clear later on this page.

Suppose we have n+1 data points D0, D1, ..., Dn and wish to fit them with a B-spline curve of degree p,
where p £ n is an input. With the technique discussed in Parameter Selection and Knot Vector
Generation, we can select a set of parameter values t0, t1, ..., tn. Note that the number of parameters is equal
to the number of data points, because parameter tk corresponds to data point Dk. From these parameters, a
knot vector of m + 1 knots is computed where m = n + p + 1. So, we have a knot vector and the degree p,
the only missing part is a set of n+1 control points. Global interpolation is a simple way of finding these
control points.

Global Curve Interpolation

Given a set of n+1 data points, D0, D2, ..., Dn and a degree p, find a B-spline
curve of degree p defined by n+1 control points that passes all data points in
the given order.

Finding a Solution

Suppose the desired interpolating B-spline curve of degree p is the following:

n
C(u ) = å N i , p (u )Pi
i=0

This B-spline curve has n + 1 unknown control points. Since parameter tk corresponds to data point Dk,
plugging tk into the above equation yields the following:

n
Dk = C(uk ) = å N i , p (uk )Pi for 0 £ k £ n
i=0

Because there are n + 1 B-spline basis functions in the above equation (N0,p(u), N1,p(u), N2,p(u), ..., and
Nn,p(u) and n +1 parameters (t0, t1, t2, .., and tn), plugging these tk's into the Ni,p(u)'s yields (n+1)2 values.
Click here to learn how to compute these coefficients. These values can be organized into a (n+1)×(n+1)
matrix N in which the k-th row contains the values of N0,p(u), N1,p(u), N2,p(u), ..., and Ni,p(u) evaluated at tk
as shown below:

19
 é N 0, p (t0 ) N1, p (t0 ) N 2, p (t0 ) L N n , p (t0 ) ù
ê N (t ) N1, p (t1 ) N 2, p (t1 ) L N n , p (t1 ) úú
N=ê
0, p 1

ê M M M O M ú
ê ú
ëê N 0, p (tn ) N1, p (tn ) N 2, p (tn ) L N n , p (tn ) ûú

Let us also collect vectors Dk and Pi into two matrices D and P as follows:

é d 01 d 02 d 03 L d 0 s ù é p01 p02 p03 L p0 s ù

êd d12 d13 L d1s úú êp p12 p13 L p1s úú
D = ê 11 P = ê 11
ê M M M O M ú ê M M M O M ú
ê ú ê ú
ë d n1 dn 2 d n 3 L d ns û ë pn1 pn 2 pn 3 L pns û

Here, we treat data point Dk as a vector in s-dimensional space (i.e., Dk = [ dk1, ....., dks]) and appears on the
k-th row of matrix D. Similarly, we treat Pi as a vector in s-dimensional space (i.e., Pi = [ pi1, ....., pis]) and
appears on the i-th row of matrix P. Note that in three-dimensional space, we have s = 3, and in a plane we
have s = 2. Note also that D and P are both (n+1)×s matrices. Then, it is not difficult to verify that the
relation of dk's and ti's is transformed into the following simpler form:

D=N×P

Since matrix D contains the input data points and matrix N are obtained by evaluating B-spline basis
functions at the given parameters, D and N both are known and the only unknown is matrix P. Since the
above is simply a system of linear equations with unknown P, solving for P yields the control points and the
desired B-spline interpolation curve becomes available. Therefore, the interpolation problem is solved!

An Algorithm

You might have something in doubt because matrices D and P are not column matrices as discussed in your
numerical methods and linear algebra textbooks. This actually causes no problem at all. We can solve this
system column by column. More precisely, let the i-th column of D and the i-the column of P be di and pi,
respectively. From the above system of linear equations, we have the following:
di=N×Pi
Then, solving for pi from N and di gives the i-th column of P. Do this for every i in the range of 0 and h,
and we will have a complete P. Hence, all control points are computed. As you might have noticed, this is
very inefficient. Fortunately, many numerical libraries provide us with a linear system solver that is capable
of solving the equation D = N.P efficiently. Thus, fitting a B-spline curve to a set of n+1 data points is not
very difficult. It boils down to the solution of a system of linear equations. The following algorithm
summarizes the required steps:
Input: n+1 data points D0, ..... Dn and a degree p
Output: A B-spline curve of degree p that contains all data points in the given order
Algorithm:
Select a method for computing a set of n+1 parameters t0, ..., tn;
As a by-product, we also receive a knot vector U;

20
 for i := 0 to n do
for j := 0 to n do
Evaluate Nj,p(ti) into row i and column j of matrix N;
/* matrix N is available */
for i := 0 to n do
Place data point Di on row i of matrix D;
/* matrix D is constructed */
Use a linear system solver to solve for P from D = N.P
Row i of P is control point Pi;
Control points P0, ..., Pn, knot vector U, and degree p determine an interpolating B-spline
curve;
Figure (a) below is an example. There are 9 data points in black. The computed control points are in blue.
The small blue dots on the interpolating curve are points corresponding to the knots which are computed
using the chord length method. In this case, these "knot points" are quite close to the data points and the
control polygon also follows the data polygon closely, although it is not always the case. In Figure (b)
below, the data polygon and the control polygon are very different.

(a) (b)

It is important to point out that matrix N is totally positive and banded with a semi-bandwidth less than p
(i.e., Ni,p(tk) = 0 if |i - k| ³ p) if the knots are computed by averaging consecutive p parameters. This was
proved by de Boor in 1978. See Knot Vector Generation for the details. This means the system of linear
equations obtained from the above algorithm can be solved using Gaussian elimination without pivoting.

The Impact of Parameters and Knots

In general, the impact of the selected parameters and knots cannot be predicted easily. However, we can
safely say that if the chord length distribution is about the same, all four parameter selection methods
should perform similarly. Moreover, the universal method should perform similar to the uniform method
because the maximums of B-spline basis functions with uniform knots are distributed quite uniformly.
Similarly, the centripetal method should work similar to the chord length method because the former is an

21
extension to the latter. However, it is not always the case when the distribution of chord lengths change
wildly.
The following four curves are obtained using four parameter selection methods and the degree of the
interpolating B-spline curve is 3. The uniform method generates a cusp, the chord length method forces the
curve wiggle through data points wildly, the centripetal method resemblances the chord length method but
performs better, and the universal method follows the data polygon closely (better than the uniform method)
but produces a small loop.

uniform chord length

centripetal universal

How about the relation between parameters and knots? The following diagram shows the parameters and
knots of all four methods. As you can see, the parameters and knots obtained from the universal method is
more evenly distributed than those of the chord length method and the centripetal method. Moreover, the
parameters and knots obtained from the centripetal method stretch the shorter (resp, longer) chords longer
(resp., shorter) and hence are more evenly distributed. Therefore, the longer curve segments obtained by the

22
 chord length method become shorter in the centripetal method and the curve does not wiggle wildly through
data points.

The Impact of Degree

The impact of degree to the shape of the interpolating B-spline curve is also difficult to predict. One can
easily observe, from the following images, that the uniformly spaced method and universal method usually
follow long chords very well. On the other hand, these two methods have problems with short chords.
Because the parameters are equally or almost equally spaced, the interpolating curves have to stretch a little
longer for shorter chords. As a result, we see peaks and loops. This situation gets worse with higher degree
curves because higher degree curves provide more freedom to wiggle.

Uniform Chord Length Centripetal Universal

Deg 2

Deg 3

23
Deg 4

Deg 5

As for the chord length method, the above figures show that it does not work very well for longer chords,
especially those followed or preceded by a number of shorter chords, for which big bulges may occur.
There is no significant impact of degree on the shape of the interpolating curves shown above. Since the
centripetal method is an extension to the chord length method, they share the same characteristics. However,
since the centripetal method has a tendency to even out the distance between two adjacent parameters, it
also share the same characteristics of the uniform and universal methods. For example, the generated
interpolating curves follow longer chords closely and loops may occur for shorter chords when degree
increases. In fact, results above show that the centripetal method and the universal method perform
similarly.

Why Is This Method Global?

This interpolation method is global even with the use of B-spline curves which satisfy the local
modification property, because changing the position of a
single data point changes the shape of the interpolating
curve completely. In the following figure, those yellow
dots are data points and one of them is moved to its new
position, marked in light blue and indicated with a red
arrow. These nine data points are fit with a B-spline curve
of degree 4 using the centripetal method. As you can see,
the resulting curve (in blue) and the original curve have
eight data points identical, and the eight curve segments
are all different. Therefore, changing the position of s
single data point changes the shape of the interpolating
curve globally!

24
Curve Global Approximation

In interpolation, the interpolating curve passes through all given data points in the given order. As discussed
on the Global Interpolation page, an interpolating curve may wiggle through all data points rather than
following the data polygon closely. The approximation technique is introduced to overcome this problem
by relaxing the strict requirement that the curve must contain all data points. Except for the first and last
data points, the curve does not have to contain any other point. To measure how well a curve can
"approximate" the given data polygon, the concept of error distance is used. The error distance is the
distance between a data point and its "corresponding" point on the curve. Thus, if the sum of these error
distances is minimized, the curve should follow the shape of the data polygon closely. An interpolating
curve is certainly such a solution because the error distance of each data point is zero. However, the
formulation to be discussed is unlikely to get an interpolating curve. A curve obtained this way is referred
to as an approximation curve.
Suppose we are given n+1 data points D0, D2, ..., Dn, and wish to find a B-spline curve that can follow the
shape of the data polygon without actually containing the data points. To do so, we need two more input:
the number of control points (i.e., h+1) and a degree (p), where n > h ³ p³ 1 must hold. Thus,
approximation is more flexible than interpolation, because we not only select a degree but also the number
of control points. The following is a summary of our problem:

Global Curve Approximation

Given a set of n+1 data points, D0, D2, ..., Dn, a degree p, and a number h, where
n > h ³ p ³ 1, find a B-spline curve of degree p defined by h+1 control points that
satisfies the following conditions:
1. this curve contains the first and last data points (i.e., D0 and Dn) and
2. this curve approximates the data polygon in the sense of least square.

With h and p in hand, we can determine a set of parameters and a knot vector. Let the parameters be t0, t1, ...,
tn. Note that the number of parameters is equal to the number of data points. Now, suppose the
approximation B-spline of degree p is
h
C(u ) = å N i , p (u )Pi
i=0

where P0, P1, ..., Ph are the h+1 unknown control points. Since we want the curve to pass the first and last
data points, we have D0 = C(0) = P0 and Dn = C(1) = Ph. Therefore, there are only h - 1 unknown control
points P1, P2, ..., Ph-1. Taking this into consideration, the curve equation becomes the following:
h
C(u ) = å N i , p (u )Pi
i=0

h -1
C(u ) = N 0, p (u )P0 + å N i , p (u )Pi + N h , p (u )Ph
i =1

The Meaning of Least Square

25
How do we measure the error distance? Because parameter tk corresponds to data point Dk, the distance
between Dk and the corresponding point of tk on the curve is |Dk - C(tk)|. Since this distance consists of a
square root which is not easy to handle, we choose to use the squared distance |Dk - C(tk)|2. Hence, the sum
of all squared error distances is
h -1
f (P1 ,..., Ph -1 ) = å Dk - C (tk )
2

k =1

Our goal is, of course, to find those control points P1, ..., Ph-1 such that the function f() is minimized. Thus,
the approximation is done in the sense of least square.

Finding a Solution

This section is going to be a little messy and requires some knowledge in linear algebra. First, let us rewrite
Dk - C(tk) into a different form:

é h -1
ù
Dk - C (tk ) = Dk - ê N 0, p (u )P0 + å N i , p (u )Pi + N h , p (u )Ph ú
ë i =1 û
h -1
= ( Dk - N 0, p (u )P0 - N h , p (u )Ph ) - å N i , p (u )Pi
i =1

In the above, D0, Dk and Dn are given, and N0,p(tk) and Nh,p(tk) can be obtained by evaluating N0,p(u) and
Nh,p(u) at tk. Click here to learn how to compute these coefficients. For convenience, let us define a new
vector Qk as:

Q k = ( Dk - N 0, p (u )P0 - N h, p (u )Ph )

Then, the sum-of-square function f() can be written as follows:

h -1 2
f (P1 ,..., Ph -1 ) = å Q k - å i =1 N i , p (tk )Pi
h -1

k =1

Next, we shall find out what the squared error distance looks like. Recall the identity x.x = |x|2. This means
the inner product of vector x with itself gives the squared length of x. Thus, the error square term can be
rewritten as

h -1 2

Q k - å N i , p (tk )Pi
i =1

æ h -1
öæ h -1
ö
= ç Q k - å N i , p (tk )Pi ÷ ç Q k - å N i , p (tk )Pi ÷
è i =1 øè i =1 ø
æ h -1 ö æ h -1 ö æ h -1 ö
= Q k Q k - 2Q k ç å N i , p (tk )Pi ÷ + ç å N i , p (tk )Pi ÷ ç å N i , p (tk )Pi ÷
è i =1 ø è i =1 ø è i =1 ø

Then, function f() becomes

h -1
é æ h -1 ö æ h -1 ö æ h -1 öù
f (P1 ,..., Ph -1 ) = å êQ k Q k - 2Q k ç å N i , p (tk )Pi ÷ + ç å N i , p (tk )Pi ÷ ç å N i , p (tk )Pi ÷ ú
k =1 ë è i =1 ø è i =1 ø è i =1 øû

How do we minimize this function? Function f() is actually an elliptic paraboloid in variables P1, ..., Ph-1.
Therefore, we can differentiate f() with respect to each pg and find the common zeros of these partial
derivatives. These zeros are the values at which function f() attends its minimum.

26
In computing the derivative with respect to Pg, note that all Qk's and Ni,p(tk) are constants (i.e., no Pk
involved) and their partial derivatives with respect any Pg must be zero. Therefore, we have
¶
( Qk × Qk ) = 0
¶Pg

Consider the second term in the summation, which is the sum of Ni,p(tk)Pi .Qk's. The derivative of each
sub-term is computed as follows:
¶ ¶P
¶Pg
( N i , p (tk )Pi × Q k ) = N i , p (tk ) i × Q k
¶Pg

The partial derivative of Pi with respect to Pg is non-zero only if i = g. Therefore, the partial derivative of
the second term is the following:

¶ æ h -1 ö
å Ni, p (tk )Pi × Q k ÷ø = N g , p (tk ) × Qk
¶Pg çè i =1

The derivative of the third term is more complicated; but it is still simple. The following uses the
multiplication rule (f.g)' = f'.g + f.g'.

¶ æ h -1 öæ h -1 ö æ h -1 ö æ h -1 ¶Pi ö
ç å N i , p (tk )Pi ÷ç å N i , p (tk )Pi ÷ = 2 ç å N i , p (tk )Pi ÷ çç å N i , p (tk ) ÷÷
¶Pg è i =1 øè i =1 ø è i =1 ø è i =1 ¶Pg ø

Since the partial derivative of Pi with respect to Pg is zero if i is not equal to g, the partial derivative of the
third term in the summation with respect to pg is:

¶ æ h -1 öæ h -1 ö æ h -1 ö
ç å
¶Pg è i =1
N i , p (tk )Pi ÷ç å N i , p (tk )Pi ÷ = 2 N g , p (tk ) ç å N i , p (tk )Pi ÷
øè i =1 ø è i =1 ø

Combining these results, the partial derivative of f() with respect to pg is

¶f æ h -1 ö
= -2 N g , p (tk )Q k + 2 N g , p (tk ) ç å N i , p (tk )Pi ÷
¶Pg è i =1 ø

Setting it to zero, we have the following:

n -1 h -1 n -1

åN
k =1
g, p (tk )å N i , p (tk )Pi = å N g , p (tk )Q k
i =1 k =1

Since we have h-1 variables, g runs from 1 to h-1 and there are h-1 such equations. Note that these
equations are linear in the unknowns Pi's. What is the coefficient of Pi? Before we go on, let us define three
matrices:

é n -1 ù
ê å N1, p (tk )Q k ú
é P1 ù ê k =1 ú é N1, p (t1 ) N 2, p (t1 ) L N h -1, p (t1 ) ù
ê n -1 ú ê N (t )
êP ú
ê å N 2, p (tk )Q k ú N 2, p (t2 ) L N h -1, p (t2 ) úú
P=ê 2 ú N=ê
1, p 2
Q = ê k =1 ú
ê M ú ê M M O M ú
ê ú ê M ú ê ú
ë Ph-1 û ê n -1 ú ëê 1, p (tn -1 )
N N 2, p (tn -1 ) L N h -1, p (tn -1 ) ûú
ê N
å h -1, p k k úûú
ëê k =1
(t )Q

27
Here, the k-th row of P is the vector Pk, the k-th row of Q is the right-hand side of the k-th equation above,
and the k-th row of N contains the values of evaluating N1,p(u), N2,p(u), ..., Nh-1,p(u) at tk. Therefore, if the
input data points are s-dimensional vectors, P, N and Q are (h-1)×s, (n-1)×(h-1) and (h-1)×s matrices,
respectively.
Now, let us rewrite the g-th linear equation
n -1 h -1 n -1

åN
k =1
g, p (tk )å N i , p (tk )Pi = å N g , p (tk )Q k
i =1 k =1

into a different form so that the coefficient of pi can easily be read off:
h -1
æ n -1 ö n -1

å çè å N
i =1 k =1
(t
g, p k ) N (t )
i, p k ÷ i
ø
P = å
k =1
N g , p (tk )Q k

Therefore, the coefficient of Pi is

n -1

åN
k =1
g, p (tk ) N i , p (tk )

If you look at matrix N, you should see that Ng,p(t1), Ng,p(t2), ..., Ng,p(tn-1) are the g-th column of N, and
Ni,p(t1), Ni,p(t2), ..., Ni,p(tn-1) are the i-th column of N. Note that the g-th column of N is the g-th row of N's
transpose matrix NT, and the coefficient of Pi is the "inner" product of the g-th row of NT and the i-th
column of N. With this observation, the system of linear equations can be rewritten as

(N N) P = Q
T

Since N and Q are known, solving this system of linear equations for P gives us the desired control points.
And, of course, we have found the solution.

An Algorithm

Although the development in the previous section is lengthy and tedious, the final result is surprisingly
simple: solving a system of linear equations just like we encountered in global interpolation. As a result, the
algorithm for global approximation is also quite simple.
Input: n+1 data points D0, D1, ..., Dn, a degree p, and the desired number of control points
h+1;
Output: A B-spline curve of degree p, defined by h+1 control points, that approximates the
given data points;
Algorithm:
Obtain a set of parameters t0, ..., tn and a knot vector U;
Let P0 = D0 and Ph = Dn;
for k := 1 to n-1 do
Compute Qk using the following formula:

Q k = ( Dk - N 0, p (u )P0 - N h, p (u )Ph )

for i := 1 to h-1 do
Compute the following and save it to the i-th row of matrix Q;
h -1

åN
i =1
i, p (tk )Q k

28
 /* matrix Q is available */
for k := 1 to n-1 do
for i := 1 to h-1 do
Compute Ni,p(tk) and save to row k and column i of N;
/* matrix N is available */
Compute M = NT.N;
Solving for P from M.P = Q;
Row i of P is control point Pi;
Control points P0, ..., Ph, knot vector U and degree p determines an approximation
B-spline curve;

The Impact of Degree and Number of Control Points

It is obvious that the data points affect the shape of the approximation curve. What is the impact of degree p
and the number of control points on the shape of the curve? The following figures show some
approximation curves for 10 data points (n=9) with various degrees and numbers of control points. Each
row gives the approximation curves of the same degree but with different number of control points, while
each column shows the same number of control points with different degree. Note that all curves use the
centripetal method for computing its parameters.

Points = 4 Points = 5 Points = 6 Points = 6

Deg 2

Deg 3

Deg 4

29
Deg 5

It is understandable that lower degree in general will not be able to approximate the data polygon well
because lower degree curves lack of flexibility. Therefore, on each column, a higher the degree yields a
better result (i.e., closer to the data polygon). By the same reason, more control points offer higher
flexibility of the approximation curve. Hence, on each row, as the number of control points increases, the
curve becomes closer to the data polygon.
Shall we use higher degree and many control points? The answer is no because the purpose of using
approximation is that we can use fewer number of control points than global interpolation. If the number of
control points is equal to the number of data points, we can just use global interpolation! As for degree, we
certainly want the smallest one as long as the produced curve can capture the shape of the data polygon.

Why Is This Method Global?

This approximation method is global because changing the position of data point causes the entire curve to
change. The yellow dots in the following image are the given data points to be approximated by a B-spline
curve of degree 3 and 5 control points (i.e., n = 7, p=3 and h
= 4). The parameters are computed with the centripetal
method. Suppose data point 4 is moved to a new position
indicated by a light blue dot. The new approximation
B-spline is the one in blue color. As you can see, except for
the first and last data points which are fixed by the
algorithm, the shapes of the original curve and the new one
are very different. Therefore, changes due to modifying
data points are global!

30
Surface Global Interpolation

Suppose we have m+1 rows and n+1 columns of data points Dij (0 £ i £ m and 0£ j£ n) and wish to find a
B-spline surface of degree (p,q) that contains all of them. Similar to the curve case, we have data points and
degrees p and q as input. To define an interpolating a B-spline surface, we need two knot vectors U and V,
one for each direction, and a set of control points. Like the curve case, the number of control points and the
number of data points are equal (i.e., (m+1)×(n+1) control points in m+1 rows and n+1 columns). With the
technique discussed in Parameters and Knot Vectors for Surfaces, we can compute two sets of
parameters sc (0 £ c £ m) in the u-direction and td (0 £ d £ n) in the v-direction by setting e and f to m and n,
respectively. We also obtain knot vectors U and V for the u- and v-direction, respectively. Therefore, what
remains to do is to find the desired control points!

Global Surface Interpolation

Given a grid of (m+1)×(n+1) data points Dij (0 £ i £ m and 0 £ j £ n) and a

degree (p, q), find a B-spline surface of degree (p,q) defined by (m+1)×(n+1)
control points that passes all data points in the given order.

Finding a Solution

Suppose the B-spline surface is given as follows:

m n
S(u , v) = åå Pij N i , p (u ) N j , q (v)
i =0 j = 0

Since it contains all data points and since parameters sc and td correspond to data point Dcd, plugging u=sc
and v=td into the surface equation yields:
m n
Dcd = S( sc , td ) = åå Pij N i , p ( sc ) N j , q (td )
i =0 j = 0

Since Ni,p(sc) is independent of the index j, we can move this term out of the summation of j:

m n
Dcd = S( sc , td ) = åå Pij N i , p ( sc ) N j , q (td )
i =0 j = 0

m æ n ö
= å N i , p ( sc ) ç å Pij N j , q (td ) ÷
i=0 è j =0 ø

If we examine the inner term, we shall see that the index i only appears in Pij. Therefore, we can define the
inner expression to be a new term as follows:
n
Qid = å Pij N j , q (td )
j =0

More precisely, if i is fixed to the same value, Qid is the point, evaluated at td, on the B-spline curve of
degree q defined by n+1 "unknown" control points on row i of the P's (i.e., Pi0, Pi1, ..., Pin). Plugging Qid
into the equation of Dcd above yields

31
 m
Dcd = å N i , p ( sc )Qid
i=0

Thus, data point Dcd is the point, evaluated at sc, of a B-spline curve of degree p defined by m+1 "unknown"
control points on column d of the Q's (i.e., Q0d, Q1d, ..., Qmd). Repeating this for every c (0 £ c £ m), the
d-th column of data points (i.e., D0d, D1d, ..., Dmd) is obtained from the d-th column of Q's (i.e., Q0d, Q1d, ...,
Qmd) and parameters s0, s1, ..., sm. Since data points D0d, D1d, ..., Dmd, degree p and parameters s0, s1, ..., sm
are known, this equation means

Find the d-th column of control points Q0d, Q1d, ..., Qmd, given degree p,
parameters s0, s1, ..., sm and the d-th column of data points D0d, D1d, ..., Dmd.

So, it is a curve interpolation problem! More precisely, the curve global interpolation can be applied to
each column of the data points to obtain the columns of control points Qcd's. Since we have n+1 columns of
data points, we shall obtain n+1 columns of Q's.
Now, we can use the same idea but apply to the equation of Qid, which is duplicated below. In this equation,
data points on row i of the Q's (i.e., Qi0, Qi1, ..., Qin) are the points on a B-spline curve, evaluated at t0, t1, ...,
tn, of degree q defined by n+1 unknown control points Pi0, Pi1, ..., Pin. Therefore, applying curve
interpolation with degree q and parameters t0, t1, ..., tn to row i of the Q's gives row i of the desired control
points!
m
Qid = å N j , q (td )Pij
i =0

Once all rows of control points are found, these control points, along with the two knot vectors and degrees
p and q define an interpolating B-spline surface of degree (p,q) of the given data points. Therefore, surface
interpolation using B-spline consists of (m+1) + (n+1) curve interpolations! The following summarizes the
required steps:

Input: (m+1)×(n+1) data points Dij and degree (p,q);

Output: A B-spline surface of degree (p,q) that contains all data points;
Algorithm:
Compute parameters in the u-direction s0, s1, ..., sm and the knot vector U;
Compute parameters in the v-direction t0, t1, ..., tn and the knot vector V;
for d := 0 to n do /* for column d */
begin /* computing "intermediate data points" Q's */
Apply curve interpolation to column d of data points (i.e., D0d, D1d, ... Dmd) using

§ degree p
§ parameters s0, s1, ..., sm
§ knot vector U

The result is column d of the "intermediate data points" Q0d, Q1d, ..., Qmd
end
for c := 0 to m do /* for row c */
begin /* computing the desired control points P's */
Apply curve interpolation to row c of the Q's (i.e., Qc0, Qc1, ... Qcn) using

32
 § degree q
§ parameters t0, t1, ..., tn
§ knot vector V

The result is row c of the desired control points Pc0, Pc1, ..., Pcn
end
The computed (m+1)×(n+1) control points Pij, degree (p,q), and knot vectors U and V define an
interpolating B-spline surface of degree (p,q) for the given data points.

Note that matrix N used in the first for loop does not change when going from column to column. Therefore,
if we naively follow the above algorithm, we would end up solving the system of equation D = N.Q n+1
times. This is, of course, not efficient. To speed up, the LU-decomposition of N should be computed before
the first for loop starts, and the interpolation in each iteration is simply a forward substitution followed by a
backward substitution. Similarly, a new matrix N will be used for the second for loop. Its
LU-decomposition should also be computed before the second for starts. Without doing so, we will perform
(m+1)+(n+1)=m+n+2 LU-decompositions when only 2 would be sufficient.

Why Is This Method Global?

Because the interpolating surface is obtained through m+n+2 global curve interpolations, it is obvious that
this interpolation technique is global. The following images show the impact of moving a data point. This
B-spline surface is obtained by using six rows (m=5) and five columns (n=4) data points with degree (3,3).
Images on the top row show the knot curves (i.e., isoparametric curves that correspond to knots), and the
data point marked by a yellow circle indicates the point to be moved. The bottom row shows the
corresponding rendered surfaces. It is obvious that after changing the position of a data point, the blue knot
curve changes its shape drastically, and moves closer to its neighboring data points. The impact propagates
to the right as well because the valley in the far right end is a little deeper as can be seen from the red knot
curves. The more dramatic change is shown in the rendered surfaces. Although the indicated data point is
moved to the right, this move creates a big bulge in the left end. Please also note that the front boundary
curve, actually a red knot curve, also changes its shape slightly. Therefore, the impact of changing the
position of a single data point can propagate to entire surface, and, hence, this is a global method.

33
 Before Move After Move

Knot
Curves

Surfaces

34
Surface Global Approximation

Suppose we wish to find a B-spline surface to approximate the (m+1)×(n+1) data points. Because an
approximation surface does not contain all given data points, like in the curve approximation case, we have
controls over the degrees and the number of control points. Therefore, in addition to the data points, the
input to this problem includes the degrees p and q in the u- and v-direction, and the number of rows e+1 and
the number of columns f+1 of control points. As in the curve case, the input values must satisfy m > e ³ p ³
1 and n > f ³ q ³ 1 to have a solution.

Global Surface Approximation

Given a grid of (m+1)×(n+1) data points Dij (0 £ i £ m and 0 £ j £ n) a

degree (p, q), and e and f satisfying m > e ³ p ³ 1 and n > f ³ q ³ 1, find a
B-spline surface of degree (p,q) defined by (e+1)×(f+1) control points Pij
(0 £ i £ e and 0 £ j £ f) that approximates the data point grid in the given
order.

Finding a Solution

Let the desired B-spline surface of degree (p,q) defined by the unknown (e+1)×(f+1) control points Pij's is
the following:
e f
S(u , v) = åå Pij N i , p (u ) N j , q (v)
i =0 j = 0

Because there are m+1 rows of data points, we need m+1 parameters in the u-direction, s0, s1, ..., sm.
Similarly, we need n+1 parameters in the v-direction, t0, t1, ..., tn. The computation of these parameters is
discussed in Parameters and Knot Vectors for Surfaces. With these parameters in hand, the point on
surface that corresponds to data point Dcd is computed as follows:
e f
S( sc , td ) = åå Pij N i , p ( sc ) N j , q (td )
i =0 j = 0

The squared error distance between Dcd and its corresponding point on surface is

2
Dcd - S( sc , td )

The sum of all squared error distances is therefore the following:

m n
f (P00 , P01 ,..., Pef ) = åå Dcd - S( sc , td )
2

c = 0 d =0

This is a function in the (e+1)×(f+1) unknown control points Pij's. Like what we did for the global
approximation, to minimize f(), we can compute the partial derivatives and set them to zero:
¶f
=0
¶Pij

35
Then, we have (e+1)×(f+1) equations whose common zeros are the desired control points. Unfortunately,
these equations are not linear and solving system of non-linear equations is a very time consuming process.
Rather than aiming for an optimal solution, we can find a reasonably good solution that does not minimize
function f().
To find a non-optimal solution, we shall employ the same technique used in Global Surface Interpolation.
More precisely, apply curve approximation to each column of data points to compute some "intermediate"
data points. In this way, for each column of m+1 data points we compute e+1 "intermediate" data points.
Since there are n+1 columns, these "intermediate" data points form a (e+1)×(n+1) grid. Then, apply curve
approximation to each row of these intermediate data points to compute the desired control points. Since
each row has n+1 "intermediate" data point and there are e+1 rows, each approximation for a row generates
f+1 desired control points, and, as a result, we will finally obtain (e+1)×(f+1) control points. The following
summarizes the required steps:

Input: (m+1)×(n+1) data points Dij, degree (p,q), and e and f for
e+1 rows and f+1 columns of control points;
Output: A B-spline surface of degree (p,q) that approximates the data points;
Algorithm:
Compute parameters in the u-direction s0, s1, ..., sm and the knot vector U;
Compute parameters in the v-direction t0, t1, ..., tn and the knot vector V;
for d := 0 to n do /* for column d of D's */
begin /* computing "intermediate data points" Q's */
Apply curve approximation to column d of data points (i.e., D0d, D1d, ... Dmd)
using

§ degree p
§ parameters s0, s1, ..., sm
§ knot vector U

The result is column d of the "intermediate data points" Q0d, Q1d, ..., Qed
/* Q's form a (e+1)×(n+1) matrix */
end
for c := 0 to e do /* for row c of Q's */
begin /* computing the desired control points P's */
Apply curve approximation to row c of the Q's (i.e., Qc0, Qc1, ... Qcn) using

§ degree q
§ parameters t0, t1, ..., tn
§ knot vector V

The result is row c of the desired control points Pc0, Pc1, ..., Pcf
/* P's form a (e+1)×(f+1) matrix */
end
The computed (e+1)×(f+1) control points Pij, degree (p,q), and knot vectors U and V define an
approximation B-spline surface of degree (p,q) for the given data points.

36
In this algorithm, n+1 curve approximations are applied to the columns of D's and then e+1 curve
approximations are applied to the rows of the "intermediate" data points. Thus, n+e+2 curve approximations
are performed. On the other hand, we could apply m+1 curve approximations to each row of the data points
to create (m+1)×(f+1) "intermediate" data points. Then, f+1 curve approximations are applied to each
column of this "intermediate" data points to obtain the final (e+1)×(f+1) control points. In this way, m+f+2
curve approximations are performed.
Since this algorithm does not minimize function f(), it is not an optimal one, although it is adequate for
many applications.
Note that matrix NT.N used in the first for loop does not change when going from column to column.
Therefore, if we naively follow the above algorithm, we would end up solving the system of equation Q =
(NT.N).P n+1 times. This is, of course, not efficient. To speed up, the LU-decomposition of NT.N should be
computed before the first for loop starts, and the approximation in each iteration is simply a forward
substitution followed by a backward substitution. Similarly, a new matrix NT.N will be used for the second
for loop. Its LU-decomposition should also be computed before the second for starts. Without doing so, we
will perform (e+1)+(n+1)=e+n+2 LU-decompositions when only 2 would be sufficient.

A Simple Comparison

The following is a grid of six rows and five columns. We shall use it to illustrate the differences between
interpolation and approximation. Images on the first column are the results before changing the indicated
data points, while images on the second depict the results of this move.
The first two rows show the results of global interpolation using the chord length method for determining
parameters and knot vectors. The interpolating surface is of degree (3,3). The red and blue curves are knots
curves (i.e., isoparametric curves at knots). The indicated data point changes its position slightly, but the
blue knot curve changes its shape drastically. The difference between the two interpolating surfaces is even
greater. Moving a data point to the right produces a large bulge in the left end of the surface. Moreover, a
deep valley below the red knot curve is generated due to this move.
The last two rows show the results of global approximation. The number of control points in u- and
v-directions are 5 and 4, respectively. It is obvious that the impact of changing the position of a data point
on the shape of the approximation surface is not so dramatic. In fact, in this example, approximation is not
so sensitive to changing the position of the indicated data point, and the approximation surfaces follow the
shape of the data point net much better than those interpolating surfaces.

37
 Before Move After Move

Int:
Curve

Int:
Surface

App:
Curve

38
App:
Surface

39