htrdr_atmosphere_compute_radiance_sw.c (17508B)
1 /* Copyright (C) 2018-2019, 2022-2025 Centre National de la Recherche Scientifique 2 * Copyright (C) 2020-2022 Institut Mines Télécom Albi-Carmaux 3 * Copyright (C) 2022-2025 Institut Pierre-Simon Laplace 4 * Copyright (C) 2022-2025 Institut de Physique du Globe de Paris 5 * Copyright (C) 2018-2025 |Méso|Star> (contact@meso-star.com) 6 * Copyright (C) 2022-2025 Observatoire de Paris 7 * Copyright (C) 2022-2025 Université de Reims Champagne-Ardenne 8 * Copyright (C) 2022-2025 Université de Versaille Saint-Quentin 9 * Copyright (C) 2018-2019, 2022-2025 Université Paul Sabatier 10 * 11 * This program is free software: you can redistribute it and/or modify 12 * it under the terms of the GNU General Public License as published by 13 * the Free Software Foundation, either version 3 of the License, or 14 * (at your option) any later version. 15 * 16 * This program is distributed in the hope that it will be useful, 17 * but WITHOUT ANY WARRANTY; without even the implied warranty of 18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19 * GNU General Public License for more details. 20 * 21 * You should have received a copy of the GNU General Public License 22 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 23 24 #include "atmosphere/htrdr_atmosphere_c.h" 25 #include "atmosphere/htrdr_atmosphere_ground.h" 26 #include "atmosphere/htrdr_atmosphere_sun.h" 27 28 #include "core/htrdr_interface.h" 29 30 #include <high_tune/htsky.h> 31 32 #include <star/s3d.h> 33 #include <star/ssf.h> 34 #include <star/ssp.h> 35 #include <star/svx.h> 36 37 #include <rsys/double2.h> 38 #include <rsys/float2.h> 39 #include <rsys/float3.h> 40 41 struct scattering_context { 42 struct ssp_rng* rng; 43 const struct htsky* sky; 44 size_t iband; /* Index of the spectral band */ 45 size_t iquad; /* Index of the quadrature point into the band */ 46 47 double Ts; /* Sampled optical thickness */ 48 double traversal_dst; /* Distance traversed along the ray */ 49 }; 50 static const struct scattering_context SCATTERING_CONTEXT_NULL = { 51 NULL, NULL, 0, 0, 0, 0 52 }; 53 54 struct transmissivity_context { 55 struct ssp_rng* rng; 56 const struct htsky* sky; 57 size_t iband; /* Index of the spectral */ 58 size_t iquad; /* Index of the quadrature point into the band */ 59 60 double Ts; /* Sampled optical thickness */ 61 double Tmin; /* Minimal optical thickness */ 62 double traversal_dst; /* Distance traversed along the ray */ 63 64 enum htsky_property prop; 65 }; 66 static const struct transmissivity_context TRANSMISSION_CONTEXT_NULL = { 67 NULL, NULL, 0, 0, 0, 0, 0, 0 68 }; 69 70 /******************************************************************************* 71 * Helper functions 72 ******************************************************************************/ 73 static int 74 scattering_hit_filter 75 (const struct svx_hit* hit, 76 const double org[3], 77 const double dir[3], 78 const double range[2], 79 void* context) 80 { 81 struct scattering_context* ctx = context; 82 double ks_max; 83 int pursue_traversal = 1; 84 ASSERT(hit && ctx && !SVX_HIT_NONE(hit) && org && dir && range); 85 (void)range; 86 87 ks_max = htsky_fetch_svx_voxel_property(ctx->sky, HTSKY_Ks, 88 HTSKY_SVX_MAX, HTSKY_CPNT_MASK_ALL, ctx->iband, ctx->iquad, &hit->voxel); 89 90 ctx->traversal_dst = hit->distance[0]; 91 92 /* Iterate until a collision occurs into the voxel or until the ray 93 * does not collide the voxel */ 94 for(;;) { 95 /* Compute tau for the current leaf */ 96 const double vox_dst = hit->distance[1] - ctx->traversal_dst; 97 const double T = vox_dst * ks_max; 98 99 /* A collision occurs behind `vox_dst' */ 100 if(ctx->Ts > T) { 101 ctx->Ts -= T; 102 ctx->traversal_dst = hit->distance[1]; 103 pursue_traversal = 1; 104 break; 105 106 /* A real/null collision occurs before `vox_dst' */ 107 } else { 108 double pos[3]; 109 double proba; 110 double ks; 111 const double collision_dst = ctx->Ts / ks_max; 112 113 /* Compute the traversed distance up to the challenged collision */ 114 ctx->traversal_dst += collision_dst; 115 ASSERT(ctx->traversal_dst >= hit->distance[0]); 116 ASSERT(ctx->traversal_dst <= hit->distance[1]); 117 118 /* Stop the ray whenever the traversal distance without any scattering 119 * event is too high. It means the maximum scattering coefficient has a 120 * very small value, and the returned radiance is null. This can only 121 * happen when the voxel has a [quasi] infinite length in the propagation 122 * direction. */ 123 if(ctx->traversal_dst > 1e9) break; 124 125 /* Compute the world space position where a collision may occur */ 126 pos[0] = org[0] + ctx->traversal_dst * dir[0]; 127 pos[1] = org[1] + ctx->traversal_dst * dir[1]; 128 pos[2] = org[2] + ctx->traversal_dst * dir[2]; 129 130 ks = htsky_fetch_raw_property(ctx->sky, HTSKY_Ks, 131 HTSKY_CPNT_MASK_ALL, ctx->iband, ctx->iquad, pos, -DBL_MAX, DBL_MAX); 132 133 /* Handle the case that ks_max is not *really* the max */ 134 proba = ks / ks_max; 135 136 if(ssp_rng_canonical(ctx->rng) < proba) {/* Collide <=> real scattering */ 137 pursue_traversal = 0; 138 break; 139 } else { /* Null collision */ 140 ctx->Ts = ssp_ran_exp(ctx->rng, 1); /* Sample a new optical thickness */ 141 } 142 } 143 } 144 return pursue_traversal; 145 } 146 147 static int 148 transmissivity_hit_filter 149 (const struct svx_hit* hit, 150 const double org[3], 151 const double dir[3], 152 const double range[2], 153 void* context) 154 { 155 struct transmissivity_context* ctx = context; 156 int comp_mask = HTSKY_CPNT_MASK_ALL; 157 double k_max; 158 double k_min; 159 int pursue_traversal = 1; 160 ASSERT(hit && ctx && !SVX_HIT_NONE(hit) && org && dir && range); 161 (void)range; 162 163 k_min = htsky_fetch_svx_voxel_property(ctx->sky, ctx->prop, 164 HTSKY_SVX_MIN, comp_mask, ctx->iband, ctx->iquad, &hit->voxel); 165 k_max = htsky_fetch_svx_voxel_property(ctx->sky, ctx->prop, 166 HTSKY_SVX_MAX, comp_mask, ctx->iband, ctx->iquad, &hit->voxel); 167 ASSERT(k_min <= k_max); 168 169 ctx->Tmin += (hit->distance[1] - hit->distance[0]) * k_min; 170 ctx->traversal_dst = hit->distance[0]; 171 172 /* Iterate until a collision occurs into the voxel or until the ray 173 * does not collide the voxel */ 174 for(;;) { 175 const double vox_dst = hit->distance[1] - ctx->traversal_dst; 176 const double Tdif = vox_dst * (k_max-k_min); 177 178 /* A collision occurs behind `vox_dst' */ 179 if(ctx->Ts > Tdif) { 180 ctx->Ts -= Tdif; 181 ctx->traversal_dst = hit->distance[1]; 182 pursue_traversal = 1; 183 break; 184 185 /* A real/null collision occurs before `vox_dst' */ 186 } else { 187 double x[3]; 188 double k; 189 double proba; 190 double collision_dst = ctx->Ts / (k_max - k_min); 191 192 /* Compute the traversed distance up to the challenged collision */ 193 ctx->traversal_dst += collision_dst; 194 ASSERT(ctx->traversal_dst >= hit->distance[0]); 195 ASSERT(ctx->traversal_dst <= hit->distance[1]); 196 197 /* Compute the world space position where a collision may occur */ 198 x[0] = org[0] + ctx->traversal_dst * dir[0]; 199 x[1] = org[1] + ctx->traversal_dst * dir[1]; 200 x[2] = org[2] + ctx->traversal_dst * dir[2]; 201 202 k = htsky_fetch_raw_property(ctx->sky, ctx->prop, 203 comp_mask, ctx->iband, ctx->iquad, x, k_min, k_max); 204 ASSERT(k >= k_min && k <= k_max); 205 206 proba = (k - k_min) / (k_max - k_min); 207 208 if(ssp_rng_canonical(ctx->rng) < proba) { /* Collide */ 209 pursue_traversal = 0; 210 break; 211 } else { /* Null collision */ 212 ctx->Ts = ssp_ran_exp(ctx->rng, 1); /* Sample a new optical thickness */ 213 } 214 } 215 } 216 return pursue_traversal; 217 } 218 219 static double 220 transmissivity 221 (struct htrdr_atmosphere* cmd, 222 struct ssp_rng* rng, 223 const enum htsky_property prop, 224 const size_t iband, 225 const size_t iquad, 226 const double pos[3], 227 const double dir[3], 228 const double range[2]) 229 { 230 struct svx_hit svx_hit; 231 struct transmissivity_context transmissivity_ctx = TRANSMISSION_CONTEXT_NULL; 232 233 ASSERT(cmd && rng && pos && dir && range); 234 235 transmissivity_ctx.rng = rng; 236 transmissivity_ctx.sky = cmd->sky; 237 transmissivity_ctx.iband = iband; 238 transmissivity_ctx.iquad = iquad; 239 transmissivity_ctx.Ts = ssp_ran_exp(rng, 1); /* Sample an optical thickness */ 240 transmissivity_ctx.prop = prop; 241 242 /* Compute the transmissivity */ 243 HTSKY(trace_ray(cmd->sky, pos, dir, range, NULL, 244 transmissivity_hit_filter, &transmissivity_ctx, iband, iquad, &svx_hit)); 245 246 if(SVX_HIT_NONE(&svx_hit)) { 247 return transmissivity_ctx.Tmin ? exp(-transmissivity_ctx.Tmin) : 1.0; 248 } else { 249 return 0; 250 } 251 } 252 253 /******************************************************************************* 254 * Local functions 255 ******************************************************************************/ 256 double 257 atmosphere_compute_radiance_sw 258 (struct htrdr_atmosphere* cmd, 259 const size_t ithread, 260 struct ssp_rng* rng, 261 const int cpnt_mask, /* Combination of enum htrdr_radiance_cpnt_flag */ 262 const double pos_in[3], 263 const double dir_in[3], 264 const double wlen, /* In nanometer */ 265 const size_t iband, 266 const size_t iquad) 267 { 268 struct s3d_hit s3d_hit = S3D_HIT_NULL; 269 struct s3d_hit s3d_hit_tmp = S3D_HIT_NULL; 270 struct s3d_hit s3d_hit_prev = S3D_HIT_NULL; 271 struct svx_hit svx_hit = SVX_HIT_NULL; 272 struct ssf_phase* phase_hg = NULL; 273 struct ssf_phase* phase_rayleigh = NULL; 274 275 double pos[3]; 276 double dir[3]; 277 double range[2]; 278 double pos_next[3]; 279 double dir_next[3]; 280 double band_bounds[2]; /* In nanometers */ 281 282 double R; 283 double r; /* Random number */ 284 double wo[3]; /* -dir */ 285 double pdf; 286 double Tr; /* Overall transmissivity */ 287 double Tr_abs; /* Absorption transmissivity */ 288 double L_sun; /* Sun radiance in W.m^-2.sr^-1 */ 289 double sun_dir[3]; 290 double ksi = 1; /* Throughput */ 291 double w = 0; /* MC weight */ 292 double g = 0; /* Asymmetry parameter of the HG phase function */ 293 294 ASSERT(cmd && rng && pos_in && dir_in); 295 ASSERT(cmd->spectral_type == HTRDR_SPECTRAL_SW 296 || cmd->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ); 297 298 /* Create the Henyey-Greenstein phase function */ 299 CHK(RES_OK == ssf_phase_create 300 (htrdr_get_thread_allocator(cmd->htrdr, ithread), 301 &ssf_phase_hg, 302 &phase_hg)); 303 304 /* Create the Rayleigh phase function */ 305 CHK(RES_OK == ssf_phase_create 306 (htrdr_get_thread_allocator(cmd->htrdr, ithread), 307 &ssf_phase_rayleigh, 308 &phase_rayleigh)); 309 310 /* Setup the phase function for this wavelength */ 311 g = htsky_fetch_per_wavelength_particle_phase_function_asymmetry_parameter 312 (cmd->sky, wlen); 313 SSF(phase_hg_setup(phase_hg, g)); 314 315 /* Fetch sun properties. Note that the sun spectral data are defined by bands 316 * that, actually are the same of the SW spectral bands defined in the 317 * default "ecrad_opt_prot.txt" file provided by the HTGOP project. */ 318 htsky_get_spectral_band_bounds(cmd->sky, iband, band_bounds); 319 ASSERT(band_bounds[0] <= wlen && wlen <= band_bounds[1]); 320 L_sun = htrdr_atmosphere_sun_get_radiance(cmd->sun, wlen); 321 d3_set(pos, pos_in); 322 d3_set(dir, dir_in); 323 324 if((cpnt_mask & ATMOSPHERE_RADIANCE_DIRECT) /* Handle direct contribution */ 325 && htrdr_atmosphere_sun_is_dir_in_solar_cone(cmd->sun, dir)) { 326 /* Check that the ray is not occluded along the submitted range */ 327 d2(range, 0, FLT_MAX); 328 HTRDR(atmosphere_ground_trace_ray 329 (cmd->ground, pos, dir, range, NULL, &s3d_hit_tmp)); 330 if(!S3D_HIT_NONE(&s3d_hit_tmp)) { 331 Tr = 0; 332 } else { 333 Tr = transmissivity 334 (cmd, rng, HTSKY_Kext, iband, iquad , pos, dir, range); 335 w = L_sun * Tr; 336 } 337 } 338 339 if((cpnt_mask & ATMOSPHERE_RADIANCE_DIFFUSE) == 0) 340 goto exit; /* Discard diffuse contribution */ 341 342 /* Radiative random walk */ 343 for(;;) { 344 struct scattering_context scattering_ctx = SCATTERING_CONTEXT_NULL; 345 struct ssf_bsdf* bsdf = NULL; 346 struct ssf_phase* phase = NULL; 347 double N[3]; 348 double bounce_reflectivity = 1; 349 double sun_dir_pdf; 350 int surface_scattering = 0; /* Define if hit a surface */ 351 int bsdf_type = 0; 352 353 /* Find the first intersection with a surface */ 354 d2(range, 0, DBL_MAX); 355 HTRDR(atmosphere_ground_trace_ray 356 (cmd->ground, pos, dir, range, &s3d_hit_prev, &s3d_hit)); 357 358 /* Sample an optical thickness */ 359 scattering_ctx.Ts = ssp_ran_exp(rng, 1); 360 361 /* Setup the remaining scattering context fields */ 362 scattering_ctx.rng = rng; 363 scattering_ctx.sky = cmd->sky; 364 scattering_ctx.iband = iband; 365 scattering_ctx.iquad = iquad; 366 367 /* Define if a scattering event occurs */ 368 d2(range, 0, s3d_hit.distance); 369 HTSKY(trace_ray(cmd->sky, pos, dir, range, NULL, 370 scattering_hit_filter, &scattering_ctx, iband, iquad, &svx_hit)); 371 372 /* No scattering and no surface reflection. Stop the radiative random walk */ 373 if(S3D_HIT_NONE(&s3d_hit) && SVX_HIT_NONE(&svx_hit)) { 374 break; 375 } 376 ASSERT(SVX_HIT_NONE(&svx_hit) 377 || ( svx_hit.distance[0] <= scattering_ctx.traversal_dst 378 && svx_hit.distance[1] >= scattering_ctx.traversal_dst)); 379 380 /* Negate the incoming dir to match the convention of the SSF library */ 381 d3_minus(wo, dir); 382 383 /* Define if the scattering occurs at a surface */ 384 surface_scattering = SVX_HIT_NONE(&svx_hit); 385 386 /* Compute the new position */ 387 pos_next[0] = pos[0] + dir[0]*scattering_ctx.traversal_dst; 388 pos_next[1] = pos[1] + dir[1]*scattering_ctx.traversal_dst; 389 pos_next[2] = pos[2] + dir[2]*scattering_ctx.traversal_dst; 390 391 /* Define the previous hit surface used to avoid self hit */ 392 s3d_hit_prev = surface_scattering ? s3d_hit : S3D_HIT_NULL; 393 394 /* Define the absorption transmissivity from the current position to the 395 * next position */ 396 d2(range, 0, scattering_ctx.traversal_dst); 397 Tr_abs = transmissivity 398 (cmd, rng, HTSKY_Ka, iband, iquad, pos, dir, range); 399 if(Tr_abs <= 0) break; 400 401 /* Sample the scattering direction */ 402 if(surface_scattering) { /* Scattering at a surface */ 403 struct htrdr_interface interf = HTRDR_INTERFACE_NULL; 404 const struct htrdr_mtl* mtl = NULL; 405 406 /* Fetch the hit interface material and build its BSDF */ 407 htrdr_atmosphere_ground_get_interface(cmd->ground, &s3d_hit, &interf); 408 mtl = htrdr_interface_fetch_hit_mtl(&interf, dir, &s3d_hit); 409 HTRDR(mtl_create_bsdf(cmd->htrdr, mtl, ithread, wlen, rng, &bsdf)); 410 411 /* Revert the normal if necessary to match the SSF convention */ 412 d3_normalize(N, d3_set_f3(N, s3d_hit.normal)); 413 if(d3_dot(N, wo) < 0) d3_minus(N, N); 414 415 /* Sample scattering direction */ 416 bounce_reflectivity = ssf_bsdf_sample 417 (bsdf, rng, wo, N, dir_next, &bsdf_type, &pdf); 418 if(!(bsdf_type & SSF_REFLECTION)) { /* Handle only reflections */ 419 bounce_reflectivity = 0; 420 } 421 422 } else { /* Scattering in a volume */ 423 double ks_particle; /* Scattering coefficient of the particles */ 424 double ks_gas; /* Scattering coefficient of the gaz */ 425 double ks; /* Overall scattering coefficient */ 426 427 ks_gas = htsky_fetch_raw_property(cmd->sky, HTSKY_Ks, 428 HTSKY_CPNT_FLAG_GAS, iband, iquad, pos_next, -DBL_MAX, DBL_MAX); 429 ks_particle = htsky_fetch_raw_property(cmd->sky, HTSKY_Ks, 430 HTSKY_CPNT_FLAG_PARTICLES, iband, iquad, pos_next, -DBL_MAX, DBL_MAX); 431 ks = ks_particle + ks_gas; 432 433 r = ssp_rng_canonical(rng); 434 if(r < ks_gas / ks) { /* Gas scattering */ 435 phase = phase_rayleigh; 436 } else { /* Cloud scattering */ 437 phase = phase_hg; 438 } 439 440 /* Sample scattering direction */ 441 ssf_phase_sample(phase, rng, wo, dir_next, NULL); 442 ssf_phase_ref_get(phase); 443 } 444 445 /* Sample the direction of the direct contribution */ 446 if(surface_scattering && (bsdf_type & SSF_SPECULAR)) { 447 if(!htrdr_atmosphere_sun_is_dir_in_solar_cone(cmd->sun, dir_next)) { 448 R = 0; /* No direct lightning */ 449 } else { 450 sun_dir[0] = dir_next[0]; 451 sun_dir[1] = dir_next[1]; 452 sun_dir[2] = dir_next[2]; 453 R = d3_dot(N, sun_dir)<0/* Below the ground*/ ? 0 : bounce_reflectivity; 454 } 455 sun_dir_pdf = 1.0; 456 } else { 457 /* Sample a sun direction */ 458 sun_dir_pdf = htrdr_atmosphere_sun_sample_direction 459 (cmd->sun, rng, sun_dir); 460 if(surface_scattering) { 461 R = d3_dot(N, sun_dir) < 0/* Below the ground */ 462 ? 0 : ssf_bsdf_eval(bsdf, wo, N, sun_dir) * d3_dot(N, sun_dir); 463 } else { 464 R = ssf_phase_eval(phase, wo, sun_dir); 465 } 466 } 467 468 /* The direct contribution to the scattering point is not null so we need 469 * to compute the transmissivity from sun to scatt pt */ 470 if(R <= 0) { 471 Tr = 0; 472 } else { 473 /* Check that the sun is visible from the new position */ 474 d2(range, 0, FLT_MAX); 475 HTRDR(atmosphere_ground_trace_ray 476 (cmd->ground, pos_next, sun_dir, range, &s3d_hit_prev, &s3d_hit_tmp)); 477 478 /* Compute the sun transmissivity */ 479 if(!S3D_HIT_NONE(&s3d_hit_tmp)) { 480 Tr = 0; 481 } else { 482 Tr = transmissivity 483 (cmd, rng, HTSKY_Kext, iband, iquad, pos_next, sun_dir, range); 484 } 485 } 486 487 /* Release the scattering function */ 488 if(surface_scattering) { 489 SSF(bsdf_ref_put(bsdf)); 490 } else { 491 SSF(phase_ref_put(phase)); 492 } 493 494 /* Update the MC weight */ 495 ksi *= Tr_abs; 496 w += ksi * L_sun * Tr * R / sun_dir_pdf; 497 498 /* Russian roulette wrt surface scattering */ 499 if(surface_scattering && ssp_rng_canonical(rng) >= bounce_reflectivity) 500 break; 501 502 /* Setup the next random walk state */ 503 d3_set(pos, pos_next); 504 d3_set(dir, dir_next); 505 } 506 507 exit: 508 SSF(phase_ref_put(phase_hg)); 509 SSF(phase_ref_put(phase_rayleigh)); 510 return w; 511 } 512