// Copyright (c) 2011 Public Software Group e. V., Berlin, Germany
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.

/*

hdrcompose is a simple command line utility to create HDR images.

Version: 0.2

Dependencies:
  ImageMagick  available at http://www.imagemagick.org/
  libtiff      available at http://www.libtiff.org/

Compile with:
gcc -O3 -ffast-math -march=native -I/usr/local/include -L/usr/local/lib -lm -ltiff -o hdrcompose hdrcompose.c `MagickWand-config --cflags --cppflags --ldflags --libs`

Replace /usr/local/include and /usr/local/lib as neccessary for libtiff.


Usage: hdrcompose <curve> <file1> <ev1> [<file2> <ev2> ...] <output_tiff_file>

Supported curves: "linear", "srgb", "adobe"
If you are unsure about which response curve to use, then try "srgb".

Example:
hdrcompose srgb short-exposure.jpg -2 medium.jpg 0 long-exposure.jpg 2 hdr.tiff

*/


#include <stdlib.h>
#include <stdint.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <wand/MagickWand.h>
#include <tiffio.h>

#define limiter(x) ( \
    ((x) > 1.0) ? 0.75 : \
    ((x) > 0.5) ? (-(x)*(x)+2*(x)-0.25) : (x) \
  )

#define scaled_limiter(x, f) ( limiter((x) * (f)) / (f) )

double remove_gamma_linear(double v) {
  return v;
}

double remove_gamma_srgb(double v) {
  if (v <= 0.03928) return v / 12.92;
  else return pow((v + 0.055) / 1.055, 2.4);
}

double remove_gamma_adobe(double v) {
  return pow(v, 2.19921875);
}

typedef double (*remove_gamma_fptr)(double);

