Numeric problem with ray transformation

Question

I am performing ray triangle intersection by transforming the ray into the object space for efficiency. I am using the same technique from this question (The correct way to transform a ray with a matrix?).

My question has to do with the numeric stability of the transformed ray. This is my code

Matrix4f matModel= this->getWorldTransform();
Matrix4f matModelInverted = matModel.Inverted();
Matrix4f matInverseNormal = matModel.Transposed();

rayOrigin = matModelInverted.Transform(rayOrigin);
rayDirection = matInverseNormal.Transform(rayDirection).Normalized();

matModel is the Object-to-World transform and it is a very simple translation matrix with this value

| 1 0 0 -4 |
| 0 1 0 -15|
| 0 0 1 -13|
| 0 0 0  1 |

matInverseNormal is the transposed inverse of the matrix used to transform the ray origin and has the value

| 1  0   0  0|
| 0  1   0  0|
| 0  0   1  0|
|-4 -15 -13 1|

My code works as it should most of the time, the problem appears around certain value where the transformed ray direction behaves erratically and the ray flip flops by as much as 180 degrees for very small changes in the input.

For input ray direction {-0.59531, 0.64250, -0.48249} I get a transformed direction {-0.59531, 0.64250, -0.48249} which is visually the correct behavior. Perturbing the values of the input ray direction to {-0.59520 0.64320 -0.48169} (a change of about 10e-3 magnitude) gets me {0.59520, -0.64320, 0.48169}, which has flipped signs on all axes and is visually the incorrect behavior.

What am I doing wrong?

Answer 1

I managed to solve the problem but I would certainly appreciate a mathematical explanation as to why it works.

I performed the calculations using 4D Vectors instead of 3D. The fourth component for the ray origin was set to 1.0 and the fourth element in the direction component was set to 0.0. Then I did matrix multiplication as usual. And from the output values I used only the 1st three components.

The code now looks like this

Vector4f o = Vector4f(rayOrigin.x, rayOrigin.y, rayOrigin.z, 1.0f);
Vector4f d = Vector4f(rayDirection.x, rayDirection.y, rayDirection.z, 0.0f);

o = matModelInverted.Transform(o);
d = matInverseNormal.Transform(d).Normalized();

Works like a charm.