# Cubic spline (Python)

## [edit] cubic spline interpolator

Spline interpolation is repetitive math, not symbolic computation, so we will use the `Numeric Python`

package.

<<cubicspline.py>>=fromnumpyimport*

We precalculate a set of cubic Bernstein bases,
starting with a linear base.
Instead of a continuous *t*,
we'll **s**tep from 0 to 256 (inclusive!) by 1/256
to generate a discrete table useful over the range [0,1].

<<lookup table>>=s = arange(257)/256.0

We need 1 − *t* as well, but that is simple:
it is the mirror image of *t*.
Rather than recalulate, we take **z** to be that same as **s**, only read backwards. (with a stride of -1)

<<lookup table>>=z = s[::-1]

Now we build up each cubic basis as a function of our existing linear terms.

<<lookup table>>=b = transpose(array((z*z*z, 3*z*z*s, 3*z*s*s, s*s*s)))

Why the factors of three? When we perform a linear interpolation, we use a linear combination of *t* and 1 − *t* as weights. The square of that generates 4 weights, and the cube 8, but multiplication is commutative, so given *t**t**t* + *t**t*(1 − *t*) + *t*(1 − *t*)*t* + ... + (1 − *t*)(1 − *t*)*t* + (1 − *t*)(1 − *t*)(1 − *t*) we can cut the work in half by grouping like terms.

Why the transposition? So that cubic spline evaluation simplifies to the dot product *B*_{3}(*t*)*c*_{0} + *B*_{2}(*t*)*c*_{1} + *B*_{1}(*t*)*c*_{2} + *B*_{0}(*t*)*c*_{3}

<<definition>>=defcubicspline(c,t):returndot(b[t],c)

*Exercise: what changes are necessary to interpolate several sets of control points in parallel?*

*Question: which symmetry axes of the cube produce the 1-3-3-1 pattern?*

## [edit] wrapping up

Finally, we test the output by generating a little Postscript,

<<cubicspline.py>>=lookup table definitionif__name__ == "__main__": cs = reshape(array([ 533,179, 533,179, 483,0, 483,0, 483,0, 483,0, 28,0, 28,0, 28,0,28,0, 28,24, 28,24, 28,24, 28,24, 81,24, 81,24, 81,24, 109,24, 117,53, 117,76, 117,76, 117,76, 117,600, 117,600, 117,600, 117,618, 99,637, 81,637, 81,637, 81,637, 28,637, 28,637, 28,637, 28,637, 28,661, 28,661, 28,661, 28,661, 259,661, 259,661, 259,661, 259,661, 259,637, 259,637, 259,637, 259,637, 206,637, 206,637, 206,637, 187,637, 172,618, 172,600, 172,600, 172,600, 172,73, 172,73, 172,73, 172,50, 184,23, 210,23, 210,23, 210,23, 330,23, 330,23, 330,23, 416,23, 495,99, 518,179, 518,179, 518,179, 533,179, 533,179, ]),(-1,4,2))for(x,y)in[cubicspline(c,16*t)forcincsfortinarange(17)]:

which, when run, yields the following output:

Download code |