CODEWITHSUNDEEP
pythonmatplotlib
ruiruiyang
ruiruiyang

Reputation: 13

rotate a set of 3d coordinates in python

I have a set of 3d coordinates of a molecule and a vector passing through its center of mass. I want to rotate the molecule coordinates and the vector to make the vector align with z-axis.
I used a script in the link below to calculate the rotation matrix from the vector to z-axis, then apply the same rotation matrix to the other 3d coordinates: Calculate rotation matrix to align two vectors in 3D space?
But this method makes my molecule "flat" (in magentas the molecule before rotation and in yellow after rotation): Front view of molecule and Side view of molecule

Does anyone know why this method doesn't work? Is it mathematically correct? Thank you!

Upvotes: 1

Views: 2358

Answers (1)

Pierre D
Pierre D

Reputation: 26201

The method in this answer of the question you linked to seems correct to me, and produces one rotation matrix (from the infinite set of rotation matrices that will align vec1 to 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

That rotation matrix is orthonormal, as it should be.

Perhaps (aka, wild guess) what's happening with your data is that the various axes have considerably different variance (perhaps different units?) In this case, you should first normalize your data before rotation. For example, say your original data is an array x with x.shape == (n, 3) and your vector is v with shape (3,):

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

Now, try to apply your rotation on x2, aligning v2 to [0,0,1].

Here is a toy example for illustration:

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

Without normalization

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)

With prior normalization

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)

Upvotes: 1

Related Questions