Simple spline interpolation with

spline_interpolation.py

import numpy as np
from scipy import ndimage

SPLINE_ORDER = 3

def b_spline(x):
  """Computes 3rd order B-spline at x."""
  x = np.fabs(x)
  b_spline = np.zeros_like(x, dtype=float)
  mask = x < 1.0
  masked_x = x[mask]
  b_spline[mask] = masked_x * masked_x * (masked_x - 2.0) * 3.0 + 4.0
  mask = (x < 2.0) & ~mask
  masked_x = x[mask]
  masked_x -= 2.0
  b_spline[mask] = -masked_x * masked_x * masked_x
  b_spline /= 6.0

  return b_spline


def mirror_bounds(indices, length):
  indices = np.array(indices)
  mask = indices < 0
  indices[mask] *= -1
  mask = indices >= 2 * length - 2
  indices[mask] %= 2 * length - 2
  mask = indices >= length
  indices[mask] = 2 * length - 2 - indices[mask]
  return indices


def reflect_bounds(indices, length):
  indices = np.array(indices)
  mask = indices < -0.5
  indices[mask] = -indices[mask] - 1
  mask = indices >= 2 * length
  indices[mask] %= 2 * length
  mask = indices >= length - 0.5
  indices[mask] = 2 * length - indices[mask] - 1
  return indices


def interpolate(x, bspline_coefficients, x_bounds_fn, coefficients_bounds_fn):
  """Interpolates at x using the passed 3rd order spline coefficients.

  Parameters
  -----------
  x : array-like
    The coordinates at which to interpolate.
  spline_coefficients : array-like
    A 1-dimensional array of spline coefficients.
  x_bounds_fn : callable
    Function taking a double array and an integer length and returning a 
  coefficients_bounds_fn : callable
    Function taking an array of integers and an integer length, and converting
    them to in-bounds indices for an array of the required length.
  """
  bspline_coefficients = np.asarray(bspline_coefficients)
  n, = bspline_coefficients.shape  # This ensures the array is 1-D
  x = x_bounds_fn(x, n)
  shape = x.shape + (SPLINE_ORDER + 1,)

  bspline_weight_coords = np.empty(shape, dtype=float)
  bspline_weight_coords[:] = np.arange(-1, SPLINE_ORDER, dtype=float)
  bspline_weight_coords += (np.floor(x) - x)[:, None]
  bspline_weights = b_spline(bspline_weight_coords)

  coefficient_indices = np.empty(shape, dtype=np.intp)
  coefficient_indices[:] = np.arange(-1, SPLINE_ORDER)
  coefficient_indices += np.floor(x).astype(np.intp)[:, None]
  coefficient_indices = coefficients_bounds_fn(coefficient_indices, n)
  coefficients = bspline_coefficients[coefficient_indices]

  return np.einsum('...i,...i', coefficients, bspline_weights)