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? enter image description here

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();
}

If your new plane is parallel to the old one, you should get $n=a$, $m=b$, $k=c$, right? — Muphrid, May 12 '14 at 20:19
That's right. I missed it. But still, I need to shift my plane. — padawan, May 12 '14 at 20:31
Well, you have one unknown $d$, so you need one equation to determine it. Are there any criteria that the shifted plane must obey? — Muphrid, May 12 '14 at 20:49
Sorry, I don't understand. Either $d$ is given to you, or there is some other condition the shifted plane must satisfy (for example, maybe you're told the shifted plane should pass through a specific point). If $d$ is not given to you and you have no other condition, you would never be able to determine $d$. — Muphrid, May 13 '14 at 14:03
Some specific data would certainly make it easier to provide a good answer! Since three non-collinear points determine a plane, perhaps it would be sufficient to indicate three points in one plane that you would like to map to three points in another. — Mark McClure, May 21 '14 at 11:59
What do you mean by "estimated coordinates" and "actual coordinates"? — Mark McClure, May 21 '14 at 14:40
@MarkMcClure actual coordinates are coordinates of points in $M$. Estimated coordinates are the coordinaes of sensors $N$ after rotation. — padawan, May 21 '14 at 14:46
If "both answers are great", then they're certainly worth at least an upvote, if not an accept or the full bounty. Also, if you want more details, can you please describe the details you're looking for in comments to the answers themselves? My answer uses your data verbatim, for example, so I'm hard pressed to know what details you're looking for. It appears that you're working numerically in some programming environment. Perhaps code would help, but in what environment? It's nice for the question to have these kinds of details, if you expect detailed answers in return. — Mark McClure, May 21 '14 at 17:43
@MarkMcClure I already upvoted them but I think there were some problems with my browser. I'm working on Java. I'm going to post my code here as well. — padawan, May 21 '14 at 17:57

score 2 · Accepted Answer · edited Apr 13 '17 at 12:21

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

No, the two planes can be parallel to each other, and $\varphi \neq 0$. This angle is a planar coordinate rotation to convert between the two (arbitrary) chosen local coordinates for each plane. — John Alexiou, May 21 '14 at 18:49
$m_i$ in the Summary section is the transformed $n_i$ right? — padawan, May 24 '14 at 11:55
That is the coordinate offset on the plane. The point on the N plane with planar coordinates $(0,0)$ has to map to a point $(u_0, v_0)$ on the M plane. — John Alexiou, May 26 '14 at 04:39

score 1 · Answer 2 · answered May 21 '14 at 15:30

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.

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

2 Answers2

Summary

Linked