Ray vs AABB (2D) reporting hits

Question

I've been working a 2D ray versus AABB implementation. I worked it out on paper, then read about a popular "branchless" algorithm. I decied to try the branchless one, but it gives me completely wrong results when one of the components is zero. I used the "naive" version and it's still reporting bad hits. (The difference is the check for a zero component vs simply allowing 1.0-over-0.0 to turn into +/- infinity)

Example: https://tavianator.com/2011/ray_box.html (adapted to C below):

#include <float.h>
#include <math.h>
#include <stdio.h>
typedef struct vec2_s { double x, y; } vec2_t;
typedef struct box_s { vec2_t min, max; } box_t;
typedef struct ray_s { vec2_t x0, n; } ray_t;

int intersection(box_t b, ray_t r) {
    double tmin = -INFINITY, tmax = INFINITY;

    if (r.n.x != 0.0) {
        double tx1 = (b.min.x - r.x0.x)/r.n.x;
        double tx2 = (b.max.x - r.x0.x)/r.n.x;

        tmin = fmax(tmin, fmin(tx1, tx2));
        tmax = fmin(tmax, fmax(tx1, tx2));
    }

    if (r.n.y != 0.0) {
        double ty1 = (b.min.y - r.x0.y)/r.n.y;
        double ty2 = (b.max.y - r.x0.y)/r.n.y;

        tmin = fmax(tmin, fmin(ty1, ty2));
        tmax = fmin(tmax, fmax(ty1, ty2));
    }

    return tmax >= tmin;
}

These don't work at all for me when one of direction components is zero. Consider the example below that used the method. This is an AABB from (-1,2) to (1,4). The ray is at (3,0) and goes up. Standard Cartesian coordinates.

This returns true that there is an intersection. Clearly there is not. It seems to be ignoring the X-axis completely.

int main() {
    box_t b;
    ray_t r;

    b.min.x = -1; b.max.x = 1;
    b.min.y =  2; b.max.y = 4;

    r.x0.x = 3;
    r.x0.y = 0;
    r.n.x = 0;
    r.n.y = 5;

    printf("Intersection? %s\n", intersection(b,r)? "yes" : "no");

}

I've worked it out on paper, and I came to the conclusion that when a component of the ray's direction vector is zero, then the test for that axis is reduced to a simple interval check. However, obviously, this requires quite a bit of "branching". Before I implement what I have on paper, can anyone explain to me what is going wrong here? I've read various implementations of the so-called "branchless" algorithm, and they all seem to do more or less the same thing and even mention that it correctly (?) handles cases when a component is zero -- but it's failing badly in my unit test. What am I missing?

Answer 1

The short answer is that the "branchy" code in the blog post is wrong, and my intuition is correct. The branchless does work correctly, as explained the blog post.