htrdr

commit 93b11041d1cdebe18f985e97b39f98bd26681f25
parent 0ab02d8283977bdd9a29e0dd8f392d0201257fde
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri, 20 Mar 2020 11:17:43 +0100

Fix the Long Wave integration weight

Divide the weight by the pdf of the sampled spectral band.

Diffstat:
M
src/htrdr.c
 | 
20
+++++++++++++++-----
M
src/htrdr.h
 | 
3
++-
M
src/htrdr_draw_radiance.c
 | 
15
++++++++++++---

3 files changed, 29 insertions(+), 9 deletions(-)
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -101,7 +101,7 @@ log_msg
 }
 
 static res_T
-dump_accum_buffer
+dump_buffer
   (struct htrdr* htrdr,
    struct htrdr_buffer* buf,
    struct htrdr_accum* time_acc, /* May be NULL */
@@ -389,6 +389,7 @@ setup_lw_cdf(struct htrdr* htrdr)
   ASSERT(htrdr && htsky_is_long_wave(htrdr->sky));
 
   nbands = htsky_get_spectral_bands_count(htrdr->sky);
+
   res = darray_double_resize(&htrdr->lw_cdf, nbands);
   if(res != RES_OK) {
     htrdr_log_err(htrdr,
@@ -397,9 +398,15 @@ setup_lw_cdf(struct htrdr* htrdr)
     goto error;
   }
 
-  /* Alias the same array by the cdf and the pdf variable to make easier the
-   * reading of the code */
-  pdf = darray_double_data_get(&htrdr->lw_cdf);
+  res = darray_double_resize(&htrdr->lw_pdf, nbands);
+  if(res != RES_OK) {
+    htrdr_log_err(htrdr,
+      "error allocating the PDF of the long wave spectral bands -- %s.\n",
+      res_to_cstr(res));
+    goto error;
+  }
+
+  pdf = darray_double_data_get(&htrdr->lw_pdf);
   cdf = darray_double_data_get(&htrdr->lw_cdf);
 
   /* Compute the *unormalized* probability to sample a long wave band */
@@ -440,6 +447,7 @@ exit:
   return res;
 error:
   darray_double_clear(&htrdr->lw_cdf);
+  darray_double_clear(&htrdr->lw_pdf);
   goto exit;
 }
 
@@ -475,6 +483,7 @@ htrdr_init
   str_init(htrdr->allocator, &htrdr->output_name);
 
   darray_double_init(htrdr->allocator, &htrdr->lw_cdf);
+  darray_double_init(htrdr->allocator, &htrdr->lw_pdf);
 
   nthreads_max = MMAX(omp_get_max_threads(), omp_get_num_procs());
   htrdr->dump_vtk = args->dump_vtk;
@@ -645,6 +654,7 @@ htrdr_release(struct htrdr* htrdr)
   }
   str_release(&htrdr->output_name);
   darray_double_release(&htrdr->lw_cdf);
+  darray_double_release(&htrdr->lw_pdf);
   logger_release(&htrdr->logger);
 }
 
@@ -678,7 +688,7 @@ htrdr_run(struct htrdr* htrdr)
       struct htrdr_accum path_time_acc = HTRDR_ACCUM_NULL;
       struct htrdr_estimate path_time;
 
-      res = dump_accum_buffer(htrdr, htrdr->buf, &path_time_acc,
+      res = dump_buffer(htrdr, htrdr->buf, &path_time_acc,
         str_cget(&htrdr->output_name), htrdr->output);
       if(res != RES_OK) goto error;
 
diff --git a/src/htrdr.h b/src/htrdr.h
@@ -56,7 +56,8 @@ struct htrdr {
   struct htsky* sky;
   double wlen_range_m[2]; /* Integration range in *meters* */
 
-  struct darray_double lw_cdf; /* CDF to sample a Long Waves band */
+  struct darray_double lw_cdf; /* CDF to sample a long wave bands */
+  struct darray_double lw_pdf; /* PDF of the long wave  bands*/
 
   size_t spp; /* #samples per pixel */
   size_t width; /* Image width */
diff --git a/src/htrdr_draw_radiance.c b/src/htrdr_draw_radiance.c
@@ -479,14 +479,19 @@ error:
 }
 
 static INLINE size_t
-sample_lw_spectral_interval(struct htrdr* htrdr, const double r)
+sample_lw_spectral_interval
+  (struct htrdr* htrdr,
+   const double r,
+   double* pdf)
 {
   const double* cdf = NULL;
   const double* find = NULL;
   double r_next = nextafter(r, DBL_MAX);
   size_t cdf_length = 0;
   size_t i;
-  ASSERT(htrdr && r >= 0 && r < 1);
+  ASSERT(htrdr && r >= 0 && r < 1 && pdf);
+  ASSERT(darray_double_size_get(&htrdr->lw_cdf)
+      == darray_double_size_get(&htrdr->lw_pdf));
 
   cdf = darray_double_cdata_get(&htrdr->lw_cdf);
   cdf_length = darray_double_size_get(&htrdr->lw_cdf);
@@ -498,6 +503,8 @@ sample_lw_spectral_interval(struct htrdr* htrdr, const double r)
 
   i = (size_t)(find - cdf);
   ASSERT(i < cdf_length && cdf[i] > r && (!i || cdf[i-1] <= r));
+  *pdf = darray_double_cdata_get(&htrdr->lw_pdf)[i];
+
   return htsky_get_spectral_band_id(htrdr->sky, i);
 }
 
@@ -629,6 +636,7 @@ draw_pixel_lw
     size_t iband;
     size_t iquad;
     double usec;
+    double band_pdf;
 
     /* Begin the registration of the time spent to in the realisation */
     time_current(&t0);
@@ -645,12 +653,13 @@ draw_pixel_lw
     r1 = ssp_rng_canonical(rng);
 
     /* Sample a spectral band and a quadrature point */
-    iband = sample_lw_spectral_interval(htrdr, r0);
+    iband = sample_lw_spectral_interval(htrdr, r0, &band_pdf);
     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);
+    weight /= band_pdf;
     ASSERT(weight >= 0);
 
     /* End the registration of the per realisation time */