Ray/AABB Intersection Incorrect

Question

I've tried to reimplement the Fast Graphics Gems Ray/AABB Intersection Method in C#:

// Based on "Fast Ray-Box Intersection" algorithm by Andrew Woo, "Graphics Gems", Academic Press, 1990
public unsafe Vector? IntersectionWith(Cuboid other) {
    const int NUM_DIMENSIONS = 3;
    Assure.Equal(NUM_DIMENSIONS, 3); // If that value is ever changed, this algorithm will need some maintenance

    const byte QUADRANT_MIN = 0;
    const byte QUADRANT_MAX = 1;
    const byte QUADRANT_BETWEEN = 2;

    // Step 1: Work out which direction from the start point to test for intersection for all 3 dimensions, and the distance
    byte* quadrants = stackalloc byte[NUM_DIMENSIONS];
    float* candidatePlanes = stackalloc float[NUM_DIMENSIONS];
    float* cuboidMinPoints = stackalloc float[NUM_DIMENSIONS];
    float* cuboidMaxPoints = stackalloc float[NUM_DIMENSIONS];
    float maxDistance = Single.NegativeInfinity;
    byte maxDistanceDimension = 0;

    bool startPointIsInsideCuboid = true;

    cuboidMinPoints[0] = other.X;
    cuboidMinPoints[1] = other.Y;
    cuboidMinPoints[2] = other.Z;
    cuboidMaxPoints[0] = other.X + other.Width;
    cuboidMaxPoints[1] = other.Y + other.Height;
    cuboidMaxPoints[2] = other.Z + other.Depth;

    for (byte i = 0; i < NUM_DIMENSIONS; ++i) {
        if (StartPoint[i] < cuboidMinPoints[i]) {
            quadrants[i] = QUADRANT_MIN;
            candidatePlanes[i] = cuboidMinPoints[i];
            startPointIsInsideCuboid = false;
        }
        else if (StartPoint[i] > cuboidMaxPoints[i]) {
            quadrants[i] = QUADRANT_MAX;
            candidatePlanes[i] = cuboidMaxPoints[i];
            startPointIsInsideCuboid = false;
        }
        else {
            quadrants[i] = QUADRANT_BETWEEN;
        }
    }

    if (startPointIsInsideCuboid) return StartPoint;

    // Step 2: Find farthest dimension from cuboid
    for (byte i = 0; i < NUM_DIMENSIONS; ++i) {
        // ReSharper disable once CompareOfFloatsByEqualityOperator Exact check is desired here: Anything other than 0f is usable
        if (quadrants[i] != QUADRANT_BETWEEN && Orientation[i] != 0f) {
            float thisDimensionDist = (candidatePlanes[i] - StartPoint[i]) / Orientation[i];
            if (thisDimensionDist > maxDistance) {
                maxDistance = thisDimensionDist;
                maxDistanceDimension = i;
            }
        }
    }

    if (maxDistance < 0f) return null;

    if (maxDistance - Length > MathUtils.FlopsErrorMargin) return null;

    float* intersectionPoint = stackalloc float[NUM_DIMENSIONS];

    for (byte i = 0; i < NUM_DIMENSIONS; ++i) {
        if (maxDistanceDimension == i) {
            intersectionPoint[i] = StartPoint[i] + maxDistance * Orientation[i];
            if (cuboidMinPoints[i] - intersectionPoint[i] > MathUtils.FlopsErrorMargin || intersectionPoint[i] - cuboidMaxPoints[i] > MathUtils.FlopsErrorMargin) return null;
        }
        else intersectionPoint[i] = candidatePlanes[i];

    }

    Vector result = new Vector(intersectionPoint[0], intersectionPoint[1], intersectionPoint[2]);
    if (!IsInfiniteLength && Vector.DistanceSquared(StartPoint, result) > Length * Length) return null;
    else return result;
}

However, although it sort of works, I'm getting incorrect results on the following part of a unit test:

Cuboid cuboid = new Cuboid(frontBottomLeft: new Vector(0f, 7.1f, 0f), width: 0f, height: 5f, depth: 0f);
Ray testRayC = new Ray(startPoint: new Vector(30f, 30f, 30f), orientation: new Vector(-1f, -1f, -1f));

Assert.AreEqual(
    null, 
    testRayC.IntersectionWith(cuboid)
);

I am expecting null from the call to testRayC.IntersectionWith(cuboid), but instead it returns a Vector(0, 12.1, 0), which is not a point on the ray at all.

---

So is it just a case of adding a final check that the calculated point is on the ray? Or (and this is what I suspect), have I made an error in transcribing the code? I have double and triple checked but didn't see anything obvious...

Answer 1

