ssol_solver.c (48210B)
1 /* Copyright (C) 2018-2026 |Meso|Star> (contact@meso-star.com) 2 * Copyright (C) 2016, 2018 CNRS 3 * 4 * This program is free software: you can redistribute it and/or modify 5 * it under the terms of the GNU General Public License as published by 6 * the Free Software Foundation, either version 3 of the License, or 7 * (at your option) any later version. 8 * 9 * This program is distributed in the hope that it will be useful, 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 * GNU General Public License for more details. 13 * 14 * You should have received a copy of the GNU General Public License 15 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 16 17 #define _POSIX_C_SOURCE 200112L /* nextafterf support */ 18 19 #include "ssol.h" 20 #include "ssol_c.h" 21 #include "ssol_atmosphere_c.h" 22 #include "ssol_device_c.h" 23 #include "ssol_estimator_c.h" 24 #include "ssol_scene_c.h" 25 #include "ssol_shape_c.h" 26 #include "ssol_object_c.h" 27 #include "ssol_sun_c.h" 28 #include "ssol_material_c.h" 29 #include "ssol_instance_c.h" 30 #include "ssol_ranst_sun_dir.h" 31 #include "ssol_ranst_sun_wl.h" 32 33 #include <rsys/float2.h> 34 #include <rsys/float3.h> 35 #include <rsys/double3.h> 36 #include <rsys/mem_allocator.h> 37 #include <rsys/ref_count.h> 38 #include <rsys/rsys.h> 39 #include <rsys/stretchy_array.h> 40 41 #include <star/ssf.h> 42 #include <star/ssp.h> 43 44 #include <limits.h> 45 #include <omp.h> 46 47 /******************************************************************************* 48 * Thread context 49 ******************************************************************************/ 50 struct thread_context { 51 struct ssp_rng* rng; 52 struct mc_data cos_factor; 53 struct mc_data absorbed_by_receivers; 54 struct mc_data shadowed; 55 struct mc_data missing; 56 struct mc_data extinguished_by_atmosphere; 57 struct mc_data other_absorbed; 58 struct htable_receiver mc_rcvs; 59 struct htable_sampled mc_samps; 60 struct darray_path paths; /* paths */ 61 size_t realisation_count; 62 }; 63 64 static void 65 thread_context_release(struct thread_context* ctx) 66 { 67 ASSERT(ctx); 68 if(ctx->rng) SSP(rng_ref_put(ctx->rng)); 69 htable_receiver_release(&ctx->mc_rcvs); 70 htable_sampled_release(&ctx->mc_samps); 71 darray_path_release(&ctx->paths); 72 } 73 74 static res_T 75 thread_context_init(struct mem_allocator* allocator, struct thread_context* ctx) 76 { 77 ASSERT(ctx); 78 memset(ctx, 0, sizeof(ctx[0])); 79 htable_receiver_init(allocator, &ctx->mc_rcvs); 80 htable_sampled_init(allocator, &ctx->mc_samps); 81 darray_path_init(allocator, &ctx->paths); 82 return RES_OK; 83 } 84 85 /* Define a copy functor only for consistency since this function will not be 86 * used */ 87 static res_T 88 thread_context_copy 89 (struct thread_context* dst, const struct thread_context* src) 90 { 91 res_T res = RES_OK; 92 ASSERT(dst && src); 93 dst->rng = src->rng; 94 dst->cos_factor = src->cos_factor; 95 dst->absorbed_by_receivers = src->absorbed_by_receivers; 96 dst->shadowed = src->shadowed; 97 dst->missing = src->missing; 98 dst->extinguished_by_atmosphere = src->extinguished_by_atmosphere; 99 dst->other_absorbed = src->other_absorbed; 100 res = htable_receiver_copy(&dst->mc_rcvs, &src->mc_rcvs); 101 if(res != RES_OK) return res; 102 res = htable_sampled_copy(&dst->mc_samps, &src->mc_samps); 103 if(res != RES_OK) return res; 104 res = darray_path_copy(&dst->paths, &src->paths); 105 if(res != RES_OK) return res; 106 return RES_OK; 107 } 108 109 static void 110 thread_context_clear(struct thread_context* ctx) 111 { 112 ASSERT(ctx); 113 if(ctx->rng) SSP(rng_ref_put(ctx->rng)); 114 htable_receiver_clear(&ctx->mc_rcvs); 115 htable_sampled_clear(&ctx->mc_samps); 116 darray_path_clear(&ctx->paths); 117 } 118 119 static res_T 120 thread_context_setup 121 (struct thread_context* ctx, 122 struct ssp_rng_proxy* rng_proxy, 123 const size_t ithread) 124 { 125 res_T res = RES_OK; 126 ASSERT(rng_proxy && ctx); 127 thread_context_clear(ctx); 128 res = ssp_rng_proxy_create_rng(rng_proxy, ithread, &ctx->rng); 129 if(res != RES_OK) goto error; 130 exit: 131 return res; 132 error: 133 thread_context_clear(ctx); 134 goto exit; 135 } 136 137 /* Declare the container of the per thread contexts */ 138 #define DARRAY_NAME thread_ctx 139 #define DARRAY_DATA struct thread_context 140 #define DARRAY_FUNCTOR_INIT thread_context_init 141 #define DARRAY_FUNCTOR_RELEASE thread_context_release 142 #define DARRAY_FUNCTOR_COPY thread_context_copy 143 #include <rsys/dynamic_array.h> 144 145 /******************************************************************************* 146 * Random walk point 147 ******************************************************************************/ 148 struct point { 149 const struct ssol_instance* inst; 150 const struct shaded_shape* sshape; 151 struct mc_sampled* mc_samp; 152 struct s3d_primitive prim; 153 double N[3]; 154 double pos[3]; 155 double dir[3]; 156 float uv[2]; 157 double wl; /* Sampled wavelength */ 158 const struct ssol_material* material; 159 /* tmp quantities to compute weights */ 160 double kabs_at_pt; 161 size_t survivor_score; 162 /* for conservation of energy check */ 163 double energy_loss; 164 /* MC weights */ 165 /* Set once */ 166 double initial_flux; /* the initial flux*/ 167 double cos_factor; /* local cos at the starting point */ 168 /* outgoing weights at previous hit */ 169 double prev_outgoing_flux; 170 double prev_outgoing_if_no_atm_loss; 171 double prev_outgoing_if_no_field_loss; 172 /* incoming weights at current hit */ 173 double incoming_flux; 174 double incoming_if_no_atm_loss; 175 double incoming_if_no_field_loss; 176 /* outgoing weights at current hit */ 177 double outgoing_flux; 178 double outgoing_if_no_atm_loss; 179 double outgoing_if_no_field_loss; 180 enum ssol_side_flag side; 181 }; 182 183 #define POINT_NULL__ { \ 184 NULL, /* Instance */ \ 185 NULL, /* Shaded shape */ \ 186 NULL, /* Primary data */ \ 187 S3D_PRIMITIVE_NULL__, /* Primitive */ \ 188 {0, 0, 0}, /* Normal */ \ 189 {0, 0, 0}, /* Position */ \ 190 {0, 0, 0}, /* Direction */ \ 191 {0, 0}, /* UV */ \ 192 0, /* Wavelength */ \ 193 NULL, /* Material */ \ 194 0, 0, /* tmp values */ \ 195 0, /* Energy loss */ \ 196 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, /* MC weights */ \ 197 SSOL_FRONT /* Side */ \ 198 } 199 static const struct point POINT_NULL = POINT_NULL__; 200 201 static FINLINE struct ssol_material* 202 point_get_material(const struct point* pt) 203 { 204 return pt->side == SSOL_FRONT ? pt->sshape->mtl_front : pt->sshape->mtl_back; 205 } 206 207 static res_T 208 point_init 209 (struct point* pt, 210 struct ssol_scene* scn, 211 struct htable_sampled* sampled, 212 struct s3d_scene_view* view_samp, 213 struct s3d_scene_view* view_rt, 214 struct ranst_sun_dir* ran_sun_dir, 215 struct ranst_sun_wl* ran_sun_wl, 216 struct ssp_rng* rng, 217 struct ssol_medium* current_medium, 218 int* is_lit) 219 { 220 struct s3d_attrib attr; 221 struct s3d_hit hit; 222 struct ray_data ray_data = RAY_DATA_NULL; 223 double N[3]; 224 double surface_sun_cos; 225 double surface_sun0_cos; 226 double sun0_sun_cos; 227 double surface_proxy_cos; 228 double cos_ratio; 229 double w0; 230 float dir[3], pos[3], range[2] = { 0, FLT_MAX }; 231 size_t id; 232 res_T res = RES_OK; 233 ASSERT(pt && scn && sampled && view_samp && view_rt); 234 ASSERT(ran_sun_dir && ran_sun_wl && rng && is_lit); 235 236 /* Sample a point into the scene view */ 237 S3D(scene_view_sample 238 (view_samp, 239 ssp_rng_canonical_float(rng), 240 ssp_rng_canonical_float(rng), 241 ssp_rng_canonical_float(rng), 242 &pt->prim, pt->uv)); 243 244 /* Retrieve the position of the sampled point */ 245 S3D(primitive_get_attrib(&pt->prim, S3D_POSITION, pt->uv, &attr)); 246 d3_set_f3(pt->pos, attr.value); 247 248 /* Retrieve the normal of the sampled point */ 249 S3D(primitive_get_attrib(&pt->prim, S3D_GEOMETRY_NORMAL, pt->uv, &attr)); 250 f3_normalize(attr.value, attr.value); 251 d3_set_f3(pt->N, attr.value); 252 253 /* Retrieve the sampled instance and shaded shape */ 254 pt->inst = *htable_instance_find(&scn->instances_samp, &pt->prim.inst_id); 255 id = *htable_shaded_shape_find 256 (&pt->inst->object->shaded_shapes_samp, &pt->prim.geom_id); 257 pt->sshape = darray_shaded_shape_cdata_get 258 (&pt->inst->object->shaded_shapes) + id; 259 260 /* Sample a sun direction */ 261 ranst_sun_dir_get(ran_sun_dir, rng, pt->dir); 262 263 /* Sample a wavelength */ 264 pt->wl = ranst_sun_wl_get(ran_sun_wl, rng); 265 266 if(pt->sshape->shape->type != SHAPE_PUNCHED) { 267 d3_set(N, pt->N); 268 } else { 269 /* For punched surface, retrieve the sampled position and normal onto the 270 * quadric surface */ 271 punched_shape_project_point 272 (pt->sshape->shape, pt->inst->transform, pt->pos, pt->pos, N); 273 } 274 275 /* Define the primitive side on which the point lies */ 276 if(d3_dot(N, pt->dir) < 0) { 277 pt->side = SSOL_FRONT; 278 } else { 279 pt->side = SSOL_BACK; 280 d3_minus(N, N); /* Force the normal to look forward dir */ 281 } 282 283 /* Initialise the Monte Carlo weight */ 284 surface_sun_cos = d3_dot(N, pt->dir); 285 surface_sun0_cos = fabs(d3_dot(scn->sun->direction, N)); 286 sun0_sun_cos = d3_dot(scn->sun->direction, pt->dir); 287 surface_proxy_cos = 288 (pt->sshape->shape->type == SHAPE_MESH) ? 1 : fabs(d3_dot(pt->N, N)); 289 cos_ratio = fabs(surface_sun_cos / (surface_proxy_cos * sun0_sun_cos)); 290 w0 = scn->sun->dni * scn->sampled_area_proxy * cos_ratio; 291 pt->cos_factor = scn->sampled_area_proxy / scn->sampled_area 292 * surface_sun0_cos / surface_proxy_cos; 293 pt->energy_loss = w0; 294 pt->initial_flux = w0; 295 pt->prev_outgoing_flux = w0; 296 pt->prev_outgoing_if_no_atm_loss = w0; 297 pt->prev_outgoing_if_no_field_loss = w0; 298 pt->survivor_score = 0; 299 d3_set(pt->N, N); 300 ASSERT(d3_dot(pt->N, pt->dir) <= 0); 301 302 /* Store sampled entity related weights */ 303 res = get_mc_sampled(sampled, pt->inst, &pt->mc_samp); 304 if(res != RES_OK) goto error; 305 pt->mc_samp->nb_samples++; 306 307 /* Define the medium in which the sampled point lies */ 308 pt->material = point_get_material(pt); 309 switch (pt->material->type) { 310 case SSOL_MATERIAL_DIELECTRIC: 311 case SSOL_MATERIAL_THIN_DIELECTRIC: 312 /* TODO: check sampled face role!!! */ 313 ssol_medium_copy(current_medium, 314 (pt->side == SSOL_FRONT) ? 315 &pt->material->in_medium : &pt->material->out_medium); 316 break; 317 case SSOL_MATERIAL_MATTE: 318 case SSOL_MATERIAL_MIRROR: 319 case SSOL_MATERIAL_VIRTUAL: 320 ssol_medium_copy(current_medium, &scn->air); 321 break; 322 default: FATAL("Unreachable code\n"); break; 323 } 324 325 /* Initialise the ray data to avoid self intersection */ 326 ray_data.scn = scn; 327 ray_data.prim_from = pt->prim; 328 ray_data.inst_from = pt->inst; 329 ray_data.sshape_from = pt->sshape; 330 ray_data.side_from = pt->side; 331 ray_data.discard_virtual_materials = 1; /* Do not intersect virtual mtl */ 332 ray_data.reversed_ray = 1; /* The ray direction is reversed */ 333 ray_data.dst = FLT_MAX; 334 335 /* pt->prim must live in RT space */ 336 f3_set_d3(pos, pt->pos); 337 ray_data.point_init_closest_point = 1; 338 S3D(shape_get_id(pt->sshape->shape->shape_rt, &ray_data.prim_from.geom_id)); 339 S3D(shape_get_id(pt->inst->shape_rt, &ray_data.prim_from.inst_id)); 340 S3D(scene_view_closest_point(view_rt, pos, FLT_MAX, &ray_data, &hit)); 341 CHK(!S3D_HIT_NONE(&hit)); 342 /* Sample and RT meshes are supposed to be identical only for SHAPE_MESH */ 343 ASSERT(pt->sshape->shape->type != SHAPE_MESH 344 || hit.distance <= (1 + f3_len(pos)) * 1e-6); 345 pt->prim = hit.prim; 346 ray_data.prim_from = pt->prim; 347 348 /* Trace a ray toward the sun to check if the sampled point is occluded */ 349 f3_minus(dir, f3_set_d3(dir, pt->dir)); 350 ray_data.point_init_closest_point = 0; 351 S3D(scene_view_trace_ray(view_rt, pos, dir, range, &ray_data, &hit)); 352 *is_lit = S3D_HIT_NONE(&hit); 353 354 exit: 355 return res; 356 error: 357 goto exit; 358 } 359 360 static FINLINE void 361 point_update_from_hit 362 (struct point* pt, 363 struct ssol_scene* scn, /* Scene into which the hit occurs */ 364 const float org[3], /* Origin of the ray that generates the hit */ 365 const float dir[3], /* Direction of the ray that generates the hit */ 366 const struct s3d_hit* hit, 367 struct ray_data* rdata) /* Ray data used to generate the hit */ 368 { 369 double tmp[3]; 370 float tmpf[3]; 371 size_t id; 372 373 /* Retrieve the hit instance and shaded shape */ 374 pt->inst = *htable_instance_find(&scn->instances_rt, &hit->prim.inst_id); 375 id = *htable_shaded_shape_find 376 (&pt->inst->object->shaded_shapes_rt, &hit->prim.geom_id); 377 pt->sshape = darray_shaded_shape_cdata_get 378 (&pt->inst->object->shaded_shapes) + id; 379 380 /* Fetch the current position and its associated normal */ 381 switch(pt->sshape->shape->type) { 382 case SHAPE_MESH: 383 d3_set_f3(pt->N, hit->normal); 384 d3_normalize(pt->N, pt->N); 385 f3_mulf(tmpf, dir, hit->distance); 386 f3_add(tmpf, org, tmpf); 387 d3_set_f3(pt->pos, tmpf); 388 break; 389 case SHAPE_PUNCHED: 390 d3_normalize(pt->N, rdata->N); 391 d3_muld(tmp, pt->dir, rdata->dst); 392 f3_set_d3(tmpf, tmp); 393 f3_add(tmpf, org, tmpf); 394 d3_set_f3(pt->pos, tmpf); 395 break; 396 default: FATAL("Unreachable code"); break; 397 } 398 399 pt->prim = hit->prim; 400 401 /* Define the primitive side on which the point lies */ 402 if(d3_dot(pt->dir, pt->N) < 0) { 403 pt->side = SSOL_FRONT; 404 } else { 405 pt->side = SSOL_BACK; 406 d3_minus(pt->N, pt->N); /* Force the normal to look forward dir */ 407 } 408 409 /* Update material */ 410 pt->material = point_get_material(pt); 411 } 412 413 static FINLINE int 414 point_is_receiver(const struct point* pt) 415 { 416 return (pt->inst->receiver_mask & (int)pt->side) != 0; 417 } 418 419 static FINLINE res_T 420 point_shade 421 (struct point* pt, 422 const struct ssol_medium* in_medium, 423 struct ssol_medium* out_medium, 424 struct ssp_rng* rng, 425 double dir[3]) 426 { 427 struct ssol_material* mtl; 428 struct ssol_surface_fragment frag; 429 struct ssf_bsdf* bsdf = NULL; 430 double propagated = 0; 431 double wi[3], N[3], pdf; 432 int type = 0; 433 res_T res; 434 ASSERT(pt && in_medium && out_medium && rng && dir); 435 436 /* TODO ensure that if `prim' was sampled, then the surface fragment setup 437 * remains valid in *all* situations, i.e. even though the point primitive 438 * comes from a sampling operation. 439 * 440 * NOTE VF: actually a fragment generated from a RT or a sampled primitive is 441 * the same. Indeed it may be inconsistent only if the two kind of primitives 442 * does not have the same set of parameters. For triangulated meshes, the RT 443 * and sampled shape are the same and thus shared the same attribs. For 444 * punched surfaces, no attrib is defined on both representation. 445 * Consequently, it seems that there is no specific work to do to ensure the 446 * `surface_fragment_setup' consistency. */ 447 surface_fragment_setup(&frag, pt->pos, pt->dir, pt->N, &pt->prim, pt->uv); 448 449 /* Shade the surface fragment */ 450 mtl = point_get_material(pt); 451 452 res = material_create_bsdf(mtl, &frag, pt->wl, in_medium, 0, &bsdf); 453 if(res != RES_OK) goto error; 454 455 /* Perturbe the normal */ 456 material_shade_normal(mtl, &frag, pt->wl, N); 457 458 /* By convention, Star-SF assumes that incoming and reflected 459 * directions point outward the surface => negate incoming dir */ 460 d3_minus(wi, pt->dir); 461 462 if(d3_dot(wi, N) <= 0) { 463 propagated = 0; 464 } else { 465 double cos_dir_Ng; 466 propagated = ssf_bsdf_sample(bsdf, rng, wi, N, dir, &type, &pdf); 467 ASSERT(0 <= propagated && propagated <= 1); 468 469 /* Due to the perturbed normal, the sampled direction may point in the 470 * wrong direction wrt the sampled BSDF component. */ 471 cos_dir_Ng = d3_dot(frag.Ng, dir); 472 if((cos_dir_Ng > 0 && (type & SSF_TRANSMISSION)) 473 || (cos_dir_Ng < 0 && (type & SSF_REFLECTION))) { 474 propagated = 0; 475 } 476 } 477 pt->kabs_at_pt = (1 - propagated); 478 pt->outgoing_flux = pt->incoming_flux * propagated; 479 pt->outgoing_if_no_atm_loss = pt->incoming_if_no_atm_loss * propagated; 480 pt->outgoing_if_no_field_loss = point_is_receiver(pt) 481 ? pt->incoming_if_no_field_loss*propagated : pt->incoming_if_no_field_loss; 482 483 if(type & SSF_TRANSMISSION) { 484 material_get_next_medium(mtl, in_medium, out_medium); 485 } else { 486 ssol_medium_copy(out_medium, in_medium); 487 } 488 489 exit: 490 if(bsdf) SSF(bsdf_ref_put(bsdf)); 491 return res; 492 error: 493 goto exit; 494 } 495 496 static FINLINE void 497 point_hit_virtual 498 (struct point* pt, 499 const struct ssol_medium* in_medium, 500 struct ssol_medium* out_medium) 501 { 502 pt->kabs_at_pt = 0; 503 pt->outgoing_flux = pt->incoming_flux; 504 pt->outgoing_if_no_atm_loss = pt->incoming_if_no_atm_loss; 505 pt->outgoing_if_no_field_loss = pt->incoming_if_no_field_loss; 506 ssol_medium_copy(out_medium, in_medium); 507 } 508 509 static FINLINE int32_t 510 point_get_id(const struct point* pt) 511 { 512 uint32_t inst_id; 513 SSOL(instance_get_id(pt->inst, &inst_id)); 514 return pt->side == SSOL_FRONT ? (int32_t)inst_id : -(int32_t)inst_id; 515 } 516 517 /******************************************************************************* 518 * Helper functions 519 ******************************************************************************/ 520 static INLINE void 521 check_energy_conservation 522 (struct ssol_scene* scn, 523 struct ssol_estimator* estimator, 524 const int64_t nrealisations) 525 { 526 struct ssol_mc_global global; 527 double dni; 528 double dni_s, pot; 529 double cos, rcv, atm, other, shadow, miss; 530 double cos_err, rcv_err, atm_err, other_err, shadow_err, miss_err; 531 double err, max_loss; 532 ASSERT(scn && estimator); 533 534 if(RES_OK != ssol_estimator_get_mc_global(estimator, &global)) return; 535 if(RES_OK != ssol_sun_get_dni(scn->sun, &dni)) return; 536 537 /* Fetch data */ 538 cos = global.cos_factor.E; 539 rcv = global.absorbed_by_receivers.E; 540 atm = global.extinguished_by_atmosphere.E; 541 other = global.other_absorbed.E; 542 shadow = global.shadowed.E; 543 miss = global.missing.E; 544 cos_err = global.cos_factor.SE; 545 rcv_err = global.absorbed_by_receivers.SE; 546 atm_err = global.extinguished_by_atmosphere.SE; 547 other_err = global.other_absorbed.SE; 548 shadow_err = global.shadowed.SE; 549 miss_err = global.missing.SE; 550 551 /* Check energy conservation */ 552 dni_s = dni * scn->sampled_area; 553 pot = cos * dni_s; 554 err = dni_s * cos_err + rcv_err + atm_err + other_err + shadow_err + miss_err; 555 max_loss = 3 * err + (double)nrealisations * pot * DBL_EPSILON; 556 if(fabs(pot - (rcv + atm + other + shadow + miss)) > max_loss) 557 FATAL("error: the energy conservation property is not verified\n"); 558 } 559 560 /* Compute an empirical length of the path segment coming from/going to the 561 * infinite, wrt the scene bounding box */ 562 static INLINE double 563 compute_infinite_path_segment_extend(struct s3d_scene_view* view) 564 { 565 float lower[3], upper[3], size[3]; 566 ASSERT(view); 567 S3D(scene_view_get_aabb(view, lower, upper)); 568 f3_sub(size, upper, lower); 569 return MMAX(size[0], MMAX(size[1], size[2])) * 0.75; 570 } 571 572 static INLINE res_T 573 path_register_and_clear 574 (struct darray_path* paths, 575 struct path* path) 576 { 577 struct path* dst_path; 578 size_t ipath; 579 res_T res = RES_OK; 580 ASSERT(paths && path); 581 582 ipath = darray_path_size_get(paths); 583 res = darray_path_resize(paths, ipath + 1); 584 if(res != RES_OK) return res; 585 586 dst_path = darray_path_data_get(paths) + ipath; 587 return path_copy_and_clear(dst_path, path); 588 } 589 590 static res_T 591 accum_mc_receivers_1side 592 (struct mc_receiver_1side* dst, 593 struct mc_receiver_1side* src) 594 { 595 struct htable_shape2mc_iterator it_shape, end_shape; 596 res_T res = RES_OK; 597 ASSERT(dst && src); 598 599 #define ACCUM_ALL { \ 600 ACCUM_WEIGHT(incoming_flux); \ 601 ACCUM_WEIGHT(incoming_if_no_atm_loss); \ 602 ACCUM_WEIGHT(incoming_lost_in_field); \ 603 ACCUM_WEIGHT(incoming_lost_in_atmosphere); \ 604 ACCUM_WEIGHT(incoming_if_no_field_loss); \ 605 ACCUM_WEIGHT(absorbed_flux); \ 606 ACCUM_WEIGHT(absorbed_if_no_atm_loss); \ 607 ACCUM_WEIGHT(absorbed_if_no_field_loss); \ 608 ACCUM_WEIGHT(absorbed_lost_in_field); \ 609 ACCUM_WEIGHT(absorbed_lost_in_atmosphere); \ 610 } (void)0 611 612 #define ACCUM_WEIGHT(Name) mc_data_accum(&dst->Name, &src->Name) 613 ACCUM_ALL; 614 #undef ACCUM_WEIGHT 615 616 /* Merge the per shape MC */ 617 htable_shape2mc_begin(&src->shape2mc, &it_shape); 618 htable_shape2mc_end(&src->shape2mc, &end_shape); 619 while(!htable_shape2mc_iterator_eq(&it_shape, &end_shape)) { 620 struct htable_prim2mc_iterator it_prim, end_prim; 621 const struct ssol_shape* shape = *htable_shape2mc_iterator_key_get(&it_shape); 622 struct mc_shape_1side* mc_shape1_src; 623 struct mc_shape_1side* mc_shape1_dst; 624 625 mc_shape1_src = htable_shape2mc_iterator_data_get(&it_shape); 626 627 res = mc_receiver_1side_get_mc_shape(dst, shape, &mc_shape1_dst); 628 if(res != RES_OK) goto error; 629 630 /* Merge the per primitive MC */ 631 htable_prim2mc_begin(&mc_shape1_src->prim2mc, &it_prim); 632 htable_prim2mc_end(&mc_shape1_src->prim2mc, &end_prim); 633 while(!htable_prim2mc_iterator_eq(&it_prim, &end_prim)) { 634 const unsigned iprim = *htable_prim2mc_iterator_key_get(&it_prim); 635 struct mc_primitive_1side* mc_prim1_src; 636 struct mc_primitive_1side* mc_prim1_dst; 637 638 mc_prim1_src = htable_prim2mc_iterator_data_get(&it_prim); 639 640 res = mc_shape_1side_get_mc_primitive(mc_shape1_dst, iprim, &mc_prim1_dst); 641 if(res != RES_OK) goto error; 642 643 #define ACCUM_WEIGHT(Name) \ 644 mc_data_accum(&mc_prim1_dst->Name, &mc_prim1_src->Name) 645 ACCUM_ALL; 646 #undef ACCUM_WEIGHT 647 648 htable_prim2mc_iterator_next(&it_prim); 649 } 650 htable_shape2mc_iterator_next(&it_shape); 651 } 652 #undef ACCUM_ALL 653 654 exit: 655 return res; 656 error: 657 goto exit; 658 } 659 660 static res_T 661 accum_mc_sampled(struct mc_sampled* dst, struct mc_sampled* src) 662 { 663 struct htable_receiver_iterator it, end; 664 struct mc_receiver mc_rcv_null; 665 res_T res = RES_OK; 666 ASSERT(dst && src); 667 668 mc_receiver_init(NULL, &mc_rcv_null); 669 670 #define ACCUM_WEIGHT(Name) mc_data_accum(&dst->Name, &src->Name) 671 ACCUM_WEIGHT(cos_factor); 672 ACCUM_WEIGHT(shadowed); 673 #undef ACCUM_WEIGHT 674 675 dst->nb_samples += src->nb_samples; 676 677 /* dst->by_receiver += src->by_receiver; */ 678 htable_receiver_begin(&src->mc_rcvs, &it); 679 htable_receiver_end(&src->mc_rcvs, &end); 680 while(!htable_receiver_iterator_eq(&it, &end)) { 681 struct mc_receiver* src_mc_rcv = htable_receiver_iterator_data_get(&it); 682 const struct ssol_instance* inst = *htable_receiver_iterator_key_get(&it); 683 struct mc_receiver* dst_mc_rcv = htable_receiver_find(&dst->mc_rcvs, &inst); 684 htable_receiver_iterator_next(&it); 685 686 if(!dst_mc_rcv) { 687 res = htable_receiver_set(&dst->mc_rcvs, &inst, &mc_rcv_null); 688 if(res != RES_OK) goto error; 689 dst_mc_rcv = htable_receiver_find(&dst->mc_rcvs, &inst); 690 } 691 692 if(inst->receiver_mask & (int)SSOL_FRONT) { 693 res = accum_mc_receivers_1side(&dst_mc_rcv->front, &src_mc_rcv->front); 694 if(res != RES_OK) goto error; 695 } 696 if(inst->receiver_mask & (int)SSOL_BACK) { 697 res = accum_mc_receivers_1side(&dst_mc_rcv->back, &src_mc_rcv->back); 698 if(res != RES_OK) goto error; 699 } 700 } 701 exit: 702 mc_receiver_release(&mc_rcv_null); 703 return res; 704 error: 705 goto exit; 706 } 707 708 static res_T 709 update_mc 710 (struct point* pt, 711 const size_t irealisation, 712 struct thread_context* thread_ctx) 713 { 714 struct mc_receiver_1side* mc_rcv1 = NULL; 715 struct mc_receiver_1side* mc_samp_x_rcv1 = NULL; 716 res_T res = RES_OK; 717 ASSERT(pt && thread_ctx && point_is_receiver(pt)); 718 719 #define ACCUM_WEIGHT(Name, W)\ 720 mc_data_add_weight(&thread_ctx->Name, irealisation, W) 721 ACCUM_WEIGHT(absorbed_by_receivers, pt->incoming_flux - pt->outgoing_flux); 722 pt->energy_loss -= (pt->incoming_flux - pt->outgoing_flux); 723 #undef ACCUM_WEIGHT 724 725 /* Per receiver MC accumulation */ 726 res = get_mc_receiver_1side(&thread_ctx->mc_rcvs, pt->inst, pt->side, &mc_rcv1); 727 if(res != RES_OK) goto error; 728 729 #define ACCUM_ALL { \ 730 ACCUM_WEIGHT(incoming_flux, pt->incoming_flux); \ 731 ACCUM_WEIGHT(incoming_if_no_atm_loss, pt->incoming_if_no_atm_loss); \ 732 ACCUM_WEIGHT(incoming_if_no_field_loss, pt->incoming_if_no_field_loss); \ 733 ACCUM_WEIGHT(incoming_lost_in_field, \ 734 pt->incoming_if_no_field_loss - pt->incoming_flux); \ 735 ACCUM_WEIGHT(incoming_lost_in_atmosphere, \ 736 pt->incoming_if_no_atm_loss - pt->incoming_flux); \ 737 ACCUM_WEIGHT(absorbed_flux, pt->incoming_flux * pt->kabs_at_pt); \ 738 ACCUM_WEIGHT(absorbed_if_no_atm_loss, \ 739 pt->incoming_if_no_atm_loss * pt->kabs_at_pt); \ 740 ACCUM_WEIGHT(absorbed_if_no_field_loss, \ 741 pt->incoming_if_no_field_loss * pt->kabs_at_pt); \ 742 ACCUM_WEIGHT(absorbed_lost_in_field, \ 743 (pt->incoming_if_no_field_loss - pt->incoming_flux) * pt->kabs_at_pt); \ 744 ACCUM_WEIGHT(absorbed_lost_in_atmosphere, \ 745 (pt->incoming_if_no_atm_loss - pt->incoming_flux) * pt->kabs_at_pt); \ 746 } (void)0 747 748 #define ACCUM_WEIGHT(Name, W) \ 749 mc_data_add_weight(&mc_rcv1->Name, irealisation, W) 750 ACCUM_ALL; 751 #undef ACCUM_WEIGHT 752 753 /* Per-sampled/receiver MC accumulation */ 754 res = mc_sampled_get_mc_receiver_1side 755 (pt->mc_samp, pt->inst, pt->side, &mc_samp_x_rcv1); 756 if(res != RES_OK) goto error; 757 758 #define ACCUM_WEIGHT(Name, W) \ 759 mc_data_add_weight(&mc_samp_x_rcv1->Name, irealisation, W) 760 ACCUM_ALL; 761 #undef ACCUM_WEIGHT 762 763 /* Per primitive receiver MC accumulation */ 764 if(pt->inst->receiver_per_primitive) { 765 struct mc_shape_1side* mc_shape1; 766 struct mc_primitive_1side* mc_prim1; 767 768 res = mc_receiver_1side_get_mc_shape(mc_rcv1, pt->sshape->shape, &mc_shape1); 769 if(res != RES_OK) goto error; 770 771 res = mc_shape_1side_get_mc_primitive(mc_shape1, pt->prim.prim_id, &mc_prim1); 772 if(res != RES_OK) goto error; 773 774 #define ACCUM_WEIGHT(Name, W) \ 775 mc_data_add_weight(&mc_prim1->Name, irealisation, W) 776 ACCUM_ALL; 777 #undef ACCUM_WEIGHT 778 } 779 #undef ACCUM_ALL 780 781 exit: 782 return res; 783 error: 784 goto exit; 785 } 786 787 static void 788 apply_factor_mc_receiver_1side 789 (struct mc_receiver_1side* rcv, 790 const size_t irealisation, 791 const double factor) 792 { 793 mc_data_apply_factor(&rcv->incoming_flux, irealisation, factor); 794 mc_data_apply_factor(&rcv->incoming_if_no_atm_loss, irealisation, factor); 795 mc_data_apply_factor(&rcv->incoming_if_no_field_loss, irealisation, factor); 796 mc_data_apply_factor(&rcv->incoming_lost_in_field, irealisation, factor); 797 mc_data_apply_factor(&rcv->incoming_lost_in_atmosphere, irealisation, factor); 798 mc_data_apply_factor(&rcv->absorbed_flux, irealisation, factor); 799 mc_data_apply_factor(&rcv->absorbed_if_no_atm_loss, irealisation, factor); 800 mc_data_apply_factor(&rcv->absorbed_if_no_field_loss, irealisation, factor); 801 mc_data_apply_factor(&rcv->absorbed_lost_in_field, irealisation, factor); 802 mc_data_apply_factor(&rcv->absorbed_lost_in_atmosphere, irealisation, factor); 803 } 804 805 static void 806 apply_factor_mc 807 (struct thread_context* thread_ctx, 808 const size_t irealisation, 809 const double factor) 810 { 811 struct htable_receiver_iterator r_it, r_end; 812 struct htable_sampled_iterator s_it, s_end; 813 814 /* Cancel global MC estimations */ 815 mc_data_apply_factor(&thread_ctx->cos_factor, irealisation, factor); 816 mc_data_apply_factor(&thread_ctx->extinguished_by_atmosphere, irealisation, factor); 817 mc_data_apply_factor(&thread_ctx->absorbed_by_receivers, irealisation, factor); 818 mc_data_apply_factor(&thread_ctx->other_absorbed, irealisation, factor); 819 mc_data_apply_factor(&thread_ctx->missing, irealisation, factor); 820 mc_data_apply_factor(&thread_ctx->shadowed, irealisation, factor); 821 822 /* Cancel receiver MC estimations */ 823 htable_receiver_begin(&thread_ctx->mc_rcvs, &r_it); 824 htable_receiver_end(&thread_ctx->mc_rcvs, &r_end); 825 while(!htable_receiver_iterator_eq(&r_it, &r_end)) { 826 struct mc_receiver* mc_rcv = htable_receiver_iterator_data_get(&r_it); 827 const struct ssol_instance* inst = *htable_receiver_iterator_key_get(&r_it); 828 829 htable_receiver_iterator_next(&r_it); 830 831 if(inst->receiver_mask & (int)SSOL_FRONT) { 832 apply_factor_mc_receiver_1side(&mc_rcv->front, irealisation, factor); 833 } 834 if(inst->receiver_mask & (int)SSOL_BACK) { 835 apply_factor_mc_receiver_1side(&mc_rcv->back, irealisation, factor); 836 } 837 } 838 /* Cancel sampled instance MC estimations */ 839 htable_sampled_begin(&thread_ctx->mc_samps, &s_it); 840 htable_sampled_end(&thread_ctx->mc_samps, &s_end); 841 while(!htable_sampled_iterator_eq(&s_it, &s_end)) { 842 struct mc_sampled* mc_samp = htable_sampled_iterator_data_get(&s_it); 843 htable_sampled_iterator_next(&s_it); 844 845 mc_data_apply_factor(&mc_samp->cos_factor, irealisation, factor); 846 mc_data_apply_factor(&mc_samp->shadowed, irealisation, factor); 847 848 /* dst->by_receiver += src->by_receiver; */ 849 htable_receiver_begin(&mc_samp->mc_rcvs, &r_it); 850 htable_receiver_end(&mc_samp->mc_rcvs, &r_end); 851 while(!htable_receiver_iterator_eq(&r_it, &r_end)) { 852 struct mc_receiver* mc_rcv = htable_receiver_iterator_data_get(&r_it); 853 const struct ssol_instance* inst = *htable_receiver_iterator_key_get(&r_it); 854 htable_receiver_iterator_next(&r_it); 855 856 if(inst->receiver_mask & (int)SSOL_FRONT) { 857 apply_factor_mc_receiver_1side(&mc_rcv->front, irealisation, factor); 858 } 859 if(inst->receiver_mask & (int)SSOL_BACK) { 860 apply_factor_mc_receiver_1side(&mc_rcv->back, irealisation, factor); 861 } 862 } 863 } 864 } 865 static void 866 cancel_mc 867 (struct thread_context* thread_ctx, 868 const size_t irealisation) 869 { 870 apply_factor_mc(thread_ctx, irealisation, 0); 871 } 872 873 static res_T 874 trace_radiative_path 875 (const size_t irealisation, /* Unique id of the realisation */ 876 struct thread_context* thread_ctx, 877 struct ssol_scene* scn, 878 struct s3d_scene_view* view_samp, 879 struct s3d_scene_view* view_rt, 880 struct ranst_sun_dir* ran_sun_dir, 881 struct ranst_sun_wl* ran_sun_wl, 882 const struct ssol_path_tracker* tracker) /* May be NULL */ 883 { 884 struct path path; 885 struct ssol_medium in_medium = SSOL_MEDIUM_VACUUM; 886 struct ssol_medium out_medium = SSOL_MEDIUM_VACUUM; 887 struct s3d_hit hit = S3D_HIT_NULL; 888 struct point pt = POINT_NULL; 889 float org[3], dir[3], range[2] = { 0, FLT_MAX }; 890 size_t depth = 0; 891 size_t roulette_interval, typical_max_depth; 892 int is_lit = 0; 893 int hit_a_receiver = 0; 894 int killed_by_roulette = 0; 895 res_T res = RES_OK; 896 ASSERT(thread_ctx && scn && view_samp && view_rt && ran_sun_dir && ran_sun_wl); 897 898 if(tracker) path_init(scn->dev->allocator, &path); 899 900 typical_max_depth = 16; /* This one could come through scn */ 901 roulette_interval = 4 * typical_max_depth; /* First roulette */ 902 903 /* Find a new starting point of the radiative random walk */ 904 res = point_init(&pt, scn, &thread_ctx->mc_samps, 905 view_samp, view_rt, ran_sun_dir, ran_sun_wl, thread_ctx->rng, 906 &in_medium, &is_lit); 907 if(res != RES_OK) goto error; 908 909 if(tracker) { 910 /* Add the first point of the starting segment */ 911 if(tracker->sun_ray_length > 0) { 912 double pos[3], wi[3]; 913 d3_minus(wi, pt.dir); 914 d3_muld(wi, wi, tracker->sun_ray_length); 915 d3_add(pos, pt.pos, wi); 916 res = path_add_vertex(&path, pos, scn->sun->dni); 917 if(res != RES_OK) goto error; 918 } 919 920 /* Register the init position onto the sampled geometry */ 921 res = path_add_vertex(&path, pt.pos, pt.initial_flux); 922 if(res != RES_OK) goto error; 923 } 924 925 #define ACCUM_WEIGHT(Res, W) mc_data_add_weight(&Res, irealisation, W) 926 927 if(!is_lit) { /* The starting point is not lit */ 928 ACCUM_WEIGHT(pt.mc_samp->shadowed, pt.initial_flux); 929 ACCUM_WEIGHT(thread_ctx->shadowed, pt.initial_flux); 930 pt.energy_loss -= pt.initial_flux; 931 if(tracker) path.type = SSOL_PATH_SHADOW; 932 } else { 933 /* Setup the ray as if it starts from the current point position in order 934 * to handle the points that start from a virtual material */ 935 f3_set_d3(org, pt.pos); 936 f3_set_d3(dir, pt.dir); 937 hit.distance = 0; /* first loop has no atmospheric extinction */ 938 939 for(;;) { /* Here we go for the radiative random walk */ 940 const int in_atm = media_ceq(&in_medium, &scn->air); 941 const int hit_receiver = point_is_receiver(&pt); 942 const int hit_virtual = pt.material->type == SSOL_MATERIAL_VIRTUAL; 943 int last_segment = 0; 944 int weight_is_zero = 0; 945 struct ray_data ray_data = RAY_DATA_NULL; 946 double trans = 1; 947 948 /* Compute medium extinction along the incoming segment. */ 949 if(hit.distance > 0) { 950 const double k_ext = ssol_data_get_value(&in_medium.extinction, pt.wl); 951 ASSERT(0 <= k_ext && k_ext <= 1); 952 if(k_ext > 0) { 953 trans = exp(-k_ext * hit.distance); 954 } 955 } 956 pt.incoming_flux = pt.prev_outgoing_flux * trans; 957 pt.incoming_if_no_atm_loss = in_atm ? 958 pt.prev_outgoing_if_no_atm_loss : pt.prev_outgoing_if_no_atm_loss * trans; 959 pt.incoming_if_no_field_loss = (!in_atm) ? 960 pt.prev_outgoing_if_no_field_loss : pt.prev_outgoing_if_no_field_loss * trans; 961 962 /* Compute interaction with material */ 963 if(hit_virtual) { 964 point_hit_virtual(&pt, &in_medium, &out_medium); 965 } else { 966 /* Modulate the point weights wrt its scattering functions and generate 967 * an outgoing direction and set out_medium accordingly */ 968 res = point_shade(&pt, &in_medium, &out_medium, thread_ctx->rng, pt.dir); 969 if(res != RES_OK) goto error; 970 } 971 972 /* If receiver update MC results */ 973 if(hit_receiver) { 974 hit_a_receiver = 1; 975 res = update_mc(&pt, irealisation, thread_ctx); 976 if(res != RES_OK) goto error; 977 } else { 978 ACCUM_WEIGHT(thread_ctx->other_absorbed, 979 pt.incoming_flux * pt.kabs_at_pt); 980 pt.energy_loss -= (pt.incoming_flux * pt.kabs_at_pt); 981 } 982 983 /* Stop the radiative random walk if no more flux */ 984 if(!pt.outgoing_flux) { 985 weight_is_zero = 1; 986 } 987 988 /* Setup new ray parameters */ 989 if(!weight_is_zero) { 990 if(hit_virtual) { 991 /* Note that for Virtual materials, the ray parameters 'org' & 'dir' 992 * are not updated to ensure that it pursues its traversal without any 993 * accuracy issue */ 994 range[0] = nextafterf(hit.distance, FLT_MAX); 995 range[1] = FLT_MAX; 996 } else { 997 f2(range, 0, FLT_MAX); 998 f3_set_d3(org, pt.pos); 999 f3_set_d3(dir, pt.dir); 1000 } 1001 1002 /* Trace the next ray */ 1003 ray_data.scn = scn; 1004 ray_data.prim_from = pt.prim; 1005 ray_data.inst_from = pt.inst; 1006 ray_data.sshape_from = pt.sshape; 1007 ray_data.side_from = pt.side; 1008 ray_data.discard_virtual_materials = 0; 1009 ray_data.reversed_ray = 0; 1010 ray_data.dst = FLT_MAX; 1011 S3D(scene_view_trace_ray(view_rt, org, dir, range, &ray_data, &hit)); 1012 if(S3D_HIT_NONE(&hit)) { /* The ray is lost! */ 1013 /* Add the point of the last path segment going to the infinite */ 1014 if(tracker && tracker->infinite_ray_length > 0) { 1015 double pos[3], wi[3]; 1016 d3_set_f3(wi, dir); 1017 d3_add(pos, pt.pos, d3_muld(wi, wi, tracker->infinite_ray_length)); 1018 res = path_add_vertex(&path, pos, pt.outgoing_flux); 1019 if (res != RES_OK) goto error; 1020 } 1021 last_segment = 1; /* Path reached its last segment */ 1022 1023 /* Check medium consistency. Note that one has to check `out_medium' - 1024 * and not `in_medium' - against the atmosphere since it is actually 1025 * the medium in which the ray was traced; at this step, `in_medium' is 1026 * still the medium of the previous path segment. */ 1027 if(!media_ceq(&out_medium, &scn->air)) { 1028 log_error(scn->dev, "Inconsistent medium description.\n"); 1029 res = RES_BAD_OP; 1030 goto error; 1031 } 1032 } 1033 } 1034 1035 /* Don't change prev_outgoing weigths nor record segment extinction until 1036 * a non-virtual material is hit or this segment is the last one. 1037 * This is because propagation is restarted from the same origin until 1038 * a non-virtual material is hit or no further hit can be found. */ 1039 if(weight_is_zero || last_segment || !hit_virtual) { 1040 const double absorbed = pt.prev_outgoing_flux - pt.incoming_flux; 1041 if(in_atm) { 1042 ACCUM_WEIGHT(thread_ctx->extinguished_by_atmosphere, absorbed); 1043 } else { 1044 ACCUM_WEIGHT(thread_ctx->other_absorbed, absorbed); 1045 } 1046 pt.energy_loss -= absorbed; 1047 1048 if(weight_is_zero || last_segment) { 1049 break; 1050 } 1051 pt.prev_outgoing_flux = pt.outgoing_flux; 1052 pt.prev_outgoing_if_no_atm_loss = pt.outgoing_if_no_atm_loss; 1053 pt.prev_outgoing_if_no_field_loss = pt.outgoing_if_no_field_loss; 1054 } 1055 1056 depth += !hit_virtual; 1057 if(depth > roulette_interval && depth % roulette_interval == 1) { 1058 /* This could be in an infinite path. To avoid to crash the app while 1059 * preserving MC weights we have to use a russian roulette: 1/2 1060 * probability to end the path now compensated by a 2x factor on 1061 * weights if the path eventually ends somewhere. We could have 1062 * written a more traditional russian roulette that relies on not 1063 * applying kabs VS setting weights=0, but this doesn't work with 1064 * kabs=0. The present code works even if weigths remain unchanged 1065 * along the path. */ 1066 double p = ssp_rng_canonical(thread_ctx->rng); 1067 if(p > 0.5) { 1068 pt.survivor_score += 1; /* This path survived once more */ 1069 roulette_interval = typical_max_depth; 1070 } else { 1071 cancel_mc(thread_ctx, irealisation); 1072 killed_by_roulette = 1; 1073 goto exit; /* break is not enough */ 1074 } 1075 } 1076 1077 /* Update the point */ 1078 point_update_from_hit(&pt, scn, org, dir, &hit, &ray_data); 1079 1080 if(tracker) { 1081 res = path_add_vertex(&path, pt.pos, pt.outgoing_flux); 1082 if (res != RES_OK) goto error; 1083 } 1084 1085 ssol_medium_copy(&in_medium, &out_medium); 1086 } 1087 /* Register the remaining flux as missing */ 1088 ACCUM_WEIGHT(thread_ctx->missing, pt.outgoing_flux); 1089 pt.energy_loss -= pt.outgoing_flux; 1090 1091 1092 if(tracker) { 1093 path.type = hit_a_receiver ? SSOL_PATH_SUCCESS : SSOL_PATH_MISSING; 1094 } 1095 } 1096 /* Now that the sample ends successfully, record MC weights */ 1097 ACCUM_WEIGHT(pt.mc_samp->cos_factor, pt.cos_factor); 1098 ACCUM_WEIGHT(thread_ctx->cos_factor, pt.cos_factor); 1099 #undef ACCUM_WEIGHT 1100 1101 /* Check conservation of energy at the realisation level */ 1102 ASSERT(((double)depth*DBL_EPSILON*10)*pt.initial_flux >= fabs(pt.energy_loss)); 1103 1104 /* this realisation accounts for many that where canceled */ 1105 if(pt.survivor_score) { 1106 const double factor = (double)(1 << pt.survivor_score); 1107 apply_factor_mc(thread_ctx, irealisation, factor); 1108 } 1109 1110 exit: 1111 if(tracker && !killed_by_roulette) { 1112 res_T tmp_res = path_register_and_clear(&thread_ctx->paths, &path); 1113 if(tmp_res != RES_OK && res == RES_OK) { 1114 res = tmp_res; 1115 goto error; 1116 } 1117 } 1118 ssol_medium_clear(&in_medium); 1119 ssol_medium_clear(&out_medium); 1120 if(tracker) path_release(&path); 1121 return res; 1122 error: 1123 if (tracker) { 1124 path.type = SSOL_PATH_ERROR; 1125 } 1126 goto exit; 1127 } 1128 1129 /******************************************************************************* 1130 * Exported functions 1131 ******************************************************************************/ 1132 res_T 1133 ssol_solve 1134 (struct ssol_scene* scn, 1135 const struct ssp_rng* rng_state, 1136 const size_t realisations_count, 1137 const size_t max_failed_count, 1138 const struct ssol_path_tracker* path_tracker, 1139 struct ssol_estimator** out_estimator) 1140 { 1141 struct htable_receiver_iterator r_it, r_end; 1142 struct htable_sampled_iterator s_it, s_end; 1143 struct s3d_scene_view* view_rt = NULL; 1144 struct s3d_scene_view* view_samp = NULL; 1145 struct ranst_sun_dir* ran_sun_dir = NULL; 1146 struct ranst_sun_wl* ran_sun_wl = NULL; 1147 struct darray_thread_ctx thread_ctxs; 1148 struct ssol_estimator* estimator = NULL; 1149 struct ssol_path_tracker tracker; 1150 struct ssp_rng_proxy* rng_proxy = NULL; 1151 int64_t nrealisations = 0; 1152 int64_t max_failures = 0; 1153 int nthreads = 0; 1154 int i = 0; 1155 ATOMIC mt_res = RES_OK; 1156 ATOMIC nfailures = 0; 1157 res_T res; 1158 ASSERT(nrealisations <= INT_MAX); 1159 1160 if(!scn || !rng_state || !realisations_count || !out_estimator) 1161 return RES_BAD_ARG; 1162 1163 darray_thread_ctx_init(scn->dev->allocator, &thread_ctxs); 1164 1165 /* CL compiler supports OpenMP parallel loop whose indices are signed. The 1166 * following line ensures that the unsigned number of realisations does not 1167 * overflow the realisation index. */ 1168 if(realisations_count > INT64_MAX || max_failed_count > INT64_MAX) { 1169 res = RES_BAD_ARG; 1170 goto error; 1171 } 1172 nrealisations = (int64_t)realisations_count; 1173 max_failures = (int64_t)max_failed_count; 1174 nthreads = (int)scn->dev->nthreads; 1175 1176 res = scene_check(scn, FUNC_NAME); 1177 if(res != RES_OK) goto error; 1178 1179 /* init air properties */ 1180 if(scn->atmosphere) 1181 ssol_data_copy(&scn->air.extinction, &scn->atmosphere->extinction); 1182 else 1183 ssol_data_copy(&scn->air.extinction, &SSOL_MEDIUM_VACUUM.extinction); 1184 1185 /* Create data structures shared by all threads */ 1186 res = scene_create_s3d_views(scn, &view_rt, &view_samp); 1187 if(res != RES_OK) goto error; 1188 res = sun_create_direction_distribution(scn->sun, &ran_sun_dir); 1189 if(res != RES_OK) goto error; 1190 res = sun_create_wavelength_distribution(scn->sun, &ran_sun_wl); 1191 if(res != RES_OK) goto error; 1192 1193 /* Create a RNG proxy from the submitted RNG state */ 1194 res = ssp_rng_proxy_create_from_rng 1195 (scn->dev->allocator, rng_state, scn->dev->nthreads, &rng_proxy); 1196 if(res != RES_OK) goto error; 1197 1198 /* Create the estimator */ 1199 res = estimator_create(scn->dev, scn, &estimator); 1200 if (res != RES_OK) goto error; 1201 1202 /* Create per thread data structures */ 1203 res = darray_thread_ctx_resize(&thread_ctxs, scn->dev->nthreads); 1204 if(res != RES_OK) goto error; 1205 FOR_EACH(i, 0, nthreads) { 1206 struct thread_context* ctx = darray_thread_ctx_data_get(&thread_ctxs)+i; 1207 res = thread_context_setup(ctx, rng_proxy, (size_t)i); 1208 if(res != RES_OK) goto error; 1209 } 1210 1211 /* Setup the path tracker */ 1212 if(path_tracker) { 1213 tracker = *path_tracker; 1214 if(tracker.sun_ray_length < 0 || tracker.infinite_ray_length < 0) { 1215 const double extend = compute_infinite_path_segment_extend(view_rt); 1216 if(tracker.sun_ray_length < 0) tracker.sun_ray_length = extend; 1217 if(tracker.infinite_ray_length < 0) tracker.infinite_ray_length = extend; 1218 } 1219 path_tracker = &tracker; 1220 } 1221 1222 /* Launch the parallel MC estimation */ 1223 #pragma omp parallel for schedule(static) 1224 for(i = 0; i < nrealisations; ++i) { 1225 struct thread_context* thread_ctx; 1226 const int ithread = omp_get_thread_num(); 1227 res_T res_local; 1228 1229 if(ATOMIC_GET(&mt_res) != RES_OK) continue; /* An error occured */ 1230 1231 /* Fetch per thread data */ 1232 thread_ctx = darray_thread_ctx_data_get(&thread_ctxs) + ithread; 1233 1234 /* Execute a MC experiment */ 1235 res_local = trace_radiative_path((size_t)i, thread_ctx, 1236 scn, view_samp, view_rt, ran_sun_dir, ran_sun_wl, path_tracker); 1237 if(res_local != RES_OK) { 1238 /* Cancel partial MC results */ 1239 cancel_mc(thread_ctx, (size_t)i); 1240 } 1241 if(res_local == RES_BAD_OP) { 1242 if(ATOMIC_INCR(&nfailures) >= max_failures) { 1243 log_error(scn->dev, "Too many unexpected radiative paths.\n"); 1244 ATOMIC_SET(&mt_res, res_local); 1245 } 1246 } else if(res_local != RES_OK) { 1247 ATOMIC_SET(&mt_res, res_local); 1248 } 1249 if(res_local != RES_OK) continue; 1250 thread_ctx->realisation_count++; 1251 } 1252 estimator->failed_count = (size_t)nfailures; 1253 1254 /* Merge per thread global MC estimations */ 1255 FOR_EACH(i, 0, nthreads) { 1256 struct thread_context* thread_ctx; 1257 thread_ctx = darray_thread_ctx_data_get(&thread_ctxs)+i; 1258 #define ACCUM_WEIGHT(Name) \ 1259 mc_data_accum(&estimator->Name, &thread_ctx->Name) 1260 ACCUM_WEIGHT(cos_factor); 1261 ACCUM_WEIGHT(absorbed_by_receivers); 1262 ACCUM_WEIGHT(shadowed); 1263 ACCUM_WEIGHT(missing); 1264 ACCUM_WEIGHT(extinguished_by_atmosphere); 1265 ACCUM_WEIGHT(other_absorbed); 1266 estimator->realisation_count += thread_ctx->realisation_count; 1267 #undef ACCUM_WEIGHT 1268 } 1269 1270 /* Merge per thread receiver MC estimations */ 1271 htable_receiver_begin(&estimator->mc_receivers, &r_it); 1272 htable_receiver_end(&estimator->mc_receivers, &r_end); 1273 while(!htable_receiver_iterator_eq(&r_it, &r_end)) { 1274 struct mc_receiver* mc_rcv = htable_receiver_iterator_data_get(&r_it); 1275 const struct ssol_instance* inst = *htable_receiver_iterator_key_get(&r_it); 1276 htable_receiver_iterator_next(&r_it); 1277 1278 FOR_EACH(i, 0, nthreads) { 1279 struct thread_context* thread_ctx; 1280 struct mc_receiver* mc_rcv_thread; 1281 1282 thread_ctx = darray_thread_ctx_data_get(&thread_ctxs) + i; 1283 mc_rcv_thread = htable_receiver_find(&thread_ctx->mc_rcvs, &inst); 1284 if(!mc_rcv_thread) continue; /* Receiver was not visited in this thread */ 1285 1286 if(inst->receiver_mask & (int)SSOL_FRONT) { 1287 res = accum_mc_receivers_1side(&mc_rcv->front, &mc_rcv_thread->front); 1288 if(res != RES_OK) goto error; 1289 } 1290 if(inst->receiver_mask & (int)SSOL_BACK) { 1291 res = accum_mc_receivers_1side(&mc_rcv->back, &mc_rcv_thread->back); 1292 if(res != RES_OK) goto error; 1293 } 1294 } 1295 } 1296 1297 /* Merge per thread sampled instance MC estimations */ 1298 htable_sampled_begin(&estimator->mc_sampled, &s_it); 1299 htable_sampled_end(&estimator->mc_sampled, &s_end); 1300 while(!htable_sampled_iterator_eq(&s_it, &s_end)) { 1301 struct mc_sampled* mc_samp = htable_sampled_iterator_data_get(&s_it); 1302 const struct ssol_instance* inst = *htable_sampled_iterator_key_get(&s_it); 1303 htable_sampled_iterator_next(&s_it); 1304 1305 FOR_EACH(i, 0, nthreads) { 1306 struct thread_context* thread_ctx; 1307 struct mc_sampled* mc_samp_thread; 1308 1309 thread_ctx = darray_thread_ctx_data_get(&thread_ctxs) + i; 1310 mc_samp_thread = htable_sampled_find(&thread_ctx->mc_samps, &inst); 1311 if(!mc_samp_thread) continue; /* Instance was not sampled in this thread */ 1312 1313 res = accum_mc_sampled(mc_samp, mc_samp_thread); 1314 if(res != RES_OK) goto error; 1315 } 1316 } 1317 1318 /* Merge per thread tracked paths */ 1319 if(path_tracker) { 1320 FOR_EACH(i, 0, nthreads) { 1321 struct thread_context* thread_ctx; 1322 size_t ipath, npaths; 1323 1324 thread_ctx = darray_thread_ctx_data_get(&thread_ctxs) + i; 1325 npaths = darray_path_size_get(&thread_ctx->paths); 1326 FOR_EACH(ipath, 0, npaths) { 1327 struct path* path; 1328 path = darray_path_data_get(&thread_ctx->paths) + ipath; 1329 res = path_register_and_clear(&estimator->paths, path); 1330 if(res != RES_OK) goto error; 1331 } 1332 } 1333 } 1334 1335 estimator->sampled_area = scn->sampled_area; 1336 1337 res = estimator_save_rng_state(estimator, rng_proxy); 1338 if(res != RES_OK) goto error; 1339 1340 if(mt_res != RES_OK) res = (res_T)mt_res; 1341 1342 #ifndef NDEBUG 1343 check_energy_conservation(scn, estimator, nrealisations); 1344 #endif 1345 1346 exit: 1347 darray_thread_ctx_release(&thread_ctxs); 1348 if(view_rt) S3D(scene_view_ref_put(view_rt)); 1349 if(view_samp) S3D(scene_view_ref_put(view_samp)); 1350 if(ran_sun_dir) ranst_sun_dir_ref_put(ran_sun_dir); 1351 if(ran_sun_wl) ranst_sun_wl_ref_put(ran_sun_wl); 1352 if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy)); 1353 if(out_estimator) *out_estimator = estimator; 1354 return res; 1355 error: 1356 if(estimator) { 1357 SSOL(estimator_ref_put(estimator)); 1358 estimator = NULL; 1359 } 1360 goto exit; 1361 } 1362