How to rotate a plane in 3-D using standard form?

Question

I have a set of points $N = \{n_1, n_2, ...\}$ on a plane $z = 0$.
And I have another set of points $M = \{m_1, m_2, ...\}$ on plane $ax + by + cz + d = 0$.
$|N| = |M|$
$\forall n_i, n_j \in N$ and $\forall m_i, m_j \in M$, $d(n_i, n_j) = d(m_i, m_j)$ where $d(i,j)$ denotes the euclidean distance between $i$ and $j$.

I want to rotate the point set $N$ in 3-D, such that $\forall n_i \in N$ and $\forall m_i \in M$, $n_i = m_i$.

First, I use three points from $M$, namely $m_i, m_j, m_k$ to find $a, b, c, d$.
$\vec{PQ} = m_j - m_i$
$\vec{PR} = m_k - m_i$
$\vec{abc} = \vec{PQ} \times \vec{PR}$
$d = \vec{m_i} \cdot \vec{abc}$

Then I rotate all the points in $N$.

New coordinates of points $N$ are on $ax+by+cz = 0$ plane. But still, they are in a different plane in 3-D space.

How can I do the proper rotation? I also want to use $d$ component of the plane and make $M$ and $N$ identical.

As you can see in the image, there is a point cluster on the right. This cluster needs to be a part of the cluster on the left. I think I have to shift the points. But how?

Here are the results of another trial:

coordinates on z=0              estimated coordinates           actual coordinates
(-436.31, -270.11, 0.00)        (40.06, -460.13, 223.60)        (987.20, 968.48, 31.88)
(-451.92, -243.67, 0.00)        (62.38, -448.82, 241.41)        (966.37, 972.10, 9.61)
(-38.95, -50.98, 0.00)          (-11.96, -61.74, 12.67)         (873.16, 533.53, 91.34)
(-404.89, -228.80, 0.00)        (49.82, -410.18, 213.43)        (960.88, 924.93, 22.94)
(-434.90, -246.23, 0.00)        (53.25, -440.94, 229.13)        (970.11, 958.16, 18.91)
(-185.93, -64.57, 0.00)         (46.30, -157.20, 109.01)        (866.53, 668.09, 31.03)
(207.76, 51.25, 0.00)           (-63.82, 159.58, -127.48)       (826.81, 276.60, 147.47)
(-136.68, -124.92, 0.00)        (-10.76, -175.15, 59.10)        (915.75, 648.40, 88.10)
(-400.62, -203.17, 0.00)        (62.73, -387.99, 217.49)        (942.80, 911.10, 10.37)
(-430.38, -244.84, 0.00)        (52.02, -437.27, 226.43)        (969.61, 953.65, 20.22)
(575.78, -92.01, 0.00)          (-312.24, 262.21, -416.82)      (971.80, 9.06, 399.19)
(-465.18, -267.26, 0.00)        (54.70, -474.63, 244.03)        (982.02, 993.08, 16.90)
(-126.75, -134.57, 0.00)        (-20.80, -176.83, 49.73)        (923.85, 643.45, 98.18)
(300.67, 15.86, 0.00)           (-126.09, 186.09, -200.32)      (862.85, 208.75, 210.59)
(-450.28, -245.68, 0.00)        (60.48, -449.42, 239.75)        (968.01, 971.45, 11.51)
(587.99, -108.94, 0.00)         (-327.52, 256.24, -429.73)      (985.43, 4.95, 414.46)
(160.43, -56.50, 0.00)          (-104.85, 49.32, -124.52)       (899.52, 360.38, 186.76)
(-146.45, -150.19, 0.00)        (-20.98, -200.24, 58.89)        (932.95, 666.88, 97.93)
(-117.88, -189.20, 0.00)        (-56.39, -213.73, 28.86)        (964.40, 657.03, 133.32)
(568.87, -109.73, 0.00)         (-319.38, 244.58, -416.93)      (983.86, 22.07, 406.06)
(-153.17, -15.13, 0.00)         (60.15, -100.22, 100.13)        (834.40, 619.94, 18.11)
(-461.96, -245.70, 0.00)        (65.73, -456.18, 247.69)        (966.74, 981.76, 6.12)
(-76.37, -185.94, 0.00)         (-73.18, -187.21, 1.50)         (966.69, 619.25, 150.68)

Edit:

Both answers are great. But terminology does not match with mine. Please post an answer that shows the steps that I need to apply $\forall n_i \in N$

