""" basic spline interpolation in a few lines of numpy.
Keywords: interpolate, resample, B-spline, Catmull-Rom, numpy, getting-started
Examples:
    t = arange( 1000+1 )
    points = f(t)  # points in 1d 2d 3d ...
    curve = dotspline1( points, times=10 )  # at 0 .1 .2 ... 999.9 1000

    surface = dotspline2( array, times=10 )  # an array of e.g. pixels

Aims: short, clear, understandable.

dotspline1 In:
    pt[ 0 1 2 ... ] scalars, or points in space, or vectors, or ...
    times, e.g. 10 => output points at at 0 .1 .2 ...
Out:
    "times" as many points, a spline curve through or near the input points.
    (Actually, len(out) is a little less than that, (len(in)-1)*times + 1 .)
    The spline is local: for example, out[ 3 to 4 ] is calculated
    from the 4 nearest input points at 2 3 4 5 .
Optional args:
    spline = "Catmull-Rom" default, or "Bspline";
        Catmull-Rom goes through the input points, Bspline does not
    border, default 1: duplicate the end points, pt[0] pt pt[-1]
        so that out[ 0 to 1 ] sees pt[0] pt[0] pt[1] pt[2].
    imatrix=M if it's constant over many calls:
        M = dotspline1matrix( times, spline )

dotspline2 is similar, over a 2d input array e.g. of pixels;
it just does dotspline1 on columns, then rows.

Notes:
    "Points" can be scalars, or vectors, or anything for which interpolation
    makes sense: audio samples in time, a butterfly's xyz and color too in time,
    pixels along a curve, rows of a matrix, whole matrices.

    There are many many kinds of splines, for curves and for surfaces,
    with a 5-foot shelf of books and papers (4M google hits x 1um ?) .
    The uniform local splines here are about the simplest.

How it works:
    split each output time = integer T + fraction t, 0 <= t < 1.
    Call the 4 nearest input points p_1 p0 p1 p2 = pt[ T-1 T T+1 T+2 ].
    M is a fixed 4x4 "blending matrix", Catmull-Rom or Bspline, see the code.
    Then out[time] =
        [p_1 p0 p1 p2] . M . [1 t t**2 t**3], the sum of 16 terms:

        [1 t t**2 t**3]
    p_1 [0 -1  2 -1]
    p0  [2  0 -5  3]
    p1  [0  1  4 -3]
    p2  [0  0 -1  1] / 2  # Catmull-Rom

    (Check this:
    t=0 -> p0
    t=1 -> p1
    t=.5 -> [p_1 p0 p1 p2] . [-1 9 9 -1] / 16
            [0 1 1 0] -> 9/8, overshoot
            [1 0 0 1] -> -1/8, negative for all-positive input.)

    Compute the 4 x times interpolation matrix  M . [1 t t**2 t**3] once
    for a block of times, e.g. 0 .1 .2 ... .9
    and dot it, slide it along, each group of 4 input points.
    (This is straightforward on paper, takes some fiddling with .T s in numpy.)

Exercises:
    Plot Catmull-Rom and Bspline through
        [0 1 0], [1 0 -1 0 ...], [1 -1 1 -1 ...] .
    Do linear and bilinear interpolation with numpy.dot .

See also:
    http://en.wikipedia.org/wiki/ B_spline Catmull-Rom Bicubic_interpolation
    scipy.interpolate.*Spline -- f2py fitpack, fast but beware
    scipy.ndimage spline_filter* zoom

"""
# speed: _dotblas.dot is pretty slow for short vecs, fast for long
# applets ? http://joningram.org/blog/2009/04/further-reading-on-splines ff


from __future__ import division
import numpy as np
dot = np.dot

__version__ = "2009-12-20 denis"
__reviews__ = ""
__author_email__ = "denis-bz-py at t-online dot de"

