Rotate $z = 0$ plane in 3D

Question

I have 100 points on $z=0$ plane. I want to rotate those points, such that they lie on any plane $P(a,b,c,d)$, preserving distances.

Hence, I need a rotation matrix.

For instance, if my points are $(2,5,0)$, $(4,4,0)$, $(1,3,0)$ and $(3,5,0)$, how do I rotate them and put them onto $x + 3y - z + 1 = 0$ plane?

This is what I have tried:
First, I have calculated the angle between the normal vector of $P$ and all the normal vectors of $x=0$, $y=0$ and $z=0$ plane $(1,0,0)$, $(0,1,0)$, $(0,0,1)$
Let us refer to those angles as $\alpha$, $\beta$ and $\theta$ respectively.
Now, I rotate my points $\alpha$ around $x-axis$, $\beta$ around $y-axis$ and $\theta$ around $z-axis$.

I think my approach is wrong. Distances are preserved but the new plane is not the plane I want.

Answer 1

To expand on my earlier comment, the approach I suggested does work and is relatively straightforward. See the plots produced by the Matlab script below.

function rotatePoints

  origin = [0 0 0];

  % Points on z = 0 plane
  pts = [[2 5 0];[4 4 0];[1 3 0];[3 5 0]];

  % Normal to z = 0 plane
  zaxis = [0 0 1];

  % Normal to x + 3y - z + 1 = 0 plane
  normal0 = [1 3 -1];
  normal = normal0/norm(normal0);

  % Angle and axis of rotation
  rot_angle = acos(dot(normal,zaxis))
  rot_axis = cross(normal, zaxis)

  % Rotation matrix
  R = rotation_matrix(rot_angle, rot_axis);

  % Rotate points
  [npts, m] = size(pts);
  pts_rotated = zeros(npts,3);
  for ii=1:npts
   pts_rotated(ii,:) = (R*pts(ii,:)')';
  end

  % Plot the two normals
  plot3([origin(1) zaxis(1)],[origin(2) zaxis(2)], [origin(3) zaxis(3)],'k'); hold on;
  plot3([origin(1) normal0(1)],[origin(2) normal0(2)], [origin(3) normal0(3)],'r');

  % Plot reference points
  for ii=1:npts
    plot3(pts(ii,1), pts(ii,2), pts(ii,3), 'ko'); hold on;
  end
  fill3(pts(:,1),pts(:,2),pts(:,3),'k');

  % Plot rotated points
  for ii=1:npts
    plot3(pts_rotated(ii,1), pts_rotated(ii,2), pts_rotated(ii,3), 'ro');
  end
  fill3(pts_rotated(:,1),pts_rotated(:,2),pts_rotated(:,3),'r','EdgeColor','r');
  grid on;
  axis equal;

  % Check distances
  for ii=1:npts-1
    xdiff = pts(ii+1,1)-pts(ii,1);
    ydiff = pts(ii+1,2)-pts(ii,2);
    zdiff = pts(ii+1,3)-pts(ii,3);
    dist(ii) = sqrt(xdiff^2+ydiff^2+zdiff^2);
  end
  for ii=1:npts-1
    xdiff = pts_rotated(ii+1,1)-pts_rotated(ii,1);
    ydiff = pts_rotated(ii+1,2)-pts_rotated(ii,2);
    zdiff = pts_rotated(ii+1,3)-pts_rotated(ii,3);
    dist_rotated(ii) = sqrt(xdiff^2+ydiff^2+zdiff^2);
  end
  dist_err = dist_rotated - dist

function [R] = rotation_matrix(angle, axis)

   axis = axis/norm(axis);
   ca = cos(angle);
   sa = sin(angle);
   I = [[1 0 0];[0 1 0];[0 0 1]];
   aa = axis'*axis
   A = [[0 axis(3) -axis(2)];[-axis(3) 0 axis(1)];[axis(2) -axis(1) 0]]
   R = (I - aa)*ca + aa + A*sa;

Answer 2

No need to compute angles anywhere.

The columns of the transformation matrix are the images of the unit vectors under said transformation. Since you didn't specify placement of your points in that plane, you can use any pair of orthogonal uint vectors in that plane. A nice trick to obtain orthogonal vectors in 3d is using the cross product.

So start with the normal vector $(1,3,-1)$ of the plane, and an arbitrarily chosen second vector $(0,0,1)$ Compute the cross products

$$v_1=\begin{pmatrix}1\\3\\-1\end{pmatrix}\times\begin{pmatrix}0\\0\\1\end{pmatrix} =\begin{pmatrix}3\\-1\\0\end{pmatrix} \\ v_2=\begin{pmatrix}1\\3\\-1\end{pmatrix}\times\begin{pmatrix}3\\-1\\0\end{pmatrix} =\begin{pmatrix}-1\\-3\\-10\end{pmatrix}$$

To preserve distances, you'll have to scale them to unit length. Then you can plug them into the columns of a transformation matrix, and choose an offset such that the image of the origin lies within your plane. For example like this:

$$\begin{pmatrix}x\\y\\z\end{pmatrix}\mapsto\begin{pmatrix} \frac{3}{\sqrt{10}} & -\frac{1}{\sqrt{110}} & \frac{1}{\sqrt{11}} \\ -\frac{1}{\sqrt{10}} & -\frac{3}{\sqrt{110}} & \frac{3}{\sqrt{11}} \\ 0 & -\frac{10}{\sqrt{110}} & -\frac{1}{\sqrt{11}} \end{pmatrix}\begin{pmatrix}x\\y\\z\end{pmatrix} +\begin{pmatrix}0\\0\\1\end{pmatrix}$$

If you really want to rotate around the axis where both your planes intersect, so that points on that axis remain fixed, you'll need a bit more work. First find the direction of the common line. It is perpendicular to both the normals of the two planes. Since I used $(0,0,1)$ as the arbitrary vector in the first step above, and that's the normal of the $z=0$ plane, you already have $v_1$ as the direction of that axis. You can get the transformation matrix like this:

$$\begin{pmatrix} \frac{3}{\sqrt{10}} & -\frac{1}{\sqrt{110}} & \frac{1}{\sqrt{11}} \\ -\frac{1}{\sqrt{10}} & -\frac{3}{\sqrt{110}} & \frac{3}{\sqrt{11}} \\ 0 & -\frac{10}{\sqrt{110}} & -\frac{1}{\sqrt{11}} \end{pmatrix}\begin{pmatrix} \frac{3}{\sqrt{10}} & \frac{1}{\sqrt{10}} & 0 \\ -\frac{1}{\sqrt{10}} & \frac{3}{\sqrt{10}} & 0 \\ 0 & 0 & 1 \end{pmatrix}^{-1}$$

The right matrix (prior to being inverted) will map the $x$ direction to the direction of the line of intersection, and it will map the $y$ direction to the direction perpendicular to that line within the $z=0$ plane. The $z$ axis is preserved, so this is a rotation around the $z$ axis. The inverse of that matrix will map the direction of the line of intersection to the $x$ direction, from where the left matrix will map it back to the line of intersection. You have to make sure that the second column of the right matrix has the correct sign, otherwise you'll end up describing a reflection (in the angle bisector plane) instead of a rotation. I did this at first, accidentially.

To get the offset right, simply plug in a point on the line of intersection (i.e. one which satisfies $x+3y+1=0$ and $z=0$), and compute the difference between the transformed and the original point. You want to subtract that difference, so that the point remains fixed.

In the end you obtain

$$\begin{pmatrix}x\\y\\z\end{pmatrix}\mapsto \frac1{110}\begin{pmatrix} -\sqrt{11} + 99 & -3 \, \sqrt{11} - 33 & 10 \, \sqrt{11} \\ -3 \, \sqrt{11} - 33 & -9 \, \sqrt{11} + 11 & 30 \, \sqrt{11} \\ -10 \, \sqrt{11} & -30 \, \sqrt{11} & -10 \, \sqrt{11} \end{pmatrix}\begin{pmatrix}x\\y\\z\end{pmatrix}- \frac1{110}\begin{pmatrix} \sqrt{11} + 11 \\ 3 \, \sqrt{11} + 33 \\ 10 \, \sqrt{11} \end{pmatrix}$$