1.3.77
 
Loading...
Searching...
No Matches
OptiX8DeviceCode.cu
Go to the documentation of this file.
1
16#include <optix.h>
17#include <optix_device.h>
18#include "OptiX8Math.h"
19#include "OptiX8LaunchParams.h"
20
21// ---------------------------------------------------------------------------
22// Launch params declared as constant memory (filled by optixLaunch)
23// ---------------------------------------------------------------------------
24extern "C" __constant__ OptiX8LaunchParams params;
25
26// ---------------------------------------------------------------------------
27// Device utility functions (adapted from RayTracing.cuh, OptiX-API-agnostic)
28// ---------------------------------------------------------------------------
29
30static __forceinline__ __device__ unsigned int lcg(unsigned int &prev) {
31 const unsigned int LCG_A = 1664525u;
32 const unsigned int LCG_C = 1013904223u;
33 prev = (LCG_A * prev + LCG_C);
34 return prev & 0x00FFFFFF;
35}
36
37static __forceinline__ __device__ float rnd(unsigned int &prev) {
38 return (float)lcg(prev) / (float)0x01000000;
39}
40
41template<unsigned int N>
42static __forceinline__ __device__ unsigned int tea(unsigned int val0, unsigned int val1) {
43 unsigned int v0 = val0, v1 = val1, s0 = 0;
44 for (unsigned int n = 0; n < N; n++) {
45 s0 += 0x9e3779b9;
46 v0 += ((v1 << 4) + 0xa341316c) ^ (v1 + s0) ^ ((v1 >> 5) + 0xc8013ea4);
47 v1 += ((v0 << 4) + 0xad90777d) ^ (v0 + s0) ^ ((v0 >> 5) + 0x7e95761e);
48 }
49 return v0;
50}
51
52__device__ __forceinline__ void atomicFloatAdd(float *address, float val) {
53 atomicAdd(address, val);
54}
55
56__device__ __forceinline__ void d_transformPoint(const float (&T)[16], float3 &v) {
57 float3 V;
58 V.x = T[0]*v.x + T[1]*v.y + T[2]*v.z + T[3];
59 V.y = T[4]*v.x + T[5]*v.y + T[6]*v.z + T[7];
60 V.z = T[8]*v.x + T[9]*v.y + T[10]*v.z + T[11];
61 v = V;
62}
63
64__device__ __forceinline__ float3 d_rotatePoint(const float3 &pos, float theta, float phi) {
65 float st = sinf(theta), ct = cosf(theta);
66 float sp = sinf(phi), cp = cosf(phi);
67 float3 tmp;
68 tmp.x = cp*ct*pos.x + (-sp)*pos.y + cp*st*pos.z;
69 tmp.y = sp*ct*pos.x + cp*pos.y + sp*st*pos.z;
70 tmp.z = -st*pos.x + ct*pos.z;
71 return tmp;
72}
73
74__device__ __forceinline__ float d_magnitude(const float3 v) {
75 return sqrtf(v.x*v.x + v.y*v.y + v.z*v.z);
76}
77
78static __forceinline__ __device__ float acos_safe(float x) {
79 return acosf(fmaxf(-1.f, fminf(1.f, x)));
80}
81
82static __forceinline__ __device__ float asin_safe(float x) {
83 return asinf(fmaxf(-1.f, fminf(1.f, x)));
84}
85
86// Translucent cover (glass/plastic) optical model: angular transmittance/reflectance/absorptance of a
87// single dielectric sheet via Fresnel reflection (two interfaces, both polarizations, internal
88// reflections summed in closed form) plus Bouguer absorption. This is the classic solar-glazing model
89// (Duffie & Beckman). Inputs: cos_theta = |cos| of the incidence angle, n = refractive index (>1),
90// KL = absorption product K*L (0 = lossless). Returns (tau, rho, alpha) with tau+rho+alpha == 1.
91static __forceinline__ __device__ float3 glass_tau_rho_alpha(float cos_theta, float n, float KL) {
92 cos_theta = fmaxf(1e-4f, fminf(1.f, cos_theta)); // guard grazing/degenerate
93 const float theta = acos_safe(cos_theta);
94 const float sin_t = sinf(theta);
95 const float sin_tr = sin_t / n; // Snell
96 const float cos_tr = sqrtf(fmaxf(0.f, 1.f - sin_tr * sin_tr));
97 const float theta_r = asin_safe(sin_tr);
98
99 // Fresnel interface reflectance for the two polarizations. As theta -> 0 the sin/tan ratios
100 // become 0/0, so fall back to the normal-incidence form r0 = ((n-1)/(n+1))^2.
101 float r_par, r_per;
102 if (theta < 1e-3f) {
103 const float r0 = ((n - 1.f) / (n + 1.f)) * ((n - 1.f) / (n + 1.f));
104 r_par = r0;
105 r_per = r0;
106 } else {
107 const float s_minus = sinf(theta_r - theta);
108 const float s_plus = sinf(theta_r + theta);
109 const float t_minus = tanf(theta_r - theta);
110 const float t_plus = tanf(theta_r + theta);
111 r_per = (s_minus * s_minus) / fmaxf(1e-12f, s_plus * s_plus); // perpendicular
112 r_par = (t_minus * t_minus) / fmaxf(1e-12f, t_plus * t_plus); // parallel
113 }
114
115 // Bouguer single-pass absorption through the angled path length L/cos_tr.
116 const float tau_a = (KL > 0.f) ? expf(-KL / fmaxf(1e-4f, cos_tr)) : 1.f;
117
118 // Per-polarization sheet transmittance/reflectance with internal reflections summed in closed form.
119 float tau = 0.f, rho = 0.f;
120 #pragma unroll
121 for (int pol = 0; pol < 2; pol++) {
122 const float r = (pol == 0) ? r_per : r_par;
123 const float denom = fmaxf(1e-6f, 1.f - (r * tau_a) * (r * tau_a));
124 const float tau_i = tau_a * (1.f - r) * (1.f - r) / denom;
125 const float rho_i = r * (1.f + tau_a * tau_i);
126 tau += 0.5f * tau_i;
127 rho += 0.5f * rho_i;
128 }
129 tau = fmaxf(0.f, fminf(1.f, tau));
130 rho = fmaxf(0.f, fminf(1.f, rho));
131 float alpha = 1.f - tau - rho;
132 if (alpha < 0.f) { alpha = 0.f; } // numerical guard
133 return make_float3(tau, rho, alpha);
134}
135
136// Initialize the per-band translucent-cover transmittance accumulator to 1 (no attenuation yet).
137static __forceinline__ __device__ void initCoverTransmittance(PerRayData &prd) {
138 #pragma unroll
139 for (int i = 0; i < HELIOS_MAX_RADIATION_BANDS; i++) {
140 prd.cover_transmittance[i] = 1.f;
141 }
142}
143
144// Evaluates the diffuse angular distribution (fd factor).
145// Priority 1: Power-law (Harrison & Coombes) if extinction > 0
146// Priority 2: Prague sky model if sky_params.w > 0
147// Priority 3: Isotropic (returns 1.0)
148static __device__ float evaluateDiffuseAngularDistribution(const float3 &ray_dir, const float3 &peak_dir,
149 float power_law_K, float power_law_norm,
150 const float4 &sky_params) {
151 if (power_law_K > 0.f) {
152 float psi = acos_safe(dot(peak_dir, ray_dir));
153 psi = fmaxf(psi, M_PI / 180.f);
154 return powf(psi, -power_law_K) * power_law_norm;
155 }
156 if (sky_params.w > 0.f) {
157 float gamma = acos_safe(dot(ray_dir, peak_dir)) * 180.f / M_PI;
158 float cos_theta = fmaxf(ray_dir.z, 0.f);
159 float pattern = (1.f + sky_params.x * expf(-gamma / sky_params.y))
160 * (1.f + (sky_params.z - 1.f) * (1.f - cos_theta));
161 return pattern * sky_params.w * M_PI;
162 }
163 return 1.f; // isotropic
164}
165
166// Load transform matrix for primitive at global position pos
167__device__ __forceinline__ void loadTransformMatrix(uint32_t pos, float (&T)[16]) {
168 for (int i = 0; i < 16; i++) {
169 T[i] = params.transform_matrix[pos * 16 + i];
170 }
171}
172
173// Sample texture mask at UV (uv_u, uv_v) in [0,1]^2 for mask at index msk_id.
174// Returns true if the texel is opaque (intersection should be reported),
175// false if transparent (reject the hit).
176// Standard texture convention: uv_v=0 maps to the top row (iy=0).
177static __forceinline__ __device__ bool sampleMask(int32_t msk_id, float uv_u, float uv_v) {
178 if (msk_id < 0) return true; // no mask → always opaque
179 const int32_t width = params.mask_sizes[msk_id * 2];
180 const int32_t height = params.mask_sizes[msk_id * 2 + 1];
181 const uint32_t offset = params.mask_offsets[msk_id];
182 int ix = (int)(floorf(float(width - 1) * uv_u));
183 int iy = (int)(floorf(float(height - 1) * (1.f - uv_v)));
184 ix = max(0, min(ix, width - 1));
185 iy = max(0, min(iy, height - 1));
186 return params.mask_data[offset + (uint32_t)(iy * width + ix)] != 0u;
187}
188
189// Build a rotation-only 4×4 transform matrix (row-major) from Euler angles (rx, ry, rz)
190// Matches OptiX 6's d_makeTransformMatrix convention.
191static __forceinline__ __device__ void d_makeTransformMatrix(float3 rotation, float (&T)[16]) {
192 float sx = sinf(rotation.x), cx = cosf(rotation.x);
193 float sy = sinf(rotation.y), cy = cosf(rotation.y);
194 float sz = sinf(rotation.z), cz = cosf(rotation.z);
195 T[0] = cz * cy; T[1] = cz * sy * sx - sz * cx; T[2] = cz * sy * cx + sz * sx; T[3] = 0.f;
196 T[4] = sz * cy; T[5] = sz * sy * sx + cz * cx; T[6] = sz * sy * cx - cz * sx; T[7] = 0.f;
197 T[8] = -sy; T[9] = cy * sx; T[10] = cy * cx; T[11] = 0.f;
198 T[12] = 0.f; T[13] = 0.f; T[14] = 0.f; T[15] = 1.f;
199}
200
201// Sample a point uniformly on the unit disk in the xy-plane (z=0).
202// Based on Suffern (2007) "Ray Tracing from the Ground Up", Ch. 6 concentric mapping.
203static __forceinline__ __device__ void d_sampleDisk(uint32_t &seed, float3 &sample) {
204 float Rx = rnd(seed), Ry = rnd(seed);
205 float sp_x = -1.f + 2.f * Rx;
206 float sp_y = -1.f + 2.f * Ry;
207 float r, p;
208 if (sp_x > -sp_y) {
209 if (sp_x > sp_y) { r = sp_x; p = sp_y / sp_x; }
210 else { r = sp_y; p = 2.f - sp_x / sp_y; }
211 } else {
212 if (sp_x < sp_y) { r = -sp_x; p = 4.f + sp_y / sp_x; }
213 else { r = -sp_y; p = (sp_y != 0.f) ? 6.f - sp_x / sp_y : 0.f; }
214 }
215 p *= 0.25f * M_PI;
216 sample = make_float3(r * cosf(p), r * sinf(p), 0.f);
217}
218
219// Sample a point uniformly on the unit square [-0.5,0.5]^2 in the xy-plane (z=0).
220static __forceinline__ __device__ void d_sampleSquare(uint32_t &seed, float3 &sample) {
221 sample = make_float3(-0.5f + rnd(seed), -0.5f + rnd(seed), 0.f);
222}
223
224// Invert a 4×4 row-major matrix (used for rect/disk source intersection tests).
225static __forceinline__ __device__ void d_invertMatrix(const float (&m)[16], float (&minv)[16]) {
226 float inv[16];
227 inv[0] = m[5]*m[10]*m[15] - m[5]*m[11]*m[14] - m[9]*m[6]*m[15] + m[9]*m[7]*m[14] + m[13]*m[6]*m[11] - m[13]*m[7]*m[10];
228 inv[4] = -m[4]*m[10]*m[15] + m[4]*m[11]*m[14] + m[8]*m[6]*m[15] - m[8]*m[7]*m[14] - m[12]*m[6]*m[11] + m[12]*m[7]*m[10];
229 inv[8] = m[4]*m[9]*m[15] - m[4]*m[11]*m[13] - m[8]*m[5]*m[15] + m[8]*m[7]*m[13] + m[12]*m[5]*m[11] - m[12]*m[7]*m[9];
230 inv[12] = -m[4]*m[9]*m[14] + m[4]*m[10]*m[13] + m[8]*m[5]*m[14] - m[8]*m[6]*m[13] - m[12]*m[5]*m[10] + m[12]*m[6]*m[9];
231 inv[1] = -m[1]*m[10]*m[15] + m[1]*m[11]*m[14] + m[9]*m[2]*m[15] - m[9]*m[3]*m[14] - m[13]*m[2]*m[11] + m[13]*m[3]*m[10];
232 inv[5] = m[0]*m[10]*m[15] - m[0]*m[11]*m[14] - m[8]*m[2]*m[15] + m[8]*m[3]*m[14] + m[12]*m[2]*m[11] - m[12]*m[3]*m[10];
233 inv[9] = -m[0]*m[9]*m[15] + m[0]*m[11]*m[13] + m[8]*m[1]*m[15] - m[8]*m[3]*m[13] - m[12]*m[1]*m[11] + m[12]*m[3]*m[9];
234 inv[13] = m[0]*m[9]*m[14] - m[0]*m[10]*m[13] - m[8]*m[1]*m[14] + m[8]*m[2]*m[13] + m[12]*m[1]*m[10] - m[12]*m[2]*m[9];
235 inv[2] = m[1]*m[6]*m[15] - m[1]*m[7]*m[14] - m[5]*m[2]*m[15] + m[5]*m[3]*m[14] + m[13]*m[2]*m[7] - m[13]*m[3]*m[6];
236 inv[6] = -m[0]*m[6]*m[15] + m[0]*m[7]*m[14] + m[4]*m[2]*m[15] - m[4]*m[3]*m[14] - m[12]*m[2]*m[7] + m[12]*m[3]*m[6];
237 inv[10] = m[0]*m[5]*m[15] - m[0]*m[7]*m[13] - m[4]*m[1]*m[15] + m[4]*m[3]*m[13] + m[12]*m[1]*m[7] - m[12]*m[3]*m[5];
238 inv[14] = -m[0]*m[5]*m[14] + m[0]*m[6]*m[13] + m[4]*m[1]*m[14] - m[4]*m[2]*m[13] - m[12]*m[1]*m[6] + m[12]*m[2]*m[5];
239 inv[3] = -m[1]*m[6]*m[11] + m[1]*m[7]*m[10] + m[5]*m[2]*m[11] - m[5]*m[3]*m[10] - m[9]*m[2]*m[7] + m[9]*m[3]*m[6];
240 inv[7] = m[0]*m[6]*m[11] - m[0]*m[7]*m[10] - m[4]*m[2]*m[11] + m[4]*m[3]*m[10] + m[8]*m[2]*m[7] - m[8]*m[3]*m[6];
241 inv[11] = -m[0]*m[5]*m[11] + m[0]*m[7]*m[9] + m[4]*m[1]*m[11] - m[4]*m[3]*m[9] - m[8]*m[1]*m[7] + m[8]*m[3]*m[5];
242 inv[15] = m[0]*m[5]*m[10] - m[0]*m[6]*m[9] - m[4]*m[1]*m[10] + m[4]*m[2]*m[9] + m[8]*m[1]*m[6] - m[8]*m[2]*m[5];
243 float det = m[0]*inv[0] + m[1]*inv[4] + m[2]*inv[8] + m[3]*inv[12];
244 det = 1.0f / det;
245 for (int i = 0; i < 16; i++) minv[i] = inv[i] * det;
246}
247
248// Test if ray hits a sphere source (any intersection in front of origin)
249static __forceinline__ __device__ bool d_raySphereIntersect(const float3 &ray_origin, const float3 &ray_direction,
250 const float3 &sphere_center, float sphere_radius) {
251 const float3 oc = make_float3(ray_origin.x - sphere_center.x, ray_origin.y - sphere_center.y,
252 ray_origin.z - sphere_center.z);
253 const float b = dot(oc, ray_direction);
254 const float c = dot(oc, oc) - sphere_radius * sphere_radius;
255 const float disc = b * b - c;
256 if (disc < 0.0f) return false;
257 return (-b - sqrtf(disc)) > 0.0f;
258}
259
260// Test if ray hits the front face of a rectangular source
261static __forceinline__ __device__ bool d_rayRectangleIntersect(const float3 &ray_origin, const float3 &ray_direction,
262 const float3 &rect_center, float rect_width, float rect_length,
263 const float3 &rect_rotation, float &out_cos_angle) {
264 float transform[16];
265 d_makeTransformMatrix(rect_rotation, transform);
266 const float3 normal = make_float3(transform[2], transform[6], transform[10]);
267 const float denom = dot(ray_direction, normal);
268 if (denom >= -1e-6f) return false;
269 const float3 oc = make_float3(rect_center.x - ray_origin.x, rect_center.y - ray_origin.y,
270 rect_center.z - ray_origin.z);
271 const float t = dot(oc, normal) / denom;
272 if (t <= 0.0f) return false;
273 float3 hit = make_float3(ray_origin.x + t * ray_direction.x - rect_center.x,
274 ray_origin.y + t * ray_direction.y - rect_center.y,
275 ray_origin.z + t * ray_direction.z - rect_center.z);
276 float inv_t[16];
277 d_invertMatrix(transform, inv_t);
278 d_transformPoint(inv_t, hit);
279 if (fabsf(hit.x) > rect_width * 0.5f || fabsf(hit.y) > rect_length * 0.5f) return false;
280 out_cos_angle = -denom;
281 return true;
282}
283
284// Test if ray hits the front face of a disk source
285static __forceinline__ __device__ bool d_rayDiskIntersect(const float3 &ray_origin, const float3 &ray_direction,
286 const float3 &disk_center, float disk_radius,
287 const float3 &disk_rotation, float &out_cos_angle) {
288 float transform[16];
289 d_makeTransformMatrix(disk_rotation, transform);
290 const float3 normal = make_float3(transform[2], transform[6], transform[10]);
291 const float denom = dot(ray_direction, normal);
292 if (denom >= -1e-6f) return false;
293 const float3 oc = make_float3(disk_center.x - ray_origin.x, disk_center.y - ray_origin.y,
294 disk_center.z - ray_origin.z);
295 const float t = dot(oc, normal) / denom;
296 if (t <= 0.0f) return false;
297 const float3 hit = make_float3(ray_origin.x + t * ray_direction.x - disk_center.x,
298 ray_origin.y + t * ray_direction.y - disk_center.y,
299 ray_origin.z + t * ray_direction.z - disk_center.z);
300 if (dot(hit, hit) > disk_radius * disk_radius) return false;
301 out_cos_angle = -denom;
302 return true;
303}
304
305// ---------------------------------------------------------------------------
306// PerRayData accessor (uses getPRD() from OptiX8LaunchParams.h)
307// ---------------------------------------------------------------------------
308
309// getPRD() is defined in OptiX8LaunchParams.h (guarded by #ifdef __CUDACC__)
310
311// ---------------------------------------------------------------------------
312// Intersection dispatch: handles all primitive types via primitive_type lookup.
313// All hitgroup SBT records reference this single entry point.
314// prim_idx = optixGetPrimitiveIndex() is the global AABB index (= global pos).
315// UUID is read from params.primitive_uuid[prim_idx] (global pos → UUID array).
316// Triangle and patch vertices are derived from the transform matrix using
317// canonical-space vertices that match those used in buildAABBs() on the host.
318// ---------------------------------------------------------------------------
319
320extern "C" __global__ void __intersection__patch() {
321 const uint32_t prim_idx = optixGetPrimitiveIndex();
322 const uint32_t pos = prim_idx; // global AABB index == global primitive position
323 const uint32_t ptype = params.primitive_type[pos];
324
325 const float3 ray_origin = optixGetWorldRayOrigin();
326 const float3 ray_direction = optixGetWorldRayDirection();
327 const float t_min = optixGetRayTmin();
328 const float t_max = optixGetRayTmax();
329
330 const uint32_t uuid = params.primitive_uuid[prim_idx];
331
332 if (ptype == 5) {
333 // ---- Bbox face: planar quad intersection (periodic boundary wall) ----
334 // Each bbox has 4 world-space vertices stored at bbox_vertices[bbox_local * 4 + v].
335 const uint32_t bbox_local = uuid - params.bbox_UUID_base;
336 const float3 v0 = params.bbox_vertices[bbox_local * 4 + 0];
337 const float3 v1 = params.bbox_vertices[bbox_local * 4 + 1];
338 const float3 v2 = params.bbox_vertices[bbox_local * 4 + 2];
339 const float3 v3 = params.bbox_vertices[bbox_local * 4 + 3];
340
341 // Quad normal from two edges
342 float3 e1 = make_float3(v1.x - v0.x, v1.y - v0.y, v1.z - v0.z);
343 float3 e3 = make_float3(v3.x - v0.x, v3.y - v0.y, v3.z - v0.z);
344 float3 n = cross(e1, e3);
345 float nd = dot(ray_direction, n);
346 if (fabsf(nd) < 1e-8f) return; // ray parallel to face
347
348 float t = dot(make_float3(v0.x - ray_origin.x, v0.y - ray_origin.y, v0.z - ray_origin.z), n) / nd;
349 if (t < t_min || t > t_max) return;
350
351 // Hit point must lie within the bounding box of the quad vertices
352 float3 hit = make_float3(ray_origin.x + t * ray_direction.x,
353 ray_origin.y + t * ray_direction.y,
354 ray_origin.z + t * ray_direction.z);
355 const float slack = 1e-4f;
356 float mnx = fminf(fminf(v0.x, v1.x), fminf(v2.x, v3.x)) - slack;
357 float mxx = fmaxf(fmaxf(v0.x, v1.x), fmaxf(v2.x, v3.x)) + slack;
358 float mny = fminf(fminf(v0.y, v1.y), fminf(v2.y, v3.y)) - slack;
359 float mxy = fmaxf(fmaxf(v0.y, v1.y), fmaxf(v2.y, v3.y)) + slack;
360 float mnz = fminf(fminf(v0.z, v1.z), fminf(v2.z, v3.z)) - slack;
361 float mxz = fmaxf(fmaxf(v0.z, v1.z), fmaxf(v2.z, v3.z)) + slack;
362 if (hit.x < mnx || hit.x > mxx || hit.y < mny || hit.y > mxy ||
363 hit.z < mnz || hit.z > mxz) return;
364
365 optixReportIntersection(t, 0, uuid, 0u);
366 return;
367 }
368
369 if (ptype != 0 && ptype != 1 && ptype != 3) return; // only patch, triangle, tile
370
371 float T[16];
372 loadTransformMatrix(pos, T);
373
374 if (ptype == 0 || ptype == 3) {
375 // ---- Patch: rectangle in canonical [-0.5, 0.5]^2 space ----
376 // Normal = third column of rotation part of T
377 float3 normal = make_float3(T[2], T[6], T[10]);
378 float nd = dot(ray_direction, normal);
379 if (fabsf(nd) < 1e-8f) return; // ray parallel to patch plane
380
381 // Patch centroid is the translation column of T
382 float3 patch_origin = make_float3(T[3], T[7], T[11]);
383 float t = dot(patch_origin - ray_origin, normal) / nd;
384 if (t < t_min || t > t_max) return;
385
386 // Project hit point into patch-local 2D coordinates
387 float3 hit_local = ray_origin + t * ray_direction - patch_origin;
388 float3 local_x = make_float3(T[0], T[4], T[8]);
389 float3 local_y = make_float3(T[1], T[5], T[9]);
390 float lx2 = dot(local_x, local_x);
391 float ly2 = dot(local_y, local_y);
392 if (lx2 < 1e-12f || ly2 < 1e-12f) return;
393 float u = dot(hit_local, local_x) / lx2;
394 float v = dot(hit_local, local_y) / ly2;
395 if (u < -0.5f || u > 0.5f || v < -0.5f || v > 0.5f) return;
396
397 // Texture mask check
398 const int32_t msk_id = params.mask_IDs[pos];
399 if (msk_id >= 0) {
400 float uv_u, uv_v;
401 if (params.uv_IDs[pos] >= 0) {
402 // Custom UV: bilinear from stored corner UVs
403 float2 uv0 = params.uv_data[pos * 4 + 0]; // UV at (u=0, v=0) corner
404 float2 uv1 = params.uv_data[pos * 4 + 1]; // UV at (u=1, v=0) corner
405 float2 uv2 = params.uv_data[pos * 4 + 2]; // UV at (u=0, v=1) corner
406 float du = uv1.x - uv0.x;
407 float dv = uv2.y - uv0.y;
408 uv_u = uv0.x + (u + 0.5f) * du;
409 uv_v = uv0.y + (v + 0.5f) * dv;
410 } else {
411 // Parametric UV: remap from [-0.5, 0.5] to [0, 1]
412 uv_u = u + 0.5f;
413 uv_v = v + 0.5f;
414 }
415 if (!sampleMask(msk_id, uv_u, uv_v)) return;
416 }
417
418 // Report the intersection for BOTH faces. A one-sided primitive (twosided_flag==0) still
419 // occludes a ray hitting its back face — matching the canonical OptiX 6 and Vulkan backends.
420 // The hit face is recorded in face_attr (1 = front/top, 0 = back/bottom) for downstream
421 // energy/face bookkeeping; one-sidedness is enforced at ray launch (raygen), not here.
422 uint32_t face_attr = (nd < 0.f) ? 1u : 0u;
423 optixReportIntersection(t, 0, uuid, face_attr);
424
425 } else {
426 // ---- Triangle: canonical vertices (0,0,0), (0,1,0), (1,1,0) ----
427 // World-space vertices via T (consistent with buildAABBs canonical vertices)
428 const float3 v0 = make_float3(T[3], T[7], T[11]);
429 const float3 v1 = make_float3(T[1] + T[3], T[5] + T[7], T[9] + T[11]);
430 const float3 v2 = make_float3(T[0] + T[1] + T[3], T[4] + T[5] + T[7], T[8] + T[9] + T[11]);
431
432 // Shirley's ray-triangle intersection (Möller–Trumbore style)
433 float a = v0.x - v1.x, b = v0.x - v2.x, c = ray_direction.x, d = v0.x - ray_origin.x;
434 float e = v0.y - v1.y, f = v0.y - v2.y, g = ray_direction.y, h = v0.y - ray_origin.y;
435 float i = v0.z - v1.z, j = v0.z - v2.z, k = ray_direction.z, l = v0.z - ray_origin.z;
436
437 float m = f * k - g * j, n = h * k - g * l, p = f * l - h * j;
438 float q = g * i - e * k, s = e * j - f * i;
439
440 float tri_denom = a * m + b * q + c * s;
441 if (fabsf(tri_denom) < 1e-12f) return;
442 float inv_denom = 1.f / tri_denom;
443
444 float e1 = d * m - b * n - c * p;
445 float beta = e1 * inv_denom;
446 if (beta < 0.f) return;
447
448 float r = e * l - h * i;
449 float e2 = a * n + d * q + c * r;
450 float gamma = e2 * inv_denom;
451 if (gamma < 0.f || beta + gamma > 1.f) return;
452
453 float e3 = a * p - b * r + d * s;
454 float t = e3 * inv_denom;
455 if (t < t_min || t > t_max) return;
456
457 // Texture mask check (using barycentric UV)
458 const int32_t msk_id = params.mask_IDs[pos];
459 if (msk_id >= 0) {
460 float uv_u, uv_v;
461 if (params.uv_IDs[pos] >= 0) {
462 // Custom UV: interpolate from stored per-vertex UV
463 float2 uv0 = params.uv_data[pos * 4 + 0];
464 float2 uv1 = params.uv_data[pos * 4 + 1];
465 float2 uv2 = params.uv_data[pos * 4 + 2];
466 // beta = weight at v1, gamma = weight at v2, (1-beta-gamma) = weight at v0
467 float2 uv = make_float2(uv0.x + beta * (uv1.x - uv0.x) + gamma * (uv2.x - uv0.x),
468 uv0.y + beta * (uv1.y - uv0.y) + gamma * (uv2.y - uv0.y));
469 uv_u = uv.x;
470 uv_v = 1.f - uv.y; // Y-flip to match OptiX 6 convention
471 } else {
472 // Parametric UV: use barycentric coordinates directly
473 uv_u = beta + gamma; // along u-axis (v1 is at u=1)
474 uv_v = gamma; // along v-axis (v2 is at v=1)
475 }
476 if (!sampleMask(msk_id, uv_u, uv_v)) return;
477 }
478
479 // Face from cross product normal vs ray direction
480 float3 edge0 = make_float3(v1.x - v0.x, v1.y - v0.y, v1.z - v0.z);
481 float3 edge1 = make_float3(v2.x - v0.x, v2.y - v0.y, v2.z - v0.z);
482 float3 tri_nrm = make_float3(edge0.y * edge1.z - edge0.z * edge1.y,
483 edge0.z * edge1.x - edge0.x * edge1.z,
484 edge0.x * edge1.y - edge0.y * edge1.x);
485 // Report the intersection for BOTH faces (see note in the patch branch above): one-sided
486 // primitives still occlude back-face hits, matching the canonical OptiX 6 / Vulkan backends.
487 uint32_t face_attr = (dot(ray_direction, tri_nrm) < 0.f) ? 1u : 0u;
488 optixReportIntersection(t, 0, uuid, face_attr);
489 }
490}
491
492extern "C" __global__ void __intersection__disk() {
493 // Disk intersection is not yet implemented. Disk primitives are rejected at
494 // the host level in updateGeometry(), so this program is never invoked.
495}
496
497extern "C" __global__ void __intersection__tile() {
498 // Tile intersection is handled by __intersection__patch (ptype == 3).
499}
500
501extern "C" __global__ void __intersection__voxel() {
502 // Voxel intersection is not yet implemented. Voxel primitives are rejected at
503 // the host level in updateGeometry(), so this program is never invoked.
504}
505
506extern "C" __global__ void __intersection__bbox() {
507 // Bbox intersection is not yet implemented in this backend.
508}
509
510// ---------------------------------------------------------------------------
511// Miss programs
512// ---------------------------------------------------------------------------
513
514extern "C" __global__ void __miss__direct() {
515 PerRayData *prd = getPayloadPRD();
516
517 const uint32_t origin_position = params.primitive_positions[prd->origin_UUID];
518 const uint32_t Nprims = params.Nprimitives;
519 const uint32_t Nbands_global = params.Nbands_global;
520 const uint32_t Nbands_launch = params.Nbands_launch;
521
522 int b = -1;
523 for (uint32_t b_global = 0; b_global < Nbands_global; b_global++) {
524 if (!params.band_launch_flag[b_global]) continue;
525 b++;
526
527 // radiation_in layout: [prim * Nbands_launch + band_launch]
528 const uint32_t ind_origin = origin_position * Nbands_launch + (uint32_t)b;
529
530 // rho/tau layout: [source * Nprims * Nbands_global + prim * Nbands_global + band_global]
531 const uint32_t radprop_ind = prd->source_ID * Nprims * Nbands_global
532 + origin_position * Nbands_global
533 + b_global;
534 const float t_rho = params.rho[radprop_ind];
535 const float t_tau = params.tau[radprop_ind];
536
537 // source_fluxes layout: [source * Nbands_launch + band_launch]
538 const uint32_t flux_idx = prd->source_ID * Nbands_launch + (uint32_t)b;
539 const float source_flux = params.source_fluxes[flux_idx];
540
541 // Attenuation from any translucent covers (glass/plastic) the ray passed through, per band.
542 const float cover_tau = (b < HELIOS_MAX_RADIATION_BANDS) ? prd->cover_transmittance[b] : 1.f;
543
544 const double strength = prd->strength * (double)source_flux * (double)cover_tau;
545 const float absorption = (float)(strength * (1.0 - t_rho - t_tau));
546
547 atomicFloatAdd(&params.radiation_in[ind_origin], absorption);
548
549 if (t_rho > 0.f || t_tau > 0.f) {
550 if (prd->face) {
551 atomicFloatAdd(&params.scatter_buff_top[ind_origin], (float)(strength * t_rho));
552 atomicFloatAdd(&params.scatter_buff_bottom[ind_origin], (float)(strength * t_tau));
553 } else {
554 atomicFloatAdd(&params.scatter_buff_bottom[ind_origin], (float)(strength * t_rho));
555 atomicFloatAdd(&params.scatter_buff_top[ind_origin], (float)(strength * t_tau));
556 }
557 }
558
559 // Camera-weighted scatter: mirrors scatter_buff but uses rho_cam/tau_cam
560 if (params.Ncameras > 0 && params.rho_cam && params.scatter_buff_top_cam) {
561 const uint32_t Ncameras = params.Ncameras;
562 const uint32_t cam_id = params.camera_ID;
563 const uint32_t rc_idx = prd->source_ID * Nprims * Nbands_global * Ncameras
564 + origin_position * Nbands_global * Ncameras
565 + b_global * Ncameras + cam_id;
566 const float t_rho_cam = params.rho_cam[rc_idx];
567 const float t_tau_cam = params.tau_cam ? params.tau_cam[rc_idx] : 0.f;
568 if ((t_rho_cam > 0.f || t_tau_cam > 0.f) && strength > 0.0) {
569 if (prd->face) {
570 atomicFloatAdd(&params.scatter_buff_top_cam[ind_origin], (float)(strength * t_rho_cam));
571 atomicFloatAdd(&params.scatter_buff_bottom_cam[ind_origin], (float)(strength * t_tau_cam));
572 } else {
573 atomicFloatAdd(&params.scatter_buff_bottom_cam[ind_origin], (float)(strength * t_rho_cam));
574 atomicFloatAdd(&params.scatter_buff_top_cam[ind_origin], (float)(strength * t_tau_cam));
575 }
576 }
577 }
578
579 // Accumulate incident radiation for specular for ALL cameras (per source, camera-weighted).
580 // Direct rays are launched once (not per camera), so we must populate every camera's
581 // slot in radiation_specular here so each camera's closest-hit can read its own data.
582 if (params.radiation_specular && params.source_fluxes_cam && strength > 0.0) {
583 for (uint32_t cam = 0; cam < params.Ncameras; cam++) {
584 // source_fluxes_cam layout: [source][band][camera] (full 3D buffer uploaded in updateSources)
585 const uint32_t weight_idx = prd->source_ID * Nbands_launch * params.Ncameras
586 + (uint32_t)b * params.Ncameras + cam;
587 const float camera_weight = params.source_fluxes_cam[weight_idx];
588 // radiation_specular layout: [source][camera][primitive][band]
589 const uint32_t ind_specular = prd->source_ID * params.Ncameras * Nprims * Nbands_launch
590 + cam * Nprims * Nbands_launch
591 + origin_position * Nbands_launch + (uint32_t)b;
592 atomicFloatAdd(&params.radiation_specular[ind_specular], (float)(strength * camera_weight));
593 }
594 }
595 }
596}
597
598extern "C" __global__ void __miss__diffuse() {
599 PerRayData *prd = getPayloadPRD();
600
601 const uint32_t origin_position = params.primitive_positions[prd->origin_UUID];
602 if (origin_position == UINT_MAX) return;
603
604 const uint32_t Nprims = params.Nprimitives;
605 const uint32_t Nbands_global = params.Nbands_global;
606 const uint32_t Nbands_launch = params.Nbands_launch;
607
608 if (params.diffuse_flux == nullptr) {
609 printf("ERROR (OptiX8 __miss__diffuse): diffuse_flux is null. "
610 "Call updateDiffuseRadiation() before launchDiffuseRays().\n");
611 __trap();
612 }
613
614 const float3 ray_dir = optixGetWorldRayDirection();
615
616 int b = -1;
617 for (uint32_t b_global = 0; b_global < Nbands_global; b_global++) {
618 if (!params.band_launch_flag[b_global]) continue;
619 b++;
620
621 if (params.diffuse_flux[b] <= 0.f) continue;
622
623 const float4 sky_p = params.sky_radiance_params ? params.sky_radiance_params[b] : make_float4(0.f, 0.f, 0.f, 0.f);
624 const float3 peak_d = params.diffuse_peak_dir ? params.diffuse_peak_dir[b] : make_float3(0.f, 0.f, 1.f);
625 const float power_K = params.diffuse_extinction ? params.diffuse_extinction[b] : 0.f;
626 const float power_N = params.diffuse_dist_norm ? params.diffuse_dist_norm[b] : 1.f;
627
628 const float fd = evaluateDiffuseAngularDistribution(ray_dir, peak_d, power_K, power_N, sky_p);
629 // Attenuation from any translucent covers (glass/plastic) the sky ray passed through, per band.
630 const float cover_tau = (b < HELIOS_MAX_RADIATION_BANDS) ? prd->cover_transmittance[b] : 1.f;
631 const float strength = fd * params.diffuse_flux[b] * (float)prd->strength * cover_tau;
632
633 const uint32_t ind_origin = origin_position * Nbands_launch + (uint32_t)b;
634 const uint32_t radprop_ind = prd->source_ID * Nprims * Nbands_global
635 + origin_position * Nbands_global + b_global;
636 const float t_rho = params.rho[radprop_ind];
637 const float t_tau = params.tau[radprop_ind];
638
639 atomicFloatAdd(&params.radiation_in[ind_origin], strength * (1.f - t_rho - t_tau));
640
641 if (t_rho > 0.f || t_tau > 0.f) {
642 if (prd->face) { // top-face origin
643 atomicFloatAdd(&params.scatter_buff_top[ind_origin], strength * t_rho);
644 atomicFloatAdd(&params.scatter_buff_bottom[ind_origin], strength * t_tau);
645 } else { // bottom-face origin
646 atomicFloatAdd(&params.scatter_buff_bottom[ind_origin], strength * t_rho);
647 atomicFloatAdd(&params.scatter_buff_top[ind_origin], strength * t_tau);
648 }
649 }
650
651 // Camera-weighted scatter: mirrors scatter_buff but uses rho_cam/tau_cam
652 if (params.Ncameras > 0 && params.rho_cam && params.scatter_buff_top_cam) {
653 const uint32_t Ncameras = params.Ncameras;
654 const uint32_t cam_id = params.camera_ID;
655 const uint32_t rc_idx = prd->source_ID * Nprims * Nbands_global * Ncameras
656 + origin_position * Nbands_global * Ncameras
657 + b_global * Ncameras + cam_id;
658 const float t_rho_cam = params.rho_cam[rc_idx];
659 const float t_tau_cam = params.tau_cam ? params.tau_cam[rc_idx] : 0.f;
660 if ((t_rho_cam > 0.f || t_tau_cam > 0.f) && strength > 0.f) {
661 if (prd->face) {
662 atomicFloatAdd(&params.scatter_buff_top_cam[ind_origin], strength * t_rho_cam);
663 atomicFloatAdd(&params.scatter_buff_bottom_cam[ind_origin], strength * t_tau_cam);
664 } else {
665 atomicFloatAdd(&params.scatter_buff_bottom_cam[ind_origin], strength * t_rho_cam);
666 atomicFloatAdd(&params.scatter_buff_top_cam[ind_origin], strength * t_tau_cam);
667 }
668 }
669 }
670 }
671}
672
673extern "C" __global__ void __miss__camera() {
674 PerRayData *prd = getPayloadPRD();
675 const uint32_t pixel_idx = prd->origin_UUID;
676 const uint32_t Nbands_l = params.Nbands_launch;
677 const float3 ray_origin = optixGetWorldRayOrigin();
678 const float3 ray_dir = optixGetWorldRayDirection();
679
680 for (uint32_t b = 0; b < Nbands_l; b++) {
681 float radiance = 0.0f;
682
683 for (uint32_t s = 0; s < params.Nsources; s++) {
684 const float flux = params.source_fluxes[s * Nbands_l + b];
685 if (flux <= 0.0f) continue;
686
687 const uint32_t stype = params.source_types[s];
688 if (stype == 0 || stype == 2) {
689 // Collimated / sun-sphere: treat as solar disk
690 if (params.solar_disk_radiance && params.solar_disk_radiance[b] > 0.0f &&
691 dot(ray_dir, params.sun_direction) >= params.solar_disk_cos_angle) {
692 radiance += params.solar_disk_radiance[b];
693 }
694 } else if (stype == 1) {
695 // Sphere
696 if (d_raySphereIntersect(ray_origin, ray_dir, params.source_positions[s],
697 params.source_widths[s].x * 0.5f)) {
698 const float area = 4.0f * M_PI * params.source_widths[s].x * 0.5f * params.source_widths[s].x * 0.5f;
699 radiance += (flux / area) / M_PI;
700 }
701 } else if (stype == 3) {
702 // Rectangle
703 float cos_angle;
704 if (d_rayRectangleIntersect(ray_origin, ray_dir, params.source_positions[s],
705 params.source_widths[s].x, params.source_widths[s].y,
706 params.source_rotations[s], cos_angle)) {
707 const float area = params.source_widths[s].x * params.source_widths[s].y;
708 radiance += (flux / area) * cos_angle / M_PI;
709 }
710 } else if (stype == 4) {
711 // Disk
712 float cos_angle;
713 if (d_rayDiskIntersect(ray_origin, ray_dir, params.source_positions[s],
714 params.source_widths[s].x, params.source_rotations[s], cos_angle)) {
715 const float area = M_PI * params.source_widths[s].x * params.source_widths[s].x;
716 radiance += (flux / area) * cos_angle / M_PI;
717 }
718 }
719 }
720
721 // Sky radiance fallback
722 if (radiance <= 0.0f && params.camera_sky_radiance && params.camera_sky_radiance[b] > 0.0f) {
723 const float4 sky_p = params.sky_radiance_params ? params.sky_radiance_params[b]
724 : make_float4(0.f, 0.f, 0.f, 0.f);
725 radiance = params.camera_sky_radiance[b] *
726 evaluateDiffuseAngularDistribution(ray_dir, params.sun_direction, 0.0f, 1.0f, sky_p);
727 }
728
729 // Isotropic sky emission/longwave: when band has emission enabled, the user-set diffuse_flux
730 // represents hemispherical sky thermal/longwave irradiance. Convert to isotropic radiance.
731 if (params.band_emission_flag && params.band_emission_flag[b] != 0u &&
732 params.camera_diffuse_flux && params.camera_diffuse_flux[b] > 0.0f) {
733 radiance += params.camera_diffuse_flux[b] / M_PI;
734 }
735
736 if (radiance > 0.0f) {
737 atomicFloatAdd(&params.radiation_in_camera[pixel_idx * Nbands_l + b],
738 radiance * (float)prd->strength);
739 }
740 }
741}
742
743extern "C" __global__ void __miss__pixel_label() {
744 PerRayData *prd = getPayloadPRD();
745 if (params.camera_pixel_depth) {
746 params.camera_pixel_depth[prd->origin_UUID] = -1.0f;
747 }
748}
749
750// ---------------------------------------------------------------------------
751// Closest-hit: direct radiation
752// ---------------------------------------------------------------------------
753
754// Shared helper: populate PerRayData with periodic boundary wrapping info.
755// Called from closesthit programs when a bbox (type-5) face is hit.
756// hit_uuid identifies which bbox face was hit (bbox_local = hit_uuid - bbox_UUID_base).
757// Bbox face ordering in RadiationModel.cpp:
758// x-only: 0=x-min, 1=x-max
759// y-only: 0=y-min, 1=y-max
760// xy: 0=x-min, 1=x-max, 2=y-min, 3=y-max
761// Vertex 0 of each face is always the "min" corner, so bbox_vertices[b*4+0] gives
762// the face's position (x or y coordinate) without floating-point tolerance issues.
763static __forceinline__ __device__ void handlePeriodicBoundaryHit(PerRayData *prd, uint32_t hit_uuid) {
764 const float t_hit = optixGetRayTmax();
765 const float3 ray_orig = optixGetWorldRayOrigin();
766 const float3 ray_dir = optixGetWorldRayDirection();
767 const float3 hit_pos = make_float3(ray_orig.x + t_hit * ray_dir.x,
768 ray_orig.y + t_hit * ray_dir.y,
769 ray_orig.z + t_hit * ray_dir.z);
770
771 const uint32_t bbox_local = hit_uuid - params.bbox_UUID_base;
772 prd->periodic_hit = hit_pos;
773
774 if (params.periodic_flag.x == 1 && bbox_local < 2) {
775 // x-faces: bbox 0 = x-min, bbox 1 = x-max
776 // vertex 0 of each face gives the face's x coordinate
777 const float xmin = params.bbox_vertices[0 * 4].x; // x-min face, any vertex .x = xmin
778 const float xmax = params.bbox_vertices[1 * 4].x; // x-max face, any vertex .x = xmax
779 const float width = xmax - xmin;
780 prd->periodic_hit.x += (bbox_local == 0) ? width : -width;
781 } else {
782 // y-faces: if x-periodic exists they are bboxes 2,3; otherwise 0,1
783 const uint32_t y_base = (params.periodic_flag.x == 1) ? 2u : 0u;
784 const float ymin = params.bbox_vertices[y_base * 4].y; // y-min face vertex 0 .y = ymin
785 const float ymax = params.bbox_vertices[(y_base + 1) * 4].y; // y-max face vertex 0 .y = ymax
786 const float width = ymax - ymin;
787 prd->periodic_hit.y += (bbox_local == y_base) ? width : -width;
788 }
789}
790
791// ---------------------------------------------------------------------------
792// Translucent cover (glass/plastic) support
793// ---------------------------------------------------------------------------
794
795// Compute |cos(theta)| between the incoming ray and the surface normal of the hit primitive.
796// Mirrors the normal construction in __intersection__patch (3rd transform column for patch/tile,
797// edge cross-product for triangle). Returns a value in (0, 1]; degenerate cases return 1.
798static __forceinline__ __device__ float coverCosTheta(uint32_t hit_position, const float3 &ray_dir) {
799 const uint32_t ptype = params.primitive_type[hit_position];
800 float T[16];
801 loadTransformMatrix(hit_position, T);
802 float3 normal;
803 if (ptype == 1) { // triangle
804 const float3 v0 = make_float3(T[3], T[7], T[11]);
805 const float3 v1 = make_float3(T[1] + T[3], T[5] + T[7], T[9] + T[11]);
806 const float3 v2 = make_float3(T[0] + T[1] + T[3], T[4] + T[5] + T[7], T[8] + T[9] + T[11]);
807 const float3 e0 = make_float3(v1.x - v0.x, v1.y - v0.y, v1.z - v0.z);
808 const float3 e1 = make_float3(v2.x - v0.x, v2.y - v0.y, v2.z - v0.z);
809 normal = cross(e0, e1);
810 } else { // patch or tile: normal is the 3rd column of the rotation part
811 normal = make_float3(T[2], T[6], T[10]);
812 }
813 const float nmag = d_magnitude(normal);
814 if (nmag < 1e-12f) return 1.f;
815 return fminf(1.f, fabsf(dot(ray_dir, normal)) / nmag);
816}
817
818// Shared any-hit body for translucent covers. For a glass hit it accumulates the per-band
819// transmittance (so the ray continues attenuated toward the source/sky), deposits the cover's own
820// reflected and absorbed energy into its buffers, and ignores the intersection so traversal
821// continues. Non-glass hits are left to the normal closest-hit (occlusion) path.
822// incoming_flux_for_band(b, b_global): the irradiance the ray delivers to the cover for launch-band
823// b BEFORE this cover's attenuation (i.e. prd->strength * external_flux * prior cover_transmittance).
824static __forceinline__ __device__ void coverAnyHitBody(bool is_direct) {
825 const uint32_t hit_uuid = optixGetAttribute_0();
826 const uint32_t hit_position = params.primitive_positions[hit_uuid];
827 if (hit_position == UINT_MAX) return;
828 if (params.is_glass == nullptr) return; // glass model not configured this launch
829
830 PerRayData *prd = getPayloadPRD();
831
832 // Never let a cover attenuate a ray it launched itself.
833 if (hit_uuid == prd->origin_UUID) { optixIgnoreIntersection(); return; }
834
835 const uint32_t Nprims = params.Nprimitives;
836 const uint32_t Nbands_global = params.Nbands_global;
837 const uint32_t Nbands_launch = params.Nbands_launch;
838
839 // A primitive is treated as a translucent cover (pass-through) only if it is glass in EVERY launched
840 // band. A single ray carries all launched bands, so a primitive that is glass in some launched bands
841 // but opaque in others cannot be both transmitted and blocked — we treat it as a normal (opaque)
842 // occluder for the whole ray (accept the hit). This matches the OptiX 6 and Vulkan backends.
843 {
844 int b_check = -1;
845 for (uint32_t b_global = 0; b_global < Nbands_global; b_global++) {
846 if (!params.band_launch_flag[b_global]) continue;
847 b_check++;
848 const uint32_t mat_ind = prd->source_ID * Nprims * Nbands_global + hit_position * Nbands_global + b_global;
849 if (params.is_glass[mat_ind] == 0) return; // opaque in this launched band → block the ray
850 }
851 }
852
853 const bool face_top = (optixGetAttribute_1() == 1u);
854 const float3 ray_dir = optixGetWorldRayDirection();
855 const float cos_theta = coverCosTheta(hit_position, ray_dir);
856
857 int b = -1;
858 for (uint32_t b_global = 0; b_global < Nbands_global; b_global++) {
859 if (!params.band_launch_flag[b_global]) continue;
860 b++;
861 if (b >= HELIOS_MAX_RADIATION_BANDS) break; // host guarantees this never trips
862
863 const uint32_t mat_ind = prd->source_ID * Nprims * Nbands_global
864 + hit_position * Nbands_global + b_global;
865
866 const float n = params.glass_n[mat_ind];
867 const float KL = params.glass_KL[mat_ind];
868 const float3 tra = glass_tau_rho_alpha(cos_theta, n, KL); // (tau, rho, alpha)
869
870 // Irradiance delivered to the cover for this band, before this cover's attenuation.
871 double ext_flux = 0.0;
872 if (is_direct) {
873 const uint32_t flux_idx = prd->source_ID * Nbands_launch + (uint32_t)b;
874 ext_flux = (double)params.source_fluxes[flux_idx];
875 } else {
876 // Diffuse: the cover's reflected/absorbed bookkeeping uses the isotropic sky flux as an
877 // approximation of the per-ray incident irradiance (the full miss-program angular
878 // distribution is not reproduced here). The transmittance accumulation below is exact and
879 // flux-independent; only this cover ρ/α deposit is approximate for anisotropic skies.
880 ext_flux = (params.diffuse_flux != nullptr) ? (double)params.diffuse_flux[b] : 0.0;
881 }
882 const double incoming = prd->strength * ext_flux * (double)prd->cover_transmittance[b];
883
884 // radiation_in / scatter buffers are indexed [prim * Nbands_launch + band_launch].
885 const uint32_t cov_ind = hit_position * Nbands_launch + (uint32_t)b;
886
887 if (incoming > 0.0) {
888 atomicFloatAdd(&params.radiation_in[cov_ind], (float)(incoming * (double)tra.z)); // absorbed
889 // Reflection leaves from the hit face; transmission continues out the far face.
890 if (face_top) {
891 atomicFloatAdd(&params.scatter_buff_top[cov_ind], (float)(incoming * (double)tra.y));
892 } else {
893 atomicFloatAdd(&params.scatter_buff_bottom[cov_ind], (float)(incoming * (double)tra.y));
894 }
895 }
896
897 // Attenuate the through-beam for this band.
898 prd->cover_transmittance[b] *= tra.x;
899 }
900
901 // All launched bands are glass: pass the (attenuated) ray through toward the source/sky.
902 optixIgnoreIntersection();
903}
904
905extern "C" __global__ void __anyhit__direct() {
906 coverAnyHitBody(true);
907}
908
909extern "C" __global__ void __anyhit__diffuse() {
910 coverAnyHitBody(false);
911}
912
913extern "C" __global__ void __closesthit__direct() {
914 // A direct ray hit an obstacle — it is simply blocked.
915 // Exception: bbox (type-5) hits trigger periodic boundary wrapping.
916 if (params.periodic_flag.x == 0 && params.periodic_flag.y == 0) return;
917
918 const uint32_t hit_uuid = optixGetAttribute_0();
919 if (params.bbox_UUID_base == 0xFFFFFFFFu || hit_uuid < params.bbox_UUID_base) return;
920
921 PerRayData *prd = getPayloadPRD();
922 prd->hit_periodic_boundary = true;
923 handlePeriodicBoundaryHit(prd, hit_uuid);
924}
925
926// ---------------------------------------------------------------------------
927// Closest-hit: diffuse radiation
928// ---------------------------------------------------------------------------
929
930extern "C" __global__ void __closesthit__diffuse() {
931 const uint32_t hit_uuid = optixGetAttribute_0();
932 const bool face_top = (optixGetAttribute_1() == 1u);
933
934 PerRayData *prd = getPayloadPRD();
935
936 // Periodic boundary: bbox (type-5) hit — wrap ray and return without depositing energy
937 if ((params.periodic_flag.x == 1 || params.periodic_flag.y == 1) &&
938 params.bbox_UUID_base != 0xFFFFFFFFu && hit_uuid >= params.bbox_UUID_base) {
939 prd->hit_periodic_boundary = true;
940 handlePeriodicBoundaryHit(prd, hit_uuid);
941 return;
942 }
943
944 const uint32_t origin_position = params.primitive_positions[prd->origin_UUID];
945 const uint32_t hit_position = params.primitive_positions[hit_uuid];
946 if (origin_position == UINT_MAX || hit_position == UINT_MAX) return;
947
948 const uint32_t Nprims = params.Nprimitives;
949 const uint32_t Nbands_global = params.Nbands_global;
950 const uint32_t Nbands_launch = params.Nbands_launch;
951
952 int b = -1;
953 for (uint32_t b_global = 0; b_global < Nbands_global; b_global++) {
954 if (!params.band_launch_flag[b_global]) continue;
955 b++;
956
957 const uint32_t ind_origin = origin_position * Nbands_launch + (uint32_t)b;
958 const uint32_t ind_hit = hit_position * Nbands_launch + (uint32_t)b;
959
960 const double strength = face_top ? params.radiation_out_top[ind_hit] * prd->strength
961 : params.radiation_out_bottom[ind_hit] * prd->strength;
962 if (strength == 0.0) continue;
963
964 const uint32_t radprop_ind = prd->source_ID * Nprims * Nbands_global
965 + origin_position * Nbands_global + b_global;
966 const float t_rho = params.rho[radprop_ind];
967 const float t_tau = params.tau[radprop_ind];
968
969 atomicFloatAdd(&params.radiation_in[ind_origin], (float)(strength * (1.0 - t_rho - t_tau)));
970
971 if (t_rho > 0.f || t_tau > 0.f) {
972 if (prd->face) { // top-face origin
973 atomicFloatAdd(&params.scatter_buff_top[ind_origin], (float)(strength * t_rho));
974 atomicFloatAdd(&params.scatter_buff_bottom[ind_origin], (float)(strength * t_tau));
975 } else { // bottom-face origin
976 atomicFloatAdd(&params.scatter_buff_bottom[ind_origin], (float)(strength * t_rho));
977 atomicFloatAdd(&params.scatter_buff_top[ind_origin], (float)(strength * t_tau));
978 }
979 }
980
981 // Camera-weighted scatter: mirrors scatter_buff but uses rho_cam/tau_cam
982 if (params.Ncameras > 0 && params.rho_cam && params.scatter_buff_top_cam) {
983 const uint32_t Ncameras = params.Ncameras;
984 const uint32_t cam_id = params.camera_ID;
985 const uint32_t rc_idx = prd->source_ID * Nprims * Nbands_global * Ncameras
986 + origin_position * Nbands_global * Ncameras
987 + b_global * Ncameras + cam_id;
988 const float t_rho_cam = params.rho_cam[rc_idx];
989 const float t_tau_cam = params.tau_cam ? params.tau_cam[rc_idx] : 0.f;
990 if ((t_rho_cam > 0.f || t_tau_cam > 0.f) && strength > 0.0) {
991 if (prd->face) {
992 atomicFloatAdd(&params.scatter_buff_top_cam[ind_origin], (float)(strength * t_rho_cam));
993 atomicFloatAdd(&params.scatter_buff_bottom_cam[ind_origin], (float)(strength * t_tau_cam));
994 } else {
995 atomicFloatAdd(&params.scatter_buff_bottom_cam[ind_origin], (float)(strength * t_rho_cam));
996 atomicFloatAdd(&params.scatter_buff_top_cam[ind_origin], (float)(strength * t_tau_cam));
997 }
998 }
999 }
1000 }
1001}
1002
1003// ---------------------------------------------------------------------------
1004// Closest-hit: camera
1005// ---------------------------------------------------------------------------
1006
1007extern "C" __global__ void __closesthit__camera() {
1008 const uint32_t hit_uuid = optixGetAttribute_0();
1009 PerRayData *prd = getPayloadPRD();
1010 const uint32_t hit_position = params.primitive_positions[hit_uuid];
1011
1012 if (hit_position == 0xFFFFFFFFu) return; // invalid primitive
1013
1014 // Periodic boundary: treat as transparent wall and re-launch
1015 if ((params.periodic_flag.x != 0.f || params.periodic_flag.y != 0.f) &&
1016 params.primitive_type[hit_position] == 5) {
1017 handlePeriodicBoundaryHit(prd, hit_uuid);
1018 prd->hit_periodic_boundary = true;
1019 return;
1020 }
1021
1022 const uint32_t pixel_index = prd->origin_UUID;
1023 const uint32_t Nbands_l = params.Nbands_launch;
1024 const uint32_t Nprims = params.Nprimitives;
1025 const float t_hit = optixGetRayTmax();
1026 const float3 ray_origin = optixGetWorldRayOrigin();
1027 const float3 ray_direction = optixGetWorldRayDirection();
1028
1029 // Use face attribute from intersection program (correct for both patches and triangles)
1030 const bool face_top = (optixGetAttribute_1() == 1u);
1031
1032 // Compute surface normal for specular reflection (needed for Blinn-Phong half-vector)
1033 float T[16];
1034 loadTransformMatrix(hit_position, T);
1035 float3 n0 = make_float3(0.f, 0.f, 0.f); d_transformPoint(T, n0);
1036 float3 n1 = make_float3(1.f, 0.f, 0.f); d_transformPoint(T, n1);
1037 float3 n2 = make_float3(0.f, 1.f, 0.f); d_transformPoint(T, n2);
1038 float3 normal = normalize(cross(n1 - n0, n2 - n0));
1039 // Ensure normal points toward the camera (consistent with face_top)
1040 if (face_top != (dot(normal, ray_direction) < 0.f)) {
1041 normal = make_float3(-normal.x, -normal.y, -normal.z);
1042 }
1043
1044 for (uint32_t b = 0; b < Nbands_l; b++) {
1045 // Radiance from hit surface (outgoing flux / pi = radiance)
1046 const uint32_t ind_hit = hit_position * Nbands_l + b;
1047 float strength = (float)prd->strength *
1048 (float)(face_top ? params.radiation_out_top[ind_hit]
1049 : params.radiation_out_bottom[ind_hit]);
1050
1051 // Check sources visible between camera origin and hit point
1052 for (uint32_t s = 0; s < params.Nsources; s++) {
1053 const float flux = params.source_fluxes[s * Nbands_l + b];
1054 if (flux <= 0.0f) continue;
1055
1056 const uint32_t stype = params.source_types[s];
1057 float source_radiance = 0.0f;
1058
1059 if (stype == 1) {
1060 // Sphere
1061 const float radius = params.source_widths[s].x * 0.5f;
1062 const float3 oc = make_float3(ray_origin.x - params.source_positions[s].x,
1063 ray_origin.y - params.source_positions[s].y,
1064 ray_origin.z - params.source_positions[s].z);
1065 const float bd = dot(oc, ray_direction);
1066 const float cd = dot(oc, oc) - radius * radius;
1067 const float disc = bd * bd - cd;
1068 if (disc >= 0.0f) {
1069 const float t_sphere = -bd - sqrtf(disc);
1070 if (t_sphere > 0.0f && t_sphere < t_hit) {
1071 const float area = 4.0f * M_PI * radius * radius;
1072 source_radiance = (flux / area) / M_PI;
1073 }
1074 }
1075 } else if (stype == 3) {
1076 // Rectangle
1077 float trans[16];
1078 d_makeTransformMatrix(params.source_rotations[s], trans);
1079 const float3 snormal = make_float3(trans[2], trans[6], trans[10]);
1080 const float denom = dot(ray_direction, snormal);
1081 if (denom < -1e-6f) {
1082 const float3 oc = make_float3(params.source_positions[s].x - ray_origin.x,
1083 params.source_positions[s].y - ray_origin.y,
1084 params.source_positions[s].z - ray_origin.z);
1085 const float t_r = dot(oc, snormal) / denom;
1086 if (t_r > 0.0f && t_r < t_hit) {
1087 float3 hp = make_float3(ray_origin.x + t_r * ray_direction.x - params.source_positions[s].x,
1088 ray_origin.y + t_r * ray_direction.y - params.source_positions[s].y,
1089 ray_origin.z + t_r * ray_direction.z - params.source_positions[s].z);
1090 float inv_t[16];
1091 d_invertMatrix(trans, inv_t);
1092 d_transformPoint(inv_t, hp);
1093 if (fabsf(hp.x) <= params.source_widths[s].x * 0.5f &&
1094 fabsf(hp.y) <= params.source_widths[s].y * 0.5f) {
1095 const float area = params.source_widths[s].x * params.source_widths[s].y;
1096 source_radiance = (flux / area) * (-denom) / M_PI;
1097 }
1098 }
1099 }
1100 } else if (stype == 4) {
1101 // Disk
1102 float trans[16];
1103 d_makeTransformMatrix(params.source_rotations[s], trans);
1104 const float3 snormal = make_float3(trans[2], trans[6], trans[10]);
1105 const float denom = dot(ray_direction, snormal);
1106 if (denom < -1e-6f) {
1107 const float3 oc = make_float3(params.source_positions[s].x - ray_origin.x,
1108 params.source_positions[s].y - ray_origin.y,
1109 params.source_positions[s].z - ray_origin.z);
1110 const float t_d = dot(oc, snormal) / denom;
1111 if (t_d > 0.0f && t_d < t_hit) {
1112 const float3 hp = make_float3(ray_origin.x + t_d * ray_direction.x - params.source_positions[s].x,
1113 ray_origin.y + t_d * ray_direction.y - params.source_positions[s].y,
1114 ray_origin.z + t_d * ray_direction.z - params.source_positions[s].z);
1115 const float radius = params.source_widths[s].x;
1116 if (dot(hp, hp) <= radius * radius) {
1117 const float area = M_PI * radius * radius;
1118 source_radiance = (flux / area) * (-denom) / M_PI;
1119 }
1120 }
1121 }
1122 }
1123
1124 if (source_radiance > 0.0f) {
1125 strength += source_radiance * (float)prd->strength;
1126 }
1127 }
1128
1129 // Specular contribution (only if enabled and on iteration 0)
1130 float strength_spec = 0.0f;
1131 if (params.specular_reflection_enabled > 0 &&
1132 params.specular_exponent && params.specular_exponent[hit_position] > 0.f &&
1133 params.scattering_iteration == 0 &&
1134 params.radiation_specular) {
1135 for (uint32_t rr = 0; rr < params.Nsources; rr++) {
1136 const uint32_t ind_spec = rr * params.Ncameras * Nprims * Nbands_l
1137 + params.camera_ID * Nprims * Nbands_l
1138 + hit_position * Nbands_l + b;
1139 const float spec = params.radiation_specular[ind_spec] * 0.25f;
1140 if (spec > 0.0f) {
1141 float3 light_dir;
1142 if (params.source_types[rr] == 0 || params.source_types[rr] == 2) {
1143 light_dir = normalize(params.source_positions[rr]);
1144 } else {
1145 const float3 hp = make_float3(ray_origin.x + t_hit * ray_direction.x,
1146 ray_origin.y + t_hit * ray_direction.y,
1147 ray_origin.z + t_hit * ray_direction.z);
1148 light_dir = normalize(make_float3(params.source_positions[rr].x - hp.x,
1149 params.source_positions[rr].y - hp.y,
1150 params.source_positions[rr].z - hp.z));
1151 }
1152 const float3 spec_dir = normalize(light_dir - ray_direction);
1153 const float exponent = params.specular_exponent[hit_position];
1154 float scale_coeff = 1.0f;
1155 if (params.specular_reflection_enabled == 2 && params.specular_scale) {
1156 scale_coeff = params.specular_scale[hit_position];
1157 }
1158 const float cos_spec = fmaxf(0.f, dot(spec_dir, normal));
1159 strength_spec += spec * scale_coeff
1160 * powf(cos_spec, exponent) * (exponent + 2.f)
1161 / ((float)params.launch_dim_x * 2.f * M_PI);
1162 }
1163 }
1164 }
1165
1166 // Accumulate into camera radiation buffer: [pixel][band]
1167 atomicFloatAdd(&params.radiation_in_camera[pixel_index * Nbands_l + b],
1168 (strength + strength_spec) / M_PI);
1169 }
1170}
1171
1172// ---------------------------------------------------------------------------
1173// Closest-hit: pixel label
1174// ---------------------------------------------------------------------------
1175
1176extern "C" __global__ void __closesthit__pixel_label() {
1177 const uint32_t hit_uuid = optixGetAttribute_0();
1178 PerRayData *prd = getPayloadPRD();
1179 const uint32_t origin_UUID = prd->origin_UUID;
1180 const uint32_t hit_position = params.primitive_positions[hit_uuid];
1181
1182 // Periodic boundary: treat as transparent wall and re-launch
1183 if ((params.periodic_flag.x != 0.f || params.periodic_flag.y != 0.f) &&
1184 hit_position != 0xFFFFFFFFu && params.primitive_type[hit_position] == 5) {
1185 handlePeriodicBoundaryHit(prd, hit_uuid);
1186 prd->hit_periodic_boundary = true;
1187 return;
1188 }
1189
1190 // Store UUID+1 (0 is reserved for sky/miss)
1191 if (params.camera_pixel_label) {
1192 params.camera_pixel_label[origin_UUID] = hit_uuid + 1u;
1193 }
1194
1195 // Depth: project ray parameter along camera view direction
1196 if (params.camera_pixel_depth) {
1197 const float t_hit = optixGetRayTmax() + (float)prd->strength; // strength=0 for pixel label
1198 const float3 ray_dir = optixGetWorldRayDirection();
1199 const float3 cam_dir = d_rotatePoint(make_float3(1.f, 0.f, 0.f),
1200 -0.5f * M_PI + params.camera_direction.x,
1201 0.5f * M_PI - params.camera_direction.y);
1202 params.camera_pixel_depth[origin_UUID] = fabsf(dot(cam_dir, ray_dir)) * t_hit;
1203 }
1204}
1205
1206// ---------------------------------------------------------------------------
1207// Raygen: direct rays
1208// ---------------------------------------------------------------------------
1209
1210extern "C" __global__ void __raygen__direct() {
1211 // 3D launch: x = x-strat index [0, dim_x), y = y-strat index [0, dim_y), z = prim
1212 const uint3 idx = optixGetLaunchIndex();
1213 const uint32_t xi = idx.x;
1214 const uint32_t yi = idx.y;
1215 const uint32_t prim_local = idx.z;
1216
1217 const uint32_t dim_x = params.launch_dim_x;
1218 const uint32_t dim_y = params.launch_dim_y;
1219 const uint32_t Nrays = dim_x * dim_y;
1220 const uint32_t prim_pos = params.launch_offset + prim_local;
1221
1222 if (prim_pos >= params.Nprimitives) return;
1223
1224 const uint32_t ptype = params.primitive_type[prim_pos];
1225 const int32_t NX = params.object_subdivisions[prim_pos * 2];
1226 const int32_t NY = params.object_subdivisions[prim_pos * 2 + 1];
1227
1228 float T[16];
1229 loadTransformMatrix(prim_pos, T);
1230
1231 const uint32_t linear_idx = xi + dim_x * yi;
1232 uint32_t seed = tea<16>(linear_idx + Nrays * prim_local, params.random_seed);
1233
1234 for (int jj = 0; jj < NY; jj++) {
1235 for (int ii = 0; ii < NX; ii++) {
1236
1237 // UUID for this sub-patch
1238 const uint32_t UUID = params.primitiveID[prim_pos] + (uint32_t)(jj * NX + ii);
1239
1240 const float Rx = rnd(seed);
1241 const float Ry = rnd(seed);
1242
1243 float3 sp;
1244 float3 normal;
1245
1246 if (ptype == 0 || ptype == 3) { // Patch or Tile: canonical space [-0.5, 0.5]^2
1247 const float dx = 1.0f / float(NX);
1248 const float dy = 1.0f / float(NY);
1249 sp.x = -0.5f + ii * dx + (float(xi) + Rx) * dx / float(dim_x);
1250 sp.y = -0.5f + jj * dy + (float(yi) + Ry) * dy / float(dim_y);
1251 sp.z = 0.f;
1252
1253 // Origin mask rejection sampling: only launch from opaque regions (matches OptiX 6 raygen behavior)
1254 const int32_t msk_orig = params.mask_IDs[prim_pos];
1255 if (msk_orig >= 0) {
1256 for (int attempt = 0; attempt < 10; attempt++) {
1257 float uv_u, uv_v;
1258 if (params.uv_IDs[prim_pos] >= 0) {
1259 float2 uv0 = params.uv_data[prim_pos * 4 + 0];
1260 float2 uv1 = params.uv_data[prim_pos * 4 + 1];
1261 float2 uv2 = params.uv_data[prim_pos * 4 + 2];
1262 uv_u = uv0.x + (sp.x + 0.5f) * (uv1.x - uv0.x);
1263 uv_v = uv0.y + (sp.y + 0.5f) * (uv2.y - uv0.y);
1264 } else {
1265 uv_u = sp.x + 0.5f;
1266 uv_v = sp.y + 0.5f;
1267 }
1268 if (sampleMask(msk_orig, uv_u, uv_v)) break;
1269 sp.x = -0.5f + (ii + rnd(seed)) * dx;
1270 sp.y = -0.5f + (jj + rnd(seed)) * dy;
1271 }
1272 }
1273
1274 float3 v0 = make_float3(0.f, 0.f, 0.f); d_transformPoint(T, v0);
1275 float3 v1 = make_float3(1.f, 0.f, 0.f); d_transformPoint(T, v1);
1276 float3 v2 = make_float3(0.f, 1.f, 0.f); d_transformPoint(T, v2);
1277 normal = normalize(cross(v1 - v0, v2 - v0));
1278
1279 } else if (ptype == 1) { // Triangle: canonical space (0,0,0)-(0,1,0)-(1,1,0)
1280 if (Rx < Ry) { sp.x = Rx; sp.y = Ry; }
1281 else { sp.x = Ry; sp.y = Rx; }
1282 sp.z = 0.f;
1283
1284 // Origin mask rejection sampling
1285 const int32_t msk_orig_t = params.mask_IDs[prim_pos];
1286 if (msk_orig_t >= 0) {
1287 for (int attempt = 0; attempt < 10; attempt++) {
1288 float uv_u, uv_v;
1289 if (params.uv_IDs[prim_pos] >= 0) {
1290 float2 uv0 = params.uv_data[prim_pos * 4 + 0];
1291 float2 uv1 = params.uv_data[prim_pos * 4 + 1];
1292 float2 uv2 = params.uv_data[prim_pos * 4 + 2];
1293 const float beta = sp.y - sp.x; // weight at v1
1294 const float gamma = sp.x; // weight at v2
1295 float2 uv = make_float2(
1296 uv0.x + beta * (uv1.x - uv0.x) + gamma * (uv2.x - uv0.x),
1297 uv0.y + beta * (uv1.y - uv0.y) + gamma * (uv2.y - uv0.y));
1298 uv_u = uv.x;
1299 uv_v = 1.f - uv.y; // Y-flip (matches intersection convention)
1300 } else {
1301 uv_u = sp.y; // = beta + gamma
1302 uv_v = sp.x; // = gamma
1303 }
1304 if (sampleMask(msk_orig_t, uv_u, uv_v)) break;
1305 float Rx2 = rnd(seed), Ry2 = rnd(seed);
1306 if (Rx2 < Ry2) { sp.x = Rx2; sp.y = Ry2; }
1307 else { sp.x = Ry2; sp.y = Rx2; }
1308 }
1309 }
1310
1311 float3 v0 = make_float3(0.f, 0.f, 0.f); d_transformPoint(T, v0);
1312 float3 v1 = make_float3(0.f, 1.f, 0.f); d_transformPoint(T, v1);
1313 float3 v2 = make_float3(1.f, 1.f, 0.f); d_transformPoint(T, v2);
1314 normal = normalize(cross(v1 - v0, v2 - v0));
1315
1316 } else {
1317 // Unsupported primitive type — should have been caught by updateGeometry().
1318 printf("ERROR (OptiX8 __raygen__direct): unsupported primitive type %u at index %u\n",
1319 ptype, prim_pos);
1320 __trap();
1321 }
1322
1323 // Transform sample point to world space
1324 float3 ray_origin = sp;
1325 d_transformPoint(T, ray_origin);
1326
1327 // Send a ray toward each source
1328 for (uint32_t rr = 0; rr < params.Nsources; rr++) {
1329
1330 const uint32_t src_type = params.source_types[rr];
1331
1332 float3 ray_direction;
1333 float ray_tmax;
1334 double strength;
1335
1336 if (src_type == 0) { // Collimated source
1337 ray_direction = normalize(params.source_positions[rr]);
1338 ray_tmax = 1e38f;
1339 strength = (1.0 / double(dim_x * dim_y)) * (double)fabsf(dot(normal, ray_direction));
1340
1341 } else if (src_type == 1 || src_type == 2) { // Sphere source (type 1 = point, type 2 = sphere)
1342 float theta_s = acos_safe(1.f - 2.f * rnd(seed));
1343 float phi_s = rnd(seed) * 2.f * M_PI;
1344 float3 sphere_pt = make_float3(0.5f * params.source_widths[rr].x * sinf(theta_s) * cosf(phi_s),
1345 0.5f * params.source_widths[rr].x * sinf(theta_s) * sinf(phi_s),
1346 0.5f * params.source_widths[rr].x * cosf(theta_s));
1347 ray_direction = sphere_pt + params.source_positions[rr] - ray_origin;
1348 ray_tmax = d_magnitude(ray_direction);
1349 ray_direction = normalize(ray_direction);
1350
1351 // Integrate over sphere surface for strength
1352 strength = 0.0;
1353 const uint32_t N = 10;
1354 for (uint32_t j = 0; j < N; j++) {
1355 for (uint32_t i = 0; i < N; i++) {
1356 float theta = acos_safe(1.f - 2.f * (float(i) + 0.5f) / float(N));
1357 float phi = (float(j) + 0.5f) * 2.f * M_PI / float(N);
1358 float3 ldir = make_float3(sinf(theta)*cosf(phi), sinf(theta)*sinf(phi), cosf(theta));
1359 if (dot(ldir, ray_direction) < 0.f) {
1360 strength += (1.0 / double(dim_x * dim_y)) * (double)fabsf(dot(normal, ray_direction))
1361 * (double)fabsf(dot(ldir, ray_direction))
1362 / ((double)ray_tmax * (double)ray_tmax)
1363 / double(N * N)
1364 * (double)params.source_widths[rr].x * (double)params.source_widths[rr].x;
1365 }
1366 }
1367 }
1368
1369 } else if (src_type == 3) { // Rectangle source
1370 float light_transform[16];
1371 float3 rot3 = params.source_rotations[rr];
1372 d_makeTransformMatrix(rot3, light_transform);
1373
1374 float3 square_pt;
1375 d_sampleSquare(seed, square_pt);
1376 square_pt = make_float3(params.source_widths[rr].x * square_pt.x, params.source_widths[rr].y * square_pt.y, square_pt.z);
1377 d_transformPoint(light_transform, square_pt);
1378
1379 float3 light_dir = make_float3(0.f, 0.f, 1.f);
1380 d_transformPoint(light_transform, light_dir);
1381
1382 ray_direction = square_pt + params.source_positions[rr] - ray_origin;
1383 if (dot(ray_direction, light_dir) > 0.f) continue; // don't emit from back of source
1384
1385 ray_tmax = d_magnitude(ray_direction);
1386 ray_direction = normalize(ray_direction);
1387 strength = (1.0 / double(dim_x * dim_y))
1388 * (double)fabsf(dot(normal, ray_direction))
1389 * (double)fabsf(dot(light_dir, ray_direction))
1390 / ((double)ray_tmax * (double)ray_tmax)
1391 * (double)params.source_widths[rr].x * (double)params.source_widths[rr].y
1392 / M_PI;
1393
1394 } else if (src_type == 4) { // Disk source
1395 float light_transform[16];
1396 float3 rot3 = params.source_rotations[rr];
1397 d_makeTransformMatrix(rot3, light_transform);
1398
1399 float3 disk_pt;
1400 d_sampleDisk(seed, disk_pt);
1401 d_transformPoint(light_transform, disk_pt);
1402
1403 float3 light_dir = make_float3(0.f, 0.f, 1.f);
1404 d_transformPoint(light_transform, light_dir);
1405
1406 ray_direction = params.source_widths[rr].x * disk_pt + params.source_positions[rr] - ray_origin;
1407 if (dot(ray_direction, light_dir) > 0.f) continue; // don't emit from back of source
1408
1409 ray_tmax = d_magnitude(ray_direction);
1410 ray_direction = normalize(ray_direction);
1411 strength = (1.0 / double(dim_x * dim_y))
1412 * (double)fabsf(dot(normal, ray_direction))
1413 * (double)fabsf(dot(light_dir, ray_direction))
1414 / ((double)ray_tmax * (double)ray_tmax)
1415 * (double)params.source_widths[rr].x * (double)params.source_widths[rr].x;
1416
1417 } else {
1418 continue; // Unknown source type
1419 }
1420
1421 PerRayData prd;
1422 prd.seed = seed;
1423 prd.origin_UUID = UUID;
1424 prd.source_ID = (unsigned char)rr;
1425 prd.hit_periodic_boundary = false;
1426 prd.face = (dot(ray_direction, normal) > 0.f);
1427 initCoverTransmittance(prd); // translucent-cover attenuation starts at 1 (no covers crossed)
1428
1429 // Strength set above per source type
1430 prd.strength = strength;
1431
1432 // Only fire from the face pointing toward source (or two-sided)
1433 const int8_t tsf = params.twosided_flag[prim_pos];
1434 if (!prd.face && tsf == 0) continue;
1435 if (tsf == 3) continue; // reserved flag — skip
1436
1437 uint32_t u0, u1;
1438 float3 current_origin = ray_origin;
1439
1440 for (int wrap = 0; wrap < 10; ++wrap) {
1441 packPointer(&prd, u0, u1);
1442 optixTrace(
1443 params.traversable,
1444 current_origin,
1445 ray_direction,
1446 1e-4f, // tmin
1447 ray_tmax, // tmax
1448 0.f, // time
1449 OptixVisibilityMask(255),
1450 OPTIX_RAY_FLAG_NONE,
1451 0, // SBT offset (direct hit group = 0)
1452 0, // SBT stride
1453 0, // miss SBT index (direct miss = 0)
1454 u0, u1
1455 );
1456 if (!prd.hit_periodic_boundary) break;
1457 current_origin = prd.periodic_hit;
1458 prd.hit_periodic_boundary = false;
1459 }
1460
1461 seed = prd.seed;
1462 }
1463 }
1464 }
1465}
1466
1467// ---------------------------------------------------------------------------
1468// Raygen: diffuse rays
1469// ---------------------------------------------------------------------------
1470
1471extern "C" __global__ void __raygen__diffuse() {
1472 // 3D launch: x=theta_idx, y=phi_idx, z=prim_local
1473 const uint3 idx = optixGetLaunchIndex();
1474 const uint32_t theta_idx = idx.x;
1475 const uint32_t phi_idx = idx.y;
1476 const uint32_t prim_local = idx.z;
1477
1478 const uint32_t dim_x = params.launch_dim_x;
1479 const uint32_t dim_y = params.launch_dim_y;
1480 const uint32_t dimxy = dim_x * dim_y;
1481
1482 const uint32_t prim_pos = params.launch_offset + prim_local;
1483 if (prim_pos >= params.Nprimitives) return;
1484
1485 // Skip bottom-face launch for one-sided primitives
1486 if (params.launch_face == 0 && params.twosided_flag[prim_pos] == 0) return;
1487
1488 const uint32_t ptype = params.primitive_type[prim_pos];
1489 const int32_t NX = params.object_subdivisions[prim_pos * 2];
1490 const int32_t NY = params.object_subdivisions[prim_pos * 2 + 1];
1491
1492 float T[16];
1493 loadTransformMatrix(prim_pos, T);
1494
1495 // Seed once per ray index
1496 const uint32_t linear_idx = theta_idx + dim_x * phi_idx;
1497 uint32_t seed = tea<16>(linear_idx + dimxy * prim_local, params.random_seed);
1498
1499 // Stratified cosine-weighted hemisphere sampling
1500 const float Rt = (theta_idx + rnd(seed)) / float(dim_x);
1501 const float Rp = (phi_idx + rnd(seed)) / float(dim_y);
1502 const float t = asin_safe(sqrtf(Rt));
1503 const float p = 2.f * M_PI * Rp;
1504 float3 ray_dir_canonical;
1505 ray_dir_canonical.x = sinf(t) * cosf(p);
1506 ray_dir_canonical.y = sinf(t) * sinf(p);
1507 ray_dir_canonical.z = cosf(t);
1508
1509 for (int jj = 0; jj < NY; jj++) {
1510 for (int ii = 0; ii < NX; ii++) {
1511
1512 const uint32_t UUID = params.primitiveID[prim_pos] + (uint32_t)(jj * NX + ii);
1513
1514 const float Rx = rnd(seed);
1515 const float Ry = rnd(seed);
1516
1517 float3 sp;
1518 float3 normal;
1519
1520 if (ptype == 0 || ptype == 3) { // Patch or Tile
1521 const float dx = 1.f / float(NX);
1522 const float dy = 1.f / float(NY);
1523 sp.x = -0.5f + (ii + Rx) * dx;
1524 sp.y = -0.5f + (jj + Ry) * dy;
1525 sp.z = 0.f;
1526
1527 // Origin mask rejection sampling
1528 const int32_t msk_orig_d = params.mask_IDs[prim_pos];
1529 if (msk_orig_d >= 0) {
1530 for (int attempt = 0; attempt < 10; attempt++) {
1531 float uv_u, uv_v;
1532 if (params.uv_IDs[prim_pos] >= 0) {
1533 float2 uv0 = params.uv_data[prim_pos * 4 + 0];
1534 float2 uv1 = params.uv_data[prim_pos * 4 + 1];
1535 float2 uv2 = params.uv_data[prim_pos * 4 + 2];
1536 uv_u = uv0.x + (sp.x + 0.5f) * (uv1.x - uv0.x);
1537 uv_v = uv0.y + (sp.y + 0.5f) * (uv2.y - uv0.y);
1538 } else {
1539 uv_u = sp.x + 0.5f;
1540 uv_v = sp.y + 0.5f;
1541 }
1542 if (sampleMask(msk_orig_d, uv_u, uv_v)) break;
1543 sp.x = -0.5f + (ii + rnd(seed)) * dx;
1544 sp.y = -0.5f + (jj + rnd(seed)) * dy;
1545 }
1546 }
1547
1548 float3 v0 = make_float3(0.f, 0.f, 0.f); d_transformPoint(T, v0);
1549 float3 v1 = make_float3(1.f, 0.f, 0.f); d_transformPoint(T, v1);
1550 float3 v2 = make_float3(0.f, 1.f, 0.f); d_transformPoint(T, v2);
1551 normal = normalize(cross(v1 - v0, v2 - v0));
1552
1553 } else if (ptype == 1) { // Triangle
1554 if (Rx < Ry) { sp.x = Rx; sp.y = Ry; }
1555 else { sp.x = Ry; sp.y = Rx; }
1556 sp.z = 0.f;
1557
1558 // Origin mask rejection sampling
1559 const int32_t msk_orig_dt = params.mask_IDs[prim_pos];
1560 if (msk_orig_dt >= 0) {
1561 for (int attempt = 0; attempt < 10; attempt++) {
1562 float uv_u, uv_v;
1563 if (params.uv_IDs[prim_pos] >= 0) {
1564 float2 uv0 = params.uv_data[prim_pos * 4 + 0];
1565 float2 uv1 = params.uv_data[prim_pos * 4 + 1];
1566 float2 uv2 = params.uv_data[prim_pos * 4 + 2];
1567 const float beta = sp.y - sp.x;
1568 const float gamma = sp.x;
1569 float2 uv = make_float2(
1570 uv0.x + beta * (uv1.x - uv0.x) + gamma * (uv2.x - uv0.x),
1571 uv0.y + beta * (uv1.y - uv0.y) + gamma * (uv2.y - uv0.y));
1572 uv_u = uv.x;
1573 uv_v = 1.f - uv.y;
1574 } else {
1575 uv_u = sp.y;
1576 uv_v = sp.x;
1577 }
1578 if (sampleMask(msk_orig_dt, uv_u, uv_v)) break;
1579 float Rx2 = rnd(seed), Ry2 = rnd(seed);
1580 if (Rx2 < Ry2) { sp.x = Rx2; sp.y = Ry2; }
1581 else { sp.x = Ry2; sp.y = Rx2; }
1582 }
1583 }
1584
1585 float3 v0 = make_float3(0.f, 0.f, 0.f); d_transformPoint(T, v0);
1586 float3 v1 = make_float3(0.f, 1.f, 0.f); d_transformPoint(T, v1);
1587 float3 v2 = make_float3(1.f, 1.f, 0.f); d_transformPoint(T, v2);
1588 normal = normalize(cross(v1 - v0, v2 - v0));
1589
1590 } else {
1591 // Unsupported primitive type — should have been caught by updateGeometry().
1592 printf("ERROR (OptiX8 __raygen__diffuse): unsupported primitive type %u at index %u\n",
1593 ptype, prim_pos);
1594 __trap();
1595 }
1596
1597 // Rotate hemisphere direction by primitive normal orientation
1598 float3 ray_dir = d_rotatePoint(ray_dir_canonical,
1599 acos_safe(normal.z),
1600 atan2f(normal.y, normal.x));
1601
1602 // Transform origin point to world space
1603 float3 ray_origin = sp;
1604 d_transformPoint(T, ray_origin);
1605
1606 PerRayData prd;
1607 prd.seed = seed;
1608 prd.origin_UUID = UUID;
1609 prd.source_ID = 0;
1610 prd.hit_periodic_boundary = false;
1611 prd.strength = 1.0 / double(dimxy);
1612 initCoverTransmittance(prd); // translucent-cover attenuation starts at 1 (no covers crossed)
1613
1614 uint32_t u0, u1;
1615
1616 if (params.launch_face == 1 && params.twosided_flag[prim_pos] != 3) {
1617 prd.face = true;
1618 float3 cur_origin = ray_origin;
1619 for (int wrap = 0; wrap < 10; ++wrap) {
1620 packPointer(&prd, u0, u1);
1621 optixTrace(
1622 params.traversable,
1623 cur_origin, ray_dir,
1624 1e-4f, 1e38f, 0.f,
1625 OptixVisibilityMask(255),
1626 OPTIX_RAY_FLAG_NONE,
1627 1, 0, 1, // SBT offset=1 (diffuse hit), stride=0, miss index=1 (diffuse miss)
1628 u0, u1
1629 );
1630 if (!prd.hit_periodic_boundary) break;
1631 cur_origin = prd.periodic_hit;
1632 prd.hit_periodic_boundary = false;
1633 }
1634 } else if (params.launch_face == 0 && params.twosided_flag[prim_pos] == 1) {
1635 prd.face = false;
1636 float3 neg_dir = make_float3(-ray_dir.x, -ray_dir.y, -ray_dir.z);
1637 float3 cur_origin = ray_origin;
1638 for (int wrap = 0; wrap < 10; ++wrap) {
1639 packPointer(&prd, u0, u1);
1640 optixTrace(
1641 params.traversable,
1642 cur_origin, neg_dir,
1643 1e-4f, 1e38f, 0.f,
1644 OptixVisibilityMask(255),
1645 OPTIX_RAY_FLAG_NONE,
1646 1, 0, 1,
1647 u0, u1
1648 );
1649 if (!prd.hit_periodic_boundary) break;
1650 cur_origin = prd.periodic_hit;
1651 prd.hit_periodic_boundary = false;
1652 }
1653 }
1654
1655 seed = prd.seed;
1656 }
1657 }
1658}
1659
1660// ---------------------------------------------------------------------------
1661// Raygen: camera rays
1662// 3D launch: x=ray_within_pixel [0,anti_samples), y=tile_column, z=tile_row
1663// ---------------------------------------------------------------------------
1664
1665extern "C" __global__ void __raygen__camera() {
1666 const uint3 idx = optixGetLaunchIndex();
1667 const uint32_t ray_idx = idx.x; // sample index within pixel
1668 const uint32_t col = idx.y; // tile column
1669 const uint32_t row = idx.z; // tile row
1670
1671 const uint32_t dim_x = params.launch_dim_x; // antialiasing_samples
1672 const uint32_t dim_y = params.launch_dim_y; // tile_width
1673
1674 // Global pixel coordinates
1675 const uint32_t ii = (uint32_t)params.camera_pixel_offset.x + col;
1676 const uint32_t jj = (uint32_t)params.camera_pixel_offset.y + row;
1677
1678 // Linear pixel index in the full image
1679 const uint32_t pixel_index = jj * (uint32_t)params.camera_resolution_full.x + ii;
1680
1681 // Seed: unique per (ray_idx, col, row)
1682 const uint32_t linear_idx = dim_x * col + ray_idx;
1683 uint32_t seed = tea<16>(linear_idx + dim_x * dim_y * row, params.random_seed);
1684
1685 const float Rx = rnd(seed);
1686 const float Ry = rnd(seed);
1687
1688 // Map sub-pixel sample to view-space point on viewplane
1689 // sp.x = viewplane distance, sp.y/sp.z = horizontal/vertical offsets
1690 const float multiplier = 1.0f / params.FOV_aspect_ratio;
1691 float3 sp;
1692 sp.y = -0.5f + ((float)ii + Rx) / (float)params.camera_resolution_full.x;
1693 sp.z = ( 0.5f - ((float)jj + Ry) / (float)params.camera_resolution_full.y) * multiplier;
1694 sp.x = params.camera_viewplane_length;
1695
1696 // Focal point on focus plane
1697 const float3 p = make_float3(
1698 params.camera_focal_length,
1699 sp.y / params.camera_viewplane_length * params.camera_focal_length,
1700 sp.z / params.camera_viewplane_length * params.camera_focal_length);
1701
1702 // Sample lens (pinhole if lens_diameter == 0)
1703 float3 ray_origin = make_float3(0.f, 0.f, 0.f);
1704 if (params.camera_lens_diameter > 0.f) {
1705 float3 disk_sample;
1706 d_sampleDisk(seed, disk_sample);
1707 ray_origin = make_float3(0.f, 0.5f * disk_sample.x * params.camera_lens_diameter,
1708 0.5f * disk_sample.y * params.camera_lens_diameter);
1709 }
1710
1711 float3 ray_direction = make_float3(p.x - ray_origin.x, p.y - ray_origin.y, p.z - ray_origin.z);
1712
1713 // Rotate into world space
1714 const float theta = -0.5f * M_PI + params.camera_direction.x;
1715 const float phi = 0.5f * M_PI - params.camera_direction.y;
1716 ray_origin = d_rotatePoint(ray_origin, theta, phi) + params.camera_position;
1717 ray_direction = d_rotatePoint(ray_direction, theta, phi);
1718 ray_direction = ray_direction * (1.0f / d_magnitude(ray_direction));
1719
1720 PerRayData prd;
1721 prd.strength = 1.0f / (float)dim_x;
1722 prd.origin_UUID = pixel_index;
1723 prd.face = true;
1724 prd.source_ID = 0;
1725 prd.seed = seed;
1726 prd.hit_periodic_boundary = false;
1727
1728 uint32_t p0, p1;
1729 packPointer(&prd, p0, p1);
1730
1731 const float t_min = 1e-5f;
1732 const float t_max = 1e30f;
1733
1734 float3 cur_origin = ray_origin;
1735 for (int wrap = 0; wrap < 10; wrap++) {
1736 prd.hit_periodic_boundary = false;
1737 optixTrace(params.traversable, cur_origin, ray_direction,
1738 t_min, t_max, 0.f,
1739 OptixVisibilityMask(255), OPTIX_RAY_FLAG_NONE,
1740 2u, 0u, 2u, // SBT: offset=2 (camera hit), stride=0, miss=2
1741 p0, p1);
1742 if (!prd.hit_periodic_boundary) break;
1743 cur_origin = prd.periodic_hit;
1744 }
1745}
1746
1747// ---------------------------------------------------------------------------
1748// Raygen: pixel label rays
1749// 3D launch: x=1 (always), y=tile_column, z=tile_row
1750// ---------------------------------------------------------------------------
1751
1752extern "C" __global__ void __raygen__pixel_label() {
1753 const uint3 idx = optixGetLaunchIndex();
1754 const uint32_t col = idx.y;
1755 const uint32_t row = idx.z;
1756
1757 const uint32_t dim_y = params.launch_dim_y; // tile_width
1758
1759 // Global pixel coordinates
1760 const uint32_t ii = (uint32_t)params.camera_pixel_offset.x + col;
1761 const uint32_t jj = (uint32_t)params.camera_pixel_offset.y + row;
1762 const uint32_t pixel_index = jj * (uint32_t)params.camera_resolution_full.x + ii;
1763
1764 uint32_t seed = tea<16>(dim_y * row + col, params.random_seed);
1765
1766 // Center of pixel, no antialiasing jitter
1767 const float multiplier = 1.0f / params.FOV_aspect_ratio;
1768 float3 sp;
1769 sp.y = -0.5f + ((float)ii + 0.5f) / (float)params.camera_resolution_full.x;
1770 sp.z = ( 0.5f - ((float)jj + 0.5f) / (float)params.camera_resolution_full.y) * multiplier;
1771 sp.x = params.camera_viewplane_length;
1772
1773 const float3 p = make_float3(
1774 params.camera_focal_length,
1775 sp.y / params.camera_viewplane_length * params.camera_focal_length,
1776 sp.z / params.camera_viewplane_length * params.camera_focal_length);
1777
1778 float3 ray_origin = make_float3(0.f, 0.f, 0.f);
1779 float3 ray_direction = p;
1780
1781 const float theta = -0.5f * M_PI + params.camera_direction.x;
1782 const float phi = 0.5f * M_PI - params.camera_direction.y;
1783 ray_origin = d_rotatePoint(ray_origin, theta, phi) + params.camera_position;
1784 ray_direction = d_rotatePoint(ray_direction, theta, phi);
1785 ray_direction = ray_direction * (1.0f / d_magnitude(ray_direction));
1786
1787 PerRayData prd;
1788 prd.strength = 0.f; // used as distance offset (always 0 for pixel label)
1789 prd.origin_UUID = pixel_index;
1790 prd.face = true;
1791 prd.source_ID = 0;
1792 prd.seed = seed;
1793 prd.hit_periodic_boundary = false;
1794
1795 uint32_t p0, p1;
1796 packPointer(&prd, p0, p1);
1797
1798 const float t_min = 1e-5f;
1799 const float t_max = 1e30f;
1800
1801 float3 cur_origin = ray_origin;
1802 for (int wrap = 0; wrap < 10; wrap++) {
1803 prd.hit_periodic_boundary = false;
1804 optixTrace(params.traversable, cur_origin, ray_direction,
1805 t_min, t_max, 0.f,
1806 OptixVisibilityMask(255), OPTIX_RAY_FLAG_NONE,
1807 3u, 0u, 3u, // SBT: offset=3 (pixel label hit), stride=0, miss=3
1808 p0, p1);
1809 if (!prd.hit_periodic_boundary) break;
1810 cur_origin = prd.periodic_hit;
1811 }
1812}