commit 62babff860e9dd45960eeb023ec255d8e10ec3d9
parent 5cfeac159892a245957b0f9d16d47a6714e87f07
Author: Arfaux Anthony <anthony.arfaux@obspm.fr>
Date: Fri, 16 May 2025 12:06:34 +0200
Decomposition of the radiative budget
The realization function for the radiative budget has been modified in
such a way that its direct and diffuse components have been
differenciated. The direct component of the radiative budget only has to
sample the solar cone, while the diffuse component still samples the
unit sphere. This will reduce the variance associated to the total in
the case the external source is restricted to a very small solid angle.
Diffstat:
5 files changed, 78 insertions(+), 182 deletions(-)
diff --git a/src/planets/htrdr_planets_c.h b/src/planets/htrdr_planets_c.h
@@ -77,6 +77,15 @@ struct planets_voxel_radiative_budget {
HTRDR_ACCUM_NULL__ \
}
+enum planets_radiance_cpnt_flag {
+ PLANETS_RADIANCE_CPNT_DIRECT = BIT(0),
+ PLANETS_RADIANCE_CPNT_DIFFUSE = BIT(1),
+ PLANETS_RADIANCE_CPNT_NONE = 0,
+ PLANETS_RADIANCE_CPNT_ALL =
+ PLANETS_RADIANCE_CPNT_DIRECT
+ | PLANETS_RADIANCE_CPNT_DIFFUSE
+};
+
struct planets_compute_radiance_args {
struct ssp_rng* rng;
size_t ithread; /* Index of the thread executing the function */
@@ -87,8 +96,11 @@ struct planets_compute_radiance_args {
double wlen; /* In nm */
size_t iband; /* Spectral band index */
size_t iquad; /* Quadrature point */
+
+ int component; /* Combination of planets_radiance_cpnt_flag */
};
-#define PLANETS_COMPUTE_RADIANCE_ARGS_NULL__ {NULL, 0, {0,0,0}, {0,0,0}, 0, 0, 0}
+#define PLANETS_COMPUTE_RADIANCE_ARGS_NULL__ \
+ {NULL, 0, {0,0,0}, {0,0,0}, 0, 0, 0, PLANETS_RADIANCE_CPNT_ALL}
static const struct planets_compute_radiance_args
PLANETS_COMPUTE_RADIANCE_ARGS_NULL = PLANETS_COMPUTE_RADIANCE_ARGS_NULL__;
diff --git a/src/planets/htrdr_planets_compute_radiance.c b/src/planets/htrdr_planets_compute_radiance.c
@@ -504,7 +504,6 @@ volume_scattering
ssf_phase_sample(phase, args->rng, wo, sc_dir, NULL);
-
/* In short wave, manage the contribution of the external source */
switch(cmd->spectral_domain.type) {
case HTRDR_SPECTRAL_LW:
@@ -604,16 +603,34 @@ planets_compute_radiance
double L = 0; /* Radiance in W/m²/sr/m */
size_t nsc = 0; /* Number of surface or volume scatterings (for debug) */
int longwave = 0;
+ int shortwave = 0;
+ int direct = 0;
+ int diffuse = 0;
ASSERT(cmd && check_planets_compute_radiance_args(cmd, args) == RES_OK);
d3_set(pos, args->path_org);
d3_set(dir, args->path_dir);
+
longwave = cmd->spectral_domain.type == HTRDR_SPECTRAL_LW;
+ shortwave = !longwave;
- if(!longwave && htrdr_planets_source_is_targeted(cmd->source, pos, dir)) {
+ /* In shortwave define which components are enabled */
+ if(shortwave) {
+ direct = (args->component & PLANETS_RADIANCE_CPNT_DIRECT) != 0;
+ diffuse = (args->component & PLANETS_RADIANCE_CPNT_DIFFUSE) != 0;
+ }
+
+ /* Handle direct shortwave contribution */
+ if(shortwave
+ && direct
+ && htrdr_planets_source_is_targeted(cmd->source, pos, dir)) {
L = direct_contribution(cmd, args, pos, dir, NULL); /* In W/m²/sr/m */
}
+ /* Nothing left to do: if only the diffuse component is required
+ * in the SW, the diffuse component should not be computed */
+ if(shortwave && !diffuse) return L;
+
for(;;) {
double ev_pos[3];
double sc_dir[3];
diff --git a/src/planets/htrdr_planets_solve_volrad_budget.c b/src/planets/htrdr_planets_solve_volrad_budget.c
@@ -35,6 +35,13 @@
#include <rsys/clock_time.h>
+enum weight_type {
+ DIRECT, /* 0 */
+ DIFFUSE, /* 1 */
+ TOTAL, /* 2 */
+ WEIGHTS_COUNT /* 3 */
+};
+
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -183,11 +190,12 @@ get_ka
return ka;
}
-static double
+static void
realisation
(struct htrdr_planets* cmd,
const struct htrdr_solve_item_args* args,
- const struct smsh_desc* volrad_mesh_desc)
+ const struct smsh_desc* volrad_mesh_desc,
+ double weights[WEIGHTS_COUNT]) /* [W/m^3] */
{
struct planets_compute_radiance_args rad_args =
PLANETS_COMPUTE_RADIANCE_ARGS_NULL;
@@ -200,161 +208,69 @@ realisation
size_t iquad = 0; /* Quadrature point */
/* Spatial & angular integration */
+ double dir_src[3] = {0,0,0}; /* Direction toward the source */
double dir[3] = {0,0,0};
double pos[3] = {0,0,0};
+ double dir_src_pdf = 0;
+ double dir_pdf = 0;
double S = 0; /* Source [W/m^2/sr/m] */
- double L = 0; /* Radiance [W/m^2/sr/m] */
+ double L_direct = 0; /* Direct radiance [W/m^2/sr/m] */
+ double L_diffuse = 0; /* Diffuse radiance [W/m^2/sr/m] */
double ka = 0; /* Absorption coefficient */
- /* Monte Carlo weight */
- double weight = 0; /* weight [W/m^3] */
-
/* Preconditions */
ASSERT(cmd && args && volrad_mesh_desc);
ASSERT(cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_VOLUMIC_RADIATIVE_BUDGET);
- spectral_sampling(cmd, args, &wlen, &wlen_pdf_nm, &iband, &iquad);
- position_sampling(args, volrad_mesh_desc, pos);
- ssp_ran_sphere_uniform(args->rng, dir, NULL/*pdf*/);
-
- /* Compute the radiance in W/m^2/sr/m */
- d3_set(rad_args.path_org, pos);
- d3_set(rad_args.path_dir, dir);
- rad_args.rng = args->rng;
- rad_args.ithread = args->ithread;
- rad_args.wlen = wlen; /* [nm] */
- rad_args.iband = iband;
- rad_args.iquad = iquad;
- L = planets_compute_radiance(cmd, &rad_args); /* [W/m^2/sr/m] */
-
- S = get_source(cmd, pos, wlen); /* [W/m^2/sr/m] */
-
- ka = get_ka(cmd, pos, iband, iquad);
- wlen_pdf_m = wlen_pdf_nm * 1.e9; /* Transform pdf from nm^-1 to m^-1 */
-
- weight = ka * (L - S) * 4*PI / wlen_pdf_m; /* [W/m^3] */
- return weight;
-}
-
-static double
-realisation_cone
- (struct htrdr_planets* cmd,
- const struct htrdr_solve_item_args* args,
- const struct smsh_desc* volrad_mesh_desc)
-{
- struct planets_compute_radiance_args rad_args =
- PLANETS_COMPUTE_RADIANCE_ARGS_NULL;
-
- /* Spectral integration */
- double wlen = 0; /* Wavelength [nm] */
- double wlen_pdf_nm = 0; /* Wavelength pdf [nm^-1] */
- double wlen_pdf_m = 0; /* Wavelength pdf [m^-1] */
- size_t iband = 0; /* Spectral band */
- size_t iquad = 0; /* Quadrature point */
-
- /* Spatial & angular integration */
- double dir[3] = {0,0,0};
- double pos[3] = {0,0,0};
-
- double S = 0; /* Source [W/m^2/sr/m] */
- double L = 0; /* Radiance [W/m^2/sr/m] */
- double ka = 0; /* Absorption coefficient */
-
- /* Monte Carlo weight */
- double weight = 0; /* weight [W/m^3] */
-
- /* For sampling in cone */
- const struct htrdr_planets_source* source;
- double dir_pdf=0;
- source = cmd->source;
-
- /* Preconditions */
- ASSERT(cmd && args && volrad_mesh_desc);
- ASSERT(cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_VOLUMIC_RADIATIVE_BUDGET);
- ASSERT(source);
+ /* Initialise the weights */
+ memset(weights, 0, sizeof(double)*WEIGHTS_COUNT);
spectral_sampling(cmd, args, &wlen, &wlen_pdf_nm, &iband, &iquad);
position_sampling(args, volrad_mesh_desc, pos);
- dir_pdf = htrdr_planets_source_sample_direction(source, args->rng, pos, dir);
+ ssp_ran_sphere_uniform(args->rng, dir, &dir_pdf);
/* Compute the radiance in W/m^2/sr/m */
d3_set(rad_args.path_org, pos);
- d3_set(rad_args.path_dir, dir);
rad_args.rng = args->rng;
rad_args.ithread = args->ithread;
rad_args.wlen = wlen; /* [nm] */
rad_args.iband = iband;
rad_args.iquad = iquad;
- L = planets_compute_radiance(cmd, &rad_args); /* [W/m^2/sr/m] */
-
- S = get_source(cmd, pos, wlen); /* [W/m^2/sr/m] */
-
- ka = get_ka(cmd, pos, iband, iquad);
- wlen_pdf_m = wlen_pdf_nm * 1.e9; /* Transform pdf from nm^-1 to m^-1 */
-
- weight = ka * (L - S) / (wlen_pdf_m * dir_pdf); /* [W/m^3] */
- return weight;
-}
-
-static double
-realisation_out_of_cone
- (struct htrdr_planets* cmd,
- const struct htrdr_solve_item_args* args,
- const struct smsh_desc* volrad_mesh_desc)
-{
- struct planets_compute_radiance_args rad_args =
- PLANETS_COMPUTE_RADIANCE_ARGS_NULL;
-
- /* Spectral integration */
- double wlen = 0; /* Wavelength [nm] */
- double wlen_pdf_nm = 0; /* Wavelength pdf [nm^-1] */
- double wlen_pdf_m = 0; /* Wavelength pdf [m^-1] */
- size_t iband = 0; /* Spectral band */
- size_t iquad = 0; /* Quadrature point */
- /* Spatial & angular integration */
- double dir[3] = {0,0,0};
- double pos[3] = {0,0,0};
+ if(cmd->spectral_domain.type == HTRDR_SPECTRAL_LW) {
+ /* In the longwave (radiation due to the medium), simply sample a radiative
+ * path for the sampled direction and position: the radiance is considered
+ * as purely diffuse. */
+ d3_set(rad_args.path_dir, dir);
+ L_direct = 0.0;
+ L_diffuse = planets_compute_radiance(cmd, &rad_args); /* [W/m^2/sr/m] */
- double S = 0; /* Source [W/m^2/sr/m] */
- double L = 0; /* Radiance [W/m^2/sr/m] */
- double ka = 0; /* Absorption coefficient */
-
- /* Monte Carlo weight */
- double weight = 0; /* weight [W/m^3] */
-
- /* For sampling out of cone cone */
- const struct htrdr_planets_source* source;
- double dir_pdf=0;
- source = cmd->source;
+ } else {
+ /* In the so-called shortwave region (actually, the radiation due the
+ * external source) is decomposed in its direct and diffuse components */
- /* Preconditions */
- ASSERT(cmd && args && volrad_mesh_desc);
- ASSERT(cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_VOLUMIC_RADIATIVE_BUDGET);
- ASSERT(source);
+ dir_src_pdf = htrdr_planets_source_sample_direction
+ (cmd->source, args->rng, pos, dir_src);
- spectral_sampling(cmd, args, &wlen, &wlen_pdf_nm, &iband, &iquad);
- position_sampling(args, volrad_mesh_desc, pos);
- dir_pdf = htrdr_planets_source_sample_direction_out(source, args->rng, pos, dir);
+ d3_set(rad_args.path_dir, dir_src);
+ rad_args.component = PLANETS_RADIANCE_CPNT_DIRECT;
+ L_direct = planets_compute_radiance(cmd, &rad_args); /* [W/m^2/sr/m] */
- /* Compute the radiance in W/m^2/sr/m */
- d3_set(rad_args.path_org, pos);
- d3_set(rad_args.path_dir, dir);
- rad_args.rng = args->rng;
- rad_args.ithread = args->ithread;
- rad_args.wlen = wlen; /* [nm] */
- rad_args.iband = iband;
- rad_args.iquad = iquad;
- L = planets_compute_radiance(cmd, &rad_args); /* [W/m^2/sr/m] */
+ d3_set(rad_args.path_dir, dir);
+ rad_args.component = PLANETS_RADIANCE_CPNT_DIFFUSE;
+ L_diffuse = planets_compute_radiance(cmd, &rad_args); /* [W/m^2/sr/m] */
+ }
S = get_source(cmd, pos, wlen); /* [W/m^2/sr/m] */
ka = get_ka(cmd, pos, iband, iquad);
wlen_pdf_m = wlen_pdf_nm * 1.e9; /* Transform pdf from nm^-1 to m^-1 */
- weight = ka * (L - S) / (wlen_pdf_m * dir_pdf); /* [W/m^3] */
- return weight;
+ /* Calculate the weights [W/m^3] */
+ weights[DIRECT] = ka * (L_direct - S) / (wlen_pdf_m * dir_src_pdf);
+ weights[DIFFUSE] = ka * (L_diffuse - S) / (wlen_pdf_m * dir_pdf);
+ weights[TOTAL] = weights[DIRECT] + weights[DIFFUSE];
}
static void
@@ -390,30 +306,22 @@ solve_volumic_radiative_budget
struct time t0, t1;
/* Monte Carlo weights */
- double w = 0; /* [W/m^3] */
+ double w[WEIGHTS_COUNT] = {0}; /* [W/m^3] */
double usec = 0; /* [us] */
/* Start of realisation time recording */
time_current(&t0);
- if (cmd->spectral_domain.type == HTRDR_SPECTRAL_LW) {
- /* Run the realisation */
- w = realisation(cmd, args, &volrad_mesh_desc);
- }
- else{
- /* Run the realisation within source cone */
- w = realisation_cone(cmd, args, &volrad_mesh_desc);
- /* Run the realisation out of source cone */
- w += realisation_out_of_cone(cmd, args, &volrad_mesh_desc);
- }
+ /* Run the realisation */
+ realisation(cmd, args, &volrad_mesh_desc, w);
/* Stop time recording */
time_sub(&t0, time_current(&t1), &t0);
usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
/* Update the volumic radiance budget accumulator */
- budget.sum_weights += w;
- budget.sum_weights_sqr += w*w;
+ budget.sum_weights += w[TOTAL];
+ budget.sum_weights_sqr += w[TOTAL]*w[TOTAL];
budget.nweights += 1;
/* Update the per realisation time accumulator */
diff --git a/src/planets/htrdr_planets_source.c b/src/planets/htrdr_planets_source.c
@@ -296,40 +296,6 @@ htrdr_planets_source_sample_direction
return pdf;
}
-double
-htrdr_planets_source_sample_direction_out
- (const struct htrdr_planets_source* source,
- struct ssp_rng* rng,
- const double pos[3],
- double dir[3])
-{
- double main_dir[3];
- double half_angle; /* In radians */
- double cos_half_angle;
- double dst; /* In m */
- double pdf;
- ASSERT(source && rng && pos && dir);
-
- /* compute the direction of `pos' toward the center of the source */
- d3_sub(main_dir, source->position, pos);
-
- /* Normalize the direction and keep the distance from `pos' to the center of
- * the source */
- dst = d3_normalize(main_dir, main_dir);
- CHK(dst > source->radius);
-
- /* Sample the source according to its solid angle,
- * i.e. 2*PI*(1 - cos(half_angle)) */
- half_angle = asin(source->radius/dst);
- cos_half_angle = - cos(half_angle);
- main_dir[0] = -main_dir[0];
- main_dir[1] = -main_dir[1];
- main_dir[2] = -main_dir[2];
- ssp_ran_sphere_cap_uniform(rng, main_dir, cos_half_angle, dir, &pdf);
-
- return pdf;
-}
-
double /* In W/m²/sr/m */
htrdr_planets_source_get_radiance
(const struct htrdr_planets_source* source,
diff --git a/src/planets/htrdr_planets_source.h b/src/planets/htrdr_planets_source.h
@@ -64,13 +64,6 @@ htrdr_planets_source_sample_direction
const double pos[3], /* Position from which direction is sampled */
double dir[3]);
-extern LOCAL_SYM double
-htrdr_planets_source_sample_direction_out
- (const struct htrdr_planets_source* source,
- struct ssp_rng* rng,
- const double pos[3], /* Position from which direction is sampled */
- double dir[3]);
-
extern LOCAL_SYM double /* In W/m²/sr/m */
htrdr_planets_source_get_radiance
(const struct htrdr_planets_source* source,
