在python中旋转一组三维坐标

Question

我有一组分子的三维坐标和一个通过它的质心的矢量。我想旋转分子坐标和矢量使矢量与z轴对齐。

我在下面的链接中使用一个脚本来计算从向量到z轴的旋转矩阵，然后将相同的旋转矩阵应用到其他3d坐标：计算旋转矩阵来对齐三维空间中的两个矢量？。

但这种方法使我的分子“扁平”(旋转前的分子和旋转后的黄色)：分子前视镜和分子侧视。

有人知道为什么这种方法不起作用吗？在数学上是正确的吗？谢谢!

Answer 1

在这个答案中，您所链接的问题的方法在我看来是正确的，并生成一个旋转矩阵(从无限个旋转矩阵集合中将vec1与vec2对齐)：

def rotation_matrix_from_vectors(vec1, vec2):
    """ Find the rotation matrix that aligns vec1 to vec2
    :param vec1: A 3d "source" vector
    :param vec2: A 3d "destination" vector
    :return mat: A transform matrix (3x3) which when applied to vec1, aligns it with vec2.
    """
    a, b = (vec1 / np.linalg.norm(vec1)).reshape(3), (vec2 / np.linalg.norm(vec2)).reshape(3)
    v = np.cross(a, b)
    c = np.dot(a, b)
    s = np.linalg.norm(v)
    kmat = np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]])
    rotation_matrix = np.eye(3) + kmat + kmat.dot(kmat) * ((1 - c) / (s ** 2))
    return rotation_matrix

这个旋转矩阵应该是正交的。

--可能是(也就是猜测)--您的数据发生了什么变化:不同的轴有相当大的差异(可能是不同的单位？)在这种情况下，首先应该在旋转之前对数据进行规范化。例如，假设您的原始数据是一个带有x.shape == (n, 3)的数组x.shape == (n, 3)，而您的向量是带有形状(3,)的v

u, s = x.mean(0), x.std(0)
x2 = (x - u) / s
v2 = (v - u) / s

现在，尝试在x2上应用您的旋转，将v2与[0,0,1]对齐。

下面是一个用于说明的玩具示例：

n = 100
x = np.c_[
    np.random.normal(0, 100, n),
    np.random.normal(0, 1, n),
    np.random.normal(4, 3, n),
]
v = np.array([1,2,3])
x = np.r_[x, v[None, :]]  # adding v into x so we can visualize it easily

无归一化的

A = rotation_matrix_from_vectors(np.array(v), np.array((0,0,1)))
y = x @ A.T

fig, axes = plt.subplots(nrows=2, ncols=2)
for ax, (a, b) in zip(np.ravel(axes), combinations(range(3), 2)):
    ax.plot(y[:, a], y[:, b], '.')
    ax.plot(y[-1, a], y[-1, b], 'ro')
    ax.set_xlabel(a)
    ax.set_ylabel(b)
axes[1][1].set_visible(False)

具有先验归一化的

u, s = x.mean(0), x.std(0)
x2 = (x - u) / s
v2 = (v - u) / s

A = rotation_matrix_from_vectors(np.array(v2), np.array((0,0,1)))
y = x2 @ A.T

fig, axes = plt.subplots(nrows=2, ncols=2)
for ax, (a, b) in zip(np.ravel(axes), combinations(range(3), 2)):
    ax.plot(y[:, a], y[:, b], '.')
    ax.plot(y[-1, a], y[-1, b], 'ro')
    ax.set_xlabel(a)
    ax.set_ylabel(b)
axes[1][1].set_visible(False)