使用numpy/scipy对3D体积进行转换和重采样

7

更新:

我创建了一个有良好文档的ipython笔记本。 如果你只是想要代码,请看第一个答案。

问题

我有一个40x40x40的灰度值体积。 这需要旋转/移位/剪切。

这里有一些有用的齐次变换集合:http://www.lfd.uci.edu/~gohlke/code/transformations.py.html

我需要把体积中的每个体素视为一对(位置向量,值)。 然后我会对每个坐标进行变换并从变换向量集中采样新值。

采样似乎相当困难,我很高兴找到了这个: https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.ndimage.affine_transform.html#scipy.ndimage.affine_transform

给定的矩阵和偏移量用于通过仿射变换找到输出中每个点在输入中的对应坐标。使用插值法确定这些坐标处的输入值。在输入边界之外的点按给定模式填充。

听起来太棒了。

但使用方法非常棘手。 这里 有人正在使用该代码旋转图像。 他的旋转矩阵是2x2,所以不是在齐次坐标下的。 我尝试将一个在齐次坐标下的平移矩阵(2D)传递给该函数:

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)

哪些会导致错误:

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

显然,这不能用于齐次坐标系。如何使用它来移动数据?

而且这仅适用于2D,对于3D,我甚至无法旋转体积:

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)

错误信息是相同的:仿射矩阵行数有误 可以使用这种方法来转换我的体积吗?
我找到了一组辅助方法,它们提供了移位和旋转,但没有剪切: https://docs.scipy.org/doc/scipy-0.14.0/reference/ndimage.html 但我更喜欢使用自定义变换矩阵。

请注意,自1.0版本以来,“affine_transform”支持齐次坐标。 - Juan
1个回答

4
我发现了另一个选项:map_coordinates 使用numpy可以生成坐标的网格,并将它们重塑/堆叠以形成位置向量。这些向量被转换并转换回网格坐标格式。最后,使用map_coordinates解决了采样问题。
我认为这是一个常见的问题,我创建了一个ipython笔记本,逐步解释了一切:

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

仍有一个问题:坐标的顺序很奇怪。你需要以一种不直观的方式重新排序网格,可能是我的代码中的错误。
请注意,这种坐标的重新排序会影响转换的轴。如果你想围绕x轴旋转某物,对应的向量不是(1,0,0),而是(0,1,0),真的很奇怪。
但它确实起作用,我认为原则是清楚的。

我认为奇怪的行为是由于(r,c,z)约定是(y,x,z),因此,绕x坐标旋转的相应向量是(0,1,0)。 - Ariel
什么是 (r,c,z) 约定?它的来源是什么?@Ariel - palimboa
它代表r=row,c=col和其他维度,然而,空间x、y对应于x=col和y=row。这类似于numpy meshgrid输出的笛卡尔('xy',默认)或矩阵('ij')索引。 - Ariel

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接