htrdr_planets_draw_map.c (15068B)
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 "planets/htrdr_planets_c.h" 25 #include "planets/htrdr_planets_source.h" 26 27 #include "core/htrdr.h" 28 #include "core/htrdr_accum.h" 29 #include "core/htrdr_buffer.h" 30 #include "core/htrdr_draw_map.h" 31 #include "core/htrdr_log.h" 32 #include "core/htrdr_ran_wlen_cie_xyz.h" 33 #include "core/htrdr_ran_wlen_discrete.h" 34 #include "core/htrdr_ran_wlen_planck.h" 35 36 #include <rad-net/rnatm.h> 37 #include <star/scam.h> 38 #include <star/ssp.h> 39 40 #include <rsys/clock_time.h> 41 42 /******************************************************************************* 43 * Helper functions 44 ******************************************************************************/ 45 static void 46 draw_pixel_xwave 47 (struct htrdr* htrdr, 48 const struct htrdr_draw_pixel_args* args, 49 void* data) 50 { 51 struct planets_compute_radiance_args rad_args = 52 PLANETS_COMPUTE_RADIANCE_ARGS_NULL; 53 54 struct htrdr_accum radiance; 55 struct htrdr_accum time; 56 struct htrdr_planets* cmd; 57 struct planets_pixel_xwave* pixel = data; 58 size_t isamp = 0; 59 ASSERT(htrdr && htrdr_draw_pixel_args_check(args) && data); 60 (void)htrdr; 61 62 cmd = args->context; 63 ASSERT(cmd); 64 ASSERT(cmd->spectral_domain.type == HTRDR_SPECTRAL_SW 65 || cmd->spectral_domain.type == HTRDR_SPECTRAL_LW); 66 ASSERT(cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_IMAGE); 67 68 /* Reset accumulators */ 69 radiance = HTRDR_ACCUM_NULL; 70 time = HTRDR_ACCUM_NULL; 71 72 FOR_EACH(isamp, 0, args->spp) { 73 struct time t0, t1; 74 struct scam_sample sample = SCAM_SAMPLE_NULL; 75 struct scam_ray ray = SCAM_RAY_NULL; 76 double weight; 77 double r0, r1, r2; 78 double wlen[2]; /* Sampled wavelength */ 79 double pdf; 80 size_t iband[2]; 81 size_t iquad; 82 double usec; 83 84 /* Begin the registration of the time spent to in the realisation */ 85 time_current(&t0); 86 87 /* Sample a position into the pixel, in the normalized image plane */ 88 sample.film[0] = (double)args->pixel_coord[0]+ssp_rng_canonical(args->rng); 89 sample.film[1] = (double)args->pixel_coord[1]+ssp_rng_canonical(args->rng); 90 sample.film[0] *= args->pixel_normalized_size[0]; 91 sample.film[1] *= args->pixel_normalized_size[1]; 92 sample.lens[0] = ssp_rng_canonical(args->rng); 93 sample.lens[1] = ssp_rng_canonical(args->rng); 94 95 /* Generate a camera ray */ 96 scam_generate_ray(cmd->camera, &sample, &ray); 97 98 r0 = ssp_rng_canonical(args->rng); 99 r1 = ssp_rng_canonical(args->rng); 100 r2 = ssp_rng_canonical(args->rng); 101 102 /* Sample a wavelength */ 103 switch(cmd->spectral_domain.type) { 104 case HTRDR_SPECTRAL_LW: 105 wlen[0] = htrdr_ran_wlen_planck_sample(cmd->planck, r0, r1, &pdf); 106 break; 107 case HTRDR_SPECTRAL_SW: 108 if(htrdr_planets_source_does_radiance_vary_spectrally(cmd->source)) { 109 wlen[0] = htrdr_ran_wlen_discrete_sample(cmd->discrete, r0, r1, &pdf); 110 } else { 111 wlen[0] = htrdr_ran_wlen_planck_sample(cmd->planck, r0, r1, &pdf); 112 } 113 break; 114 default: FATAL("Unreachable code\n"); break; 115 116 } 117 wlen[1] = wlen[0]; 118 pdf *= 1.e9; /* Transform the pdf from nm⁻¹ to m⁻¹ */ 119 120 /* Find the band the wavelength belongs to and sample a quadrature point */ 121 RNATM(find_bands(cmd->atmosphere, wlen, iband)); 122 RNATM(band_sample_quad_pt(cmd->atmosphere, r2, iband[0], &iquad)); 123 ASSERT(iband[0] == iband[1]); 124 125 /* Compute the radiance in W/m²/sr/m */ 126 d3_set(rad_args.path_org, ray.org); 127 d3_set(rad_args.path_dir, ray.dir); 128 rad_args.rng = args->rng; 129 rad_args.ithread = args->ithread; 130 rad_args.wlen = wlen[0]; 131 rad_args.iband = iband[0]; 132 rad_args.iquad = iquad; 133 weight = planets_compute_radiance(cmd, &rad_args); 134 ASSERT(weight >= 0); 135 136 weight /= pdf; /* In W/m²/sr */ 137 138 /* End of time recording per realisation */ 139 time_sub(&t0, time_current(&t1), &t0); 140 usec = (double)time_val(&t0, TIME_NSEC) * 0.001; 141 142 /* Update the radiance accumulator */ 143 radiance.sum_weights += weight; 144 radiance.sum_weights_sqr += weight*weight; 145 radiance.nweights += 1; 146 147 /* Update the per realisation time accumulator */ 148 time.sum_weights += usec; 149 time.sum_weights_sqr += usec*usec; 150 time.nweights += 1; 151 } 152 153 /* Flush pixel data */ 154 pixel->radiance = radiance; 155 pixel->time = time; 156 157 if(cmd->spectral_domain.type == HTRDR_SPECTRAL_SW) { 158 pixel->radiance_temperature.E = 0; 159 pixel->radiance_temperature.SE = 0; 160 } else { 161 struct htrdr_estimate L; 162 163 /* Wavelength in meters */ 164 const double wmin_m = cmd->spectral_domain.wlen_range[0] * 1.e-9; 165 const double wmax_m = cmd->spectral_domain.wlen_range[1] * 1.e-9; 166 167 /* Brightness temperatures in W/m²/sr */ 168 double T, Tmin, Tmax; 169 170 htrdr_accum_get_estimation(&radiance, &L); 171 172 T = htrdr_radiance_temperature(cmd->htrdr, wmin_m, wmax_m, L.E); 173 Tmin = htrdr_radiance_temperature(cmd->htrdr, wmin_m, wmax_m, L.E - L.SE); 174 Tmax = htrdr_radiance_temperature(cmd->htrdr, wmin_m, wmax_m, L.E + L.SE); 175 176 pixel->radiance_temperature.E = T; 177 pixel->radiance_temperature.SE = Tmax - Tmin; 178 } 179 } 180 181 static void 182 draw_pixel_image 183 (struct htrdr* htrdr, 184 const struct htrdr_draw_pixel_args* args, 185 void* data) 186 { 187 struct planets_compute_radiance_args rad_args = 188 PLANETS_COMPUTE_RADIANCE_ARGS_NULL; 189 190 struct htrdr_accum XYZ[3]; /* X, Y, and Z */ 191 struct htrdr_accum time; 192 struct htrdr_planets* cmd; 193 struct planets_pixel_image* pixel = data; 194 size_t ichannel; 195 ASSERT(htrdr && htrdr_draw_pixel_args_check(args) && data); 196 (void)htrdr; 197 198 cmd = args->context; 199 ASSERT(cmd); 200 ASSERT(cmd->spectral_domain.type == HTRDR_SPECTRAL_SW_CIE_XYZ); 201 ASSERT(cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_IMAGE); 202 203 /* Reset accumulators */ 204 XYZ[0] = HTRDR_ACCUM_NULL; 205 XYZ[1] = HTRDR_ACCUM_NULL; 206 XYZ[2] = HTRDR_ACCUM_NULL; 207 time = HTRDR_ACCUM_NULL; 208 209 FOR_EACH(ichannel, 0, 3) { 210 size_t isamp; 211 212 FOR_EACH(isamp, 0, args->spp) { 213 struct time t0, t1; 214 struct scam_sample sample = SCAM_SAMPLE_NULL; 215 struct scam_ray ray = SCAM_RAY_NULL; 216 double weight; 217 double r0, r1, r2; 218 double wlen[2]; /* Sampled wavelength into the spectral band */ 219 double pdf; 220 size_t iband[2]; /* Sampled spectral band */ 221 size_t iquad; /* Sampled quadrature point into the spectral band */ 222 double usec; 223 224 /* Begin the registration of the time spent to in the realisation */ 225 time_current(&t0); 226 227 /* Sample a position into the pixel, in the normalized image plane */ 228 sample.film[0] = (double)args->pixel_coord[0]+ssp_rng_canonical(args->rng); 229 sample.film[1] = (double)args->pixel_coord[1]+ssp_rng_canonical(args->rng); 230 sample.film[0] *= args->pixel_normalized_size[0]; 231 sample.film[1] *= args->pixel_normalized_size[1]; 232 sample.lens[0] = ssp_rng_canonical(args->rng); 233 sample.lens[1] = ssp_rng_canonical(args->rng); 234 235 /* Generate a camera ray */ 236 SCAM(generate_ray(cmd->camera, &sample, &ray)); 237 238 r0 = ssp_rng_canonical(args->rng); 239 r1 = ssp_rng_canonical(args->rng); 240 r2 = ssp_rng_canonical(args->rng); 241 242 /* Sample a wavelength according to the CIE XYZ color space */ 243 switch(ichannel) { 244 case 0: 245 wlen[0] = htrdr_ran_wlen_cie_xyz_sample_X(cmd->cie, r0, r1, &pdf); 246 break; 247 case 1: 248 wlen[0] = htrdr_ran_wlen_cie_xyz_sample_Y(cmd->cie, r0, r1, &pdf); 249 break; 250 case 2: 251 wlen[0] = htrdr_ran_wlen_cie_xyz_sample_Z(cmd->cie, r0, r1, &pdf); 252 break; 253 default: FATAL("Unreachable code.\n"); break; 254 } 255 pdf *= 1.e9; /* Transform the pdf from nm⁻¹ to m⁻¹ */ 256 257 /* Find the band the wavelength belongs to and sample a quadrature point */ 258 wlen[1] = wlen[0]; 259 RNATM(find_bands(cmd->atmosphere, wlen, iband)); 260 RNATM(band_sample_quad_pt(cmd->atmosphere, r2, iband[0], &iquad)); 261 ASSERT(iband[0] == iband[1]); 262 263 /* Compute the radiance in W/m²/sr/m */ 264 d3_set(rad_args.path_org, ray.org); 265 d3_set(rad_args.path_dir, ray.dir); 266 rad_args.rng = args->rng; 267 rad_args.ithread = args->ithread; 268 rad_args.wlen = wlen[0]; 269 rad_args.iband = iband[0]; 270 rad_args.iquad = iquad; 271 weight = planets_compute_radiance(cmd, &rad_args); 272 ASSERT(weight >= 0); 273 274 weight /= pdf; /* In W/m²/sr */ 275 276 /* End of time recording per realisation */ 277 time_sub(&t0, time_current(&t1), &t0); 278 usec = (double)time_val(&t0, TIME_NSEC) * 0.001; 279 280 /* Update pixel channel accumulators */ 281 XYZ[ichannel].sum_weights += weight; 282 XYZ[ichannel].sum_weights_sqr += weight*weight; 283 XYZ[ichannel].nweights += 1; 284 285 /* Update the per realisation time accumulator */ 286 time.sum_weights += usec; 287 time.sum_weights_sqr += usec*usec; 288 time.nweights += 1; 289 } 290 } 291 292 /* Flush pixel data */ 293 htrdr_accum_get_estimation(XYZ+0, &pixel->X); 294 htrdr_accum_get_estimation(XYZ+1, &pixel->Y); 295 htrdr_accum_get_estimation(XYZ+2, &pixel->Z); 296 pixel->time = time; 297 } 298 299 300 static INLINE void 301 write_accum 302 (const struct htrdr_accum* acc, /* Accum to dump */ 303 struct htrdr_accum* out_acc, /* May be NULL */ 304 FILE* stream) 305 { 306 ASSERT(acc && stream); 307 308 if(acc->nweights == 0) { 309 fprintf(stream, "0 0 "); 310 } else { 311 struct htrdr_estimate estimate = HTRDR_ESTIMATE_NULL; 312 313 htrdr_accum_get_estimation(acc, &estimate); 314 fprintf(stream, "%g %g ", estimate.E, estimate.SE); 315 316 if(out_acc) { 317 out_acc->sum_weights += acc->sum_weights; 318 out_acc->sum_weights_sqr += acc->sum_weights_sqr; 319 out_acc->nweights += acc->nweights; 320 } 321 } 322 } 323 324 static INLINE void 325 write_pixel_image 326 (const struct planets_pixel_image* pix, 327 struct htrdr_accum* time_acc, /* May be NULL */ 328 FILE* stream) 329 { 330 ASSERT(pix && stream); 331 fprintf(stream, "%g %g ", pix->X.E, pix->X.SE); 332 fprintf(stream, "%g %g ", pix->Y.E, pix->Y.SE); 333 fprintf(stream, "%g %g ", pix->Z.E, pix->Z.SE); 334 write_accum(&pix->time, time_acc, stream); 335 fprintf(stream, "\n"); 336 } 337 338 static INLINE void 339 write_pixel_xwave 340 (const struct planets_pixel_xwave* pix, 341 struct htrdr_accum* radiance_acc, /* May be NULL */ 342 struct htrdr_accum* time_acc, /* May be NULL */ 343 FILE* stream) 344 { 345 ASSERT(pix && stream); 346 fprintf(stream, "%g %g ", 347 pix->radiance_temperature.E, 348 pix->radiance_temperature.SE); 349 write_accum(&pix->radiance, radiance_acc, stream); 350 fprintf(stream, "0 0 "); 351 write_accum(&pix->time, time_acc, stream); 352 fprintf(stream, "\n"); 353 } 354 355 static res_T 356 write_buffer 357 (struct htrdr_planets* cmd, 358 struct htrdr_buffer* buf, 359 struct htrdr_accum* radiance_acc, /* May be NULL */ 360 struct htrdr_accum* time_acc, 361 FILE* stream) 362 { 363 struct htrdr_pixel_format pixfmt; 364 struct htrdr_buffer_layout layout; 365 size_t x, y; 366 ASSERT(cmd && buf && time_acc && stream); 367 ASSERT(cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_IMAGE); 368 369 planets_get_pixel_format(cmd, &pixfmt); 370 htrdr_buffer_get_layout(buf, &layout); 371 CHK(pixfmt.size == layout.elmt_size); 372 373 fprintf(stream, "%lu %lu\n", layout.width, layout.height); 374 *time_acc = HTRDR_ACCUM_NULL; 375 376 FOR_EACH(y, 0, layout.height) { 377 FOR_EACH(x, 0, layout.width) { 378 void* pix_raw = htrdr_buffer_at(buf, x, y); 379 ASSERT(IS_ALIGNED(pix_raw, pixfmt.alignment)); 380 381 switch(cmd->spectral_domain.type) { 382 case HTRDR_SPECTRAL_LW: 383 case HTRDR_SPECTRAL_SW: 384 write_pixel_xwave(pix_raw, radiance_acc, time_acc, stream); 385 break; 386 case HTRDR_SPECTRAL_SW_CIE_XYZ: 387 write_pixel_image(pix_raw, time_acc, stream); 388 break; 389 default: FATAL("Unreachable code\n"); break; 390 } 391 } 392 } 393 return RES_OK; 394 } 395 396 /******************************************************************************* 397 * Local functions 398 ******************************************************************************/ 399 res_T 400 planets_draw_map(struct htrdr_planets* cmd) 401 { 402 struct htrdr_draw_map_args args = HTRDR_DRAW_MAP_ARGS_NULL; 403 struct htrdr_estimate path_time = HTRDR_ESTIMATE_NULL; 404 struct htrdr_accum path_time_acc = HTRDR_ACCUM_NULL; 405 struct htrdr_accum radiance_acc = HTRDR_ACCUM_NULL; 406 res_T res = RES_OK; 407 ASSERT(cmd && cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_IMAGE); 408 409 args.buffer_layout = cmd->buf_layout; 410 args.spp = cmd->spp; 411 args.context = cmd; 412 switch(cmd->spectral_domain.type) { 413 case HTRDR_SPECTRAL_LW: 414 case HTRDR_SPECTRAL_SW: 415 args.draw_pixel = draw_pixel_xwave; 416 break; 417 case HTRDR_SPECTRAL_SW_CIE_XYZ: 418 args.draw_pixel = draw_pixel_image; 419 break; 420 default: FATAL("Unreachable code\n"); break; 421 } 422 423 res = htrdr_draw_map(cmd->htrdr, &args, cmd->buf); 424 if(res != RES_OK) goto error; 425 426 /* Nothing more to do on non master processes */ 427 if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit; 428 429 /* Write output image */ 430 res = write_buffer(cmd, cmd->buf, &radiance_acc, &path_time_acc, cmd->output); 431 if(res != RES_OK) goto error; 432 433 CHK(fflush(cmd->output) == 0); 434 435 /* Log time per realisation */ 436 htrdr_accum_get_estimation(&path_time_acc, &path_time); 437 htrdr_log(cmd->htrdr, "Time per radiative path (in µs): %g +/- %g\n", 438 path_time.E, path_time.SE); 439 440 /* Log measured radiance on the whole image */ 441 if(cmd->spectral_domain.type == HTRDR_SPECTRAL_LW 442 || cmd->spectral_domain.type == HTRDR_SPECTRAL_SW) { 443 struct htrdr_estimate L; 444 double omega; /* Solid angle of the camera */ 445 enum scam_type type = SCAM_NONE; 446 447 htrdr_accum_get_estimation(&radiance_acc, &L); 448 htrdr_log(cmd->htrdr, "Radiance in W/m²/sr: %g +/- %g\n", L.E, L.SE); 449 450 SCAM(get_type(cmd->camera, &type)); 451 if(type == SCAM_PERSPECTIVE) { 452 SCAM(perspective_get_solid_angle(cmd->camera, &omega)); 453 454 htrdr_log(cmd->htrdr, "Radiance in W/m² (solid angle = %g sr): %g +/- %g\n", 455 omega, L.E*omega, L.SE*omega); 456 } 457 } 458 459 exit: 460 return res; 461 error: 462 goto exit; 463 }