Edit2: Here is my Java code:

double[] rotate(double[] point, double[] currentEquation, double[] targetEquation)
{

    double[] currentNormal = new double[]{currentEquation[0], currentEquation[1], currentEquation[2]};
    double[] targetNormal = new double[]{targetEquation[0], targetEquation[1], targetEquation[2]};
    targetNormal = normalize(targetNormal);
    double angle = angleBetween(currentNormal, targetNormal);
    double[] axis = cross(targetNormal, currentNormal);
    double[][] R = getRotationMatrix(axis, angle);
    return rotated; 
}

double[][] getRotationMatrix(double[] axis, double angle)
{
    axis = normalize(axis);

    double cA = (float)Math.cos(angle);
    double sA = (float)Math.sin(angle);
    Matrix I = Matrix.identity(3, 3);
    Matrix a = new Matrix(axis, 3);


    Matrix aT = a.transpose();

    Matrix a2 = a.times(aT);

    double[][] B = 
        {
            {0, axis[2], -1*axis[1]},
            {-1*axis[2], 0, axis[0]},
            {axis[1], -1*axis[0], 0}
        };
    Matrix A = new Matrix(B);

    Matrix R = I.minus(a2);


    R = R.times(cA);
    R = R.plus(a2);
    R = R.plus(A.times(sA));
    return R.getArray();
}

Answer 1

You have two planes $W_N = (0,0,1,0)$ and $W_M = (a,b,c,d)$. For each plane, a coordinate system can be defined with origin on the point closest to the world origin $(0,0,0)$, and directions such that the local $z$ coordinate is out of plane.

The origin of the N plane is $\vec{r}_N = (0,0,0)$ and its coordinate directions $$E_N = \begin{bmatrix} 1 & 0& 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} $$

The origin and orientation of the M plane can be defined as $$\begin{align} \vec{r}_M & = \begin{pmatrix} \frac{-a d}{a^2+b^2+c^2} \\\frac{-b d}{a^2+b^2+c^2} \\ \frac{-c d}{a^2+b^2+c^2} \end{pmatrix} & E_M & = \begin{bmatrix} \frac{\sqrt{b^2+c^2}}{\sqrt{a^2+b^2+c^2}} & 0 & \frac{a}{\sqrt{a^2+b^2+c^2}} \\ \frac{-a b}{\sqrt{a^2+b^2+c^2}\sqrt{b^2+c^2}} & \frac{c}{\sqrt{b^2+c^2}} & \frac{b}{\sqrt{a^2+b^2+c^2}} \\ \frac{-a c}{\sqrt{a^2+b^2+c^2}\sqrt{b^2+c^2}} & \frac{-b}{\sqrt{b^2+c^2}} & \frac{c}{\sqrt{a^2+b^2+c^2}} \end{bmatrix} \end{align}$$

You can check that the local z axis is along the normal of the plane $\hat{k}_M = \frac{ \vec{r}_M}{ | \vec{r}_M | } = {\rm unit}(a,b,c)$

Also you can check that if the columns of $E$ are orthogonal and the last column is along the normal of the plane.

So now you take each point on N convert it to local coordinates $(u_i^N,v_i^N,0)$ such that the local to global transformation is $\vec{n}_i = \vec{r}_N + E\;(u_i^N,v_i^N,0)$ with $$(u_i^N,v_i^N,0) = E_N^\top ( \vec{n}_i - \vec{r}_N )$$ which is trivial. The $^\top$ denotes matrix transpose which is also the inverse for an orthonormal matrix like $E$.

On the M plane we do the world to local coordinate transformation with $$(u_i^M,v_i^M,0) = E_M^\top ( \vec{m}_i - \vec{r}_M )$$ which is not trivial. Or to go from local to world coordinates $\vec{m}_i = \vec{r}_N + E \;(u_i^M,v_i^M,0)$

Now if you make the identity transference with $(u_i^M,v_i^M,0) = (u_i^N,v_i^N,0)$ then the points on N will be mapped on M in the same arrangement relative to the plane origin (closest point to world origin) and orientation.

If you want to match the points in 3D with existing data, then you need to devise an affine transformation such that with $$ \begin{pmatrix} u_i^M \\ v_i^M \end{pmatrix} = \begin{pmatrix} u_0 \\ v_0 \end{pmatrix} + \begin{bmatrix} \cos \varphi & -\sin\varphi \\ \sin\varphi & \cos\varphi \end{bmatrix} \begin{pmatrix} u_i^N \\ v_i^N \end{pmatrix}$$ you get matching of $\vec{m}_i = \vec{r}_M + E_M\;(u_i^M,v_i^M,0)$.

