/* This file is an image processing operation for GEGL
 *
 * GEGL is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 3 of the License, or (at your option) any later version.
 *
 * GEGL is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with GEGL; if not, see <http://www.gnu.org/licenses/>.
 *
 * Copyright        2010 Danny Robson      <danny@blubinc.net>
 * (pfscalibration) 2006 Grzegorz Krawczyk <gkrawczyk@users.sourceforge.net>
 */

#include "config.h"
#include <glib/gi18n-lib.h>

#ifdef GEGL_PROPERTIES

property_string(exposures, _("Exposure values"), "")
    description (_("Relative brightness of each exposure in EV"))

property_int (steps, _("Discretization bits"), 13)
    description (_("Log2 of source's discretization steps"))
    value_range (8, 32)

property_double (sigma, _("Weight sigma"), 8.0f)
    description (_("Weight distribution sigma controlling response contributions"))
    value_range (0.0f, 32.0f)

#else

/*#define DEBUG_LINEARIZE*/
/*#define DEBUG_SAVE_CURVES*/

#include <gegl-plugin.h>
struct _GeglOp
{
  GeglOperationFilter parent_instance;
  gpointer            properties;
};

typedef struct
{
  GeglOperationFilterClass parent_class;
} GeglOpClass;


#define GEGL_OP_C_SOURCE exp-combine.c
#define GEGL_OP_NAME     exp_combine
#include "gegl-op.h"
GEGL_DEFINE_DYNAMIC_OPERATION(GEGL_TYPE_OPERATION_FILTER)

#include <errno.h>
#include <math.h>
#include <stdio.h>

#include "gegl-debug.h"
#include "graph/gegl-node-private.h"
#include "graph/gegl-pad.h"

static const gchar *PAD_FORMAT = "R'G'B' float";
static const gchar *EXP_PREFIX = "exposure-";

/* maximum iterations after algorithm accepts local minima */
static const guint MAXIT = 500;

/* maximum accepted error */
static const gfloat MAX_DELTA  = 1e-5f;
static const gfloat MIN_WEIGHT = 1e-3f;

/* number of pixels used if downscaling for curve estimation */
static const guint  MAX_SCALED_PIXELS = 1000 * 1000;

typedef enum
{
  PIXELS_ACTIVE,    /* Must be lowest valid ID, zero */
  PIXELS_FULL,
  PIXELS_SCALED,

  NUM_PIXEL_BUCKETS
} pixel_bucket;


typedef struct _exposure
{
  /* Immediately higher and lower exposures */
  struct _exposure *hi;
  struct _exposure *lo;

  gfloat *pixels[NUM_PIXEL_BUCKETS];

  /* 'Duration' of exposure. May only be relatively correct against the
   * other exposures in sequence.
   */
  gfloat  ti;
} exposure;


static void
gegl_expcombine_exposures_set_active (GSList       *exposures,
                                      pixel_bucket  bucket)
{
  g_return_if_fail (bucket > PIXELS_ACTIVE && bucket < NUM_PIXEL_BUCKETS);

  for (; exposures; exposures = exposures->next)
    {
      exposure *e              = exposures->data;
      e->pixels[PIXELS_ACTIVE] = e->pixels[bucket];
    }
}


/* Comparator for struct exposure: order by ti */
static int
gegl_expcombine_exposure_cmp (gconstpointer _a,
                              gconstpointer _b)
{
  const exposure *a = _a,
                 *b = _b;

  if (a->ti > b->ti)
    return  1;
  else if (a->ti < b->ti)
    return -1;
  else
    return  0;
}


static exposure*
gegl_expcombine_new_exposure (void)
{
  exposure *e = g_new (exposure, 1);
  e->hi = e->lo = e;

  memset (e->pixels, 0, sizeof (e->pixels));
  e->ti = NAN;

  return e;
}


static void
gegl_expcombine_destroy_exposure (exposure *e)
{
  guint i, j;

  /* Unlink ourselves from the next hi and lo exposures */
  g_return_if_fail (e->lo);
  g_return_if_fail (e->hi);

  e->lo->hi = (e->hi == e) ? e->lo : e->hi;
  e->hi->lo = (e->lo == e) ? e->hi : e->lo;

  /* Free each pixel bucket. Each time a non-null is found, free it then null
   * any references to the same pixels in later buckets.
   */
  for (i = PIXELS_ACTIVE + 1; i < NUM_PIXEL_BUCKETS; ++i)
    {
      if (e->pixels[i])
        {
          g_free (e->pixels[i]);

          for (j = i + 1; j < NUM_PIXEL_BUCKETS; ++j)
              if (e->pixels[j] == e->pixels[i])
                  e->pixels[j] = NULL;
        }
    }

  g_free (e);
}


/* Initialise a response curve to linear. The absolute values are probably
 * irrelevant as the curve should be normalised later.
 */
static void
gegl_expcombine_response_linear (gfloat *response,
                                 guint   steps)
{
  guint i;
  for (i = 0; i < steps; ++i)
    response[i] = i / 255.0f;
}


