sgf_realisation.h (11055B)
1 /* Copyright (C) 2021, 2024 |Meso|Star> (contact@meso-star.com) 2 * Copyright (C) 2015-2018 EDF S.A., France (syrthes-support@edf.fr) 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 /* 18 * Generate the gebhart_radiative_path_<2|3>d realisation functions with 19 * respect to the SGF_DIMENSIONALITY macro 20 */ 21 22 #ifndef SGF_DIMENSIONALITY 23 24 #ifndef SGF_REALISATION_H 25 #define SGF_REALISATION_H 26 27 #include "sgf_scene_c.h" 28 29 #include <rsys/float2.h> 30 #include <rsys/float3.h> 31 #include <rsys/hash_table.h> 32 33 #include <star/s2d.h> 34 #include <star/s3d.h> 35 #include <star/ssp.h> 36 37 /* How many self intersections are tested on the same ray before an error 38 * occurs */ 39 #define NSELF_HITS_MAX 10 40 41 /* Monte Carlo estimator */ 42 struct accum { 43 double radiative_flux; 44 double sqr_radiative_flux; 45 }; 46 47 /* Declare the htable_bounce data structure */ 48 #define HTABLE_NAME bounce 49 #define HTABLE_KEY unsigned /* Primitive ID */ 50 #define HTABLE_DATA double /* Weight */ 51 #include <rsys/hash_table.h> 52 53 /* Type of the realisation function. Return RES_BAD_OP on numerical issue. i.e. 54 * the radiative path is skipped */ 55 typedef res_T 56 (*gebhart_radiative_path_T) 57 (struct sgf_device* dev, 58 struct accum* accums, 59 struct accum* accum_infinity, 60 struct accum* accum_medium, 61 struct ssp_rng* rng, 62 struct htable_bounce* path, /* Store the path */ 63 const double absorption_coef, /* In m^-1 */ 64 const size_t ispectral_band, 65 const size_t primitive_id, 66 struct sgf_scene* scn); 67 68 static FINLINE void 69 accum_weight(struct accum* accum, const double weight) 70 { 71 ASSERT(accum); 72 accum->radiative_flux += weight; 73 accum->sqr_radiative_flux += weight * weight; 74 } 75 76 /******************************************************************************* 77 * 2D helper functions 78 ******************************************************************************/ 79 static FINLINE void 80 primitive_get_normal_2d(struct s2d_primitive* prim, float normal[3]) 81 { 82 struct s2d_attrib attr; 83 const float s = 0; 84 ASSERT(prim && normal); 85 S2D(primitive_get_attrib(prim, S2D_GEOMETRY_NORMAL, s, &attr)); 86 ASSERT(attr.type == S2D_FLOAT2); 87 f2_normalize(normal, attr.value); 88 normal[2] = 0.f; 89 } 90 91 static FINLINE void 92 primitive_sample_position_2d 93 (struct s2d_primitive* prim, 94 struct ssp_rng* rng, 95 float position[3]) 96 { 97 struct s2d_attrib attr; 98 float s; 99 ASSERT(prim && position); 100 S2D(primitive_sample(prim, ssp_rng_canonical_float(rng), &s)); 101 S2D(primitive_get_attrib(prim, S2D_POSITION, s, &attr)); 102 ASSERT(attr.type == S2D_FLOAT2); 103 f2_set(position, attr.value); 104 position[2] = 0.f; 105 } 106 107 static FINLINE void 108 hit_get_normal_2d(struct s2d_hit* hit, float normal[3]) 109 { 110 ASSERT(hit && normal); 111 f2_normalize(normal, hit->normal); 112 normal[2] = 0.f; 113 } 114 115 /******************************************************************************* 116 * 3D helper functions 117 ******************************************************************************/ 118 static FINLINE void 119 primitive_get_normal_3d(struct s3d_primitive* prim, float normal[3]) 120 { 121 struct s3d_attrib attr; 122 const float st[2] = { 0.f, 0.f }; 123 ASSERT(prim && normal); 124 S3D(primitive_get_attrib(prim, S3D_GEOMETRY_NORMAL, st, &attr)); 125 ASSERT(attr.type == S3D_FLOAT3); 126 f3_normalize(normal, attr.value); 127 } 128 129 static FINLINE void 130 primitive_sample_position_3d 131 (struct s3d_primitive* prim, 132 struct ssp_rng* rng, 133 float position[3]) 134 { 135 struct s3d_attrib attr; 136 float st[2]; 137 float u, v; 138 ASSERT(prim && position); 139 140 u = ssp_rng_canonical_float(rng); 141 v = ssp_rng_canonical_float(rng); 142 S3D(primitive_sample(prim, u, v, st)); 143 S3D(primitive_get_attrib(prim, S3D_POSITION, st, &attr)); 144 ASSERT(attr.type == S3D_FLOAT3); 145 f3_set(position, attr.value); 146 } 147 148 static FINLINE void 149 hit_get_normal_3d(struct s3d_hit* hit, float normal[3]) 150 { 151 ASSERT(hit && normal); 152 f3_normalize(normal, hit->normal); 153 } 154 155 #endif /* SGF_REALISATION_H */ 156 157 #else /* !SGF_DIMENSIONALITY */ 158 159 #if SGF_DIMENSIONALITY == 2 160 /* Wrap */ 161 #define sXd(Name) CONCAT(s2d_, Name) 162 #define SXD(Name) CONCAT(S2D_, Name) 163 #define Xd(Name) CONCAT(Name, _2d) 164 165 #elif SGF_DIMENSIONALITY == 3 166 /* Wrap */ 167 #define sXd(Name) CONCAT(s3d_, Name) 168 #define SXD(Name) CONCAT(S3D_, Name) 169 #define Xd(Name) CONCAT(Name, _3d) 170 171 #else 172 #error Unexpected dimensionility 173 #endif 174 175 static res_T 176 Xd(gebhart_radiative_path) 177 (struct sgf_device* dev, 178 struct accum* accums, 179 struct accum* accum_infinity, 180 struct accum* accum_medium, 181 struct ssp_rng* rng, 182 struct htable_bounce* path, 183 const double absorption_coef, /* m^-1 */ 184 const size_t ispectral_band, 185 const size_t primitive_id, 186 struct sgf_scene* scn) 187 { 188 struct sXd(scene_view)* view; 189 struct sXd(hit) hit; 190 struct sXd(primitive) prim; 191 struct htable_bounce_iterator it, end; 192 const double trans_min = 1.e-8; /* Minimum transmissivity threshold */ 193 double proba_reflec_spec; 194 double transmissivity, emissivity, reflectivity, specularity; 195 double medium_transmissivity; 196 #ifndef NDEBUG 197 double radiative_flux = 0.0; 198 #endif 199 double infinite_radiative_flux = 0.0; 200 double medium_radiative_flux = 0.0; 201 float vec0[3]; /* Temporary vector */ 202 float normal[3]; /* Geometric normal */ 203 float pos[3]; /* Radiative path position */ 204 float dir[3]; /* Radiative path direction */ 205 float range[2]; /* Traced ray range */ 206 res_T res = RES_OK; 207 ASSERT(accums && rng && scn); 208 ASSERT(absorption_coef < 0 || accum_medium); 209 ASSERT(absorption_coef >= 0 || accum_infinity); 210 211 #if SGF_DIMENSIONALITY == 2 212 view = scn->geometry.s2d.view; 213 #else 214 view = scn->geometry.s3d.view; 215 #endif 216 htable_bounce_clear(path); 217 218 /* Discard faces with no emissivity */ 219 emissivity = scene_get_emissivity(scn, primitive_id, ispectral_band); 220 if(emissivity <= 0.0) 221 return RES_OK; 222 223 /* Retrieve the sXd_scn primitive */ 224 sXd(scene_view_get_primitive)(view, (unsigned)primitive_id, &prim); 225 /* Get the geometric normal of the input primitive */ 226 Xd(primitive_get_normal)(&prim, normal); 227 /* Uniformly sample prim to define the origin of the radiative path */ 228 Xd(primitive_sample_position)(&prim, rng, pos); 229 /* Cosine weighted sampling of the radiative path direction around `normal' */ 230 ssp_ran_hemisphere_cos_float(rng, normal, dir, NULL); 231 232 transmissivity = 1.0; 233 for(;;) { /* Here we go */ 234 struct sXd(primitive) prim_from = prim; 235 range[0] = FLT_MIN, range[1] = FLT_MAX; 236 237 #if SGF_DIMENSIONALITY == 2 238 s2d_scene_view_trace_ray_3d(view, pos, dir, range, &prim_from, &hit); 239 #else 240 s3d_scene_view_trace_ray(view, pos, dir, range, &prim_from, &hit); 241 #endif 242 243 /* Handle medium absorption */ 244 if(absorption_coef >= 0) { 245 if(SXD(HIT_NONE)(&hit)) { /* The ray shoulnd't be outside the volume */ 246 log_error(dev, 247 "The radiative random walk goes to the infinity while the submitted geometry \n" 248 "should surround a close medium. This may be due to numerical issues or to an\n" 249 "invalid geometry definition.\n" 250 "Ray origin: %g, %g, %g\n", 251 SPLIT3(pos)); 252 return RES_BAD_OP; 253 } else { 254 double weight, ka; 255 256 /* Check the consistency of the medium description, i.e. the absorption 257 * coefficient must be the same when it is fetched from any primitive 258 * surrounding the current enclosure */ 259 ka = scene_get_absorption(scn, primitive_id, ispectral_band); 260 if(ka != absorption_coef) { 261 log_error(dev, "Inconsistent medium description.\n"); 262 return RES_BAD_ARG; 263 } 264 265 medium_transmissivity = exp(-absorption_coef*hit.distance); 266 weight = transmissivity * (1-medium_transmissivity); 267 transmissivity *= medium_transmissivity; 268 medium_radiative_flux += weight; 269 } 270 } else if(SXD(HIT_NONE)(&hit)) { /* The ray is outside the volume */ 271 infinite_radiative_flux = transmissivity; 272 break; 273 } 274 275 /* Retrieve the hit position and the hit primitive id */ 276 f3_mulf(vec0, dir, hit.distance); 277 f3_add(pos, vec0, pos); 278 prim = hit.prim; 279 280 /* Fetch material property */ 281 emissivity = scene_get_emissivity(scn, prim.scene_prim_id, ispectral_band); 282 specularity = scene_get_specularity(scn, prim.scene_prim_id, ispectral_band); 283 reflectivity = scene_get_reflectivity(scn, prim.scene_prim_id, ispectral_band); 284 if(transmissivity > trans_min) { 285 const double weight = transmissivity * emissivity; 286 double* pweight = htable_bounce_find(path, &prim.scene_prim_id); 287 if(pweight) { 288 *pweight += weight; 289 } else { 290 res = htable_bounce_set(path, &prim.scene_prim_id, &weight); 291 if(res != RES_OK) return res; 292 } 293 transmissivity = transmissivity * (1.0 - emissivity); 294 #ifndef NDEBUG 295 radiative_flux += weight; 296 #endif 297 } else { 298 /* Russian roulette */ 299 if(ssp_rng_canonical(rng) < emissivity) { 300 double* pweight = htable_bounce_find(path, &prim.scene_prim_id); 301 if(pweight) { 302 *pweight += transmissivity; 303 } else { 304 res = htable_bounce_set(path, &prim.scene_prim_id, &transmissivity); 305 if(res != RES_OK) return res; 306 } 307 #ifndef NDEBUG 308 radiative_flux += transmissivity; 309 #endif 310 break; 311 } 312 } 313 314 if(reflectivity <= 0.0) break; 315 316 proba_reflec_spec = specularity / reflectivity; 317 Xd(hit_get_normal)(&hit, normal); 318 319 if(ssp_rng_canonical(rng) >= proba_reflec_spec) { /* Diffuse reflection */ 320 ssp_ran_hemisphere_cos_float(rng, normal, dir, NULL); 321 ASSERT(f3_dot(normal, dir) > 0); 322 } else { /* Specular reflection */ 323 const float tmp = -2.f * f3_dot(dir, normal); 324 f3_mulf(vec0, normal, tmp); 325 f3_add(dir, dir, vec0); 326 f3_normalize(dir, dir); 327 } 328 } 329 330 /* Flush MC weights */ 331 htable_bounce_begin(path, &it); 332 htable_bounce_end(path, &end); 333 while(!htable_bounce_iterator_eq(&it, &end)) { 334 const size_t iprim = *htable_bounce_iterator_key_get(&it); 335 const double weight = *htable_bounce_iterator_data_get(&it); 336 accum_weight(accums + iprim, weight); 337 htable_bounce_iterator_next(&it); 338 } 339 if(medium_radiative_flux) { 340 accum_weight(accum_medium, medium_radiative_flux); 341 } 342 if(infinite_radiative_flux) { 343 accum_weight(accum_infinity, infinite_radiative_flux); 344 } 345 346 #if !defined(NDEBUG) 347 /* Check the energy conservation property */ 348 ASSERT(eq_eps 349 (radiative_flux + infinite_radiative_flux + medium_radiative_flux, 350 1.0, 1.e6)); 351 #endif 352 return RES_OK; 353 } 354 355 #undef sXd 356 #undef SXD 357 #undef Xd 358 #undef SGF_DIMENSIONALITY 359 360 #endif /* !SGF_DIMENSIONALITY */