1/* Template (included once per instantiation after pathtrace/t_math.h; see pathtrace/real.h): the
2 * nearest hit of a ray in (EPS, tmax) against the scene's quads, spheres and triangles (the triangles
3 * behind their bounding sphere). Returns t, and the hit id in *id (0 none, 1 + i quad i, PT_ID_SPH + i
4 * sphere i, PT_ID_TRI + i triangle i). */
5#include "pathtrace/real.h"
6#include "pathtrace/scene.h"
7
8PT_API float PT_FN(intersect)(const uint32_t *sc, const pt_layout *l, pt_v o, pt_v d, float tmax, uint32_t *id_out) {
9    float tmin = PT_K(PT_C_EPS), one = PT_K(PT_C_ONE), zero = PT_K(PT_C_ZERO);
10    float best = tmax;
11    uint32_t id = 0;
12    for (uint32_t i = 0; i < l->nquad; i++) {
13        size_t b = l->quad + PT_QUAD_WORDS * (size_t)i;
14        pt_v n = PT_FN(ldv)(sc, b + 9);
15        float den = PT_FN(dot)(n, d);
16        if (PT_GT(den, PT_NEG(PT_K(PT_C_PAR_EPS))) && PT_LT(den, PT_K(PT_C_PAR_EPS))) continue;
17        pt_v p = PT_FN(ldv)(sc, b);
18        float t = PT_DIV(PT_FN(dot)(n, PT_FN(vsub)(p, o)), den);
19        if (!(PT_GT(t, tmin) && PT_LT(t, best))) continue;
20        pt_v w = PT_FN(vsub)(PT_FN(vadd)(o, PT_FN(vscale)(d, t)), p);
21        float a = PT_MUL(PT_FN(dot)(w, PT_FN(ldv)(sc, b + 3)), PT_FN(ldr)(sc, b + 12));
22        if (PT_LT(a, zero) || PT_GT(a, one)) continue;
23        float c = PT_MUL(PT_FN(dot)(w, PT_FN(ldv)(sc, b + 6)), PT_FN(ldr)(sc, b + 13));
24        if (PT_LT(c, zero) || PT_GT(c, one)) continue;
25        best = t;
26        id = 1 + i;
27    }
28    for (uint32_t i = 0; i < l->nsph; i++) {
29        size_t b = l->sph + PT_SPH_WORDS * (size_t)i;
30        float r = PT_FN(ldr)(sc, b + 3);
31        pt_v oc = PT_FN(vsub)(o, PT_FN(ldv)(sc, b));
32        float bq = PT_FN(dot)(oc, d);
33        float cq = PT_SUB(PT_FN(dot)(oc, oc), PT_MUL(r, r));
34        float disc = PT_SUB(PT_MUL(bq, bq), cq);
35        if (PT_LT(disc, zero)) continue;
36        float s = PT_FN(sqrt)(disc);
37        float t = PT_SUB(PT_NEG(bq), s);
38        if (!PT_GT(t, tmin)) t = PT_SUB(s, bq);
39        if (PT_GT(t, tmin) && PT_LT(t, best)) {
40            best = t;
41            id = PT_ID_SPH + i;
42        }
43    }
44    if (l->ntri > 0) {
45        /* The mesh's bounding sphere first. */
46        float mr = PT_FN(ldr)(sc, PT_H_MESH + 3);
47        pt_v oc = PT_FN(vsub)(o, PT_FN(ldv)(sc, PT_H_MESH));
48        float bq = PT_FN(dot)(oc, d);
49        float cq = PT_SUB(PT_FN(dot)(oc, oc), PT_MUL(mr, mr));
50        float disc = PT_SUB(PT_MUL(bq, bq), cq);
51        int near = !PT_LT(disc, zero);
52        if (near) {
53            float s = PT_FN(sqrt)(disc);
54            near = PT_GT(PT_SUB(s, bq), tmin) && PT_LT(PT_SUB(PT_NEG(bq), s), best);
55        }
56        if (near) {
57            for (uint32_t i = 0; i < l->ntri; i++) {
58                size_t b = l->tri + PT_TRI_WORDS * (size_t)i;
59                pt_v e1 = PT_FN(ldv)(sc, b + 3);
60                pt_v e2 = PT_FN(ldv)(sc, b + 6);
61                pt_v pv = PT_FN(cross)(d, e2);
62                float det = PT_FN(dot)(e1, pv);
63                if (PT_GT(det, PT_NEG(PT_K(PT_C_DET_EPS))) && PT_LT(det, PT_K(PT_C_DET_EPS))) continue;
64                float inv = PT_DIV(one, det);
65                pt_v tv = PT_FN(vsub)(o, PT_FN(ldv)(sc, b));
66                float u = PT_MUL(PT_FN(dot)(tv, pv), inv);
67                if (PT_LT(u, zero) || PT_GT(u, one)) continue;
68                pt_v qv = PT_FN(cross)(tv, e1);
69                float v = PT_MUL(PT_FN(dot)(d, qv), inv);
70                if (PT_LT(v, zero) || PT_GT(PT_ADD(u, v), one)) continue;
71                float t = PT_MUL(PT_FN(dot)(e2, qv), inv);
72                if (PT_GT(t, tmin) && PT_LT(t, best)) {
73                    best = t;
74                    id = PT_ID_TRI + i;
75                }
76            }
77        }
78    }
79    *id_out = id;
80    return best;
81}