Created
March 8, 2018 12:35
-
-
Save jaimefrio/79ee7a9f8d7ff6ca10ae72633ec68c80 to your computer and use it in GitHub Desktop.
Simple spline interpolation with
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment