6

I have a set of 3d coordinates that was generated using meshgrid(). I want to be able to rotate these about the 3 axes.

I tried unraveling the meshgrid and doing a rotation on each point but the meshgrid is large and I run out of memory.

This question addresses this in 2d with einsum(), but I can't figure out the string format when extending it to 3d.

I have read several other pages about einsum() and its format string but haven't been able to figure it out.

EDIT:

I call my meshgrid axes X, Y, and Z, each is of shape (213, 48, 37). Also, the actual memory error came when I tried to put the results back into a meshgrid.

When I attempted to 'unravel' it to do point by point rotation I used the following function:

def mg2coords(X, Y, Z):
    return np.vstack([X.ravel(), Y.ravel(), Z.ravel()]).T

I looped over the result with the following:

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)

After the rotation I will be using the points to interpolate onto.

4

1 回答 1

6

使用您的定义：

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维度是第一个，数组的 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代替，即使用.ijrot.T

于 2015-08-04T21:30:04.310 回答