Rotating arbitrary plane to be z-axis-aligned

Question

Correct me if I'm wrong, but...

Given the normal of an arbitrary source plane, and the normal of the plane after applying a desired rotation:

Vector3F sourceNormal = (x, y, z).normalize()
Vector3F desiredNormal = (0, 0, 1).normalize()

1) We can find the "axis of rotation" through the cross-product of the two normals

Vector3F rotationAxis = Vector3F.cross(sourceNormal, desiredNormal).normalize()

2) We can find the "angle of rotation" through the arc-cosine of the dot product of the two normals.

// Thanks nico - it was there in my project source, but it was omitted here.
float theta = Math.acos(Vector3F.dot(sourceNormal, desiredNormal))

3) We can apply the rotation to a set of points in order to orient the source plane to our desired plane.

float[] rotationMatrix = new float[16];

// X component
rotationMatrix[5] = rotationMatrix[10] = (float)Math.cos(theta);
rotationMatrix[9] = (float)Math.sin(theta);
rotationMatrix[6] = -rotationMatrix[9];

// Y component
rotationMatrix[0] = rotationMatrix[10] = (float)Math.cos(theta);
rotationMatrix[2] = (float)Math.sin(theta);
rotationMatrix[8] = -rotationMatrix[2];

// Z component
rotationMatrix[0] = rotationMatrix[5] = (float)Math.cos(theta);
rotationMatrix[1] = (float)Math.sin(theta);
rotationMatrix[4] = -rotationMatrix[1];

for(Point3F point : polygon)
{
    float x = pt.getX();
    float y = pt.getY();
    float z = pt.getZ();
    float[] xs = new float[3];
    float[] ys = new float[3];
    float[] zs = new float[3];
    for(int j = 0; j < 3; ++j)
    {
    xs[j] = rotationMatrix[j] * x;
    ys[j] = rotationMatrix[j + 4] * y;
    zs[j] = rotationMatrix[j + 8] * z;
    }
    x = 0; y = 0; z = 0;
    for(int j = 0; j < 3; ++j)
    {
    x += xs[j];
    y += ys[j];
    z += zs[j];
    }
    pt.set(x, y, z);
}

My output is incorrect.

In Points:

(-56.00, 72.01, 48.02)
(-48.00, 72.01, 48.02)
(-48.00, 86.01, 24.02)
(-56.00, 86.01, 24.02)

Out Points:

(-124.960010, -88.105451, 24.185812)
(-107.108590, -88.105451, 24.185812)
(-107.108590, -105.237051, 12.0929052)
(-124.960010, -105.237051, 12.0929052)

If I had to guess, I'd say that I am applying the rotations to the points incorrectly...perhaps I've interpreted the rotation matrix found in this article ( http://en.wikipedia.org/wiki/Rotation_matrix ) incorrectly?

Thanks for any input.

...assuming that this is the correct way to set up the rotation matrix, the output is still incorrect:

Vector3F axis = Vector3F.cross(sourceNormal, desiredNormal).normalize();
float angle = (float) Math.acos(p.normal.dot(new Vector3F(0, 0, 1)));
float s = (float)Math.sin(angle);
float c = (float)Math.cos(angle);
float x = axis.getX(), y = axis.getY(), z = axis.getZ();
float[] matrix = new float[16];
matrix[0] = x * x * (1 - c) + c;
matrix[1] = x * y * (1 - c) - (z * s);
matrix[2] = x * z * (1 - c) + (y * s);

matrix[4] = y * x * (1 - c) + (z * s);
matrix[5] = y * y * (1 - c) + c;
matrix[6] = y * z * (1 - c) - (x * s);

matrix[8] = x * z * (1 - c) - (y * s);
matrix[9] = y * z * (1 - c) + (x * s);
matrix[10] = z * z * (1 - c) + c;


float nx = x * matrix[0] + y * matrix[1] + z * matrix[2];
float ny = x * matrix[4] + y * matrix[5] + z * matrix[6];
float nz = x * matrix[8] + y * matrix[9] + z * matrix[10];

In: (-56.00, 56.01, -16.02)
In: (-48.00, 56.01, -16.02)
In: (-48.00, 72.01, -8.02)
In: (-56.00, 72.01, -8.02)
Out: (-51.270340, 46.887921, -25.5108342)
Out: (-43.9460070, 46.887921, -25.5108342)
Out: (-43.9460070, 62.798761, -21.5554182)
Out: (-51.270340, 62.798761, -21.5554182)

Answer 1

Your rotation matrix is wrong.

The rotations need to be combined per the formula in Section 10 of the Wikipedia article you've linked to ("Rotation matrix from axis and angle").

A simplified version of this matrix can also be found in the manual page for the glRotate function.

NB: unless you're repeating the rotation over and over, there's no need to take acos(dotp) only to take the cos and sin of it again. Use the dot product directly as cos(theta), and use the relationship sin(theta)^2 = 1 - cos(theta)^2 to work out sin(theta).

Answer 2

Ultimately, there were many issues. Alnitak solved the vast majority of them, and I've just found the last one.

In short, the following were wrong:

1) Matrix Creation

2) Matrix Multiplication

3) Floating Point Truncation

Even if it might be a hacky solution, it's fine for my purposes. The working snip is below.

Vector3F axis = new Vector3F(p.normal);
// Seems to work when adding precision...but NOT with (0, 0, 1)
Vector3F target = new Vector3F(0.000000000000000001f, 0.000000000000001f, 1.0f);
axis.cross(target);
axis.normalize();
float angle = (float) Math.acos(p.normal.dot(new Vector3F(0, 0, 1)));
float s = (float) Math.sin(angle);
float c = (float) Math.cos(angle);
float x = axis.getX(), y = axis.getY(), z = axis.getZ();
float[] matrix = new float[16];
matrix[0] = x * x * (1 - c) + c;
matrix[1] = x * y * (1 - c) - (z * s);
matrix[2] = x * z * (1 - c) + (y * s);
matrix[4] = y * x * (1 - c) + (z * s);
matrix[5] = y * y * (1 - c) + c;
matrix[6] = y * z * (1 - c) - (x * s);
matrix[8] = x * z * (1 - c) - (y * s);
matrix[9] = y * z * (1 - c) + (x * s);
matrix[10] = z * z * (1 - c) + c;

float nx = x * matrix[0] + y * matrix[1] + z * matrix[2];
float ny = x * matrix[4] + y * matrix[5] + z * matrix[6];
float nz = x * matrix[8] + y * matrix[9] + z * matrix[10];

// triangulate (now 2D) polygon
// Rotate the points back
angle = -angle;
s = (float) Math.sin(angle);
c = (float) Math.cos(angle);
matrix[0] = x * x * (1 - c) + c;
matrix[1] = x * y * (1 - c) - (z * s);
matrix[2] = x * z * (1 - c) + (y * s);
matrix[4] = y * x * (1 - c) + (z * s);
matrix[5] = y * y * (1 - c) + c;
matrix[6] = y * z * (1 - c) - (x * s);
matrix[8] = x * z * (1 - c) - (y * s);
matrix[9] = y * z * (1 - c) + (x * s);
matrix[10] = z * z * (1 - c) + c;
nx = x * matrix[0] + y * matrix[1] + z * matrix[2];
ny = x * matrix[4] + y * matrix[5] + z * matrix[6];
nz = x * matrix[8] + y * matrix[9] + z * matrix[10];
// All's well that ends well