htpp

commit fbe414e922dbab69858c489f8c6fdb82451b25ed
parent f217010d2574b20cf993cc4732f4b3d31d5707f4
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue,  2 Oct 2018 16:23:32 +0200

Full rewrite of the color correction process

Diffstat:
M
src/htpp.c
 | 
73
+++++++++++++++++++++++++++++++------------------------------------------

1 file changed, 31 insertions(+), 42 deletions(-)
diff --git a/src/htpp.c b/src/htpp.c
@@ -128,6 +128,16 @@ error:
   goto exit;
 }
 
+static int
+cmp_dbl(const void* a, const void* b)
+{
+  const double* dbl0 = a;
+  const double* dbl1 = b;
+  if(*dbl0 < *dbl1) return-1;
+  if(*dbl0 > *dbl1) return 1;
+  return 0;
+}
+
 static res_T
 open_output_stream(const char* path, const int force_overwrite, FILE** stream)
 {
@@ -203,7 +213,7 @@ error:
 
 /* http://filmicworlds.com/blog/filmic-tonemapping-operators/ */
 static double*
-filmic_tone_mapping(double pixel[3])
+filmic_tone_mapping(double pixel[3], const double exposure, const double Ymax)
 {
   const double A = 0.15;
   const double B = 0.50;
@@ -211,14 +221,14 @@ filmic_tone_mapping(double pixel[3])
   const double D = 0.20;
   const double E = 0.02;
   const double F = 0.30;
-  const double W = 11.2;
+  const double W = Ymax;
   #define TONE_MAP(X) ((((X)*(A*(X)+C*B)+D*E)/((X)*(A*(X)+B)+D*F))-E/F)
   const double white_scale = TONE_MAP(W);
   ASSERT(pixel);
 
-  pixel[0] = TONE_MAP(pixel[0]) * white_scale;
-  pixel[1] = TONE_MAP(pixel[1]) * white_scale;
-  pixel[2] = TONE_MAP(pixel[2]) * white_scale;
+  pixel[0] = TONE_MAP(pixel[0]) / white_scale * exposure;
+  pixel[1] = TONE_MAP(pixel[1]) / white_scale * exposure;
+  pixel[2] = TONE_MAP(pixel[2]) / white_scale * exposure;
   #undef TONE_MAP
   return pixel;
 }
@@ -402,49 +412,33 @@ error:
 static double
 compute_img_normalization_factor(const struct img* img)
 {
-  size_t histo[16] = {0};
-  const size_t histo_len = sizeof(histo) / sizeof(histo[0]);
-  double Yscale;
+  double* Y = NULL;
   double Ymax;
-  double Ydelta;
-  double rcp_npixels;
   size_t i, x, y;
   ASSERT(img);
 
-  /* Use nextarfter on the Luminance range to prevent a normalized luminance
-   * equal to 1, i.e. to prevent an histogram index equal to histo_len. */
-  Yscale = 1.0 / nextafter(img->Yrange[1] - img->Yrange[0], DBL_MAX);
+  CHK(Y = mem_alloc(sizeof(*Y)*img->width*img->height));
 
-  /* Build the histo */
+  /* Copy the pixel luminance in the Y array */
+  i = 0;
   FOR_EACH(y, 0, img->height) {
     const double* row = (double*)(img->pixels + y*img->pitch);
     FOR_EACH(x, 0, img->width) {
       const double* pixel = row + x*3;
-      const double Y = pixel[1];
-      const double Ynormalized = (Y - img->Yrange[0]) * Yscale;
-      const size_t ientry = (size_t)((double)histo_len * Ynormalized);
-      ASSERT(ientry < histo_len);
-      ++histo[ientry];
+      Y[i] = pixel[1];
+      i++;
     }
   }
 
-  rcp_npixels = 1.0 / (double)(img->width * img->height);
-
-  /* Find the first histogram entry whose cumulative contains N% of the pixels.
-   * N is defined empirically to 90% */
-  FOR_EACH(i, 0, histo_len) {
-    if(i > 0) histo[i] += histo[i-1]; /* Compute the cumulative */
-    if((double)histo[i] * rcp_npixels > 0.9) break;
-  }
-
-  /* Define the width of an histogram entry */
-  Ydelta = (img->Yrange[1] - img->Yrange[0]) / (double)histo_len;
+  /* Sort the Pixel luminance in ascending order */
+  qsort(Y, img->width*img->height, sizeof(*Y), cmp_dbl);
 
-  /* Compute the upper bound luminance value of the found histogram entry */
-  Ymax = img->Yrange[0]
-       + (double)(i+1/*Upper bound of the histogram entry*/)*Ydelta;
+  /* Arbitrarly use the luminance at 99.5% of the overal pixel at the maxium
+   * luminance of the image */
+  i = (size_t)((double)(img->width*img->height)*0.995);
+  Ymax = Y[i];
 
-  /* Use the previous value as the normalization factor */
+  mem_rm(Y);
   return Ymax;
 }
 
@@ -460,7 +454,7 @@ main(int argc, char** argv)
   const char* stream_in_name = "stdin";
   struct img img = IMG_NULL;
   struct args args = ARGS_DEFAULT;
-  double Yscale;
+  double Ymax;
   const int nprocs = omp_get_num_procs();
   int err = 0;
   int64_t i;
@@ -484,7 +478,7 @@ main(int argc, char** argv)
   res = img_load(&img, args.uncertainties, stream_in, stream_in_name);
   if(res != RES_OK) goto error;
 
-  Yscale = args.exposure / compute_img_normalization_factor(&img);
+  Ymax = compute_img_normalization_factor(&img);
 
   omp_set_num_threads(nprocs);
   /* Convert input HDR XYZ image in RGB LDR image */
@@ -495,14 +489,9 @@ main(int argc, char** argv)
     double* row = (double*)(img.pixels + img.pitch*y);
     double* pixel = row + x*3;
 
-    /* Scale the pixel value */
-    pixel[0] = pixel[0] * Yscale;
-    pixel[1] = pixel[1] * Yscale;
-    pixel[2] = pixel[2] * Yscale;
-
+    filmic_tone_mapping(pixel, args.exposure, Ymax); /* Tone map the RGB pixel */
     if(!args.uncertainties) {
       XYZ_to_RGB(pixel); /* Convert in RGB color space */
-      filmic_tone_mapping(pixel); /* Tone map the RGB pixel */
       RGB_to_sRGB(pixel); /* Convert in sRGB color space (i.e. gamma correction) */
     }
   }