## n is the degree.  If n == 3, the b-spline is cubic.
## m := len( P ) - 1
## len( knots ) == len( P ) + n + 1
## t is the parameter along the curve.  t is in [0,1].

def b_j_n( knots, j, n, t ):
    if n == 0:
        if knots[j] <= t and t < knots[j+1]: return 1.
        else: return 0.
    
    lhs_coeff_numer = t - knots[j]
    lhs_coeff_denom = knots[j+n] - knots[j]
    if lhs_coeff_numer != 0. and lhs_coeff_denom != 0.:
        lhs = ( lhs_coeff_numer / lhs_coeff_denom ) * b_j_n( knots, j, n-1, t )
    else:
        lhs = 0.
    
    rhs_coeff_numer = knots[j+n+1] - t
    rhs_coeff_denom = knots[j+n+1] - knots[j+1]
    if rhs_coeff_numer != 0. and rhs_coeff_denom != 0.:
        rhs = ( rhs_coeff_numer / rhs_coeff_denom ) * b_j_n( knots, j+1, n-1, t )
    else:
        rhs = 0.
    
    return lhs + rhs

def S( knots, n, P, t ):
    
    m = len( P ) - 1
    assert len( knots ) == len( P ) + n + 1
    
    from numpy import zeros
    result = zeros( len( P[0] ) )
    for i in xrange( len(P) ):
        result += P[i] * b_j_n( knots, i, n, t )
    
    return result

def uniform_knots( num_pts, degree ):
    '''
    A helper function for Uniform B-Splines.
    Returns an evenly spaced knot vector for passing to S().
    '''
    assert num_pts > 1
    
    from numpy import linspace
    return linspace( 0, 1, num_pts + degree + 1 )

def magic_uniform_knots( num_pts, degree ):
    '''
    Like uniform knots, but returns knots shifted and scaled
    such that t in [0,1] is evaluatable.
    '''
    
    assert num_pts + 1 > degree
    assert degree > 1
    
    knots = uniform_knots( num_pts, degree )
    
    n = degree
    
    k0 = n
    k1 = len( knots )-1 - n
    
    shift = -knots[k0]
    scale = 1. / ( knots[k1] - knots[k0] )
    
    knots += shift
    knots *= scale
    
    return knots


def Uniform_B_Spline( P, n ):
    '''
    Returns a function object that can be called with parameters between 0 and 1
    to evaluate the B-Spline with control points 'P' and degree 'n'.
    '''
    
    if len(P) <= n: return lambda t: P[ len(P) // 2 ]
    
    knots = magic_uniform_knots( len( P ), n )
    #knots = uniform_knots( len( P ), n )
    return lambda t: S( knots, n, P, t )


def S_3_i( P, i, t ):
    '''
    Returns the evaluation of a cubic b-spline S evaluated at index 'i', 1 <= i < len(P) - 2,
    and parameter t, 0 <= t <= 1.
    '''
    
    assert i - 1 >= 0
    assert i + 2 < len(P)
    
    from numpy import array, dot
    A = array( [[ -1, 3, -3, 1 ], [ 3, -6, 3, 0 ], [ -3, 0, 3, 0 ], [ 1, 4, 1, 0 ]] )
    ps = P[i-1:i+3]
    ts = array( [ t**3, t**2, t, 1 ] )
    
    return 1./6. * dot( ts, dot( A, ps ) )

def S_3( P, t ):
    '''
    Returns the evaluation of a cubic b-spline S evaluated at t, 0 <= t <= 1.
    '''
    
    from math import floor
    
    a = 1
    b = len(P) - 2
    
    texpand = a + t*(b-a)
    i = min( int( floor( texpand ) ), len(P) - 3 )
    tt = texpand - i
    
    return S_3_i( P, i, tt )

def Uniform_Cubic_B_Spline( P ):
    '''
    Returns a function object that can be called with parameters between 0 and 1
    to evaluate the uniform, cubic B-Spline with control points 'P'.
    '''
    
    if len(P) < 4: return lambda t: P[ len(P)//2 ]
    
    return lambda t: S_3( P, t )


def test1():
    print '====== test1() ======'
    
    from numpy import meshgrid, linspace
    P = meshgrid([0,0],range(10))[1]
    
    degree = 3
    bspline = Uniform_B_Spline( P, degree )
    V = 5
    for val in linspace( 0, 1, V ):
        print '%s: %s' % (val, bspline( val ))

def test2():
    print '====== test2() ======'
    
    from numpy import meshgrid, linspace
    P = meshgrid([0,0],range(10))[1]
    
    bspline = Uniform_Cubic_B_Spline( P )
    V = 5
    for val in linspace( 0, 1, V ):
        print '%s: %s' % (val, bspline( val ))

def main():
    test1()
    test2()

if __name__ == '__main__': main()