htpp

htrdr-image post-processing
git clone git://git.meso-star.fr/htpp.git
Log | Files | Refs | README | LICENSE

htpp.c (27860B)


      1 /* Copyright (C) 2018-2020, 2023 |Méso|Star> (contact@meso-star.com)
      2  * Copyright (C) 2018, 2019 Centre National de la Recherche Scientifique
      3  * Copyright (C) 2018, 2019 Université Paul Sabatier
      4  *
      5  * This program is free software: you can redistribute it and/or modify
      6  * it under the terms of the GNU General Public License as published by
      7  * the Free Software Foundation, either version 3 of the License, or
      8  * (at your option) any later version.
      9  *
     10  * This program is distributed in the hope that it will be useful,
     11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     13  * GNU General Public License for more details.
     14  *
     15  * You should have received a copy of the GNU General Public License
     16  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     17 
     18 #define _POSIX_C_SOURCE 200112L /* getopt/close support */
     19 
     20 #include "htpp_version.h"
     21 
     22 #include <rsys/cstr.h>
     23 #include <rsys/double33.h>
     24 #include <rsys/mem_allocator.h>
     25 #include <rsys/rsys.h>
     26 #include <rsys/text_reader.h>
     27 
     28 #include <star/scmap.h>
     29 
     30 #include <errno.h>
     31 #include <fcntl.h> /* open */
     32 #include <omp.h>
     33 #include <string.h>
     34 #include <sys/stat.h> /* S_IRUSR & S_IWUSR */
     35 #include <unistd.h> /* getopt & close functions */
     36 
     37 enum pixcpnt {
     38   PIXCPNT_X,
     39   PIXCPNT_R = PIXCPNT_X,
     40   PIXCPNT_X_STDERR,
     41   PIXCPNT_Y,
     42   PIXCPNT_G = PIXCPNT_Y,
     43   PIXCPNT_Y_STDERR,
     44   PIXCPNT_Z,
     45   PIXCPNT_B = PIXCPNT_Z,
     46   PIXCPNT_Z_STDERR,
     47   PIXCPNT_TIME,
     48   PIXCPNT_TIME_STDERR,
     49   PIXCPNTS_COUNT__
     50 };
     51 
     52 enum pp_type {
     53   PP_IMAGE,
     54   PP_MAP
     55 };
     56 
     57 struct args {
     58   const char* input;
     59   const char* output;
     60 
     61   enum pp_type pp_type;
     62 
     63   struct image_opt {
     64     double exposure; /* In [0, inf) */
     65     double white; /* In ]0, inf) */
     66   } image;
     67 
     68   struct map_opt {
     69     const struct scmap_palette* palette;
     70     unsigned pixcpnt; /* In [0, PIXCPNTS_COUNT__[ */
     71     double range[2];
     72     int gnuplot;
     73   } map;
     74 
     75   int verbose;
     76   int force_overwrite;
     77   int nthreads;
     78   int quit;
     79 };
     80 #define ARGS_DEFAULT__ {                                                       \
     81   NULL,                     /* Input */                                        \
     82   NULL,                     /* Output */                                       \
     83   PP_IMAGE,                 /* Post process type */                            \
     84   {                                                                            \
     85     1.0,                    /* Image exposure */                               \
     86     -1.0,                   /* Image white scale */                            \
     87   }, {                                                                         \
     88     &scmap_palette_inferno, /* Map palette */                                  \
     89     0,                      /* Map channel */                                  \
     90     {DBL_MAX,-DBL_MAX},     /* Range */                                        \
     91     0                       /* Gnuplot */                                      \
     92   },                                                                           \
     93   0,                        /* Verbosity level */                              \
     94   0,                        /* Force overwrite? */                             \
     95   INT_MAX,                  /* #threads */                                     \
     96   0                         /* Quit? */                                        \
     97 }
     98 
     99 struct img {
    100   char* pixels; /* row majored pixels */
    101   size_t width;
    102   size_t height;
    103   size_t pitch; /* #bytes of a row */
    104 
    105   /* Ranges of the loaded value */
    106   double ranges[PIXCPNTS_COUNT__][2];
    107 };
    108 #define IMG_NULL__ { NULL, 0, 0, 0, {{ 0, 0}} }
    109 static const struct img IMG_NULL = IMG_NULL__;
    110 
    111 /*******************************************************************************
    112  * Helper functions
    113  ******************************************************************************/
    114 static void
    115 usage(void)
    116 {
    117   printf(
    118 "usage: htpp [-fhVv] [-i image_option[:image_option ...]]\n"
    119 "            [-m map_option[:map_option ...]] [-o output]\n"
    120 "            [-t threads_count] [input]\n");
    121 }
    122 
    123 static res_T
    124 parse_multiple_options
    125   (struct args* args,
    126    const char* str,
    127    res_T (*parse_option)(struct args* args, const char* str))
    128 {
    129   char buf[512];
    130   char* tk;
    131   char* ctx;
    132   res_T res = RES_OK;
    133   ASSERT(args && str);
    134 
    135   if(strlen(str) >= sizeof(buf) - 1/*NULL char*/) {
    136     fprintf(stderr, "Could not duplicate the option string `%s'.\n", str);
    137     res = RES_MEM_ERR;
    138     goto error;
    139   }
    140   strncpy(buf, str, sizeof(buf));
    141 
    142   tk = strtok_r(buf, ":", &ctx);
    143   do {
    144     res = parse_option(args, tk);
    145     if(res != RES_OK) goto error;
    146     tk = strtok_r(NULL, ":", &ctx);
    147   } while(tk);
    148 
    149 exit:
    150   return res;
    151 error:
    152   goto exit;
    153 }
    154 
    155 static res_T
    156 parse_img_option(struct args* args, const char* str)
    157 {
    158   char buf[128];
    159   char* key;
    160   char* val;
    161   char* tk_ctx;
    162   const struct args args_default = ARGS_DEFAULT__;
    163   res_T res = RES_OK;
    164 
    165   if(strlen(str) >= sizeof(buf) - 1/*NULL char*/) {
    166     fprintf(stderr,
    167       "Could not duplicate the image options string `%s'.\n", str);
    168     res = RES_MEM_ERR;
    169     goto error;
    170   }
    171   strncpy(buf, str, sizeof(buf));
    172 
    173   key = strtok_r(buf, "=", &tk_ctx);
    174   val = strtok_r(NULL, "", &tk_ctx);
    175 
    176   if(!strcmp(key, "default")) {
    177     if(val) {
    178       fprintf(stderr, "Unexpected value to the image option `%s'.\n", key);
    179       res = RES_BAD_ARG;
    180       goto error;
    181     }
    182     args->image = args_default.image;
    183 
    184   } else {
    185 
    186     if(!val) {
    187       fprintf(stderr, "Missing value to the image option `%s'.\n", key);
    188       res = RES_BAD_ARG;
    189       goto error;
    190     }
    191 
    192     if(!strcmp(key, "exposure")) {
    193       res = cstr_to_double(val, &args->image.exposure);
    194       if(res != RES_OK) goto error;
    195       if(args->image.exposure < 0) {
    196         fprintf(stderr, "Invalid image exposure %g.\n", args->image.exposure);
    197         res = RES_BAD_ARG;
    198         goto error;
    199       }
    200     } else if(!strcmp(key, "white")) {
    201       res = cstr_to_double(val, &args->image.white);
    202       if(res != RES_OK) goto error;
    203       if(args->image.exposure < 0) {
    204         fprintf(stderr, "Invalid image white scale %g.\n", args->image.white);
    205         res = RES_BAD_ARG;
    206         goto error;
    207       }
    208     } else {
    209       fprintf(stderr, "Invalid image option `%s'.\n", key);
    210       res = RES_BAD_ARG;
    211       goto error;
    212     }
    213   }
    214 
    215 exit:
    216   return res;
    217 error:
    218   goto exit;
    219 }
    220 
    221 static res_T
    222 parse_map_option(struct args* args, const char* str)
    223 {
    224   char buf[128];
    225   char* key;
    226   char* val;
    227   char* tk_ctx;
    228   const struct args args_default = ARGS_DEFAULT__;
    229   res_T res = RES_OK;
    230 
    231   if(strlen(str) >= sizeof(buf) - 1/*NULL char*/) {
    232     fprintf(stderr,
    233       "Could not duplicate the map options string `%s'.\n", str);
    234     res = RES_MEM_ERR;
    235     goto error;
    236   }
    237   strncpy(buf, str, sizeof(buf));
    238 
    239   key = strtok_r(buf, "=", &tk_ctx);
    240   val = strtok_r(NULL, "", &tk_ctx);
    241 
    242   if(!strcmp(key, "default")) {
    243     if(val) {
    244       fprintf(stderr, "Unexpected value to the map option `%s'.\n", key);
    245       res = RES_BAD_ARG;
    246       goto error;
    247     }
    248     args->map = args_default.map;
    249   } else if(!strcmp(key, "gnuplot")) {
    250     args->map.gnuplot = 1;
    251   } else {
    252     if(!val) {
    253       fprintf(stderr, "Missing value to the map option `%s'.\n", key);
    254       res = RES_BAD_ARG;
    255       goto error;
    256     }
    257 
    258     if(!strcmp(key, "pixcpnt")) {
    259       res = cstr_to_uint(val, &args->map.pixcpnt);
    260       if(res != RES_OK) goto error;
    261       if(args->map.pixcpnt >= PIXCPNTS_COUNT__) {
    262         fprintf(stderr, "Invalid pixel component `%u'.\n", args->map.pixcpnt);
    263         res = RES_BAD_ARG;
    264         goto error;
    265       }
    266     } else if(!strcmp(key, "palette")) {
    267       args->map.palette = scmap_get_builtin_palette(val);
    268       if(!args->map.palette) {
    269         fprintf(stderr, "Invalid palette `%s'.\n", val);
    270         res = RES_BAD_ARG;
    271         goto error;
    272       }
    273     } else if(!strcmp(key, "range")) {
    274       size_t len;
    275       res = cstr_to_list_double(val, ',', args->map.range, &len, 2);
    276       if(res != RES_OK) goto error;
    277     } else {
    278       fprintf(stderr, "Invalid map option `%s'.\n", key);
    279       res = RES_BAD_ARG;
    280       goto error;
    281     }
    282   }
    283 
    284 exit:
    285   return res;
    286 error:
    287   goto exit;
    288 }
    289 
    290 static void
    291 args_release(struct args* args)
    292 {
    293   const struct args args_default = ARGS_DEFAULT__;
    294   ASSERT(args);
    295   *args = args_default;
    296 }
    297 
    298 static res_T
    299 args_init(struct args* args, const int argc, char** argv)
    300 {
    301   int opt;
    302   res_T res = RES_OK;
    303   ASSERT(args && argc && argv);
    304 
    305   /* Begin the optstring by ':' to make silent getopt */
    306   while((opt = getopt(argc, argv, "fhi:m:o:t:vV")) != -1) {
    307     switch(opt) {
    308       case 'f': args->force_overwrite = 1; break;
    309       case 'h':
    310         usage();
    311         args_release(args);
    312         args->quit = 1;
    313         goto exit;
    314       case 'i':
    315         args->pp_type = PP_IMAGE;
    316         res = parse_multiple_options(args, optarg, parse_img_option);
    317         break;
    318       case 'm':
    319         args->pp_type = PP_MAP;
    320         res = parse_multiple_options(args, optarg, parse_map_option);
    321         break;
    322       case 'o': args->output = optarg; break;
    323       case 't':
    324         res = cstr_to_int(optarg, &args->nthreads);
    325         if(res == RES_OK && args->nthreads <= 0) res = RES_BAD_ARG;
    326         break;
    327       case 'v': args->verbose = 1; break;
    328       case 'V':
    329         printf("High-Tune: Post-Process %d.%d.%d\n",
    330           HTPP_VERSION_MAJOR,
    331           HTPP_VERSION_MINOR,
    332           HTPP_VERSION_PATCH);
    333         args->quit = 1;
    334         goto exit;
    335       default: res = RES_BAD_ARG; break;
    336     }
    337     if(res != RES_OK) {
    338       if(optarg) {
    339         fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n",
    340           argv[0], optarg, opt);
    341       }
    342       goto error;
    343     }
    344   }
    345 
    346   if(optind < argc) {
    347     args->input = argv[optind];
    348   }
    349 
    350   args->nthreads = MMIN(omp_get_num_procs(), args->nthreads);
    351 
    352 exit:
    353   optind = 1;
    354   return res;
    355 error:
    356   usage();
    357   args_release(args);
    358   goto exit;
    359 }
    360 
    361 static int
    362 cmp_dbl(const void* a, const void* b)
    363 {
    364   const double* dbl0 = a;
    365   const double* dbl1 = b;
    366   if(*dbl0 < *dbl1) return-1;
    367   if(*dbl0 > *dbl1) return 1;
    368   return 0;
    369 }
    370 
    371 static res_T
    372 open_output_stream(const char* path, const int force_overwrite, FILE** stream)
    373 {
    374   int fd = -1;
    375   FILE* fp = NULL;
    376   res_T res = RES_OK;
    377   ASSERT(path && stream);
    378 
    379   if(force_overwrite) {
    380     fp = fopen(path, "w");
    381     if(!fp) {
    382       fprintf(stderr, "Could not open the output file '%s' -- %s.\n",
    383         path, strerror(errno));
    384       goto error;
    385     }
    386   } else {
    387     fd = open(path, O_CREAT|O_WRONLY|O_EXCL|O_TRUNC, S_IRUSR|S_IWUSR);
    388     if(fd >= 0) {
    389       fp = fdopen(fd, "w");
    390       if(fp == NULL) {
    391         fprintf(stderr, "Could not open the output file '%s' -- %s.\n",
    392           path, strerror(errno));
    393         goto error;
    394       }
    395     } else if(errno == EEXIST) {
    396       fprintf(stderr,
    397         "The output file '%s' already exists. Use -f to overwrite it.\n", path);
    398       goto error;
    399     } else {
    400       fprintf(stderr,
    401         "Unexpected error while opening the output file `'%s' -- %s.\n",
    402          path, strerror(errno));
    403       goto error;
    404     }
    405   }
    406 
    407 exit:
    408   *stream = fp;
    409   return res;
    410 error:
    411   res = RES_IO_ERR;
    412   if(fp) {
    413     CHK(fclose(fp) == 0);
    414     fp = NULL;
    415   } else if(fd >= 0) {
    416     CHK(close(fd) == 0);
    417   }
    418   goto exit;
    419 }
    420 
    421 static res_T
    422 open_input_stream(const char* path, FILE** stream)
    423 {
    424   FILE* fp = NULL;
    425   res_T res = RES_OK;
    426   ASSERT(path && stream);
    427 
    428   fp = fopen(path, "r");
    429   if(!fp) {
    430     fprintf(stderr, "Could not open the input file '%s' -- %s.\n",
    431       path, strerror(errno));
    432     res = RES_IO_ERR;
    433     goto error;
    434   }
    435 
    436 exit:
    437   *stream = fp;
    438   return res;
    439 error:
    440   if(fp) fclose(fp);
    441   goto exit;
    442 }
    443 
    444 /* http://filmicworlds.com/blog/filmic-tonemapping-operators/ */
    445 static double*
    446 filmic_tone_mapping
    447   (double pixel[PIXCPNTS_COUNT__],
    448    const double exposure,
    449    const double Ymax)
    450 {
    451   const double A = 0.15;
    452   const double B = 0.50;
    453   const double C = 0.10;
    454   const double D = 0.20;
    455   const double E = 0.02;
    456   const double F = 0.30;
    457   const double W = Ymax;
    458   #define TONE_MAP(X) ((((X)*(A*(X)+C*B)+D*E)/((X)*(A*(X)+B)+D*F))-E/F)
    459   const double white_scale = TONE_MAP(W*exposure);
    460   ASSERT(pixel);
    461 
    462   pixel[PIXCPNT_X] = TONE_MAP(pixel[PIXCPNT_X]*exposure) / white_scale;
    463   pixel[PIXCPNT_Y] = TONE_MAP(pixel[PIXCPNT_Y]*exposure) / white_scale;
    464   pixel[PIXCPNT_Z] = TONE_MAP(pixel[PIXCPNT_Z]*exposure) / white_scale;
    465   #undef TONE_MAP
    466   return pixel;
    467 }
    468 
    469 static double*
    470 XYZ_to_sRGB(double pixel[PIXCPNTS_COUNT__])
    471 {
    472   #define D65_x 0.31271
    473   #define D65_y 0.32902
    474   const double D65[3] = {D65_x / D65_y, 1.0, (1 - D65_x - D65_y) / D65_y};
    475   const double mat[9] = { /* XYZ to sRGB matrix */
    476      3.2404542, -0.9692660,  0.0556434,
    477     -1.5371385,  1.8760108, -0.2040259,
    478     -0.4985314,  0.0415560,  1.0572252
    479   };
    480   double XYZ[3], sRGB[3], XYZ_D65[3];
    481   ASSERT(pixel);
    482 
    483   XYZ[0] = pixel[PIXCPNT_X];
    484   XYZ[1] = pixel[PIXCPNT_Y];
    485   XYZ[2] = pixel[PIXCPNT_Z];
    486 
    487   d3_mul(XYZ_D65, XYZ, D65);
    488   d33_muld3(sRGB, mat, XYZ_D65);
    489   sRGB[0] = MMAX(sRGB[0], 0);
    490   sRGB[1] = MMAX(sRGB[1], 0);
    491   sRGB[2] = MMAX(sRGB[2], 0);
    492 
    493   pixel[PIXCPNT_R] = sRGB[0];
    494   pixel[PIXCPNT_G] = sRGB[1];
    495   pixel[PIXCPNT_B] = sRGB[2];
    496   return pixel;
    497 }
    498 
    499 static double
    500 sRGB_gamma_correct(double value)
    501 {
    502   if(value <= 0.0031308) {
    503     return value * 12.92;
    504   } else {
    505     return 1.055 * pow(value, 1.0/2.4) - 0.055;
    506   }
    507 }
    508 
    509 static void
    510 img_release(struct img* img)
    511 {
    512   ASSERT(img);
    513   if(img->pixels) mem_rm(img->pixels);
    514 }
    515 
    516 static res_T
    517 img_load
    518   (struct img* img,
    519    FILE* stream,
    520    const char* stream_name)
    521 {
    522   struct txtrdr* txtrdr = NULL;
    523   size_t icpnt;
    524   size_t x, y;
    525   unsigned resolution[2];
    526   res_T res = RES_OK;
    527   ASSERT(img && stream);
    528 
    529   *img = IMG_NULL;
    530 
    531   res = txtrdr_stream(NULL, stream, stream_name, '#', &txtrdr);
    532   if(res != RES_OK) {
    533     fprintf(stderr, "%s: could not setup the text reader -- %s.\n",
    534       txtrdr_get_name(txtrdr), res_to_cstr(res));
    535     goto error;
    536   }
    537 
    538   res = txtrdr_read_line(txtrdr);
    539   if(res != RES_OK) {
    540     fprintf(stderr, "%s: could not read the image resolution -- %s.\n",
    541       txtrdr_get_name(txtrdr), res_to_cstr(res));
    542     goto error;
    543   }
    544 
    545   if(!txtrdr_get_line(txtrdr)) {
    546     fprintf(stderr, "%s: missing image definition.\n", txtrdr_get_name(txtrdr));
    547     res = RES_BAD_ARG;
    548     goto error;
    549   }
    550 
    551   res = cstr_to_list_uint(txtrdr_get_cline(txtrdr), ' ', resolution, NULL, 2);
    552   if(res != RES_OK) {
    553     fprintf(stderr, "%s: invalid image resolution `%s' -- %s\n",
    554       txtrdr_get_name(txtrdr), txtrdr_get_cline(txtrdr), res_to_cstr(res)) ;
    555     goto error;
    556   }
    557 
    558   img->width = (size_t)resolution[0];
    559   img->height = (size_t)resolution[1];
    560   img->pitch = ALIGN_SIZE(sizeof(double[PIXCPNTS_COUNT__])*img->width, 16u);
    561   img->pixels = mem_calloc(img->height, img->pitch);
    562   if(!img->pixels) {
    563     fprintf(stderr, "Could not allocate the image pixels.\n");
    564     res = RES_MEM_ERR;
    565     goto error;
    566   }
    567 
    568   /* Reset the range of each pixel components */
    569   FOR_EACH(icpnt, 0, PIXCPNTS_COUNT__) {
    570     img->ranges[icpnt][0] = DBL_MAX;
    571     img->ranges[icpnt][1] = -DBL_MAX;
    572   }
    573 
    574   /* Read pixel components */
    575   FOR_EACH(y, 0, img->height) {
    576     double* row = (double*)(img->pixels + y*img->pitch);
    577     FOR_EACH(x, 0, img->width) {
    578       double* pixel = row + x*PIXCPNTS_COUNT__;
    579       size_t ncpnts;
    580 
    581       res = txtrdr_read_line(txtrdr);
    582       if(res != RES_OK) {
    583         fprintf(stderr,
    584           "%s: could not read the components of the pixel (%lu, %lu) -- %s.\n",
    585           txtrdr_get_name(txtrdr), (unsigned long)x, (unsigned long)y,
    586           res_to_cstr(res));
    587         goto error;
    588       }
    589 
    590       res = cstr_to_list_double(txtrdr_get_cline(txtrdr), ' ', pixel, &ncpnts, 8);
    591       if(res != RES_OK || ncpnts < 6) {
    592         fprintf(stderr, "%s: invalid components for the (%lu, %lu) pixel.\n",
    593           txtrdr_get_name(txtrdr), (unsigned long)x, (unsigned long)y);
    594         res = RES_BAD_ARG;
    595         goto error;
    596       }
    597 
    598       FOR_EACH(icpnt, 0, ncpnts) {
    599         img->ranges[icpnt][0] = MMIN(img->ranges[icpnt][0], pixel[icpnt]);
    600         img->ranges[icpnt][1] = MMAX(img->ranges[icpnt][1], pixel[icpnt]);
    601       }
    602     }
    603   }
    604 exit:
    605   if(txtrdr) txtrdr_ref_put(txtrdr);
    606   return res;
    607 error:
    608   img_release(img);
    609   goto exit;
    610 }
    611 
    612 static res_T
    613 img_write_ppm(const struct img* img, FILE* stream, const char* stream_name)
    614 {
    615   size_t x, y;
    616   int i;
    617   res_T res = RES_OK;
    618   ASSERT(img && stream && stream_name);
    619 
    620   /* Write LDR image */
    621   i = fprintf(stream, "P3 %lu %lu\n255\n",
    622     (unsigned long)img->width, (unsigned long)img->height);
    623   if(i < 0) {
    624     fprintf(stderr, "%s: could not write the PPM header.\n", stream_name);
    625     res = RES_IO_ERR;
    626     goto error;
    627   }
    628   FOR_EACH(y, 0, img->height) {
    629     const double* row = (double*)(img->pixels + y*img->pitch);
    630     FOR_EACH(x, 0, img->width) {
    631       const double* pixel = row + x*PIXCPNTS_COUNT__;
    632       i = fprintf(stream, "%i %i %i\n",
    633         (uint8_t)(CLAMP(pixel[PIXCPNT_X], 0.0, 1.0) * 255.0 + 0.5/*Round*/),
    634         (uint8_t)(CLAMP(pixel[PIXCPNT_Y], 0.0, 1.0) * 255.0 + 0.5/*Round*/),
    635         (uint8_t)(CLAMP(pixel[PIXCPNT_Z], 0.0, 1.0) * 255.0 + 0.5/*Round*/));
    636       if(i < 0) {
    637         fprintf(stderr, "%s: could not write the (%lu, %lu) pixel.\n",
    638           stream_name, (unsigned long)x, (unsigned long)y);
    639         res = RES_IO_ERR;
    640         goto error;
    641       }
    642     }
    643   }
    644 
    645 exit:
    646   return res;
    647 error:
    648   goto exit;
    649 }
    650 
    651 static res_T
    652 img_write_gnuplot
    653   (const struct img* img,
    654    const struct args* args,
    655    FILE* stream,
    656    const char* stream_name)
    657 {
    658   double cbox_width = 0.8;
    659   double cbox_height = 0.08;
    660   double cbox_tmargin = 0.02;
    661   size_t icol;
    662   size_t x, y;
    663   res_T res = RES_OK;
    664   ASSERT(img && args && stream && stream_name);
    665 
    666   #define CHKWR(FPrintf) {                                                     \
    667     const int i = FPrintf;                                                     \
    668     if(i < 0) {                                                                \
    669       fprintf(stderr, "%s: could not write the gnuplot map.\n", stream_name);  \
    670       res = RES_IO_ERR;                                                        \
    671       goto error;                                                              \
    672     }                                                                          \
    673   } (void)0
    674   CHKWR(fprintf(stream, "unset xtics\n"));
    675   CHKWR(fprintf(stream, "unset ytics\n"));
    676   CHKWR(fprintf(stream, "unset key\n"));
    677   CHKWR(fprintf(stream, "unset colorbox\n"));
    678   CHKWR(fprintf(stream, "unset origin\n"));
    679   CHKWR(fprintf(stream, "unset border\n"));
    680   CHKWR(fprintf(stream, "unset title\n"));
    681   CHKWR(fprintf(stream, "unset margin\n"));
    682   CHKWR(fprintf(stream, "set margins 0,0,0,0\n"));
    683   CHKWR(fprintf(stream, "set xrange[0:%lu]\n", (unsigned long)img->width-1));
    684   CHKWR(fprintf(stream, "set yrange[%lu:0]\n", (unsigned long)img->height-1));
    685   if(args->map.range[0] < args->map.range[1]) {
    686     CHKWR(fprintf(stream, "set cbrange[%g:%g]\n",
    687       args->map.range[0], args->map.range[1]));
    688   }
    689   CHKWR(fprintf(stream, "set terminal png size %lu,%lu*(1+%g+%g)\n",
    690     (long unsigned)img->width, (unsigned long)img->height,
    691     cbox_height, cbox_tmargin));
    692   CHKWR(fprintf(stream, "set origin 0, %g\n", (cbox_height+cbox_tmargin)*0.5));
    693   CHKWR(fprintf(stream, "set size ratio %g\n",
    694     (double)img->height/(double)img->width));
    695   CHKWR(fprintf(stream, "set colorbox horiz user origin %g,%g size %g,%g\n",
    696     (1.0-cbox_width)*0.5, cbox_height*0.5, cbox_width, cbox_height*0.5));
    697 
    698   CHKWR(fprintf(stream, "set palette defined (\\\n"));
    699   FOR_EACH(icol, 0, args->map.palette->ncolors) {
    700     double col[3];
    701     args->map.palette->get_color(icol, col, args->map.palette->context);
    702     CHKWR(fprintf(stream, "  %lu %g %g %g",
    703       (unsigned long)icol, col[0], col[1], col[2]));
    704     if(icol < args->map.palette->ncolors-1) {
    705       CHKWR(fprintf(stream, ",\\\n"));
    706     } else {
    707       CHKWR(fprintf(stream, ")\n"));
    708     }
    709   }
    710 
    711   CHKWR(fprintf(stream, "$map2 << EOD\n"));
    712   FOR_EACH(y, 0, img->height) {
    713     FOR_EACH(x, 0, img->width) {
    714       double* row = (double*)(img->pixels + img->pitch*y);
    715       double* pixel = row + x*PIXCPNTS_COUNT__;
    716       CHKWR(fprintf(stream, "%lu %lu %g\n",
    717         (unsigned long)y, (unsigned long)x, pixel[args->map.pixcpnt]));
    718     }
    719     if(y != img->height-1) {
    720       CHKWR(fprintf(stream, "\n"));
    721     }
    722   }
    723   CHKWR(fprintf(stream, "EOD\n"));
    724   CHKWR(fprintf(stream, "plot '$map2' using 2:1:3 with image pixels\n"));
    725   #undef CHKWR
    726 
    727 exit:
    728   return res;
    729 error:
    730   goto exit;
    731 }
    732 
    733 static double
    734 compute_XYZ_normalization_factor(const struct img* img)
    735 {
    736   double* Y = NULL;
    737   double Ymax;
    738   size_t i, x, y;
    739   ASSERT(img);
    740 
    741   CHK((Y = mem_alloc(sizeof(*Y)*img->width*img->height))!= NULL);
    742 
    743   /* Copy the pixel luminance in the Y array */
    744   i = 0;
    745   FOR_EACH(y, 0, img->height) {
    746     const double* row = (double*)(img->pixels + y*img->pitch);
    747     FOR_EACH(x, 0, img->width) {
    748       const double* pixel = row + x*PIXCPNTS_COUNT__;
    749       Y[i] = pixel[PIXCPNT_Y];
    750       i++;
    751     }
    752   }
    753 
    754   /* Sort the Pixel luminance in ascending order */
    755   qsort(Y, img->width*img->height, sizeof(*Y), cmp_dbl);
    756 
    757   /* Arbitrarly use the luminance at 99.5% of the overal pixel at the maxium
    758    * luminance of the image */
    759   i = (size_t)((double)(img->width*img->height)*0.995);
    760   Ymax = Y[i];
    761 
    762   mem_rm(Y);
    763   return Ymax;
    764 }
    765 
    766 static uint8_t
    767 rgb_to_c256(const uint8_t rgb[3])
    768 {
    769   uint8_t c256 = 0;
    770   ASSERT(rgb);
    771 
    772   if(rgb[0] == rgb[1] && rgb[1] == rgb[2]) {
    773     int c = rgb[0];
    774     if(c < 4) {
    775       c256 = 16; /* Grey 0 */
    776     } else if(c >= 243) {
    777       c256 = 231; /* Grey 100 */
    778     } else {
    779       c256 = (uint8_t)(232 + (c-8+5) / 10);
    780     }
    781   } else {
    782     int r, g, b;
    783     r = rgb[0] < 48 ? 0 : (rgb[0] < 95 ? 1 : 1 + (rgb[0] - 95 + 20) / 40);
    784     g = rgb[1] < 48 ? 0 : (rgb[1] < 95 ? 1 : 1 + (rgb[1] - 95 + 20) / 40);
    785     b = rgb[2] < 48 ? 0 : (rgb[2] < 95 ? 1 : 1 + (rgb[2] - 95 + 20) / 40);
    786     c256 = (uint8_t)(36*r + 6*g + b + 16);
    787   }
    788   return c256;
    789 }
    790 
    791 static void
    792 print_color_map(const struct scmap* scmap, const double range[2])
    793 {
    794   const double ransz = range[1] - range[0];
    795   const int map_length = 65;
    796   const int map_quarter = map_length / 4;
    797   const int label_length = map_length / 4;
    798   int i;
    799   ASSERT(range && range[0] <= range[1]);
    800 
    801   FOR_EACH(i, 0, map_length) {
    802     const double u = (double)i / (double)(map_length-1);
    803     double color[3] = {0,0,0};
    804     uint8_t rgb[3];
    805     uint8_t c256;
    806     SCMAP(fetch_color(scmap, u, SCMAP_FILTER_LINEAR, color));
    807 
    808     rgb[0] = (uint8_t)(CLAMP(color[0], 0, 1) * 255. + 0.5/*round*/);
    809     rgb[1] = (uint8_t)(CLAMP(color[1], 0, 1) * 255. + 0.5/*round*/);
    810     rgb[2] = (uint8_t)(CLAMP(color[2], 0, 1) * 255. + 0.5/*round*/);
    811     c256 = rgb_to_c256(rgb);
    812     if(i == 0 * map_quarter
    813     || i == 1 * map_quarter
    814     || i == 2 * map_quarter
    815     || i == 3 * map_quarter
    816     || i == 4 * map_quarter) {
    817       fprintf(stderr, "\x1b[0m|");
    818     } else {
    819       fprintf(stderr, "\x1b[48;5;%dm ", c256);
    820     }
    821   }
    822   fprintf(stderr, "\n");
    823   fprintf(stderr, "%-*.5g", label_length, range[0]);
    824   fprintf(stderr, "%-*.5g", label_length, 0.25 * ransz + range[0]);
    825   fprintf(stderr, "%-*.5g", label_length, 0.50 * ransz + range[0]);
    826   fprintf(stderr, "%-*.5g", label_length, 0.75 * ransz + range[0]);
    827   fprintf(stderr, "%-*.5g", label_length, range[1]);
    828   fprintf(stderr, "\n");
    829 }
    830 
    831 static res_T
    832 pp_map(struct img* img, const struct args* args)
    833 {
    834   struct scmap* scmap = NULL;
    835   int64_t i;
    836   double range[2];
    837   double ransz;
    838   res_T res = RES_OK;
    839   ASSERT(img && args && args->pp_type == PP_MAP);
    840 
    841   res = scmap_create(NULL, NULL, 1, args->map.palette, &scmap);
    842   if(res != RES_OK) {
    843     fprintf(stderr, "Could not create the color map -- %s.\n",
    844       res_to_cstr(res));
    845     goto error;
    846   }
    847 
    848   if(args->map.range[0] < args->map.range[1]) {
    849     /* The range is fixed */
    850     range[0] = args->map.range[0];
    851     range[1] = args->map.range[1];
    852   } else {
    853     /* The range is defined from the loaded data  */
    854     range[0] = img->ranges[args->map.pixcpnt][0];
    855     range[1] = img->ranges[args->map.pixcpnt][1];
    856   }
    857   ransz = range[1] - range[0];
    858 
    859   omp_set_num_threads(args->nthreads);
    860   #pragma omp parallel for
    861   for(i=0; i < (int64_t)(img->width*img->height); ++i) {
    862     const size_t y = (size_t)i / img->width;
    863     const size_t x = (size_t)i % img->width;
    864     double* row = (double*)(img->pixels + img->pitch*y);
    865     double* pixel = row + x*PIXCPNTS_COUNT__;
    866     if(ransz == 0) {
    867       pixel[PIXCPNT_R] = 0;
    868       pixel[PIXCPNT_G] = 0;
    869       pixel[PIXCPNT_B] = 0;
    870     } else {
    871       const double val = CLAMP((pixel[args->map.pixcpnt] - range[0])/ransz, 0, 1);
    872       double color[3] = {0,0,0};
    873 
    874       SCMAP(fetch_color(scmap, val, SCMAP_FILTER_LINEAR, color));
    875       pixel[PIXCPNT_R] = color[0];
    876       pixel[PIXCPNT_G] = color[1];
    877       pixel[PIXCPNT_B] = color[2];
    878     }
    879   }
    880 
    881   if(args->verbose) {
    882     print_color_map(scmap, range);
    883   }
    884 
    885 exit:
    886   if(scmap) SCMAP(ref_put(scmap));
    887   return res;
    888 error:
    889   goto exit;
    890 }
    891 
    892 static res_T
    893 pp_image(struct img* img, const struct args* args)
    894 {
    895   int64_t i;
    896   double Ymax;
    897   res_T res = RES_OK;
    898   ASSERT(img && args && args->pp_type == PP_IMAGE);
    899 
    900   if(args->image.white> 0) {
    901     Ymax = args->image.white;
    902   } else {
    903     Ymax = compute_XYZ_normalization_factor(img);
    904     if(args->verbose) {
    905       fprintf(stderr, "White scale = %g\n", Ymax);
    906     }
    907   }
    908 
    909   omp_set_num_threads(args->nthreads);
    910   #pragma omp parallel for
    911   for(i=0; i < (int64_t)(img->width*img->height); ++i) {
    912     const size_t y = (size_t)i / img->width;
    913     const size_t x = (size_t)i % img->width;
    914     double* row = (double*)(img->pixels + img->pitch*y);
    915     double* pixel = row + x*PIXCPNTS_COUNT__;
    916 
    917     filmic_tone_mapping(pixel, args->image.exposure, Ymax);
    918     XYZ_to_sRGB(pixel);
    919     pixel[PIXCPNT_R] = sRGB_gamma_correct(pixel[PIXCPNT_R]);
    920     pixel[PIXCPNT_G] = sRGB_gamma_correct(pixel[PIXCPNT_G]);
    921     pixel[PIXCPNT_B] = sRGB_gamma_correct(pixel[PIXCPNT_B]);
    922   }
    923 
    924   return res;
    925 }
    926 
    927 /*******************************************************************************
    928  * Program
    929  ******************************************************************************/
    930 int
    931 main(int argc, char** argv)
    932 {
    933   FILE* stream_out = stdout;
    934   FILE* stream_in = stdin;
    935   const char* stream_out_name = "stdout";
    936   const char* stream_in_name = "stdin";
    937   struct img img = IMG_NULL;
    938   struct args args = ARGS_DEFAULT__;
    939   int img_is_loaded = 0;
    940   int err = 0;
    941   res_T res = RES_OK;
    942 
    943   res = args_init(&args, argc, argv);
    944   if(res != RES_OK) goto error;
    945   if(args.quit) goto exit;
    946 
    947   if(args.output) {
    948     res = open_output_stream(args.output, args.force_overwrite, &stream_out);
    949     if(res != RES_OK) goto error;
    950     stream_out_name = args.output;
    951   }
    952   if(args.input) {
    953     res = open_input_stream(args.input, &stream_in);
    954     if(res != RES_OK) goto error;
    955     stream_in_name = args.input;
    956   } else if(args.verbose) {
    957     fprintf(stderr, "Read image from standard input.\n");
    958   }
    959 
    960   res = img_load(&img, stream_in, stream_in_name);
    961   if(res != RES_OK) goto error;
    962   img_is_loaded = 1;
    963 
    964   switch(args.pp_type) {
    965     case PP_IMAGE:
    966       res = pp_image(&img, &args);
    967       if(res != RES_OK) goto error;
    968       res = img_write_ppm(&img, stream_out, stream_out_name);
    969       if(res != RES_OK) goto error;
    970       break;
    971     case PP_MAP:
    972       if(args.map.gnuplot) {
    973         img_write_gnuplot(&img, &args, stream_out, stream_out_name);
    974         if(res != RES_OK) goto error;
    975       } else {
    976         res = pp_map(&img, &args);
    977         if(res != RES_OK) goto error;
    978         res = img_write_ppm(&img, stream_out, stream_out_name);
    979         if(res != RES_OK) goto error;
    980       }
    981       break;
    982     default: FATAL("Unreachable code.\n"); break;
    983   }
    984 
    985 exit:
    986   if(stream_out && stream_out != stdout) fclose(stream_out);
    987   if(stream_in && stream_in != stdin) fclose(stream_in);
    988   if(img_is_loaded) img_release(&img);
    989   args_release(&args);
    990   if(mem_allocated_size() != 0) {
    991     fprintf(stderr, "Memory leaks: %lu Bytes.\n",
    992       (unsigned long)mem_allocated_size());
    993   }
    994   return err;
    995 error:
    996   err = -1;
    997   goto exit;
    998 }
    999