commit 41046ab9f62215b2bf40b38dd6de0bd75c64d462 parent effa61f00d8558d9d98b77d26c7db8f421ef0755 Author: Vincent Forest <vincent.forest@meso-star.com> Date: Fri, 19 Feb 2021 16:57:05 +0100 Fix atmosphere compilation Diffstat:
19 files changed, 476 insertions(+), 374 deletions(-)
diff --git a/src/atmosphere/htrdr_atmosphere.c b/src/atmosphere/htrdr_atmosphere.c @@ -15,25 +15,27 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#include "htrdr_atmosphere.h" -#include "htrdr_atmosphere_args.h" -#include "htrdr_atmosphere_sun.h" - -#include "htrdr_buffer.h" -#include "htrdr_camera.h" -#include "htrdr_cie_xyz.h" -#include "htrdr_geometry.h" -#include "htrdr_log.h" -#include "htrdr_materials.h" - -#include <star/s3d.h> - -struct pixel_format { - size_t size; /* In bytes */ - size_t alignment; /* In bytes */ -}; -#define PIXEL_FORMAT_NULL__ {0, 0} -static const struct pixel_format PIXEL_FORMAT_NULL = PIXEL_FORMAT_NULL__; +#define _POSIX_C_SOURCE 200112L + +#include "atmosphere/htrdr_atmosphere.h" +#include "atmosphere/htrdr_atmosphere_c.h" +#include "atmosphere/htrdr_atmosphere_args.h" +#include "atmosphere/htrdr_atmosphere_ground.h" +#include "atmosphere/htrdr_atmosphere_sun.h" + +#include "core/htrdr_buffer.h" +#include "core/htrdr_camera.h" +#include "core/htrdr_cie_xyz.h" +#include "core/htrdr_log.h" +#include "core/htrdr_materials.h" +#include "core/htrdr_ran_wlen.h" +#include "core/htrdr_rectangle.h" + +#include <high_tune/htsky.h> + +#include <rsys/cstr.h> + +#include <math.h> /******************************************************************************* * Helper functions @@ -75,10 +77,10 @@ static size_t compute_spectral_bands_count(const struct htrdr_atmosphere* cmd) { double wlen_range[2]; - double wlen_size; - const size_t nbands; - const double band_len; - const double band_len_max; + double wlen_range_size; + size_t nbands; + double band_len; + double band_len_max; ASSERT(cmd); /* Compute size of the spectral range in nanometers */ @@ -101,7 +103,7 @@ compute_spectral_bands_count(const struct htrdr_atmosphere* cmd) if(band_len > band_len_max) { nbands = (size_t)ceil(wlen_range_size / band_len_max); } - return nbands + return nbands; } static enum htsky_spectral_type @@ -145,34 +147,6 @@ spherical_to_cartesian_dir dir[2] = sin_elevation; } -static INLINE void -get_pixel_format - (const struct htrdr_atmosphere* cmd, - struct pixel_format* fmt) -{ - switch(cmd->sensor_type) { - case HTRDR_SENSOR_RECTANGLE: - fmt->size = sizeof(struct atmosphere_pixel_flux); - fmt->alignment = ALIGNOF(struct atmosphere_pixel_flux); - break; - case HTRDR_SENSOR_CAMERA: - switch(cmd->spectral_type) { - case HTRDR_SPECTRAL_LW: - case HTRDR_SPECTRAL_SW: - fmt->size = sizeof(struct atmosphere_pixel_xwave); - fmt->alignment = ALIGNOF(struct atmosphere_pixel_xwave); - break; - case HTRDR_SPECTRAL_SW_CIE_XYZ: - fmt->size = sizeof(struct atmosphere_pixel_image); - fmt->alignment = ALIGNOF(struct atmosphere_pixel_image); - break; - default: FATAL("Unreachable code.\n"); break; - } - break; - default: FATAL("Unreachable code.\n"); break; - } -} - static res_T setup_sensor (struct htrdr_atmosphere* cmd, @@ -180,11 +154,11 @@ setup_sensor { double proj_ratio; res_T res = RES_OK; - ASSERT(htrdr && args); + ASSERT(cmd && args); cmd->sensor.type = args->sensor_type; - if(args->spectral.type == HTRDR_SPECTRAL_SW_CIE_XYZ + if(args->spectral.spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ && args->sensor_type != HTRDR_SENSOR_CAMERA) { htrdr_log_err(cmd->htrdr, "the CIE 1931 XYZ spectral integration can be used " "only with a camera sensor.\n"); @@ -197,14 +171,23 @@ setup_sensor proj_ratio = (double)args->image.definition[0] / (double)args->image.definition[1]; - res = htrdr_camera_create(htrdr, args->camera.position, - args->camera.target, args->camera.up, proj_ratio, - MDEG2RAD(args->camera.fov_y), &cmd->sensor.camera); + res = htrdr_camera_create + (cmd->htrdr, + args->sensor.camera.position, + args->sensor.camera.target, + args->sensor.camera.up, + proj_ratio, + MDEG2RAD(args->sensor.camera.fov_y), + &cmd->sensor.camera); break; case HTRDR_SENSOR_RECTANGLE: - res = htrdr_rectangle_create(htrdr, args->rectangle.size, - args->rectangle.position, args->rectangle.target, args->rectangle.up, - &cmd->sensor.rectangle); + res = htrdr_rectangle_create + (cmd->htrdr, + args->sensor.rectangle.size, + args->sensor.rectangle.position, + args->sensor.rectangle.target, + args->sensor.rectangle.up, + &cmd->sensor.rectangle); break; default: FATAL("Unreachable code.\n"); break; } @@ -227,7 +210,7 @@ dump_volumetric_acceleration_structure(struct htrdr_atmosphere* cmd) nbands = htsky_get_spectral_bands_count(cmd->sky); /* Nothing to do */ - if(htrdr->mpi_rank != 0) goto exit; + if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit; FOR_EACH(i, 0, nbands) { size_t iquad; @@ -236,136 +219,10 @@ dump_volumetric_acceleration_structure(struct htrdr_atmosphere* cmd) (cmd->sky, iband); FOR_EACH(iquad, 0, nquads) { - res = htsky_dump_cloud_vtk(htrdr->sky, iband, iquad, htrdr->output); + res = htsky_dump_cloud_vtk(cmd->sky, iband, iquad, cmd->output); if(res != RES_OK) goto error; - fprintf(htrdr->output, "---\n"); - } - } - -exit: - return res; -error: - goto exit; -} - -static INLINE void -dump_accum - (const struct htrdr_accum* acc, /* Accum to dump */ - struct htrdr_accum* out_acc, /* May be NULL */ - FILE* stream) -{ - ASSERT(acc && stream); - - if(acc->nweights == 0) { - fprintf(streamn, "0 0 "); - } else { - struct htrdr_estimate estimate = HTRDR_ESTIMATE_NULL; - - htrdr_accum_get_estimation(acc, &estimate); - fprintf("%g %g\n", estimate.E, estimate.SE); - - if(out_acc) { - out_acc->sum_weights += acc->sum_weights; - out_acc->sum_weights_sqr += acc->sum_weights_sqr; - out_acc->nweights += acc->nweights; - } - } -} - -static INLINE void -dump_pixel_flux - (const struct htrdr_pixel_flux* pix, - struct htrdr_accum* time_acc, /* May be NULL */ - struct htrdr_accum* flux_acc, /* May be NULL */ - FILE* stream) -{ - struct htrdr_estimate pix_time = HTRDR_ESTIMATE_NULL; - - ASSERT(pix && stream); - dump_accum(&pix->flux, flux_acc, stream_name, stream); - fprintf(stream, "0 0 0 0 "); - dump_accum(&pix->time, time_acc, stream_name, stream); - fprintf(stream, "\n"); -} - -static INLINE void -dump_pixel_image - (const struct htrdr_pixel_image* pix, - struct htrdr_accum* time_acc, /* May be NULL */ - FILE* stream) -{ - ASSERT(pix && stream_name && stream); - fprintf(stream, "%g %g ", pix->X.E, pix->X.SE); - fprintf(stream, "%g %g ", pix->Y.E, pix->Y.SE); - fprintf(stream, "%g %g ", pix->Z.E, pix->Z.SE); - dump_accum(pix->time, time_acc, stream); - fprintf(stream, "\n"); -} - -static INLINE void -dump_pixel_xwave - (const struct htrdr_pixel_xwave* pix, - struct htrdr_accum* time_acc, /* May be NULL */ - const char* stream_name, - FILE* stream) -{ - fprintf(stream, "%g %g %f %f 0 0 ", - pix->radiance_temperature.E, - pix->radiance_temperature.SE, - pix->radiance.E, - pix->radiance.SE); - dump_accum(pix->time, time_acc, stream); - fprintf(stream, "\n"); -} - -static res_T -dump_buffer - (struct htrdr_atmosphere* cmd, - struct htrdr_buffer* buf, - struct htrdr_accum* time_acc, /* May be NULL */ - struct htrdr_accum* flux_acc, /* May be NULL */ - const char* stream_name, - FILE* stream) -{ - struct htrdr_buffer_layout layout; - size_t pixsz, pixal; - size_t x, y; - res_T res = RES_OK; - ASSERT(cmd && buf && stream_name && stream); - (void)stream_name; - - pixsz = htrdr_spectral_type_get_pixsz(cmd->spectral_type, cmd->sensor.type); - pixal = htrdr_spectral_type_get_pixal(cmd->spectral_type, cmd->sensor.type); - - htrdr_buffer_get_layout(buf, &layout); - if(layout.elmt_size != pixsz || layout.alignment != pixal) { - htrdr_log_err(cmd->htrdr, "%s: invalid buffer layout. ", FUNC_NAME); - res = RES_BAD_ARG; - goto error; - } - - fprintf(stream, "%lu %lu\n", layout.width, layout.height); - - if(time_acc) *time_acc = HTRDR_ACCUM_NULL; - if(flux_acc) *flux_acc = HTRDR_ACCUM_NULL; - - FOR_EACH(y, 0, layout.height) { - FOR_EACH(x, 0, layout.width) { - struct htrdr_estimate pix_time = HTRDR_ESTIMATE_NULL; - const struct htrdr_accum* pix_time_acc = NULL; - - if(cmd->sensor.type == HTRDR_SENSOR_RECTANGLE) { - const struct htrdr_pixel_flux* pix = htrdr_buffer_at(buf, x, y); - dump_pixel_flux(pix, time_acc, flux_acc, stream_name, stream); - } else if(cmd->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ) { - const struct htrdr_pixel_image* pix = htrdr_buffer_at(buf, x, y); - dump_pixel_image(pix, time_acc, stream_name, stream); - } else { - const struct htrdr_pixel_xwave* pix = htrdr_buffer_at(buf, x, y); - dump_pixel_xwave(pix, time_acc, stream_name, stream); - } + fprintf(cmd->output, "---\n"); } - fprintf(stream, "\n"); } exit: @@ -381,7 +238,6 @@ atmosphere_release(ref_T* ref) struct htrdr* htrdr = NULL; ASSERT(ref); - if(cmd->s3d) S3D(device_ref_put(cmd->s3d)); if(cmd->ground) htrdr_atmosphere_ground_ref_put(cmd->ground); if(cmd->mats) htrdr_materials_ref_put(cmd->mats); if(cmd->sun) htrdr_atmosphere_sun_ref_put(cmd->sun); @@ -389,18 +245,18 @@ atmosphere_release(ref_T* ref) if(cmd->ran_wlen) htrdr_ran_wlen_ref_put(cmd->ran_wlen); if(cmd->sensor.camera) htrdr_camera_ref_put(cmd->sensor.camera); if(cmd->sensor.rectangle) htrdr_rectangle_ref_put(cmd->sensor.rectangle); - if(cmd->but) htrdr_buffer_ref_put(cmd->but); + if(cmd->buf) htrdr_buffer_ref_put(cmd->buf); if(cmd->sky) HTSKY(ref_put(cmd->sky)); if(cmd->output && cmd->output != stdout) fclose(cmd->output); - str_release(cmd->output_name); + str_release(&cmd->output_name); htrdr = cmd->htrdr; - MEM_RM(htrdr->allocator, cmd); + MEM_RM(htrdr_get_allocator(htrdr), cmd); htrdr_ref_put(htrdr); } /******************************************************************************* - * Local functions + * Exported functions ******************************************************************************/ res_T htrdr_atmosphere_create @@ -417,7 +273,7 @@ htrdr_atmosphere_create res_T res = RES_OK; ASSERT(htrdr && args && out_cmd); - cmd = MEM_CALLOC(htrdr->allocator, 1, sizeof(*cmd)); + cmd = MEM_CALLOC(htrdr_get_allocator(htrdr), 1, sizeof(*cmd)); if(!cmd) { htrdr_log_err(htrdr, "%s: could not allocate the htrdr_atmosphere data.\n", FUNC_NAME); @@ -425,29 +281,35 @@ htrdr_atmosphere_create goto error; } ref_init(&cmd->ref); - str_init(htrdr->allocator, &cmd->output_name); - cmd->dump_vtk = args->dump_vtk; + str_init(htrdr_get_allocator(htrdr), &cmd->output_name); + cmd->dump_volumetric_acceleration_structure = + args->dump_volumetric_acceleration_structure; cmd->verbose = args->verbose; cmd->spp = args->image.spp; cmd->width = args->image.definition[0]; cmd->height = args->image.definition[1]; cmd->grid_max_definition[0] = args->grid_max_definition[0]; + cmd->grid_max_definition[1] = args->grid_max_definition[1]; + cmd->grid_max_definition[2] = args->grid_max_definition[2]; + cmd->spectral_type = args->spectral.spectral_type; + cmd->ref_temperature = args->spectral.ref_temperature; + cmd->sky_mtl_name = args->sky_mtl_name; - /* Get ownership onf the htrdr structure */ + /* Get ownership on the htrdr structure */ htrdr_ref_get(htrdr); cmd->htrdr = htrdr; if(!args->filename_output) { cmd->output = stdout; output_name = "<stdout>"; - } else if(htrdr->mpi_rank != 0) { + } else if(htrdr_get_mpi_rank(htrdr) != 0) { cmd->output = NULL; output_name = "<null>"; } else { - res = htrdr_open_output_stream - (htrdr, args->output, 0/*read*/, args->force_overwriting, &cmd->output); + res = htrdr_open_output_stream(htrdr, args->filename_output, 0/*read*/, + args->force_overwriting, &cmd->output); if(res != RES_OK) goto error; - output_name = args->output; + output_name = args->filename_output; } res = str_set(&cmd->output_name, output_name); if(res != RES_OK) { @@ -457,32 +319,20 @@ htrdr_atmosphere_create goto error; } - /* Disable the Star-3D verbosity since the Embree backend prints some messages - * on stdout rather than stderr. This is annoying since stdout may be used by - * htrdr to write output data */ - res = s3d_device_create - (&htrdr->logger, htrdr->allocator, 0, &cmd->s3d); - if(res != RES_OK) { - htrdr_log_err(htrdr, - "%s: could not create the Star-3D device -- %s.\n", - FUNC_NAME, res_to_cstr(res)); - goto error; - } - /* Materials are necessary only if a ground geometry is defined */ if(args->filename_obj) { res = htrdr_materials_create(htrdr, args->filename_mtl, &cmd->mats); if(res != RES_OK) goto error; } - res = htrdr_atmosphere_ground_create(htrdr, args->filename_obj, + res = htrdr_atmosphere_ground_create(htrdr, args->filename_obj, cmd->mats, args->repeat_ground, &cmd->ground); if(res != RES_OK) goto error; res = setup_sensor(cmd, args); if(res != RES_OK) goto error; - res = htrdr_atmosphere_sun_create(htrdr, &cmd->sun); + res = htrdr_atmosphere_sun_create(cmd->htrdr, &cmd->sun); if(res != RES_OK) goto error; spherical_to_cartesian_dir (MDEG2RAD(args->sun_azimuth), MDEG2RAD(args->sun_elevation), sun_dir); @@ -496,23 +346,25 @@ htrdr_atmosphere_create htsky_args.grid_max_definition[1] = args->grid_max_definition[1]; htsky_args.grid_max_definition[2] = args->grid_max_definition[2]; htsky_args.optical_thickness = args->optical_thickness; - htsky_args.nthreads = htrdr->nthreads; + htsky_args.nthreads = (unsigned)htrdr_get_threads_count(htrdr); htsky_args.repeat_clouds = args->repeat_clouds; - htsky_args.verbose = htrdr->mpi_rank == 0 ? args->verbose : 0; - htsky_args.spectral_type = htrdr_to_sky_spectral_type(args->spectral_type); + htsky_args.verbose = htrdr_get_mpi_rank(htrdr) == 0 ? args->verbose : 0; + htsky_args.spectral_type = htrdr_to_sky_spectral_type(args->spectral.spectral_type); htsky_args.wlen_range[0] = args->spectral.wlen_range[0]; htsky_args.wlen_range[1] = args->spectral.wlen_range[1]; - res = htsky_create(&htrdr->logger, htrdr->allocator, &htsky_args, &cmd->sky); + res = htsky_create(htrdr_get_logger(htrdr), htrdr_get_allocator(htrdr), + &htsky_args, &cmd->sky); if(res != RES_OK) goto error; HTSKY(get_raw_spectral_bounds(cmd->sky, spectral_range)); - spectral_range[0] = MMAX(args->wlen_range[0], spectral_range[0]); - spectral_range[1] = MMIN(args->wlen_range[1], spectral_range[1]); - if(spectral_range[0] != args->wlen_range[0] - || spectral_range[1] != args->wlen_range[1]) { + spectral_range[0] = MMAX(args->spectral.wlen_range[0], spectral_range[0]); + spectral_range[1] = MMIN(args->spectral.wlen_range[1], spectral_range[1]); + if(spectral_range[0] != args->spectral.wlen_range[0] + || spectral_range[1] != args->spectral.wlen_range[1]) { htrdr_log_warn(htrdr, - "%s: the submitted spectral range overflowed the spectral data.\n", FUNC_NAME); + "%s: the submitted spectral range overflowed the spectral data.\n", + FUNC_NAME); } cmd->wlen_range_m[0] = spectral_range[0]*1e-9; /* Convert in meters */ @@ -539,9 +391,10 @@ htrdr_atmosphere_create /* Create the image buffer only on the master process; the image parts * rendered by the processes are gathered onto the master process. */ - if(!cmd->dump_volumetric_acceleration_structure && htrdr->mpi_rank == 0) { - struct pixel_format pixfmt = PIXEL_FORMAT_NULL; - get_pixel_format(cmd, &pixfmt); + if(!cmd->dump_volumetric_acceleration_structure + && htrdr_get_mpi_rank(htrdr) == 0) { + struct atmosphere_pixel_format pixfmt = ATMOSPHERE_PIXEL_FORMAT_NULL; + atmosphere_get_pixel_format(cmd, &pixfmt); res = htrdr_buffer_create(htrdr, args->image.definition[0], /* Width */ @@ -568,17 +421,17 @@ void htrdr_atmosphere_ref_get(struct htrdr_atmosphere* cmd) { ASSERT(cmd); - ref_get(cmd->ref); + ref_get(&cmd->ref); } void htrdr_atmosphere_ref_put(struct htrdr_atmosphere* cmd) { ASSERT(cmd); - ref_put(cmd->ref, atmosphere_release); + ref_put(&cmd->ref, atmosphere_release); } -void +res_T htrdr_atmosphere_run(struct htrdr_atmosphere* cmd) { res_T res = RES_OK; @@ -587,7 +440,7 @@ htrdr_atmosphere_run(struct htrdr_atmosphere* cmd) res = dump_volumetric_acceleration_structure(cmd); if(res != RES_OK) goto error; } else { - res = draw_map(cmd); + res = atmosphere_draw_map(cmd); if(res != RES_OK) goto error; } @@ -597,3 +450,34 @@ error: goto exit; } +/******************************************************************************* + * Local functions + ******************************************************************************/ +void +atmosphere_get_pixel_format + (const struct htrdr_atmosphere* cmd, + struct atmosphere_pixel_format* fmt) +{ + switch(cmd->sensor.type) { + case HTRDR_SENSOR_RECTANGLE: + fmt->size = sizeof(struct atmosphere_pixel_flux); + fmt->alignment = ALIGNOF(struct atmosphere_pixel_flux); + break; + case HTRDR_SENSOR_CAMERA: + switch(cmd->spectral_type) { + case HTRDR_SPECTRAL_LW: + case HTRDR_SPECTRAL_SW: + fmt->size = sizeof(struct atmosphere_pixel_xwave); + fmt->alignment = ALIGNOF(struct atmosphere_pixel_xwave); + break; + case HTRDR_SPECTRAL_SW_CIE_XYZ: + fmt->size = sizeof(struct atmosphere_pixel_image); + fmt->alignment = ALIGNOF(struct atmosphere_pixel_image); + break; + default: FATAL("Unreachable code.\n"); break; + } + break; + default: FATAL("Unreachable code.\n"); break; + } +} + diff --git a/src/atmosphere/htrdr_atmosphere_args.c b/src/atmosphere/htrdr_atmosphere_args.c @@ -15,13 +15,17 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#include "htrdr_atmosphere_args.h" +#include "atmosphere/htrdr_atmosphere_args.h" + +#include <rsys/cstr.h> + +#include <getopt.h> /******************************************************************************* * Helper functions ******************************************************************************/ static void -print_help_atmosphere(const char* cmd) +print_help(const char* cmd) { ASSERT(cmd); printf("Usage: %s [<opions>] -a GAS\n", cmd); @@ -46,7 +50,7 @@ print_help_atmosphere(const char* cmd) HTRDR_ATMOSPHERE_ARGS_DEFAULT.sun_elevation); printf( " -d dump volumetric acceleration structures to OUTPUT\n" -" and exit.\n");" +" and exit.\n"); printf( " -f overwrite the OUTPUT file if it already exists.\n"); printf( @@ -65,8 +69,8 @@ print_help_atmosphere(const char* cmd) " Its default value is `%s'.\n", HTRDR_ATMOSPHERE_ARGS_DEFAULT.sky_mtl_name); printf( -" -O CACHE filenaname of the cache file used to store/restore the - volumetric data. By default do not use any cache.\n"); +" -O CACHE filenaname of the cache file used to store/restore the\n" +" volumetric data. By default do not use any cache.\n"); printf( " -o OUTPUT file where data are written. If not defined, data are\n" " written to standard output.\n"); @@ -184,7 +188,6 @@ htrdr_atmosphere_args_init char** argv) { int opt; - int i; res_T res = RES_OK; ASSERT(args && argc && argv); @@ -204,7 +207,7 @@ htrdr_atmosphere_args_init case 'g': args->filename_obj = optarg; break; case 'h': print_help(argv[0]); - htrdr_args_release(args); + htrdr_atmosphere_args_release(args); args->quit = 1; goto exit; case 'i': @@ -213,8 +216,8 @@ htrdr_atmosphere_args_init case 'M': args->filename_mtl = optarg; break; case 'm': args->filename_mie = optarg; break; case 'n': args->sky_mtl_name = optarg; break; - case 'O': args->cache = optarg; break; - case 'o': args->output = optarg; break; + case 'O': args->filename_cache = optarg; break; + case 'o': args->filename_output = optarg; break; case 'p': args->sensor_type = HTRDR_SENSOR_RECTANGLE; res = htrdr_args_rectangle_parse(&args->sensor.rectangle, optarg); @@ -222,7 +225,7 @@ htrdr_atmosphere_args_init case 'r': args->repeat_clouds = 1; break; case 'R': args->repeat_ground = 1; break; case 's': - res = htrdr_args_spectral_parse(&args->sensor.spectral, optarg); + res = htrdr_args_spectral_parse(&args->spectral, optarg); break; case 'T': res = cstr_to_double(optarg, &args->optical_thickness); @@ -264,16 +267,16 @@ htrdr_atmosphere_args_init } /* Setup default ref temperature if necessary */ - if(args->ref_temperature <= 0) { - switch(args->spectral_type) { + if(args->spectral.ref_temperature <= 0) { + switch(args->spectral.spectral_type) { case HTRDR_SPECTRAL_LW: - args->ref_temperature = HTRDR_DEFAULT_LW_REF_TEMPERATURE; + args->spectral.ref_temperature = HTRDR_DEFAULT_LW_REF_TEMPERATURE; break; case HTRDR_SPECTRAL_SW: - args->ref_temperature = HTRDR_SUN_TEMPERATURE; + args->spectral.ref_temperature = HTRDR_SUN_TEMPERATURE; break; case HTRDR_SPECTRAL_SW_CIE_XYZ: - args->ref_temperature = -1; /* Unused */ + args->spectral.ref_temperature = -1; /* Unused */ break; default: FATAL("Unreachable code.\n"); break; } @@ -282,7 +285,14 @@ htrdr_atmosphere_args_init exit: return res; error: - htrdr_args_release(args); + htrdr_atmosphere_args_release(args); goto exit; } +void +htrdr_atmosphere_args_release(struct htrdr_atmosphere_args* args) +{ + ASSERT(args); + *args = HTRDR_ATMOSPHERE_ARGS_DEFAULT; +} + diff --git a/src/atmosphere/htrdr_atmosphere_args.h.in b/src/atmosphere/htrdr_atmosphere_args.h.in @@ -18,9 +18,13 @@ #ifndef HTRDR_ATMOSPHERE_ARGS_H #define HTRDR_ATMOSPHERE_ARGS_H -#include "htrdr_args.h" +#include "core/htrdr_args.h" +#include "core/htrdr_sensor.h" + #include <rsys/rsys.h> +#include <limits.h> /* UINT_MAX support */ + struct htrdr_atmosphere_args { const char* filename_gas; /* Path of the gas file */ const char* filename_les; /* Path of the HTCP file */ @@ -70,7 +74,7 @@ struct htrdr_atmosphere_args { NULL, /* Cache filename */ \ NULL, /* Output filename */ \ \ - @HTRDR_ARGS_DEFAULT_SKY_MTL_NAME@, /* Sky mtl name */ \ + @HTRDR_ATMOSPHERE_ARGS_DEFAULT_SKY_MTL_NAME@, /* Sky mtl name */ \ \ {HTRDR_ARGS_CAMERA_DEFAULT__}, /* Sensor */ \ HTRDR_SENSOR_CAMERA, /* Sensor type */ \ @@ -81,7 +85,7 @@ struct htrdr_atmosphere_args { 0, /* Sun azimuth */ \ 90, /* Sun elevation */ \ \ - @HTRDR_ARGS_ATMOPSHERE_DEFAULT_OPTICAL_THICKNESS_THRESHOLD@, \ + @HTRDR_ATMOSPHERE_ARGS_DEFAULT_OPTICAL_THICKNESS_THRESHOLD@, \ {UINT_MAX, UINT_MAX, UINT_MAX}, /* Maximum definition of the grid */ \ \ 0, /* Repeat clouds */ \ @@ -102,7 +106,7 @@ htrdr_atmosphere_args_init int argc, char** argv); -extern LOCAL_SYM res_T +extern LOCAL_SYM void htrdr_atmosphere_args_release (struct htrdr_atmosphere_args* args); diff --git a/src/atmosphere/htrdr_atmosphere_c.h b/src/atmosphere/htrdr_atmosphere_c.h @@ -18,11 +18,13 @@ #ifndef HTRDR_ATMOSPHERE_C_H #define HTRDR_ATMOSPHERE_C_H -#include "htrdr_accum.h" +#include "core/htrdr_accum.h" +#include "core/htrdr_sensor.h" +#include "core/htrdr_spectral.h" #include <rsys/ref_count.h> #include <rsys/rsys.h> -#include <rsys/string.h> +#include <rsys/str.h> /* Define the radiance component */ enum atmosphere_radiance_cpnt_flag { @@ -33,6 +35,14 @@ enum atmosphere_radiance_cpnt_flag { | ATMOSPHERE_RADIANCE_DIFFUSE }; +struct atmosphere_pixel_format { + size_t size; /* In bytes */ + size_t alignment; /* Power of two, in Bytes */ +}; +#define ATMOSPHERE_PIXEL_FORMAT_NULL__ {0, 0} +static const struct atmosphere_pixel_format ATMOSPHERE_PIXEL_FORMAT_NULL = + ATMOSPHERE_PIXEL_FORMAT_NULL__; + struct atmosphere_pixel_xwave { struct htrdr_estimate radiance; /* In W/m^2/sr */ struct htrdr_estimate radiance_temperature; /* In K */ @@ -79,11 +89,9 @@ struct htrdr_atmosphere_args; struct htrdr_buffer; struct htrdr_cie_xyz; struct htrdr_materials; -struct s3d_device; +struct htrdr_ran_wlen; struct htrdr_atmosphere { - struct s3d_device* s3d; - struct htrdr_atmosphere_ground* ground; struct htrdr_atmosphere_sun* sun; struct htrdr_materials* mats; @@ -115,12 +123,17 @@ struct htrdr_atmosphere { struct htrdr* htrdr; }; +extern LOCAL_SYM void +atmosphere_get_pixel_format + (const struct htrdr_atmosphere* cmd, + struct atmosphere_pixel_format* fmt); + extern LOCAL_SYM res_T atmosphere_draw_map (struct htrdr_atmosphere* cmd); /* Return the shortwave radiance in W/m^2/sr/m */ -extern LOCAL_SYM res_T +extern LOCAL_SYM double atmosphere_compute_radiance_sw (struct htrdr_atmosphere* cmd, const size_t ithread, @@ -133,7 +146,7 @@ atmosphere_compute_radiance_sw const size_t iquad); /* Return the longwave radiance in W/m^2/sr/m */ -extern LOCAL_SYM res_T +extern LOCAL_SYM double atmosphere_compute_radiance_lw (struct htrdr_atmosphere* cmd, const size_t ithread, diff --git a/src/atmosphere/htrdr_atmosphere_compute_radiance_lw.c b/src/atmosphere/htrdr_atmosphere_compute_radiance_lw.c @@ -15,11 +15,11 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#include "htrdr.h" -#include "htrdr_c.h" -#include "htrdr_interface.h" -#include "htrdr_ground.h" -#include "htrdr_solve.h" +#include "atmosphere/htrdr_atmosphere_c.h" +#include "atmosphere/htrdr_atmosphere_ground.h" + +#include "core/htrdr.h" +#include "core/htrdr_interface.h" #include <high_tune/htsky.h> @@ -132,8 +132,8 @@ hit_filter * Local functions ******************************************************************************/ double -htrdr_compute_radiance_lw - (struct htrdr* htrdr, +atmosphere_compute_radiance_lw + (struct htrdr_atmosphere* cmd, const size_t ithread, struct ssp_rng* rng, const double pos_in[3], @@ -158,13 +158,16 @@ htrdr_compute_radiance_lw double g; double w = 0; /* Weight */ - ASSERT(htrdr && rng && pos_in && dir_in && ithread < htrdr->nthreads); + ASSERT(cmd && rng && pos_in && dir_in); + ASSERT(ithread < htrdr_get_threads_count(cmd->htrdr)); /* Setup the phase function for this spectral band & quadrature point */ CHK(RES_OK == ssf_phase_create - (&htrdr->lifo_allocators[ithread], &ssf_phase_hg, &phase_hg)); + (htrdr_get_thread_allocator(cmd->htrdr, ithread), + &ssf_phase_hg, + &phase_hg)); g = htsky_fetch_per_wavelength_particle_phase_function_asymmetry_parameter - (htrdr->sky, wlen); + (cmd->sky, wlen); SSF(phase_hg_setup(phase_hg, g)); /* Initialise the random walk */ @@ -180,20 +183,20 @@ htrdr_compute_radiance_lw /* Setup the remaining fields of the hit filter context */ ctx.rng = rng; - ctx.htsky = htrdr->sky; + ctx.htsky = cmd->sky; ctx.iband = iband; ctx.iquad = iquad; /* Found the first intersection with the surface geometry */ - HTRDR(ground_trace_ray - (htrdr->ground, pos, dir, range, &s3d_hit_prev, &s3d_hit)); + HTRDR(atmosphere_ground_trace_ray + (cmd->ground, pos, dir, range, &s3d_hit_prev, &s3d_hit)); /* Fit the ray range to the surface distance along the ray */ range[0] = 0; range[1] = s3d_hit.distance; /* Trace a ray into the participating media */ - HTSKY(trace_ray(htrdr->sky, pos, dir, range, NULL, + HTSKY(trace_ray(cmd->sky, pos, dir, range, NULL, hit_filter, &ctx, iband, iquad, &svx_hit)); /* No scattering and no surface reflection. @@ -211,9 +214,9 @@ htrdr_compute_radiance_lw /* Absorption event. Stop the realisation */ if(ctx.event_type == EVENT_ABSORPTION) { ASSERT(!SVX_HIT_NONE(&svx_hit)); - temperature = htsky_fetch_temperature(htrdr->sky, pos_next); + temperature = htsky_fetch_temperature(cmd->sky, pos_next); /* weight is planck integrated over the spectral sub-interval */ - w = planck_monochromatic(wlen_m, temperature); + w = htrdr_planck_monochromatic(wlen_m, temperature); break; } @@ -238,9 +241,9 @@ htrdr_compute_radiance_lw ASSERT(!S3D_HIT_NONE(&s3d_hit)); /* Fetch the hit interface materal and build its BSDF */ - htrdr_ground_get_interface(htrdr->ground, &s3d_hit, &interf); + htrdr_atmosphere_ground_get_interface(cmd->ground, &s3d_hit, &interf); mtl = htrdr_interface_fetch_hit_mtl(&interf, dir, &s3d_hit); - HTRDR(mtl_create_bsdf(htrdr, mtl, ithread, wlen, rng, &bsdf)); + HTRDR(mtl_create_bsdf(cmd->htrdr, mtl, ithread, wlen, rng, &bsdf)); d3_normalize(N, d3_set_f3(N, s3d_hit.normal)); if(d3_dot(N, wo) < 0) d3_minus(N, N); @@ -257,7 +260,7 @@ htrdr_compute_radiance_lw if(ssp_rng_canonical(rng) >= bounce_reflectivity) { /* Absorbed at boundary */ temperature = mtl->temperature; /* Fetch mtl temperature */ /* Weight is planck integrated over the spectral sub-interval */ - w = temperature > 0 ? planck_monochromatic(wlen_m, temperature) : 0; + w = temperature>0 ? htrdr_planck_monochromatic(wlen_m, temperature) : 0; break; } s3d_hit_prev = s3d_hit; diff --git a/src/atmosphere/htrdr_atmosphere_compute_radiance_sw.c b/src/atmosphere/htrdr_atmosphere_compute_radiance_sw.c @@ -15,12 +15,11 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#include "htrdr.h" -#include "htrdr_c.h" -#include "htrdr_interface.h" -#include "htrdr_ground.h" -#include "htrdr_solve.h" -#include "htrdr_sun.h" +#include "atmosphere/htrdr_atmosphere_c.h" +#include "atmosphere/htrdr_atmosphere_ground.h" +#include "atmosphere/htrdr_atmosphere_sun.h" + +#include "core/htrdr_interface.h" #include <high_tune/htsky.h> @@ -213,7 +212,7 @@ transmissivity_hit_filter static double transmissivity - (struct htrdr* htrdr, + (struct htrdr_atmosphere* cmd, struct ssp_rng* rng, const enum htsky_property prop, const size_t iband, @@ -225,17 +224,17 @@ transmissivity struct svx_hit svx_hit; struct transmissivity_context transmissivity_ctx = TRANSMISSION_CONTEXT_NULL; - ASSERT(htrdr && rng && pos && dir && range); + ASSERT(cmd && rng && pos && dir && range); transmissivity_ctx.rng = rng; - transmissivity_ctx.sky = htrdr->sky; + transmissivity_ctx.sky = cmd->sky; transmissivity_ctx.iband = iband; transmissivity_ctx.iquad = iquad; transmissivity_ctx.Ts = ssp_ran_exp(rng, 1); /* Sample an optical thickness */ transmissivity_ctx.prop = prop; /* Compute the transmissivity */ - HTSKY(trace_ray(htrdr->sky, pos, dir, range, NULL, + HTSKY(trace_ray(cmd->sky, pos, dir, range, NULL, transmissivity_hit_filter, &transmissivity_ctx, iband, iquad, &svx_hit)); if(SVX_HIT_NONE(&svx_hit)) { @@ -292,13 +291,13 @@ atmosphere_compute_radiance_sw /* Create the Henyey-Greenstein phase function */ CHK(RES_OK == ssf_phase_create - (htrdr_get_lifo_allocator(cmd->htrdr, ithread), + (htrdr_get_thread_allocator(cmd->htrdr, ithread), &ssf_phase_hg, &phase_hg)); /* Create the Rayleigh phase function */ CHK(RES_OK == ssf_phase_create - (htrdr_get_lifo_allocator(cmd->htrdr, ithread), + (htrdr_get_thread_allocator(cmd->htrdr, ithread), &ssf_phase_rayleigh, &phase_rayleigh)); @@ -316,11 +315,12 @@ atmosphere_compute_radiance_sw d3_set(pos, pos_in); d3_set(dir, dir_in); - if((cpnt_mask & HTRDR_RADIANCE_DIRECT) /* Handle direct contribation */ + if((cpnt_mask & ATMOSPHERE_RADIANCE_DIRECT) /* Handle direct contribation */ && htrdr_atmosphere_sun_is_dir_in_solar_cone(cmd->sun, dir)) { /* Check that the ray is not occluded along the submitted range */ d2(range, 0, FLT_MAX); - HTRDR(ground_trace_ray(cmd->ground, pos, dir, range, NULL, &s3d_hit_tmp)); + HTRDR(atmosphere_ground_trace_ray + (cmd->ground, pos, dir, range, NULL, &s3d_hit_tmp)); if(!S3D_HIT_NONE(&s3d_hit_tmp)) { Tr = 0; } else { @@ -330,7 +330,7 @@ atmosphere_compute_radiance_sw } } - if((cpnt_mask & HTRDR_RADIANCE_DIFFUSE) == 0) + if((cpnt_mask & ATMOSPHERE_RADIANCE_DIFFUSE) == 0) goto exit; /* Discard diffuse contribution */ /* Radiative random walk */ @@ -346,7 +346,7 @@ atmosphere_compute_radiance_sw /* Find the first intersection with a surface */ d2(range, 0, DBL_MAX); - HTRDR(ground_trace_ray + HTRDR(atmosphere_ground_trace_ray (cmd->ground, pos, dir, range, &s3d_hit_prev, &s3d_hit)); /* Sample an optical thickness */ @@ -398,9 +398,9 @@ atmosphere_compute_radiance_sw const struct htrdr_mtl* mtl = NULL; /* Fetch the hit interface materal and build its BSDF */ - htrdr_ground_get_interface(cmd->ground, &s3d_hit, &interf); + htrdr_atmosphere_ground_get_interface(cmd->ground, &s3d_hit, &interf); mtl = htrdr_interface_fetch_hit_mtl(&interf, dir, &s3d_hit); - HTRDR(mtl_create_bsdf(cmd, mtl, ithread, wlen, rng, &bsdf)); + HTRDR(mtl_create_bsdf(cmd->htrdr, mtl, ithread, wlen, rng, &bsdf)); /* Revert the normal if necessary to match the SSF convention */ d3_normalize(N, d3_set_f3(N, s3d_hit.normal)); @@ -449,7 +449,8 @@ atmosphere_compute_radiance_sw sun_dir_pdf = 1.0; } else { /* Sample a sun direction */ - sun_dir_pdf = htrdr_sun_sample_direction(cmd->sun, rng, sun_dir); + sun_dir_pdf = htrdr_atmosphere_sun_sample_direction + (cmd->sun, rng, sun_dir); if(surface_scattering) { R = d3_dot(N, sun_dir) < 0/* Below the ground */ ? 0 : ssf_bsdf_eval(bsdf, wo, N, sun_dir) * d3_dot(N, sun_dir); @@ -465,7 +466,7 @@ atmosphere_compute_radiance_sw } else { /* Check that the sun is visible from the new position */ d2(range, 0, FLT_MAX); - HTRDR(ground_trace_ray + HTRDR(atmosphere_ground_trace_ray (cmd->ground, pos_next, sun_dir, range, &s3d_hit_prev, &s3d_hit_tmp)); /* Compute the sun transmissivity */ diff --git a/src/atmosphere/htrdr_atmosphere_draw_map.c b/src/atmosphere/htrdr_atmosphere_draw_map.c @@ -15,7 +15,29 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#include "htrdr_draw_map.h" +#include "atmosphere/htrdr_atmosphere_c.h" +#include "atmosphere/htrdr_atmosphere_ground.h" +#include "atmosphere/htrdr_atmosphere_sun.h" + +#include "core/htrdr.h" +#include "core/htrdr_buffer.h" +#include "core/htrdr_camera.h" +#include "core/htrdr_cie_xyz.h" +#include "core/htrdr_draw_map.h" +#include "core/htrdr_interface.h" +#include "core/htrdr_log.h" +#include "core/htrdr_ran_wlen.h" +#include "core/htrdr_rectangle.h" + +#include <high_tune/htsky.h> + +#include <star/s3d.h> +#include <star/ssp.h> + +#include <rsys/clock_time.h> +#include <rsys/str.h> + +#include <string.h> /******************************************************************************* * Helper functions @@ -45,7 +67,8 @@ sample_rectangle_ray htrdr_rectangle_sample_pos(rect, pix_samp, ray_org); /* Check that `ray_org' is not included into a geometry */ - HTRDR(ground_trace_ray(cmd->ground, ray_org, up_dir, range, NULL, &hit)); + HTRDR(atmosphere_ground_trace_ray + (cmd->ground, ray_org, up_dir, range, NULL, &hit)); /* Up direction is occluded. Check if the sample must be rejected, i.e. does it * lies inside a geometry? */ @@ -63,7 +86,7 @@ sample_rectangle_ray /* Fetch the hit interface and retrieve the material into which the ray was * traced */ - htrdr_ground_get_interface(cmd->ground, &hit, &interf); + htrdr_atmosphere_ground_get_interface(cmd->ground, &hit, &interf); mtl = cos_wi_N < 0 ? &interf.mtl_front : &interf.mtl_back; /* Reject the sample if the incident direction do not travel into the sky */ @@ -80,15 +103,16 @@ sample_rectangle_ray static void draw_pixel_image (struct htrdr* htrdr, - struct htrdr_draw_pixel_args* args, + const struct htrdr_draw_pixel_args* args, void* data) { struct htrdr_accum XYZ[3]; /* X, Y, and Z */ struct htrdr_accum time; struct htrdr_atmosphere* cmd; - struct htrdr_pixel_image* pixel = data; + struct atmosphere_pixel_image* pixel = data; size_t ichannel; ASSERT(htrdr && htrdr_draw_pixel_args_check(args) && data); + (void)htrdr; cmd = args->context; ASSERT(cmd); @@ -147,7 +171,7 @@ draw_pixel_image /* Compute the radiance in W/m^2/sr/m */ weight = atmosphere_compute_radiance_sw(cmd, args->ithread, args->rng, - HTRDR_RADIANCE_ALL, ray_org, ray_dir, wlen, iband, iquad); + ATMOSPHERE_RADIANCE_ALL, ray_org, ray_dir, wlen, iband, iquad); ASSERT(weight >= 0); pdf *= 1.e9; /* Transform the pdf from nm^-1 to m^-1 */ @@ -185,9 +209,10 @@ draw_pixel_flux struct htrdr_accum flux; struct htrdr_accum time; struct htrdr_atmosphere* cmd; - struct htrdr_pixel_flux* pixel = data; + struct atmosphere_pixel_flux* pixel = data; size_t isamp; ASSERT(htrdr && htrdr_draw_pixel_args_check(args) && data); + (void)htrdr; cmd = args->context; ASSERT(cmd); @@ -199,7 +224,7 @@ draw_pixel_flux flux = HTRDR_ACCUM_NULL; time = HTRDR_ACCUM_NULL; - FOR_EACH(isamp, 0, spp) { + FOR_EACH(isamp, 0, args->spp) { struct time t0, t1; double ray_org[3]; double ray_dir[3]; @@ -215,13 +240,13 @@ draw_pixel_flux /* Begin the registration of the time spent in the realisation */ time_current(&t0); - res = sample_rectangle_ray(cmd, htrdr->sensor.rectangle, args->pixel_coord, + res = sample_rectangle_ray(cmd, cmd->sensor.rectangle, args->pixel_coord, args->pixel_normalized_size, args->rng, ray_org, ray_dir); if(res != RES_OK) continue; /* Reject the current sample */ - r0 = ssp_rng_canonical(rng); - r1 = ssp_rng_canonical(rng); - r2 = ssp_rng_canonical(rng); + r0 = ssp_rng_canonical(args->rng); + r1 = ssp_rng_canonical(args->rng); + r2 = ssp_rng_canonical(args->rng); /* Sample a wavelength */ wlen = htrdr_ran_wlen_sample(cmd->ran_wlen, r0, r1, &band_pdf); @@ -244,7 +269,7 @@ draw_pixel_flux ASSERT(cmd->spectral_type == HTRDR_SPECTRAL_SW); /* Compute direct contribution if necessary */ - htrdr_sun_sample_direction(cmd->sun, rng, sun_dir); + htrdr_atmosphere_sun_sample_direction(cmd->sun, args->rng, sun_dir); htrdr_rectangle_get_normal(cmd->sensor.rectangle, N); cos_N_sun_dir = d3_dot(N, sun_dir); @@ -252,15 +277,15 @@ draw_pixel_flux L_direct = 0; } else { L_direct = atmosphere_compute_radiance_sw(cmd, args->ithread, - args->rng, HTRDR_RADIANCE_DIRECT, ray_org, sun_dir, wlen, iband, + args->rng, ATMOSPHERE_RADIANCE_DIRECT, ray_org, sun_dir, wlen, iband, iquad); } /* Compute diffuse contribution */ L_diffuse = atmosphere_compute_radiance_sw(cmd, args->ithread, args->rng, - HTRDR_RADIANCE_DIFFUSE, ray_org, ray_dir, wlen, iband, iquad); + ATMOSPHERE_RADIANCE_DIFFUSE, ray_org, ray_dir, wlen, iband, iquad); - sun_solid_angle = htrdr_sun_get_solid_angle(cmd->sun); + sun_solid_angle = htrdr_atmosphere_sun_get_solid_angle(cmd->sun); /* Compute the weight in W/m^2/m */ weight = cos_N_sun_dir * sun_solid_angle * L_direct + PI * L_diffuse; @@ -297,10 +322,12 @@ draw_pixel_xwave { struct htrdr_accum radiance; struct htrdr_accum time; - struct htrdr_pixel_xwave* pixel = data; + struct htrdr_atmosphere* cmd; + struct atmosphere_pixel_xwave* pixel = data; size_t isamp; double temp_min, temp_max; ASSERT(htrdr && htrdr_draw_pixel_args_check(args) && data); + (void)htrdr; cmd = args->context; ASSERT(cmd->sensor.type == HTRDR_SENSOR_CAMERA); @@ -311,7 +338,7 @@ draw_pixel_xwave radiance = HTRDR_ACCUM_NULL; time = HTRDR_ACCUM_NULL; - FOR_EACH(isamp, 0, spp) { + FOR_EACH(isamp, 0, args->spp) { struct time t0, t1; double pix_samp[2]; double ray_org[3]; @@ -323,7 +350,6 @@ draw_pixel_xwave size_t iquad; double usec; double band_pdf; - res_T res = RES_OK; /* Begin the registration of the time spent in the realisation */ time_current(&t0); @@ -338,9 +364,9 @@ draw_pixel_xwave * pixel sample */ htrdr_camera_ray(cmd->sensor.camera, pix_samp, ray_org, ray_dir); - r0 = ssp_rng_canonical(rng); - r1 = ssp_rng_canonical(rng); - r2 = ssp_rng_canonical(rng); + r0 = ssp_rng_canonical(args->rng); + r1 = ssp_rng_canonical(args->rng); + r2 = ssp_rng_canonical(args->rng); /* Sample a wavelength */ wlen = htrdr_ran_wlen_sample(cmd->ran_wlen, r0, r1, &band_pdf); @@ -350,14 +376,14 @@ draw_pixel_xwave iquad = htsky_spectral_band_sample_quadrature(cmd->sky, r2, iband); /* Compute the spectral radiance in W/m^2/sr/m */ - switch(htrdr->spectral_type) { + switch(cmd->spectral_type) { case HTRDR_SPECTRAL_LW: weight = atmosphere_compute_radiance_lw(cmd, args->ithread, args->rng, ray_org, ray_dir, wlen, iband, iquad); break; case HTRDR_SPECTRAL_SW: weight = atmosphere_compute_radiance_sw(cmd, args->ithread, args->rng, - HTRDR_RADIANCE_ALL, ray_org, ray_dir, wlen, iband, iquad); + ATMOSPHERE_RADIANCE_ALL, ray_org, ray_dir, wlen, iband, iquad); break; default: FATAL("Unreachable code.\n"); break; } @@ -388,7 +414,7 @@ draw_pixel_xwave /* Compute the brightness_temperature of the pixel and estimate its standard * error if the sources were in the medium (<=> longwave) */ - if(htrdr->spectral_type == HTRDR_SPECTRAL_LW) { + if(cmd->spectral_type == HTRDR_SPECTRAL_LW) { const double wlen_min = cmd->wlen_range_m[0]; const double wlen_max = cmd->wlen_range_m[1]; pixel->radiance_temperature.E = htrdr_radiance_temperature @@ -403,7 +429,7 @@ draw_pixel_xwave static INLINE void setup_draw_map_args_rectangle - (const struct htrdr_atmosphere* cmd, + (struct htrdr_atmosphere* cmd, struct htrdr_draw_map_args* args) { ASSERT(cmd && cmd->sensor.type == HTRDR_SENSOR_RECTANGLE && args); @@ -415,7 +441,7 @@ setup_draw_map_args_rectangle static INLINE void setup_draw_map_args_camera - (const struct htrdr_atmosphere* cmd, + (struct htrdr_atmosphere* cmd, struct htrdr_draw_map_args* args) { ASSERT(cmd && cmd->sensor.type == HTRDR_SENSOR_CAMERA && args); @@ -436,6 +462,118 @@ setup_draw_map_args_camera } } +static INLINE void +dump_accum + (const struct htrdr_accum* acc, /* Accum to dump */ + struct htrdr_accum* out_acc, /* May be NULL */ + FILE* stream) +{ + ASSERT(acc && stream); + + if(acc->nweights == 0) { + fprintf(stream, "0 0 "); + } else { + struct htrdr_estimate estimate = HTRDR_ESTIMATE_NULL; + + htrdr_accum_get_estimation(acc, &estimate); + fprintf(stream, "%g %g ", estimate.E, estimate.SE); + + if(out_acc) { + out_acc->sum_weights += acc->sum_weights; + out_acc->sum_weights_sqr += acc->sum_weights_sqr; + out_acc->nweights += acc->nweights; + } + } +} + +static INLINE void +dump_pixel_flux + (const struct atmosphere_pixel_flux* pix, + struct htrdr_accum* time_acc, /* May be NULL */ + struct htrdr_accum* flux_acc, /* May be NULL */ + FILE* stream) +{ + ASSERT(pix && stream); + dump_accum(&pix->flux, flux_acc, stream); + fprintf(stream, "0 0 0 0 "); + dump_accum(&pix->time, time_acc, stream); + fprintf(stream, "\n"); +} + +static INLINE void +dump_pixel_image + (const struct atmosphere_pixel_image* pix, + struct htrdr_accum* time_acc, /* May be NULL */ + FILE* stream) +{ + ASSERT(pix && stream); + fprintf(stream, "%g %g ", pix->X.E, pix->X.SE); + fprintf(stream, "%g %g ", pix->Y.E, pix->Y.SE); + fprintf(stream, "%g %g ", pix->Z.E, pix->Z.SE); + dump_accum(&pix->time, time_acc, stream); + fprintf(stream, "\n"); +} + +static INLINE void +dump_pixel_xwave + (const struct atmosphere_pixel_xwave* pix, + struct htrdr_accum* time_acc, /* May be NULL */ + FILE* stream) +{ + ASSERT(pix && stream); + fprintf(stream, "%g %g %f %f 0 0 ", + pix->radiance_temperature.E, + pix->radiance_temperature.SE, + pix->radiance.E, + pix->radiance.SE); + dump_accum(&pix->time, time_acc, stream); + fprintf(stream, "\n"); +} + +static res_T +dump_buffer + (struct htrdr_atmosphere* cmd, + struct htrdr_buffer* buf, + struct htrdr_accum* time_acc, /* May be NULL */ + struct htrdr_accum* flux_acc, /* May be NULL */ + FILE* stream) +{ + struct atmosphere_pixel_format pixfmt; + struct htrdr_buffer_layout layout; + size_t x, y; + ASSERT(cmd && buf && stream); + + atmosphere_get_pixel_format(cmd, &pixfmt); + htrdr_buffer_get_layout(buf, &layout); + CHK(pixfmt.size == layout.elmt_size); + + fprintf(stream, "%lu %lu\n", layout.width, layout.height); + + if(time_acc) *time_acc = HTRDR_ACCUM_NULL; + if(flux_acc) *flux_acc = HTRDR_ACCUM_NULL; + + FOR_EACH(y, 0, layout.height) { + FOR_EACH(x, 0, layout.width) { + void* pix_raw = htrdr_buffer_at(buf, x, y); + ASSERT(IS_ALIGNED(pix_raw, pixfmt.alignment)); + + if(cmd->sensor.type == HTRDR_SENSOR_RECTANGLE) { + const struct atmosphere_pixel_flux* pix = pix_raw; + dump_pixel_flux(pix, time_acc, flux_acc, stream); + } else if(cmd->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ) { + const struct atmosphere_pixel_image* pix = pix_raw; + dump_pixel_image(pix, time_acc, stream); + } else { + const struct atmosphere_pixel_xwave* pix = pix_raw; + dump_pixel_xwave(pix, time_acc, stream); + } + } + fprintf(stream, "\n"); + } + + return RES_OK; +} + /******************************************************************************* * Local functions ******************************************************************************/ @@ -447,6 +585,7 @@ atmosphere_draw_map(struct htrdr_atmosphere* cmd) struct htrdr_accum flux_acc = HTRDR_ACCUM_NULL; struct htrdr_estimate path_time; struct htrdr_estimate flux; + res_T res = RES_OK; ASSERT(cmd); args.spp = cmd->spp; @@ -465,11 +604,10 @@ atmosphere_draw_map(struct htrdr_atmosphere* cmd) if(res != RES_OK) goto error; /* No more to do on non master processes */ - if(htrdr_mpi_rank(cmd->htrdr) != 0) goto exit; + if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit; /* Write buffer to output */ - res = dump_buffer(cmd->htrdr, cmd->buf, &path_time_acc, &flux_acc, - str_cget(&cmd->output_name), cmd->output); + res = dump_buffer(cmd, cmd->buf, &path_time_acc, &flux_acc, cmd->output); if(res != RES_OK) goto error; htrdr_accum_get_estimation(&path_time_acc, &path_time); @@ -477,7 +615,7 @@ atmosphere_draw_map(struct htrdr_atmosphere* cmd) "Time per radiative path (in micro seconds): %g +/- %g\n", path_time.E, path_time.SE); - if(htrdr->sensor.type == HTRDR_SENSOR_RECTANGLE) { + if(cmd->sensor.type == HTRDR_SENSOR_RECTANGLE) { htrdr_accum_get_estimation(&flux_acc, &flux); htrdr_log(cmd->htrdr, "Radiative flux density (in W/(external m^2)): %g +/- %g\n", diff --git a/src/atmosphere/htrdr_atmosphere_ground.c b/src/atmosphere/htrdr_atmosphere_ground.c @@ -17,13 +17,17 @@ #define _POSIX_C_SOURCE 200112L /* strtok_r support */ -#include "htrdr.h" -#include "htrdr_interface.h" -#include "htrdr_atmosphere_ground.h" -#include "htrdr_materials.h" -#include "htrdr_slab.h" +#include "atmosphere/htrdr_atmosphere_ground.h" + +#include "core/htrdr.h" +#include "core/htrdr_interface.h" +#include "core/htrdr_log.h" +#include "core/htrdr_materials.h" +#include "core/htrdr_slab.h" #include <aw.h> +#include <star/s3d.h> + #include <rsys/clock_time.h> #include <rsys/cstr.h> #include <rsys/dynamic_array_double.h> @@ -33,8 +37,7 @@ #include <rsys/float2.h> #include <rsys/float3.h> #include <rsys/hash_table.h> - -#include <star/s3d.h> +#include <rsys/ref_count.h> #include <string.h> /* strtok_r */ @@ -162,6 +165,7 @@ trace_ground static res_T parse_shape_interface (struct htrdr* htrdr, + struct htrdr_materials* mats, const char* name, struct htrdr_interface* interf) { @@ -177,9 +181,9 @@ parse_shape_interface int has_thin = 0; int has_back = 0; res_T res = RES_OK; - ASSERT(htrdr && name && interf); + ASSERT(htrdr && mats && name && interf); - str_init(htrdr->allocator, &str); + str_init(htrdr_get_allocator(htrdr), &str); /* Locally copy the string to parse */ res = str_set(&str, name); @@ -216,8 +220,7 @@ parse_shape_interface /* Fetch the interface material if any */ if(mtl_name_thin) { - has_thin = htrdr_materials_find_mtl - (htrdr->mats, mtl_name_thin, &interf->mtl_thin); + has_thin = htrdr_materials_find_mtl(mats, mtl_name_thin, &interf->mtl_thin); if(!has_thin) { htrdr_log_err(htrdr, "Invalid interface `%s'. The interface material `%s' is unknown.\n", @@ -228,8 +231,7 @@ parse_shape_interface } /* Fetch the front material */ - has_front = htrdr_materials_find_mtl - (htrdr->mats, mtl_name_front, &interf->mtl_front); + has_front = htrdr_materials_find_mtl(mats, mtl_name_front, &interf->mtl_front); if(!has_front) { htrdr_log_err(htrdr, "Invalid interface `%s'. The front material `%s' is unknown.\n", @@ -239,8 +241,7 @@ parse_shape_interface } /* Fetch the back material */ - has_back = htrdr_materials_find_mtl - (htrdr->mats, mtl_name_back, &interf->mtl_back); + has_back = htrdr_materials_find_mtl(mats, mtl_name_back, &interf->mtl_back); if(!has_back) { htrdr_log_err(htrdr, "Invalid interface `%s'. The back material `%s' is unknown.\n", @@ -256,7 +257,6 @@ error: *interf = HTRDR_INTERFACE_NULL; goto exit; } - static res_T setup_mesh (struct htrdr* htrdr, @@ -389,7 +389,7 @@ create_s3d_shape mesh.positions = positions; mesh.indices = indices; - res = s3d_shape_create_mesh(htrdr->s3d, &shape); + res = s3d_shape_create_mesh(htrdr_get_s3d(htrdr), &shape); if(res != RES_OK) { htrdr_log_err(htrdr, "Error creating a Star-3D shape -- %s.\n", res_to_cstr(res)); @@ -429,7 +429,10 @@ error: } static res_T -setup_ground(struct htrdr_atmosphere_ground* ground, const char* obj_filename) +setup_ground + (struct htrdr_atmosphere_ground* ground, + struct htrdr_materials* mats, + const char* obj_filename) { struct aw_obj_desc desc; struct htable_vertex vertices; @@ -441,21 +444,24 @@ setup_ground(struct htrdr_atmosphere_ground* ground, const char* obj_filename) size_t iusemtl; res_T res = RES_OK; - ASSERT(obj_filename); + ASSERT(ground && mats && obj_filename); - htable_vertex_init(ground->htrdr->allocator, &vertices); - darray_double_init(ground->htrdr->allocator, &positions); - darray_size_t_init(ground->htrdr->allocator, &indices); + htable_vertex_init(htrdr_get_allocator(ground->htrdr), &vertices); + darray_double_init(htrdr_get_allocator(ground->htrdr), &positions); + darray_size_t_init(htrdr_get_allocator(ground->htrdr), &indices); - res = aw_obj_create(&ground->htrdr->logger, ground->htrdr->allocator, - ground->htrdr->verbose, &obj); + res = aw_obj_create + (htrdr_get_logger(ground->htrdr), + htrdr_get_allocator(ground->htrdr), + htrdr_get_verbosity_level(ground->htrdr), + &obj); if(res != RES_OK) { htrdr_log_err(ground->htrdr, "Could not create the obj loader -- %s.\n", res_to_cstr(res)); goto error; } - res = s3d_scene_create(ground->htrdr->s3d, &scene); + res = s3d_scene_create(htrdr_get_s3d(ground->htrdr), &scene); if(res != RES_OK) { htrdr_log_err(ground->htrdr, "Could not create the Star-3D scene -- %s.\n", res_to_cstr(res)); @@ -483,7 +489,7 @@ setup_ground(struct htrdr_atmosphere_ground* ground, const char* obj_filename) AW(obj_get_mtl(obj, iusemtl , &mtl)); - res = parse_shape_interface(ground->htrdr, mtl.name, &interf); + res = parse_shape_interface(ground->htrdr, mats, mtl.name, &interf); if(res != RES_OK) goto error; res = setup_mesh @@ -552,7 +558,7 @@ release_ground(ref_T* ref) if(ground->view) S3D(scene_view_ref_put(ground->view)); htable_interface_release(&ground->interfaces); htrdr = ground->htrdr; - MEM_RM(ground->htrdr->allocator, ground); + MEM_RM(htrdr_get_allocator(ground->htrdr), ground); htrdr_ref_put(htrdr); } @@ -563,6 +569,7 @@ res_T htrdr_atmosphere_ground_create (struct htrdr* htrdr, const char* obj_filename, /* May be NULL */ + struct htrdr_materials* mats, /* May be NULL */ const int repeat_ground, /* Infinitely repeat the ground in X and Y */ struct htrdr_atmosphere_ground** out_ground) { @@ -584,15 +591,20 @@ htrdr_atmosphere_ground_create ground->repeat = repeat_ground; f3_splat(ground->lower, (float)INF); f3_splat(ground->upper,-(float)INF); - htable_interface_init(ground->htrdr->allocator, &ground->interfaces); + htable_interface_init(htrdr_get_allocator(htrdr), &ground->interfaces); htrdr_ref_get(htrdr); ground->htrdr = htrdr; if(!obj_filename) goto exit; + if(!mats) { + res = RES_BAD_ARG; + goto error; + } + htrdr_log(ground->htrdr, "Loading ground geometry from `%s'.\n",obj_filename); time_current(&t0); - res = setup_ground(ground, obj_filename); + res = setup_ground(ground, mats, obj_filename); if(res != RES_OK) goto error; time_sub(&t0, time_current(&t1), &t0); time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf)); diff --git a/src/atmosphere/htrdr_atmosphere_ground.h b/src/atmosphere/htrdr_atmosphere_ground.h @@ -24,6 +24,7 @@ struct htrdr; struct htrdr_atmosphere_ground; struct htrdr_interface; +struct htrdr_materials; struct s3d_hit; struct ssf_bsdf; @@ -31,6 +32,7 @@ extern LOCAL_SYM res_T htrdr_atmosphere_ground_create (struct htrdr* htrdr, const char* obj_filename, /* May be NULL <=> No ground geometry */ + struct htrdr_materials* mats, /* May be NULL if no ground geometry */ const int repeat_ground, /* Infinitely repeat the ground in X and Y */ struct htrdr_atmosphere_ground** ground); diff --git a/src/atmosphere/htrdr_atmosphere_main.c b/src/atmosphere/htrdr_atmosphere_main.c @@ -15,9 +15,12 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#include "htrdr.h" -#include "htrdr_atmosphere.h" -#include "htrdr_atmosphere_args.h" +#include "atmosphere/htrdr_atmosphere.h" +#include "atmosphere/htrdr_atmosphere_args.h" + +#include "core/htrdr.h" + +#include <rsys/mem_allocator.h> int main(int argc, char** argv) @@ -33,9 +36,9 @@ main(int argc, char** argv) res = htrdr_mpi_init(argc, argv); if(res != RES_OK) goto error; - res = htrdr_atmosphere_args_init(&cmd.args, argc, argv); + res = htrdr_atmosphere_args_init(&cmd_args, argc, argv); if(res != RES_OK) goto error; - if(args.quit) goto exit; + if(cmd_args.quit) goto exit; htrdr_args.nthreads = cmd_args.nthreads; htrdr_args.verbose = cmd_args.verbose; diff --git a/src/atmosphere/htrdr_atmosphere_sun.c b/src/atmosphere/htrdr_atmosphere_sun.c @@ -15,8 +15,11 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#include "htrdr_atmosphere_c.h" -#include "htrdr_atmosphere_sun.h" +#include "atmosphere/htrdr_atmosphere_c.h" +#include "atmosphere/htrdr_atmosphere_sun.h" + +#include "core/htrdr.h" +#include "core/htrdr_log.h" #include <rsys/algorithm.h> #include <rsys/double33.h> @@ -47,7 +50,7 @@ release_sun(ref_T* ref) ASSERT(ref); sun = CONTAINER_OF(ref, struct htrdr_atmosphere_sun, ref); htrdr = sun->htrdr; - MEM_RM(htrdr_allocator(htrdr), sun); + MEM_RM(htrdr_get_allocator(htrdr), sun); htrdr_ref_put(htrdr); } diff --git a/src/atmosphere/htrdr_atmosphere_sun.h b/src/atmosphere/htrdr_atmosphere_sun.h @@ -27,7 +27,7 @@ struct ssp_rng; extern LOCAL_SYM res_T htrdr_atmosphere_sun_create - (struct htrdr_atmosphere* sun, + (struct htrdr* htrdr, struct htrdr_atmosphere_sun** out_sun); extern LOCAL_SYM void diff --git a/src/core/htrdr.c b/src/core/htrdr.c @@ -244,6 +244,7 @@ release_htrdr(ref_T* ref) struct htrdr* htrdr = CONTAINER_OF(ref, struct htrdr, ref); ASSERT(ref); + if(htrdr->s3d) S3D(device_ref_put(htrdr->s3d)); release_mpi(htrdr); if(htrdr->lifo_allocators) { size_t i; @@ -308,11 +309,11 @@ htrdr_create size_t ithread; int nthreads_max; res_T res = RES_OK; - ASSERT(args && htrdr); + ASSERT(args && out_htrdr); allocator = mem_allocator ? mem_allocator : &mem_default_allocator; - htrdr = MEM_CALLOC(allocator, 0, sizeof(*htrdr)); - if(htrdr) { + htrdr = MEM_CALLOC(allocator, 1, sizeof(*htrdr)); + if(!htrdr) { fprintf(stderr, HTRDR_LOG_ERROR_PREFIX "Could not allocate the htrdr data structure.\n"); res = RES_MEM_ERR; @@ -329,6 +330,17 @@ htrdr_create res = init_mpi(htrdr); if(res != RES_OK) goto error; + /* Disable the Star-3D verbosity since the Embree backend prints some messages + * on stdout rather than stderr. This is annoying since stdout may be used by + * htrdr to write output data */ + res = s3d_device_create(&htrdr->logger, htrdr->allocator, 0, &htrdr->s3d); + if(res != RES_OK) { + htrdr_log_err(htrdr, + "%s: could not create the Star-3D device -- %s.\n", + FUNC_NAME, res_to_cstr(res)); + goto error; + } + htrdr->lifo_allocators = MEM_CALLOC (htrdr->allocator, htrdr->nthreads, sizeof(*htrdr->lifo_allocators)); if(!htrdr->lifo_allocators) { @@ -424,6 +436,13 @@ htrdr_get_verbosity_level(const struct htrdr* htrdr) return htrdr->verbose; } +struct s3d_device* +htrdr_get_s3d(struct htrdr* htrdr) +{ + ASSERT(htrdr); + return htrdr->s3d; +} + const char* htrdr_mpi_error_string(struct htrdr* htrdr, const int mpi_err) { diff --git a/src/core/htrdr.h b/src/core/htrdr.h @@ -136,6 +136,10 @@ HTRDR_API int htrdr_get_verbosity_level (const struct htrdr* htrdr); +HTRDR_API struct s3d_device* +htrdr_get_s3d + (struct htrdr* htrdr); + HTRDR_API res_T htrdr_open_output_stream (struct htrdr* htrdr, diff --git a/src/core/htrdr_accum.h b/src/core/htrdr_accum.h @@ -19,6 +19,7 @@ #define HTRDR_ACCUM_H #include <rsys/rsys.h> +#include <rsys/math.h> /* Monte carlo accumulator */ struct htrdr_accum { diff --git a/src/core/htrdr_args.c b/src/core/htrdr_args.c @@ -24,7 +24,6 @@ #include <rsys/cstr.h> #include <rsys/double3.h> -#include <getopt.h> #include <string.h> /******************************************************************************* diff --git a/src/core/htrdr_c.h b/src/core/htrdr_c.h @@ -35,7 +35,11 @@ enum htrdr_mpi_message { HTRDR_MPI_TILE_DATA }; +struct s3d_device; + struct htrdr { + struct s3d_device* s3d; + unsigned nthreads; /* #threads of the process */ int mpi_rank; /* Rank of the process in the MPI group */ diff --git a/src/core/htrdr_draw_map.h b/src/core/htrdr_draw_map.h @@ -49,12 +49,14 @@ typedef void struct htrdr_draw_map_args { htrdr_draw_pixel_T draw_pixel; + struct ssp_rng* rng; size_t spp; /* Samples per pixel */ void* context; /* User defined data */ }; #define HTRDR_DRAW_MAP_ARGS_NULL__ { \ NULL, /* Draw pixel functor */ \ + NULL, /* Random Number Generator */ \ 0, /* #Samples per pixel */ \ NULL /* User defined data */ \ } @@ -66,7 +68,7 @@ struct htrdr; struct htrdr_buffer; static INLINE int -htrdr_draw_pixel_args_check(struct htrdr_draw_pixel_args* args) +htrdr_draw_pixel_args_check(const struct htrdr_draw_pixel_args* args) { return args && args->pixel_normalized_size[0] > 0 diff --git a/src/core/htrdr_log.h b/src/core/htrdr_log.h @@ -39,7 +39,7 @@ htrdr_log #endif ; -extern LOCAL_SYM void +HTRDR_API void htrdr_log_err (struct htrdr* htrdr, const char* msg, @@ -49,7 +49,7 @@ htrdr_log_err #endif ; -extern LOCAL_SYM void +HTRDR_API void htrdr_log_warn (struct htrdr* htrdr, const char* msg,