commit 1d182b5d7eb7d6d908b4bf638271c7c3f6184c59
parent 585067cb887fb9a0ae78ae6acb8b8bea0b96afb7
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 2 Jun 2020 15:28:31 +0200
Refactor how the spectral dimension is handled
Diffstat:
8 files changed, 283 insertions(+), 219 deletions(-)
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -102,6 +102,23 @@ log_msg
}
}
+static enum htsky_spectral_type
+htrdr_to_sky_spectral_type(const enum htrdr_spectral_type type)
+{
+ enum htsky_spectral_type spectype;
+ switch(type) {
+ case HTRDR_SPECTRAL_LW:
+ spectype = HTSKY_SPECTRAL_LW;
+ break;
+ case HTRDR_SPECTRAL_SW:
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ spectype = HTSKY_SPECTRAL_SW;
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return spectype;
+}
+
static res_T
dump_buffer
(struct htrdr* htrdr,
@@ -117,13 +134,8 @@ dump_buffer
ASSERT(htrdr && buf && stream_name && stream);
(void)stream_name;
- if(!htrdr->is_image) {
- pixsz = sizeof(struct htrdr_pixel_integ);
- pixal = ALIGNOF(struct htrdr_pixel_integ);
- } else {
- pixsz = sizeof(struct htrdr_pixel_image);
- pixal = ALIGNOF(struct htrdr_pixel_image);
- }
+ pixsz = htrdr_spectral_type_get_pixsz(htrdr->spectral_type);
+ pixal = htrdr_spectral_type_get_pixal(htrdr->spectral_type);
htrdr_buffer_get_layout(buf, &layout);
if(layout.elmt_size != pixsz || layout.alignment != pixal) {
@@ -140,8 +152,8 @@ dump_buffer
struct htrdr_estimate pix_time = HTRDR_ESTIMATE_NULL;
const struct htrdr_accum* pix_time_acc = NULL;
- if(!htrdr->is_image){
- const struct htrdr_pixel_integ* pix = htrdr_buffer_at(buf, x, y);
+ if(htrdr->spectral_type != HTRDR_SPECTRAL_SW_CIE_XYZ){
+ const struct htrdr_pixel_xwave* pix = htrdr_buffer_at(buf, x, y);
fprintf(stream, "%g %g ",
pix->radiance_temperature.E, pix->radiance_temperature.SE);
fprintf(stream, "%g %g ", pix->radiance.E, pix->radiance.SE);
@@ -414,6 +426,7 @@ htrdr_init
htrdr->grid_max_definition[0] = args->grid_max_definition[0];
htrdr->grid_max_definition[1] = args->grid_max_definition[1];
htrdr->grid_max_definition[2] = args->grid_max_definition[2];
+ htrdr->spectral_type = args->spectral_type;
res = init_mpi(htrdr);
if(res != RES_OK) goto error;
@@ -484,62 +497,37 @@ htrdr_init
htsky_args.nthreads = htrdr->nthreads;
htsky_args.repeat_clouds = args->repeat_clouds;
htsky_args.verbose = htrdr->mpi_rank == 0 ? args->verbose : 0;
- /* should the sky load short or long wave data ? */
- /* if longwave is degenerated => sw ; else : lw */
- if(args->wlen_lw_range[0] > args->wlen_lw_range[1]) {
- htsky_args.is_long_wave = 0 ;
- htsky_args.wlen_range[0] = args->wlen_sw_range[0];
- htsky_args.wlen_range[1] = args->wlen_sw_range[1];
- if(args->wlen_sw_range[0] > args->wlen_sw_range[1]) { /* image */
- htrdr->is_image = 1 ;
- } else {
- htrdr->is_image = 0 ;
- }
- } else {
- htsky_args.is_long_wave = 1 ;
- htsky_args.wlen_range[0] = args->wlen_lw_range[0];
- htsky_args.wlen_range[1] = args->wlen_lw_range[1];
- htrdr->is_image=0;
- }
-
+ htsky_args.spectral_type = htrdr_to_sky_spectral_type(args->spectral_type);
+ htsky_args.wlen_range[0] = args->wlen_range[0];
+ htsky_args.wlen_range[1] = args->wlen_range[1];
res = htsky_create(&htrdr->logger, htrdr->allocator, &htsky_args, &htrdr->sky);
if(res != RES_OK) goto error;
- if(htrdr->is_image) {
- const double* range = HTRDR_CIE_XYZ_RANGE_DEFAULT;
- size_t n;
+ htrdr->wlen_range_m[0] = args->wlen_range[0]*1e-9; /* Convert in meters */
+ htrdr->wlen_range_m[1] = args->wlen_range[1]*1e-9; /* Convert in meters */
- n = (size_t)(range[1] - range[0]);
- res = htrdr_cie_xyz_create(htrdr, range, n, &htrdr->cie);
+ if(htrdr->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ) {
+ size_t n;
+ n = (size_t)(args->wlen_range[1] - args->wlen_range[0]);
+ res = htrdr_cie_xyz_create(htrdr, args->wlen_range, n, &htrdr->cie);
if(res != RES_OK) goto error;
} else {
- if(htsky_is_long_wave(htrdr->sky)) { /* Long wave random variate */
- const double Tref=290 ; /* In Kelvin */
- size_t n;
-
- htrdr->wlen_range_m[0] = args->wlen_lw_range[0]*1e-9; /* Convert in meters */
- htrdr->wlen_range_m[1] = args->wlen_lw_range[1]*1e-9; /* Convert in meters */
- ASSERT(htrdr->wlen_range_m[0] <= htrdr->wlen_range_m[1]);
- n = (size_t)(args->wlen_lw_range[1] - args->wlen_lw_range[0]);
-
- res = htrdr_ran_wlen_create
- (htrdr, args->wlen_lw_range, n, Tref, &htrdr->ran_wlen);
- if(res != RES_OK) goto error;
+ double Tref; /* In Kelvin */
+ size_t n;
- } else {
- const double Tref=5778 ; /* Tsun In Kelvin */
- size_t n;
+ switch(htrdr->spectral_type) {
+ case HTRDR_SPECTRAL_LW: Tref = 290; break;
+ case HTRDR_SPECTRAL_SW: Tref = 5778; /* Tsun */ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
- htrdr->wlen_range_m[0] = args->wlen_sw_range[0]*1e-9; /* Convert in meters */
- htrdr->wlen_range_m[1] = args->wlen_sw_range[1]*1e-9; /* Convert in meters */
- ASSERT(htrdr->wlen_range_m[0] <= htrdr->wlen_range_m[1]);
- n = (size_t)(args->wlen_sw_range[1] - args->wlen_sw_range[0]);
+ ASSERT(htrdr->wlen_range_m[0] <= htrdr->wlen_range_m[1]);
+ n = (size_t)(args->wlen_range[1] - args->wlen_range[0]);
- res = htrdr_ran_wlen_create
- (htrdr, args->wlen_sw_range, n, Tref, &htrdr->ran_wlen);
- if(res != RES_OK) goto error;
- }
+ res = htrdr_ran_wlen_create
+ (htrdr, args->wlen_range, n, Tref, &htrdr->ran_wlen);
+ if(res != RES_OK) goto error;
}
htrdr->lifo_allocators = MEM_CALLOC
@@ -566,16 +554,9 @@ htrdr_init
/* Create the image buffer only on the master process; the image parts
* rendered by the processes are gathered onto the master process. */
if(!htrdr->dump_vtk && htrdr->mpi_rank == 0) {
- size_t pixsz = 0; /* sizeof(pixel) */
- size_t pixal = 0; /* alignof(pixel) */
+ const size_t pixsz = htrdr_spectral_type_get_pixsz(htrdr->spectral_type);
+ const size_t pixal = htrdr_spectral_type_get_pixal(htrdr->spectral_type);
- if(!htrdr->is_image) {
- pixsz = sizeof(struct htrdr_pixel_integ);
- pixal = ALIGNOF(struct htrdr_pixel_integ);
- } else {
- pixsz = sizeof(struct htrdr_pixel_image);
- pixal = ALIGNOF(struct htrdr_pixel_image);
- }
res = htrdr_buffer_create(htrdr,
args->image.definition[0], /* Width */
args->image.definition[1], /* Height */
diff --git a/src/htrdr.h b/src/htrdr.h
@@ -30,6 +30,12 @@
#define HTRDR(Func) htrdr_ ## Func
#endif
+enum htrdr_spectral_type {
+ HTRDR_SPECTRAL_LW, /* Longwave */
+ HTRDR_SPECTRAL_SW, /* Shortwave */
+ HTRDR_SPECTRAL_SW_CIE_XYZ /* Shortwave wrt the CIE XYZ tristimulus */
+};
+
/* Forward declarations */
struct htsky;
struct htrdr_args;
@@ -56,8 +62,8 @@ struct htrdr {
struct htrdr_buffer* buf;
struct htsky* sky;
+ enum htrdr_spectral_type spectral_type;
double wlen_range_m[2]; /* Integration range in *meters* */
- int is_image ; /* XYZ Image or spectral integration? */
size_t spp; /* #samples per pixel */
size_t width; /* Image width */
diff --git a/src/htrdr_args.c b/src/htrdr_args.c
@@ -61,12 +61,6 @@ print_help(const char* cmd)
" -i <image> define the image to compute. Refer to the htrdr man\n"
" page for the list of image options\n");
printf(
-" -l WLEN_MIN,WLEN_MAX\n"
-" switch in infrared rendering for longwave in\n"
-" [WLEN_MIN, WLEN_MAX], in nanometers. By default, the\n"
-" rendering is performed for the visible part of the\n"
-" spectrum in [380, 780] nanometers.\n");
- printf(
" -M MATERIALS file listing the ground materials.\n");
printf(
" -m MIE file of Mie's data.\n");
@@ -81,11 +75,9 @@ print_help(const char* cmd)
printf(
" -r infinitely repeat the clouds along the X and Y axis.\n");
printf(
-" -s WLEN_MIN,WLEN_MAX\n"
-" switch in solar rendering for the short waves in\n"
-" [WLEN_MIN, WLEN_MAX], in nanometers. By default, the\n"
-" rendering is performed for the visible part of the\n"
-" spectrum in [380, 780] nanometers with CIE.\n");
+" -s <spectral> define the type and range of the spectral\n"
+" integration. Refer to the htrdr man page for the list\n"
+" of spectral options\n");
printf(
" -T THRESHOLD optical thickness used as threshold during the octree\n"
" building. By default its value is `%g'.\n",
@@ -223,10 +215,11 @@ parse_camera_parameter(struct htrdr_args* args, const char* str)
char* val;
char* ctx;
res_T res = RES_OK;
+ ASSERT(args);
if(strlen(str) >= sizeof(buf) -1/*NULL char*/) {
fprintf(stderr,
- "Could not duplicate the rectangle option string `%s'.\n", str);
+ "Could not duplicate the camera option string `%s'.\n", str);
res = RES_MEM_ERR;
goto error;
}
@@ -268,6 +261,82 @@ error:
}
static res_T
+parse_spectral_range(const char* str, double wlen_range[2])
+{
+ double range[2];
+ size_t len;
+ res_T res = RES_OK;
+ ASSERT(wlen_range && str);
+
+ res = cstr_to_list_double(str, ',', range, &len, 2);
+ if(res == RES_OK && len != 2) res = RES_BAD_ARG;
+ if(res == RES_OK && range[0] > range[1]) res = RES_BAD_ARG;
+ if(res != RES_OK) {
+ fprintf(stderr, "Invalid spectral range `%s'.\n", str);
+ goto error;
+ }
+ wlen_range[0] = range[0];
+ wlen_range[1] = range[1];
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+parse_spectral_parameter(struct htrdr_args* args, const char* str)
+{
+ char buf[128];
+ char* key;
+ char* val;
+ char* ctx;
+ res_T res = RES_OK;
+ ASSERT(args);
+
+ if(strlen(str) >= sizeof(buf) -1/*NULL char*/) {
+ fprintf(stderr,
+ "Could not duplicate the spectral option string `%s'.\n", str);
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ strncpy(buf, str, sizeof(buf));
+
+ key = strtok_r(buf, "=", &ctx);
+ val = strtok_r(buf, "", &ctx);
+
+ if(!strcmp(key, "cie_xyz")) {
+ args->spectral_type = HTRDR_SPECTRAL_SW_CIE_XYZ;
+ args->wlen_range[0] = HTRDR_CIE_XYZ_RANGE_DEFAULT[0];
+ args->wlen_range[1] = HTRDR_CIE_XYZ_RANGE_DEFAULT[1];
+ } else {
+ if(!val) {
+ fprintf(stderr, "Missing value to the spectral option `%s'.\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ if(!strcmp(key, "sw")) {
+ args->spectral_type = HTRDR_SPECTRAL_SW;
+ res = parse_spectral_range(val, args->wlen_range);
+ if(res != RES_OK) goto error;
+ } else if(!strcmp(val, "lw")) {
+ args->spectral_type = HTRDR_SPECTRAL_LW;
+ res = parse_spectral_range(val, args->wlen_range);
+ if(res != RES_OK) goto error;
+ } else {
+ fprintf(stderr, "Invalid spectral parameter `%s'.\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
parse_multiple_parameters
(struct htrdr_args* args,
const char* str,
@@ -371,53 +440,6 @@ error:
goto exit;
}
-static res_T
-parse_spectral_range(const char* str, double wlen_range[2])
-{
- double range[2];
- size_t len;
- res_T res = RES_OK;
- ASSERT(wlen_range && str);
-
- res = cstr_to_list_double(str, ',', range, &len, 2);
- if(res == RES_OK && len != 2) res = RES_BAD_ARG;
- if(res == RES_OK && range[0] > range[1]) res = RES_BAD_ARG;
- if(res != RES_OK) {
- fprintf(stderr, "Invalid spectral range `%s'.\n", str);
- goto error;
- }
-
- wlen_range[0] = range[0];
- wlen_range[1] = range[1];
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-static res_T
-parse_lw_range(struct htrdr_args* args, const char* str)
-{
- res_T res = RES_OK;
- ASSERT(args && str);
-
- res = parse_spectral_range(str, args->wlen_lw_range) ;
-
- return res ;
-}
-
-static res_T
-parse_sw_range(struct htrdr_args* args, const char* str)
-{
- res_T res = RES_OK;
- ASSERT(args && str);
-
- res = parse_spectral_range(str, args->wlen_sw_range) ;
-
- return res ;
-}
-
/*******************************************************************************
* Local functions
******************************************************************************/
@@ -442,7 +464,7 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
}
}
- while((opt = getopt(argc, argv, "a:C:c:D:dfg:hi:l:M:m:O:o:Rrs:T:t:V:v")) != -1) {
+ while((opt = getopt(argc, argv, "a:C:c:D:dfg:hi:M:m:O:o:Rrs:T:t:V:v")) != -1) {
switch(opt) {
case 'a': args->filename_gas = optarg; break;
case 'C':
@@ -463,9 +485,6 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
res = parse_multiple_parameters
(args, optarg, parse_image_parameter);
break;
- case 'l':
- res = parse_lw_range(args, optarg);
- break;
case 'M': args->filename_mtl = optarg; break;
case 'm': args->filename_mie = optarg; break;
case 'O': args->cache = optarg; break;
@@ -473,7 +492,8 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
case 'r': args->repeat_clouds = 1; break;
case 'R': args->repeat_ground = 1; break;
case 's':
- res = parse_sw_range(args, optarg);
+ res = parse_multiple_parameters
+ (args, optarg, parse_spectral_parameter);
break;
case 'T':
res = cstr_to_double(optarg, &args->optical_thickness);
diff --git a/src/htrdr_args.h.in b/src/htrdr_args.h.in
@@ -16,7 +16,8 @@
#ifndef HTRDR_ARGS_H
#define HTRDR_ARGS_H
-#include "htrdr_ground.h"
+#include "htrdr.h"
+#include "htrdr_cie_xyz.h"
#include <float.h>
#include <limits.h>
@@ -55,8 +56,8 @@ struct htrdr_args {
double optical_thickness; /* Threshold used during octree building */
unsigned grid_max_definition[3]; /* Maximum definition of the grid */
- double wlen_sw_range[2]; /* Spectral range of integration in nm (solar) */
- double wlen_lw_range[2]; /* Spectral range of integration in nm (thermal) */
+ enum htrdr_spectral_type spectral_type;
+ double wlen_range[2]; /* Spectral range of integration in nm */
unsigned nthreads; /* Hint on the number of threads to use */
int force_overwriting;
@@ -93,8 +94,8 @@ struct htrdr_args {
90, /* Sun elevation */ \
@HTRDR_ARGS_DEFAULT_OPTICAL_THICKNESS_THRESHOLD@, /* Optical thickness */ \
{UINT_MAX, UINT_MAX, UINT_MAX}, /* Maximum definition of the grid */ \
- {DBL_MAX, -DBL_MAX}, /* Shortwave range. Degenerated <=> CIE sw */ \
- {DBL_MAX, -DBL_MAX}, /* Long wave range. Degenerated <=> CIE sw */ \
+ HTRDR_SPECTRAL_SW_CIE_XYZ, \
+ HTRDR_CIE_XYZ_RANGE_DEFAULT__, /* Spectral range */ \
(unsigned)~0, /* #threads */ \
0, /* Force overwriting */ \
0, /* dump VTK */ \
diff --git a/src/htrdr_cie_xyz.h b/src/htrdr_cie_xyz.h
@@ -22,7 +22,9 @@ struct htrdr;
struct htrdr_cie_xyz;
/* Wavelength boundaries of the CIE XYZ color space in nanometers */
-static const double HTRDR_CIE_XYZ_RANGE_DEFAULT[2] = {380, 780};
+#define HTRDR_CIE_XYZ_RANGE_DEFAULT__ {380, 780}
+static const double HTRDR_CIE_XYZ_RANGE_DEFAULT[2] =
+ HTRDR_CIE_XYZ_RANGE_DEFAULT__;
extern LOCAL_SYM res_T
htrdr_cie_xyz_create
diff --git a/src/htrdr_compute_radiance_sw.c b/src/htrdr_compute_radiance_sw.c
@@ -275,13 +275,12 @@ htrdr_compute_radiance_sw
double L_sun; /* Sun radiance in W.m^-2.sr^-1 */
double sun_dir[3];
double ksi = 1; /* Throughput */
- double temperature=5778;
- double wlen_m = wlen * 1.e-9;
double w = 0; /* MC weight */
double g = 0; /* Asymmetry parameter of the HG phase function */
ASSERT(htrdr && rng && pos_in && dir_in && ithread < htrdr->nthreads);
- ASSERT(!htsky_is_long_wave(htrdr->sky));
+ ASSERT(htrdr->spectral_type == HTRDR_SPECTRAL_SW
+ || htrdr->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ);
CHK(RES_OK == ssf_phase_create
(&htrdr->lifo_allocators[ithread], &ssf_phase_hg, &phase_hg));
@@ -299,11 +298,7 @@ htrdr_compute_radiance_sw
htsky_get_spectral_band_bounds(htrdr->sky, iband, band_bounds);
ASSERT(band_bounds[0] <= wlen && wlen <= band_bounds[1]);
sun_solid_angle = htrdr_sun_get_solid_angle(htrdr->sun);
- if (htrdr->is_image) {
- L_sun = htrdr_sun_get_radiance(htrdr->sun, wlen);
- } else {
- L_sun = planck_monochromatic(wlen_m, temperature);
- }
+ L_sun = htrdr_sun_get_radiance(htrdr->sun, wlen);
d3_set(pos, pos_in);
d3_set(dir, dir_in);
diff --git a/src/htrdr_draw_radiance.c b/src/htrdr_draw_radiance.c
@@ -46,12 +46,12 @@
STATIC_ASSERT(IS_POW2(TILE_SIZE), TILE_SIZE_must_be_a_power_of_2);
enum pixel_format {
- PIXEL_INTEG,
+ PIXEL_XWAVE,
PIXEL_IMAGE
};
union pixel {
- struct htrdr_pixel_integ integ;
+ struct htrdr_pixel_xwave xwave;
struct htrdr_pixel_image image;
};
@@ -101,6 +101,19 @@ morton2D_encode(const uint16_t u16)
return u32;
}
+static INLINE enum pixel_format
+spectral_type_to_pixfmt(const enum htrdr_spectral_type spectral_type)
+{
+ enum pixel_format pixfmt;
+ switch(spectral_type) {
+ case HTRDR_SPECTRAL_LW: pixfmt = PIXEL_XWAVE; break;
+ case HTRDR_SPECTRAL_SW: pixfmt = PIXEL_XWAVE; break;
+ case HTRDR_SPECTRAL_SW_CIE_XYZ: pixfmt = PIXEL_IMAGE; break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return pixfmt;
+}
+
static FINLINE struct tile*
tile_create(struct mem_allocator* allocator, const enum pixel_format fmt)
{
@@ -155,9 +168,11 @@ tile_at
ASSERT(tile && x < TILE_SIZE && y < TILE_SIZE);
return tile->data.pixels + (y*TILE_SIZE + x);
}
+
static void
write_tile_data
- (struct htrdr_buffer* buf,
+ (struct htrdr* htrdr,
+ struct htrdr_buffer* buf,
const struct tile_data* tile_data)
{
struct htrdr_buffer_layout layout = HTRDR_BUFFER_LAYOUT_NULL;
@@ -165,13 +180,12 @@ write_tile_data
size_t irow_tile;
size_t ncols_tile, nrows_tile;
char* buf_mem;
- ASSERT(buf && tile_data);
+ ASSERT(htrdr && buf && tile_data);
+ (void)htrdr;
htrdr_buffer_get_layout(buf, &layout);
buf_mem = htrdr_buffer_get_data(buf);
- ASSERT(layout.elmt_size == (tile_data->format == PIXEL_INTEG
- ? sizeof(struct htrdr_pixel_integ)
- : sizeof(struct htrdr_pixel_image)));
+ ASSERT(layout.elmt_size==htrdr_spectral_type_get_pixsz(htrdr->spectral_type));
/* Compute the row/column of the tile origin into the buffer */
icol = tile_data->x * (size_t)TILE_SIZE;
@@ -190,8 +204,8 @@ write_tile_data
FOR_EACH(x, 0, ncols_tile) {
switch(tile_data->format) {
- case PIXEL_INTEG:
- ((struct htrdr_pixel_integ*)buf_col)[x] = tile_row[x].integ;
+ case PIXEL_XWAVE:
+ ((struct htrdr_pixel_xwave*)buf_col)[x] = tile_row[x].xwave;
break;
case PIXEL_IMAGE:
((struct htrdr_pixel_image*)buf_col)[x] = tile_row[x].image;
@@ -445,17 +459,18 @@ mpi_gather_tiles
LIST_FOR_EACH(node, tiles) {
struct tile* t = CONTAINER_OF(node, struct tile, node);
- write_tile_data(buf, &t->data);
+ write_tile_data(htrdr, buf, &t->data);
++itile;
}
if(itile != ntiles) {
+ enum pixel_format pixfmt;
ASSERT(htrdr->mpi_nprocs > 1);
/* Create a temporary tile to receive the tile data computed by the
* concurrent MPI processes */
- tile = tile_create
- (htrdr->allocator, htrdr->is_image ? PIXEL_IMAGE : PIXEL_INTEG) ;
+ pixfmt = spectral_type_to_pixfmt(htrdr->spectral_type);
+ tile = tile_create(htrdr->allocator, pixfmt);
if(!tile) {
res = RES_MEM_ERR;
htrdr_log_err(htrdr,
@@ -468,7 +483,7 @@ mpi_gather_tiles
FOR_EACH(itile, itile, ntiles) {
MPI(Recv(&tile->data, (int)msg_sz, MPI_CHAR, MPI_ANY_SOURCE,
HTRDR_MPI_TILE_DATA, MPI_COMM_WORLD, MPI_STATUS_IGNORE));
- write_tile_data(buf, &tile->data);
+ write_tile_data(htrdr, buf, &tile->data);
}
}
}
@@ -495,7 +510,7 @@ draw_pixel_image
struct htrdr_accum time;
size_t ichannel;
ASSERT(ipix && ipix && pix_sz && cam && rng && pixel);
- ASSERT(!htsky_is_long_wave(htrdr->sky));
+ ASSERT(htrdr->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ);
/* Reset accumulators */
XYZ[0] = HTRDR_ACCUM_NULL;
@@ -576,8 +591,39 @@ draw_pixel_image
pixel->time = time;
}
+
+static INLINE double
+radiance_temperature
+ (struct htrdr* htrdr,
+ const double radiance) /* In W/m^2/sr */
+{
+ double temperature = 0;
+ double radiance_avg = radiance;
+ res_T res = RES_OK;
+ ASSERT(htrdr && radiance >= 0);
+
+ /* From integrated radiance to average radiance in W/m^2/sr/m */
+ if(htrdr->wlen_range_m[0] != htrdr->wlen_range_m[1]) { /* !monochromatic */
+ radiance_avg /= (htrdr->wlen_range_m[1] - htrdr->wlen_range_m[0]);
+ }
+
+ res = brightness_temperature
+ (htrdr,
+ htrdr->wlen_range_m[0],
+ htrdr->wlen_range_m[1],
+ radiance_avg,
+ &temperature);
+ if(res != RES_OK) {
+ htrdr_log_warn(htrdr,
+ "Could not compute the brightness temperature for the radiance %g.\n",
+ radiance_avg);
+ temperature = 0;
+ }
+ return temperature;
+}
+
static void
-draw_pixel_integ
+draw_pixel_xwave
(struct htrdr* htrdr,
const size_t ithread,
const size_t ipix[2],
@@ -585,14 +631,15 @@ draw_pixel_integ
const struct htrdr_camera* cam,
const size_t spp,
struct ssp_rng* rng,
- struct htrdr_pixel_integ* pixel)
+ struct htrdr_pixel_xwave* pixel)
{
struct htrdr_accum radiance;
struct htrdr_accum time;
size_t isamp;
double temp_min, temp_max;
ASSERT(ipix && ipix && pix_sz && cam && rng && pixel);
- ASSERT(!(htrdr->is_image));
+ ASSERT(htrdr->spectral_type == HTRDR_SPECTRAL_LW
+ || htrdr->spectral_type == HTRDR_SPECTRAL_SW);
/* Reset the pixel accumulators */
radiance = HTRDR_ACCUM_NULL;
@@ -633,25 +680,23 @@ draw_pixel_integ
iband = htsky_find_spectral_band(htrdr->sky, wlen);
iquad = htsky_spectral_band_sample_quadrature(htrdr->sky, r2, iband);
- /* Compute the integrated luminance in W/m^2/sr/m */
- if(htsky_is_long_wave(htrdr->sky)) {
- weight = htrdr_compute_radiance_lw(htrdr, ithread, rng, ray_org, ray_dir,
- wlen, iband, iquad);
- } else {
- weight = htrdr_compute_radiance_sw(htrdr, ithread, rng, ray_org, ray_dir,
- wlen, iband, iquad);
+ /* Compute the spectral radiance in W/m^2/sr/m */
+ switch(htrdr->spectral_type) {
+ case HTRDR_SPECTRAL_LW:
+ weight = htrdr_compute_radiance_lw(htrdr, ithread, rng, ray_org,
+ ray_dir, wlen, iband, iquad);
+ break;
+ case HTRDR_SPECTRAL_SW:
+ weight = htrdr_compute_radiance_sw(htrdr, ithread, rng, ray_org,
+ ray_dir, wlen, iband, iquad);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
}
+ ASSERT(weight >= 0);
/* Importance sampling: correct weight with pdf */
weight /= band_pdf; /* In W/m^2/sr */
- /* From integrated radiance to average radiance in W/m^2/sr/m */
- if(htrdr->wlen_range_m[0] != htrdr->wlen_range_m[1]) {
- /* Is not monochromatic */
- weight /= (htrdr->wlen_range_m[1] - htrdr->wlen_range_m[0]) ;
- }
- ASSERT(weight >= 0);
-
/* End the registration of the per realisation time */
time_sub(&t0, time_current(&t1), &t0);
usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
@@ -674,32 +719,13 @@ draw_pixel_integ
pixel->time = time;
/* Compute the brightness_temperature of the pixel and estimate its standard
- * error if the sources were in the medium (i.e., is_long_wave) */
- if(htsky_is_long_wave(htrdr->sky)) {
- #define BRIGHTNESS_TEMPERATURE(Radiance, Temperature) { \
- res_T res = brightness_temperature \
- (htrdr, \
- htrdr->wlen_range_m[0], \
- htrdr->wlen_range_m[1], \
- (Radiance), \
- &(Temperature)); \
- if(res != RES_OK) { \
- htrdr_log_warn(htrdr, \
- "Could not compute the brightness temperature for the radiance %g.\n", \
- (Radiance)); \
- (Temperature) = 0; \
- } \
- } (void)0
- BRIGHTNESS_TEMPERATURE(pixel->radiance.E, pixel->radiance_temperature.E);
- BRIGHTNESS_TEMPERATURE(pixel->radiance.E - pixel->radiance.SE, temp_min);
- BRIGHTNESS_TEMPERATURE(pixel->radiance.E + pixel->radiance.SE, temp_max);
+ * error if the sources were in the medium (<=> longwave) */
+ if(htrdr->spectral_type == HTRDR_SPECTRAL_LW) {
+ pixel->radiance_temperature.E = radiance_temperature(htrdr, pixel->radiance.E);
+ temp_min = radiance_temperature(htrdr, pixel->radiance.E - pixel->radiance.SE);
+ temp_max = radiance_temperature(htrdr, pixel->radiance.E + pixel->radiance.SE);
pixel->radiance_temperature.SE = temp_max - temp_min;
- #undef BRIGHTNESS_TEMPERATURE
}
-
- /* Transform the pixel radiance from W/m^2/sr/m to W/m^/sr/nm */
- pixel->radiance.E *= 1.0e-9;
- pixel->radiance.SE *= 1.0e-9;
}
static res_T
@@ -741,10 +767,15 @@ draw_tile
ipix[1] = tile_org[1] + ipix_tile[1];
/* Draw the pixel */
- if(!(htrdr->is_image)) {
- draw_pixel_integ(htrdr, ithread, ipix, pix_sz, cam, spp, rng, &pixel->integ);
- } else {
- draw_pixel_image(htrdr, ithread, ipix, pix_sz, cam, spp, rng, &pixel->image);
+ switch(htrdr->spectral_type) {
+ case HTRDR_SPECTRAL_LW:
+ case HTRDR_SPECTRAL_SW:
+ draw_pixel_xwave(htrdr, ithread, ipix, pix_sz, cam, spp, rng, &pixel->xwave);
+ break;
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ draw_pixel_image(htrdr, ithread, ipix, pix_sz, cam, spp, rng, &pixel->image);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
}
}
return RES_OK;
@@ -790,7 +821,7 @@ draw_image
htrdr->mpi_working_procs[htrdr->mpi_rank] = 0;
--htrdr->mpi_nworking_procs;
- pixfmt = htrdr->is_image ? PIXEL_IMAGE : PIXEL_INTEG ;
+ pixfmt = spectral_type_to_pixfmt(htrdr->spectral_type);
omp_set_num_threads((int)nthreads);
#pragma omp parallel
@@ -944,19 +975,13 @@ htrdr_draw_radiance
proc_work_init(htrdr->allocator, &work);
if(htrdr->mpi_rank == 0) {
- size_t pixsz, pixal;
+ const size_t pixsz = htrdr_spectral_type_get_pixsz(htrdr->spectral_type);
+ const size_t pixal = htrdr_spectral_type_get_pixal(htrdr->spectral_type);
+
htrdr_buffer_get_layout(buf, &layout);
ASSERT(layout.width || layout.height || layout.elmt_size);
ASSERT(layout.width == width && layout.height == height);
- if(!(htrdr->is_image)) {
- pixsz = sizeof(struct htrdr_pixel_integ);
- pixal = ALIGNOF(struct htrdr_pixel_integ);
- } else {
- pixsz = sizeof(struct htrdr_pixel_image);
- pixal = ALIGNOF(struct htrdr_pixel_image);
- }
-
if(layout.elmt_size != pixsz || layout.alignment < pixal) {
htrdr_log_err(htrdr, "%s: invalid buffer layout.\n", FUNC_NAME);
res = RES_BAD_ARG;
diff --git a/src/htrdr_solve.h b/src/htrdr_solve.h
@@ -35,6 +35,13 @@ struct htrdr_estimate {
};
static const struct htrdr_estimate HTRDR_ESTIMATE_NULL;
+struct htrdr_pixel_xwave {
+ struct htrdr_estimate radiance; /* In W/m^2/sr */
+ struct htrdr_estimate radiance_temperature; /* In K */
+ struct htrdr_accum time; /* In microseconds */
+};
+static const struct htrdr_pixel_xwave HTRDR_PIXEL_XWAVE_NULL;
+
struct htrdr_pixel_image {
struct htrdr_estimate X; /* In W/m^2/sr */
struct htrdr_estimate Y; /* In W/m^2/sr */
@@ -43,13 +50,6 @@ struct htrdr_pixel_image {
};
static const struct htrdr_pixel_image HTRDR_PIXEL_IMAGE_NULL;
-struct htrdr_pixel_integ {
- struct htrdr_estimate radiance; /* In W/m^2/sr */
- struct htrdr_estimate radiance_temperature; /* In K */
- struct htrdr_accum time; /* In microseconds */
-};
-static const struct htrdr_pixel_integ HTRDR_PIXEL_INTEG_NULL;
-
/* Forward declarations */
struct htrdr;
struct htrdr_camera;
@@ -120,4 +120,38 @@ htrdr_accum_get_estimation
}
}
+static INLINE size_t
+htrdr_spectral_type_get_pixsz(const enum htrdr_spectral_type type)
+{
+ size_t sz = 0;
+ switch(type) {
+ case HTRDR_SPECTRAL_LW:
+ case HTRDR_SPECTRAL_SW:
+ sz = sizeof(struct htrdr_pixel_xwave);
+ break;
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ sz = sizeof(struct htrdr_pixel_image);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return sz;
+}
+
+static INLINE size_t
+htrdr_spectral_type_get_pixal(const enum htrdr_spectral_type type)
+{
+ size_t al = 0;
+ switch(type) {
+ case HTRDR_SPECTRAL_LW:
+ case HTRDR_SPECTRAL_SW:
+ al = ALIGNOF(struct htrdr_pixel_xwave);
+ break;
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ al = ALIGNOF(struct htrdr_pixel_image);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return al;
+}
+
#endif /* HTRDR_SOLVE_H */
