Python中3D数组的任意切片

Question

我将3D体积数据存储在一个3维数组中，背景值为0，体积值为1。现在，我想获得此体积的任意剖面。我在这里阅读答案：How can an almost arbitrary plane in a 3D dataset be plotted by matplotlib?

但是似乎接受的答案是错误的，它生成xoy平面的映射坐标，而不是切片坐标。那么如何获得切片平面的正确形状呢？是否有将映射的形状转换为原始形状的方法？谢谢！

Answer 1

该问题可能已过时，但看来following function中的scipy.ndimage可能解决了您的问题。

scipy.ndimage.interpolation.rotate所做的是，它将整个3d数组绕3个轴中的任意一个旋转一定角度，将存储的值插入到新的“单元格”中。它还会相应地调整（扩展）新数组的大小，并使用您指定的值填充新的空白单元格。之后，您可以像平常一样进行切片，例如：array[:,sy//2,:]。

简而言之，这是在平行于z轴的对角线上的中心切口（为简单起见：]

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

您实际上可以再想一想，并通过绕不同轴旋转几次将其扩展到任意切片。

PS。如果感觉花费了太多时间，可以通过传递order=0（默认值为order=3）来牺牲插值的质量。这将运行得更快。