Spline interpolation in 3D in python

Question

I am searching the equivalent Matlab command

Vq = interp3(X,Y,Z,V,Xq,Yq,Zq)

in Python. In Matlab I can use the method 'spline' interpolation, which I can not find in python for 3D data. There exists scipy.interpolate.griddata, but it doesn't have the option spline for 3D data.

The data I want to interpolate is a 3D matrix (51x51x51), which is regularly distributed on a 3D grid.

scipy.interpolate.Rbf may be the option, but I don't get it working:

xi = yi = zi = np.linspace(1, 132651, 132651) interp = scipy.interpolate.Rbf(xi, yi, zi, data, function='cubic')

leads to a memory error.

Edit: A minimal example of what I want (without interpolation): Matlab code

v=rand([51,51,51]);
isosurface (v, 0.3);

For simplicity, I use random data in this example. I want to make isosurface plots (in particular, Fermi surface plots). Since some structures are very small, a high grid resolution of 51x51x51 is needed.

A further comment: The data set in the matrix is independent from each other, z (or the 3rd component) is NOT a function of x and y.

Objectivism

Memory error may be related to the size of your data. Why are you using interpolate.Rbf? scipy equivalent to interp3 would more likely be griddata, docs.scipy.org/doc/scipy-0.14.0/reference/generated/…. – Savannahsavant

griddata supports (cubic) splines only in 2D (see your link). docs.scipy.org/doc/scipy/reference/generated/…, interpolate.Rbf doesn't seem to have this limitation. – Objectivism

True, haven't check the cubic. But Rbf seems to work - since you got a memory error, have you tried smaller datasets? – Savannahsavant

"I can't use scipy.interpolate.Rbf because I have 5,000+ datapoints." according to #39881247 . Since I have 132000 points, it is not an option. And in Matlab the spline interpolation is working with far higher datasets. Maybe ndimage.map_coordinates is an option ( docs.scipy.org/doc/scipy/reference/generated/… ), but I haven't understood it yet. – Objectivism

Answer 1

Spline interpolation on for 3+ dimensions can be done using scipy.interpolate.Rbf as your described. For plotting purposes you can use a smaller resolution (1000 points is a good rule of thumb), and when you want to evaluate your spline, you can interpolate on much greater than 132000 points without problem (see example below).

Can you add a Minimal, Complete, and Verifiable example for what you are trying to do in matlab? This will explain why do you need to create a grid space with a resolution of 132000 points. Also, please note, there is a curse of dimensionality. Matlab uses a cubic spline or a piecewise polynomial which can be dangerous due to overfitting. I recommend you used a more sane method for training on 51 datapoints and applying to 132000+ datapoints. This is a great example on polynomial curve fitting and model selection.

Example:

Generate data:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

%matplotlib inline
import random
# set seed to reproducible
random.seed(1)
data_size = 51
max_value_range = 132651
x = np.array([random.random()*max_value_range for p in range(0,data_size)])
y = np.array([random.random()*max_value_range for p in range(0,data_size)])
z = 2*x*x*x + np.sqrt(y)*y + random.random()
fig = plt.figure(figsize=(10,6))
ax = axes3d.Axes3D(fig)
ax.scatter3D(x,y,z, c='r')

Fit spline and interpolate

x_grid = np.linspace(0, 132651, 1000*len(x))
y_grid = np.linspace(0, 132651, 1000*len(y))
B1, B2 = np.meshgrid(x_grid, y_grid, indexing='xy')
Z = np.zeros((x.size, z.size))

import scipy as sp
import scipy.interpolate
spline = sp.interpolate.Rbf(x,y,z,function='thin_plate',smooth=5, episilon=5)

Z = spline(B1,B2)
fig = plt.figure(figsize=(10,6))
ax = axes3d.Axes3D(fig)
ax.plot_wireframe(B1, B2, Z)
ax.plot_surface(B1, B2, Z,alpha=0.2)
ax.scatter3D(x,y,z, c='r')

Fit spline on large data

predict_data_size = 132000
x_predict = np.array([random.random()*max_value_range for p in range(0,predict_data_size)])
y_predict = np.array([random.random()*max_value_range for p in range(0,predict_data_size)])
z_predict = spline(x_predict, y_predict)
fig = plt.figure(figsize=(10,6))
ax = axes3d.Axes3D(fig)
ax.plot_wireframe(B1, B2, Z)
ax.plot_surface(B1, B2, Z,alpha=0.2)
ax.scatter3D(x_predict,y_predict,z_predict, c='r')

Stucker

Thank you very much for your comment and sorry for my late response. I modified my question, to make it, hopefully, clearer. Summarised, I want to make isosurface plots. Can I do that with your answer? I think, if z_predict = spline(x_predict, y_predict) is not true - as it is the case for me- then your approach does not work, or did I get it wrong? – Objectivism

These splines are shockingly slow! – Abduct

My 8Gb RAM are not enough. I could run it with max_value_range = 13265 (ten less), x_grid = np.linspace(0, max_value_range, 1000*len(x)) (hundred less - same for y_grid) – Silva

If you don't run it with Jupyter, suppress the line %matplotlib inline and add plt.show() at the end of each cell. – Silva

yeah, I have the same problem as the comment above, ... but even with smaller values, the Jupyer notebook browser still fails with the memory. I think the example needs to have settings or arguments/code that consumes reasonable computing power? – Seclusion

Answer 2

The mess example worked great. I want to create a b spline of the data in 3d

data = np.array([[ 41. ,  57. ,  92. ],[ 39. ,  57.  , 92.4],[ 43. ,  57.  , 91.2], [ 23.,   47. , 119.6],
                 [ 27. ,  47. , 115.2], [ 25. ,  45. , 122. ], [ 25. ,  49. , 114. ],[ 29.,   49. , 109.6],
                 [ 29. ,  47. , 114.4], [ 27. ,  49. , 111.2], [ 23. ,  45. , 125.6], [ 31.,   49.,  106.8],
                 [ 25. ,  47. , 117.6], [ 39. ,  55. ,  95.6],[ 37.  , 53.  , 98.4], [ 35. ,  55. ,  96.8],
                 [ 33. ,  53. , 116.8], [ 23. ,  43. , 132.8], [ 25. ,  41. , 145.2],[ 25. ,  43.,  133.6],
                 [ 29. ,  51. , 106.4],[ 31.  , 53. , 121.2],[ 31., 51. , 104.8],[ 41.,   55.,   93.6],
                 [ 33. ,  51. , 103.6],[ 35.  , 53. ,  99.6],[ 37. ,  55. ,  96.4]])

x = data[:,0]
y = data[:,1]
z = data[:,2]

# sort data to avoid plotting problems
x, y, z = zip(*sorted(zip(x, y, z)))

x = np.array(x)
y = np.array(y)
z = np.array(z)

import scipy as sp
import scipy.interpolate
from mpl_toolkits.mplot3d import Axes3D

spline = sp.interpolate.Rbf(x,y,z,function='thin_plate',smooth=5, episilon=5)

x_grid = np.linspace(min(x),max(x), len(x))
y_grid = np.linspace(min(y),max(y), len(y))
B1, B2 = np.meshgrid(x_grid, y_grid, indexing='xy')

Z = spline(B1,B2)
fig = plt.figure(figsize=(10,6))
ax = Axes3D(fig)
ax.plot_wireframe(B1, B2, Z)
ax.plot_surface(B1, B2, Z,alpha=0.2)
ax.scatter3D(x,y,z, c='r')
plt.show()

Lovesome