/* Generate a weighting curve for the range of pixel values. Weights pixels in
 * the centre of the range more highly than the extremes.
 *
 * -sigma              weights the exponential arguments
 * -step_min,step_max  describes the region we provide weighting for
 */
static void
gegl_expcombine_weights_gauss (gfloat *weights,
                               guint   steps,
                               guint   step_min,
                               guint   step_max,
                               gfloat  sigma)
{
  gfloat step_mid1, step_mid2;
  gfloat weight;
  guint i;

  g_return_if_fail (weights);
  g_return_if_fail (step_min <= step_max);
  g_return_if_fail (step_min <  steps);
  g_return_if_fail (step_max <  steps);

  step_mid1 = step_min + (step_max - step_min) / 2.0f - 0.5f;
  step_mid2 = (step_mid1 - step_min) * (step_mid1 - step_min);

  for (i = 0; i < steps; ++i)
    {
      if (i < step_min || i > step_max)
        {
          weights[i] = 0.0f;
          continue;
        }

      /* gkrawczyk: that's not really a Gaussian, but equation is
       * taken from Robertson02 paper.
       */
      weight = exp ( -sigma                 *
                    ((gfloat)i - step_mid1) *
                    ((gfloat)i - step_mid1) /
                     step_mid2);

      /* ignore very low weights. we rely on these weights being set to zero
       * in the response estimation code, thus this should not be removed.
       */
      if (weight < MIN_WEIGHT)
        weights[i] = 0.0f;
      else
        weights[i] = weight;
    }
}


/* Normalise the response curve, in place, using the midpoint of the non-zero
 * region of the response curve. Returns the mid-point value of the curve.
 */
static gfloat
gegl_expcombine_normalize (gfloat *response,
                           guint   steps)
{
  guint  step_min, step_max, step_mid;
  guint  i;
  gfloat val_mid;

  g_return_val_if_fail (response, NAN);
  g_return_val_if_fail (steps > 0, NAN);

  /* Find the first and last non-zero values in response curve */
  for (step_min = 0;
       step_min < steps && response[step_min] == 0;
       ++step_min)
    ;
  for (step_max = steps - 1;
       step_max > 0 && response[step_max] == 0;
       --step_max)
    ;

  g_return_val_if_fail (step_max >= step_min, NAN);
  step_mid = step_min + (step_max - step_min) / 2;

  /* Find the non-zero mid-value of the response curve */
  val_mid = response[step_mid];
  if (val_mid == 0.0f)
  {
    /* find first non-zero middle response */
    while (step_mid < step_max && response[step_mid] == 0.0f)
        ++step_mid;
    val_mid = response[step_mid];
  }

  /* Normalize response curve values via the mid-value */
  g_return_val_if_fail (val_mid != 0.0f, 0.0f);

  for (i = 0; i < steps; ++i)
      response[i] /= val_mid;

  return val_mid;
}


/* Apply debevec model for response curve application. Does not suffer from
 * apparent discolouration in some areas of image reconstruction when
 * compared with robertson.
 */
