# Cubic spline (Python)

##  cubic spline interpolator

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

<<cubicspline.py>>=
from numpy import *


We precalculate a set of cubic Bernstein bases, starting with a linear base. Instead of a continuous t, we'll step 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.
$\begin{matrix} B_3 & = & (1-t)^3 \\ B_2 & = & 3(1-t)^2t \\ B_1 & = & 3(1-t)t^2 \\ B_0 & = & t^3 \end{matrix}$

<<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 ttt + tt(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 B3(t)c0 + B2(t)c1 + B1(t)c2 + B0(t)c3

<<definition>>=
def cubicspline(c,t): return dot(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?

##  wrapping up

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

<<cubicspline.py>>=
lookup table
definition
if __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))

print "%!"
print "20 200 translate"
print ".8 .8 scale"
print "%d %d moveto" % tuple(cs[0][0])

for (x,y) in [cubicspline(c,16*t) for c in cs for t in arange(17)]:
print x,y,"lineto"

print "stroke showpage"


which, when run, yields the following output:

hijacker
hijacker
hijacker
hijacker