Gridbased multivariate cubic interpolation

Question

Say I am interested in

y = f(n, m, o)

I create discretized grids for n, m, o:

import numpy as np
nGrid = np.arange(0,10)
mGrid = np.arange(100,110)
oGrid = np.arange(1000,1010)
N,M,O = np.meshgrid(nGrid, mGrid, oGrid, indexing='ij')

For the sake of simplicity, let's assume f(n,m,o) = n+m+p:

Y = N+M+O

So, I have now a grid wise approximation of f inside Y. From my look into the documentation,

• I need to look at multivariate interpolation
• scipy.interpolate.Griddata looks great, but it is inteded for unstructured grids. That means I have to translate all my results into long arrays - this does not take advantage of me actually having structured grids
• The other methods listed under structured do not support cubic and multivariate (higher than 2 dimensions) inputs at the same time.

Given the variables that I already created, what is a good way to interpolate on this matter?

Update

This is the implementation of the first answer:

import numpy as np
nGrid = np.arange(0,10)
mGrid = np.arange(100,110)
oGrid = np.arange(1000,1010)
N,M,O = np.meshgrid(nGrid, mGrid, oGrid, indexing='ij')
Y = N+M+O

n,m,o = 5, 103, 1007

m_i = np.argmin(np.abs(np.floor(m)-mGrid))
m_f = m - m_i

n_i = np.argmin(np.abs(np.floor(n)-nGrid))
n_f = n - n_i

o_i = np.argmin(np.abs(np.floor(o)-oGrid))
o_f = o - o_i

A = Y[m_i-1:m_i+3, n_i-1:n_i+3, o_i-1:o_i+3]

# cubic smoothing polynome in 1D
# Catmull-Rom style
def v(p0, p1, p2, p3, f):
    return 0.5 * (p0 * (-f**3 + 2*f**2 - f) +
                  p1 * (3*f**3 - 5*f**2 + 2) +
                  p2 * (-3*f**3 + 4*f**2 + f) +
                  p3 * (f**3 - f**2))



B = v(A[0], A[1], A[2], A[3], m_f)
C = v(B[0], B[1], B[2], B[3], n_f)
D = v(C[0], C[1], C[2], C[3], o_f)

# D is the final interpolated array

print m+n+o, D

Unfortunately, f(n,m,o) = 1115, while D = 2215. Since there is no link to the methodology, it is difficult for me to understand what's exactly happening and why the approximation is so far off.

Answer 1

If you do not find a good one, roll your own!

But be warned, the tricubic interpolation is a bit tricky. The following describes a simple tricubic interpolation. There are much more complicated and faster methods available (matrix-based) depending on your need. One should also be aware of the fact that there are several cubic extrapolation methods, I will just use the one which is probably most common (Catmull-Rom).