static int
gegl_expcombine_apply_debevec  (gfloat              *hdr,
                                const guint          components,
                                GSList              *imgs,
                                gfloat             **response,
                                const gfloat        *weighting,
                                const guint          steps,
                                const GeglRectangle *extent)
{
  guint  step_min, step_max;
  guint  num_imgs    = g_slist_length (imgs);
  guint  saturated   = 0;
  guint  pixel_count = (guint)(extent->width * extent->height);
  guint  i, j;

  g_return_val_if_fail (hdr,                        G_MAXINT);
  g_return_val_if_fail (g_slist_length (imgs) > 0,  G_MAXINT);
  g_return_val_if_fail (response,                   G_MAXINT);
  g_return_val_if_fail (weighting,                  G_MAXINT);
  g_return_val_if_fail (steps > 0,                  G_MAXINT);
  g_return_val_if_fail (extent,                     G_MAXINT);
  g_return_val_if_fail (extent->width  > 0,         G_MAXINT);
  g_return_val_if_fail (extent->height > 0,         G_MAXINT);

  /* anti saturation: calculate trusted camera output range */
  for (step_min = 0, i = step_min; i <    steps; ++i)
    if (weighting[i] > 0)
      {
        step_min = i;
        break;
      }
  for (step_max = steps - 1, i = step_max; i > step_min; --i)
    if (weighting[i] > 0)
      {
        step_max = i;
        break;
      }
  g_return_val_if_fail (step_max >= step_min, G_MAXINT);

  for (j = 0; j < pixel_count; ++j)
    {
      gfloat  sum[3]  = { 0.0f, 0.0f, 0.0f },
              div     = 0.0f;
      gfloat  ti_max  = G_MINFLOAT,
              ti_min  = G_MAXFLOAT;
      gfloat  average;
      guint   white_step[3] = {0,0,0}, black_step[3] = {0,0,0};

      /* all exposures for each pixel */
      for (i = 0; i < num_imgs; ++i)
        {
          exposure *exp_i;
          guint     step[3], step_hi[3], step_lo[3];

          exp_i    = g_slist_nth_data (imgs, i);
          step[0]  = exp_i->pixels[PIXELS_ACTIVE][0 + j * components];
          step[1]  = exp_i->pixels[PIXELS_ACTIVE][1 + j * components];
          step[2]  = exp_i->pixels[PIXELS_ACTIVE][2 + j * components];
          g_return_val_if_fail (step[0] < steps && step[1] < steps && step[2] < steps, G_MAXINT);

          /* anti saturation: observe minimum exposure time at which saturated
           * value is present, and maximum exp time at which black value is
           * present
           */
          if ((step[0] > step_max ||
               step[1] > step_max ||
               step[2] > step_max) && exp_i->ti < ti_min)
            {
              white_step[0] = MIN (step[0], steps);
              white_step[1] = MIN (step[1], steps);
              white_step[2] = MIN (step[2], steps);
              ti_min        = exp_i->ti;
            }

          if ((step[0] < step_min ||
               step[1] < step_min ||
               step[2] < step_min) && exp_i->ti > ti_max)
            {
              /* This could be protected by a max (0, step), but we're using
               * unsigned ints so it would be worthless currently.
               */
              black_step[0] = step[0];
              black_step[1] = step[1];
              black_step[2] = step[2];
              ti_max        = exp_i->ti;
            }

          /* anti ghosting: monotonous increase in time should result in
           * monotonous increase in intensity; make forward and backward check,
           * ignore value if condition not satisfied
           */
          step_lo[0] = exp_i->lo->pixels[PIXELS_ACTIVE][0 + j * components];
          step_lo[1] = exp_i->lo->pixels[PIXELS_ACTIVE][1 + j * components];
          step_lo[2] = exp_i->lo->pixels[PIXELS_ACTIVE][2 + j * components];
          g_return_val_if_fail (step_lo[0] < steps &&
                                step_lo[1] < steps &&
                                step_lo[2] < steps, G_MAXINT);

          step_hi[0] = exp_i->hi->pixels[PIXELS_ACTIVE][0 + j * components];
          step_hi[1] = exp_i->hi->pixels[PIXELS_ACTIVE][1 + j * components];
          step_hi[2] = exp_i->hi->pixels[PIXELS_ACTIVE][2 + j * components];
          g_return_val_if_fail (step_hi[0] < steps &&
                                step_hi[1] < steps &&
                                step_hi[2] < steps, G_MAXINT);

          if (step_lo[0] > step[0] ||
              step_lo[1] > step[1] ||
              step_lo[2] > step[2])
            {
              white_step[0] = step[0];
              white_step[1] = step[1];
              white_step[2] = step[2];
              /* TODO: This is not an fminf in luminance. Should it be? */
              /* ti_min        = fminf (exp_i->ti, ti_min); */
              continue;
            }

          if (step_hi[0] < step[0] ||
              step_hi[1] < step[1] ||
              step_hi[2] < step[2])
            {
              black_step[0] = step[0];
              black_step[1] = step[1];
              black_step[2] = step[2];
              /* TODO: This is not an fminf in luminance. Should it be? */
              /* ti_max        = fmaxf (exp_i->ti, ti_max); */
              continue;
            }

          /* Weight the addition into each channel by the average weight of
           * the channels and the exposure duration. This is the key
           * difference to robertson.
           */
          average  = weighting [step[0]] +
                     weighting [step[1]] +
                     weighting [step[2]];
          average /= 3.0f;

          sum[0] += response[0][step[0]] * average / exp_i->ti;
          sum[1] += response[1][step[1]] * average / exp_i->ti;
          sum[2] += response[2][step[2]] * average / exp_i->ti;

          div += average;
        }

      g_return_val_if_fail (sum[0] >= 0.0f, G_MAXINT);
      g_return_val_if_fail (sum[1] >= 0.0f, G_MAXINT);
      g_return_val_if_fail (sum[2] >= 0.0f, G_MAXINT);
      g_return_val_if_fail (   div >= 0.0f, G_MAXINT);
      g_return_val_if_fail (ti_max <= ti_min, G_MAXINT);

      /* anti saturation: if a meaningful representation of pixel was not
       * found, replace it with information from observed data
       */
      if (div == 0.0f)
          ++saturated;

      if (div == 0.0f && ti_max != G_MINFLOAT)
        {
          sum[0] = response[0][black_step[0]] / ti_max;
          sum[1] = response[1][black_step[1]] / ti_max;
          sum[2] = response[2][black_step[2]] / ti_max;
          div = 1.0f;
        }
      else if (div == 0.0f && ti_min != G_MAXFLOAT)
        {
          sum[0] = response[0][white_step[0]] / ti_min;
          sum[1] = response[1][white_step[1]] / ti_min;
          sum[2] = response[2][white_step[2]] / ti_min;
          div = 1.0f;
        }

      if (div == 0.0f)
        {
          hdr[0 + j * components] = 0.0f;
          hdr[1 + j * components] = 0.0f;
          hdr[2 + j * components] = 0.0f;
        }
      else
        {
          hdr[0 + j * components] = sum[0] / div;
          hdr[1 + j * components] = sum[1] / div;
          hdr[2 + j * components] = sum[2] / div;
        }
    }

  return saturated;
}


