4
\$\begingroup\$

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:

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

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

enter image description here

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

enter image description here

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

\$\endgroup\$
13
  • \$\begingroup\$ If your goal is that m[i] = n[i] for all i, well then copy the array, given that the lengths are equal as well. \$\endgroup\$ May 25, 2014 at 15:16
  • \$\begingroup\$ Oh, I think I see what you are trying. You actually don't have the target points. Well, then the best way to do this is to use coordinate transformation. Tomorrow I'm having an exam. If no one answered until then, I will try to formulate an answer. Basically the trick is to create a matrix that converts from one basis to another. en.wikipedia.org/wiki/Basis_(linear_algebra) \$\endgroup\$ May 25, 2014 at 15:20
  • \$\begingroup\$ @MartijnCourteaux I need the transformation/rotation matrix as well. \$\endgroup\$
    – padawan
    May 25, 2014 at 15:20
  • \$\begingroup\$ Can you maybe rewrite the question to say precisely what data you know and what you don’t? Right now it looks like you know a, b, c, d, but it seems that you actually don’t? Also in your 2D representation, the b distance isn’t actually preserved, is that an error? \$\endgroup\$ May 29, 2014 at 12:33
  • \$\begingroup\$ @SamHocevar I know a, b, c, d but there are infinitely many numbers of transformations of points from z=0 to ax+by+cz+d=0. I know three points on second plane as well. That reduces the last degree of freedom. \$\endgroup\$
    – padawan
    May 29, 2014 at 15:24

3 Answers 3

1
\$\begingroup\$

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.

enter image description here

Here you have the check:

enter image description here

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;
}
\$\endgroup\$
9
  • \$\begingroup\$ That code did the trick. Last modification I did is to shift the points on z=0 to positive side of the plane (both x and y coordinates are greater than zero). \$\endgroup\$
    – padawan
    May 30, 2014 at 12:06
  • \$\begingroup\$ I think those errors are not because of the method. \$\endgroup\$
    – padawan
    May 30, 2014 at 12:13
  • \$\begingroup\$ But again: what are those black lines? \$\endgroup\$ May 30, 2014 at 12:20
  • \$\begingroup\$ Those black lines are where the points should be. Blue points are where they are calculated. \$\endgroup\$
    – padawan
    May 30, 2014 at 12:28
  • \$\begingroup\$ And how do you know where they should be? You're trying to find their position? \$\endgroup\$ May 30, 2014 at 12:30
4
\$\begingroup\$

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));
}
\$\endgroup\$
15
  • \$\begingroup\$ a, b, c and d are the components of target plane, right? \$\endgroup\$
    – padawan
    May 29, 2014 at 11:46
  • \$\begingroup\$ Yes, they’re from your question: “I have another set of points double[] m on the plane ax + by + cz + d = 0”. \$\endgroup\$ May 29, 2014 at 11:52
  • \$\begingroup\$ Look at his requirements. He says he has 3 points on the original plane, and 3 points on the target plane. The transformation matrix should match these 3 points. By only fixing the plane, there is the degree of freedom I talked about in my old answer (still online) and the comments. \$\endgroup\$ May 29, 2014 at 11:53
  • \$\begingroup\$ @MartijnCourteaux But the question doesn’t mention these three points, and they appear nowhere in the code… The question is “I want to turn a plane into another plane, preserving the pairwise distances”. \$\endgroup\$ May 29, 2014 at 12:30
  • \$\begingroup\$ Yes, I know. It's a bit vague. But I explained it first by just rotating it correctly, but he wanted to get rid of the degree of freedom. \$\endgroup\$ May 29, 2014 at 12:34
3
+50
\$\begingroup\$

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:

Source and Dest points

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

Translate Matrices

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:

Centered 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:

Solve lin comb

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:

Image of X unit

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.

Image of Z unit

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:

Images and normalize

(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:

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:

Total matrix


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.

enter image description 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. enter image description 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.

\$\endgroup\$
17
  • \$\begingroup\$ May the force be with you in your exams, my noble friend. But there is no randomness. Why would you suggest that? \$\endgroup\$
    – padawan
    May 25, 2014 at 15:42
  • \$\begingroup\$ You can rotate all the points in the target plane around it's normal vector, and they will still be in the same plane. That is your degree of freedom. \$\endgroup\$ May 25, 2014 at 15:46
  • \$\begingroup\$ What if the orientation is the same? I mean, when you rotate the initial plane, it is parallel to the target plane? \$\endgroup\$
    – padawan
    May 25, 2014 at 15:48
  • \$\begingroup\$ No of course they are not parallel anymore. Parallel planes have the same normal vector. So a parallel plane to z=0 would be z=d. You really need to find a way to specify how you are going get rid of your degree of freedom. Think about it. If you find something, I will be glad to do the calculations in Maple for you :) \$\endgroup\$ May 25, 2014 at 15:51
  • \$\begingroup\$ Yes, I see what you are trying. But this preservers distance. Every rotation preserves distance. Another way to clarify your degree of freedom is to see that, before you rotate your plane into space, that you rotate all your points in z=0 around the z axis, so your points are still in z=0, but rotate around the origin. \$\endgroup\$ May 25, 2014 at 16:14

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .