rotate a set of 3d coordinates in python

Question

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!

could you please provide some sample data (just a few rows), as well as the coordinates mean and std-deviation? — Pierre D, Dec 31 '21 at 16:53
@PierreD Hi, I'm new to Stack overflow and I don't seem to have enough "reputation" to edit my question as there is already an edit pending, do you know if there is an alternative way to provide data? I haven't found any answers online. The mean of x, y, z coordinates are 85.3097, 25.8701, 5.7173 and std are 9.8511, 8.1084, 13.4442 — ruiruiyang, Jan 04 '22 at 11:20

Pierre D · Accepted Answer · 2021-12-31T16:52:00.010

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)

Thank you Pierre, you are right, my data and my vector have different units, after change one of it, my problem is solved! thanks again — ruiruiyang, Dec 31 '21 at 17:03
Glad to hear! Please consider [What should I do when someone answers](https://stackoverflow.com/help/someone-answers). — Pierre D, Dec 31 '21 at 18:13

rotate a set of 3d coordinates in python

1 Answers1