How to rotate a set of points on z = 0 plane in 3-D, preserving pairwise distances?

Question

I have a set of points double n[] on the plane z = 0.
And I have another set of points double[] m on the plane ax + by + cz + d = 0.
Length of n is equal to length of m.
Also, euclidean distance between n[i] and n[j] is equal to euclidean distance between m[i] and m[j].

I want to rotate n[] in 3-D, such that for all i, n[i] = m[i] would be true.

In other words, I want to turn a plane into another plane, preserving the pairwise distances.

Here's my code in java. But it does not help so much:

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

This is what I get. The point set on the right side is actually part of a point set on the left side. But they are on another plane.

Here's a 2-D representation of what I try to do:

There are two lines. The line on the bottom is the line I have. The line on the top is the target line. The distances are preserved (a, b and c).

Edit:
I have tried both methods written in answers. They both fail (I guess).

Method of Martijn Courteaux

   public static double[][] getRotationMatrix(double[] v0, double[] v1, double[] v2, double[] u0, double[] u1, double[] u2)
        {

            RealMatrix M1 = new Array2DRowRealMatrix(new double[][]{
                    {1,0,0,-1*v0[0]},
                    {0,1,0,-1*v0[1]},
                    {0,0,1,0},
                    {0,0,0,1}
            });

            RealMatrix M2 = new Array2DRowRealMatrix(new double[][]{
                    {1,0,0,-1*u0[0]},
                    {0,1,0,-1*u0[1]},
                    {0,0,1,-1*u0[2]},
                    {0,0,0,1}
            });

            Vector3D imX = new Vector3D((v0[1] - v1[1])*(u2[0] - u0[0]) - (v0[1] - v2[1])*(u1[0] - u0[0]), 
                    (v0[1] - v1[1])*(u2[1] - u0[1]) - (v0[1] - v2[1])*(u1[1] - u0[1]),
                    (v0[1] - v1[1])*(u2[2] - u0[2]) - (v0[1] - v2[1])*(u1[2] - u0[2])
            ).scalarMultiply(1/((v0[0]*v1[1])-(v0[0]*v2[1])-(v1[0]*v0[1])+(v1[0]*v2[1])+(v2[0]*v0[1])-(v2[0]*v1[1])));

            Vector3D imZ = new Vector3D(findEquation(u0, u1, u2));
            Vector3D imY = Vector3D.crossProduct(imZ, imX);

            double[] imXn = imX.normalize().toArray();
            double[] imYn = imY.normalize().toArray();
            double[] imZn = imZ.normalize().toArray();
            RealMatrix M = new Array2DRowRealMatrix(new double[][]{
                    {imXn[0], imXn[1], imXn[2], 0},
                    {imYn[0], imYn[1], imYn[2], 0},
                    {imZn[0], imZn[1], imZn[2], 0},
                    {0, 0, 0, 1}
            });
            RealMatrix rotationMatrix = MatrixUtils.inverse(M2).multiply(M).multiply(M1);
            return rotationMatrix.getData(); 
        }

Method of Sam Hocevar

static double[][] makeMatrix(double[] p1, double[] p2, double[] p3)
    {
        double[] v1 = normalize(difference(p2,p1));
        double[] v2 = normalize(cross(difference(p3,p1), difference(p2,p1)));
        double[] v3 = cross(v1, v2);
        double[][] M = { { v1[0], v2[0], v3[0], p1[0] },
                         { v1[1], v2[1], v3[1], p1[1] },
                         { v1[2], v2[2], v3[2], p1[2] },
                         {   0.0,   0.0,   0.0,   1.0 } };
        return M;
    }

    static double[][] createTransform(double[] A, double[] B, double[] C,
                         double[] P, double[] Q, double[] R)
    {
        RealMatrix c = new Array2DRowRealMatrix(makeMatrix(A,B,C));
        RealMatrix t = new Array2DRowRealMatrix(makeMatrix(P,Q,R));
        return MatrixUtils.inverse(c).multiply(t).getData();
    }

The blue points are the calculated points. The black lines indicate the offset from the real position.

Answer 1

New post, clean and well. I tried to implement Sam's answer, but slightly different, and it works. This approach is simpler than all the others, simply because it relies on you to make sure that your preserve the distance. If you don't, it will still match the points, but without distance and angles being preserved.

Here you have the check:

In code this would be (based on Sam's snippet):

static double[][] makeMatrix(double[] p1, double[] p2, double[] p3)
{
    double[] v1 = difference(p2,p1);
    double[] v2 = difference(p3,p1);
    double[] v3 = cross(v1, v2);
    double[][] M = { { v1[0], v2[0], v3[0], p1[0] },
                     { v1[1], v2[1], v3[1], p1[1] },
                     { v1[2], v2[2], v3[2], p1[2] },
                     {   0.0,   0.0,   0.0,   1.0 } };
    return M;
}

static double[][] createTransform(double[] A, double[] B, double[] C,
                     double[] P, double[] Q, double[] R)
{
    RealMatrix c = new Array2DRowRealMatrix(makeMatrix(A,B,C));
    RealMatrix t = new Array2DRowRealMatrix(makeMatrix(P,Q,R));
    return t.multiply(MatrixUtils.inverse(c)).getData();
}

static double[][] transformPoints(double[][] points, double[][] matrix)
{
    double[][] tfd = new double[points.length][];
    for (int i = 0; i < points.length; ++i)
    {
        tfd[i] = multiply(matrix, points[i]);
    }
    return tfd;
}

static double[] multiply(double[][] matrix, double[] point)
{
     double[] temp = new double[4];
     for (int i = 0; i < 3; ++i) temp[i] = point[i];
     temp[3] = 1.0;

     double[] res = new double[3];
     for (int i = 0; i < 3; ++i)
     {
         for (int j = 0; j < 4; ++j)
         {
             res[i] += temp[j] * matrix[i][j];
         }
     }

     return res;
}

Answer 2

If I understand correctly, you have three points A, B, C, and three points P, Q, R and you would like to know the affine transform (i.e. preserving distances) that transforms the first set into the second set.

Finding the transform is straightforward; you just need to create two orthonormal bases from your sets of points, fill matrixes with the basis vectors, and combine them.

Here is one way to do it (assuming you got some method to subtract two vectors):

double[][] makeMatrix(double[] p1, double[] p2, double[] p3)
{
    double[] v1 = normalize(p2 - p1);
    double[] v2 = normalize(cross(p3 - p1, v1));
    double[] v3 = cross(v1, v2);
    double[][] M = { { v1[0], v2[0], v3[0], p1[0] },
                     { v1[1], v2[1], v3[1], p1[1] },
                     { v1[2], v2[2], v3[2], p1[2] },
                     {   0.0,   0.0,   0.0,   1.0 } };
    return M;
}

mat4 createTransform(double[] A, double[] B, double[] C,
                     double[] P, double[] Q, double[] R)
{
    return mul(inverse(makeMatrix(A, B, C)),
               makeMatrix(P, Q, R));
}

Answer 3

Update:

This was much harder to solve than I thought. However, I think I managed it to fix it, after trying 5 different approaches. So, the requirements really are that you have 3 source points and 3 corresponding destination points. And that you are sure that the distance between the points is conserved. And finally that the source points define a plane.

(Note that all indices in Maple start at 1 instead of 0)

Your 3 source and destination points are represented like this:

First step is to create matrices that translate these points to the origin. The matrixes will be:

Here M1 will put the source points to the origin and M2 will do the same for the dest points, both using the first of the three points. Now we have two times a zero vector and two source vectors and two dest vectors:

The names refer to Source Points Centered (spC) and Destination Points Centered (dpC).

Next, we are going to find the linear combination of those two source vectors that produce the unit x vector (1, 0, 0). So we are looking for a solution for:

alpha * spC[0] + beta * spC[1] = (1, 0, 0)

Maple can solve this for us:

So, now, we have found a representation for the unit vector in function of the source vectors. We know that after applying the final linear transformation, the image of the x unit vector will be a linear combination of the destination vectors, again with the same coefficients we just found: alpha and beta. So, let's compute the image of the x unit vector after our linear transform:

So, at this point, we already know what the image of the x unit vector will be. We can also compute the normal vector of our three destination points. This will be the image of our z unit vector. This normal vector is given by the cross product of the two destination vectors.

So, now we know already two images of the unit vectors x and z. The image of the unit vector y can be found by taking into account that we want an orthonormal transformation. Just take the crossproduct of these two image vectors, to get the image of the y unit vector. Once we found the images, we have to normalize them of course by dividing by their length:

(Note: the requirement was that the distance was conserved. That means that normalizing the image of x would not be necessary.)

Now we can create the linear component of our transformation matrix by putting imXn, imYn and imZn in that order as column in a matrix:

You can see I homogenized the matrix to an affine matrix by adding zeroes at the right and bottom and a 1 at right bottom position of the matrix. To clarify:

a b c              a b c 0
d e f    becomes   d e f 0
g h i              g h i 0
                   0 0 0 1

Now we can compute the total matrix! The transformation is basically:

1. Translate the source points to the origin
2. Apply the linear transform
3. Translate the points to the dest position

This gives this matrix in total:

---

Old answer:

Well, this was actually a great exercise for my finals tomorrow. Here you have the technical part. If you want to, I can add explanation later:

Large image here.

This is done with Maple. So as you can see, if we let the normal vector of plane z = 0, which is <0, 0, 1> get through the transformation matrix, we have the expected normal vector for your ax + by + cz + d = 0 plane. If you assume that a^2 + b^2 + c^2 = 1, then it's even nicer.

Note that you have one degree of freedom: randomVector can be random and determines the rotation, because you didn't specify anything on how to rotate this thing.

If you make the randomVector less random, you can get a cleaner matrix, but I'm not sure if this will work for all kinds of target planes, because it can occur that the CrossProduct turns out to be zero, if the normalVector and the randomVector can be parallel:

Large image here.

Also note that you might want to add a translational component afterwards to make the plane match the equation ax + by + cz + d = 0, because if you transform all your points with this matrix, they are going to be in the plane ax + by + cz = 0 right now.