Search code examples
pythonscipyrotationpcapoint-clouds

Align PCA component with cartesian axis with rotation

I'm trying to rotate my point cloud such that the least significant PCA component is aligned with z-axis but with little success.

I first calculate the PCA components

U, S, Vt = np.linalg.svd(vertices - vertices.mean(axis=0), full_matrices=False)

but then I have trouble constructing rotation matrix, I tried scipy.spatial.transform.Rotation with from_rotvec() method and I'm not sure what I'm doing wrong since the results don't look as I would expect.

angles = np.arctan2(Vt[:, 2], np.array([0, 0, 1]))
rot = scipy.spatial.transform.Rotation.from_rotvec(angles)
new_vertices = np.dot(vertices, rot.T)
Solution

  • I generated an example data as follows

    import numpy as np;
import matplotlib.pyplot as plt

vertices = np.random.randn(10000, 2) / 2
vertices[:, 0] *= 3
vertices[:, 1] += vertices[:, 0] * 0.5;

    vc = vertices - vertices.mean(axis=0)
U, S, Vt = np.linalg.svd(vc)
vr = vc @ Vt.T

plt.figure(figsize=(10, 5))
plt.subplot(1,2,1)
plt.title('original vertices')
plt.scatter(vc[:, 0], vc[:, 1], alpha=0.1), plt.xlim([-6, 6]), plt.ylim([-6, 6])
plt.subplot(1,2,2)
plt.title('rotated vertices')
plt.scatter(vr[:, 0], vr[:, 1], alpha=0.1), plt.xlim([-6, 6]), plt.ylim([-6, 6])

    aligned points plot

    Basically X = U[:, :2] @ np.diag(S) @ Vt,

    np.allclose(U[:, :2] @ np.diag(S) @ Vt, vc)

    U is orthogonal, and S just scale the columns of U, and Vt applies the rotation. If we multiply both sides of the equation by inv(Vt) = Vt.T we get the aligned points.