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).

Why don't you simply project all the points on the plane and then rotate everything so that you can just draw the points using their xz (or xy) coordinates? — Jasper Bekkers, Jun 21 '09 at 14:30
That, or attach the camera directly to the plane such that it's looking directly at it. — Jasper Bekkers, Jun 21 '09 at 14:31
The camera idea is a good one. You move the camera out from the plane from the origin some distance along the normal. Then you point the camera at the origin. Of course, this assumes that you have camera projection code, which can be done with with rotations or vectors. — Nosredna, Jun 21 '09 at 20:34
Good suggestion. Sadly enough, I don't have camera code there and I don't want to introduce any operations I do not fully control. — ypnos, Jun 21 '09 at 21:18

Nosredna · Accepted Answer · 2009-06-21T20:22:17.427

If you have a plane, you have a normal vector and an origin. I wouldn't do any "rotations" at all. You're just a few vector operations away from your answer.

Let's call your plane's normal vector the new z axis.
You can generate the new y axis by crossing the old x axis with the new z axis (your plane's normal).
Generate the new x axis by crossing the new z with the new y.
Make all your new axis vectors into unit vectors (length 1).
For every point you have, create a vector that's from your new origin to the point (vector subtraction of point - plane_origin). Just dot with the new x and new y unit vectors and you get a pair (x,y) you can plot!

If you have cross and dot product functions already, this is just a few lines of code. I know it works because most of the 3D videogames I wrote worked this way.

Tricks:

Pay attention to which directions your vectors are pointing. If they point the wrong way, negate the resultant vector or change the order of the cross product.
You have trouble if your plane's normal is exactly the same as your original x axis.

Fantastic answer! This fixes a lot of problems I've been having (dealing with rotational errors via trig) this whole time. I think I need to take a linear algebra course. — Kaiged, Aug 17 '12 at 04:28

score 1 · Answer 2 · answered Jun 21 '09 at 14:27

How about:

Decompose the normal vector into a vector in the XY-plane and a Z vector. Then apply a rotation around the Z-axis to line up the XY vector with one of the axis. Then find the dot product of the the normal with the Z-axis, and rotate along which ever of X,Y you lined up with.

The idea is to line the normal vector up with Z, and by doing that your plane is now the XY plane.

score 1 · Answer 3 · answered Mar 16 '21 at 15:41

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]))

score 0 · Answer 4 · answered Jun 21 '09 at 15:13

Although there were other interesting responses, this is the solution we figured out while waiting for answers:

function roti = magic_cosini(n)
    b = acos(n(2) / sqrt(n(1)*n(1) + n(2)*n(2)));
    bwinkel = b * 360 / 2 / pi;
    if (n(1) >= 0)
        rotb = [cos(-b) -sin(-b) 0; sin(-b) cos(-b) 0; 0 0 1];
    else
        rotb = [cos(-b) sin(-b) 0; -sin(-b) cos(-b) 0; 0 0 1];
    end
    n2 = n * rotb;
    a = acos(n2(3) / sqrt(n2(2)*n2(2) + n2(3)*n2(3)));
    awinkel = a * 360 / 2 / pi;
    rota = [1 0 0; 0 cos(-a) -sin(-a); 0 sin(-a) cos(-a)];
    roti = rotb * rota;

(It's returning a hopefully correct double rotation matrix)

The flaw we had before and fixed here was to esp. deal with the sign of the X component, which was not covered in the cosine computations. This made us rotate in the wrong direction once (rotating with 180° - angle).

I hope I will also find time to try Nosredna's solution! It's always good to avoid trigonometry.

I wanted to mention that any time you go the trig route rather than the vector route and you find yourself using atan, see if atan2 will work for you. It takes both x and y as arguments instead of taking a single argument. It puts the answer in the correct quadrant so you can avoid the usual conditional mess. http://en.wikipedia.org/wiki/Atan2 — Nosredna, Jun 21 '09 at 20:28
Yes, atan2 is awesome and it saved me lots of time before. A good hint. — ypnos, Jun 21 '09 at 21:19

Rotate normal vector onto axis plane

4 Answers4