htrdr_combustion_draw_map.c (10383B)
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 "combustion/htrdr_combustion_c.h" 25 26 #include "core/htrdr_accum.h" 27 #include "core/htrdr_draw_map.h" 28 #include "core/htrdr_log.h" 29 #include "core/htrdr_rectangle.h" 30 31 #include <star/scam.h> 32 #include <star/ssp.h> 33 34 #include <rsys/clock_time.h> 35 36 /******************************************************************************* 37 * Helper functions 38 ******************************************************************************/ 39 static void 40 draw_pixel_image 41 (struct htrdr* htrdr, 42 const struct htrdr_draw_pixel_args* args, 43 void* data) 44 { 45 struct htrdr_accum radiance = HTRDR_ACCUM_NULL; 46 struct htrdr_accum time = HTRDR_ACCUM_NULL; 47 struct htrdr_combustion* cmd = NULL; 48 struct combustion_pixel_image* pixel = NULL; 49 size_t isamp; 50 res_T res = RES_OK; 51 ASSERT(htrdr && htrdr_draw_pixel_args_check(args) && data); 52 (void)htrdr; 53 54 cmd = args->context; 55 pixel = data; 56 57 FOR_EACH(isamp, 0, args->spp) { 58 struct time t0, t1; 59 struct scam_sample sample = SCAM_SAMPLE_NULL; 60 struct scam_ray ray = SCAM_RAY_NULL; 61 double weight; 62 double usec; 63 64 /* Begin the registration of the time spent in the realisation */ 65 time_current(&t0); 66 67 /* Sample a position into the pixel, in the normalized image plane */ 68 sample.film[0] = (double)args->pixel_coord[0]+ssp_rng_canonical(args->rng); 69 sample.film[1] = (double)args->pixel_coord[1]+ssp_rng_canonical(args->rng); 70 sample.film[0] *= args->pixel_normalized_size[0]; 71 sample.film[1] *= args->pixel_normalized_size[1]; 72 sample.lens[0] = ssp_rng_canonical(args->rng); 73 sample.lens[1] = ssp_rng_canonical(args->rng); 74 75 /* Generate a camera ray */ 76 scam_generate_ray(cmd->camera, &sample, &ray); 77 78 /* Backward trace the path */ 79 res = combustion_compute_radiance_sw(cmd, args->ithread, args->rng, 80 ray.org, ray.dir, &weight); 81 if(res != RES_OK) continue; /* Reject the path */ 82 83 /* End the registration of the per realisation time */ 84 time_sub(&t0, time_current(&t1), &t0); 85 usec = (double)time_val(&t0, TIME_NSEC) * 0.001; 86 87 /* Update the pixel accumulator */ 88 radiance.sum_weights += weight; 89 radiance.sum_weights_sqr += weight*weight; 90 radiance.nweights += 1; 91 92 /* Update the pixel accumulator of per realisation time */ 93 time.sum_weights += usec; 94 time.sum_weights_sqr += usec*usec; 95 time.nweights += 1; 96 } 97 98 /* Compute the estimation of the pixel radiance */ 99 htrdr_accum_get_estimation(&radiance, &pixel->radiance); 100 101 /* Save the per realisation time */ 102 pixel->time = time; 103 } 104 105 static void 106 draw_pixel_flux 107 (struct htrdr* htrdr, 108 const struct htrdr_draw_pixel_args* args, 109 void* data) 110 { 111 struct htrdr_accum flux = HTRDR_ACCUM_NULL; 112 struct htrdr_accum time = HTRDR_ACCUM_NULL; 113 struct htrdr_combustion* cmd = NULL; 114 struct combustion_pixel_flux* pixel = NULL; 115 size_t isamp; 116 res_T res = RES_OK; 117 ASSERT(htrdr && htrdr_draw_pixel_args_check(args) && data); 118 (void)htrdr; 119 120 cmd = args->context; 121 pixel = data; 122 123 FOR_EACH(isamp, 0, args->spp) { 124 struct time t0, t1; 125 double pix_samp[2]; 126 double ray_org[3]; 127 double ray_dir[3]; 128 double normal[3]; 129 double weight; 130 double usec; 131 132 /* Begin the registration of the time spent in the realisation */ 133 time_current(&t0); 134 135 /* Sample a position into the pixel, in the normalized image plane */ 136 pix_samp[0] = (double)args->pixel_coord[0] + ssp_rng_canonical(args->rng); 137 pix_samp[1] = (double)args->pixel_coord[1] + ssp_rng_canonical(args->rng); 138 pix_samp[0] *= args->pixel_normalized_size[0]; 139 pix_samp[1] *= args->pixel_normalized_size[1]; 140 141 /* Retrieve the world space position of pix_samp */ 142 htrdr_rectangle_sample_pos(cmd->flux_map, pix_samp, ray_org); 143 144 /* Sample a ray direction */ 145 htrdr_rectangle_get_normal(cmd->flux_map, normal); 146 ssp_ran_hemisphere_cos(args->rng, normal, ray_dir, NULL); 147 148 /* Backward trace the path */ 149 res = combustion_compute_radiance_sw(cmd, args->ithread, 150 args->rng, ray_org, ray_dir, &weight); 151 if(res != RES_OK) continue; /* Reject the path */ 152 weight *= PI; /* Transform form W/m^2/sr to W/m^2 */ 153 154 /* End the registration of the per realisation time */ 155 time_sub(&t0, time_current(&t1), &t0); 156 usec = (double)time_val(&t0, TIME_NSEC) * 0.001; 157 158 /* Update the pixel accumulator */ 159 flux.sum_weights += weight; 160 flux.sum_weights_sqr += weight*weight; 161 flux.nweights += 1; 162 163 /* Update the pixel accumulator of per realisation time */ 164 time.sum_weights += usec; 165 time.sum_weights_sqr += usec*usec; 166 time.nweights += 1; 167 } 168 169 /* Write the accumulators */ 170 pixel->flux = flux; 171 pixel->time = time; 172 } 173 174 static INLINE void 175 dump_accum 176 (const struct htrdr_accum* acc, /* Accum to dump */ 177 struct htrdr_accum* out_acc, /* May be NULL */ 178 FILE* stream) 179 { 180 ASSERT(acc && stream); 181 182 if(acc->nweights == 0) { 183 fprintf(stream, "0 0 "); 184 } else { 185 struct htrdr_estimate estimate = HTRDR_ESTIMATE_NULL; 186 187 htrdr_accum_get_estimation(acc, &estimate); 188 fprintf(stream, "%g %g ", estimate.E, estimate.SE); 189 190 if(out_acc) { 191 out_acc->sum_weights += acc->sum_weights; 192 out_acc->sum_weights_sqr += acc->sum_weights_sqr; 193 out_acc->nweights += acc->nweights; 194 } 195 } 196 } 197 198 static void 199 dump_pixel_flux 200 (const struct combustion_pixel_flux* pix, 201 struct htrdr_accum* time_acc, /* May be NULL */ 202 struct htrdr_accum* flux_acc, /* May be NULL */ 203 FILE* stream) 204 { 205 ASSERT(pix && stream); 206 dump_accum(&pix->flux, flux_acc, stream); 207 fprintf(stream, "0 0 0 0 "); 208 dump_accum(&pix->time, time_acc, stream); 209 fprintf(stream, "\n"); 210 } 211 212 static void 213 dump_pixel_image 214 (const struct combustion_pixel_image* pix, 215 struct htrdr_accum* time_acc, /* May be NULL */ 216 FILE* stream) 217 { 218 ASSERT(pix && stream); 219 fprintf(stream, "%g %g ", pix->radiance.E, pix->radiance.SE); 220 fprintf(stream, "0 0 0 0 "); 221 dump_accum(&pix->time, time_acc, stream); 222 fprintf(stream, "\n"); 223 } 224 225 static res_T 226 dump_buffer 227 (struct htrdr_combustion* cmd, 228 struct htrdr_buffer* buf, 229 struct htrdr_accum* time_acc, /* May be NULL */ 230 struct htrdr_accum* flux_acc, /* May be NULL */ 231 FILE* stream) 232 { 233 struct htrdr_pixel_format pixfmt; 234 struct htrdr_buffer_layout layout; 235 size_t x, y; 236 ASSERT(cmd && buf && stream); 237 238 combustion_get_pixel_format(cmd, &pixfmt); 239 htrdr_buffer_get_layout(buf, &layout); 240 CHK(pixfmt.size == layout.elmt_size); 241 242 fprintf(stream, "%lu %lu\n", layout.width, layout.height); 243 244 if(time_acc) *time_acc = HTRDR_ACCUM_NULL; 245 246 FOR_EACH(y, 0, layout.height) { 247 FOR_EACH(x, 0, layout.width) { 248 void* pix_raw = htrdr_buffer_at(buf, x, y); 249 ASSERT(IS_ALIGNED(pix_raw, pixfmt.alignment)); 250 if(cmd->output_type == HTRDR_COMBUSTION_ARGS_OUTPUT_FLUX_MAP) { 251 const struct combustion_pixel_flux* pix = pix_raw; 252 dump_pixel_flux(pix, time_acc, flux_acc, stream); 253 } else if(cmd->output_type == HTRDR_COMBUSTION_ARGS_OUTPUT_IMAGE) { 254 const struct combustion_pixel_image* pix = pix_raw; 255 dump_pixel_image(pix, time_acc, stream); 256 } else { 257 FATAL("Unreachable code.\n"); 258 } 259 } 260 fprintf(stream, "\n"); 261 } 262 return RES_OK; 263 } 264 265 /******************************************************************************* 266 * Local functions 267 ******************************************************************************/ 268 res_T 269 combustion_draw_map(struct htrdr_combustion* cmd) 270 { 271 struct htrdr_draw_map_args args = HTRDR_DRAW_MAP_ARGS_NULL; 272 struct htrdr_accum path_time_acc = HTRDR_ACCUM_NULL; 273 struct htrdr_accum flux_acc = HTRDR_ACCUM_NULL; 274 struct htrdr_estimate path_time = HTRDR_ESTIMATE_NULL; 275 struct htrdr_estimate flux = HTRDR_ESTIMATE_NULL; 276 res_T res = RES_OK; 277 ASSERT(cmd); 278 279 /* Setup the draw map input arguments */ 280 switch(cmd->output_type) { 281 case HTRDR_COMBUSTION_ARGS_OUTPUT_IMAGE: 282 args.draw_pixel = draw_pixel_image; 283 break; 284 case HTRDR_COMBUSTION_ARGS_OUTPUT_FLUX_MAP: 285 args.draw_pixel = draw_pixel_flux; 286 break; 287 default: FATAL("Unreachable code.\n"); break; 288 } 289 args.buffer_layout = cmd->buf_layout; 290 args.spp = cmd->spp; 291 args.context = cmd; 292 293 res = htrdr_draw_map(cmd->htrdr, &args, cmd->buf); 294 if(res != RES_OK) goto error; 295 296 /* Nothing more to do on non master processes */ 297 if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit; 298 299 /* Write buffer to output */ 300 res = dump_buffer(cmd, cmd->buf, &path_time_acc, &flux_acc, cmd->output); 301 if(res != RES_OK) goto error; 302 303 htrdr_accum_get_estimation(&path_time_acc, &path_time); 304 htrdr_log(cmd->htrdr, 305 "Time per radiative path (in micro seconds): %g +/- %g\n", 306 path_time.E, path_time.SE); 307 308 if(cmd->output_type == HTRDR_COMBUSTION_ARGS_OUTPUT_FLUX_MAP) { 309 double map_size[2]; 310 double area; 311 312 htrdr_accum_get_estimation(&flux_acc, &flux); 313 htrdr_log(cmd->htrdr, 314 "Radiative flux density (in W/m^2): %g +/- %g\n", 315 flux.E, flux.SE); 316 317 htrdr_rectangle_get_size(cmd->flux_map, map_size); 318 area = map_size[0] * map_size[1]; 319 htrdr_log(cmd->htrdr, 320 "Radiative flux (in W): %g +/- %g\n", 321 flux.E*area, flux.SE*area); 322 } 323 324 exit: 325 return res; 326 error: 327 goto exit; 328 } 329