我有一组使用 meshgrid() 生成的 3d 坐标。我希望能够围绕 3 个轴旋转它们。
我尝试解开网格并在每个点上进行旋转,但网格很大,内存不足。
This question在 2d 中使用 einsum() 解决了这个问题,但在将它扩展到 3d 时我无法弄清楚字符串格式。
我已经阅读了其他几页有关 einsum() 及其格式字符串的内容,但一直无法弄明白。
编辑:
我将网格轴称为 X、Y 和 Z,每个轴的形状均为 (213, 48, 37)。此外,当我试图将结果放回网格时,实际的内存错误出现了。
当我试图“解开”它以逐点旋转时,我使用了以下函数:
def mg2coords(X, Y, Z):
return np.vstack([X.ravel(), Y.ravel(), Z.ravel()]).T
我用以下循环遍历结果:
def rotz(angle, point):
rad = np.radians(angle)
sin = np.sin(rad)
cos = np.cos(rad)
rot = [[cos, -sin, 0],
[sin, cos, 0],
[0, 0, 1]]
return np.dot(rot, point)
旋转后,我将使用这些点进行插值。
最佳答案
使用您的定义:
In [840]: def mg2coords(X, Y, Z):
return np.vstack([X.ravel(), Y.ravel(), Z.ravel()]).T
In [841]: def rotz(angle):
rad = np.radians(angle)
sin = np.sin(rad)
cos = np.cos(rad)
rot = [[cos, -sin, 0],
[sin, cos, 0],
[0, 0, 1]]
return np.array(rot)
# just to the rotation matrix
定义一个示例网格:
In [842]: X,Y,Z=np.meshgrid([0,1,2],[0,1,2,3],[0,1,2],indexing='ij')
In [843]: xyz=mg2coords(X,Y,Z)
逐行旋转:
In [844]: xyz1=np.array([np.dot(rot,i) for i in xyz])
等价于einsum
逐行计算:
In [845]: xyz2=np.einsum('ij,kj->ki',rot,xyz)
它们匹配:
In [846]: np.allclose(xyz2,xyz1)
Out[846]: True
或者,我可以将 3 个数组收集到一个 4d 数组中,然后使用 einsum
旋转它。这里 np.array
在开头添加了一个维度。所以 dot
和 j
维度是 1st,数组的 3d 如下:
In [871]: XYZ=np.array((X,Y,Z))
In [872]: XYZ2=np.einsum('ij,jabc->iabc',rot,XYZ)
In [873]: np.allclose(xyz2[:,0], XYZ2[0,...].ravel())
Out[873]: True
1
和 2
类似。
或者,我可以将 XYZ2
分成 3 个分量数组:
In [882]: X2,Y2,Z2 = XYZ2
In [883]: np.allclose(X2,xyz2[:,0].reshape(X.shape))
Out[883]: True
如果您想在另一个方向旋转,请使用 ji
而不是 ij
,即使用 rot.T
。
关于python - 用于网格旋转的 Numpy einsum(),我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/31816754/