commit 476b6100867ab3e4124cbd8ebd9c40d068950d24
parent 019a06f1ab1ff62aabcfdabb4636123ad1718ca8
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 19 Mar 2020 10:45:40 +0100
Add the brightness_temperature function and refactor the draw function
Diffstat:
4 files changed, 253 insertions(+), 98 deletions(-)
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -434,9 +434,11 @@ htrdr_init
const struct htrdr_args* args,
struct htrdr* htrdr)
{
+ size_t nbands, iband0, iband1;
struct htsky_args htsky_args = HTSKY_ARGS_DEFAULT;
double proj_ratio;
double sun_dir[3];
+ double wlen0[2], wlen1[2];
const char* output_name = NULL;
size_t ithread;
int nthreads_max;
@@ -554,6 +556,15 @@ htrdr_init
res = htsky_create(&htrdr->logger, htrdr->allocator, &htsky_args, &htrdr->sky);
if(res != RES_OK) goto error;
+ /* Fetch the wavelengths integration range */
+ nbands = htsky_get_spectral_bands_count(htrdr->sky);
+ iband0 = htsky_get_spectral_band_id(htrdr->sky, 0);
+ iband1 = htsky_get_spectral_band_id(htrdr->sky, nbands-1);
+ HTSKY(get_spectral_band_bounds(htrdr->sky, iband0, wlen0));
+ HTSKY(get_spectral_band_bounds(htrdr->sky, iband1, wlen1));
+ htrdr->wlen_range[0] = wlen0[0];
+ htrdr->wlen_range[1] = wlen1[1];
+
if(htsky_is_long_wave(htrdr->sky)) {
/* Define the CDF used to sample a long wave band */
res = setup_lw_cdf(htrdr);
@@ -731,7 +742,67 @@ htrdr_fflush(struct htrdr* htrdr, FILE* stream)
/*******************************************************************************
* Local functions
******************************************************************************/
-extern LOCAL_SYM res_T
+res_T
+brightness_temperature
+ (struct htrdr* htrdr,
+ const double lambda_min,
+ const double lambda_max,
+ const double radiance,
+ double* temperature)
+{
+ const size_t MAX_ITER = 100;
+ const double epsilon_T = 1e-4; /* In K */
+ const double epsilon_B = radiance * 1e-8;
+ double T, T0, T1, T2;
+ double B, B0;
+ size_t i;
+ res_T res = RES_OK;
+ ASSERT(temperature && lambda_min <= lambda_max);
+
+ /* Search for a brightness temperature whose radiance is greater than or
+ * equal to the estimated radiance */
+ T2 = 200;
+ FOR_EACH(i, 0, MAX_ITER) {
+ const double B2 = planck(lambda_min, lambda_max, T2);
+ if(B2 >= radiance) break;
+ T2 *= 2;
+ }
+ if(i >= MAX_ITER) { res = RES_BAD_OP; goto error; }
+
+ B0 = T0 = T1 = 0;
+ FOR_EACH(i, 0, MAX_ITER) {
+ T = (T1+T2)*0.5;
+ B = planck(lambda_min, lambda_max, T);
+
+ if(B < radiance) {
+ T1 = T;
+ } else {
+ T2 = T;
+ }
+
+ if(fabs(T-T0) < epsilon_T || fabs(B-B0) < epsilon_B)
+ break;
+
+ T0 = T;
+ B0 = B;
+ }
+ if(i >= MAX_ITER) { res = RES_BAD_OP; goto error; }
+
+ *temperature = T;
+
+exit:
+ return res;
+error:
+ htrdr_log_err(htrdr,
+ "Could not compute the brightness temperature for the radiance %g "
+ "estimated over [%g, %g] nanometers.\n",
+ radiance,
+ lambda_min*1e9,
+ lambda_max*1e9);
+ goto exit;
+}
+
+res_T
open_output_stream
(struct htrdr* htrdr,
const char* filename,
diff --git a/src/htrdr.h b/src/htrdr.h
@@ -54,6 +54,7 @@ struct htrdr {
struct htrdr_buffer* buf;
struct htsky* sky;
+ double wlen_range[2]; /* Integration range in nanometers */
struct darray_double lw_cdf; /* CDF to sample a Long Waves band */
diff --git a/src/htrdr_c.h b/src/htrdr_c.h
@@ -35,6 +35,7 @@ enum htrdr_mpi_message {
enum htrdr_estimate {
HTRDR_ESTIMATE_X,
+ HTRDR_ESTIMATE_RADIANCE = HTRDR_ESTIMATE_X,
HTRDR_ESTIMATE_Y,
HTRDR_ESTIMATE_Z,
HTRDR_ESTIMATE_TIME, /* Time per realisation */
@@ -159,9 +160,18 @@ planck
const double BOLTZMANN_CONSTANT = 5.6696e-8; /* W/m^2/K^4 */
ASSERT(lambda_min < lambda_max && temperature > 0);
return blackbody_fraction(lambda_min, lambda_max, temperature)
- * BOLTZMANN_CONSTANT * T4;
+ * BOLTZMANN_CONSTANT * T4 / PI; /* In W/m^2/sr */
}
+extern LOCAL_SYM res_T
+brightness_temperature
+ (struct htrdr* htrdr,
+ const double lambda_min, /* In meter */
+ const double lambda_max, /* In meter */
+ /* Integrated over [lambda_min, lambda_max]. In W/m^2/sr */
+ const double radiance,
+ double* temperature);
+
extern LOCAL_SYM res_T
open_output_stream
(struct htrdr* htrdr,
diff --git a/src/htrdr_draw_radiance.c b/src/htrdr_draw_radiance.c
@@ -473,6 +473,170 @@ sample_lw_spectral_interval(struct htrdr* htrdr, const double r)
return htsky_get_spectral_band_id(htrdr->sky, i);
}
+static void
+draw_pixel_sw
+ (struct htrdr* htrdr,
+ const size_t ithread,
+ const size_t ipix[2],
+ const double pix_sz[2], /* Size of a pixel in the normalized image plane */
+ const struct htrdr_camera* cam,
+ const size_t spp,
+ struct ssp_rng* rng,
+ struct htrdr_accum* pix_accums)
+{
+ size_t ichannel;
+ ASSERT(ipix && ipix && pix_sz && cam && rng && pix_accums);
+ ASSERT(!htsky_is_long_wave(htrdr->sky));
+
+ /* Reset the pixel accumulators */
+ pix_accums[HTRDR_ESTIMATE_X] = HTRDR_ACCUM_NULL;
+ pix_accums[HTRDR_ESTIMATE_Y] = HTRDR_ACCUM_NULL;
+ pix_accums[HTRDR_ESTIMATE_Z] = HTRDR_ACCUM_NULL;
+ pix_accums[HTRDR_ESTIMATE_TIME] = HTRDR_ACCUM_NULL;
+
+ FOR_EACH(ichannel, 0, 3) {
+ /* Check that the X, Y and Z estimates are stored in accumulators 0, 1 et
+ * 2, respectively */
+ STATIC_ASSERT
+ ( HTRDR_ESTIMATE_X == 0
+ && HTRDR_ESTIMATE_Y == 1
+ && HTRDR_ESTIMATE_Z == 2,
+ Unexpected_htrdr_estimate_enumerate);
+ size_t isamp;
+
+ FOR_EACH(isamp, 0, spp) {
+ struct time t0, t1;
+ double pix_samp[2];
+ double ray_org[3];
+ double ray_dir[3];
+ double weight;
+ double r0, r1;
+ size_t iband;
+ size_t iquad;
+ double usec;
+
+ /* Begin the registration of the time spent to in the realisation */
+ time_current(&t0);
+
+ /* Sample a position into the pixel, in the normalized image plane */
+ pix_samp[0] = ((double)ipix[0] + ssp_rng_canonical(rng)) * pix_sz[0];
+ pix_samp[1] = ((double)ipix[1] + ssp_rng_canonical(rng)) * pix_sz[1];
+
+ /* Generate a ray starting from the pinhole camera and passing through the
+ * pixel sample */
+ htrdr_camera_ray(cam, pix_samp, ray_org, ray_dir);
+
+ r0 = ssp_rng_canonical(rng);
+ r1 = ssp_rng_canonical(rng);
+
+ /* Sample a spectral band and a quadrature point */
+ switch(ichannel) {
+ case 0:
+ htsky_sample_sw_spectral_data_CIE_1931_X
+ (htrdr->sky, r0, r1, &iband, &iquad);
+ break;
+ case 1:
+ htsky_sample_sw_spectral_data_CIE_1931_Y
+ (htrdr->sky, r0, r1, &iband, &iquad);
+ break;
+ case 2:
+ htsky_sample_sw_spectral_data_CIE_1931_Z
+ (htrdr->sky, r0, r1, &iband, &iquad);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+
+ /* Compute the luminance */
+ weight = htrdr_compute_radiance_sw
+ (htrdr, ithread, rng, ray_org, ray_dir, iband, iquad);
+ 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;
+
+ /* Update the pixel accumulator of the current channel */
+ pix_accums[ichannel].sum_weights += weight;
+ pix_accums[ichannel].sum_weights_sqr += weight*weight;
+ pix_accums[ichannel].nweights += 1;
+
+ /* Update the pixel accumulator of per realisation time */
+ pix_accums[HTRDR_ESTIMATE_TIME].sum_weights += usec;
+ pix_accums[HTRDR_ESTIMATE_TIME].sum_weights_sqr += usec*usec;
+ pix_accums[HTRDR_ESTIMATE_TIME].nweights += 1;
+ }
+ }
+}
+
+static void
+draw_pixel_lw
+ (struct htrdr* htrdr,
+ const size_t ithread,
+ const size_t ipix[2],
+ const double pix_sz[2], /* Size of a pixel in the normalized image plane */
+ const struct htrdr_camera* cam,
+ const size_t spp,
+ struct ssp_rng* rng,
+ struct htrdr_accum* pix_accums)
+{
+ size_t isamp;
+ ASSERT(ipix && ipix && pix_sz && cam && rng && pix_accums);
+ ASSERT(htsky_is_long_wave(htrdr->sky));
+
+ /* Reset the pixel accumulators */
+ pix_accums[HTRDR_ESTIMATE_RADIANCE] = HTRDR_ACCUM_NULL;
+ pix_accums[HTRDR_ESTIMATE_TIME] = HTRDR_ACCUM_NULL;
+
+ FOR_EACH(isamp, 0, spp) {
+ struct time t0, t1;
+ double pix_samp[2];
+ double ray_org[3];
+ double ray_dir[3];
+ double weight;
+ double r0, r1;
+ size_t iband;
+ size_t iquad;
+ double usec;
+
+ /* Begin the registration of the time spent to in the realisation */
+ time_current(&t0);
+
+ /* Sample a position into the pixel, in the normalized image plane */
+ pix_samp[0] = ((double)ipix[0] + ssp_rng_canonical(rng)) * pix_sz[0];
+ pix_samp[1] = ((double)ipix[1] + ssp_rng_canonical(rng)) * pix_sz[1];
+
+ /* Generate a ray starting from the pinhole camera and passing through the
+ * pixel sample */
+ htrdr_camera_ray(cam, pix_samp, ray_org, ray_dir);
+
+ r0 = ssp_rng_canonical(rng);
+ r1 = ssp_rng_canonical(rng);
+
+ /* Sample a spectral band and a quadrature point */
+ iband = sample_lw_spectral_interval(htrdr, r0);
+ iquad = htsky_spectral_band_sample_quadrature(htrdr->sky, r1, iband);
+
+ /* Compute the luminance */
+ weight = htrdr_compute_radiance_lw
+ (htrdr, ithread, rng, ray_org, ray_dir, iband, iquad);
+ 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;
+
+ /* Update the pixel accumulator of the current channel */
+ pix_accums[HTRDR_ESTIMATE_RADIANCE].sum_weights += weight;
+ pix_accums[HTRDR_ESTIMATE_RADIANCE].sum_weights_sqr += weight*weight;
+ pix_accums[HTRDR_ESTIMATE_RADIANCE].nweights += 1;
+
+ /* Update the pixel accumulator of per realisation time */
+ pix_accums[HTRDR_ESTIMATE_TIME].sum_weights += usec;
+ pix_accums[HTRDR_ESTIMATE_TIME].sum_weights_sqr += usec*usec;
+ pix_accums[HTRDR_ESTIMATE_TIME].nweights += 1;
+ }
+}
+
static res_T
draw_tile
(struct htrdr* htrdr,
@@ -486,7 +650,6 @@ draw_tile
struct ssp_rng* rng,
struct tile* tile)
{
- size_t nchannels;
size_t npixels;
size_t mcode; /* Morton code of tile pixel */
ASSERT(htrdr && tile_org && tile_sz && pix_sz && cam && spp && tile);
@@ -495,14 +658,10 @@ draw_tile
npixels = round_up_pow2(MMAX(tile_sz[0], tile_sz[1]));
npixels *= npixels;
- /* Define how many channels to handle */
- nchannels = htsky_is_long_wave(htrdr->sky) ? 1 : 3;
-
FOR_EACH(mcode, 0, npixels) {
struct htrdr_accum* pix_accums;
size_t ipix_tile[2]; /* Pixel coord in the tile */
size_t ipix[2]; /* Pixel coord in the buffer */
- size_t ichannel;
ipix_tile[0] = morton2D_decode((uint32_t)(mcode>>0));
if(ipix_tile[0] >= tile_sz[0]) continue; /* Pixel is out of tile */
@@ -512,101 +671,15 @@ draw_tile
/* Fetch and reset the pixel accumulator */
pix_accums = tile_at(tile, ipix_tile[0], ipix_tile[1]);
- /* Reset the pixel accumulators */
- pix_accums[HTRDR_ESTIMATE_X] = HTRDR_ACCUM_NULL;
- pix_accums[HTRDR_ESTIMATE_Y] = HTRDR_ACCUM_NULL;
- pix_accums[HTRDR_ESTIMATE_Z] = HTRDR_ACCUM_NULL;
- pix_accums[HTRDR_ESTIMATE_TIME] = HTRDR_ACCUM_NULL;
-
/* Compute the pixel coordinate */
ipix[0] = tile_org[0] + ipix_tile[0];
ipix[1] = tile_org[1] + ipix_tile[1];
- FOR_EACH(ichannel, 0, nchannels) {
- /* Check that the X, Y and Z estimates are stored in accumulators 0, 1 et
- * 2, respectively */
- STATIC_ASSERT
- ( HTRDR_ESTIMATE_X == 0
- && HTRDR_ESTIMATE_Y == 1
- && HTRDR_ESTIMATE_Z == 2,
- Unexpected_htrdr_estimate_enumerate);
- size_t isamp;
-
- FOR_EACH(isamp, 0, spp) {
- struct time t0, t1;
- double pix_samp[2];
- double ray_org[3];
- double ray_dir[3];
- double weight;
- double r0, r1;
- size_t iband;
- size_t iquad;
- double usec;
-
- /* Begin the registration of the time spent to in the realisation */
- time_current(&t0);
-
- /* Sample a position into the pixel, in the normalized image plane */
- pix_samp[0] = ((double)ipix[0] + ssp_rng_canonical(rng)) * pix_sz[0];
- pix_samp[1] = ((double)ipix[1] + ssp_rng_canonical(rng)) * pix_sz[1];
-
- /* Generate a ray starting from the pinhole camera and passing through the
- * pixel sample */
- htrdr_camera_ray(cam, pix_samp, ray_org, ray_dir);
-
- r0 = ssp_rng_canonical(rng);
- r1 = ssp_rng_canonical(rng);
-
- if(htsky_is_long_wave(htrdr->sky)) {
- /* Sample a spectral band and a quadrature point */
- iband = sample_lw_spectral_interval(htrdr, r0);
- iquad = htsky_spectral_band_sample_quadrature(htrdr->sky, r1, iband);
- /* Compute the luminance */
- weight = htrdr_compute_radiance_lw
- (htrdr, ithread, rng, ray_org, ray_dir, iband, iquad);
- } else {
- /* Sample a spectral band and a quadrature point */
- switch(ichannel) {
- case 0:
- htsky_sample_sw_spectral_data_CIE_1931_X
- (htrdr->sky, r0, r1, &iband, &iquad);
- break;
- case 1:
- htsky_sample_sw_spectral_data_CIE_1931_Y
- (htrdr->sky, r0, r1, &iband, &iquad);
- break;
- case 2:
- htsky_sample_sw_spectral_data_CIE_1931_Z
- (htrdr->sky, r0, r1, &iband, &iquad);
- break;
- default: FATAL("Unreachable code.\n"); break;
- }
- /* Compute the luminance */
- weight = htrdr_compute_radiance_sw
- (htrdr, ithread, rng, ray_org, ray_dir, iband, iquad);
- }
- 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;
-
- /* Update the pixel accumulator of the current channel */
- pix_accums[ichannel].sum_weights += weight;
- pix_accums[ichannel].sum_weights_sqr += weight*weight;
- pix_accums[ichannel].nweights += 1;
-
- /* Update the pixel accumulator of per realisation time */
- pix_accums[HTRDR_ESTIMATE_TIME].sum_weights += usec;
- pix_accums[HTRDR_ESTIMATE_TIME].sum_weights_sqr += usec*usec;
- pix_accums[HTRDR_ESTIMATE_TIME].nweights += 1;
- }
- }
-
- /* Fill up the remaining channels with the estimate of the first one */
- while(ichannel < 3) {
- pix_accums[ichannel] = pix_accums[0];
- ++ichannel;
+ /* Draw the pixel */
+ if(htsky_is_long_wave(htrdr->sky)) {
+ draw_pixel_lw(htrdr, ithread, ipix, pix_sz, cam, spp, rng, pix_accums);
+ } else {
+ draw_pixel_sw(htrdr, ithread, ipix, pix_sz, cam, spp, rng, pix_accums);
}
}
return RES_OK;
