Ray-plane intersection: inaccurate results - rounding errors?

Question

I've set up a demo with simple 3d first-person demo using C++ and OpenGL, and it seems to work reasonably well. My goal is this: when the user points the camera at a plane and clicks the left mouse button, I want to draw the intersection of a ray pointing in the direction the camera is facing from the player's position with that plane.

So, I start off with two Vectors, Vector position and Vector rotation, where Vector is a pretty standard three-dimensional vector class:

class Vector
{
  public:
    GLfloat x, y, z;

    Vector() {};

    Vector(GLfloat x, GLfloat y, GLfloat z)
    {
      this->x = x;
      this->y = y;
      this->z = z;
    }

    GLfloat dot(const Vector &vector) const
    {
      return x * vector.x + y * vector.y + z * vector.z;
    }

    ... etc ...

And Plane p, with Plane being a simple struct storing the normal of the plane and d. I copied this struct directly from the book "Real-Time Collision Detection," by Christer Ericson:

struct Plane 
{
  Vector n; // Plane normal. Points x on the plane satisfy Dot(n,x) = d
  float d; // d = dot(n,p) for a given point p on the plane
};

To start, I take position as the start of the ray, which I call a. I use that point and rotation to find the end of the ray, b. Then I use an algorithm for finding the intersection of a ray and a plane from that same book. I've actually implemented the same method myself, but I'm using the code from the book directly here just to make sure I didn't mess anything up:

void pickPoint()
{
    const float length = 100.0f;

    // Points a and b
    Vector a = State::position;
    Vector b = a;

    // Find point b of directed line ab
    Vector radians(Math::rad(State::rotation.x), Math::rad(State::rotation.y), 0);
    const float lengthYZ = Math::cos(radians.x) * length;

    b.y -= Math::sin(radians.x) * length;
    b.x += Math::sin(radians.y) * lengthYZ;
    b.z -= Math::cos(radians.y) * lengthYZ;

    // Compute the t value for the directed line ab intersecting the plane
    Vector ab = b - a;

    GLfloat t = (p.d - p.n.dot(a)) / p.n.dot(ab);

    printf("Plane normal: %f, %f, %f\n", p.n.x, p.n.y, p.n.z);
    printf("Plane value d: %f\n", p.d);
    printf("Rotation (degrees): %f, %f, %f\n", State::rotation.x, State::rotation.y, State::rotation.z);
    printf("Rotation (radians): %f, %f, %f\n", radians.x, radians.y, radians.z);
    printf("Point a: %f, %f, %f\n", a.x, a.y, a.z);
    printf("Point b: %f, %f, %f\n", b.x, b.y, b.z);
    printf("Expected length of ray: %f\n", length);
    printf("Actual length of ray: %f\n", ab.length());
    printf("Value t: %f\n", t);

    // If t in [0..1] compute and return intersection point
    if(t >= 0.0f && t <= 1.0f) 
    {
        point = a + t * ab;
        printf("Intersection: %f, %f, %f\n", point.x, point.y, point.z);
    }
    // Else no intersection
    else
    {
        printf("No intersection found\n");
    }

    printf("\n\n");
}

When I render this point with OpenGL, it looks to be pretty close to the where the intersection of the ray and the plane would be. But from printing out the actual values, I discovered that for specific positions and rotations, the intersection point can be off by up to 0.000004. Here's an example of where the intersection is inaccurate - I know the intersection point is NOT on the plane because its Y value should be 0, not 0.000002. I could also sub it back into the plane equation and get an inequality:

Plane normal: 0.000000, 1.000000, 0.000000
Plane value d: 0.000000
Rotation (degrees): 70.100044, 1.899823, 0.000000
Rotation (radians): 1.223477, 0.033158, 0.000000
Point a: 20.818802, 27.240383, 15.124892
Point b: 21.947229, -66.788452, -18.894285
Expected length of ray: 100.000000
Actual length of ray: 100.000000
Value t: 0.289702
Intersection: 21.145710, 0.000002, 5.269455

Now, I know floating-point numbers are just approximations of real numbers, so I'm guessing this inaccuracy is just the effect of floating-point rounding, though it's possible I made a mistake somewhere else in the code. I know the intersection is off only by an extremely small amount, but I still care about it because I'm planning to use these points to define vertices of a model or level by snapping them to an arbitrarily-oriented grid, so I actually want those points to be ON that grid, even if they're slightly inaccurate. This might be a misguided approach - I don't really know.

So my question is: is this inaccuracy just floating-point rounding at work, or did I make a mistake somewhere else?

If it is just floating-point rounding, is there any way to deal with it? I've tried rounding the values of the rotation and position vectors in various ways, which obviously results in a less accurate intersection point, but I still sometimes get intersections that aren't on the plane. I did read an answer to a similar question (Is this plane-ray intersection code correct?) that mentions keeping the dimensions large, but I'm not sure exactly what that means.

Sorry if this question has been asked before - I searched, but I didn't see anything that was quite what I'm having trouble with. Thanks!

Answer 1

Your math seems correct and this definitely looks like a rounding error. I have a strong feeling that it is this line:

GLfloat t = (p.d - p.n.dot(a)) / p.n.dot(ab);

That said I don't see any other method to compute t. You could maybe verify if you are losing precision by using "%.12f" (or more) in your printf statements. Another way to pinpoint the culprit is to try doing your t computation step by step and printing the results along the way to see if you are losing precision somewhere.

Did you try using double precision floating point, if precision really matters to you that much?