/**
 * @brief Create HDR image by applying response curve to given images
 * taken with different exposures
 *
 * @param hdr [out] HDR image
 * @param imgs      list of scene exposures
 * @param response  camera response function (array size of `steps')
 * @param weighting weighting function for camera output values (array size of `steps')
 * @param steps     number of camera output levels
 * @return          number of saturated pixels in the HDR image (0: all OK)
 */
static int
gegl_expcombine_apply_response (gfloat              *hdr,
                                const guint          offset,
                                const guint          components,
                                GSList              *imgs,
                                const gfloat        *response,
                                const gfloat        *weighting,
                                const guint          steps,
                                const GeglRectangle *extent)
{
  guint  step_min, step_max;
  guint  num_imgs  = g_slist_length (imgs);
  guint  saturated = 0;
  guint  pixel_count = (guint)(extent->width * extent->height);
  guint  i, j;

  g_return_val_if_fail (hdr,                        G_MAXINT);
  g_return_val_if_fail (g_slist_length (imgs) > 0,  G_MAXINT);
  g_return_val_if_fail (response,                   G_MAXINT);
  g_return_val_if_fail (weighting,                  G_MAXINT);
  g_return_val_if_fail (steps > 0,                  G_MAXINT);
  g_return_val_if_fail (extent,                     G_MAXINT);
  g_return_val_if_fail (extent->width  > 0,         G_MAXINT);
  g_return_val_if_fail (extent->height > 0,         G_MAXINT);

  /* anti saturation: calculate trusted camera output range */
  for (step_min = 0, i = step_min; i < steps; ++i)
    {
      if (weighting[i] > 0)
        {
          step_min = i;
          break;
        }
    }
  for (step_max = steps - 1, i = step_max; i > step_min; --i)
    {
      if (weighting[i] > 0)
        {
          step_max = i;
          break;
        }
    }

  g_return_val_if_fail (step_max >= step_min, G_MAXINT);

  for (j = 0; j < pixel_count; ++j)
    {
      gfloat  sum    = 0.0f,
              div    = 0.0f;
      gfloat  ti_max = G_MINFLOAT,
              ti_min = G_MAXFLOAT;

      /* all exposures for each pixel */
      for (i = 0; i < num_imgs; ++i)
        {
          exposure *exp_i;
          guint     step, step_hi, step_lo;

          exp_i = g_slist_nth_data (imgs, i);
          step  = exp_i->pixels[PIXELS_ACTIVE][offset + j * components];
          g_return_val_if_fail (step < steps, G_MAXINT);

          /* anti saturation: observe minimum exposure time at which saturated
           * value is present, and maximum exp time at which black value is
           * present
           */
          if (step > step_max)
              ti_min = fminf (ti_min, exp_i->ti);
          if (step < step_min)
              ti_max = fmaxf (ti_max, exp_i->ti);

          /* anti ghosting: monotonous increase in time should result in
           * monotonous increase in intensity; make forward and backward check,
           * ignore value if condition not satisfied
           */
          step_lo = exp_i->lo->pixels[PIXELS_ACTIVE][offset + j * components];
          step_hi = exp_i->hi->pixels[PIXELS_ACTIVE][offset + j * components];

          if (step_lo > step || step_hi < step)
               continue;

          sum += weighting[step] * exp_i->ti * response[step];
          div += weighting[step] * exp_i->ti * exp_i->ti;
        }

      g_return_val_if_fail (sum >= 0.0f, G_MAXINT);
      g_return_val_if_fail (div >= 0.0f, G_MAXINT);
      g_return_val_if_fail (ti_max <= ti_min, G_MAXINT);

      /* anti saturation: if a meaningful representation of pixel was not
       * found, replace it with information from observed data
       */
      if (div == 0.0f)
          ++saturated;

      if (div == 0.0f && ti_max != G_MINFLOAT)
        {
          sum = response[step_min];
          div = ti_max;
        }

      if (div == 0.0f && ti_min != G_MAXFLOAT)
        {
          sum = response[step_max];
          div = ti_min;
        }

      if (div != 0.0f)
          hdr[offset + j * components] = sum / div;
      else
          hdr[offset + j * components] = 0.0f;
    }

  return saturated;
}


/**
 * @brief Calculate camera response using Robertson02 algorithm
 *
 * @param luminance [out] estimated luminance values
 * @param imgs            list of scene exposures
 * @param response  [out] array to put response function
 * @param weighting       weights
 * @param steps           max camera output (no of discrete steps)
 * @return                number of saturated pixels in the HDR image (0: all OK)
 */

