2

I have a 3D volume data stored in an 3-dimension array with background value 0 and volume value 1. Now I want to get arbitrary section planes of this volume. I read the answer here:How can an almost arbitrary plane in a 3D dataset be plotted by matplotlib?

But it seems that the accepted answer was wrong, it generates mapped coordinates of xoy plane but not the slice coordinates. So how to get the correct shape of the slice plane? Are there any methods of transforming the mapped shape to original shape? Thanks!

Silver0427
  • 63
  • 5
  • Do you mean "slice" in the python sense (using subscripts, for example to select a 2D subarray of a 3d array parallel to one the axes), or slice as in an interpolation of volumetric data along a plane (where the plane is not necessarily parallele to any axes)? The first of these is a lot simpler than the second – Alex I Apr 30 '18 at 09:06
  • @Alex Sorry I'm sorry I didn't express my problem well. Actually the problem should be the second situation, an interpolation method is necessary. – Silver0427 Apr 30 '18 at 14:21

1 Answers1

1

The question might be outdated, but looks like the following function from scipy.ndimage might solve your issue.

What scipy.ndimage.interpolation.rotate does is it rotates the full 3d array around either of the 3 axes by a certain angle interpolating the stored values to the new "cells". It also resizes (extends) the new array accordingly, filling the new empty cells by a value you specify. After that you can take slice as you would normally do, say: array[:,sy//2,:].

So in short, here's a central cut around a diagonal parallel to the z-axis (for simplicity):

sz, sy, sx = array.shape
array[:,sy//2,:] # this is a cut in the x-z plane ... 
                 # ... passing through the middle of the y-axis

# rotating by 45 degrees around the z axis ... 
# ... `(2,1)` is `(x,y)` which defines the rotation plane
array_rotated = scipy.ndimage.interpolation.rotate(array, angle=45, axes=(2,1))

sz, sy, sx = array_rotated.shape
# now you'll notice that `sz` is the same ... 
# ... but `sx` and `sy` have increased because the diagonal is longer

array_rotated[:,sy//2,:] # this slice is now in the new x'-z' plane ...
                         # ... but in terms of the original array ...
                         # ... it passes through the diagonal of x-y along z

You can actually think a bit further and extend this to an arbitrary slice by rotating it around different axes several times.

PS. If you feel it takes too much time, you can sacrifice the quality of interpolation by passing the order=0 (default is order=3). This will run faster.

hayk
  • 388
  • 1
  • 5
  • 20