When looking up a point (m,n,o) in a 3D array M, its surrounding for the tricubic interpolation consists of a 4x4x4-point (or 3x3x3-cell) grid. The points are defined on each axis by: floor(a)-1, floor(a),floor(a)+1, floor(a)+2 etc. for each axis. (Think of a Rubik's cube where point (m,n,o) is in the invisible center cubelet.)

The tricubic interpolation result is then a weighted average of these 64 points. The weights depend on the distance of the point from the grid points.

Let us define:

m_i = np.floor(m).astype('int')    # integer part
m_f = m - m_i                      # fraction

And similar definitions for n and o. Now we need to operate on the array:

A = M[m_i-1:m_i+3, n_i-1:n_i+3, o_i-1:o_i+3]

Our point is now at (1+m_f, 1+n_f, 1+o_f). (Note that this assumes the matrix to be infinite. Otherwise edge effects have to be taken into account.)

The 64 coefficients of the array can be calculated with some clever matrix maths, but here it is enough to know that the interpolation is associative, i.e. we can do it one axis at a time:

# cubic smoothing polynome in 1D
# Catmull-Rom style
def v(p0, p1, p2, p3, f):
    return 0.5 * (p0 * (-f**3 + 2*f**2 - f) +
                  p1 * (3*f**3 - 5*f**2 + 2) +
                  p2 * (-3*f**3 + 4*f**2 + f) +
                  p3 * (f**3 - f**2))


# interpolate axis-by-axis
B = v(A[0], A[1], A[2], A[3], m_f)
C = v(B[0], B[1], B[2], B[3], n_f)
D = v(C[0], C[1], C[2], C[3], o_f)

# D is the final interpolated value

So, the interpolation is calculated one axis at a time.

---

Let us make this into something a bit more practical. First, we need to do some address calculation from our own coordinate space into the array coordinates. If m-axis is defined:

m0: position of the first grid plane on the axis
m1: position of the last grid plane on the axis
n_m: number of samples

we can calculate the equivalent array position from the "real world" position m:

m_arr = (n_m - 1) * (m -  m0) / (m1 - m0)

The -1 comes from the fencepost phenomenon (if we have 10 samples, the distance between the first and the last is 9).

Now we can put everything together into a function:

# cubic smoothing polynome in 1D
# Catmull-Rom style
def interp_catmull_rom(p0, p1, p2, p3, f):
    return 0.5 * (p0 * (-f**3 + 2*f**2 - f) +
                  p1 * (3*f**3 - 5*f**2 + 2) +
                  p2 * (-3*f**3 + 4*f**2 + f) +
                  p3 * (f**3 - f**2))

# linear interpolation
# only needs p1 and p2, the rest are there for compatibility
def interp_linear(p0, p1, p2, p3, f):
    return (1-f)*p1 + f*p2


# don't interpolate, use the nearest point
# only needs p1 and p2, the rest are there for compatibility
def interp_nearest(p0, p1, p2, p3, f):
    if f > .5:
        return p2
    else:
        return p1

# 3D interpolation done axis-by-axis
def tri_interp(M, f, m, n, o, m0, m1, n0, n1, o0, o1):
    # M: 3D array
    # f: interpolation function to use
    # m,n,o: coordinates where to interpolate
    # m0: real world minimum for m
    # m1: real world maximum for m
    # n0,n1,o0,o0: as with m0 and m1

    # calculate the array coordinates
    m_arr = (M.shape[0] - 1) * (m - m0) / (m1 - m0)
    n_arr = (M.shape[1] - 1) * (n - n0) / (n1 - n0)
    o_arr = (M.shape[2] - 1) * (o - o0) / (o1 - o0)

    # if we are out of our data box, return a nan
    if m_arr < 0 or m_arr > M.shape[0] or \
       n_arr < 0 or n_arr > M.shape[1] or \
       o_arr < 0 or o_arr > M.shape[2]:
        return np.nan

    # calculate the integer parts and fractions
    m_i = np.floor(m_arr).astype('int')
    n_i = np.floor(n_arr).astype('int')
    o_i = np.floor(o_arr).astype('int')
    m_f = m_arr - m_i
    n_f = n_arr - n_i
    o_f = o_arr - o_i

    # edge effects may be nasty, we may need elements outside of the array
    # there may be more efficient ways to avoid it, but we'll create lists of
    # coordinates:

    n_coords, m_coords, o_coords = np.mgrid[m_i-1:m_i+3, n_i-1:n_i+3, o_i-1:o_i+3]

    # these coordinate arrays are clipped so that we are in the available data
    # for example, point (-1,3,7) will use the point (0,3,7) instead
    m_coords = m_coords.clip(0, M.shape[0]-1)
    n_coords = n_coords.clip(0, M.shape[1]-1)
    o_coords = o_coords.clip(0, M.shape[2]-1)

    # get the 4x4x4 cube:
    A = M[m_coords, n_coords, o_coords]

    # interpolate along the first axis (3D to 2D)
    B = f(A[0], A[1], A[2], A[3], m_f)

    # interpolate along the second axis (2D to 1D)
    C = f(B[0], B[1], B[2], B[3], n_f)

    # interpolate along the third axis (1D to scalar)
    D = f(C[0], C[1], C[2], C[3], o_f)

    return D

Now this function makes it possible to use any interpolation method which can be carried axis-by-axis for four neighbouring points. Cubic, linear and "nearest" are now shown.

---

But does it work? Let us test it by generating some random data and a plane intersecting the data. (The latter being surprisingly difficult.)

# random data
M = np.random.random((10, 12, 16))
m0,m1 = -10.,10.
n0,n1 = -7.,5.
o0,o1 = -4.,5.

# create a grid (grid from -15..15,-15..15 in n-m plane)
gr = mgrid[-15:15.01:.1, -15:15.01:.1]

# create two perpendicular unit vectors (forming a plane)
# vn: normal vector of the plane
# vp: some vector, whose projection determines one direction
# v0: unit vector on the plane (perpendicular to vn and vp)
# v1: unit vector on the plane (perpendicular to vn and v0)
vn = np.array([-.2, .3, 1])
vp = np.array([0, -1, 0])
v1 = np.cross(vn, vp)
v2 = np.cross(vn, v1)
v1 /= numpy.linalg.norm(v1)
v2 /= numpy.linalg.norm(v2)

# the grid and the unit vectors define the 3d points on the plane
gr3d = gr[0][:,:,None] * v1 + gr[1][:,:,None] * v2

# now we can fetch the points at grid points, just must flatten and reshape back
res = [ tri_interp(M, interp_catmull_rom, p[0], p[1], p[2], m0,m1,n0,n1,o0,o1) for p in gr3d.reshape(-1,3) ]
res = np.array(res).reshape((gr3d.shape[0], gr3d.shape[1]))

So, we have a plane which goes through our brick of data. Now we can actually see the difference between three different interpolation methods by choosing interp_linear, interp_nearest, or interp_catmull_rom as the 1D interpolation function. All pictures below are in the same scale both dimensionally and color-wise.

1. Nearest: Take the nearest known point in the grid (quick'n'dirty)

2. Linear: Take a 2x2x2 cube and use linear interpolation to calculate the result; the results are always between existing data points, no edge trickery is needed. However, the first differentials are not smooth, which can be seen in the image and is sometimes a problem in signal processing.

3. Cubic: Take a 4x4x4 cube and use cubic interpolation to calculate the result. The result has continuous first differentials, i.e. the result is "smooth" in all directions. There may be edge effects due to the large sample volume. (The interpolation algorithm repeats the edgemost voxels, but they may be mirrored, taken as some constant, etc.)

---

The differences in the nature of these interpolation methods may be even better seen if we take a line segment through the block:

 # line segment from (-12,-12,-12) to (12,12,12)
 v = np.array([np.linspace(-12,12,1000)]*3).T

 res = [ tri_interp(M, interp_catmull_rom, p[0], p[1], p[2], m0,m1,n0,n1,o0,o1) for p in v ]

For the different interpolation functions, we have the following values along the line segment:

The blockiness of the "nearest" and the kinks in the linear interpolation are clearly visible. It may look surprising that the "linear" interpolation does not give linear results. It can be best understood by thinking a 2D interpolation with a square with corner values 1,0,1,0 (i.e. diagonally opposite corners have the same values). There is no way to interpolate that so that all sections would be linear.

---

The algorithm shown here is a very trivial algorithm which can be made faster with several tricks. Especially if several points are needed in one cell the matrix-based approaches are much faster. Unfortunately, there are no very fast methods, and the matrix-based methods are more complicated to explain.

Answer 2

One way I found is ndimage.map_coordinates. It's a bit unfortunate it's hidden and not linked in the interpolate module, but it works amazingly well and offers many orders of approximation.

The two big downsides:

• it is a bit fishy at the grid ends, mode='nearest' will make it pick the border value, while it otherwise would assume values of 0 for the points outside of the grid.
• it is not really efficient.

I'm still open for better suggestions, this is just to report my own status quo - and perhaps this will help someone at some time.

import numpy as np
nGrid = np.arange(0,10)
mGrid = np.arange(100,110)
oGrid = np.arange(1000,1010)
N,M,O = np.meshgrid(nGrid, mGrid, oGrid, indexing='ij')
Y = N+M+O
n,m,o = 0.01, 103, 1007
from scipy import ndimage
data = Y
# subtract lower end of grid
coords = np.array([[n-0,m-100,o-1000]])
ndimage.map_coordinates(data, coords.T, order=2, mode='nearest')