static int
gegl_expcombine_get_response (gfloat              *hdr,
                              const guint          offset,
                              const guint          components,
                              GSList              *imgs,
                              gfloat              *response,
                              const gfloat        *weighting,
                              guint                steps,
                              const GeglRectangle *extent)
{
  gfloat  *response_old;
  gfloat   delta, delta_old;

  gulong  *card;
  gfloat  *sum;

  guint    i, j, hits;
  guint    iterations;
  guint    pixel_count = (guint)(extent->height * extent->width);
  gulong   saturated  = 0;
  gboolean converged;

  g_return_val_if_fail (hdr,                        G_MAXINT);
  g_return_val_if_fail (g_slist_length (imgs) > 1,  G_MAXINT);
  g_return_val_if_fail (response,                   G_MAXINT);
  g_return_val_if_fail (weighting,                  G_MAXINT);
  g_return_val_if_fail (steps > 0,                  G_MAXINT);
  g_return_val_if_fail (extent,                     G_MAXINT);
  g_return_val_if_fail (extent->width  > 0,         G_MAXINT);
  g_return_val_if_fail (extent->height > 0,         G_MAXINT);

  response_old = g_new (gfloat, steps);

  /* 0. Initialization */
  gegl_expcombine_normalize (response, steps);
  for (i = 0; i < steps; ++i)
      response_old[i] = response[i];

  saturated = gegl_expcombine_apply_response (hdr, offset, components, imgs,
                                              response, weighting, steps,
                                              extent);

  converged  = FALSE;
  card       = g_new (gulong, steps);
  sum        = g_new (gfloat, steps);
  iterations = 0;
  delta_old  = 0.0f;

  /* Optimization process */
  while (!converged)
  {
    GSList *cursor = imgs;
    /* 1. Minimize with respect to response */
    for (i = 0; i < steps; ++i)
      {
        card[i] = 0;
        sum [i] = 0.0f;
      }

    for (cursor = imgs; cursor; cursor = cursor->next)
      {
        exposure *e = cursor->data;

        for (j = 0; j < pixel_count; ++j)
          {
            guint step = e->pixels[PIXELS_ACTIVE][offset + j * components];
            if (step < steps)
              {
                sum[step] += e->ti * hdr[offset + j * components];
                ++card[step];
              }
            else
              g_warning ("robertson02: m out of range: %u", step);
          }
      }

    for (i = 0; i < steps; ++i)
      {
        if (card[i] != 0)
          response[i] = sum[i] / card[i];
        else
          response[i] = 0.0f;
      }

    /* 2. Apply new response */
    gegl_expcombine_normalize (response, steps);
    saturated = gegl_expcombine_apply_response (hdr, offset, components, imgs,
                                                response, weighting, steps,
                                                extent);

    /* 3. Check stopping condition */
    delta = 0.0f;
    hits  = 0;

    for (i = 0; i < steps; ++i)
      {
        g_return_val_if_fail (response[i] >= 0, G_MAXINT);

        if (response[i] != 0.0f)
          {
            gdouble diff     = response[i] - response_old[i];
            delta           += diff * diff;
            response_old[i]  = response[i];
            ++hits;
          }
      }
    delta /= hits;


    GEGL_NOTE (GEGL_DEBUG_PROCESS,
               "exp-combine convergence, #%u delta=%f (coverage: %u%%)",
               iterations,
               delta,
               (guint)(100.0f * (gfloat)hits / steps));
    if (delta < MAX_DELTA)
      converged = TRUE;
    else if (isnan (delta) || (iterations > MAXIT && delta_old < delta))
      {
        g_warning ("exp-combine failed to converge. too noisy data in range.");
        break;
      }

    delta_old = delta;
    ++iterations;
  }

  g_free (response_old);
  g_free (card);
  g_free (sum);

  return saturated;
}


/* All exposure pads have the common prefix, followed by a numeric
 * identifier.
 */
static inline gboolean
gegl_expcombine_is_exposure_padname (const gchar *p)
{
  return g_str_has_prefix (p, EXP_PREFIX);
}


static inline gboolean
gegl_expcombine_is_exposure_pad (GeglPad *p)
{
  return gegl_expcombine_is_exposure_padname (gegl_pad_get_name (p));
}


/* XXX: We create a large quantity of pads for all the exposures, hoping we
 * create sufficient numbers and they don't utilise too many resources. This
 * is a bad hack, but works for the time being...
 */
static void
gegl_expcombine_attach (GeglOperation *operation)
{
  GParamSpec  *pspec;
  gchar       padname[16];
  guint       i;

  pspec = g_param_spec_object ("output",
                               "output",
                               "Output buffer",
                               GEGL_TYPE_BUFFER,
                               G_PARAM_READWRITE |
                               GEGL_PARAM_PAD_OUTPUT);

  gegl_operation_create_pad (operation, pspec);
  g_param_spec_sink (pspec);

  for (i = 0; i <= 99; ++i)
    {
      snprintf (padname, G_N_ELEMENTS (padname), "exposure_%u", i);

      pspec = g_param_spec_object (padname,
                                   padname,
                                   "Exposure input.",
                                   GEGL_TYPE_BUFFER,
                                   G_PARAM_READWRITE |
                                   GEGL_PARAM_PAD_INPUT);

      gegl_operation_create_pad (operation, pspec);
      g_param_spec_sink (pspec);
    }
}



static void
gegl_expcombine_prepare (GeglOperation *operation)
{
  GSList *inputs = gegl_node_get_input_pads (operation->node);

  /* Set all the pads output formats */
  for (; inputs; inputs = inputs->next)
    {
      GeglPad *pad = inputs->data;

      gegl_pad_set_format (pad, babl_format (PAD_FORMAT));
    }

  gegl_operation_set_format (operation, "output", babl_format (PAD_FORMAT));
}


