Rotate normal vector onto axis plane

Question

I have a set of data points in 3D space which apparently all fall onto a specific plane. I use PCA to compute the plane parameters. The 3rd component of PCA gives me the normal vector of the plane (weakest component).

What I want to do next is to transform all the points onto said plane and look at it in 2D.

My idea was to do the following:

• Find a center point (average point) on the plane
• Substract it from all data points to arrange them around the origin
• Rotate the normal so that it becomes (0,0,-1)
• Apply this rotation to all data points
• Use orthogonal projection (basically, skip z axis)

Now I'm stuck at finding the right rotation operation. I tried working with acos or atan and setting up two rotation matrices. Seems both methods (using acos, using atan) give me the wrong result. Perhaps you can help me out here!

Matlab code follows:

b = atan(n(1) / n(2));
rotb = [cos(b) -sin(b) 0; sin(b) cos(b) 0; 0 0 1];
n2 = n * rotb;
a = atan(n(1) / n(3));
rota = [cos(a) 0 sin(a); 0 1 0; -sin(a) 0 cos(a)];
n3 = n2 * rotaows:

I expect n2 to have y component of zero. However that fails already for the vector (-0.6367, 0.7697, 0.0467).

Answer 1

Here's the accepted answer made in Python:

import numpy as np

def rotate(points, normal, around):
  # Let's rotate the points such that the normal is the new Z axis
  # Following https://stackoverflow.com/questions/1023948/rotate-normal-vector-onto-axis-plane
  old_x_axis = np.array([1, 0, 0])

  z_axis = normal
  y_axis = np.cross(old_x_axis, z_axis)
  x_axis = np.cross(z_axis, y_axis)
  
  axis = np.stack([x_axis, y_axis, z_axis])

  return np.dot(points - around, axis.T)

points = np.array([
    [0, 1, 1],
    [0, 1, 0.2],
    [1, 0, -7]
])

v1 = points[1] - points[0]
v2 = points[2] - points[0]

normal = np.cross(v1, v2)
print("rotated points", rotate(points, normal, points[0]))