Best way to rotate a 3D grid (nxnxn) of values in Python with interpolation?

Question

If I have a nxnxn grid of values, say 32x32x32, and I want to rotate this cube grid of values by some rotation angle in either the x, y, or z axes, and interpolate missing values, what would be the best way to go about doing this without using any existing algorithms from packages (such as Scipy)?

I'm familiar with applying a 3D rotation matrix to a 3D grid of points when it's represented as a [n, 3] matrix, but I'm not sure how to go about applying a rotation when the representation is given in its 3D form as nxnxn.

I found a prior Stack Overflow post about this topic, but it uses three for loops for its approach, which doesn't really scale in terms of speed. Is there a more vectorized approach that can accomplish a similar task?

Thanks in advance!

Answer 1

One way I could think of would look like this:

1. reshape nxnxn matrix to an array containing n-dimensional points
2. apply rotation on this array
3. reshape array back to nxnxn

Here is some code:

import numpy as np

#just a way to create some nxnxn matrix
n = 4
a = np.arange(n)
b = np.array([a]*n)
mat = np.array([b]*n)

#creating an array containg n-dimensional points
flat_mat = mat.reshape((int(mat.size/n),n))

#just a random matrix we will use as a rotation
rot = np.eye(n) + 2

#apply the rotation on each n-dimensional point
result = np.array([rot.dot(x) for x in flat_mat])
#return to original shape
result=result.reshape((n,n,n))
print(result)