sphor_compute_mvrea.c (26617B)
1 /* Copyright (C) 2024-2026 Centre National de la Recherche Scientifique 2 * Copyright (C) 2024-2026 Clermont Auvergne INP 3 * Copyright (C) 2024-2026 INSA Lyon 4 * Copyright (C) 2024-2026 Institut Mines Télécom Albi-Carmaux 5 * Copyright (C) 2024-2026 Institut National Polytechnique de Toulouse 6 * Copyright (C) 2024-2026 |Méso|Star> (contact@meso-star.com) 7 * Copyright (C) 2024-2026 PhotonLyX (info@photonlyx.com) 8 * Copyright (C) 2024-2026 Université de Lorraine 9 * Copyright (C) 2024-2026 Université Paul Sabatier 10 * Copyright (C) 2024-2026 Université Toulouse - Jean Jaurès 11 * 12 * This program is free software: you can redistribute it and/or modify 13 * it under the terms of the GNU General Public License as published by 14 * the Free Software Foundation, either version 3 of the License, or 15 * (at your option) any later version. 16 * 17 * This program is distributed in the hope that it will be useful, 18 * but WITHOUT ANY WARRANTY; without even the implied warranty of 19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 20 * GNU General Public License for more details. 21 * 22 * You should have received a copy of the GNU General Public License 23 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 24 25 #include "sphor_c.h" 26 #include "sphor_accum.h" 27 #include "sphor_compute_mvrea.h" 28 #include "sphor_interface.h" 29 #include "sphor_ran_bsdf.h" 30 #include "sphor_ran_source.h" 31 32 #include <rsys/double3.h> 33 #include <rsys/mem_allocator.h> 34 #include <rsys/rsys.h> 35 #include <star/sphin.h> 36 #include <star/ssp.h> 37 #include <star/ssf.h> 38 39 #include <omp.h> 40 #include <limits.h> 41 42 struct sphin_prop_rad; 43 struct sphin_volume; 44 struct ssf_bsdf; 45 struct ssp_rng; 46 47 /* Syntactic sugar */ 48 #define ALL_COMPONENTS INVALID_ID 49 #define ALL_SENSORS INVALID_ID 50 51 /******************************************************************************* 52 * Helper functions 53 ******************************************************************************/ 54 static res_T 55 setup_MVREA_surfaces_accum 56 (struct sphor* sphor, 57 struct htable_accum2id* accum2id, 58 struct darray_accum* accums) 59 { 60 size_t accum_id = 0; 61 res_T res = RES_OK; 62 63 ASSERT(NULL != sphor); 64 ASSERT(NULL != accum2id); 65 ASSERT(NULL != accums); 66 67 res = register_accum 68 (sphor->allocator, "LOSSES", "\0", "\0", accums); 69 if (RES_OK != res) { goto error; } 70 accum_id = darray_accum_size_get(accums) - 1; 71 72 /* Set corresponding entry in the accum2id hash table */ 73 struct accum_key accum_key = ACCUM_KEY_NULL; 74 accum_key.sensor_type = SPHOR_SENSOR_SURFACE; 75 accum_key.sensor_id = INVALID_ID; 76 accum_key.component_id = INVALID_ID; 77 res = htable_accum2id_set(accum2id, &accum_key, &accum_id); 78 if (RES_OK != res) { goto error; } 79 80 exit: 81 return res; 82 error: 83 goto exit; 84 } 85 86 static res_T 87 setup_MVREA_per_surface_accums 88 (struct sphor* sphor, 89 struct htable_accum2id* accum2id, 90 struct darray_accum* accums) 91 { 92 char* surface_name = NULL; 93 size_t accum_id = 0; 94 size_t surface_count = 0; 95 size_t i_surface = 0; 96 struct sphin_config* config = sphor->config; 97 struct sphin_surface* surface = NULL; 98 struct sphin_sensor_surface* sensor = NULL; 99 res_T res = RES_OK; 100 101 ASSERT(NULL != sphor); 102 ASSERT(NULL != accum2id); 103 ASSERT(NULL != accums); 104 105 res = sphin_config_get_surface_count(config, &surface_count); 106 if (RES_OK != res) { goto error; } 107 108 FOR_EACH(i_surface, 0, surface_count){ 109 res = sphin_config_get_surface(config, i_surface, &surface); 110 if (RES_OK != res) { goto error; } 111 112 res = sphin_surface_get_sensor(surface, &sensor); 113 if (RES_OK != res) { goto error; } 114 if (NULL != sensor) { 115 res = sphin_surface_get_name(surface, &surface_name); 116 if (RES_OK != res) { goto error; } 117 118 res = register_accum 119 (sphor->allocator, "LOSSES", surface_name, "\0", accums); 120 if (RES_OK != res) { goto error; } 121 accum_id = darray_accum_size_get(accums) - 1; 122 123 /* Set corresponding entry in the accum2id hash table */ 124 struct accum_key accum_key = ACCUM_KEY_NULL; 125 accum_key.sensor_type = SPHOR_SENSOR_SURFACE; 126 accum_key.sensor_id = i_surface; 127 res = htable_accum2id_set(accum2id, &accum_key, &accum_id); 128 if (RES_OK != res) { goto error; } 129 } 130 } 131 132 exit: 133 return res; 134 error: 135 goto exit; 136 } 137 138 static res_T 139 setup_MVREA_volumes_accum 140 (struct sphor* sphor, 141 struct htable_accum2id* accum2id, 142 struct darray_accum* accums) 143 { 144 size_t accum_id = 0; 145 res_T res = RES_OK; 146 147 ASSERT(NULL != sphor); 148 ASSERT(NULL != accum2id); 149 ASSERT(NULL != accums); 150 151 res = register_accum 152 (sphor->allocator, "MVREA", "\0", "\0", accums); 153 if (RES_OK != res) { goto error; } 154 accum_id = darray_accum_size_get(accums) - 1; 155 156 /* Set corresponding entry in the accum2id hash table */ 157 struct accum_key accum_key = ACCUM_KEY_NULL; 158 accum_key.sensor_type = SPHOR_SENSOR_VOLUME; 159 accum_key.sensor_id = INVALID_ID; 160 accum_key.component_id = INVALID_ID; 161 res = htable_accum2id_set(accum2id, &accum_key, &accum_id); 162 163 exit: 164 return res; 165 error: 166 goto exit; 167 } 168 169 static res_T 170 setup_MVREA_per_volume_accums 171 (struct sphor* sphor, 172 struct htable_accum2id* accum2id, 173 struct darray_accum* accums) 174 { 175 char* volume_name = NULL; 176 size_t accum_id = 0; 177 size_t volume_count = 0; 178 size_t i_volume = 0; 179 struct sphin_config* config = sphor->config; 180 struct sphin_volume* volume = NULL; 181 struct sphin_sensor_volume* sensor = NULL; 182 res_T res = RES_OK; 183 184 ASSERT(NULL != sphor); 185 ASSERT(NULL != accum2id); 186 ASSERT(NULL != accums); 187 188 res = sphin_config_get_volume_count(config, &volume_count); 189 if (RES_OK != res) { goto error; } 190 191 FOR_EACH(i_volume, 0, volume_count){ 192 res = sphin_config_get_volume(config, i_volume, &volume); 193 if (RES_OK != res) { goto error; } 194 195 res = sphin_volume_get_sensor(volume, &sensor); 196 if (RES_OK != res) { goto error; } 197 if (NULL != sensor) { 198 res = sphin_volume_get_name(volume, &volume_name); 199 if (RES_OK != res) { goto error; } 200 201 res = register_accum 202 (sphor->allocator, "MVREA", volume_name, "\0", accums); 203 if (RES_OK != res) { goto error; } 204 accum_id = darray_accum_size_get(accums) - 1; 205 206 /* Set corresponding entry in the accum2id hash table */ 207 struct accum_key accum_key = ACCUM_KEY_NULL; 208 accum_key.sensor_type = SPHOR_SENSOR_VOLUME; 209 accum_key.sensor_id = i_volume; 210 accum_key.component_id = INVALID_ID; 211 res = htable_accum2id_set(accum2id, &accum_key, &accum_id); 212 if (RES_OK != res) { goto error; } 213 214 } 215 } 216 217 exit: 218 return res; 219 error: 220 goto exit; 221 } 222 223 static res_T 224 setup_MVREA_per_volume_component_accums 225 (struct sphor* sphor, 226 struct htable_accum2id *accum2id, 227 struct darray_accum* accums) 228 { 229 char* volume_name = NULL; 230 char* prop_rad_name = NULL; 231 size_t accum_id = 0; 232 size_t volume_count = 0; 233 size_t prop_rad_count = 0; 234 size_t i_volume = 0, j_proprad = 0; 235 struct sphin_config* config = sphor->config; 236 struct sphin_volume* volume = NULL; 237 struct sphin_sensor_volume* sensor = NULL; 238 struct sphin_prop_rad* prop_rad = NULL; 239 res_T res = RES_OK; 240 241 ASSERT(NULL != sphor); 242 ASSERT(NULL != accum2id); 243 ASSERT(NULL != accums); 244 245 res = sphin_config_get_volume_count(config, &volume_count); 246 if (RES_OK != res) { goto error; } 247 248 FOR_EACH(i_volume, 0, volume_count){ 249 res = sphin_config_get_volume(config, i_volume, &volume); 250 if (RES_OK != res) { goto error; } 251 252 res = sphin_volume_get_sensor(volume, &sensor); 253 if (RES_OK != res) { goto error; } 254 if (NULL != sensor) { 255 res = sphin_volume_get_name(volume, &volume_name); 256 if (RES_OK != res) { goto error; } 257 258 res = sphin_volume_get_prop_rad_count(volume, &prop_rad_count); 259 if (RES_OK != res) { goto error; } 260 261 FOR_EACH(j_proprad, 0, prop_rad_count){ 262 263 res = sphin_volume_get_prop_rad(volume, j_proprad, &prop_rad); 264 if (RES_OK != res) { goto error; } 265 266 res = sphin_prop_rad_get_name(prop_rad, &prop_rad_name); 267 if (RES_OK != res) { goto error; } 268 269 res = register_accum 270 (sphor->allocator, "MVREA", volume_name, prop_rad_name, accums); 271 if (RES_OK != res) { goto error; } 272 accum_id = darray_accum_size_get(accums) - 1; 273 274 /* Set corresponding entry in the accum2id hash table */ 275 struct accum_key accum_key = ACCUM_KEY_NULL; 276 accum_key.sensor_type = SPHOR_SENSOR_VOLUME; 277 accum_key.sensor_id = i_volume; 278 accum_key.component_id = j_proprad; 279 res = htable_accum2id_set(accum2id, &accum_key, &accum_id); 280 if (RES_OK != res) { goto error; } 281 282 } 283 } 284 } 285 286 exit: 287 return res; 288 error: 289 goto exit; 290 } 291 292 static res_T 293 setup_MVREA_accums 294 (struct sphor *sphor, 295 struct htable_accum2id *accum2id, 296 struct darray_accum *accums) 297 { 298 res_T res = RES_OK; 299 300 ASSERT(NULL != sphor); 301 ASSERT(NULL != accums); 302 303 /* Total MVREA in the scene */ 304 res = setup_MVREA_volumes_accum(sphor, accum2id, accums); 305 if (RES_OK != res) { goto error; } 306 307 /* One accum per sensor volume; all components of the volume contribute 308 * to the accum */ 309 res = setup_MVREA_per_volume_accums(sphor, accum2id, accums); 310 if (RES_OK != res) { goto error; } 311 312 /* One accum per component per sensor volume */ 313 res = setup_MVREA_per_volume_component_accums(sphor, accum2id, accums); 314 if (RES_OK != res) { goto error; } 315 316 /* Total losses in the scene */ 317 res = setup_MVREA_surfaces_accum(sphor, accum2id, accums); 318 if (RES_OK != res) { goto error; } 319 320 /* One accum per sensor surface */ 321 res = setup_MVREA_per_surface_accums(sphor, accum2id, accums); 322 if (RES_OK != res) { goto error; } 323 exit: 324 return res; 325 error: 326 goto exit; 327 } 328 329 static res_T 330 prop_rad_compute_ka 331 (struct sphin_prop_rad* prop_rad, 332 double wavelength, 333 double* ka) 334 { 335 struct sphin_scatterer* scatterer = NULL; 336 struct sphin_spectral_property* abs_cross_sec = NULL; 337 double concentration = 0; 338 double sigma_a = 0; 339 res_T res = RES_OK; 340 341 res = sphin_prop_rad_get_scatterer(prop_rad, &scatterer); 342 if (RES_OK != res) { goto error; } 343 344 res = sphin_scatterer_get_concentration(scatterer, &concentration); 345 if (RES_OK != res) { goto error; } 346 347 res = sphin_scatterer_get_abs_cross_sec(scatterer, &abs_cross_sec); 348 if (RES_OK != res) { goto error; } 349 350 res = sphin_spectral_property_interpolate_at_wavelength 351 (abs_cross_sec, wavelength, SPHIN_INTERPOLATION_LINEAR, &sigma_a); 352 if (RES_OK != res) { goto error; } 353 354 *ka = concentration * sigma_a; 355 exit: 356 return res; 357 error: 358 goto exit; 359 } 360 361 static res_T 362 volume_compute_total_ka 363 (struct sphor* sphor, 364 const struct sphin_volume* volume, /* May be NULL => out_ka = 0 */ 365 double wavelength, 366 double* out_ka) 367 { 368 size_t prop_rad_count = 0; 369 size_t j = 0; 370 double ka = 0; 371 struct sphin_prop_rad* prop_rad = NULL; 372 res_T res = RES_OK; 373 374 ASSERT(NULL != sphor); 375 ASSERT(NULL != out_ka); 376 377 (void)sphor; 378 379 /* No volume defined in the interface side */ 380 if (NULL == volume) { 381 ka = 0.; 382 goto exit; 383 } 384 385 /* Get ka corresponding to the current wavelength */ 386 res = sphin_volume_get_prop_rad_count 387 ((struct sphin_volume*)volume, &prop_rad_count); 388 if (RES_OK != res) { goto error; } 389 390 FOR_EACH(j, 0, prop_rad_count) { 391 double ka_i = 0; 392 res = sphin_volume_get_prop_rad((struct sphin_volume*)volume, j, &prop_rad); 393 if (RES_OK != res) { goto error; } 394 395 res = prop_rad_compute_ka(prop_rad, wavelength, &ka_i); 396 if (RES_OK != res) { goto error; } 397 398 ka += ka_i; 399 } 400 401 exit: 402 *out_ka = ka; 403 return res; 404 error: 405 goto exit; 406 } 407 408 static struct sphin_volume* 409 volume_get_from_primitive 410 (struct sphor* sphor, 411 const struct primitive* primitive) 412 { 413 struct sphin_volume* volume = NULL; 414 size_t volume_id = INVALID_ID; 415 416 ASSERT(NULL != sphor); 417 ASSERT(NULL != primitive); 418 419 volume_id = primitive->interface->volumes[primitive->side]; 420 421 if (INVALID_ID != volume_id) { 422 SPHIN(config_get_volume(sphor->config, volume_id, &volume)); 423 } 424 425 return volume; 426 } 427 428 static res_T 429 MVREA_update_volume_weights 430 (struct sphor* sphor, 431 struct intersection* intersection, 432 double wavelength, 433 struct htable_accum2id* accum2id, 434 struct darray_accum* accums) 435 { 436 struct accum_key accum_key = ACCUM_KEY_NULL; 437 size_t prop_rad_count = 0; 438 size_t volume_id = SIZE_MAX; 439 size_t i_proprad = 0; 440 double response_function = 0; 441 double vol_size = 0; 442 double tot_pow = 0; 443 double ka = 0; 444 double mc_weight = 0; 445 struct sphin_volume* volume = NULL; 446 struct sphin_sensor_volume* sensor_volume = NULL; 447 res_T res = RES_OK; 448 449 ASSERT(NULL != sphor); 450 ASSERT(NULL != intersection); 451 ASSERT(NULL != accum2id); 452 ASSERT(NULL != accums); 453 454 volume = volume_get_from_primitive(sphor, &intersection->position.primitive); 455 456 /* Get the total three-dimensional space occupied by the volume */ 457 res = sphin_volume_compute_total_size(volume, &vol_size/*[m^3]*/); 458 if (RES_OK != res) { goto error; } 459 res = sphin_config_get_total_surface_power(sphor->config, &tot_pow); 460 if (RES_OK != res) { goto error; } 461 462 mc_weight = tot_pow / vol_size; 463 464 res = sphin_volume_get_sensor(volume, &sensor_volume); 465 if (RES_OK != res) { goto error; } 466 467 res = sphin_sensor_volume_get_response_function 468 (sensor_volume, &response_function); 469 if (RES_OK != res) { goto error; } 470 res = sphin_volume_get_prop_rad_count(volume, &prop_rad_count); 471 if (RES_OK != res) { goto error; } 472 473 volume_id = intersection->position.primitive.interface->volumes 474 [intersection->position.primitive.side]; 475 476 /* Update total accumulator for all volumes in the scene*/ 477 accum_key.sensor_type = SPHOR_SENSOR_VOLUME; 478 accum_key.sensor_id = ALL_SENSORS; 479 accum_key.component_id = ALL_COMPONENTS; 480 update_accum(accum2id, &accum_key, mc_weight, accums); 481 482 /* Update total accumulator in the volume */ 483 accum_key.sensor_type = SPHOR_SENSOR_VOLUME; 484 accum_key.sensor_id = volume_id; 485 accum_key.component_id = ALL_COMPONENTS; 486 update_accum(accum2id, &accum_key, mc_weight, accums); 487 488 /* Increment the weight of each chemical species in the volume by 489 * ( ka_i / ka_tot ) response_function */ 490 res = volume_compute_total_ka(sphor, volume, wavelength, &ka); 491 if (RES_OK != res) { goto error; } 492 FOR_EACH(i_proprad, 0, prop_rad_count) { 493 double ka_i = 0; 494 struct sphin_prop_rad* prop_rad = NULL; 495 res = sphin_volume_get_prop_rad(volume, i_proprad, &prop_rad); 496 if (RES_OK != res) { goto error; } 497 res = prop_rad_compute_ka(prop_rad, wavelength, &ka_i); 498 if (RES_OK != res) { goto error; } 499 accum_key.component_id = i_proprad; 500 501 mc_weight = ka_i / ka * tot_pow / vol_size; 502 update_accum(accum2id, &accum_key, mc_weight, accums); 503 } 504 505 exit: 506 return res; 507 error: 508 goto exit; 509 } 510 511 static res_T 512 MVREA_update_surface_weights 513 (struct sphor* sphor, 514 struct intersection* intersection, 515 double wavelength, 516 struct htable_accum2id* accum2id, 517 struct darray_accum* accums) 518 { 519 struct accum_key accum_key = ACCUM_KEY_NULL; 520 size_t surface_id = SIZE_MAX; 521 size_t i = 0; 522 double response_function = 0; 523 double tot_pow = 0; 524 double mc_weight = 0; 525 double surf_area = 0; 526 struct sphin_surface* surface = NULL; 527 struct sphin_sensor_surface* sensor_surface = NULL; 528 struct primitive prim = PRIMITIVE_NULL; 529 res_T res = RES_OK; 530 531 ASSERT(NULL != sphor); 532 ASSERT(NULL != intersection); 533 ASSERT(NULL != accum2id); 534 ASSERT(NULL != accums); 535 536 (void)wavelength; 537 538 prim = intersection->position.primitive; 539 540 /* Retrieve the sensor and the surface_id corresponding to the sensor */ 541 FOR_EACH(i, 0, (size_t)prim.interface->surface_count[(size_t)prim.side]) { 542 surface_id = (size_t)prim.interface->surfaces[prim.side][i]; 543 res = sphin_config_get_surface(sphor->config, surface_id, &surface); 544 if (RES_OK != res) { goto error; } 545 res = sphin_surface_get_sensor(surface, &sensor_surface); 546 if (RES_OK != res) { goto error; } 547 548 if (NULL != sensor_surface) { break; } 549 } 550 551 if (NULL == sensor_surface) { res = RES_BAD_ARG; goto error; } 552 553 res = sphin_surface_compute_total_area(surface, &surf_area); 554 if (RES_OK != res) { goto error; } 555 res = sphin_config_get_total_surface_power(sphor->config, &tot_pow); 556 if (RES_OK != res) { goto error; } 557 res = sphin_sensor_surface_get_response_function 558 (sensor_surface, &response_function); 559 if (RES_OK != res) { goto error; } 560 561 mc_weight = response_function * tot_pow / surf_area; 562 563 /* Update total accumulator for all surfaces losses in the scene*/ 564 accum_key.sensor_type = SPHOR_SENSOR_SURFACE; 565 accum_key.sensor_id = ALL_SENSORS; 566 update_accum(accum2id, &accum_key, mc_weight, accums); 567 568 /* Update total accumulator in the surface */ 569 accum_key.sensor_type = SPHOR_SENSOR_SURFACE; 570 accum_key.sensor_id = surface_id; 571 update_accum(accum2id, &accum_key, mc_weight, accums); 572 573 exit: 574 return res; 575 error: 576 goto exit; 577 } 578 579 static res_T 580 MVREA_write_outputs 581 (struct sphor* sphor, 582 const struct darray_accum* accums, 583 size_t nfailures) 584 { 585 size_t samples = 0; 586 size_t i = 0; 587 struct accum* accum; 588 res_T res = RES_OK; 589 590 ASSERT(NULL != sphor); 591 ASSERT(NULL != accums); 592 593 samples = sphor->samples - nfailures; 594 595 FOR_EACH(i, 0, darray_accum_size_get(accums)){ 596 accum = (struct accum*)darray_accum_cdata_get(accums) + i; 597 accum->n_realizations = samples; 598 write_accum_estim(accum, sphor->stream); 599 } 600 601 return res; 602 } 603 604 static res_T 605 compute_MVREA_realization 606 (struct sphor* sphor, 607 struct ssp_rng* rng, 608 struct htable_accum2id* accum2id, 609 struct darray_accum* accums) 610 { 611 /* Star-Phor Input */ 612 struct sphin_sensor_volume* sensor_volume = NULL; 613 struct sphin_sensor_surface* sensor_surface = NULL; 614 struct sphin_volume* volume = NULL; 615 616 /* The sampled source */ 617 struct source_view* source_view = NULL; 618 619 /* Ray */ 620 enum intersection_ray_interaction_type intersection_ray_interaction_type = 621 INTERSECTION_RAY_INTERACTION_NONE__; 622 struct intersection intersection = INTERSECTION_NULL; 623 struct ray ray = RAY_DEFAULT; 624 625 /* Miscellaneous */ 626 double dir[3] = {0}; /* Sampled direction during path sampling */ 627 double ka = 0; 628 double free_path = 0; 629 double wavelength = 0; 630 int stop = 0; 631 632 /* For error management */ 633 res_T res = RES_OK; 634 #define CALL(Function) { \ 635 res = Function; \ 636 if (RES_OK != res) goto error; \ 637 } (void)0 638 639 /* Pre-conditions */ 640 ASSERT(NULL != sphor); 641 ASSERT(NULL != rng); 642 ASSERT(NULL != accum2id); 643 ASSERT(NULL != accums); 644 645 /* Sample a random position from a source in the whole scene */ 646 CALL(sample_source_position(sphor, 647 /* in: */ rng, 648 /* out: */ &source_view, &ray.origin_on_prim, ray.origin)); 649 650 /* Sample a direction according to the source direction distribution */ 651 CALL(source_surface_sample_direction(sphor, 652 /* in: */ rng, &ray.origin_on_prim, 653 /* out: */ ray.direction)); 654 655 /* Sample a wavelength according to the source emission spectrum */ 656 CALL(source_sample_wavelength 657 (sphor, rng, source_view, &wavelength)); 658 659 /* Retrieve the volume defined by the sampled primitive */ 660 volume = volume_get_from_primitive(sphor, &ray.origin_on_prim.primitive); 661 662 /* Retrieve the ka of the present volume */ 663 CALL(volume_compute_total_ka(sphor, volume, wavelength, &ka)); 664 665 /* While no absorption takes place and an interface is found */ 666 while(!stop) { 667 /* Sample the free path length to absorption from an exponential 668 * distribution with rate parameter ka */ 669 free_path = ssp_ran_exp(rng, ka); 670 671 /* Trace a ray in the scene from the sampled primitive in the sampled 672 * direction and retrive the hit distance */ 673 CALL(trace_ray(sphor, &ray, &intersection)); 674 675 /* If there is no intersection in the sampled direction from the sampled 676 * position, the photon is lost and it counts zero in the weight. */ 677 if (INTERSECTION_NONE(&intersection)) { break; } 678 679 /* If the distance between the ray origin and the intersection is bigger 680 * than the absorption free path length, an absorption takes place (check 681 * if the volume is a sensor and return the weight accordingly) */ 682 if (free_path < intersection.distance) { 683 CALL(sphin_volume_get_sensor(volume, &sensor_volume)); 684 685 /* Check if the volume is a sensor */ 686 if (NULL == sensor_volume) { 687 /* If the volume is not a sensor, the absorbed photon does not count in 688 * the weight. The weight is implictly 0 */ 689 } else { 690 /* If the volume is a sensor, compute the weight using the response 691 * function */ 692 CALL(MVREA_update_volume_weights(sphor, 693 /* in */ &intersection, wavelength, accum2id, 694 /* out */ accums)); 695 } 696 697 /* Absorption by a volume => stop the path */ 698 stop = 1; 699 } else { 700 /* Else, the distance between the ray origin and the intersection is 701 * smaller than the absorption free path length. The photon intersects a 702 * primitive in the scene view. Sample the type of interaction that the 703 * photon has with the interface as well as the new direction sampled */ 704 CALL(sample_interface_ray_interaction(sphor, 705 /* in: */ rng, &ray, &intersection, wavelength, 706 /* out: */ dir, &intersection_ray_interaction_type)); 707 708 if (INTERSECTION_RAY_INTERACTION_ABSORPTION == 709 intersection_ray_interaction_type) { 710 CALL(primitive_get_sensor_surface 711 (sphor, &intersection.position.primitive, &sensor_surface)); 712 713 /* If the volume is a sensor, compute the weight using the response 714 * function */ 715 if (NULL != sensor_surface) { 716 CALL(MVREA_update_surface_weights(sphor, 717 /* in */ &intersection, wavelength, accum2id, 718 /* out */ accums)); 719 } 720 721 /* Absorption by a surface => stop the path */ 722 stop = 1; 723 } 724 else if (INTERSECTION_RAY_INTERACTION_TRANSMISSION == 725 intersection_ray_interaction_type) { 726 /* Continue the path on the other side of the primitive.*/ 727 intersection.position.primitive.side = 728 !intersection.position.primitive.side; 729 730 CALL(ray_update(/* in/out: */ &ray, /* in: */ &intersection, dir)); 731 732 /* Update the current medium */ 733 volume = volume_get_from_primitive 734 (sphor, &intersection.position.primitive); 735 736 CALL(volume_compute_total_ka(sphor, volume, wavelength, &ka)); 737 } 738 else if (INTERSECTION_RAY_INTERACTION_REFLECTION == 739 intersection_ray_interaction_type) { 740 CALL(ray_update(/* in/out: */ &ray, /* in: */ &intersection, dir)); 741 } 742 } 743 } 744 745 #undef CALL 746 747 exit: 748 return res; 749 error: 750 goto exit; 751 } 752 753 /******************************************************************************* 754 * Local functions 755 ******************************************************************************/ 756 res_T 757 sphor_compute_MVREA 758 (struct sphor* sphor) 759 { 760 res_T res = RES_OK; 761 size_t nthreads = 0; 762 size_t samples = 0; 763 size_t nfailures = 0; 764 size_t i = 0; /* iterator */ 765 struct ssp_rng_proxy *rng_proxy = NULL; 766 struct ssp_rng **rngs = NULL; 767 struct darray_accum accums; 768 struct htable_accum2id accum2id; 769 struct darray_accum_list accums_threads; 770 771 samples = sphor->samples; 772 nthreads = sphor->nthreads; 773 if (UINT_MAX == nthreads) { /* use all threads available if not in args */ 774 nthreads = (size_t)omp_get_num_procs(); 775 } 776 777 /* Initialize the global accumulator and the hash table that allows to 778 * retrieve the id of a given accumulator with an accum_key */ 779 darray_accum_init(sphor->allocator, &accums); 780 htable_accum2id_init(sphor->allocator, &accum2id); 781 782 /* Create and initialize dynamic array o accums (one per observable) */ 783 res = setup_MVREA_accums(sphor, &accum2id, &accums); 784 if (RES_OK != res) { goto error; } 785 786 /* Create an array containing nthreads pointers to different accums lists */ 787 darray_accum_list_init(sphor->allocator, &accums_threads); 788 res = darray_accum_list_resize(&accums_threads, nthreads); 789 if (RES_OK != res) { goto error; } 790 791 /* Mem allocation of the accummulator list of each one of the threads 792 * Each thread handles a copy of the global accum list just created */ 793 FOR_EACH(i, 0, nthreads) { 794 struct darray_accum* accum_thread; 795 accum_thread = darray_accum_list_data_get(&accums_threads) + i; 796 res = darray_accum_copy(accum_thread, &accums); 797 if (RES_OK != res) { goto error; } 798 } 799 800 /* Create of the proxy generator RNG_MT19937_64 (Mersenne Twister) */ 801 res = ssp_rng_proxy_create(NULL, SSP_RNG_MT19937_64, nthreads, &rng_proxy); 802 if (res != RES_OK) { goto error; } 803 804 /* Create an array containing nthreads pointers to different rngs */ 805 rngs = mem_calloc(nthreads, sizeof(*rngs)); 806 if (NULL == rngs) { res = RES_MEM_ERR; goto error; } 807 808 /* Set one generator per thread using the proxy generator */ 809 FOR_EACH(i, 0, nthreads) { 810 res = ssp_rng_proxy_create_rng(rng_proxy, i, &rngs[i]); 811 if (RES_OK != res) { goto error; } 812 } 813 omp_set_num_threads((int)nthreads); 814 815 /* Realizations loop */ 816 nfailures = 0; 817 #pragma omp parallel for schedule(static) 818 for(i=0; i<samples; i++) { 819 const int ithread = omp_get_thread_num(); 820 res_T res_local = RES_OK; 821 struct darray_accum* accum_thread; 822 accum_thread = darray_accum_list_data_get(&accums_threads) + ithread; 823 824 /* Ignore the rest of the for loop if there is a fatal (not failure) error 825 * in the previous realizations */ 826 if (RES_OK != res) continue; 827 828 res_local = compute_MVREA_realization 829 (sphor, rngs[ithread], &accum2id, accum_thread); 830 if (RES_OK != res_local) { 831 /*Protect res and nfailures from concurrent write accesses*/ 832 #pragma omp critical 833 switch(res_local) { 834 /* Fatal, all remaining realizations will be skipped */ 835 case RES_UNKNOWN_ERR: 836 res = res_local; 837 break; 838 default: 839 nfailures += 1; 840 break; 841 } 842 } 843 } 844 if (res != RES_OK) { goto error; } 845 846 /* Aggregate accums from the different threads in the main accumulator */ 847 sum_accum_lists(&accums_threads, &accums); 848 849 if (0 < nfailures) { 850 WARN(sphor, 851 "%lu effective realizations and %lu failures out of %lu samples.\n", 852 samples-nfailures, nfailures, samples); 853 } 854 855 res = MVREA_write_outputs(sphor, &accums, nfailures); 856 if (res != RES_OK) { goto error; } 857 858 exit: 859 /* Free the proxy generator and each one of the thread generators */ 860 if (NULL != rng_proxy) ssp_rng_proxy_ref_put(rng_proxy); 861 if (NULL != rngs) { 862 FOR_EACH(i, 0, nthreads) { 863 if (NULL != rngs[i]) ssp_rng_ref_put(rngs[i]); 864 } 865 mem_rm(rngs); 866 } 867 /* Free local data structures accums and accum2id */ 868 darray_accum_release(&accums); 869 darray_accum_list_release(&accums_threads); 870 htable_accum2id_release(&accum2id); 871 872 return res; 873 error: 874 goto exit; 875 }