commit 266c5aed225c36bd514834979fbfca674433dacf
parent 30e8f59edd9b7d188f9bb14f369659abead4494d
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 29 May 2020 16:36:46 +0200
Fix the sun radiance and the MC weight in shortwave
Diffstat:
5 files changed, 59 insertions(+), 21 deletions(-)
diff --git a/src/htrdr_cie_xyz.c b/src/htrdr_cie_xyz.c
@@ -33,6 +33,9 @@ struct htrdr_cie_xyz {
struct darray_double cdf_X;
struct darray_double cdf_Y;
struct darray_double cdf_Z;
+ double rcp_integral_X;
+ double rcp_integral_Y;
+ double rcp_integral_Z;
double range[2]; /* Boundaries of the handled CIE XYZ color space */
double band_len; /* Length in nanometers of a band */
@@ -176,7 +179,7 @@ setup_cie_xyz
enum { X, Y, Z }; /* Helper constant */
double* pdf[3] = {NULL, NULL, NULL};
double* cdf[3] = {NULL, NULL, NULL};
- double sum[3] = {0,0,0};
+ double sum[3] = {0, 0, 0};
size_t i;
res_T res = RES_OK;
@@ -204,7 +207,6 @@ setup_cie_xyz
#undef SETUP_STIMULUS
/* Compute the *unormalized* pdf of the tristimulus */
-
FOR_EACH(i, 0, nbands) {
const double lambda_lo = cie->range[0] + (double)i * cie->band_len;
const double lambda_hi = MMIN(lambda_lo + cie->band_len, cie->range[1]);
@@ -218,6 +220,15 @@ setup_cie_xyz
sum[Y] += pdf[Y][i];
sum[Z] += pdf[Z][i];
}
+ #define CHK_SUM(Sum, Range, Fit) \
+ ASSERT(eq_eps(Sum, trapezoidal_integration(Range[0], Range[1], Fit), 1.e-3))
+ CHK_SUM(sum[X], cie->range, fit_x_bar_1931);
+ CHK_SUM(sum[Y], cie->range, fit_y_bar_1931);
+ CHK_SUM(sum[Z], cie->range, fit_z_bar_1931);
+ #undef CHK_SUM
+ cie->rcp_integral_X = 1.0 / sum[X];
+ cie->rcp_integral_Y = 1.0 / sum[Y];
+ cie->rcp_integral_Z = 1.0 / sum[Z];
FOR_EACH(i, 0, nbands) {
/* Normalize the pdf */
@@ -343,29 +354,38 @@ double
htrdr_cie_xyz_sample_X
(struct htrdr_cie_xyz* cie,
const double r0,
- const double r1)
+ const double r1,
+ double* pdf)
{
- return sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_X),
+ const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_X),
darray_double_size_get(&cie->cdf_X), fit_x_bar_1931, r0, r1);
+ if(pdf) *pdf = cie->rcp_integral_X;
+ return wlen;
}
double
htrdr_cie_xyz_sample_Y
(struct htrdr_cie_xyz* cie,
const double r0,
- const double r1)
+ const double r1,
+ double* pdf)
{
- return sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Y),
+ const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Y),
darray_double_size_get(&cie->cdf_Y), fit_y_bar_1931, r0, r1);
+ if(pdf) *pdf = cie->rcp_integral_Y;
+ return wlen;
}
double
htrdr_cie_xyz_sample_Z
(struct htrdr_cie_xyz* cie,
const double r0,
- const double r1)
+ const double r1,
+ double* pdf)
{
- return sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Z),
+ const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Z),
darray_double_size_get(&cie->cdf_Z), fit_z_bar_1931, r0, r1);
+ if(pdf) *pdf = cie->rcp_integral_Z;
+ return wlen;
}
diff --git a/src/htrdr_cie_xyz.h b/src/htrdr_cie_xyz.h
@@ -43,19 +43,22 @@ htrdr_cie_xyz_ref_put
extern LOCAL_SYM double
htrdr_cie_xyz_sample_X
(struct htrdr_cie_xyz* cie,
- const double r0, const double r1); /* Canonical numbers in [0, 1[ */
+ const double r0, const double r1, /* Canonical numbers in [0, 1[ */
+ double* pdf); /* In nm^-1. May be NULL */
/* Return a wavelength in nanometer */
extern LOCAL_SYM double
htrdr_cie_xyz_sample_Y
(struct htrdr_cie_xyz* cie,
- const double r0, const double r1); /* Canonical number in [0, 1[ */
+ const double r0, const double r1, /* Canonical number in [0, 1[ */
+ double* pdf); /* In nm^-1. May be NULL */
/* Return a wavelength in nanometer */
extern LOCAL_SYM double
htrdr_cie_xyz_sample_Z
(struct htrdr_cie_xyz* cie,
- const double r0, const double r1); /* Canonical number in [0, 1[ */
+ const double r0, const double r1, /* Canonical number in [0, 1[ */
+ double* pdf); /* In nm^-1. May be NULL */
#endif /* HTRDR_cie_xyz_H */
diff --git a/src/htrdr_draw_radiance.c b/src/htrdr_draw_radiance.c
@@ -514,6 +514,7 @@ draw_pixel_sw
double weight;
double r0, r1, r2;
double wlen; /* Sampled wavelength into the spectral band */
+ double pdf;
size_t iband; /* Sampled spectral band */
size_t iquad; /* Sampled quadrature point into the spectral band */
double usec;
@@ -535,20 +536,23 @@ draw_pixel_sw
/* Sample a spectral band and a quadrature point */
switch(ichannel) {
- case 0: wlen = htrdr_cie_xyz_sample_X(htrdr->cie, r0, r1); break;
- case 1: wlen = htrdr_cie_xyz_sample_Y(htrdr->cie, r0, r1); break;
- case 2: wlen = htrdr_cie_xyz_sample_Z(htrdr->cie, r0, r1); break;
+ case 0: wlen = htrdr_cie_xyz_sample_X(htrdr->cie, r0, r1, &pdf); break;
+ case 1: wlen = htrdr_cie_xyz_sample_Y(htrdr->cie, r0, r1, &pdf); break;
+ case 2: wlen = htrdr_cie_xyz_sample_Z(htrdr->cie, r0, r1, &pdf); break;
default: FATAL("Unreachable code.\n"); break;
}
iband = htsky_find_spectral_band(htrdr->sky, wlen);
iquad = htsky_spectral_band_sample_quadrature(htrdr->sky, r2, iband);
- /* Compute the luminance */
+ /* Compute the radiance in W/m^2/sr/m */
weight = htrdr_compute_radiance_sw
(htrdr, ithread, rng, ray_org, ray_dir, wlen, iband, iquad);
ASSERT(weight >= 0);
+ pdf *= 1.e9; /* Transform the pdf fro nm^-1 to m^-1 */
+ weight /= pdf; /* In W/m^2/sr */
+
/* End the registration of the per realisation time */
time_sub(&t0, time_current(&t1), &t0);
usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
diff --git a/src/htrdr_sun.c b/src/htrdr_sun.c
@@ -27,7 +27,7 @@
#include <star/ssp.h>
struct htrdr_sun {
- /* Short wave radiance in W.m^-2.sr^-1, for each spectral interval */
+ /* Short wave radiance in W/m^2/sr/m, for each spectral interval */
struct darray_double radiances_sw;
/* Short wave spectral interval boundaries, in cm^-1 */
@@ -37,6 +37,7 @@ struct htrdr_sun {
double cos_half_angle;
double solid_angle; /* In sr; solid_angle = 2*PI*(1 - cos(half_angle)) */
double frame[9];
+ double temperature; /* In K */
ref_T ref;
struct htrdr* htrdr;
@@ -87,13 +88,14 @@ htrdr_sun_create(struct htrdr* htrdr, struct htrdr_sun** out_sun)
if(!sun) {
htrdr_log_err(htrdr, "could not allocate the sun data structure.\n");
res = RES_MEM_ERR;
- goto error;
-}
+ goto error;
+ }
ref_init(&sun->ref);
sun->htrdr = htrdr;
darray_double_init(htrdr->allocator, &sun->radiances_sw);
darray_double_init(htrdr->allocator, &sun->wavenumbers_sw);
sun->half_angle = 4.6524e-3;
+ sun->temperature = 5778;
sun->cos_half_angle = cos(sun->half_angle);
sun->solid_angle = 2*PI*(1-sun->cos_half_angle);
d33_basis(sun->frame, main_dir);
@@ -112,8 +114,17 @@ htrdr_sun_create(struct htrdr* htrdr, struct htrdr_sun** out_sun)
}
FOR_EACH(i, 0, darray_double_size_get(&sun->radiances_sw)) {
- /* Convert the incoming flux in radiance */
- darray_double_data_get(&sun->radiances_sw)[i] = incoming_flux_sw[i] / PI;
+ double band_bounds_m[2];
+ double band_len_m;
+
+ /* Retrieve the boundaries of the spectral band (in meters) */
+ band_bounds_m[0] = wavenumber_to_wavelength(wavenumbers_sw[i+1])*1.e-9;
+ band_bounds_m[1] = wavenumber_to_wavelength(wavenumbers_sw[i+0])*1.e-9;
+ band_len_m = band_bounds_m[1] - band_bounds_m[0];
+
+ /* Convert the incoming flux in radiance (in W/m^2/sr/m) */
+ darray_double_data_get(&sun->radiances_sw)[i] =
+ incoming_flux_sw[i] / (band_len_m * sun->solid_angle);
}
FOR_EACH(i, 0, darray_double_size_get(&sun->wavenumbers_sw)) {
darray_double_data_get(&sun->wavenumbers_sw)[i] = wavenumbers_sw[i];
diff --git a/src/htrdr_sun.h b/src/htrdr_sun.h
@@ -54,7 +54,7 @@ extern LOCAL_SYM double
htrdr_sun_get_solid_angle
(const struct htrdr_sun* sun);
-extern LOCAL_SYM double /* W.sr^-1.m^-2 */
+extern LOCAL_SYM double /* W/m^2/sr/m */
htrdr_sun_get_radiance
(const struct htrdr_sun* sun,
const double wavelength);
