==> include/pathtrace/t_sample.h <==1/* Template (included once per instantiation after t_math.h and t_intersect.h; see pathtrace/real.h):
2 * one camera sample's radiance for pixel (x, y), sample s, and a pixel's pass (its samples
3 * p S .. p S + S - 1 summed in order from +0). */
4#include "pathtrace/real.h"
5#include "pathtrace/rng.h"
6#include "pathtrace/scene.h"
7
8/* A 24-bit integer to [0, 1). */
9static inline float PT_FN(u01)(uint32_t h) { return PT_MUL(PT_FROM_U32(h >> 8), PT_K(PT_C_U24)); }
10
11/* The hit's geometric normal (unit) and material. */
12static inline pt_v PT_FN(surface)(const uint32_t *sc, const pt_layout *l, uint32_t id, pt_v h, uint32_t *mat) {
13 if (id >= PT_ID_TRI) {
14 size_t b = l->tri + PT_TRI_WORDS * (size_t)(id - PT_ID_TRI);
15 *mat = sc[b + 12];
16 return PT_FN(ldv)(sc, b + 9);
17 }
18 if (id >= PT_ID_SPH) {
19 size_t b = l->sph + PT_SPH_WORDS * (size_t)(id - PT_ID_SPH);
20 *mat = sc[b + 5];
21 return PT_FN(vscale)(PT_FN(vsub)(h, PT_FN(ldv)(sc, b)), PT_FN(ldr)(sc, b + 4));
22 }
23 size_t b = l->quad + PT_QUAD_WORDS * (size_t)(id - 1);
24 *mat = sc[b + 15];
25 return PT_FN(ldv)(sc, b + 9);
26}
27
28/* An orthonormal basis around unit n (Duff et al. 2017), applied to (x, y, z). */
29static inline pt_v PT_FN(onb)(pt_v n, float x, float y, float z) {
30 float one = PT_K(PT_C_ONE);
31 float sign = (pt_to_bits(n.z) >> 31) != 0 ? PT_NEG(one) : one;
32 float a = PT_DIV(PT_NEG(one), PT_ADD(sign, n.z));
33 float b = PT_MUL(PT_MUL(n.x, n.y), a);
34 pt_v t = PT_FN(v3)(PT_ADD(one, PT_MUL(PT_MUL(PT_MUL(sign, n.x), n.x), a)), PT_MUL(sign, b), PT_MUL(PT_NEG(sign), n.x)); /* Rust's -sign * n.x is (-sign) * n.x */
35 pt_v s = PT_FN(v3)(b, PT_ADD(sign, PT_MUL(PT_MUL(n.y, n.y), a)), PT_NEG(n.y));
36 return PT_FN(v3)(PT_ADD(PT_ADD(PT_MUL(t.x, x), PT_MUL(s.x, y)), PT_MUL(n.x, z)),
37 PT_ADD(PT_ADD(PT_MUL(t.y, x), PT_MUL(s.y, y)), PT_MUL(n.y, z)),
38 PT_ADD(PT_ADD(PT_MUL(t.z, x), PT_MUL(s.z, y)), PT_MUL(n.z, z)));
39}
40
41static inline pt_v PT_FN(reflect)(pt_v d, pt_v n) {
42 return PT_FN(vsub)(d, PT_FN(vscale)(n, PT_MUL(PT_K(PT_C_TWO), PT_FN(dot)(d, n))));
43}
44
45/* A NaN sample counts 0 (a bit test), and a large one `cl` (fireflies). */
46static inline float PT_FN(clamp_sample)(float x, float cl) {
47 if ((pt_to_bits(x) & 0x7fffffffu) > 0x7f800000u) return PT_K(PT_C_ZERO);
48 return PT_GT(x, cl) ? cl : x;
49}
50
51/* One camera sample's radiance: pixel (x, y), sample s. */
52PT_API pt_v PT_FN(sample)(const uint32_t *sc, const pt_layout *l, uint32_t x, uint32_t y, uint32_t s) {
53 float zero = PT_K(PT_C_ZERO), one = PT_K(PT_C_ONE);
54 uint32_t w = sc[PT_H_WIDTH], hgt = sc[PT_H_HEIGHT];
55 uint32_t key = pt_sample_key(sc[PT_H_SEED], y * w + x, s);
56 float fx = PT_DIV(PT_ADD(PT_FROM_U32(x), PT_FN(u01)(pt_rnd(key, 0))), PT_FROM_U32(w));
57 float fy = PT_DIV(PT_ADD(PT_FROM_U32(y), PT_FN(u01)(pt_rnd(key, 1))), PT_FROM_U32(hgt));
58 pt_v eye = PT_FN(ldv)(sc, PT_H_EYE);
59 pt_v dir = PT_FN(vsub)(PT_FN(vadd)(PT_FN(vadd)(PT_FN(ldv)(sc, PT_H_CORNER), PT_FN(vscale)(PT_FN(ldv)(sc, PT_H_HORIZ), fx)),
60 PT_FN(vscale)(PT_FN(ldv)(sc, PT_H_VERT), fy)),
61 eye);
62 pt_v d = PT_FN(normalize)(dir);
63 pt_v o = eye;
64 pt_v tp = PT_FN(v3)(one, one, one);
65 pt_v acc = PT_FN(zero3)();
66 /* Emission counts on camera rays and after specular bounces; after a diffuse bounce the NEE
67 * light's was already counted by its light sample. */
68 int spec = 1;
69 uint32_t depth = sc[PT_H_DEPTH];
70 for (uint32_t b = 0; b < depth; b++) {
71 uint32_t id;
72 float t = PT_FN(intersect)(sc, l, o, d, PT_K(PT_C_BIG), &id);
73 if (id == 0) break;
74 pt_v h = PT_FN(vadd)(o, PT_FN(vscale)(d, t));
75 uint32_t mat;
76 pt_v n = PT_FN(surface)(sc, l, id, h, &mat);
77 int front = PT_LT(PT_FN(dot)(d, n), zero);
78 pt_v ns = front ? n : PT_FN(vneg)(n);
79 size_t mb = l->mat + PT_MAT_WORDS * (size_t)mat;
80 uint32_t kind = sc[mb];
81 pt_v albedo = PT_FN(ldv)(sc, mb + 1);
82 uint32_t dim = 2 + 8 * b;
83 if (kind == PT_EMIT) {
84 if (front && (spec || id != l->light + 1u)) acc = PT_FN(vadd)(acc, PT_FN(vmul)(tp, PT_FN(ldv)(sc, mb + 4)));
85 break;
86 }
87 if (kind == PT_DIFFUSE) {
88 if (l->light != PT_NO_LIGHT) {
89 size_t lb = l->quad + PT_QUAD_WORDS * (size_t)l->light;
90 pt_v lp = PT_FN(vadd)(PT_FN(vadd)(PT_FN(ldv)(sc, lb), PT_FN(vscale)(PT_FN(ldv)(sc, lb + 3), PT_FN(u01)(pt_rnd(key, dim)))),
91 PT_FN(vscale)(PT_FN(ldv)(sc, lb + 6), PT_FN(u01)(pt_rnd(key, dim + 1))));
92 pt_v wv = PT_FN(vsub)(lp, h);
93 float dist2 = PT_FN(dot)(wv, wv);
94 float dist = PT_FN(sqrt)(dist2);
95 pt_v wl = PT_FN(vscale)(wv, PT_DIV(one, dist));
96 float cs = PT_FN(dot)(ns, wl);
97 float cl = PT_NEG(PT_FN(dot)(PT_FN(ldv)(sc, lb + 9), wl));
98 if (PT_GT(cs, zero) && PT_GT(cl, zero)) {
99 uint32_t blocker;
100 PT_FN(intersect)(sc, l, PT_FN(vadd)(h, PT_FN(vscale)(ns, PT_K(PT_C_EPS))), wl,
101 PT_SUB(dist, PT_K(PT_C_SHADOW_EPS)), &blocker);
102 if (blocker == 0) {
103 pt_v le = PT_FN(ldv)(sc, l->mat + PT_MAT_WORDS * (size_t)sc[lb + 15] + 4);
104 float g = PT_DIV(PT_MUL(PT_MUL(PT_MUL(cs, cl), PT_FN(ldr)(sc, lb + 14)), PT_K(PT_C_INV_PI)), dist2);
105 acc = PT_FN(vadd)(acc, PT_FN(vscale)(PT_FN(vmul)(PT_FN(vmul)(tp, albedo), le), g));
106 }
107 }
108 }
109 /* A cosine-weighted bounce: the albedo is the whole weight. */
110 float u1 = PT_FN(u01)(pt_rnd(key, dim + 2));
111 float r = PT_FN(sqrt)(u1);
112 float c, sn;
113 PT_FN(cis)(pt_rnd(key, dim + 3), &c, &sn);
114 d = PT_FN(normalize)(PT_FN(onb)(ns, PT_MUL(r, c), PT_MUL(r, sn), PT_FN(sqrt)(PT_SUB(one, u1))));
115 o = PT_FN(vadd)(h, PT_FN(vscale)(ns, PT_K(PT_C_EPS)));
116 spec = 0;
117 } else if (kind == PT_METAL) {
118 float z = PT_SUB(one, PT_MUL(PT_K(PT_C_TWO), PT_FN(u01)(pt_rnd(key, dim + 2))));
119 float r = PT_FN(sqrt)(PT_SUB(one, PT_MUL(z, z)));
120 float c, sn;
121 PT_FN(cis)(pt_rnd(key, dim + 3), &c, &sn);
122 float fuzz = PT_FN(ldr)(sc, mb + 7);
123 d = PT_FN(normalize)(PT_FN(vadd)(PT_FN(reflect)(d, ns), PT_FN(vscale)(PT_FN(v3)(PT_MUL(r, c), PT_MUL(r, sn), z), fuzz)));
124 if (!PT_GT(PT_FN(dot)(d, ns), zero)) break;
125 o = PT_FN(vadd)(h, PT_FN(vscale)(ns, PT_K(PT_C_EPS)));
126 spec = 1;
127 } else {
128 float ior = PT_FN(ldr)(sc, mb + 7);
129 float ci = PT_NEG(PT_FN(dot)(d, ns));
130 float eta = front ? PT_DIV(one, ior) : ior;
131 float sin2 = PT_MUL(PT_MUL(eta, eta), PT_SUB(one, PT_MUL(ci, ci)));
132 float r0 = PT_DIV(PT_SUB(one, ior), PT_ADD(one, ior));
133 r0 = PT_MUL(r0, r0);
134 float m = PT_SUB(one, ci);
135 float fr = PT_ADD(r0, PT_MUL(PT_SUB(one, r0), PT_MUL(PT_MUL(PT_MUL(PT_MUL(m, m), m), m), m)));
136 if (PT_GT(sin2, one) || PT_LT(PT_FN(u01)(pt_rnd(key, dim + 4)), fr)) {
137 d = PT_FN(reflect)(d, ns);
138 o = PT_FN(vadd)(h, PT_FN(vscale)(ns, PT_K(PT_C_EPS)));
139 } else {
140 float ct = PT_FN(sqrt)(PT_SUB(one, sin2));
141 d = PT_FN(normalize)(PT_FN(vadd)(PT_FN(vscale)(d, eta), PT_FN(vscale)(ns, PT_SUB(PT_MUL(eta, ci), ct))));
142 o = PT_FN(vsub)(h, PT_FN(vscale)(ns, PT_K(PT_C_EPS)));
143 }
144 spec = 1;
145 }
146 tp = PT_FN(vmul)(tp, albedo);
147 if (b >= PT_RR_DEPTH) {
148 float p = tp.x;
149 if (PT_GT(tp.y, p)) p = tp.y;
150 if (PT_GT(tp.z, p)) p = tp.z;
151 if (PT_LT(p, PT_K(PT_C_RR_LO))) p = PT_K(PT_C_RR_LO);
152 if (PT_GT(p, PT_K(PT_C_RR_HI))) p = PT_K(PT_C_RR_HI);
153 if (!PT_LT(PT_FN(u01)(pt_rnd(key, dim + 5)), p)) break;
154 tp = PT_FN(v3)(PT_DIV(tp.x, p), PT_DIV(tp.y, p), PT_DIV(tp.z, p));
155 }
156 }
157 float cl = PT_FN(ldr)(sc, PT_H_CLAMP);
158 return PT_FN(v3)(PT_FN(clamp_sample)(acc.x, cl), PT_FN(clamp_sample)(acc.y, cl), PT_FN(clamp_sample)(acc.z, cl));
159}
160
161/* Pass p of pixel (x, y): its samples p S .. p S + S - 1 summed in order from +0. */
162PT_API pt_v PT_FN(pixel_pass)(const uint32_t *sc, const pt_layout *l, uint32_t x, uint32_t y, uint32_t p) {
163 uint32_t spp = sc[PT_H_SPP];
164 pt_v sum = PT_FN(zero3)();
165 for (uint32_t s = p * spp; s < (p + 1) * spp; s++) sum = PT_FN(vadd)(sum, PT_FN(sample)(sc, l, x, y, s));
166 return sum;
167}