So with picking a few points you find the origin and angle (and using the method posted here) of your affine transformation.

For example given two points on N $(u_1^N,v_1^N)$, $(u_2^N,v_2^N)$ and their corresponding points on M $(u_1^M,v_1^M)$, $(u_2^M,v_2^M)$ then their difference is

$$\begin{pmatrix} u_2^M \\ v_2^M \end{pmatrix} - \begin{pmatrix} u_1^M \\ v_1^M \end{pmatrix} = \begin{bmatrix} \cos \varphi & -\sin\varphi \\ \sin\varphi & \cos\varphi \end{bmatrix} \left( \begin{pmatrix} u_2^N \\ v_2^N \end{pmatrix}-\begin{pmatrix} u_1^N \\ v_1^N \end{pmatrix} \right) $$

The angle is found using the tan half angle rule as

$$ \varphi =- 2\, \tan^{-1} \left( \frac{(v_2^M-v_1^M)-(v_2^N-v_1^N)}{(u_2^M-u_1^M)-(u_2^N-u_1^N)} \right) $$

The translation is thus $(u_0,v_0) = (u_1^M,v_1^M) - (u_1^N \cos\varphi - v_1^N \sin\varphi, u_1^N \sin\varphi + v_1^N \cos\varphi)$

Summary

In the end for each point $\vec{n}_i$ in the N plane you transfer it to the M plane by

$$ \vec{m}_i = \vec{r}_M + E_M \left(\begin{pmatrix} u_0 \\ v_0 \\ 0\end{pmatrix} + \begin{bmatrix} \cos \varphi & -\sin\varphi &0 \\ \sin\varphi & \cos\varphi &0\\0&0&1 \end{bmatrix} \left( E_N^\top (\vec{n}_i - \vec{r}_N) \right ) \right) $$

Answer 2

First, your data simply doesn't satisfy your requirements perfectly, though it is close. This is not surprising given the (presumably) numerical estimates that you're working with. As a result, though, you simply cannot expect to find an isometry mapping the set $N$ to the set $M$. The best we can hope for is to find a function that comes close, perhaps by using a least squares estimate.

I guess we're looking for a function of the form $$ f(x,y) = \left(\begin{matrix} a_{1,1} & a_{1,2} \\ a_{2,1} & a_{2,2} \\ a_{3,1} & a_{3,3} \end{matrix}\right) \left(\begin{matrix} x \\ y \end{matrix}\right) + \left(\begin{matrix} s_1 \\ s_2 \\ s_3 \end{matrix}\right). $$ Our objective is to find constants that minimizes $$ \sum_i \left\|f(xn_i,yn_i) - (xm_i,ym_i,zm_i)\right\|, $$ where $(xn_i,yn_i)$ is a point in the input set $N$ and $(xm_i,ym_i,zm_i)$ is the corresponding point in the target set $M$. Fortunately, your data appears to be aligned in this fashion.

Implementing this plan in Mathematica, I found

$$ f(x,y) = \left(\begin{matrix} -0.449872 & 0.578342 \\ 0.578331 & 0.769321 \\ -0.680546 & 0.271442 \end{matrix}\right) \left(\begin{matrix} x \\ y \end{matrix}\right) + \left(\begin{matrix} 0.000242141 \\ 0.00258449 \\ -0.00249402 \end{matrix}\right), $$ with a minimum value of the norm of $0.134204$. Note that the matrix is close to orthogonal, as expected. The norms of the columns are $0.999997$ and $1.00001$ and their dot product is $0.0000133517$.

Alternatively, you can provide constraints to guarantee an orthogonal matrix but, then your fit might not be quite as good. Taking this approach, I got

$$ f(x,y) = \left(\begin{matrix} -0.449871 & 0.578337 \\ 0.578332 & 0.76931 \\ -0.680546 & 0.271443 \end{matrix}\right) \left(\begin{matrix} x \\ y \end{matrix}\right) + \left(\begin{matrix} -0.000238491 \\ 0.000989366 \\ -0.00246728 \end{matrix}\right), $$ with an error in the norm of $0.135831$.

I could provide Mathematica code, if you like, but (judging from your picture) you might be using another environment.