htrdr_atmosphere.c (16797B)
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 #define _POSIX_C_SOURCE 200112L 25 26 #include "atmosphere/htrdr_atmosphere.h" 27 #include "atmosphere/htrdr_atmosphere_c.h" 28 #include "atmosphere/htrdr_atmosphere_args.h" 29 #include "atmosphere/htrdr_atmosphere_ground.h" 30 #include "atmosphere/htrdr_atmosphere_sun.h" 31 32 #include "core/htrdr_buffer.h" 33 #include "core/htrdr_log.h" 34 #include "core/htrdr_materials.h" 35 #include "core/htrdr_ran_wlen_cie_xyz.h" 36 #include "core/htrdr_ran_wlen_planck.h" 37 #include "core/htrdr_rectangle.h" 38 39 #include <high_tune/htsky.h> 40 41 #include <star/scam.h> 42 43 #include <rsys/cstr.h> 44 #include <rsys/double3.h> 45 46 #include <math.h> 47 48 /******************************************************************************* 49 * Helper functions 50 ******************************************************************************/ 51 /* Compute the number of fixed size bands used to discretized the spectral 52 * range */ 53 static size_t 54 compute_spectral_bands_count(const struct htrdr_atmosphere* cmd) 55 { 56 double wlen_range[2]; 57 double wlen_range_size; 58 size_t nbands; 59 ASSERT(cmd); 60 61 /* Compute size of the spectral range in nanometers */ 62 wlen_range[0] = cmd->wlen_range_m[0]*1.e9; 63 wlen_range[1] = cmd->wlen_range_m[1]*1.e9; 64 wlen_range_size = wlen_range[1] - wlen_range[0]; 65 66 /* Define as many intervals as wavelengths count in the spectral range */ 67 nbands = (size_t)rint(wlen_range_size); 68 69 return nbands; 70 } 71 72 static enum htsky_spectral_type 73 htrdr_to_sky_spectral_type(const enum htrdr_spectral_type type) 74 { 75 enum htsky_spectral_type spectype; 76 switch(type) { 77 case HTRDR_SPECTRAL_LW: 78 spectype = HTSKY_SPECTRAL_LW; 79 break; 80 case HTRDR_SPECTRAL_SW: 81 case HTRDR_SPECTRAL_SW_CIE_XYZ: 82 spectype = HTSKY_SPECTRAL_SW; 83 break; 84 default: FATAL("Unreachable code.\n"); break; 85 } 86 return spectype; 87 } 88 89 static INLINE void 90 spherical_to_cartesian_dir 91 (const double azimuth, /* In radians */ 92 const double elevation, /* In radians */ 93 double dir[3]) 94 { 95 double cos_azimuth; 96 double sin_azimuth; 97 double cos_elevation; 98 double sin_elevation; 99 ASSERT(azimuth >= 0 && azimuth < 2*PI); 100 ASSERT(elevation >= 0 && elevation <= PI/2.0); 101 ASSERT(dir); 102 103 cos_azimuth = cos(azimuth); 104 sin_azimuth = sin(azimuth); 105 cos_elevation = cos(elevation); 106 sin_elevation = sin(elevation); 107 108 dir[0] = cos_elevation * cos_azimuth; 109 dir[1] = cos_elevation * sin_azimuth; 110 dir[2] = sin_elevation; 111 } 112 113 static res_T 114 setup_camera_orthographic 115 (struct htrdr_atmosphere* cmd, 116 const struct htrdr_atmosphere_args* args) 117 { 118 struct scam_orthographic_args cam_args = SCAM_ORTHOGRAPHIC_ARGS_DEFAULT; 119 res_T res = RES_OK; 120 121 ASSERT(cmd && args); 122 ASSERT(cmd->output_type == HTRDR_ATMOSPHERE_ARGS_OUTPUT_IMAGE); 123 ASSERT(args->cam_type == HTRDR_ARGS_CAMERA_ORTHOGRAPHIC); 124 ASSERT(htrdr_args_camera_orthographic_check(&args->cam_ortho) == RES_OK); 125 ASSERT(htrdr_args_image_check(&args->image) == RES_OK); 126 127 d3_set(cam_args.position, args->cam_ortho.position); 128 d3_set(cam_args.target, args->cam_ortho.target); 129 d3_set(cam_args.up, args->cam_ortho.up); 130 cam_args.height = args->cam_ortho.height; 131 cam_args.aspect_ratio = 132 (double)args->image.definition[0] 133 / (double)args->image.definition[1]; 134 135 res = scam_create_orthographic 136 (htrdr_get_logger(cmd->htrdr), 137 htrdr_get_allocator(cmd->htrdr), 138 htrdr_get_verbosity_level(cmd->htrdr), 139 &cam_args, 140 &cmd->camera); 141 if(res != RES_OK) goto error; 142 143 exit: 144 return res; 145 error: 146 goto exit; 147 } 148 149 static res_T 150 setup_camera_perspective 151 (struct htrdr_atmosphere* cmd, 152 const struct htrdr_atmosphere_args* args) 153 { 154 struct scam_perspective_args cam_args = SCAM_PERSPECTIVE_ARGS_DEFAULT; 155 res_T res = RES_OK; 156 157 ASSERT(cmd && args); 158 ASSERT(cmd->output_type == HTRDR_ATMOSPHERE_ARGS_OUTPUT_IMAGE); 159 ASSERT(args->cam_type == HTRDR_ARGS_CAMERA_PERSPECTIVE); 160 ASSERT(htrdr_args_camera_perspective_check(&args->cam_persp) == RES_OK); 161 ASSERT(htrdr_args_image_check(&args->image) == RES_OK); 162 163 d3_set(cam_args.position, args->cam_persp.position); 164 d3_set(cam_args.target, args->cam_persp.target); 165 d3_set(cam_args.up, args->cam_persp.up); 166 cam_args.aspect_ratio = 167 (double)args->image.definition[0] 168 / (double)args->image.definition[1]; 169 cam_args.field_of_view = MDEG2RAD(args->cam_persp.fov_y); 170 cam_args.lens_radius = args->cam_persp.lens_radius; 171 cam_args.focal_distance = args->cam_persp.focal_dst; 172 173 res = scam_create_perspective 174 (htrdr_get_logger(cmd->htrdr), 175 htrdr_get_allocator(cmd->htrdr), 176 htrdr_get_verbosity_level(cmd->htrdr), 177 &cam_args, 178 &cmd->camera); 179 if(res != RES_OK) goto error; 180 181 exit: 182 return res; 183 error: 184 goto exit; 185 } 186 187 static res_T 188 setup_camera 189 (struct htrdr_atmosphere* cmd, 190 const struct htrdr_atmosphere_args* args) 191 { 192 res_T res = RES_OK; 193 ASSERT(cmd->output_type == HTRDR_ATMOSPHERE_ARGS_OUTPUT_IMAGE); 194 switch(args->cam_type) { 195 case HTRDR_ARGS_CAMERA_ORTHOGRAPHIC: 196 res = setup_camera_orthographic(cmd, args); 197 break; 198 case HTRDR_ARGS_CAMERA_PERSPECTIVE: 199 res = setup_camera_perspective(cmd, args); 200 break; 201 default: FATAL("Unreachable code.\n"); break; 202 } 203 return res; 204 } 205 206 static res_T 207 setup_flux_map 208 (struct htrdr_atmosphere* cmd, 209 const struct htrdr_atmosphere_args* args) 210 { 211 ASSERT(cmd && args); 212 ASSERT(cmd->output_type == HTRDR_ATMOSPHERE_ARGS_OUTPUT_FLUX_MAP); 213 214 if(args->spectral.spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ) { 215 htrdr_log_err(cmd->htrdr, 216 "the CIE 1931 XYZ spectral integration can be used only with a camera" 217 "sensor.\n"); 218 return RES_BAD_ARG; 219 } 220 221 return htrdr_rectangle_create 222 (cmd->htrdr, 223 args->flux_map.size, 224 args->flux_map.position, 225 args->flux_map.target, 226 args->flux_map.up, 227 &cmd->flux_map); 228 } 229 230 static res_T 231 setup_sensor 232 (struct htrdr_atmosphere* cmd, 233 const struct htrdr_atmosphere_args* args) 234 { 235 res_T res = RES_OK; 236 switch(cmd->output_type) { 237 case HTRDR_ATMOSPHERE_ARGS_OUTPUT_FLUX_MAP: 238 res = setup_flux_map(cmd, args); 239 break; 240 case HTRDR_ATMOSPHERE_ARGS_OUTPUT_IMAGE: 241 res = setup_camera(cmd, args); 242 break; 243 default: /* Nothing to do */ break; 244 } 245 return res; 246 } 247 248 static res_T 249 dump_volumetric_acceleration_structure(struct htrdr_atmosphere* cmd) 250 { 251 size_t nbands; 252 size_t i; 253 res_T res = RES_OK; 254 ASSERT(cmd); 255 256 nbands = htsky_get_spectral_bands_count(cmd->sky); 257 258 /* Nothing to do */ 259 if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit; 260 261 FOR_EACH(i, 0, nbands) { 262 size_t iquad; 263 const size_t iband = htsky_get_spectral_band_id(cmd->sky, i); 264 const size_t nquads = htsky_get_spectral_band_quadrature_length 265 (cmd->sky, iband); 266 267 FOR_EACH(iquad, 0, nquads) { 268 res = htsky_dump_cloud_vtk(cmd->sky, iband, iquad, cmd->output); 269 if(res != RES_OK) goto error; 270 fprintf(cmd->output, "---\n"); 271 } 272 } 273 274 exit: 275 return res; 276 error: 277 goto exit; 278 } 279 280 static void 281 atmosphere_release(ref_T* ref) 282 { 283 struct htrdr_atmosphere* cmd = CONTAINER_OF(ref, struct htrdr_atmosphere, ref); 284 struct htrdr* htrdr = NULL; 285 ASSERT(ref); 286 287 if(cmd->ground) htrdr_atmosphere_ground_ref_put(cmd->ground); 288 if(cmd->mats) htrdr_materials_ref_put(cmd->mats); 289 if(cmd->sun) htrdr_atmosphere_sun_ref_put(cmd->sun); 290 if(cmd->cie) htrdr_ran_wlen_cie_xyz_ref_put(cmd->cie); 291 if(cmd->planck) htrdr_ran_wlen_planck_ref_put(cmd->planck); 292 if(cmd->camera) SCAM(ref_put(cmd->camera)); 293 if(cmd->flux_map) htrdr_rectangle_ref_put(cmd->flux_map); 294 if(cmd->buf) htrdr_buffer_ref_put(cmd->buf); 295 if(cmd->sky) HTSKY(ref_put(cmd->sky)); 296 if(cmd->output && cmd->output != stdout) fclose(cmd->output); 297 str_release(&cmd->output_name); 298 299 htrdr = cmd->htrdr; 300 MEM_RM(htrdr_get_allocator(htrdr), cmd); 301 htrdr_ref_put(htrdr); 302 } 303 304 /******************************************************************************* 305 * Exported functions 306 ******************************************************************************/ 307 res_T 308 htrdr_atmosphere_create 309 (struct htrdr* htrdr, 310 const struct htrdr_atmosphere_args* args, 311 struct htrdr_atmosphere** out_cmd) 312 { 313 struct htrdr_atmosphere* cmd = NULL; 314 struct htsky_args htsky_args = HTSKY_ARGS_DEFAULT; 315 double sun_dir[3]; 316 double spectral_range[2]; 317 const char* output_name = NULL; 318 size_t nintervals; /* #bands used to discretized the spectral curve */ 319 res_T res = RES_OK; 320 ASSERT(htrdr && args && out_cmd); 321 322 cmd = MEM_CALLOC(htrdr_get_allocator(htrdr), 1, sizeof(*cmd)); 323 if(!cmd) { 324 htrdr_log_err(htrdr, 325 "%s: could not allocate the htrdr_atmosphere data.\n", FUNC_NAME); 326 res = RES_MEM_ERR; 327 goto error; 328 } 329 ref_init(&cmd->ref); 330 str_init(htrdr_get_allocator(htrdr), &cmd->output_name); 331 cmd->output_type = args->output_type; 332 cmd->verbose = args->verbose; 333 cmd->spp = args->image.spp; 334 cmd->width = args->image.definition[0]; 335 cmd->height = args->image.definition[1]; 336 cmd->grid_max_definition[0] = args->grid_max_definition[0]; 337 cmd->grid_max_definition[1] = args->grid_max_definition[1]; 338 cmd->grid_max_definition[2] = args->grid_max_definition[2]; 339 cmd->spectral_type = args->spectral.spectral_type; 340 cmd->ref_temperature = args->spectral.ref_temperature; 341 cmd->sky_mtl_name = args->sky_mtl_name; 342 343 /* Get ownership on the htrdr structure */ 344 htrdr_ref_get(htrdr); 345 cmd->htrdr = htrdr; 346 347 if(!args->filename_output) { 348 cmd->output = stdout; 349 output_name = "<stdout>"; 350 } else if(htrdr_get_mpi_rank(htrdr) != 0) { 351 cmd->output = NULL; 352 output_name = "<null>"; 353 } else { 354 res = htrdr_open_output_stream(htrdr, args->filename_output, 0/*read*/, 355 args->force_overwriting, &cmd->output); 356 if(res != RES_OK) goto error; 357 output_name = args->filename_output; 358 } 359 res = str_set(&cmd->output_name, output_name); 360 if(res != RES_OK) { 361 htrdr_log_err(htrdr, 362 "%s: could not store the name of the output stream `%s' -- %s.\n", 363 FUNC_NAME, output_name, res_to_cstr(res)); 364 goto error; 365 } 366 367 /* Materials are necessary only if a ground geometry is defined */ 368 if(args->filename_obj) { 369 res = htrdr_materials_create(htrdr, args->filename_mtl, &cmd->mats); 370 if(res != RES_OK) goto error; 371 } 372 373 res = htrdr_atmosphere_ground_create(htrdr, args->filename_obj, cmd->mats, 374 args->repeat_ground, &cmd->ground); 375 if(res != RES_OK) goto error; 376 377 res = setup_sensor(cmd, args); 378 if(res != RES_OK) goto error; 379 380 res = htrdr_atmosphere_sun_create(cmd->htrdr, &cmd->sun); 381 if(res != RES_OK) goto error; 382 spherical_to_cartesian_dir 383 (MDEG2RAD(args->sun_azimuth), MDEG2RAD(args->sun_elevation), sun_dir); 384 htrdr_atmosphere_sun_set_direction(cmd->sun, sun_dir); 385 386 htsky_args.htcp_filename = args->filename_les; 387 htsky_args.htgop_filename = args->filename_gas; 388 htsky_args.htmie_filename = args->filename_mie; 389 htsky_args.cache_filename = args->filename_cache; 390 htsky_args.grid_max_definition[0] = args->grid_max_definition[0]; 391 htsky_args.grid_max_definition[1] = args->grid_max_definition[1]; 392 htsky_args.grid_max_definition[2] = args->grid_max_definition[2]; 393 htsky_args.optical_thickness = args->optical_thickness; 394 htsky_args.nthreads = (unsigned)htrdr_get_threads_count(htrdr); 395 htsky_args.repeat_clouds = args->repeat_clouds; 396 htsky_args.verbose = htrdr_get_mpi_rank(htrdr) == 0 ? args->verbose : 0; 397 htsky_args.spectral_type = htrdr_to_sky_spectral_type(args->spectral.spectral_type); 398 htsky_args.wlen_range[0] = args->spectral.wlen_range[0]; 399 htsky_args.wlen_range[1] = args->spectral.wlen_range[1]; 400 res = htsky_create(htrdr_get_logger(htrdr), htrdr_get_allocator(htrdr), 401 &htsky_args, &cmd->sky); 402 if(res != RES_OK) goto error; 403 404 HTSKY(get_raw_spectral_bounds(cmd->sky, spectral_range)); 405 406 spectral_range[0] = MMAX(args->spectral.wlen_range[0], spectral_range[0]); 407 spectral_range[1] = MMIN(args->spectral.wlen_range[1], spectral_range[1]); 408 if(spectral_range[0] != args->spectral.wlen_range[0] 409 || spectral_range[1] != args->spectral.wlen_range[1]) { 410 htrdr_log_warn(htrdr, 411 "%s: the submitted spectral range overflowed the spectral data.\n", 412 FUNC_NAME); 413 } 414 415 cmd->wlen_range_m[0] = spectral_range[0]*1e-9; /* Convert in meters */ 416 cmd->wlen_range_m[1] = spectral_range[1]*1e-9; /* Convert in meters */ 417 418 /* Compute the number of fixed sized bands used to accelerate the sampling of 419 * spectral data */ 420 nintervals = compute_spectral_bands_count(cmd); 421 422 if(cmd->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ) { 423 res = htrdr_ran_wlen_cie_xyz_create 424 (htrdr, spectral_range, nintervals, &cmd->cie); 425 if(res != RES_OK) goto error; 426 } else { 427 if(cmd->ref_temperature <= 0) { 428 htrdr_log_err(htrdr, "%s: invalid reference temperature %g K.\n", 429 FUNC_NAME, cmd->ref_temperature); 430 res = RES_BAD_ARG; 431 goto error; 432 } 433 res = htrdr_ran_wlen_planck_create 434 (htrdr, spectral_range, nintervals, cmd->ref_temperature, &cmd->planck); 435 if(res != RES_OK) goto error; 436 } 437 438 if(cmd->output_type != HTRDR_ATMOSPHERE_ARGS_OUTPUT_OCTREES) { 439 struct htrdr_pixel_format pixfmt = HTRDR_PIXEL_FORMAT_NULL; 440 atmosphere_get_pixel_format(cmd, &pixfmt); 441 442 /* Setup the buffer layout */ 443 cmd->buf_layout.width = args->image.definition[0]; 444 cmd->buf_layout.height = args->image.definition[1]; 445 cmd->buf_layout.pitch = args->image.definition[0] * pixfmt.size; 446 cmd->buf_layout.elmt_size = pixfmt.size; 447 cmd->buf_layout.alignment = pixfmt.alignment; 448 449 /* Create the image buffer only on the master process; the image parts 450 * rendered by the others processes are gathered onto the master process */ 451 if(htrdr_get_mpi_rank(htrdr) == 0) { 452 res = htrdr_buffer_create(htrdr, &cmd->buf_layout, &cmd->buf); 453 if(res != RES_OK) goto error; 454 } 455 } 456 457 exit: 458 *out_cmd = cmd; 459 return res; 460 error: 461 if(cmd) { 462 htrdr_atmosphere_ref_put(cmd); 463 cmd = NULL; 464 } 465 goto exit; 466 } 467 468 void 469 htrdr_atmosphere_ref_get(struct htrdr_atmosphere* cmd) 470 { 471 ASSERT(cmd); 472 ref_get(&cmd->ref); 473 } 474 475 void 476 htrdr_atmosphere_ref_put(struct htrdr_atmosphere* cmd) 477 { 478 ASSERT(cmd); 479 ref_put(&cmd->ref, atmosphere_release); 480 } 481 482 res_T 483 htrdr_atmosphere_run(struct htrdr_atmosphere* cmd) 484 { 485 res_T res = RES_OK; 486 switch(cmd->output_type) { 487 case HTRDR_ATMOSPHERE_ARGS_OUTPUT_IMAGE: 488 case HTRDR_ATMOSPHERE_ARGS_OUTPUT_FLUX_MAP: 489 res = atmosphere_draw_map(cmd); 490 break; 491 case HTRDR_ATMOSPHERE_ARGS_OUTPUT_OCTREES: 492 res = dump_volumetric_acceleration_structure(cmd); 493 break; 494 default: FATAL("Unreachable code.\n"); break; 495 } 496 return res; 497 } 498 499 /******************************************************************************* 500 * Local functions 501 ******************************************************************************/ 502 void 503 atmosphere_get_pixel_format 504 (const struct htrdr_atmosphere* cmd, 505 struct htrdr_pixel_format* fmt) 506 { 507 ASSERT(cmd && fmt); 508 switch(cmd->output_type) { 509 case HTRDR_ATMOSPHERE_ARGS_OUTPUT_FLUX_MAP: 510 fmt->size = sizeof(struct atmosphere_pixel_flux); 511 fmt->alignment = ALIGNOF(struct atmosphere_pixel_flux); 512 break; 513 case HTRDR_ATMOSPHERE_ARGS_OUTPUT_IMAGE: 514 switch(cmd->spectral_type) { 515 case HTRDR_SPECTRAL_LW: 516 case HTRDR_SPECTRAL_SW: 517 fmt->size = sizeof(struct atmosphere_pixel_xwave); 518 fmt->alignment = ALIGNOF(struct atmosphere_pixel_xwave); 519 break; 520 case HTRDR_SPECTRAL_SW_CIE_XYZ: 521 fmt->size = sizeof(struct atmosphere_pixel_image); 522 fmt->alignment = ALIGNOF(struct atmosphere_pixel_image); 523 break; 524 default: FATAL("Unreachable code.\n"); break; 525 } 526 break; 527 default: FATAL("Unreachable code.\n"); break; 528 } 529 } 530