Transforming and Resampling a 3D volume with numpy/scipy

Question

UPDATE:

I created a well documented ipython notebook. If you just want the code, look at the first answer.

Question

I've got a 40x40x40 volume of greyscale values. This needs to be rotated/shifted/sheared.

Here is a useful collection of homogeneous transformations: http://www.lfd.uci.edu/~gohlke/code/transformations.py.html

I need to treat every voxel in my volume like a pair of (position vector, value). Then I would transform the position and sample new values for each coordinate from the set of transformed vectors.

The sampling seems rather difficult, and I was glad to find this: https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.ndimage.affine_transform.html#scipy.ndimage.affine_transform

The given matrix and offset are used to find for each point in the output the corresponding coordinates in the input by an affine transformation. The value of the input at those coordinates is determined by spline interpolation of the requested order. Points outside the boundaries of the input are filled according to the given mode.

Sounds perfect.

But the usage is very tricky. Here someone is using that code for rotating an image. His rotation matrix is 2x2, so that's not in homogenous coordinates. I tried passing a translation matrix in homogenous coordinates (2D) to the function:

dim =10
arr=np.zeros((dim,dim))
arr[0,0]=1
mat=np.array([[1,0,1],[0,1,0],[0,0,1]])
out3=scipy.ndimage.affine_transform(arr,mat)
print("out3: ",out3)

Which produces an error:

Traceback (most recent call last):
  File "C:/Users/212590884/PycharmProjects/3DAugmentation/main.py", line 32, in <module>
    out3=scipy.ndimage.affine_transform(arr,mat)
  File "C:\Users\212590884\AppData\Local\Continuum\Anaconda2\lib\site-packages\scipy\ndimage\interpolation.py", line 417, in affine_transform
    raise RuntimeError('affine matrix has wrong number of rows')
RuntimeError: affine matrix has wrong number of rows

Apparently this doesn't work with homogeneous coordinates. How can I use this to shift the data?

And this was just in 2D, in 3D I can't even rotate the volume:

dim =10
arr=np.zeros((dim,dim,dim))
arr[0,0]=1

angle=10/180*np.pi
c=np.cos(angle)
s=np.sin(angle)
mat=np.array([[c,-s,0,0],[s,c,0,0],[0,0,1,0],[0,0,0,1]])
out3=scipy.ndimage.affine_transform(arr,mat)
print("out3: ",out3)

The error message is the same: affine matrix has wrong number of rows

Is it possible to use this method to transform my volume ?

I found a collection of helper methods, they offer shift and rotate but not shear: https://docs.scipy.org/doc/scipy-0.14.0/reference/ndimage.html

But I would prefer to use a custom transformation matrix.

Answer 1

I've found another option: map_coordinates

With numpy it is possible to generate a meshgrid of coordinates, then reshape/stack them to form position vectors. These vectors are transformed and converted back into the meshgrid coordinate format. Finally with map_coordinates the sampling problem is solved.

I think that this is a common problem and have created an ipython notebook which explains everything step by step:

http://nbviewer.jupyter.org/gist/lhk/f05ee20b5a826e4c8b9bb3e528348688

There is still one problem: The order of the coordinates is strange. You need to reorder the meshgrids in an unintuitive way. Could be a bug in my code.

Please be aware that this reordering of coordinates influences the axes of the transformations. If you want to rotate something around the x axis, the corresponding vector is not (1,0,0) but (0,1,0), it's really strange.

But it works, and I think the principle is clear.