1.3.64
 
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
230 if (dot(ray_direction, normal) > 0) {
231 prd.face = 1;
232 } else {
233 prd.face = 0;
234 }
235
236 if ((prd.face == 1 || twosided_flag[objID] == 1) && twosided_flag[objID] != 3) {
237
238 for (int wrap = 0; wrap < 10; ++wrap) {
239 rtTrace(top_object, ray, prd);
240
241 if (!prd.hit_periodic_boundary)
242 break; // real hit or miss → done
243
244 ray.origin = prd.periodic_hit;
245 prd.hit_periodic_boundary = false;
246 }
247 }
248 }
249 }
250 }
251}
252
253RT_PROGRAM void diffuse_raygen() {
254
255 uint dimx = launch_dim.x * launch_dim.y;
256 uint indx = launch_dim.x * launch_index.y + launch_index.x;
257
258 PerRayData prd;
259 prd.seed = tea<16>(indx + dimx * launch_index.z, random_seed);
260
261 uint objID = launch_offset + launch_index.z;
262
263 if (launch_face == 0 && twosided_flag[objID] == 0) { // skip the launch if from the bottom face and twosided_flag = 0
264 return;
265 }
266
267 uint pID = primitiveID[objID];
268 uint ptype = primitive_type[objID];
269 int puvID = uvID[objID];
270
271 // transformation matrix
272 float m_trans[16];
273 for (uint i = 0; i < 16; i++) {
274 m_trans[i] = transform_matrix[optix::make_uint2(i, objID)];
275 }
276
277 float3 sp, normal;
278
279 // looping over sub-patches
280 int NX = object_subdivisions[objID].x;
281 int NY = object_subdivisions[objID].y;
282 for (int jj = 0; jj < NY; jj++) {
283 for (int ii = 0; ii < NX; ii++) {
284
285 uint UUID = pID + jj * NX + ii;
286
287 // two random samples [0,1]
288 float Rx = rnd(prd.seed);
289 float Ry = rnd(prd.seed);
290
291 if (ptype == 0 || ptype == 3) { // Patch or Tile
292
293 // calculate rectangle normal vector (world coordinates)
294 float3 s0 = make_float3(0, 0, 0);
295 float3 s1 = make_float3(1, 0, 0);
296 float3 s2 = make_float3(0, 1, 0);
297 d_transformPoint(m_trans, s0);
298 d_transformPoint(m_trans, s1);
299 d_transformPoint(m_trans, s2);
300
301 normal = normalize(cross(s1 - s0, s2 - s0));
302
303 float dx = 1.f / float(NX);
304 float dy = 1.f / float(NY);
305
306 // Map sample to rectangle [-0.5,0.5] [-0.5,0.5]
307 sp.x = -0.5f + (ii + Rx) * dx;
308 sp.y = -0.5f + (jj + Ry) * dy;
309 sp.z = 0.f;
310
311 int ID = maskID[objID];
312 if (ID >= 0) { // has texture transparency
313
314 d_sampleTexture_patch(sp, optix::make_int2(ii, jj), optix::make_float2(dx, dy), prd, ID, puvID);
315 }
316
317 } else if (ptype == 1) { // Triangle
318
319 // Map sample to triangle with vertices (0,0,0), (0,1,0), (1,1,0)
320 if (Rx < Ry) {
321 sp.x = Rx;
322 sp.y = Ry;
323 } else {
324 sp.x = Ry;
325 sp.y = Rx;
326 }
327 sp.z = 0;
328
329 // calculate triangle normal vector (world coordinates)
330 float3 v0 = make_float3(0, 0, 0);
331 d_transformPoint(m_trans, v0);
332 float3 v1 = make_float3(0, 1, 0);
333 d_transformPoint(m_trans, v1);
334 float3 v2 = make_float3(1, 1, 0);
335 d_transformPoint(m_trans, v2);
336
337 normal = normalize(cross(v1 - v0, v2 - v0));
338
339 int ID = maskID[objID];
340 if (ID >= 0) { // has texture transparency
341
342 d_sampleTexture_triangle(sp, v0, v1, v2, prd, m_trans, ID, puvID);
343 }
344
345 } else if (ptype == 2) { // Disk
346
347 // Map Sample to disk - from Suffern (2007) "Ray tracing fom the ground up" Chap. 6
348
349 // first map sample point to rectangle [-1,1] [-1,1]
350 sp.x = -1.f + 2.f * Rx;
351 sp.y = -1.f + 2.f * Ry;
352
353 float r, p;
354 if (sp.x > -sp.y) {
355 if (sp.x > sp.y) {
356 r = sp.x;
357 p = sp.y / sp.x;
358 } else {
359 r = sp.y;
360 p = 2.f - sp.x / sp.y;
361 }
362 } else {
363 if (sp.x < sp.y) {
364 r = -sp.x;
365 p = 4.f + sp.y / sp.x;
366 } else {
367 r = -sp.y;
368 if (sp.y != 0.f) { // avoid division by zero at origin
369 p = 6.f - sp.x / sp.y;
370 } else {
371 p = 0.f;
372 }
373 }
374 }
375 p *= 0.25f * M_PI;
376
377 // find x,y point on unit disk
378 sp.x = r * cosf(p);
379 sp.y = r * sinf(p);
380 sp.z = 0.f;
381
382 // calculate disk normal vector (world coordinates)
383 float3 v0 = make_float3(0, 0, 0);
384 d_transformPoint(m_trans, v0);
385 float3 v1 = make_float3(1, 0, 0);
386 d_transformPoint(m_trans, v1);
387 float3 v2 = make_float3(0, 1, 0);
388 d_transformPoint(m_trans, v2);
389 normal = normalize(cross(v1 - v0, v2 - v0));
390
391 } else if (ptype == 4) { // Voxel
392
393 // Map sample to cube [-0.5,0.5] [-0.5,0.5] [-0.5,0.5]
394 sp.x = -0.5f + Rx;
395 sp.y = -0.5f + Ry;
396 sp.z = -0.5f + rnd(prd.seed);
397 }
398
399 // Choose random hemispherical direction - map samples to hemisphere (from Suffern (2007) "Ray tracing fom the ground up" Chap. 6)
400
401 float Rt;
402 float Rp;
403
404 Rt = (launch_index.x + rnd(prd.seed)) / float(launch_dim.x);
405 Rp = (launch_index.y + rnd(prd.seed)) / float(launch_dim.y);
406
407 float t;
408 if (ptype == 4) { // voxel
409 t = acos_safe(1.f - Rt);
410 } else { // other
411 t = asin_safe(sqrtf(Rt));
412 }
413 float p = 2.f * M_PI * Rp;
414
415 float3 ray_direction;
416 ray_direction.x = sin(t) * cos(p);
417 ray_direction.y = sin(t) * sin(p);
418 ray_direction.z = cos(t);
419
420 float3 ray_origin;
421 optix::Ray ray;
422
423 if (ptype == 4) { // voxel
424
425 prd.strength = 0.5f / float(dimx);
426 prd.origin_UUID = UUID;
427 prd.face = 0;
428 prd.source_ID = 0;
429 prd.hit_periodic_boundary = false;
430
431 ray_origin = sp;
432 d_transformPoint(m_trans, ray_origin);
433
434 ray = optix::make_Ray(ray_origin, ray_direction, diffuse_ray_type, 1e-5, RT_DEFAULT_MAX);
435 rtTrace(top_object, ray, prd);
436
437 ray = optix::make_Ray(ray_origin, -ray_direction, diffuse_ray_type, 1e-5, RT_DEFAULT_MAX);
438 rtTrace(top_object, ray, prd);
439
440 } else { // not a voxel
441
442 ray_direction = d_rotatePoint(ray_direction, acos_safe(normal.z), atan2(normal.y, normal.x));
443
444 prd.strength = 1.f / float(dimx);
445
446 prd.origin_UUID = UUID;
447 prd.source_ID = 0;
448 prd.hit_periodic_boundary = false;
449
450 // ---- "top" surface launch -------
451 ray_origin = sp;
452 d_transformPoint(m_trans, ray_origin);
453
454 if (launch_face == 1 && twosided_flag[objID] != 3) {
455
456 ray = optix::make_Ray(ray_origin, ray_direction, diffuse_ray_type, 1e-5, RT_DEFAULT_MAX);
457
458 prd.face = 1;
459
460 for (int wrap = 0; wrap < 10; ++wrap) {
461 rtTrace(top_object, ray, prd);
462
463 if (!prd.hit_periodic_boundary)
464 break; // real hit or miss → done
465
466 ray.origin = prd.periodic_hit;
467 prd.hit_periodic_boundary = false;
468 }
469
470 // ---- "bottom" surface launch -------
471 } else if (launch_face == 0 && twosided_flag[objID] == 1) {
472
473 ray_direction = -ray_direction;
474 ray = optix::make_Ray(ray_origin, ray_direction, diffuse_ray_type, 1e-5, RT_DEFAULT_MAX);
475
476 prd.face = 0;
477
478 for (int wrap = 0; wrap < 10; ++wrap) {
479 rtTrace(top_object, ray, prd);
480
481 if (!prd.hit_periodic_boundary)
482 break; // real hit or miss → done
483
484 ray.origin = prd.periodic_hit;
485 prd.hit_periodic_boundary = false;
486 }
487 }
488 // else: Skip ray trace for bottom face of one-sided primitives
489 }
490 }
491 }
492}
493
494RT_PROGRAM void camera_raygen() {
495
496 uint dimx = launch_dim.x * launch_dim.y; // x number of ray, y width, z length
497 uint indx = launch_dim.x * launch_index.y + launch_index.x;
498
499 // Use full camera resolution (not tile size) for correct pixel calculations
500 optix::int2 camera_resolution = camera_resolution_full;
501
502 PerRayData prd;
503 prd.seed = tea<16>(indx + dimx * launch_index.z, random_seed);
504
505 float3 sp;
506
507 // Calculate global pixel coordinates including tile offsets
508 uint ii = camera_pixel_offset_x + launch_index.y; // global x-pixel
509 uint jj = camera_pixel_offset_y + launch_index.z; // global y-pixel
510 size_t origin_ID = jj * camera_resolution_full.x + ii; // global pixel index
511
512
513 // distortion
514 // float PPointsRatiox =1.052f;
515 // float PPointsRatioy =0.999f;
516 // float sensorxscale = 1.0054;
517 // float focalxy = 710;
518 // double x =(float(ii)-camera_resolution.x/2 * PPointsRatiox)/focalxy*sensorxscale;// / focalxy; cam_res.y = 712
519 // double y = (float(jj)-camera_resolution.y/2 * PPointsRatioy )/focalxy; /// focalxy; cam_res.x = 1072
520 // double r2 = x*x + y*y;
521 // double distCoeffs[4] = {-0.3535674,0.17298, 0, 0};
522 // double ii_d = x * (1+ distCoeffs[0] * r2 + distCoeffs[1] * r2 * r2) + 2 * distCoeffs[2] * x * y + distCoeffs[3] * (r2 + 2 * x * x);
523 // double jj_d = y * (1+ distCoeffs[0] * r2 + distCoeffs[1] * r2 * r2) + 2 * distCoeffs[3] * x * y + distCoeffs[2] * (r2 + 2 * y * y);
524 // ii_d = ii_d*focalxy+float(camera_resolution.x)/2 * PPointsRatiox;
525 // jj_d = jj_d*focalxy+float(camera_resolution.y)/2 * PPointsRatioy;
526
527 // *** sample a point on the pixel (view direction coordinate aligned) *** //
528
529 float PPointsRatiox = 1.f;
530 float PPointsRatioy = 1.f;
531 float Rx = rnd(prd.seed);
532 float Ry = rnd(prd.seed);
533
534 // Map sample to pixel
535 // Calculate VFOV scaling factor directly from aspect ratio
536 // Since FOV_aspect_ratio = HFOV/VFOV, then tan(half_VFOV) = tan(half_HFOV)/FOV_aspect_ratio
537 // The multiplier = tan(half_VFOV)/tan(half_HFOV) = 1/FOV_aspect_ratio
538 float multiplier = 1.0f / FOV_aspect_ratio;
539 sp.y = (-0.5f * PPointsRatioy + (ii + Rx) / float(camera_resolution.x));
540 sp.z = (0.5f * PPointsRatiox - (jj + Ry) / float(camera_resolution.y)) * multiplier;
541 sp.x = camera_viewplane_length;
542
543 // *** Determine point 'p' on focal plane that passes through the lens center (0,0) and pixel sample (view direction coordinate aligned) *** //
544 // Note: camera_focal_length is the focal plane distance (working distance), not the lens optical focal length
545
546 float3 p = make_float3(camera_focal_length, sp.y / camera_viewplane_length * camera_focal_length, sp.z / camera_viewplane_length * camera_focal_length);
547
548 // *** Sample point on lens (view direction coordinate aligned) *** //
549
550 float3 ray_origin = make_float3(0, 0, 0);
551 if (camera_lens_diameter > 0) {
552 float3 disk_sample;
553 d_sampleDisk(prd.seed, disk_sample);
554 ray_origin = make_float3(0.f, 0.5f * disk_sample.x * camera_lens_diameter, 0.5f * disk_sample.y * camera_lens_diameter);
555 }
556
557 //*** ray direction is line from lens sample to p ***//
558
559 float3 ray_direction = p - ray_origin;
560
561 //*** rotate ray origin and direction into the direction of the camera view *** //
562
563 ray_origin = d_rotatePoint(ray_origin, -0.5 * M_PI + camera_direction.x, 0.5f * M_PI - camera_direction.y) + camera_position;
564
565 ray_direction = d_rotatePoint(ray_direction, -0.5 * M_PI + camera_direction.x, 0.5f * M_PI - camera_direction.y);
566 ray_direction /= d_magnitude(ray_direction);
567
568 optix::Ray ray;
569
570 prd.strength = 1.f / float(launch_dim.x);
571
572 prd.origin_UUID = origin_ID;
573 prd.face = 1;
574 prd.source_ID = 0;
575 prd.hit_periodic_boundary = false;
576
577 ray = optix::make_Ray(ray_origin, ray_direction, camera_ray_type, 1e-5, RT_DEFAULT_MAX);
578
579 for (int wrap = 0; wrap < 10; ++wrap) {
580 rtTrace(top_object, ray, prd);
581
582 if (!prd.hit_periodic_boundary)
583 break; // real hit or miss → done
584
585 ray.origin = prd.periodic_hit;
586 prd.hit_periodic_boundary = false;
587 }
588}
589
590RT_PROGRAM void pixel_label_raygen() {
591
592 uint indx = launch_dim.y * launch_index.z + launch_index.y;
593
594 // Use full camera resolution (not tile size) for correct pixel calculations
595 optix::int2 camera_resolution = camera_resolution_full;
596
597 PerRayData prd;
598 prd.seed = tea<16>(indx, random_seed);
599
600 float3 sp;
601
602 // Calculate global pixel coordinates including tile offsets
603 uint ii = camera_pixel_offset_x + launch_index.y; // global x-pixel
604
605 uint jj = camera_pixel_offset_y + launch_index.z; // global y-pixel
606
607 size_t origin_ID = jj * camera_resolution_full.x + ii; // global pixel index
608
609 // Map sample to center of pixel
610 // Calculate VFOV scaling factor directly from aspect ratio
611 // Since FOV_aspect_ratio = HFOV/VFOV, then tan(half_VFOV) = tan(half_HFOV)/FOV_aspect_ratio
612 // The multiplier = tan(half_VFOV)/tan(half_HFOV) = 1/FOV_aspect_ratio
613 float multiplier = 1.0f / FOV_aspect_ratio;
614 sp.y = (-0.5f + (ii + 0.5f) / float(camera_resolution.x));
615 sp.z = (0.5f - (jj + 0.5f) / float(camera_resolution.y)) * multiplier;
616 sp.x = camera_viewplane_length;
617
618
619 // *** Determine point 'p' on focal plane that passes through the lens center (0,0) and pixel sample (view direction coordinate aligned) *** //
620 // Note: camera_focal_length is the focal plane distance (working distance), not the lens optical focal length
621
622 float3 p = make_float3(camera_focal_length, sp.y / camera_viewplane_length * camera_focal_length, sp.z / camera_viewplane_length * camera_focal_length);
623
624 // *** Ray is launched from center of lens *** //
625
626 float3 ray_origin = make_float3(0, 0, 0);
627
628 //*** ray direction is line from ray origin to p ***//
629
630 float3 ray_direction = p;
631
632 //*** rotate ray origin and direction into the direction of the camera view *** //
633
634 ray_origin = d_rotatePoint(ray_origin, -0.5 * M_PI + camera_direction.x, 0.5f * M_PI - camera_direction.y) + camera_position;
635
636 ray_direction = d_rotatePoint(ray_direction, -0.5 * M_PI + camera_direction.x, 0.5f * M_PI - camera_direction.y);
637 ray_direction /= d_magnitude(ray_direction);
638
639 optix::Ray ray;
640
641 prd.strength = 0.f;
642
643 prd.origin_UUID = origin_ID;
644 prd.face = 1;
645 prd.source_ID = 0;
646 prd.hit_periodic_boundary = false;
647
648 ray = optix::make_Ray(ray_origin, ray_direction, pixel_label_ray_type, 1e-5, RT_DEFAULT_MAX);
649
650 for (int wrap = 0; wrap < 10; ++wrap) {
651 rtTrace(top_object, ray, prd);
652
653 if (!prd.hit_periodic_boundary)
654 break; // real hit or miss → done
655
656 ray.origin = prd.periodic_hit;
657 prd.hit_periodic_boundary = false;
658 }
659}