/* Order pads beginning with EXP_PREFIX by their trailing exposure value.
 * Non-exposure pads come first, then exposure pads without an EV, then
 * ordered exposure pads.
 */
static int
gegl_expcombine_pad_cmp (gconstpointer _a, gconstpointer _b)
{
  const gchar *a = gegl_pad_get_name (GEGL_PAD (_a)),
              *b = gegl_pad_get_name (GEGL_PAD (_b));
  guint64      x, y;

  if (!g_str_has_prefix (b, EXP_PREFIX)) return  1;
  if (!g_str_has_prefix (a, EXP_PREFIX)) return -1;

  a = strrchr (a, '-');
  b = strrchr (b, '-');

  g_return_val_if_fail (b, -1);
  g_return_val_if_fail (a, -1);

  y = g_ascii_strtoll (b + 1, NULL, 10);
  if (errno) return  1;
  x = g_ascii_strtoll (a + 1, NULL, 10);
  if (errno) return -1;

  if (x < y) return -1;
  if (x > y) return  1;
  return 0;
}


/* Extract a list of exposure pixel and metadata from the operation context.
 * The data is not guaranteed to be ordered within the list. May return NULL
 * for error.
 */
static GSList *
gegl_expcombine_get_exposures (GeglOperation        *operation,
                               GeglOperationContext *context,
                               const GeglRectangle  *full_roi,
                               GeglRectangle        *scaled_roi)
{
  GeglProperties    *o          = GEGL_PROPERTIES (operation);
  gchar         *ev_cursor  = o->exposures;
  GSList        *exposures  = NULL,
                *inputs,
                *cursor;
  guint          components = babl_format_get_n_components (babl_format (PAD_FORMAT));
  gfloat         scale;

  g_return_val_if_fail (operation, NULL);
  g_return_val_if_fail (context, NULL);
  g_return_val_if_fail (full_roi, NULL);
  g_return_val_if_fail (!gegl_rectangle_is_empty (full_roi), NULL);
  g_return_val_if_fail (scaled_roi, NULL);

  /* Calculate the target dimensions for the downscaled buffer used in the
   * curve recovery. Does not enlarge the buffer.
   * sqrt does not cause issues with rectangular buffers as it's merely a
   * scaling factor for later.
   */
  {
    scale       = MAX_SCALED_PIXELS /
                  (gfloat)(full_roi->width * full_roi->height);
    scale       = sqrtf (scale);
    *scaled_roi = *full_roi;

    if (scale > 1.0f)
      {
        scale = 1.0f;
      }
    else
      {
        scaled_roi->width  = full_roi->width  * scale;
        scaled_roi->height = full_roi->height * scale;
      }
  }

  inputs = g_slist_copy (gegl_node_get_input_pads (operation->node));
  inputs = g_slist_sort (inputs, gegl_expcombine_pad_cmp);

  /* Read in each of the exposure buffers */
  for (cursor = inputs; cursor; cursor = cursor->next)
    {
      GeglBuffer *buffer;
      GeglPad    *pad = cursor->data;
      exposure   *e;

      /* Weed out non-exposure pads */
      if (!gegl_expcombine_is_exposure_pad (pad))
        {
          if (strcmp (gegl_pad_get_name (pad), "aux"))
              g_warning ("Unexpected pad name '%s' in exp-combine",
                         gegl_pad_get_name (pad));
          continue;
        }

      /* Add exposure to list */
      buffer = (GeglBuffer*)gegl_operation_context_dup_object (context,
                                                   gegl_pad_get_name (pad));
      if (!buffer)
          continue;

      e                = gegl_expcombine_new_exposure ();
      e->pixels[PIXELS_FULL]   = g_new (gfloat, full_roi->width  *
                                                full_roi->height *
                                                components);
      gegl_buffer_get (buffer, full_roi, 1.0, babl_format (PAD_FORMAT),
                       e->pixels[PIXELS_FULL], GEGL_AUTO_ROWSTRIDE,
                       GEGL_ABYSS_NONE);

      g_return_val_if_fail (scale <= 1.0f, NULL);
      if (scale == 1.0f)
          e->pixels[PIXELS_SCALED] = e->pixels[PIXELS_FULL];
      else
        {
          e->pixels[PIXELS_SCALED] = g_new (gfloat,
                                            (scaled_roi->width  *
                                             scaled_roi->height *
                                             components));
          gegl_buffer_get (buffer, scaled_roi, scale, babl_format (PAD_FORMAT),
                           e->pixels[PIXELS_SCALED], GEGL_AUTO_ROWSTRIDE,
                           GEGL_ABYSS_NONE);
        }

      e->pixels[PIXELS_ACTIVE] = e->pixels[PIXELS_FULL];
      g_object_unref (buffer);

      /* Read the exposure time: relate APEX brightness value only as a
       * function of exposure time that is assume aperture = 1 and
       * sensitivity = 1
       */
      e->ti = g_ascii_strtod (ev_cursor, &ev_cursor);
      e->ti = 1.0f / powf (2.0f, e->ti);

      /* Krawczyk: absolute calibration: this magic multiplier is a result
       * of my research in internet plus a bit of tweaking with a luminance
       * meter. tested with canon 350d, so might be a bit off for other
       * cameras.
       *   e->ti /= 1.0592f * 11.4f / 3.125f;
       */

      /* For the sake of sanity, we scale the times to produce an image which
       * is much closer to the 'normal' 0-1 range. This prevents components
       * from going massively outside of our gamut.
       *
       * In particular this aids visual inspection and use of standard pixel
       * value comparison functions in testing.
       */
      e->ti *= 5e4f;

      if (errno)
        {
          g_warning ("Invalid exposure values specified for exp-combine");
          g_slist_free_full (exposures, (GDestroyNotify) gegl_expcombine_destroy_exposure);
          g_slist_free (inputs);
          gegl_expcombine_destroy_exposure (e);

          return NULL;
        }

      exposures = g_slist_prepend (exposures, e);
    }

  /* Link each exposure's high and low sibling pointers. Simplifies
   * anti-ghosting computations later.
   */
  exposures = g_slist_sort (exposures, gegl_expcombine_exposure_cmp);
  {
    exposure *prev = exposures->data;
    for (cursor = exposures; cursor; cursor = cursor->next)
      {
        exposure *e = cursor->data;

        prev->hi = e;
        e->lo    = prev;
        prev     = e;
      }
  }

  if (exposures == NULL)
      g_warning ("exp-combine did not find any suitable input pads");
  if (*ev_cursor != '\0')
      g_warning ("too many supplied EVs for input pads");

  g_slist_free (inputs);
  return exposures;
}


