htrdr

Solving radiative transfer in heterogeneous media
git clone git://git.meso-star.fr/htrdr.git
Log | Files | Refs | README | LICENSE

htrdr_combustion.c (21306B)


      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.h"
     25 #include "combustion/htrdr_combustion_args.h"
     26 #include "combustion/htrdr_combustion_c.h"
     27 #include "combustion/htrdr_combustion_laser.h"
     28 
     29 #include "core/htrdr.h"
     30 #include "core/htrdr_log.h"
     31 #include "core/htrdr_geometry.h"
     32 #include "core/htrdr_materials.h"
     33 #include "core/htrdr_rectangle.h"
     34 
     35 #include <astoria/atrstm.h>
     36 
     37 #include <star/scam.h>
     38 #include <star/ssf.h>
     39 
     40 #include <rsys/cstr.h>
     41 #include <rsys/double3.h>
     42 #include <rsys/mem_allocator.h>
     43 
     44 /*******************************************************************************
     45  * Helper functions
     46  ******************************************************************************/
     47 static void
     48 release_phase_functions(struct htrdr_combustion* cmd)
     49 {
     50   size_t i;
     51   ASSERT(cmd);
     52 
     53   if(!cmd->phase_functions) return; /* Nothing to release */
     54 
     55   FOR_EACH(i, 0, htrdr_get_threads_count(cmd->htrdr)) {
     56     if(cmd->phase_functions[i]) {
     57       SSF(phase_ref_put(cmd->phase_functions[i]));
     58     }
     59   }
     60   MEM_RM(htrdr_get_allocator(cmd->htrdr), cmd->phase_functions);
     61   cmd->phase_functions = NULL;
     62 }
     63 
     64 static res_T
     65 setup_output
     66   (struct htrdr_combustion* cmd,
     67    const struct htrdr_combustion_args* args)
     68 {
     69   const char* output_name = NULL;
     70   res_T res = RES_OK;
     71   ASSERT(cmd && args);
     72 
     73   if(htrdr_get_mpi_rank(cmd->htrdr) != 0) {
     74     /* No output stream on non master processes */
     75     cmd->output = NULL;
     76     output_name = "<null>";
     77 
     78   } else if(!args->path_output) {
     79     /* Write results to standard output when no destination file is defined */
     80     cmd->output = stdout;
     81     output_name = "<stdout>";
     82 
     83   } else {
     84     /* Open the output stream */
     85     res = htrdr_open_output_stream(cmd->htrdr, args->path_output,
     86       0/*read*/, args->force_overwriting, &cmd->output);
     87     if(res != RES_OK) goto error;
     88     output_name = args->path_output;
     89   }
     90 
     91   /* Setup the output name */
     92   res = str_set(&cmd->output_name, output_name);
     93   if(res != RES_OK) {
     94     htrdr_log_err(cmd->htrdr,
     95       "Could not store the name of the output stream `%s' -- %s.\n",
     96       output_name, res_to_cstr(res));
     97     goto error;
     98   }
     99 
    100 exit:
    101   return res;
    102 error:
    103   str_clear(&cmd->output_name);
    104   if(cmd->output && cmd->output != stdout) {
    105     CHK(fclose(cmd->output) == 0);
    106     cmd->output = NULL;
    107   }
    108   goto exit;
    109 }
    110 
    111 static res_T
    112 setup_simd
    113   (struct htrdr_combustion* cmd,
    114    const struct htrdr_combustion_args* args)
    115 {
    116   struct ssf_info ssf_info = SSF_INFO_NULL;
    117   ASSERT(cmd && args);
    118 
    119   cmd->rdgfa_simd = SSF_SIMD_NONE;
    120 
    121   if(args->phase_func_type != HTRDR_COMBUSTION_ARGS_PHASE_FUNC_RDGFA
    122   || args->use_simd == 0)
    123     return RES_OK; /* Nothing to do */
    124 
    125   /* Check SIMD support for the RDG-FA phase function */
    126   ssf_get_info(&ssf_info);
    127   if(ssf_info.simd_256) {
    128     htrdr_log(cmd->htrdr,
    129       "Use the SIMD-256 instruction set for the RDG-FA phase function.\n");
    130     cmd->rdgfa_simd = SSF_SIMD_256;
    131   } else if(ssf_info.simd_128) {
    132     htrdr_log(cmd->htrdr,
    133       "Use the SIMD-128 instruction set for the RDG-FA phase function.\n");
    134     cmd->rdgfa_simd = SSF_SIMD_128;
    135   } else {
    136     htrdr_log_warn(cmd->htrdr,
    137       "Cannot use SIMD for the RDG-FA phase function: the "
    138       "Star-ScatteringFunction library was compiled without SIMD support.\n");
    139     cmd->rdgfa_simd = SSF_SIMD_NONE;
    140   }
    141   return RES_OK;
    142 }
    143 
    144 static res_T
    145 setup_geometry
    146   (struct htrdr_combustion* cmd,
    147    const struct htrdr_combustion_args* args)
    148 {
    149   res_T res = RES_OK;
    150 
    151   if(!args->geom.path_obj) goto exit;
    152   ASSERT(args->geom.path_mats);
    153 
    154   res = htrdr_materials_create(cmd->htrdr, args->geom.path_mats, &cmd->mats);
    155   if(res != RES_OK) goto error;
    156 
    157   res = htrdr_geometry_create
    158     (cmd->htrdr, args->geom.path_obj, cmd->mats, &cmd->geom);
    159   if(res != RES_OK) goto error;
    160 
    161 exit:
    162   return res;
    163 error:
    164   if(cmd->mats) { htrdr_materials_ref_put(cmd->mats); cmd->mats = NULL; }
    165   if(cmd->geom) { htrdr_geometry_ref_put(cmd->geom); cmd->geom = NULL; }
    166   goto exit;
    167 }
    168 
    169 static res_T
    170 setup_camera_orthographic
    171   (struct htrdr_combustion* cmd,
    172    const struct htrdr_combustion_args* args)
    173 {
    174   struct scam_orthographic_args cam_args = SCAM_ORTHOGRAPHIC_ARGS_DEFAULT;
    175   ASSERT(cmd && args && args->image.definition[0] && args->image.definition[1]);
    176   ASSERT(cmd->output_type == HTRDR_COMBUSTION_ARGS_OUTPUT_IMAGE);
    177   ASSERT(args->cam_type == HTRDR_ARGS_CAMERA_ORTHOGRAPHIC);
    178 
    179   d3_set(cam_args.position, args->cam_ortho.position);
    180   d3_set(cam_args.target, args->cam_ortho.target);
    181   d3_set(cam_args.up, args->cam_ortho.up);
    182   cam_args.height = args->cam_ortho.height;
    183   cam_args.aspect_ratio =
    184     (double)args->image.definition[0]
    185   / (double)args->image.definition[1];
    186 
    187   return scam_create_orthographic
    188     (htrdr_get_logger(cmd->htrdr),
    189      htrdr_get_allocator(cmd->htrdr),
    190      htrdr_get_verbosity_level(cmd->htrdr),
    191      &cam_args,
    192      &cmd->camera);
    193 }
    194 
    195 static res_T
    196 setup_camera_perspective
    197   (struct htrdr_combustion* cmd,
    198    const struct htrdr_combustion_args* args)
    199 {
    200   struct scam_perspective_args cam_args = SCAM_PERSPECTIVE_ARGS_DEFAULT;
    201   ASSERT(cmd && args && args->image.definition[0] && args->image.definition[1]);
    202   ASSERT(cmd->output_type == HTRDR_COMBUSTION_ARGS_OUTPUT_IMAGE);
    203   ASSERT(args->cam_type == HTRDR_ARGS_CAMERA_PERSPECTIVE);
    204 
    205   d3_set(cam_args.position, args->cam_persp.position);
    206   d3_set(cam_args.target, args->cam_persp.target);
    207   d3_set(cam_args.up, args->cam_persp.up);
    208   cam_args.aspect_ratio =
    209     (double)args->image.definition[0]
    210   / (double)args->image.definition[1];
    211   cam_args.field_of_view = MDEG2RAD(args->cam_persp.fov_y);
    212   cam_args.lens_radius = args->cam_persp.lens_radius;
    213   cam_args.focal_distance = args->cam_persp.focal_dst;
    214 
    215   return scam_create_perspective
    216     (htrdr_get_logger(cmd->htrdr),
    217      htrdr_get_allocator(cmd->htrdr),
    218      htrdr_get_verbosity_level(cmd->htrdr),
    219      &cam_args,
    220      &cmd->camera);
    221 }
    222 
    223 static res_T
    224 setup_camera
    225   (struct htrdr_combustion* cmd,
    226    const struct htrdr_combustion_args* args)
    227 {
    228   res_T res = RES_OK;
    229   ASSERT(cmd->output_type == HTRDR_COMBUSTION_ARGS_OUTPUT_IMAGE);
    230   switch(args->cam_type) {
    231     case HTRDR_ARGS_CAMERA_ORTHOGRAPHIC:
    232       res = setup_camera_orthographic(cmd, args);
    233       break;
    234     case HTRDR_ARGS_CAMERA_PERSPECTIVE:
    235       res = setup_camera_perspective(cmd, args);
    236       break;
    237     default: FATAL("Unreachable code.\n"); break;
    238   }
    239   return res;
    240 }
    241 
    242 static res_T
    243 setup_flux_map
    244   (struct htrdr_combustion* cmd,
    245    const struct htrdr_combustion_args* args)
    246 {
    247   ASSERT(cmd && args);
    248   ASSERT(cmd->output_type == HTRDR_COMBUSTION_ARGS_OUTPUT_FLUX_MAP);
    249   return htrdr_rectangle_create
    250     (cmd->htrdr,
    251      args->flux_map.size,
    252      args->flux_map.position,
    253      args->flux_map.target,
    254      args->flux_map.up,
    255      &cmd->flux_map);
    256 }
    257 
    258 static res_T
    259 setup_sensor
    260   (struct htrdr_combustion* cmd,
    261    const struct htrdr_combustion_args* args)
    262 {
    263   res_T res = RES_OK;
    264   switch(cmd->output_type) {
    265     case HTRDR_COMBUSTION_ARGS_OUTPUT_FLUX_MAP:
    266       res = setup_flux_map(cmd, args);
    267       break;
    268     case HTRDR_COMBUSTION_ARGS_OUTPUT_IMAGE:
    269       res = setup_camera(cmd, args);
    270       break;
    271     default: /* Nothing to do */ break;
    272   }
    273   return res;
    274 }
    275 
    276 static res_T
    277 setup_laser
    278   (struct htrdr_combustion* cmd,
    279    const struct htrdr_combustion_args* args)
    280 {
    281   struct htrdr_combustion_laser_create_args laser_args =
    282     HTRDR_COMBUSTION_LASER_CREATE_ARGS_DEFAULT;
    283   ASSERT(cmd && args);
    284   cmd->wavelength = args->wavelength;
    285   laser_args.surface = args->laser;
    286   laser_args.wavelength = args->wavelength;
    287   laser_args.flux_density = args->laser_flux_density;
    288   return htrdr_combustion_laser_create(cmd->htrdr, &laser_args, &cmd->laser);
    289 }
    290 
    291 static res_T
    292 setup_phase_functions(struct htrdr_combustion* cmd)
    293 {
    294   struct mem_allocator* allocator = NULL;
    295   const struct ssf_phase_type* phase_type = NULL;
    296   size_t nthreads;
    297   size_t i;
    298   res_T res = RES_OK;
    299   ASSERT(cmd);
    300 
    301   nthreads = htrdr_get_threads_count(cmd->htrdr);
    302   allocator = htrdr_get_allocator(cmd->htrdr);
    303 
    304   switch(cmd->phase_func_type) {
    305     case HTRDR_COMBUSTION_ARGS_PHASE_FUNC_ISOTROPIC:
    306       htrdr_log(cmd->htrdr, "Use an isotropic phase function.\n");
    307       phase_type = &ssf_phase_hg;
    308       break;
    309     case HTRDR_COMBUSTION_ARGS_PHASE_FUNC_RDGFA:
    310       htrdr_log(cmd->htrdr, "Use the RDG-FA phase function.\n");
    311       phase_type = &ssf_phase_rdgfa;
    312       break;
    313     default: FATAL("Unreachable code.\n"); break;
    314   }
    315 
    316   /* Allocate the list of per thread phase function */
    317   cmd->phase_functions = MEM_CALLOC
    318     (allocator, nthreads, sizeof(*cmd->phase_functions));
    319   if(!cmd->phase_functions) {
    320     htrdr_log_err(cmd->htrdr,
    321       "Could not allocate the per thread RDG-FA phase function.\n");
    322     res = RES_MEM_ERR;
    323     goto error;
    324   }
    325 
    326   /* Create the per thread phase function */
    327   FOR_EACH(i, 0, nthreads) {
    328     res = ssf_phase_create(allocator, phase_type, cmd->phase_functions+i);
    329     if(res != RES_OK) {
    330       htrdr_log_err(cmd->htrdr,
    331         "Could not create the phase function for the thread %lu -- %s.\n",
    332         (unsigned long)i, res_to_cstr(res));
    333       goto error;
    334     }
    335   }
    336 
    337 exit:
    338   return res;
    339 error:
    340   release_phase_functions(cmd);
    341   goto exit;
    342 }
    343 
    344 static res_T
    345 setup_buffer
    346   (struct htrdr_combustion* cmd,
    347    const struct htrdr_combustion_args* args)
    348 {
    349   struct htrdr_pixel_format pixfmt = HTRDR_PIXEL_FORMAT_NULL;
    350   res_T res = RES_OK;
    351   ASSERT(cmd && args);
    352 
    353   if(cmd->output_type != HTRDR_COMBUSTION_ARGS_OUTPUT_FLUX_MAP
    354   && cmd->output_type != HTRDR_COMBUSTION_ARGS_OUTPUT_IMAGE)
    355     goto exit;
    356 
    357   combustion_get_pixel_format(cmd, &pixfmt);
    358 
    359   /* Setup the buffer layout */
    360   cmd->buf_layout.width = args->image.definition[0];
    361   cmd->buf_layout.height = args->image.definition[1];
    362   cmd->buf_layout.pitch = args->image.definition[0] * pixfmt.size;
    363   cmd->buf_layout.elmt_size = pixfmt.size;
    364   cmd->buf_layout.alignment = pixfmt.alignment;
    365 
    366   /* Create the image buffer only on the master process; the image parts
    367    * rendered by the others processes are gathered onto the master process */
    368   if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit;
    369 
    370   res = htrdr_buffer_create(cmd->htrdr, &cmd->buf_layout, &cmd->buf);
    371   if(res != RES_OK) goto error;
    372 
    373 exit:
    374   return res;
    375 error:
    376   if(cmd->buf) { htrdr_buffer_ref_put(cmd->buf); cmd->buf = NULL; }
    377   goto exit;
    378 }
    379 
    380 static res_T
    381 setup_medium
    382   (struct htrdr_combustion* cmd,
    383    const struct htrdr_combustion_args* args)
    384 {
    385   struct atrstm_args atrstm_args = ATRSTM_ARGS_DEFAULT;
    386   res_T res = RES_OK;
    387   ASSERT(cmd && args);
    388 
    389   /* Setup the semi-transaprent medium arguments */
    390   atrstm_args.sth_filename = args->path_tetra;
    391   atrstm_args.atrtp_filename = args->path_therm_props;
    392   atrstm_args.atrri_filename = args->path_refract_ids;
    393   atrstm_args.cache_filename = args->path_cache;
    394   atrstm_args.spectral_type = ATRSTM_SPECTRAL_SW;
    395   atrstm_args.wlen_range[0] = args->wavelength;
    396   atrstm_args.wlen_range[1] = args->wavelength;
    397   atrstm_args.fractal_prefactor = args->fractal_prefactor;
    398   atrstm_args.fractal_dimension = args->fractal_dimension;
    399   atrstm_args.optical_thickness = args->optical_thickness;
    400   atrstm_args.precompute_normals = args->precompute_normals;
    401   atrstm_args.use_simd = args->use_simd;
    402   atrstm_args.nthreads = args->nthreads;
    403   atrstm_args.verbose = args->verbose;
    404 
    405   switch(args->grid.type) {
    406     case HTRDR_COMBUSTION_ARGS_GRID_DEFINITION_AUTO:
    407       atrstm_args.auto_grid_definition = 1;
    408       atrstm_args.auto_grid_definition_hint = args->grid.definition.hint;
    409       break;
    410     case HTRDR_COMBUSTION_ARGS_GRID_DEFINITION_FIXED:
    411       atrstm_args.auto_grid_definition = 0;
    412       atrstm_args.grid_max_definition[0] = args->grid.definition.fixed[0];
    413       atrstm_args.grid_max_definition[1] = args->grid.definition.fixed[1];
    414       atrstm_args.grid_max_definition[2] = args->grid.definition.fixed[2];
    415       break;
    416     default: FATAL("Unreachable code.\n"); break;
    417   }
    418 
    419   /* Here we go! Create the semi-transparent medium */
    420   res = atrstm_create
    421     (htrdr_get_logger(cmd->htrdr),
    422      htrdr_get_allocator(cmd->htrdr),
    423      &atrstm_args,
    424      &cmd->medium);
    425   if(res != RES_OK) goto error;
    426 
    427 exit:
    428   return res;
    429 error:
    430   if(cmd->medium) { ATRSTM(ref_put(cmd->medium)); cmd->medium = NULL; }
    431   goto exit;
    432 }
    433 
    434 static res_T
    435 dump_volumetric_acceleration_structure(struct htrdr_combustion* cmd)
    436 {
    437   struct atrstm_dump_svx_octree_args args = ATRSTM_DUMP_SVX_OCTREE_ARGS_DEFAULT;
    438   res_T res = RES_OK;
    439   ASSERT(cmd);
    440 
    441   /* Nothing to do on non master process */
    442   if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit;
    443 
    444   htrdr_log(cmd->htrdr, "Write volumetric acceleration structure to '%s'.\n",
    445     str_cget(&cmd->output_name));
    446 
    447   res = atrstm_dump_svx_octree(cmd->medium, &args, cmd->output);
    448   if(res != RES_OK) goto error;
    449 
    450 exit:
    451   return res;
    452 error:
    453   goto exit;
    454 }
    455 
    456 static double
    457 compute_laser_mesh_extent(const struct htrdr_combustion* cmd)
    458 {
    459   double mdm_upp[3];
    460   double mdm_low[3];
    461   double laser_dir[3];
    462   double laser_pos[3];
    463   double t[2];
    464   int max_axis;
    465   ASSERT(cmd);
    466 
    467   /* Retrieve the medium axis aligned bounding box */
    468   atrstm_get_aabb(cmd->medium, mdm_low, mdm_upp);
    469 
    470   /* Retrieve laser parameters */
    471   htrdr_combustion_laser_get_position(cmd->laser, laser_pos);
    472   htrdr_combustion_laser_get_direction(cmd->laser, laser_dir);
    473 
    474   /* Compute the dominant axis of the laser direction */
    475   max_axis =
    476      fabs(laser_dir[0]) > fabs(laser_dir[1])
    477   ? (fabs(laser_dir[0]) > fabs(laser_dir[2]) ? 0 : 2)
    478   : (fabs(laser_dir[1]) > fabs(laser_dir[2]) ? 1 : 2);
    479 
    480   /* Define the intersection of the laser along its dominant axis with the
    481    * medium bounds along this axis */
    482   t[0] = (mdm_low[max_axis] - laser_pos[max_axis]) / laser_dir[max_axis];
    483   t[1] = (mdm_upp[max_axis] - laser_pos[max_axis]) / laser_dir[max_axis];
    484   if(t[0] > t[1]) SWAP(double, t[0], t[1]);
    485 
    486   /* Use the far intersection distance as the extent of the laser mesh */
    487   return t[1];
    488 }
    489 
    490 static res_T
    491 dump_laser_sheet(const struct htrdr_combustion* cmd)
    492 {
    493   struct htrdr_combustion_laser_mesh laser_mesh;
    494   double extent;
    495   unsigned i;
    496   res_T res = RES_OK;
    497   ASSERT(cmd);
    498 
    499   htrdr_log(cmd->htrdr, "Write laser sheet to '%s'.\n",
    500     str_cget(&cmd->output_name));
    501 
    502   /* Compute the extent of the geometry that will represent the laser sheet */
    503   extent = compute_laser_mesh_extent(cmd);
    504 
    505   /* Retreive the mesh of the laser sheet */
    506   htrdr_combustion_laser_get_mesh(cmd->laser, extent, &laser_mesh);
    507 
    508   #define FPRINTF(Fmt, Args) {                                                 \
    509     const int err = fprintf(cmd->output, Fmt COMMA_##Args LIST_##Args);        \
    510     if(err < 0) {                                                              \
    511       htrdr_log_err(cmd->htrdr, "Error writing data to `%s'.\n",               \
    512         str_cget(&cmd->output_name));                                          \
    513       res = RES_IO_ERR;                                                        \
    514       goto error;                                                              \
    515     }                                                                          \
    516   } (void)0
    517 
    518   /* Write header */
    519   FPRINTF("# vtk DataFile Version 2.0\n", ARG0());
    520   FPRINTF("Laser sheet\n", ARG0());
    521   FPRINTF("ASCII\n", ARG0());
    522   FPRINTF("DATASET POLYDATA\n", ARG0());
    523 
    524   /* Write the vertices */
    525   FPRINTF("POINTS %u double\n", ARG1(laser_mesh.nvertices));
    526   FOR_EACH(i, 0, laser_mesh.nvertices) {
    527     FPRINTF("%g %g %g\n", ARG3
    528       (laser_mesh.vertices[i*3+0],
    529        laser_mesh.vertices[i*3+1],
    530        laser_mesh.vertices[i*3+2]));
    531   }
    532 
    533   /* Write the triangles */
    534   FPRINTF("POLYGONS %u %u\n",ARG2
    535     (laser_mesh.ntriangles,
    536      laser_mesh.ntriangles*4));
    537   FOR_EACH(i, 0, laser_mesh.ntriangles) {
    538     FPRINTF("3 %u %u %u\n", ARG3
    539       (laser_mesh.triangles[i*3+0],
    540        laser_mesh.triangles[i*3+1],
    541        laser_mesh.triangles[i*3+2]));
    542   }
    543 
    544   /* Write flux density */
    545   FPRINTF("CELL_DATA %u\n", ARG1(laser_mesh.ntriangles));
    546   FPRINTF("SCALARS Flux_density double 1\n", ARG0());
    547   FPRINTF("LOOKUP_TABLE default\n", ARG0());
    548   FOR_EACH(i, 0, laser_mesh.ntriangles) {
    549     FPRINTF("%g\n", ARG1(htrdr_combustion_laser_get_flux_density(cmd->laser)));
    550   }
    551   #undef FPRINTF
    552 
    553 exit:
    554   return res;
    555 error:
    556   goto exit;
    557 }
    558 
    559 static void
    560 combustion_release(ref_T* ref)
    561 {
    562   struct htrdr_combustion* cmd = CONTAINER_OF(ref, struct htrdr_combustion, ref);
    563   struct htrdr* htrdr = NULL;
    564   ASSERT(ref);
    565 
    566   if(cmd->geom) htrdr_geometry_ref_put(cmd->geom);
    567   if(cmd->mats) htrdr_materials_ref_put(cmd->mats);
    568   if(cmd->medium) ATRSTM(ref_put(cmd->medium));
    569   if(cmd->camera) SCAM(ref_put(cmd->camera));
    570   if(cmd->flux_map) htrdr_rectangle_ref_put(cmd->flux_map);
    571   if(cmd->laser) htrdr_combustion_laser_ref_put(cmd->laser);
    572   if(cmd->buf) htrdr_buffer_ref_put(cmd->buf);
    573   if(cmd->output && cmd->output != stdout) CHK(fclose(cmd->output) == 0);
    574   release_phase_functions(cmd);
    575   str_release(&cmd->output_name);
    576 
    577   htrdr = cmd->htrdr;
    578   MEM_RM(htrdr_get_allocator(htrdr), cmd);
    579   htrdr_ref_put(htrdr);
    580 }
    581 
    582 /*******************************************************************************
    583  * Exported functions
    584  ******************************************************************************/
    585 res_T
    586 htrdr_combustion_create
    587   (struct htrdr* htrdr,
    588    const struct htrdr_combustion_args* args,
    589    struct htrdr_combustion** out_cmd)
    590 {
    591   struct htrdr_combustion* cmd = NULL;
    592   res_T res = RES_OK;
    593   ASSERT(htrdr && args && out_cmd);
    594 
    595   cmd = MEM_CALLOC(htrdr_get_allocator(htrdr), 1, sizeof(*cmd));
    596   if(!cmd) {
    597     htrdr_log_err(htrdr, "Could not allocate the htrdr_combustion data.\n");
    598     res = RES_MEM_ERR;
    599     goto error;
    600   }
    601   ref_init(&cmd->ref);
    602   str_init(htrdr_get_allocator(htrdr), &cmd->output_name);
    603 
    604   /* Get the ownership on the htrdr structure */
    605   htrdr_ref_get(htrdr);
    606   cmd->htrdr = htrdr;
    607 
    608   cmd->spp = args->image.spp;
    609   cmd->output_type = args->output_type;
    610   cmd->phase_func_type = args->phase_func_type;
    611 
    612   res = setup_output(cmd, args);
    613   if(res != RES_OK) goto error;
    614   res = setup_phase_functions(cmd);
    615   if(res != RES_OK) goto error;
    616   res = setup_simd(cmd, args);
    617   if(res != RES_OK) goto error;
    618   res = setup_geometry(cmd, args);
    619   if(res != RES_OK) goto error;
    620   res = setup_sensor(cmd, args);
    621   if(res != RES_OK) goto error;
    622   res = setup_laser(cmd, args);
    623   if(res != RES_OK) goto error;
    624   res = setup_buffer(cmd, args);
    625   if(res != RES_OK) goto error;
    626   res = setup_medium(cmd, args);
    627   if(res != RES_OK) goto error;
    628 
    629 exit:
    630   *out_cmd = cmd;
    631   return res;
    632 error:
    633   if(cmd) {
    634     htrdr_combustion_ref_put(cmd);
    635     cmd = NULL;
    636   }
    637   goto exit;
    638 }
    639 
    640 void
    641 htrdr_combustion_ref_get(struct htrdr_combustion* cmd)
    642 {
    643   ASSERT(cmd);
    644   ref_get(&cmd->ref);
    645 }
    646 
    647 void
    648 htrdr_combustion_ref_put(struct htrdr_combustion* cmd)
    649 {
    650   ASSERT(cmd);
    651   ref_put(&cmd->ref, combustion_release);
    652 }
    653 
    654 res_T
    655 htrdr_combustion_run(struct htrdr_combustion* cmd)
    656 {
    657   res_T res = RES_OK;
    658   ASSERT(cmd);
    659 
    660   switch(cmd->output_type) {
    661     case HTRDR_COMBUSTION_ARGS_OUTPUT_FLUX_MAP:
    662       res = combustion_draw_map(cmd);
    663       break;
    664     case HTRDR_COMBUSTION_ARGS_OUTPUT_IMAGE:
    665       res = combustion_draw_map(cmd);
    666       break;
    667     case HTRDR_COMBUSTION_ARGS_OUTPUT_LASER_SHEET:
    668       res = dump_laser_sheet(cmd);
    669       break;
    670     case HTRDR_COMBUSTION_ARGS_OUTPUT_OCTREES:
    671       res = dump_volumetric_acceleration_structure(cmd);
    672       break;
    673     default: FATAL("Unreachable code.\n"); break;
    674   }
    675   if(res != RES_OK) {
    676     goto error;
    677   }
    678 
    679 exit:
    680   return res;
    681 error:
    682   goto exit;
    683 }
    684 
    685 /*******************************************************************************
    686  * Local functions
    687  ******************************************************************************/
    688 void
    689 combustion_get_pixel_format
    690   (const struct htrdr_combustion* cmd,
    691    struct htrdr_pixel_format* fmt)
    692 {
    693   ASSERT(cmd && fmt);
    694   (void)cmd;
    695   switch(cmd->output_type) {
    696     case HTRDR_COMBUSTION_ARGS_OUTPUT_FLUX_MAP:
    697       fmt->size = sizeof(struct combustion_pixel_flux);
    698       fmt->alignment = ALIGNOF(struct combustion_pixel_flux);
    699       break;
    700     case HTRDR_COMBUSTION_ARGS_OUTPUT_IMAGE:
    701       fmt->size = sizeof(struct combustion_pixel_image);
    702       fmt->alignment = ALIGNOF(struct combustion_pixel_image);
    703       break;
    704     default: FATAL("Unreachable code.\n"); break;
    705   }
    706 }