61 p[k]=0.5*y2[k]

62 h1=h2

63 f[n]=y[n]+h1/w[n]

64 for k in range(n):

65 j2=k+1

66 h=g[k]

67 g[k]=(f[j2]-f[k])*h-(y2[j2]+2.*y2[k])/(6.*h)

68 g[n]=(f[n]-f[n1])*h+(2.*y2[n]+y2[n1])/(6.*h)

69 curv=0. # compute curvature

70 for k in range(n1): curv=curv+q[k]*(f[k+1]-f[k])

71 curv=-6.*curv

72 print ™curvature = ™, curv, ™ ™,

73 chisq=0. # compute chisquare sum

74 for k in range(np1): chisq=chisq+((y[k]-f[k])/sig[k])**2

75 print ™chisquare = ™, chisq

76 return [[x,f,g],[curv,chisq]]

Comments

This program is adapted from Sp¨th (1973) and contains an explicit linear solver; it does not

a

require the package linalg. A trial weight must be given; the program returns both the curvature

and the chi-square sum. The weight should be adjusted such that the chi-square sum is statistically

acceptable, i.e., near the number of data points.

548 Splines for everything

19.7 B-splines

(p)

functions11 Ni (t)

B-spline are local smooth functions of a parameter t (p

is the order of the spline), which are used to construct smooth B-spline

curves in one- or multidimensional space by linear superposition, given a set

of control points (also called de Boor points) di at the nodes ti :

(p)

F (t) = Ni (t)di . (19.62)

i

Each B-spline function is a polynomial with local support, i.e., it has (posi-

tive) values only on a limited local interval and is zero outside this interval.

The support depends on the order p of the spline and is equal to the interval

[ti , ti+p ]. The higher the order p, the smoother the curve, but the larger its

support. If the control points are vectors in 3D space, a curve in 3D space re-

sults; the control points may also be real or complex scalars that reconstruct

a (real or complex) function f (x). B-spline surfaces result from a product

of spline functions with two parameters u and v and a two-dimensional set

of control points:

(p) (p)

F (u, v) = Ni (u)Nj (v)dij . (19.63)

i,j

The B-spline functions are normalized if they integrate to unity.

(p)

The normalized B-spline function Ni is de¬ned recursively:

(1)

p = 1: Ni (t) = 1 for ti ¤ t < ti+1 ; 0 elsewhere, (19.64)

ti+p ’ t

t ’ ti

(p) (p’1) (p’1)

p ≥ 2: Ni (t) = Ni (t) + Ni+1 (t).(19.65)

ti+p’1 ’ ti ti+p ’ ti+1

Most common are the cardinal B-splines, which have equidistant intervals

ti+1 ’ ti , and especially with distance 1 (“normalized nodes”). The nodes

(or knots) are then a set of consecutive integers. They produce uniform

B-spline curves. These cardinal B-spline functions are thus de¬ned as

(1)

p = 1: Ni (t) = 1 for i ¤ t < i + 1; 0 elsewhere, (19.66)

11 See Engeln-M¨ llges and Uhlig (1996). B-splines are also introduced in Essmann et al. (1995).

u

The name derives from B for “basis”. Although in principle already known by Laplace and

Euler, B-spline functions and their applications in approximation theory were developed in the

1960s and 1970s. The basic mathematical theory can be found in Schumaker (1981), but a more

practical guide is a published lecture series by Schoenberg (1973). A book by de Boor (1978)

provides Fortran programs. A relatively more modern description can be found in Chapter 4

of Chui (1992). There is no universally agreed notation for B-splines functions. Our notation

N is most common for the normalized cardinal spline functions, although the meaning of sub-

and superscripts may di¬er. Schoenberg (1973) denotes non-normalized B-spline functions by

Q and central B-spline functions, which are shifted to have the origin at the maximum of the

function, by M . The terms knots and nodes are both used for the control points. One should

distinguish between B-splines (curves or surfaces) and B-spline functions. The former are linear

combinations of the latter and could be called B-spline series.

19.7 B-splines 549

t ’ i (p’1) i + p ’ t (p’1)

(p)

p ≥ 2: Ni (t) = Ni (t) + Ni+1 (t). (19.67)

p’1 p’1

Clearly Ni is simply translated over i units with respect to N0 :

(p) (p) (p)

Ni (t) = N0 (t ’ i); support of Ni (t) : i ¤ t < i + p. (19.68)

We may omit the subscript 0 and thus de¬ne the basic B-spline functions

(p)

N (p) (t) for the support 0 ¤ t < p, from which all other functions Ni (t) are

derived by translation, as

p = 1: N (1) (t) = 1 for 0 ¤ t < 1; 0 elsewhere, (19.69)

p ’ t (p’1)

t

p ≥ 2: N (p) (t) = (t ’ 1).

N (p’1) (t) + N (19.70)

p’1 p’1

The functions are symmetric around t = p/2, where they attain a maximum.

Some texts de¬ne centralized B-spline functions by taking the origin at p/2;

for odd p this has the disadvantage that the knots do not occur at integer

arguments. The B-spline functions de¬ned here are also called forward B-

spline functions (Silliman, 1974).

Figure 19.9 shows the ¬rst ¬ve normalized cardinal splines for i = 0. The

functions N (p) (t) do not only integrate to 1, but the sum of their values at

the nodes equals unity as well:

i+p

(p)

Ni (j) = 1. (19.71)

j=i

Thus they, so to say, distribute the number 1 among points close to the

center p/2.

An interesting property of cardinal B-spline functions is that they repre-

sent the distribution generated by the sum of p independent random numbers

·, each homogeneously distributed in the interval 0 ¤ · < 1. As p becomes

larger, according to the central limit theorem, this distribution tends to a

normal distribution with mean p/2 and variance p/12. Therefore

(t ’ p/2)2

1 p

(t) = √ exp ’

(p)

, σ2 = .

lim N (19.72)

2σ 2 12

σ 2π

p’∞

This property is a direct consequence of a general convolution theorem for

B-spline functions: each N (p) is the convolution of N (p’q) (1 ¤ q < p) with

N (q) . For the special case q = 1 N (q) is the unit pulse function and the

convolution becomes another recurrent relation

1 t

—N (t ’ t ) dt =

(p) (p’1) (1) (p’1)

N (p’1) (t ) dt .

N (t) = (N )(t) = N

0 t’1

(19.73)

550 Splines for everything

1

(1)

p=1 N0

1