int main(int argc, char **argv) {
  int input_file_count = 0;
  int width, height;
  char **filenames;
  double *ev_factors;
  float *output_image;
  MagickWand *input_image;
  MagickBooleanType magick_status;
  float *input_row_buffer;
  int i, j;
  size_t x, y;
  remove_gamma_fptr remove_gamma;
  TIFF* tif;
  if (argc < 5 || (argc-3) % 2) {
    char *command;
    if (argc >= 1) command = argv[0];
    else command = "executable";
    fprintf(stderr, "Usage: %s <curve> <file1> <ev1> [<file2> <ev2> ...] <output_tiff_file>\n", command);
    fprintf(stderr, "Supported curves: \"linear\", \"srgb\", \"adobe\"\n");
    exit(1);
  }
  if      (!strcasecmp(argv[1], "linear")) remove_gamma = remove_gamma_linear;
  else if (!strcasecmp(argv[1], "srgb"  )) remove_gamma = remove_gamma_srgb;
  else if (!strcasecmp(argv[1], "adobe" )) remove_gamma = remove_gamma_adobe;
  else {
    fprintf(stderr, "Curve \"%s\" is not supported.\n", argv[1]);
    fprintf(stderr, "Supported curves: \"linear\", \"srgb\", \"adobe\"\n");
    exit(1);
  }
  input_file_count = (argc-3) / 2;
  filenames = malloc(input_file_count * sizeof(char *));
  if (!filenames) {
    fprintf(stderr, "Failed to allocate memory.\n");
    exit(1);
  }
  ev_factors = malloc(input_file_count * sizeof(double));
  if (!ev_factors) {
    fprintf(stderr, "Failed to allocate memory.\n");
    exit(1);
  }
  for (i=0; i<input_file_count; i++) {
    char *ev_string;
    double ev_value;
    char *endptr;
    ev_string = argv[2+2*i+0];
    filenames[i] = ev_string;
    ev_value = strtod(argv[2+2*i+1], &endptr);
    if (*endptr) {
      fprintf(stderr, "Invalid floating point number: \"%s\"\n", ev_string);
      exit(1);
    }
    ev_factors[i] = pow(2.0, ev_value);
  }
  for (i=0; i<input_file_count-1; i++) {
    for (j=i+1; j<input_file_count; j++) {
      char *tmp1;
      double tmp2;
      if (ev_factors[i] > ev_factors[j]) {
        tmp1=filenames[i];  filenames[i] =filenames[j];  filenames[j] =tmp1;
        tmp2=ev_factors[i]; ev_factors[i]=ev_factors[j]; ev_factors[j]=tmp2;
      }
    }
  }
  MagickWandGenesis();
  for (i=0; i<input_file_count; i++) {
    int current_width, current_height;
    input_image = NewMagickWand();
    magick_status = MagickReadImage(input_image, filenames[i]);
    if (magick_status == MagickFalse) {
      char *description;
      ExceptionType severity;
      description = MagickGetException(input_image, &severity);
      fprintf(stderr, "Could not read file \"%s\": %s\n", filenames[i], description);
      exit(1);
    }
    current_width  = MagickGetImageWidth(input_image);
    current_height = MagickGetImageHeight(input_image);
    // Avoid integer overflow: 13377^2 * 3 * sizeof(float) < 2^31-1
    if (current_width > 13377 || current_height > 13377) {
      fprintf(stderr, "Image is too big.\n");
    }
    if (i == 0) {
      width  = current_width;
      height = current_height;
      input_row_buffer = malloc(3 * width * sizeof(float));
      if (!input_row_buffer) {
        fprintf(stderr, "Failed to allocate memory.\n");
        exit(1);
      }
      output_image = malloc(3 * width * height * sizeof(float));
      if (!output_image) {
        fprintf(stderr, "Failed to allocate memory.\n");
        exit(1);
      }
      for (x=0; x<3*width*height; x++) output_image[x] = 0.0;
    } else {
      if (current_width != width || current_height != height) {
        fprintf(stderr, "File dimensions are not matching.\n");
        exit(1);
      }
    }
    for (y=0; y<height; y++) {
      magick_status = MagickExportImagePixels(
        input_image,
        0, y, width, 1,
        "RGB", FloatPixel, input_row_buffer
      );
      if (magick_status == MagickFalse) {
        char *description;
        ExceptionType severity;
        description = MagickGetException(input_image, &severity);
        fprintf(stderr, "Error in MagickExportImagePixels call: %s\n", description);
        exit(1);
      }
      for (x=0; x<3*width; x++) {
        double value;
        value = remove_gamma(input_row_buffer[x]) / ev_factors[i];
        if (i) value = scaled_limiter(value, ev_factors[i]);
        if (i < input_file_count-1) {
          value -= scaled_limiter(value, ev_factors[i+1]);
        }
        output_image[3*width*y+x] += value;
      }
    }
    DestroyMagickWand(input_image);
  }
  MagickWandTerminus();
  free(input_row_buffer);
  free(ev_factors);
  free(filenames);
  tif = TIFFOpen(argv[argc-1], "w");
  if (!tif) exit(1);
  if (!TIFFSetField(tif, TIFFTAG_IMAGEWIDTH, (uint32_t)width)      ) exit(1);
  if (!TIFFSetField(tif, TIFFTAG_IMAGELENGTH, (uint32_t)height)    ) exit(1);
  if (!TIFFSetField(tif, TIFFTAG_ROWSPERSTRIP, (uint32_t)height)   ) exit(1);
  if (!TIFFSetField(tif, TIFFTAG_BITSPERSAMPLE,
    (uint16_t)8*(uint16_t)sizeof(float))                           ) exit(1);
  if (!TIFFSetField(tif, TIFFTAG_COMPRESSION, COMPRESSION_DEFLATE) ) exit(1);
  if (!TIFFSetField(tif, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB)     ) exit(1);
  if (!TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT, SAMPLEFORMAT_IEEEFP)) exit(1);
  if (!TIFFSetField(tif, TIFFTAG_SAMPLESPERPIXEL, (uint16_t)3)     ) exit(1);
  if (!TIFFSetField(tif, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG)) exit(1);
  if (TIFFWriteEncodedStrip(
    tif, 0, output_image,
    3 * width * height * sizeof(float)
  ) < 0) exit(1);
  TIFFClose(tif);
  free(output_image);
  return 0;
}