The problem in your code is when you do if (maxDistanceDimension == i) {. The original code checks if (whichPlane != i) {. I don't have your data structures, but a fix should look like:

        for (byte i = 0; i < NUM_DIMENSIONS; ++i)
        {
            if (maxDistanceDimension != i)
            {
                intersectionPoint[i] = StartPoint[i] + maxDistance * Orientation[i];
                if (intersectionPoint[i] < cuboidMinPoints[i] - MathUtils.FlopsErrorMargin || intersectionPoint[i] > cuboidMaxPoints[i] + MathUtils.FlopsErrorMargin)
                    return null;
            }
            else
            {
                intersectionPoint[i] = candidatePlanes[i];
            }
        }

Next, the following isn't in the original code. What is this for?

        if (maxDistance - Length > MathUtils.FlopsErrorMargin)
            return null;

If you are trying to check if the hit is within the extent of the ray, this may be a bug. Given that your Orientation does not appear to be normalized, maxDistance is not necessarily in units of length. This may not matter in the original algorithm, but if you are going to check maxDistance against some other length you need to normalize Orientation (make it dimensionless) so that

                thisDimensionDist = (candidatePlanes[i] - StartPoint[i]) / Orientation[i];

will have units of length.

Incidentally, in the original I think the following is wrong:

if(inside)  {
    coord = origin;
    return (TRUE);
}

Assuming this code is c and not c++, this simply sets the the coord pointer to have the same reference as the origin pointer, which will have no effect on the caller. This issue doesn't apply to your version, however.

Also, in the course of looking at this, I made a more literal c# transcription of the algorithm here:

public static class RayXCuboid
{
    enum HitQuadrant
    {
        Right = 0,
        Left = 1,
        Middle = 2,
    }

    const int Dimension = 3;

    [Conditional("DEBUG")]
    static void AssertValidArguments<TDoubleList>(params TDoubleList[] args) where TDoubleList : IList<double>
    {
        Debug.Assert(Dimension == 3);
        foreach (var list in args)
            Debug.Assert(list != null && list.Count == Dimension);
    }

    public static bool HitBoundingBox<TDoubleList>(TDoubleList minB, TDoubleList maxB, TDoubleList origin, TDoubleList dir, TDoubleList coord) where TDoubleList : IList<double>
    {
        AssertValidArguments(minB, maxB, origin, dir, coord);

        HitQuadrant[] quadrant = new HitQuadrant[Dimension];
        double[] maxT = new double[Dimension];
        double[] candidatePlane = new double[Dimension];

        /* Find candidate planes; this loop can be avoided if
        rays cast all from the eye(assume perpsective view) */
        bool inside = true;
        for (int i = 0; i < Dimension; i++)
            if (origin[i] < minB[i])
            {
                quadrant[i] = HitQuadrant.Left;
                candidatePlane[i] = minB[i];
                inside = false;
            }
            else if (origin[i] > maxB[i])
            {
                quadrant[i] = HitQuadrant.Right;
                candidatePlane[i] = maxB[i];
                inside = false;
            }
            else
            {
                quadrant[i] = HitQuadrant.Middle;
            }

        /* Ray origin inside bounding box */
        if (inside)
        {
            CopyTo(origin, coord);
            return true;
        }

        /* Calculate T distances to candidate planes */
        for (int i = 0; i < Dimension; i++)
            if (quadrant[i] != HitQuadrant.Middle && dir[i] != 0.0)
                maxT[i] = (candidatePlane[i] - origin[i]) / dir[i];
            else
                maxT[i] = -1.0;

        /* Get largest of the maxT's for final choice of intersection */
        int whichPlane = 0;
        for (int i = 1; i < Dimension; i++)
            if (maxT[whichPlane] < maxT[i])
                whichPlane = i;

        /* Check final candidate actually inside box */
        if (maxT[whichPlane] < 0.0)
        {
            FillWithDefault(coord);
            return false;
        }

        for (int i = 0; i < Dimension; i++)
            if (whichPlane != i)
            {
                coord[i] = origin[i] + maxT[whichPlane] * dir[i];
                if (coord[i] < minB[i] || coord[i] > maxB[i])
                {
                    FillWithDefault(coord);
                    return false;
                }
            }
            else
            {
                coord[i] = candidatePlane[i];
            }
        return true;                /* ray hits box */
    }

    static void FillWithDefault<T>(IList<T> list)
    {
        for (int i = 0; i < list.Count; i++)
            list[i] = default(T);
    }

    static void CopyTo<T>(IList<T> from, IList<T> to)
    {
        int arrayIndex = 0;
        foreach (var item in from)
            to[arrayIndex++] = item;
    }
}