#ifdef DEBUG_SAVE_CURVES
static void
gegl_expcombine_save_curves (const gfloat *response,
                             guint         steps,
                             const gchar  *path)
{
  guint i;
  FILE *fd = fopen (path, "w");

  for (i = 0; i < steps; ++i)
    fprintf (fd, "%u %f\n", i, response[i]);
  fclose (fd);

  return;
}
#endif


static gboolean
gegl_expcombine_process (GeglOperation        *operation,
                         GeglOperationContext *context,
                         const gchar          *output_pad,
                         const GeglRectangle  *full_roi,
                         gint                  level)
{
  GeglProperties *o           = GEGL_PROPERTIES (operation);
  GeglBuffer *output      = gegl_operation_context_get_target (context,
                                                               output_pad);

  GeglRectangle   scaled_roi;
  GSList         *cursor;
  GSList         *exposures   = gegl_expcombine_get_exposures (operation,
                                                               context,
                                                               full_roi,
                                                               &scaled_roi);

  guint   i;
  guint   steps      = 1 << o->steps,
          step_max   = 0,
          step_min   = steps - 1;

  guint   components = babl_format_get_n_components (babl_format (PAD_FORMAT));
  gfloat *hdr        = g_new (gfloat, full_roi->width  *
                                      full_roi->height *
                                      components);
  gfloat *weights     =   g_new (gfloat, steps);
  gfloat *response[3] = { g_new (gfloat, steps),
                          g_new (gfloat, steps),
                          g_new (gfloat, steps) };

  guint   saturated  = 0;
  gfloat  over       = 0.0f,
          under      = 0.0f;

  g_return_val_if_fail (output,    FALSE);
  g_return_val_if_fail (exposures, FALSE);

  g_return_val_if_fail (steps      > 0, FALSE);
  g_return_val_if_fail (components > 0, FALSE);

  /* Find the highest and lowest valid steps in the output. Remap from the
   * 'normal' 0.0-1.0 floats to the range of integer steps.
   */
  for (cursor = exposures; cursor; cursor = cursor->next)
    {
      exposure *e = cursor->data;

      for (i = 0; i < full_roi->width * full_roi->height * components; ++i)
        {
          /* Clamp the values we receive to [0.0, 1.0) and record the
           * magnitude of this over/underflow.
           */
          if (e->pixels[PIXELS_FULL][i] <= 0.0f)
            {
              under += fabs (e->pixels[PIXELS_FULL][i]);
              e->pixels[PIXELS_FULL][i] = 0.0f;
            }
          else if (e->pixels[PIXELS_FULL][i] > 1.0f)
            {
              over  += e->pixels[PIXELS_FULL][i] - 1.0f;
              e->pixels[PIXELS_FULL][i] = 1.0f;
            }

          e->pixels[PIXELS_FULL][i] *= (steps - 1);
          if (e->pixels[PIXELS_FULL][i] > 0.0f)
            {
              step_max = MAX (step_max, e->pixels[PIXELS_FULL][i]);
              step_min = MIN (step_min, e->pixels[PIXELS_FULL][i]);
            }
        }

      if (e->pixels[PIXELS_SCALED] == e->pixels[PIXELS_FULL])
          continue;

      for (i = 0; i < scaled_roi.width * scaled_roi.height * components; ++i)
        {
          e->pixels[PIXELS_SCALED][i]  = CLAMP (e->pixels[PIXELS_SCALED][i], 0.0f, 1.0f);
          e->pixels[PIXELS_SCALED][i] *= (steps - 1);
        }

    }

  g_return_val_if_fail (step_max >= step_min, FALSE);
  if (under || over)
      g_warning ("Unexpected pixel bounds. "
                 "Overflow error: %f, underflow error: %f",
                 over  / components,
                 under / components);

  /* Initialise response curves and weights for estimation. Estimate and
   * retrieve the converted output. Apply the response curves to the input
   * using the debevec model rather than using the robertson output.
   */
  gegl_expcombine_weights_gauss (weights, steps, step_min, step_max, o->sigma);
  gegl_expcombine_exposures_set_active (exposures, PIXELS_SCALED);

  for (i = 0; i < components; ++i)
    {
      gegl_expcombine_response_linear (response[i], steps);
      saturated += gegl_expcombine_get_response (hdr, i, components, exposures,
                                                 response[i], weights, steps,
                                                 &scaled_roi);
    }

#ifdef DEBUG_SAVE_CURVES
  gegl_expcombine_save_curves (response[0], steps, "response_r");
  gegl_expcombine_save_curves (response[1], steps, "response_g");
  gegl_expcombine_save_curves (response[2], steps, "response_b");
#endif

  gegl_expcombine_exposures_set_active (exposures, PIXELS_FULL);
  gegl_expcombine_apply_debevec (hdr, components, exposures, response, weights,
                            steps, full_roi);

  g_return_val_if_fail (G_N_ELEMENTS (response) == 3, FALSE);
  g_free (response[0]);
  g_free (response[1]);
  g_free (response[2]);
  g_free (weights);

  if (saturated)
    {
      GEGL_NOTE (GEGL_DEBUG_PROCESS,
                 "exp-combine discovered %u saturated pixels (%f%%), "
                 "results may be suboptimal", saturated,
                 100.0f * saturated / full_roi->width / full_roi->height);
    }

  /* Rescale the output so we can more easily visually inspect the output
   * using LDR outputs.
   */
#ifdef DEBUG_LINEARIZE
  {
    gfloat max = G_MINFLOAT;
    guint  i;
    for (i = 0; i < full_roi->width * full_roi->height * components; ++i)
        max = MAX (max, hdr[i]);

    for (i = 0; i < full_roi->width * full_roi->height * components; ++i)
        hdr[i] /= max;
  }
#endif

  /* Save the HDR components to the output buffer. */
  gegl_buffer_set (output, full_roi, 0, babl_format (PAD_FORMAT), hdr,
                   GEGL_AUTO_ROWSTRIDE);

  /* Cleanup */
  g_free (hdr);
  g_slist_free_full (exposures, (GDestroyNotify) gegl_expcombine_destroy_exposure);

  return TRUE;
}


