1.3.77
 
Loading...
Searching...
No Matches
rayGeneration.cu
Go to the documentation of this file.
1
16#include <optix_world.h>
17#include "RayTracing.cuh"
18
19using namespace optix;
20
21RT_PROGRAM void direct_raygen() {
22
23 uint Nrays = launch_dim.x * launch_dim.y;
24 uint ray_index = launch_dim.x * launch_index.y + launch_index.x;
25
26 PerRayData prd;
27 prd.seed = tea<16>(ray_index + Nrays * launch_index.z, random_seed);
28
29 uint objID = launch_offset + launch_index.z;
30
31 uint pID = primitiveID[objID];
32 uint ptype = primitive_type[objID];
33 int puvID = uvID[objID];
34
35 float3 sp;
36
37 float3 normal;
38
39 // transformation matrix
40 float m_trans[16];
41 for (uint i = 0; i < 16; i++) {
42 m_trans[i] = transform_matrix[optix::make_uint2(i, objID)];
43 }
44
45 // looping over sub-patches
46 int NX = object_subdivisions[objID].x;
47 int NY = object_subdivisions[objID].y;
48 for (int jj = 0; jj < NY; jj++) {
49 for (int ii = 0; ii < NX; ii++) {
50
51 uint UUID = pID + jj * NX + ii;
52
53 // two random samples [0,1]
54 float Rx = rnd(prd.seed);
55 float Ry = rnd(prd.seed);
56
57 if (ptype == 0 || ptype == 3) { // Patch or Tile
58
59 uint Nx = launch_dim.x;
60 uint Ny = launch_dim.y;
61 float dx = 1.f / float(NX);
62 float dy = 1.f / float(NY);
63
64 // Map sample to rectangle [-0.5,0.5] [-0.5,0.5]
65 sp.x = -0.5f + ii * dx + float(launch_index.x) * dx / float(Nx) + Rx * dx / float(Nx);
66 sp.y = -0.5f + jj * dy + float(launch_index.y) * dy / float(Ny) + Ry * dy / float(Ny);
67 sp.z = 0.f;
68
69 int ID = maskID[objID];
70 // FIX: Use objID (position) instead of UUID for primitive_solid_fraction access
71 if (ID >= 0 && primitive_solid_fraction[objID] > 0.f && primitive_solid_fraction[objID] < 1.f) { // has texture transparency
72
73 d_sampleTexture_patch(sp, optix::make_int2(ii, jj), optix::make_float2(dx, dy), prd, ID, puvID);
74 }
75
76 // calculate rectangle normal vector (world coordinates)
77 float3 v0 = make_float3(0, 0, 0);
78 d_transformPoint(m_trans, v0);
79 float3 v1 = make_float3(1, 0, 0);
80 d_transformPoint(m_trans, v1);
81 float3 v2 = make_float3(0, 1, 0);
82 d_transformPoint(m_trans, v2);
83
84 normal = normalize(cross(v1 - v0, v2 - v0));
85
86 } else if (ptype == 1) { // Triangle
87
88 // Map sample to triangle with vertices (0,0,0), (0,1,0), (1,1,0)
89 if (Rx < Ry) {
90 sp.x = Rx;
91 sp.y = Ry;
92 } else {
93 sp.x = Ry;
94 sp.y = Rx;
95 }
96 sp.z = 0;
97
98 // calculate triangle normal vector (world coordinates)
99 float3 v0 = make_float3(0, 0, 0);
100 d_transformPoint(m_trans, v0);
101 float3 v1 = make_float3(0, 1, 0);
102 d_transformPoint(m_trans, v1);
103 float3 v2 = make_float3(1, 1, 0);
104 d_transformPoint(m_trans, v2);
105
106 normal = normalize(cross(v1 - v0, v2 - v0));
107
108 int ID = maskID[objID];
109 // FIX: Use objID (position) instead of UUID for primitive_solid_fraction access
110 if (ID >= 0 && primitive_solid_fraction[objID] > 0.f && primitive_solid_fraction[objID] < 1.f) { // has texture transparency
111
112 d_sampleTexture_triangle(sp, v0, v1, v2, prd, m_trans, ID, puvID);
113 }
114
115 } else if (ptype == 2) { // Disk
116
117 d_sampleDisk(prd.seed, sp);
118
119 // calculate disk normal vector (world coordinates)
120 float3 v0 = make_float3(0, 0, 0);
121 d_transformPoint(m_trans, v0);
122 float3 v1 = make_float3(1, 0, 0);
123 d_transformPoint(m_trans, v1);
124 float3 v2 = make_float3(0, 1, 0);
125 d_transformPoint(m_trans, v2);
126
127 normal = normalize(cross(v1 - v0, v2 - v0));
128
129 } else if (ptype == 4) { // Voxel
130
131 float Rz = rnd(prd.seed);
132 }
133
134 // translate the ray to the location of the primitive
135
136 float3 ray_origin = sp;
137 d_transformPoint(m_trans, ray_origin);
138
139 // Send a ray toward each source
140 for (int rr = 0; rr < Nsources; rr++) {
141
142 // set the ray direction
143 float3 ray_direction;
144 float ray_magnitude;
145 if (source_types[rr] == 0) { // collimated source
146 ray_direction = normalize(source_positions[rr]);
147 ray_magnitude = RT_DEFAULT_MAX;
148 prd.strength = 1. / double(launch_dim.x * launch_dim.y) * fabs(dot(normal, ray_direction));
149 } else if (source_types[rr] == 1 || source_types[rr] == 2) { // sphere source
150
151 // sample point on surface of sphere
152 float theta_s = acos_safe(1.f - 2.f * rnd(prd.seed));
153 float phi_s = rnd(prd.seed) * 2.f * M_PI;
154 float3 sphere_point = 0.5 * source_widths[rr].x * make_float3(sin(theta_s) * cos(phi_s), sin(theta_s) * sin(phi_s), cos(theta_s));
155
156 ray_direction = sphere_point + source_positions[rr] - ray_origin;
157
158 ray_magnitude = d_magnitude(ray_direction);
159 ray_direction = normalize(ray_direction);
160 prd.strength = 0.f;
161 uint N = 10;
162 for (uint j = 0; j < N; j++) {
163 for (uint i = 0; i < N; i++) {
164 float theta = acos_safe(1.f - 2.f * (float(i) + 0.5f) / float(N));
165 float phi = (float(j) + 0.5f) * 2.f * M_PI / float(N);
166 float3 light_direction = make_float3(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
167 if (dot(light_direction, ray_direction) < 0) {
168 prd.strength +=
169 1. / double(launch_dim.x * launch_dim.y) * fabs(dot(normal, ray_direction)) * fabs(dot(light_direction, ray_direction)) / (ray_magnitude * ray_magnitude) / (N * N) * source_widths[rr].x * source_widths[rr].x;
170 }
171 }
172 }
173
174 } else if (source_types[rr] == 3) { // rectangle source
175
176 // transformation matrix
177 float light_transform[16];
178 d_makeTransformMatrix(source_rotations[rr], light_transform);
179
180 // sample point on surface of disk
181 float3 square_point;
182 d_sampleSquare(prd.seed, square_point);
183 square_point = make_float3(source_widths[rr].x * square_point.x, source_widths[rr].y * square_point.y, square_point.z);
184 d_transformPoint(light_transform, square_point);
185
186 float3 light_direction = make_float3(0, 0, 1);
187 d_transformPoint(light_transform, light_direction);
188
189 ray_direction = square_point + source_positions[rr] - ray_origin;
190
191 if (dot(ray_direction, light_direction) > 0.f) { // don't emit from back side of light source (note that ray goes toward the source, so the dot produce is negative when light is pointed at primitive)
192 continue;
193 }
194
195 ray_magnitude = d_magnitude(ray_direction);
196 ray_direction = normalize(ray_direction);
197 prd.strength = 1. / double(launch_dim.x * launch_dim.y) * fabs(dot(normal, ray_direction)) * fabs(dot(light_direction, ray_direction)) / (ray_magnitude * ray_magnitude) * source_widths[rr].x * source_widths[rr].y / M_PI;
198
199 } else if (source_types[rr] == 4) { // disk source
200
201 // transformation matrix
202 float light_transform[16];
203 d_makeTransformMatrix(source_rotations[rr], light_transform);
204
205 // sample point on surface of disk
206 float3 disk_point;
207 d_sampleDisk(prd.seed, disk_point);
208 d_transformPoint(light_transform, disk_point);
209
210 float3 light_direction = make_float3(0, 0, 1);
211 d_transformPoint(light_transform, light_direction);
212
213 ray_direction = source_widths[rr].x * disk_point + source_positions[rr] - ray_origin;
214
215 if (dot(ray_direction, light_direction) > 0.f) { // don't emit from back side of light source (note that ray goes toward the source, so the dot produce is negative when light is pointed at primitive)
216 continue;
217 }
218
219 ray_magnitude = d_magnitude(ray_direction);
220 ray_direction = normalize(ray_direction);
221 prd.strength = 1. / double(launch_dim.x * launch_dim.y) * fabs(dot(normal, ray_direction)) * fabs(dot(light_direction, ray_direction)) / (ray_magnitude * ray_magnitude) * source_widths[rr].x * source_widths[rr].x;
222 }
223
224 optix::Ray ray = optix::make_Ray(ray_origin, ray_direction, direct_ray_type, 1e-4, ray_magnitude);
225
226 prd.origin_UUID = UUID;
227 prd.source_ID = rr;
228 prd.hit_periodic_boundary = false;
229 initCoverTransmittance(prd); // translucent-cover attenuation starts at 1 (no covers crossed)
230
231 if (dot(ray_direction, normal) > 0) {
232 prd.face = 1;
233 } else {
234 prd.face = 0;
235 }
236
237 if ((prd.face == 1 || twosided_flag[objID] == 1) && twosided_flag[objID] != 3) {
238
239 for (int wrap = 0; wrap < 10; ++wrap) {
240 rtTrace(top_object, ray, prd);
241
242 if (!prd.hit_periodic_boundary)
243 break; // real hit or miss → done
244
245 ray.origin = prd.periodic_hit;
246 prd.hit_periodic_boundary = false;
247 }
248 }
249 }
250 }
251 }
252}
253
254RT_PROGRAM void diffuse_raygen() {
255
256 uint dimx = launch_dim.x * launch_dim.y;
257 uint indx = launch_dim.x * launch_index.y + launch_index.x;
258
259 PerRayData prd;
260 prd.seed = tea<16>(indx + dimx * launch_index.z, random_seed);
261
262 uint objID = launch_offset + launch_index.z;
263
264 if (launch_face == 0 && twosided_flag[objID] == 0) { // skip the launch if from the bottom face and twosided_flag = 0
265 return;
266 }
267
268 uint pID = primitiveID[objID];
269 uint ptype = primitive_type[objID];
270 int puvID = uvID[objID];
271
272 // transformation matrix
273 float m_trans[16];
274 for (uint i = 0; i < 16; i++) {
275 m_trans[i] = transform_matrix[optix::make_uint2(i, objID)];
276 }
277
278 float3 sp, normal;
279
280 // looping over sub-patches
281 int NX = object_subdivisions[objID].x;
282 int NY = object_subdivisions[objID].y;
283 for (int jj = 0; jj < NY; jj++) {
284 for (int ii = 0; ii < NX; ii++) {
285
286 uint UUID = pID + jj * NX + ii;
287
288 // two random samples [0,1]
289 float Rx = rnd(prd.seed);
290 float Ry = rnd(prd.seed);
291
292 if (ptype == 0 || ptype == 3) { // Patch or Tile
293
294 // calculate rectangle normal vector (world coordinates)
295 float3 s0 = make_float3(0, 0, 0);
296 float3 s1 = make_float3(1, 0, 0);
297 float3 s2 = make_float3(0, 1, 0);
298 d_transformPoint(m_trans, s0);
299 d_transformPoint(m_trans, s1);
300 d_transformPoint(m_trans, s2);
301
302 normal = normalize(cross(s1 - s0, s2 - s0));
303
304 float dx = 1.f / float(NX);
305 float dy = 1.f / float(NY);
306
307 // Map sample to rectangle [-0.5,0.5] [-0.5,0.5]
308 sp.x = -0.5f + (ii + Rx) * dx;
309 sp.y = -0.5f + (jj + Ry) * dy;
310 sp.z = 0.f;
311
312 int ID = maskID[objID];
313 if (ID >= 0) { // has texture transparency
314
315 d_sampleTexture_patch(sp, optix::make_int2(ii, jj), optix::make_float2(dx, dy), prd, ID, puvID);
316 }
317
318 } else if (ptype == 1) { // Triangle
319
320 // Map sample to triangle with vertices (0,0,0), (0,1,0), (1,1,0)
321 if (Rx < Ry) {
322 sp.x = Rx;
323 sp.y = Ry;
324 } else {
325 sp.x = Ry;
326 sp.y = Rx;
327 }
328 sp.z = 0;
329
330 // calculate triangle normal vector (world coordinates)
331 float3 v0 = make_float3(0, 0, 0);
332 d_transformPoint(m_trans, v0);
333 float3 v1 = make_float3(0, 1, 0);
334 d_transformPoint(m_trans, v1);
335 float3 v2 = make_float3(1, 1, 0);
336 d_transformPoint(m_trans, v2);
337
338 normal = normalize(cross(v1 - v0, v2 - v0));
339
340 int ID = maskID[objID];
341 if (ID >= 0) { // has texture transparency
342
343 d_sampleTexture_triangle(sp, v0, v1, v2, prd, m_trans, ID, puvID);
344 }
345
346 } else if (ptype == 2) { // Disk
347
348 // Map Sample to disk - from Suffern (2007) "Ray tracing fom the ground up" Chap. 6
349
350 // first map sample point to rectangle [-1,1] [-1,1]
351 sp.x = -1.f + 2.f * Rx;
352 sp.y = -1.f + 2.f * Ry;
353
354 float r, p;
355 if (sp.x > -sp.y) {
356 if (sp.x > sp.y) {
357 r = sp.x;
358 p = sp.y / sp.x;
359 } else {
360 r = sp.y;
361 p = 2.f - sp.x / sp.y;
362 }
363 } else {
364 if (sp.x < sp.y) {
365 r = -sp.x;
366 p = 4.f + sp.y / sp.x;
367 } else {
368 r = -sp.y;
369 if (sp.y != 0.f) { // avoid division by zero at origin
370 p = 6.f - sp.x / sp.y;
371 } else {
372 p = 0.f;
373 }
374 }
375 }
376 p *= 0.25f * M_PI;
377
378 // find x,y point on unit disk
379 sp.x = r * cosf(p);
380 sp.y = r * sinf(p);
381 sp.z = 0.f;
382
383 // calculate disk normal vector (world coordinates)
384 float3 v0 = make_float3(0, 0, 0);
385 d_transformPoint(m_trans, v0);
386 float3 v1 = make_float3(1, 0, 0);
387 d_transformPoint(m_trans, v1);
388 float3 v2 = make_float3(0, 1, 0);
389 d_transformPoint(m_trans, v2);
390 normal = normalize(cross(v1 - v0, v2 - v0));
391
392 } else if (ptype == 4) { // Voxel
393
394 // Map sample to cube [-0.5,0.5] [-0.5,0.5] [-0.5,0.5]
395 sp.x = -0.5f + Rx;
396 sp.y = -0.5f + Ry;
397 sp.z = -0.5f + rnd(prd.seed);
398 }
399
400 // Choose random hemispherical direction - map samples to hemisphere (from Suffern (2007) "Ray tracing fom the ground up" Chap. 6)
401
402 float Rt;
403 float Rp;
404
405 Rt = (launch_index.x + rnd(prd.seed)) / float(launch_dim.x);
406 Rp = (launch_index.y + rnd(prd.seed)) / float(launch_dim.y);
407
408 float t;
409 if (ptype == 4) { // voxel
410 t = acos_safe(1.f - Rt);
411 } else { // other
412 t = asin_safe(sqrtf(Rt));
413 }
414 float p = 2.f * M_PI * Rp;
415
416 float3 ray_direction;
417 ray_direction.x = sin(t) * cos(p);
418 ray_direction.y = sin(t) * sin(p);
419 ray_direction.z = cos(t);
420
421 float3 ray_origin;
422 optix::Ray ray;
423
424 if (ptype == 4) { // voxel
425
426 prd.strength = 0.5f / float(dimx);
427 prd.origin_UUID = UUID;
428 prd.face = 0;
429 prd.source_ID = 0;
430 prd.hit_periodic_boundary = false;
431
432 ray_origin = sp;
433 d_transformPoint(m_trans, ray_origin);
434
435 ray = optix::make_Ray(ray_origin, ray_direction, diffuse_ray_type, 1e-5, RT_DEFAULT_MAX);
436 initCoverTransmittance(prd);
437 rtTrace(top_object, ray, prd);
438
439 ray = optix::make_Ray(ray_origin, -ray_direction, diffuse_ray_type, 1e-5, RT_DEFAULT_MAX);
440 initCoverTransmittance(prd);
441 rtTrace(top_object, ray, prd);
442
443 } else { // not a voxel
444
445 ray_direction = d_rotatePoint(ray_direction, acos_safe(normal.z), atan2(normal.y, normal.x));
446
447 prd.strength = 1.f / float(dimx);
448
449 prd.origin_UUID = UUID;
450 prd.source_ID = 0;
451 prd.hit_periodic_boundary = false;
452 initCoverTransmittance(prd); // translucent-cover attenuation starts at 1 (no covers crossed)
453
454 // ---- "top" surface launch -------
455 ray_origin = sp;
456 d_transformPoint(m_trans, ray_origin);
457
458 if (launch_face == 1 && twosided_flag[objID] != 3) {
459
460 ray = optix::make_Ray(ray_origin, ray_direction, diffuse_ray_type, 1e-5, RT_DEFAULT_MAX);
461
462 prd.face = 1;
463
464 for (int wrap = 0; wrap < 10; ++wrap) {
465 rtTrace(top_object, ray, prd);
466
467 if (!prd.hit_periodic_boundary)
468 break; // real hit or miss → done
469
470 ray.origin = prd.periodic_hit;
471 prd.hit_periodic_boundary = false;
472 }
473
474 // ---- "bottom" surface launch -------
475 } else if (launch_face == 0 && twosided_flag[objID] == 1) {
476
477 ray_direction = -ray_direction;
478 ray = optix::make_Ray(ray_origin, ray_direction, diffuse_ray_type, 1e-5, RT_DEFAULT_MAX);
479
480 prd.face = 0;
481
482 for (int wrap = 0; wrap < 10; ++wrap) {
483 rtTrace(top_object, ray, prd);
484
485 if (!prd.hit_periodic_boundary)
486 break; // real hit or miss → done
487
488 ray.origin = prd.periodic_hit;
489 prd.hit_periodic_boundary = false;
490 }
491 }
492 // else: Skip ray trace for bottom face of one-sided primitives
493 }
494 }
495 }
496}
497
498RT_PROGRAM void camera_raygen() {
499
500 uint dimx = launch_dim.x * launch_dim.y; // x number of ray, y width, z length
501 uint indx = launch_dim.x * launch_index.y + launch_index.x;
502
503 // Use full camera resolution (not tile size) for correct pixel calculations
504 optix::int2 camera_resolution = camera_resolution_full;
505
506 PerRayData prd;
507 prd.seed = tea<16>(indx + dimx * launch_index.z, random_seed);
508
509 float3 sp;
510
511 // Calculate global pixel coordinates including tile offsets
512 uint ii = camera_pixel_offset_x + launch_index.y; // global x-pixel
513 uint jj = camera_pixel_offset_y + launch_index.z; // global y-pixel
514 size_t origin_ID = jj * camera_resolution_full.x + ii; // global pixel index
515
516
517 // distortion
518 // float PPointsRatiox =1.052f;
519 // float PPointsRatioy =0.999f;
520 // float sensorxscale = 1.0054;
521 // float focalxy = 710;
522 // double x =(float(ii)-camera_resolution.x/2 * PPointsRatiox)/focalxy*sensorxscale;// / focalxy; cam_res.y = 712
523 // double y = (float(jj)-camera_resolution.y/2 * PPointsRatioy )/focalxy; /// focalxy; cam_res.x = 1072
524 // double r2 = x*x + y*y;
525 // double distCoeffs[4] = {-0.3535674,0.17298, 0, 0};
526 // double ii_d = x * (1+ distCoeffs[0] * r2 + distCoeffs[1] * r2 * r2) + 2 * distCoeffs[2] * x * y + distCoeffs[3] * (r2 + 2 * x * x);
527 // double jj_d = y * (1+ distCoeffs[0] * r2 + distCoeffs[1] * r2 * r2) + 2 * distCoeffs[3] * x * y + distCoeffs[2] * (r2 + 2 * y * y);
528 // ii_d = ii_d*focalxy+float(camera_resolution.x)/2 * PPointsRatiox;
529 // jj_d = jj_d*focalxy+float(camera_resolution.y)/2 * PPointsRatioy;
530
531 // *** sample a point on the pixel (view direction coordinate aligned) *** //
532
533 float PPointsRatiox = 1.f;
534 float PPointsRatioy = 1.f;
535 float Rx = rnd(prd.seed);
536 float Ry = rnd(prd.seed);
537
538 // Map sample to pixel
539 // Calculate VFOV scaling factor directly from aspect ratio
540 // Since FOV_aspect_ratio = HFOV/VFOV, then tan(half_VFOV) = tan(half_HFOV)/FOV_aspect_ratio
541 // The multiplier = tan(half_VFOV)/tan(half_HFOV) = 1/FOV_aspect_ratio
542 float multiplier = 1.0f / FOV_aspect_ratio;
543 sp.y = (-0.5f * PPointsRatioy + (ii + Rx) / float(camera_resolution.x));
544 sp.z = (0.5f * PPointsRatiox - (jj + Ry) / float(camera_resolution.y)) * multiplier;
545 sp.x = camera_viewplane_length;
546
547 // *** Determine point 'p' on focal plane that passes through the lens center (0,0) and pixel sample (view direction coordinate aligned) *** //
548 // Note: camera_focal_length is the focal plane distance (working distance), not the lens optical focal length
549
550 float3 p = make_float3(camera_focal_length, sp.y / camera_viewplane_length * camera_focal_length, sp.z / camera_viewplane_length * camera_focal_length);
551
552 // *** Sample point on lens (view direction coordinate aligned) *** //
553
554 float3 ray_origin = make_float3(0, 0, 0);
555 if (camera_lens_diameter > 0) {
556 float3 disk_sample;
557 d_sampleDisk(prd.seed, disk_sample);
558 ray_origin = make_float3(0.f, 0.5f * disk_sample.x * camera_lens_diameter, 0.5f * disk_sample.y * camera_lens_diameter);
559 }
560
561 //*** ray direction is line from lens sample to p ***//
562
563 float3 ray_direction = p - ray_origin;
564
565 //*** rotate ray origin and direction into the direction of the camera view *** //
566
567 ray_origin = d_rotatePoint(ray_origin, -0.5 * M_PI + camera_direction.x, 0.5f * M_PI - camera_direction.y) + camera_position;
568
569 ray_direction = d_rotatePoint(ray_direction, -0.5 * M_PI + camera_direction.x, 0.5f * M_PI - camera_direction.y);
570 ray_direction /= d_magnitude(ray_direction);
571
572 optix::Ray ray;
573
574 prd.strength = 1.f / float(launch_dim.x);
575
576 prd.origin_UUID = origin_ID;
577 prd.face = 1;
578 prd.source_ID = 0;
579 prd.hit_periodic_boundary = false;
580
581 ray = optix::make_Ray(ray_origin, ray_direction, camera_ray_type, 1e-5, RT_DEFAULT_MAX);
582
583 for (int wrap = 0; wrap < 10; ++wrap) {
584 rtTrace(top_object, ray, prd);
585
586 if (!prd.hit_periodic_boundary)
587 break; // real hit or miss → done
588
589 ray.origin = prd.periodic_hit;
590 prd.hit_periodic_boundary = false;
591 }
592}
593
594RT_PROGRAM void pixel_label_raygen() {
595
596 uint indx = launch_dim.y * launch_index.z + launch_index.y;
597
598 // Use full camera resolution (not tile size) for correct pixel calculations
599 optix::int2 camera_resolution = camera_resolution_full;
600
601 PerRayData prd;
602 prd.seed = tea<16>(indx, random_seed);
603
604 float3 sp;
605
606 // Calculate global pixel coordinates including tile offsets
607 uint ii = camera_pixel_offset_x + launch_index.y; // global x-pixel
608
609 uint jj = camera_pixel_offset_y + launch_index.z; // global y-pixel
610
611 size_t origin_ID = jj * camera_resolution_full.x + ii; // global pixel index
612
613 // Map sample to center of pixel
614 // Calculate VFOV scaling factor directly from aspect ratio
615 // Since FOV_aspect_ratio = HFOV/VFOV, then tan(half_VFOV) = tan(half_HFOV)/FOV_aspect_ratio
616 // The multiplier = tan(half_VFOV)/tan(half_HFOV) = 1/FOV_aspect_ratio
617 float multiplier = 1.0f / FOV_aspect_ratio;
618 sp.y = (-0.5f + (ii + 0.5f) / float(camera_resolution.x));
619 sp.z = (0.5f - (jj + 0.5f) / float(camera_resolution.y)) * multiplier;
620 sp.x = camera_viewplane_length;
621
622
623 // *** Determine point 'p' on focal plane that passes through the lens center (0,0) and pixel sample (view direction coordinate aligned) *** //
624 // Note: camera_focal_length is the focal plane distance (working distance), not the lens optical focal length
625
626 float3 p = make_float3(camera_focal_length, sp.y / camera_viewplane_length * camera_focal_length, sp.z / camera_viewplane_length * camera_focal_length);
627
628 // *** Ray is launched from center of lens *** //
629
630 float3 ray_origin = make_float3(0, 0, 0);
631
632 //*** ray direction is line from ray origin to p ***//
633
634 float3 ray_direction = p;
635
636 //*** rotate ray origin and direction into the direction of the camera view *** //
637
638 ray_origin = d_rotatePoint(ray_origin, -0.5 * M_PI + camera_direction.x, 0.5f * M_PI - camera_direction.y) + camera_position;
639
640 ray_direction = d_rotatePoint(ray_direction, -0.5 * M_PI + camera_direction.x, 0.5f * M_PI - camera_direction.y);
641 ray_direction /= d_magnitude(ray_direction);
642
643 optix::Ray ray;
644
645 prd.strength = 0.f;
646
647 prd.origin_UUID = origin_ID;
648 prd.face = 1;
649 prd.source_ID = 0;
650 prd.hit_periodic_boundary = false;
651
652 ray = optix::make_Ray(ray_origin, ray_direction, pixel_label_ray_type, 1e-5, RT_DEFAULT_MAX);
653
654 for (int wrap = 0; wrap < 10; ++wrap) {
655 rtTrace(top_object, ray, prd);
656
657 if (!prd.hit_periodic_boundary)
658 break; // real hit or miss → done
659
660 ray.origin = prd.periodic_hit;
661 prd.hit_periodic_boundary = false;
662 }
663}