#...............................................................................
def dotspline1( x, times, spline="Catmull-Rom", border=1, imatrix=None ):
    """ simple local spline, x[ 0 1 2 ...] -> e.g. curve[ 0 .1 .2 ...] """
    if border:  # none same linear mirror
        x = np.concatenate([ np.array([x[0]]), x, np.array([x[-1]]) ])
    if imatrix is None:
        imatrix = dotspline1matrix( times, spline )  # 4 x times
    if np.isscalar(x[0]):
        dots = np.hstack([ dot( x[j:j+4], imatrix ) for j in xrange( len(x) - 3 )])
    else:
        dots = np.concatenate([ dot( x[j:j+4] .T, imatrix )
            for j in xrange( len(x) - 3 )], axis=-1 ) .T
        # assert x[0].shape == dots[0].shape
    return np.append( dots, np.array([x[-1]]), axis=0 )

blendmatrices = {
"Catmull-Rom" : np.array([  # interpolates but not convex, .5 -> [-1 9 9 -1] /16
    [0, -1, 2, -1],
    [2, 0, -5, 3],
    [0, 1, 4, -3],
    [0, 0, -1, 1]]) / 2,
"Bspline" : np.array([  # convex but doesn't interpolate, 0 -> [1 4 1] /6
   [1, -3, 3, -1],
   [4, 0, -6, 3],
   [1, 3, 3, -3],
   [0, 0, 0, 1]]) / 6
}

def dotspline1matrix( times, spline="Catmull-Rom" ):
    """ -> 4 x times interpolation matrix """
    blend = blendmatrices[spline]
    return np.array([ dot( blend, [1, t, t**2, t**3] )
        for t in np.linspace( 0, 1, num=times, endpoint=False )]) .T

#...............................................................................
def dotspline2( a, times, spline="Catmull-Rom", border=1, imatrix=None ):
    """ 2d interpolation:
        dotspline1 columns, N x N -> N x (times*N)
        then rows -> (times*N) x (times*N)
    """
    # slow for real use, but shows how separable kernels work
    # (faster: dot( 16 x times**2 imatrix, 16 points ) -- spline16, agg ?)
    # Mitchell and Netravali, "Reconstruction Filters in Computer Graphics",
    # Siggraph 1988: 2/3 CR + 1/3 Bspline

    if imatrix is None:
        imatrix = dotspline1matrix( times, spline )
    cols = dotspline1( a.T, times, border=border, imatrix=imatrix ) .T
        # print "test dotspline2 cols:", cols.shape, cols
    return dotspline1( cols, times, border=border, imatrix=imatrix )

#...............................................................................
if __name__ == "__main__":
    import sys
    N = 100
    N2 = 5
    times = 4
    spline = "Catmull-Rom"
    shape = 0  # (3,) (3,3)
    exec "\n".join( sys.argv[1:] )  # N= times= ...
    np.set_printoptions( 1, threshold=500, suppress=True )  # .1f

    def test_dotspline1( N=N, times=times, spline=spline, shape=shape ):
        print "dotspline1 N=%d times=%d spline=%s shape=%s" % (
            N, times, spline, shape )
        if shape:
            pt = np.zeros( (N,) + shape )
            pt[1] = np.ones(shape)
        else:
            pt = np.arange(N)
            pt[1] = 10
        M = dotspline1matrix( times, spline )
            # print "test dotspline1matrix:", M
        return dotspline1( pt, times, spline=spline, imatrix=M )

        # %timeit test_dotspline1( N=1000, times=10, shape=0 )
        #   10 loops, best of 3: 56 ms per loop (mac g4 ppc)
        # scipy.interpolate.UnivariateSpline N=1000 H=10: 13 msec

    def test_dotspline2( N=N2, times=times, spline=spline ):
        print "dotspline2 N=%d times=%d spline=%s" % ( N, times, spline )
        a = np.zeros(( N, N ))
        a[N//2,N//2] = 100
            # x = np.arange( N )
            # a = np.multiply( x, x[:,np.newaxis] )  # bilinear
        return dotspline2( a, times, spline=spline )

    if N > 0:
        t1 = test_dotspline1(N)
        print t1
    if N2 > 0:
        t2 = test_dotspline2(N2)
        print np.around(t2).astype(int)