在python中旋转一组3d坐标

Question

我有一组分子的 3d 坐标和一个穿过其质心的向量。我想旋转分子坐标和矢量以使矢量与 z 轴对齐。
我使用下面链接中的脚本来计算从向量到 z 轴的旋转矩阵，然后将相同的旋转矩阵应用于其他 3d 坐标： Calculate rotation matrix to align two vectors in 3D space?
但是这种方法使我的分子“变平”（旋转前为洋红色，旋转后为黄色）： Front view of molecule 和 Side view of molecule

有谁知道为什么这种方法不起作用？它在数学上正确吗？谢谢！

Answer 1

您链接到的问题this answer 中的方法对我来说似乎是正确的，并产生一个旋转矩阵（从将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 和x.shape == (n, 3)，而您的向量是v，形状为(3,)：

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)