commit 768970bf7a382cd98e41f2374fa3dac27044883d
parent d978bcea7156dda07113a67043ccf6292ba06931
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 21 Nov 2022 15:12:11 +0100
htrdr-planeto: first write of a complete lw algorithm
Diffstat:
3 files changed, 311 insertions(+), 45 deletions(-)
diff --git a/src/planeto/htrdr_planeto_c.h b/src/planeto/htrdr_planeto_c.h
@@ -113,7 +113,13 @@ planeto_get_pixel_format
/* Return the radiance in W/m²/sr/m */
extern LOCAL_SYM double
-planeto_compute_radiance
+planeto_compute_radiance_sw
+ (struct htrdr_planeto* cmd,
+ const struct planeto_compute_radiance_args* args);
+
+/* Return the radiance in W/m²/sr/m */
+extern LOCAL_SYM double
+planeto_compute_radiance_lw
(struct htrdr_planeto* cmd,
const struct planeto_compute_radiance_args* args);
diff --git a/src/planeto/htrdr_planeto_compute_radiance.c b/src/planeto/htrdr_planeto_compute_radiance.c
@@ -32,25 +32,25 @@
/* Syntactic sugar */
#define DISTANCE_NONE(Dst) ((Dst) >= FLT_MAX)
-#define SURFACE_SCATTERING(Scattering) (!S3D_HIT_NONE(&(Scattering)->hit))
+#define SURFACE_EVENT(Scattering) (!S3D_HIT_NONE(&(Scattering)->hit))
-struct scattering {
- /* Set to S3D_HIT_NULL if the scattering occurs in volume.*/
+struct event {
+ /* Set to S3D_HIT_NULL if the event occurs in volume.*/
struct s3d_hit hit;
- /* The surface normal is defined only if scattering is on the surface. It is
+ /* The surface normal is defined only if event is on the surface. It is
* normalized and looks towards the incoming direction */
double normal[3];
- /* Cells in which the scattering position is located */
+ /* Cells in which the event position is located */
struct rnatm_cell_pos cells[RNATM_MAX_COMPONENTS_COUNT];
double distance; /* Distance from ray origin to scattering position */
};
-#define SCATTERING_NULL__ { \
+#define EVENT_NULL__ { \
S3D_HIT_NULL__, {0,0,0}, {RNATM_CELL_POS_NULL__}, DBL_MAX \
}
-static const struct scattering SCATTERING_NULL = SCATTERING_NULL__;
+static const struct event EVENT_NULL = EVENT_NULL__;
/* Arguments of the filtering function used to sample a position */
struct sample_distance_context {
@@ -294,21 +294,22 @@ direct_contribution
}
static void
-sample_scattering
+sample_event
(struct htrdr_planeto* cmd,
const struct planeto_compute_radiance_args* args,
+ const enum rnatm_radcoef radcoef,
const double pos[3],
const double dir[3],
const struct s3d_hit* hit_from,
- struct scattering* sc)
+ struct event* evt)
{
struct rngrd_trace_ray_args rt = RNGRD_TRACE_RAY_ARGS_DEFAULT;
struct s3d_hit hit;
double range[2];
double dst;
- ASSERT(cmd && args && pos && dir && hit_from && sc);
+ ASSERT(cmd && args && pos && dir && hit_from && evt);
- *sc = SCATTERING_NULL;
+ *evt = EVENT_NULL;
/* Look for a surface intersection */
d3_set(rt.ray_org, pos);
@@ -320,26 +321,26 @@ sample_scattering
/* Look for an atmospheric collision */
range[0] = 0;
range[1] = hit.distance;
- dst = sample_distance(cmd, args, sc->cells, RNATM_RADCOEF_Ks, pos, dir, range);
+ dst = sample_distance(cmd, args, evt->cells, radcoef, pos, dir, range);
- /* Scattering in volume */
+ /* Event occurs in volume */
if(!DISTANCE_NONE(dst)) {
- sc->distance = dst;
- sc->hit = S3D_HIT_NULL;
+ evt->distance = dst;
+ evt->hit = S3D_HIT_NULL;
- /* Surface scattering */
+ /* Event is on surface */
} else if(!S3D_HIT_NONE(&hit)) {
/* Normalize the normal and ensure that it points to `dir' */
- d3_normalize(sc->normal, d3_set_f3(sc->normal, hit.normal));
- if(d3_dot(sc->normal, dir) > 0) d3_minus(sc->normal, sc->normal);
+ d3_normalize(evt->normal, d3_set_f3(evt->normal, hit.normal));
+ if(d3_dot(evt->normal, dir) > 0) d3_minus(evt->normal, evt->normal);
- sc->distance = hit.distance;
- sc->hit = hit;
+ evt->distance = hit.distance;
+ evt->hit = hit;
- /* No scattering */
+ /* No event */
} else {
- sc->distance = INF;
- sc->hit = S3D_HIT_NULL;
+ evt->distance = INF;
+ evt->hit = S3D_HIT_NULL;
}
}
@@ -400,12 +401,13 @@ create_phase_fn
return phase;
}
-/* Return the direct contribution at the scattering position */
+/* In shortwave, return the contribution of the external source at the bounce
+ * position. In longwave, simply return 0 */
static double
-surface_scattering
+surface_bounce
(struct htrdr_planeto* cmd,
const struct planeto_compute_radiance_args* args,
- const struct scattering* sc,
+ const struct event* sc,
const double sc_pos[3], /* Scattering position */
const double in_dir[3], /* Incident direction */
double sc_dir[3], /* Sampled scattering direction */
@@ -431,6 +433,10 @@ surface_scattering
/* Fully absorbs transmissions */
if(mask & SSF_TRANSMISSION) reflectivity = 0;
+ /* No external source in longwave */
+ if(cmd->spectral_domain.spectral_type == HTRDR_SPECTRAL_LW)
+ goto exit;
+
/* Calculate direct contribution for specular reflection */
if((mask & SSF_SPECULAR)
&& (mask & SSF_REFLECTION)
@@ -460,18 +466,19 @@ surface_scattering
}
}
+exit:
SSF(bsdf_ref_put(bsdf));
-
*out_refl = reflectivity;
return L;
}
-/* Return the direct contribution at the scattering position */
+/* In shortwave, return the contribution at the scattering position of the
+ * external source. Returns 0 in long wave */
static INLINE double
volume_scattering
(struct htrdr_planeto* cmd,
const struct planeto_compute_radiance_args* args,
- const struct scattering* sc,
+ const struct event* sc,
const double sc_pos[3], /* Scattering position */
const double in_dir[3], /* Incident direction */
double sc_dir[3]) /* Sampled scattering direction */
@@ -493,20 +500,92 @@ volume_scattering
/* Sample a direction toward the source */
pdf = htrdr_planeto_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);
- L = Ld * rho / pdf;
+ /* In short wave, manage the contribution of the external source */
+ switch(cmd->spectral_domain.spectral_type) {
+ case HTRDR_SPECTRAL_LW:
+ /* Nothing to be done */
+ break;
+
+ case HTRDR_SPECTRAL_SW:
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ /* Calculate the direct contribution at the scattering position */
+ Ld = direct_contribution(cmd, args, sc_pos, wi, NULL);
+ rho = ssf_phase_eval(phase, wo, wi);
+ L = Ld * rho / pdf;
+ break;
+
+ default: FATAL("Unreachable code.\n"); break;
+ }
SSF(phase_ref_put(phase));
return L;
}
+static INLINE enum rnatm_radcoef
+sample_volume_event_type
+ (const struct htrdr_planeto* cmd,
+ const struct planeto_compute_radiance_args* args,
+ struct event* evt)
+{
+ struct rnatm_get_radcoef_args get_k_args = RNATM_GET_RADCOEF_ARGS_NULL;
+ double ka, kext;
+ double r;
+ ASSERT(cmd && args && evt);
+
+ get_k_args.cells = evt->cells;
+ get_k_args.iband = args->iband;
+ get_k_args.iquad = args->iquad;
+
+ /* Retrieve the absorption coefficient */
+ get_k_args.radcoef = RNATM_RADCOEF_Ka;
+ RNATM(get_radcoef(cmd->atmosphere, &get_k_args, &ka));
+
+ /* Retrieve the extinction coefficient */
+ get_k_args.radcoef = RNATM_RADCOEF_Kext;
+ RNATM(get_radcoef(cmd->atmosphere, &get_k_args, &kext));
+
+ r = ssp_rng_canonical(args->rng);
+ if(r < ka / kext) {
+ return RNATM_RADCOEF_Ka; /* Absorption */
+ } else {
+ return RNATM_RADCOEF_Ks; /* Scattering */
+ }
+}
+
+static INLINE double
+get_temperature
+ (const struct htrdr_planeto* cmd,
+ const struct event* evt)
+{
+ double T = 0;
+ ASSERT(cmd && evt);
+
+ if(!SURFACE_EVENT(evt)) {
+ const struct rnatm_cell_pos* cell = NULL;
+
+ /* Get gas temperature */
+ cell = &RNATM_GET_COMPONENT_CELL(evt->cells, RNATM_GAS);
+ RNATM(cell_get_gas_temperature(cmd->atmosphere, cell, &T));
+
+ } else {
+ struct rngrd_get_temperature_args temp_args = RNGRD_GET_TEMPERATURE_ARGS_NULL;
+
+ /* Get ground temperature */
+ temp_args.prim = evt->hit.prim;
+ temp_args.barycentric_coords[0] = evt->hit.uv[0];
+ temp_args.barycentric_coords[1] = evt->hit.uv[1];
+ temp_args.barycentric_coords[2] = 1 - evt->hit.uv[0] - evt->hit.uv[1];
+ RNGRD(get_temperature(cmd->ground, &temp_args, &T));
+
+ }
+ return T;
+}
+
/*******************************************************************************
* Local functions
******************************************************************************/
-double
-planeto_compute_radiance
+double /* Radiance in W/m²/sr/m */
+planeto_compute_radiance_sw
(struct htrdr_planeto* cmd,
const struct planeto_compute_radiance_args* args)
{
@@ -526,11 +605,12 @@ planeto_compute_radiance
}
for(;;) {
- struct scattering sc;
+ struct event sc;
double sc_pos[3];
double sc_dir[3];
- sample_scattering(cmd, args, pos, dir, &hit_from, &sc);
+ /* Sample a scattering event */
+ sample_event(cmd, args, RNATM_RADCOEF_Ks, pos, dir, &hit_from, &sc);
/* No scattering. Stop the path */
if(DISTANCE_NONE(sc.distance)) break;
@@ -546,7 +626,7 @@ planeto_compute_radiance
sc_pos[2] = pos[2] + dir[2] * sc.distance;
/* Scattering in volume */
- if(!SURFACE_SCATTERING(&sc)) {
+ if(!SURFACE_EVENT(&sc)) {
double Ls; /* Direct contribution at the scattering position */
Ls = volume_scattering(cmd, args, &sc, sc_pos, dir, sc_dir);
L += Tr_abs * Ls;
@@ -554,11 +634,11 @@ planeto_compute_radiance
/* Reset surface intersection */
hit_from = S3D_HIT_NULL;
- /* Surface scattering */
+ /* Surface bounce */
} else {
- double Ls; /* Direct contribution at the scattering position */
+ double Ls; /* Direct contribution at the bounce position */
double refl; /* Surface reflectivity */
- Ls = surface_scattering(cmd, args, &sc, sc_pos, dir, sc_dir, &refl);
+ Ls = surface_bounce(cmd, args, &sc, sc_pos, dir, sc_dir, &refl);
L += Tr_abs * Ls;
/* Russian roulette wrt surface reflectivity */
@@ -577,3 +657,68 @@ planeto_compute_radiance
return L;
}
+
+double /* Radiance in W/m²/sr/m */
+planeto_compute_radiance_lw
+ (struct htrdr_planeto* cmd,
+ const struct planeto_compute_radiance_args* args)
+{
+ struct s3d_hit hit_from = S3D_HIT_NULL;
+ struct event evt;
+ double pos[3];
+ double dir[3];
+ double wlen_m = 0; /* Wavelength in meters */
+ double T = 0; /* Temperature in K */
+ double L = 0; /* Radiance in W/m²/sr/m */
+ ASSERT(cmd && check_planeto_compute_radiance_args(cmd, args) == RES_OK);
+
+ d3_set(pos, args->path_org);
+ d3_set(dir, args->path_dir);
+
+ for(;;) {
+ double next_pos[3];
+ double sc_dir[3];
+
+ sample_event(cmd, args, RNATM_RADCOEF_Kext, pos, dir, &hit_from, &evt);
+
+ /* No collision on surface or in volume. Stop the path */
+ if(DISTANCE_NONE(evt.distance)) break;
+
+ /* Compute the event position */
+ next_pos[0] = pos[0] + dir[0] * evt.distance;
+ next_pos[1] = pos[1] + dir[1] * evt.distance;
+ next_pos[2] = pos[2] + dir[2] * evt.distance;
+
+ /* Surface bounce */
+ if(SURFACE_EVENT(&evt)) {
+ double refl; /* Surface reflectivity */
+ surface_bounce(cmd, args, &evt, next_pos, dir, sc_dir, &refl);
+
+ /* Russian roulette wrt surface reflectivity */
+ if(ssp_rng_canonical(args->rng) >= refl) break;
+
+ /* Save current intersected surface to avoid self-intersection when
+ * sampling next extinction event */
+ hit_from = evt.hit;
+
+ /* Volume scattering */
+ } else if(sample_volume_event_type(cmd, args, &evt) == RNATM_RADCOEF_Ks) {
+ volume_scattering(cmd, args, &evt, next_pos, dir, sc_dir);
+ hit_from = S3D_HIT_NULL; /* Reset surface intersection */
+
+ /* Absorption */
+ } else {
+ break;
+ }
+
+ d3_set(pos, next_pos);
+ d3_set(dir, sc_dir);
+ }
+
+ /* Compute the returned weight */
+ wlen_m = args->wlen * 1.e-9;
+ T = get_temperature(cmd, &evt);
+ L = htrdr_planck_monochromatic(wlen_m, T);
+
+ return L;
+}
diff --git a/src/planeto/htrdr_planeto_draw_map.c b/src/planeto/htrdr_planeto_draw_map.c
@@ -23,6 +23,7 @@
#include "core/htrdr_draw_map.h"
#include "core/htrdr_log.h"
#include "core/htrdr_ran_wlen_cie_xyz.h"
+#include "core/htrdr_ran_wlen_planck.h"
#include <rad-net/rnatm.h>
#include <star/scam.h>
@@ -39,8 +40,122 @@ draw_pixel_xwave
const struct htrdr_draw_pixel_args* args,
void* data)
{
- /* TODO */
- (void)htrdr, (void)args, (void)data;
+ struct planeto_compute_radiance_args rad_args =
+ PLANETO_COMPUTE_RADIANCE_ARGS_NULL;
+
+ struct htrdr_estimate L;
+ struct htrdr_accum radiance;
+ struct htrdr_accum time;
+ struct htrdr_planeto* cmd;
+ struct planeto_pixel_xwave* pixel = data;
+ size_t isamp = 0;
+ ASSERT(htrdr && htrdr_draw_pixel_args_check(args) && data);
+ (void)htrdr;
+
+ cmd = args->context;
+ ASSERT(cmd);
+ ASSERT(cmd->spectral_domain.spectral_type == HTRDR_SPECTRAL_SW
+ || cmd->spectral_domain.spectral_type == HTRDR_SPECTRAL_LW);
+ ASSERT(cmd->output_type == HTRDR_PLANETO_ARGS_OUTPUT_IMAGE);
+
+ /* Reset accumulators */
+ radiance = HTRDR_ACCUM_NULL;
+ time = HTRDR_ACCUM_NULL;
+
+ FOR_EACH(isamp, 0, args->spp) {
+ struct time t0, t1;
+ struct scam_sample sample = SCAM_SAMPLE_NULL;
+ struct scam_ray ray = SCAM_RAY_NULL;
+ double weight;
+ double r0, r1, r2;
+ double wlen[2]; /* Sampled wavelength */
+ double pdf;
+ size_t iband[2];
+ 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 */
+ sample.film[0] = (double)args->pixel_coord[0]+ssp_rng_canonical(args->rng);
+ sample.film[1] = (double)args->pixel_coord[1]+ssp_rng_canonical(args->rng);
+ sample.film[0] *= args->pixel_normalized_size[0];
+ sample.film[1] *= args->pixel_normalized_size[1];
+ sample.lens[0] = ssp_rng_canonical(args->rng);
+ sample.lens[1] = ssp_rng_canonical(args->rng);
+
+ /* Generate a camera ray */
+ scam_generate_ray(cmd->camera, &sample, &ray);
+
+ r0 = ssp_rng_canonical(args->rng);
+ r1 = ssp_rng_canonical(args->rng);
+ r2 = ssp_rng_canonical(args->rng);
+
+ /* Sample a wavelength */
+ wlen[0] = wlen[1] = htrdr_ran_wlen_planck_sample(cmd->planck, r0, r1, &pdf);
+ pdf *= 1.e9; /* Transform the pdf from nm⁻¹ to m⁻¹ */
+
+ /* Find the band the wavelength belongs to and sample a quadrature point */
+ RNATM(find_bands(cmd->atmosphere, wlen, iband));
+ RNATM(band_sample_quad_pt(cmd->atmosphere, r2, iband[0], &iquad));
+ ASSERT(iband[0] == iband[1]);
+
+ /* Compute the radiance in W/m²/sr/m */
+ d3_set(rad_args.path_org, ray.org);
+ d3_set(rad_args.path_dir, ray.dir);
+ rad_args.rng = args->rng;
+ rad_args.ithread = args->ithread;
+ rad_args.wlen = wlen[0];
+ rad_args.iband = iband[0];
+ rad_args.iquad = iquad;
+ switch(cmd->spectral_domain.spectral_type) {
+ case HTRDR_SPECTRAL_LW:
+ weight = planeto_compute_radiance_lw(cmd, &rad_args);
+ break;
+ case HTRDR_SPECTRAL_SW:
+ weight = planeto_compute_radiance_sw(cmd, &rad_args);
+ break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+ ASSERT(weight >= 0);
+
+ weight /= pdf; /* In W/m²/sr */
+
+ /* End of time recording per realisation */
+ time_sub(&t0, time_current(&t1), &t0);
+ usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
+
+ /* Update the radiance accumulator */
+ radiance.sum_weights += weight;
+ radiance.sum_weights_sqr += weight*weight;
+ radiance.nweights += 1;
+
+ /* Update the per realisation time accumulator */
+ time.sum_weights += usec;
+ time.sum_weights_sqr += usec*usec;
+ time.nweights += 1;
+ }
+
+ /* Flush pixel data */
+ htrdr_accum_get_estimation(&radiance, &L);
+ pixel->radiance = L;
+ pixel->time = time;
+
+ if(cmd->spectral_domain.spectral_type == HTRDR_SPECTRAL_LW) {
+ /* Wavelength in meters */
+ const double wmin_m = cmd->spectral_domain.wlen_range[0] * 1.e-9;
+ const double wmax_m = cmd->spectral_domain.wlen_range[1] * 1.e-9;
+
+ /* Brightness temperatures in W/m²/sr */
+ double T, Tmin, Tmax;
+ T = htrdr_radiance_temperature(cmd->htrdr, wmin_m, wmax_m, L.E);
+ Tmin = htrdr_radiance_temperature(cmd->htrdr, wmin_m, wmax_m, L.E - L.SE);
+ Tmax = htrdr_radiance_temperature(cmd->htrdr, wmin_m, wmax_m, L.E + L.SE);
+
+ pixel->radiance_temperature.E = T;
+ pixel->radiance_temperature.SE = Tmax - Tmin;
+ }
}
static void
@@ -133,7 +248,7 @@ draw_pixel_image
rad_args.wlen = wlen[0];
rad_args.iband = iband[0];
rad_args.iquad = iquad;
- weight = planeto_compute_radiance(cmd, &rad_args);
+ weight = planeto_compute_radiance_sw(cmd, &rad_args);
ASSERT(weight >= 0);
weight /= pdf; /* In W/m²/sr */
