Ray/bounding box intersection test program

ray_box.c

#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>

typedef struct {
    double origin[3];
    double dir_inv[3];
} ray;

typedef struct {
    double min[3];
    double max[3];
} box;

static inline double min(double x, double y) {
#if SSE
    return x < y ? x : y;
#elif IEEE
    return y != y ? x : (x < y ? x : y);
#elif JAVA
    return x != x ? x : (x < y ? x : y);
#endif
}

static inline double max(double x, double y) {
#if SSE
    return x > y ? x : y;
#elif IEEE
    return y != y ? x : (x > y ? x : y);
#elif JAVA
    return x != x ? x : (x > y ? x : y);
#endif
}

bool intersection(box b, ray r) {
    double tmin = -INFINITY, tmax = INFINITY;

    for (int i = 0; i < 3; ++i) {
        double t1 = (b.min[i] - r.origin[i])*r.dir_inv[i];
        double t2 = (b.max[i] - r.origin[i])*r.dir_inv[i];

#if SSE
        tmin = max(tmin, min(min(t1, t2), tmax));
        tmax = min(tmax, max(max(t1, t2), tmin));
#else
        tmin = max(tmin, min(t1, t2));
        tmax = min(tmax, max(t1, t2));
#endif
    }

    return tmax > max(tmin, 0.0);
}

void check_reflections(ray r, bool expected, int* count, int* total) {
    box b = {{-1.0, -1.0, -1.0}, {1.0, 1.0, 1.0}};

    for (int i = 0; i < 1 << 3; ++i) {
        ray rr;
        for (int j = 0; j < 3; ++j) {
            if (i & (1 << j)) {
                rr.origin[j] = -r.origin[j];
                rr.dir_inv[j] = -r.dir_inv[j];
            } else {
                rr.origin[j] = r.origin[j];
                rr.dir_inv[j] = r.dir_inv[j];
            }
        }

        if (intersection(b, rr) != expected) {
            printf("(%g, %g, %g) + t*(%g, %g, %g):     \t%s\n",
                   rr.origin[0], rr.origin[1], rr.origin[2],
                   rr.dir_inv[0], rr.dir_inv[1], rr.dir_inv[2],
                   !expected ? "true" : "false");
        } else {
            ++*count;
        }

        ++*total;
    }
}

void check_permutations(ray r, bool expected, int* count, int* total) {
    static const int permutations[][3] = {
        {0, 1, 2},
        {0, 2, 1},
        {1, 0, 2},
        {1, 2, 0},
        {2, 0, 1},
        {2, 1, 0},
    };

    for (int i = 0; i < 6; ++i) {
        ray rr;
        for (int j = 0; j < 3; ++j) {
            rr.origin[j] = r.origin[permutations[i][j]];
            rr.dir_inv[j] = r.dir_inv[permutations[i][j]];
        }
        check_reflections(rr, expected, count, total);
    }
}

void check_scaled(double factor, bool expected, int* count, int* total) {
    ray r;

    r.origin[1] = factor;
    r.origin[2] = -2.0;
    r.dir_inv[2] = 1.0;

    for (int i = -1; i <= 1; ++i) {
        r.origin[0] = factor*i;

        for (int j = -1; j <= 1; j += 2) {
            r.dir_inv[1] = j*INFINITY;

            for (int k = 0; k <= 1 - abs(i); ++k) {
                r.dir_inv[0] = 8.0/k;
                check_permutations(r, expected, count, total);

                r.dir_inv[0] = -r.dir_inv[0];
                check_permutations(r, expected, count, total);
            }
        }
    }
}

int main() {
    int count = 0, total = 0;

    check_scaled(1.01, false, &count, &total);
    check_scaled(1.00, false, &count, &total);
    check_scaled(0.99, true, &count, &total);
    printf("%d/%d\n", count, total);

    return count == total ? EXIT_SUCCESS : EXIT_FAILURE;
}