ssol_estimator.c (12995B)
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 #include "ssol.h" 18 #include "ssol_c.h" 19 #include "ssol_scene_c.h" 20 #include "ssol_estimator_c.h" 21 #include "ssol_device_c.h" 22 #include "ssol_instance_c.h" 23 24 #include <rsys/double3.h> 25 #include <rsys/mem_allocator.h> 26 #include <rsys/ref_count.h> 27 #include <rsys/rsys.h> 28 29 #include <star/ssp.h> 30 31 #include <math.h> 32 33 /******************************************************************************* 34 * Helper functions 35 ******************************************************************************/ 36 static res_T 37 create_mc_receivers 38 (struct ssol_estimator* estimator, 39 struct ssol_scene* scene) 40 { 41 struct htable_instance_iterator it, end; 42 struct mc_receiver mc_rcv_null; 43 struct mc_sampled mc_samp_null; 44 res_T res = RES_OK; 45 ASSERT(scene && estimator); 46 47 htable_instance_begin(&scene->instances_rt, &it); 48 htable_instance_end(&scene->instances_rt, &end); 49 50 mc_receiver_init(estimator->dev->allocator, &mc_rcv_null); 51 mc_sampled_init(estimator->dev->allocator, &mc_samp_null); 52 53 while(!htable_instance_iterator_eq(&it, &end)) { 54 const struct ssol_instance* inst = *htable_instance_iterator_data_get(&it); 55 htable_instance_iterator_next(&it); 56 57 if(inst->receiver_mask) { 58 res = htable_receiver_set(&estimator->mc_receivers, &inst, &mc_rcv_null); 59 if(res != RES_OK) goto error; 60 } 61 if(inst->sample) { 62 res = htable_sampled_set(&estimator->mc_sampled, &inst, &mc_samp_null); 63 if(res != RES_OK) goto error; 64 } 65 } 66 exit: 67 mc_receiver_release(&mc_rcv_null); 68 mc_sampled_release(&mc_samp_null); 69 return res; 70 error: 71 htable_receiver_clear(&estimator->mc_receivers); 72 htable_sampled_clear(&estimator->mc_sampled); 73 goto exit; 74 } 75 76 static void 77 estimator_release(ref_T* ref) 78 { 79 struct ssol_device* dev; 80 struct ssol_estimator* estimator = 81 CONTAINER_OF(ref, struct ssol_estimator, ref); 82 ASSERT(ref); 83 dev = estimator->dev; 84 htable_receiver_release(&estimator->mc_receivers); 85 htable_sampled_release(&estimator->mc_sampled); 86 darray_path_release(&estimator->paths); 87 if(estimator->rng) SSP(rng_ref_put(estimator->rng)); 88 ASSERT(dev && dev->allocator); 89 MEM_RM(dev->allocator, estimator); 90 SSOL(device_ref_put(dev)); 91 } 92 93 /******************************************************************************* 94 * Exported function 95 ******************************************************************************/ 96 res_T 97 ssol_estimator_ref_get(struct ssol_estimator* estimator) 98 { 99 if(!estimator) return RES_BAD_ARG; 100 ref_get(&estimator->ref); 101 return RES_OK; 102 } 103 104 res_T 105 ssol_estimator_ref_put(struct ssol_estimator* estimator) 106 { 107 if(!estimator) return RES_BAD_ARG; 108 ref_put(&estimator->ref, estimator_release); 109 return RES_OK; 110 } 111 112 res_T 113 ssol_estimator_get_mc_global 114 (struct ssol_estimator* estimator, 115 struct ssol_mc_global* global) 116 { 117 if(!estimator || !global) return RES_BAD_ARG; 118 #define SETUP_MC_RESULT(Name) { \ 119 const double N = (double)estimator->realisation_count; \ 120 struct mc_data* data = &estimator->Name; \ 121 double weight, sqr_weight; \ 122 mc_data_get(data, &weight, &sqr_weight); \ 123 global->Name.E = weight / N; \ 124 global->Name.V = sqr_weight / N - global->Name.E*global->Name.E; \ 125 global->Name.V = global->Name.V > 0 ? global->Name.V : 0; \ 126 global->Name.SE = sqrt(global->Name.V / N); \ 127 } (void)0 128 SETUP_MC_RESULT(cos_factor); 129 SETUP_MC_RESULT(absorbed_by_receivers); 130 SETUP_MC_RESULT(shadowed); 131 SETUP_MC_RESULT(missing); 132 SETUP_MC_RESULT(extinguished_by_atmosphere); 133 SETUP_MC_RESULT(other_absorbed); 134 #undef SETUP_MC_RESULT 135 return RES_OK; 136 } 137 138 res_T 139 ssol_estimator_get_mc_sampled_x_receiver 140 (struct ssol_estimator* estimator, 141 const struct ssol_instance* samp_instance, 142 const struct ssol_instance* recv_instance, 143 const enum ssol_side_flag side, 144 struct ssol_mc_receiver* rcv) 145 { 146 struct mc_sampled* mc_samp = NULL; 147 struct mc_receiver* mc_rcv = NULL; 148 struct mc_receiver_1side* mc_rcv1 = NULL; 149 150 if(!estimator || !samp_instance || !recv_instance || !rcv 151 || (side != SSOL_BACK && side != SSOL_FRONT) 152 || !samp_instance->sample 153 || !(recv_instance->receiver_mask & (int)side)) 154 return RES_BAD_ARG; 155 156 memset(rcv, 0, sizeof(rcv[0])); 157 158 mc_samp = htable_sampled_find(&estimator->mc_sampled, &samp_instance); 159 if(!mc_samp) { 160 /* The sampled instance has no MC estimation */ 161 return RES_BAD_ARG; 162 } 163 164 mc_rcv = htable_receiver_find(&mc_samp->mc_rcvs, &recv_instance); 165 if(!mc_rcv) { 166 /* No radiative path starting from the sampled instance reaches the receiver 167 * instance. */ 168 return RES_OK; 169 } 170 171 mc_rcv1 = side == SSOL_FRONT ? &mc_rcv->front : &mc_rcv->back; 172 #define SETUP_MC_RESULT(Name) { \ 173 const double N = (double)estimator->realisation_count; \ 174 struct mc_data* data = &mc_rcv1->Name; \ 175 double weight, sqr_weight; \ 176 mc_data_get(data, &weight, &sqr_weight); \ 177 rcv->Name.E = weight / N; \ 178 rcv->Name.V = sqr_weight / N - rcv->Name.E*rcv->Name.E; \ 179 rcv->Name.V = rcv->Name.V > 0 ? rcv->Name.V : 0; \ 180 rcv->Name.SE = sqrt(rcv->Name.V / N); \ 181 } (void)0 182 SETUP_MC_RESULT(incoming_flux); 183 SETUP_MC_RESULT(incoming_if_no_atm_loss); 184 SETUP_MC_RESULT(incoming_if_no_field_loss); 185 SETUP_MC_RESULT(incoming_lost_in_field); 186 SETUP_MC_RESULT(incoming_lost_in_atmosphere); 187 SETUP_MC_RESULT(absorbed_flux); 188 SETUP_MC_RESULT(absorbed_if_no_atm_loss); 189 SETUP_MC_RESULT(absorbed_if_no_field_loss); 190 SETUP_MC_RESULT(absorbed_lost_in_field); 191 SETUP_MC_RESULT(absorbed_lost_in_atmosphere); 192 #undef SETUP_MC_RESULT 193 rcv->mc__ = mc_rcv1; 194 rcv->N__ = mc_samp->nb_samples; 195 return RES_OK; 196 } 197 198 res_T 199 ssol_estimator_get_realisation_count 200 (const struct ssol_estimator* estimator, size_t* count) 201 { 202 if (!estimator || !count) return RES_BAD_ARG; 203 *count = estimator->realisation_count; 204 return RES_OK; 205 } 206 207 res_T 208 ssol_estimator_get_failed_count 209 (const struct ssol_estimator* estimator, size_t* count) 210 { 211 if (!estimator || !count) return RES_BAD_ARG; 212 *count = estimator->failed_count; 213 return RES_OK; 214 } 215 216 res_T 217 ssol_estimator_get_sampled_area 218 (const struct ssol_estimator* estimator, double* area) 219 { 220 if(!estimator || !area) return RES_BAD_ARG; 221 *area = estimator->sampled_area; 222 return RES_OK; 223 } 224 225 res_T 226 ssol_estimator_get_sampled_count 227 (const struct ssol_estimator* estimator, size_t* count) 228 { 229 if (!estimator || !count) return RES_BAD_ARG; 230 *count = htable_sampled_size_get(&estimator->mc_sampled); 231 return RES_OK; 232 } 233 234 res_T 235 ssol_estimator_get_mc_sampled 236 (struct ssol_estimator* estimator, 237 const struct ssol_instance* samp_instance, 238 struct ssol_mc_sampled* sampled) 239 { 240 struct mc_sampled* mc = NULL; 241 if (!estimator || !samp_instance || !sampled) return RES_BAD_ARG; 242 mc = htable_sampled_find(&estimator->mc_sampled, &samp_instance); 243 if(!mc) return RES_BAD_ARG; 244 sampled->nb_samples = mc->nb_samples; 245 #define SETUP_MC_RESULT(Name, Count) { \ 246 const double N = (double)(Count); \ 247 struct mc_data* data = &mc->Name; \ 248 double weight, sqr_weight; \ 249 mc_data_get(data, &weight, &sqr_weight); \ 250 sampled->Name.E = weight / N; \ 251 sampled->Name.V = sqr_weight/N - sampled->Name.E*sampled->Name.E; \ 252 sampled->Name.V = sampled->Name.V > 0 ? sampled->Name.V : 0; \ 253 sampled->Name.SE = sqrt(sampled->Name.V / N); \ 254 } (void)0 255 SETUP_MC_RESULT(cos_factor, sampled->nb_samples); 256 SETUP_MC_RESULT(shadowed, estimator->realisation_count); 257 #undef SETUP_MC_RESULT 258 return RES_OK; 259 } 260 261 res_T 262 ssol_estimator_get_rng_state 263 (const struct ssol_estimator* estimator, const struct ssp_rng** rng_state) 264 { 265 if(!estimator || !rng_state) return RES_BAD_ARG; 266 *rng_state = estimator->rng; 267 return RES_OK; 268 } 269 270 res_T 271 ssol_estimator_get_tracked_paths_count 272 (const struct ssol_estimator* estimator, size_t* npaths) 273 { 274 if(!estimator || !npaths) return RES_BAD_ARG; 275 *npaths = darray_path_size_get(&estimator->paths); 276 return RES_OK; 277 } 278 279 res_T 280 ssol_estimator_get_tracked_path 281 (const struct ssol_estimator* estimator, 282 const size_t ipath, 283 struct ssol_path* path) 284 { 285 if(!estimator || ipath >= darray_path_size_get(&estimator->paths) || !path) 286 return RES_BAD_ARG; 287 path->path__ = darray_path_cdata_get(&estimator->paths) + ipath; 288 return RES_OK; 289 } 290 291 res_T 292 ssol_path_get_vertices_count(const struct ssol_path* path, size_t* nvertices) 293 { 294 const struct path* p; 295 if(!path || !nvertices) return RES_BAD_ARG; 296 p = path->path__; 297 *nvertices = darray_path_vertex_size_get(&p->vertices); 298 return RES_OK; 299 } 300 301 res_T 302 ssol_path_get_vertex 303 (const struct ssol_path* path, 304 const size_t ivertex, 305 struct ssol_path_vertex* vertex) 306 { 307 const struct path* p; 308 if(!path || !vertex) return RES_BAD_ARG; 309 p = path->path__; 310 if(ivertex >= darray_path_vertex_size_get(&p->vertices)) return RES_BAD_ARG; 311 *vertex = darray_path_vertex_cdata_get(&p->vertices)[ivertex]; 312 return RES_OK; 313 } 314 315 res_T 316 ssol_path_get_type(const struct ssol_path* path, enum ssol_path_type* type) 317 { 318 ASSERT(path && type); 319 if(!path || !type) return RES_BAD_ARG; 320 *type = ((struct path*)path->path__)->type; 321 return RES_OK; 322 } 323 324 /******************************************************************************* 325 * Local function 326 ******************************************************************************/ 327 res_T 328 estimator_create 329 (struct ssol_device* dev, 330 struct ssol_scene* scene, 331 struct ssol_estimator** out_estimator) 332 { 333 struct ssol_estimator* estimator = NULL; 334 res_T res = RES_OK; 335 336 if (!dev || !scene || !out_estimator) { 337 res = RES_BAD_ARG; 338 goto error; 339 } 340 341 estimator = MEM_CALLOC(dev->allocator, 1, sizeof(struct ssol_estimator)); 342 if(!estimator) { 343 res = RES_MEM_ERR; 344 goto error; 345 } 346 347 htable_receiver_init(dev->allocator, &estimator->mc_receivers); 348 htable_sampled_init(dev->allocator, &estimator->mc_sampled); 349 darray_path_init(dev->allocator, &estimator->paths); 350 SSOL(device_ref_get(dev)); 351 estimator->dev = dev; 352 ref_init(&estimator->ref); 353 354 res = create_mc_receivers(estimator, scene); 355 if(res != RES_OK) goto error; 356 357 exit: 358 if(out_estimator) *out_estimator = estimator; 359 return res; 360 error: 361 if(estimator) { 362 SSOL(estimator_ref_put(estimator)); 363 estimator = NULL; 364 } 365 goto exit; 366 } 367 368 res_T 369 estimator_save_rng_state 370 (struct ssol_estimator* estimator, const struct ssp_rng_proxy* proxy) 371 { 372 enum ssp_rng_type rng_type; 373 FILE* stream = NULL; 374 res_T res = RES_OK; 375 ASSERT(estimator && proxy); 376 377 /* Release the previous rng state */ 378 if(estimator->rng) { 379 SSP(rng_ref_put(estimator->rng)); 380 estimator->rng = NULL; 381 } 382 383 stream = tmpfile(); 384 if(!stream) { 385 log_error(estimator->dev, 386 "Could not open a stream to store the proxy RNG state.\n"); 387 res = RES_IO_ERR; 388 goto error; 389 } 390 391 SSP(rng_proxy_get_type(proxy, &rng_type)); 392 res = ssp_rng_create(estimator->dev->allocator, rng_type, &estimator->rng); 393 if(res != RES_OK) { 394 log_error(estimator->dev, 395 "Could not create the RNG to save the proxy RNG state.\n"); 396 goto error; 397 } 398 399 res = ssp_rng_proxy_write(proxy, stream); 400 if(res != RES_OK) { 401 log_error(estimator->dev, "Could not serialize the proxy RNG state.\n"); 402 goto error; 403 } 404 405 rewind(stream); 406 res = ssp_rng_read(estimator->rng, stream); 407 if(res != RES_OK) { 408 log_error(estimator->dev, "Could not save the proxy RNG state.\n"); 409 goto error; 410 } 411 412 exit: 413 if(stream) fclose(stream); 414 return res; 415 error: 416 if(estimator->rng) SSP(rng_ref_put(estimator->rng)); 417 goto exit; 418 } 419