Cubic spline (Python)

From LiteratePrograms

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.

<<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 B₃(t)c₀ + B₂(t)c₁ + B₁(t)c₂ + B₀(t)c₃

<<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:

Download code