/* Grab each source's dimensions, and update a running bounding box with the
 * union of itself and the current source.
 * TODO: properly support inputs of varying dimensions
 */
static GeglRectangle
gegl_expcombine_get_bounding_box (GeglOperation *operation)
{
  GeglRectangle  result = {0,0,0,0};
  GSList        *inputs = gegl_node_get_input_pads (operation->node);

  for (; inputs; inputs = inputs->next)
    {
      GeglPad             *pad = inputs->data;
      const GeglRectangle *newrect;

      if (!gegl_expcombine_is_exposure_pad (pad))
        continue;

      /* Get the source bounds and update with the union */
      newrect = gegl_operation_source_get_bounding_box (operation,
                                                        gegl_pad_get_name (pad));
      if (!newrect)
        continue;

      if (!gegl_rectangle_is_empty (&result) &&
          !gegl_rectangle_equal (newrect, &result))
        {
          g_warning ("expcombine inputs are of varying dimensions");
        }
      gegl_rectangle_bounding_box (&result, newrect, &result);
    }

  if (gegl_rectangle_is_empty (&result))
      g_warning ("Bounding box for exp-combine should not be empty");
  return result;
}


static GeglRectangle
gegl_expcombine_get_cached_region (GeglOperation       *operation,
                                   const GeglRectangle *roi)
{
  return gegl_expcombine_get_bounding_box (operation);
}


/* We have to use all the input to produce any output */
static GeglRectangle
gegl_expcombine_get_required_for_output (GeglOperation       *operation,
                                         const gchar         *input_pad,
                                         const GeglRectangle *output_roi)
{
  return gegl_expcombine_get_bounding_box (operation);
}

static void
gegl_op_class_init (GeglOpClass *klass)
{
  GeglOperationClass       *operation_class;

  operation_class = GEGL_OPERATION_CLASS (klass);

  operation_class->attach            = gegl_expcombine_attach;
  operation_class->process           = gegl_expcombine_process;
  operation_class->get_bounding_box  = gegl_expcombine_get_bounding_box;
  operation_class->get_cached_region = gegl_expcombine_get_cached_region;

  operation_class->prepare     = gegl_expcombine_prepare;
  operation_class->get_required_for_output = gegl_expcombine_get_required_for_output;

  gegl_operation_class_set_keys (operation_class,
  "name"       , "gegl:exp-combine",
  "title",       _("Combine Exposures"),
  "categories" , "compositors",
  "description",
      _("Combine multiple scene exposures into one high dynamic range image."),
      NULL);
}


#endif

