I have following code used to interpolate 3D volume data.
Y, X, Z = np.shape(volume)
xs = np.arange(0, X)
ys = np.arange(0, Y)
zs = np.arange(0, Z)
points = list(zip(np.ravel(result[:, :, :, 1]), np.ravel(result[:, :, :, 0]), np.ravel(result[:, :, :, 2])))
interp = interpolate.RegularGridInterpolator((ys, xs, zs), volume,
bounds_error=False, fill_value=0, method='linear')
new_volume = interp(points)
new_volume = np.reshape(new_volume, (Y, X, Z))
This code takes about 37 seconds to execute on 512x512x110 volume (about 29 millions of points), which results in more than one microsecond per voxel (which is unacceptable amount of time for me - what is more it uses 4 cores). Call new_volume=interp(points)
takes about 80% of the prodecure time and the list creation almost whole remaining time.
Is there any simple (or even more complex) way to make this computation faster? Or is there any good Python library, which provides faster interpolation? My volume and points change in every call to this prodecure.
Here is slightly modified version of your cython
solution:
import numpy as np
cimport numpy as np
from libc.math cimport floor
from cython cimport boundscheck, wraparound, nonecheck, cdivision
DTYPE = np.float
ctypedef np.float_t DTYPE_t
@boundscheck(False)
@wraparound(False)
@nonecheck(False)
def interp3D(DTYPE_t[:,:,::1] v, DTYPE_t[:,:,::1] xs, DTYPE_t[:,:,::1] ys, DTYPE_t[:,:,::1] zs):
cdef int X, Y, Z
X,Y,Z = v.shape[0], v.shape[1], v.shape[2]
cdef np.ndarray[DTYPE_t, ndim=3] interpolated = np.zeros((X, Y, Z), dtype=DTYPE)
_interp3D(&v[0,0,0], &xs[0,0,0], &ys[0,0,0], &zs[0,0,0], &interpolated[0,0,0], X, Y, Z)
return interpolated
@cdivision(True)
cdef inline void _interp3D(DTYPE_t *v, DTYPE_t *x_points, DTYPE_t *y_points, DTYPE_t *z_points,
DTYPE_t *result, int X, int Y, int Z):
cdef:
int i, x0, x1, y0, y1, z0, z1, dim
DTYPE_t x, y, z, xd, yd, zd, c00, c01, c10, c11, c0, c1, c
dim = X*Y*Z
for i in range(dim):
x = x_points[i]
y = y_points[i]
z = z_points[i]
x0 = <int>floor(x)
x1 = x0 + 1
y0 = <int>floor(y)
y1 = y0 + 1
z0 = <int>floor(z)
z1 = z0 + 1
xd = (x-x0)/(x1-x0)
yd = (y-y0)/(y1-y0)
zd = (z-z0)/(z1-z0)
if x0 >= 0 and y0 >= 0 and z0 >= 0:
c00 = v[Y*Z*x0+Z*y0+z0]*(1-xd) + v[Y*Z*x1+Z*y0+z0]*xd
c01 = v[Y*Z*x0+Z*y0+z1]*(1-xd) + v[Y*Z*x1+Z*y0+z1]*xd
c10 = v[Y*Z*x0+Z*y1+z0]*(1-xd) + v[Y*Z*x1+Z*y1+z0]*xd
c11 = v[Y*Z*x0+Z*y1+z1]*(1-xd) + v[Y*Z*x1+Z*y1+z1]*xd
c0 = c00*(1-yd) + c10*yd
c1 = c01*(1-yd) + c11*yd
c = c0*(1-zd) + c1*zd
else:
c = 0
result[i] = c
The results are still identical to yours. With a random grid data of 60x60x60
I obtain the following timings:
SciPy's solution: 982ms
Your cython solution: 24.7ms
Above modified cython solution: 8.17ms
So its nearly 4 times faster than your cython
solution. Note that
@boundscheck(False)
.range
instead of prange
in your for loop
. Hope this helps.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With