commit 5cfeac159892a245957b0f9d16d47a6714e87f07
parent 2172c05c7483399eb5e732ca18e8304213b14773
Author: Arfaux Anthony <anthony.arfaux@obspm.fr>
Date: Fri, 16 May 2025 14:43:44 +0200
Improve volumic radiative budget
Reduce variance by separating sampling of directions within the sun cone
from sampling of directions outside it.
Diffstat:
6 files changed, 180 insertions(+), 10 deletions(-)
diff --git a/Makefile b/Makefile
@@ -28,13 +28,13 @@ include config.mk
default install uninstall lint clean:
@$(MAKE) -fMakefile.core $@
- @if [ "$(ATMOSPHERE)" == "ENABLE" ]; then $(MAKE) -fMakefile.atmosphere $@; fi
- @if [ "$(COMBUSTION)" == "ENABLE" ]; then $(MAKE) -fMakefile.combustion $@; fi
- @if [ "$(PLANETS)" == ENABLE ]; then $(MAKE) -fMakefile.planets $@; fi
+ @if [ "$(ATMOSPHERE)" = "ENABLE" ]; then $(MAKE) -fMakefile.atmosphere $@; fi
+ @if [ "$(COMBUSTION)" = "ENABLE" ]; then $(MAKE) -fMakefile.combustion $@; fi
+ @if [ "$(PLANETS)" = ENABLE ]; then $(MAKE) -fMakefile.planets $@; fi
test:
- @if [ "$(COMBUSTION)" == "ENABLE" ]; then $(MAKE) -fMakefile.combustion $@; fi
- @if [ "$(PLANETS)" == ENABLE ]; then $(MAKE) -fMakefile.planets $@; fi
+ @if [ "$(COMBUSTION)" = "ENABLE" ]; then $(MAKE) -fMakefile.combustion $@; fi
+ @if [ "$(PLANETS)" = ENABLE ]; then $(MAKE) -fMakefile.planets $@; fi
default: htrdr
install: install_common
diff --git a/Makefile.core b/Makefile.core
@@ -147,7 +147,7 @@ htrdr-core.pc: config.mk htrdr-core.pc.in
# Miscellaneous
################################################################################
install:
- @if [ "$(LIB_TYPE)" == "SHARED" ]; then \
+ @if [ "$(LIB_TYPE)" = "SHARED" ]; then \
$(SHELL) install.sh 755 "$(DESTDIR)$(LIBPREFIX)" $(LIBNAME); \
fi
diff --git a/src/planets/htrdr_planets_compute_radiance.c b/src/planets/htrdr_planets_compute_radiance.c
@@ -504,8 +504,6 @@ volume_scattering
ssf_phase_sample(phase, args->rng, wo, sc_dir, NULL);
- /* Sample a direction toward the source */
- pdf = htrdr_planets_source_sample_direction(cmd->source, args->rng, sc_pos, wi);
/* In short wave, manage the contribution of the external source */
switch(cmd->spectral_domain.type) {
@@ -515,6 +513,9 @@ volume_scattering
case HTRDR_SPECTRAL_SW:
case HTRDR_SPECTRAL_SW_CIE_XYZ:
+
+ /* Sample a direction toward the source */
+ pdf = htrdr_planets_source_sample_direction(cmd->source, args->rng, sc_pos, wi);
/* Calculate the direct contribution at the scattering position */
Ld = direct_contribution(cmd, args, sc_pos, wi, NULL);
rho = ssf_phase_eval(phase, wo, wi);
diff --git a/src/planets/htrdr_planets_solve_volrad_budget.c b/src/planets/htrdr_planets_solve_volrad_budget.c
@@ -237,6 +237,126 @@ realisation
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);
+
+ 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);
+
+ /* 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};
+
+ 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;
+
+ /* Preconditions */
+ ASSERT(cmd && args && volrad_mesh_desc);
+ ASSERT(cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_VOLUMIC_RADIATIVE_BUDGET);
+ ASSERT(source);
+
+ 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);
+
+ /* 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 void
solve_volumic_radiative_budget
(struct htrdr* htrdr,
@@ -276,8 +396,16 @@ solve_volumic_radiative_budget
/* Start of realisation time recording */
time_current(&t0);
- /* Run the realisation */
- w = realisation(cmd, args, &volrad_mesh_desc);
+ 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);
+ }
/* Stop time recording */
time_sub(&t0, time_current(&t1), &t0);
diff --git a/src/planets/htrdr_planets_source.c b/src/planets/htrdr_planets_source.c
@@ -296,6 +296,40 @@ 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,6 +64,13 @@ 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,
