Adding upstream version 2.10.22.upstream/2.10.22upstream

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

1 files changed, 840 insertions, 0 deletions

diff --git a/plug-ins/common/contrast-retinex.c b/plug-ins/common/contrast-retinex.c
new file mode 100644
index 0000000..f48e5c9
--- /dev/null
+++ b/plug-ins/common/contrast-retinex.c
@@ -0,0 +1,840 @@
+/* GIMP - The GNU Image Manipulation Program
+ * Copyright (C) 1995 Spencer Kimball and Peter Mattis
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program 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 General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#include "config.h"
+
+#include <string.h>
+
+#include "libgimp/gimp.h"
+#include "libgimp/gimpui.h"
+
+#include "libgimp/stdplugins-intl.h"
+
+
+#define PLUG_IN_PROC        "plug-in-retinex"
+#define PLUG_IN_BINARY      "contrast-retinex"
+#define PLUG_IN_ROLE        "gimp-contrast-retinex"
+#define MAX_RETINEX_SCALES    8
+#define MIN_GAUSSIAN_SCALE   16
+#define MAX_GAUSSIAN_SCALE  250
+#define SCALE_WIDTH         150
+#define ENTRY_WIDTH           4
+
+
+typedef struct
+{
+  gint     scale;
+  gint     nscales;
+  gint     scales_mode;
+  gfloat   cvar;
+} RetinexParams;
+
+typedef enum
+{
+  filter_uniform,
+  filter_low,
+  filter_high
+} FilterMode;
+
+/*
+  Definit comment sont repartis les
+  differents filtres en fonction de
+  l'echelle (~= ecart type de la gaussienne)
+ */
+#define RETINEX_UNIFORM 0
+#define RETINEX_LOW     1
+#define RETINEX_HIGH    2
+
+static gfloat RetinexScales[MAX_RETINEX_SCALES];
+
+typedef struct
+{
+  gint    N;
+  gfloat  sigma;
+  gdouble B;
+  gdouble b[4];
+} gauss3_coefs;
+
+
+/*
+ * Declare local functions.
+ */
+static void     query                       (void);
+static void     run                         (const gchar      *name,
+                                             gint              nparams,
+                                             const GimpParam  *param,
+                                             gint             *nreturn_vals,
+                                             GimpParam       **return_vals);
+
+/* Gimp */
+static gboolean retinex_dialog              (gint32        drawable_ID);
+static void     retinex                     (gint32        drawable_ID,
+                                             GimpPreview  *preview);
+static void     retinex_preview             (gpointer      drawable_ID,
+                                             GimpPreview  *preview);
+
+static void     retinex_scales_distribution (gfloat       *scales,
+                                             gint          nscales,
+                                             gint          mode,
+                                             gint          s);
+
+static void     compute_mean_var            (gfloat       *src,
+                                             gfloat       *mean,
+                                             gfloat       *var,
+                                             gint          size,
+                                             gint          bytes);
+/*
+ * Gauss
+ */
+static void     compute_coefs3              (gauss3_coefs *c,
+                                             gfloat        sigma);
+
+static void     gausssmooth                 (gfloat       *in,
+                                             gfloat       *out,
+                                             gint          size,
+                                             gint          rowtride,
+                                             gauss3_coefs *c);
+
+/*
+ * MSRCR = MultiScale Retinex with Color Restoration
+ */
+static void     MSRCR                       (guchar       *src,
+                                             gint          width,
+                                             gint          height,
+                                             gint          bytes,
+                                             gboolean      preview_mode);
+
+
+/*
+ * Private variables.
+ */
+static RetinexParams rvals =
+{
+  240,             /* Scale */
+  3,               /* Scales */
+  RETINEX_UNIFORM, /* Echelles reparties uniformement */
+  1.2              /* A voir */
+};
+
+static GimpPlugInInfo PLUG_IN_INFO =
+{
+  NULL,  /* init_proc  */
+  NULL,  /* quit_proc  */
+  query, /* query_proc */
+  run,   /* run_proc   */
+};
+
+MAIN ()
+
+static void
+query (void)
+{
+  static const GimpParamDef args[] =
+  {
+    { GIMP_PDB_INT32,    "run-mode",    "The run mode { RUN-INTERACTIVE (0), RUN-NONINTERACTIVE (1) }"        },
+    { GIMP_PDB_IMAGE,    "image",       "Input image (unused)"                },
+    { GIMP_PDB_DRAWABLE, "drawable",    "Input drawable"                      },
+    { GIMP_PDB_INT32,    "scale",       "Biggest scale value"                 },
+    { GIMP_PDB_INT32,    "nscales",     "Number of scales"                    },
+    { GIMP_PDB_INT32,    "scales-mode", "Retinex distribution through scales" },
+    { GIMP_PDB_FLOAT,    "cvar",        "Variance value"                      }
+  };
+
+  gimp_install_procedure (PLUG_IN_PROC,
+                          N_("Enhance contrast using the Retinex method"),
+                          "The Retinex Image Enhancement Algorithm is an "
+                          "automatic image enhancement method that enhances "
+                          "a digital image in terms of dynamic range "
+                          "compression, color independence from the spectral "
+                          "distribution of the scene illuminant, and "
+                          "color/lightness rendition.",
+                          "Fabien Pelisson",
+                          "Fabien Pelisson",
+                          "2003",
+                          N_("Retine_x..."),
+                          "RGB*",
+                          GIMP_PLUGIN,
+                          G_N_ELEMENTS (args), 0,
+                          args, NULL);
+
+  gimp_plugin_menu_register (PLUG_IN_PROC, "<Image>/Colors/Tone Mapping");
+}
+
+static void
+run (const gchar      *name,
+     gint              nparams,
+     const GimpParam  *param,
+     gint             *nreturn_vals,
+     GimpParam       **return_vals)
+{
+  static GimpParam   values[1];
+  GimpRunMode        run_mode;
+  gint32             drawable_ID;
+  GimpPDBStatusType  status = GIMP_PDB_SUCCESS;
+  gint               x, y, width, height;
+
+  INIT_I18N ();
+  gegl_init (NULL, NULL);
+
+  *nreturn_vals = 1;
+  *return_vals  = values;
+
+  values[0].type          = GIMP_PDB_STATUS;
+  values[0].data.d_status = status;
+
+  run_mode    = param[0].data.d_int32;
+  drawable_ID = param[2].data.d_drawable;
+
+  if (! gimp_drawable_mask_intersect (drawable_ID,
+                                      &x, &y, &width, &height) ||
+      width  < MIN_GAUSSIAN_SCALE ||
+      height < MIN_GAUSSIAN_SCALE)
+    {
+      status = GIMP_PDB_EXECUTION_ERROR;
+      values[0].data.d_status = status;
+      return;
+    }
+
+  switch (run_mode)
+    {
+    case GIMP_RUN_INTERACTIVE:
+      /*  Possibly retrieve data  */
+      gimp_get_data (PLUG_IN_PROC, &rvals);
+
+      /*  First acquire information with a dialog  */
+      if (! retinex_dialog (drawable_ID))
+        return;
+      break;
+
+    case GIMP_RUN_NONINTERACTIVE:
+      /*  Make sure all the arguments are there!  */
+      if (nparams != 7)
+        {
+          status = GIMP_PDB_CALLING_ERROR;
+        }
+      else
+        {
+          rvals.scale        = (param[3].data.d_int32);
+          rvals.nscales      = (param[4].data.d_int32);
+          rvals.scales_mode  = (param[5].data.d_int32);
+          rvals.cvar         = (param[6].data.d_float);
+        }
+      break;
+
+    case GIMP_RUN_WITH_LAST_VALS:
+      gimp_get_data (PLUG_IN_PROC, &rvals);
+      break;
+
+    default:
+      break;
+    }
+
+  if (status == GIMP_PDB_SUCCESS &&
+      (gimp_drawable_is_rgb (drawable_ID)))
+    {
+      gimp_progress_init (_("Retinex"));
+
+      retinex (drawable_ID, NULL);
+
+      if (run_mode != GIMP_RUN_NONINTERACTIVE)
+        gimp_displays_flush ();
+
+      /*  Store data  */
+      if (run_mode == GIMP_RUN_INTERACTIVE)
+        gimp_set_data (PLUG_IN_PROC, &rvals, sizeof (RetinexParams));
+    }
+  else
+    {
+      status = GIMP_PDB_EXECUTION_ERROR;
+    }
+
+  values[0].data.d_status = status;
+}
+
+
+static gboolean
+retinex_dialog (gint32 drawable_ID)
+{
+  GtkWidget *dialog;
+  GtkWidget *main_vbox;
+  GtkWidget *preview;
+  GtkWidget *table;
+  GtkWidget *combo;
+  GtkObject *adj;
+  gboolean   run;
+
+  gimp_ui_init (PLUG_IN_BINARY, FALSE);
+
+  dialog = gimp_dialog_new (_("Retinex Image Enhancement"), PLUG_IN_ROLE,
+                            NULL, 0,
+                            gimp_standard_help_func, PLUG_IN_PROC,
+
+                            _("_Cancel"), GTK_RESPONSE_CANCEL,
+                            _("_OK"),     GTK_RESPONSE_OK,
+
+                            NULL);
+
+  gtk_dialog_set_alternative_button_order (GTK_DIALOG (dialog),
+                                           GTK_RESPONSE_OK,
+                                           GTK_RESPONSE_CANCEL,
+                                           -1);
+
+  gimp_window_set_transient (GTK_WINDOW (dialog));
+
+  main_vbox = gtk_box_new (GTK_ORIENTATION_VERTICAL, 12);
+  gtk_container_set_border_width (GTK_CONTAINER (main_vbox), 12);
+  gtk_box_pack_start (GTK_BOX (gtk_dialog_get_content_area (GTK_DIALOG (dialog))),
+                      main_vbox, TRUE, TRUE, 0);
+  gtk_widget_show (main_vbox);
+
+  preview = gimp_zoom_preview_new_from_drawable_id (drawable_ID);
+  gtk_box_pack_start (GTK_BOX (main_vbox), preview, TRUE, TRUE, 0);
+  gtk_widget_show (preview);
+
+  g_signal_connect_swapped (preview, "invalidated",
+                            G_CALLBACK (retinex_preview),
+                            GINT_TO_POINTER (drawable_ID));
+
+  table = gtk_table_new (4, 3, FALSE);
+  gtk_table_set_col_spacings (GTK_TABLE (table), 6);
+  gtk_table_set_row_spacings (GTK_TABLE (table), 6);
+  gtk_box_pack_start (GTK_BOX (main_vbox), table, FALSE, FALSE, 0);
+  gtk_widget_show (table);
+
+  combo = gimp_int_combo_box_new (_("Uniform"), filter_uniform,
+                                  _("Low"),     filter_low,
+                                  _("High"),    filter_high,
+                                  NULL);
+
+  gimp_int_combo_box_connect (GIMP_INT_COMBO_BOX (combo), rvals.scales_mode,
+                              G_CALLBACK (gimp_int_combo_box_get_active),
+                              &rvals.scales_mode);
+  g_signal_connect_swapped (combo, "changed",
+                            G_CALLBACK (gimp_preview_invalidate),
+                            preview);
+
+  gimp_table_attach_aligned (GTK_TABLE (table), 0, 0,
+                             _("_Level:"), 0.0, 0.5,
+                             combo, 2, FALSE);
+  gtk_widget_show (combo);
+
+  adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 1,
+                              _("_Scale:"), SCALE_WIDTH, ENTRY_WIDTH,
+                              rvals.scale,
+                              MIN_GAUSSIAN_SCALE, MAX_GAUSSIAN_SCALE, 1, 1, 0,
+                              TRUE, 0, 0, NULL, NULL);
+
+  g_signal_connect (adj, "value-changed",
+                    G_CALLBACK (gimp_int_adjustment_update),
+                    &rvals.scale);
+  g_signal_connect_swapped (adj, "value-changed",
+                            G_CALLBACK (gimp_preview_invalidate),
+                            preview);
+
+  adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 2,
+                              _("Scale _division:"), SCALE_WIDTH, ENTRY_WIDTH,
+                              rvals.nscales,
+                              0, MAX_RETINEX_SCALES, 1, 1, 0,
+                              TRUE, 0, 0, NULL, NULL);
+
+  g_signal_connect (adj, "value-changed",
+                    G_CALLBACK (gimp_int_adjustment_update),
+                    &rvals.nscales);
+  g_signal_connect_swapped (adj, "value-changed",
+                            G_CALLBACK (gimp_preview_invalidate),
+                            preview);
+
+  adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 3,
+                              _("Dy_namic:"), SCALE_WIDTH, ENTRY_WIDTH,
+                              rvals.cvar, 0, 4, 0.1, 0.1, 1,
+                              TRUE, 0, 0, NULL, NULL);
+
+  g_signal_connect (adj, "value-changed",
+                    G_CALLBACK (gimp_float_adjustment_update),
+                    &rvals.cvar);
+  g_signal_connect_swapped (adj, "value-changed",
+                            G_CALLBACK (gimp_preview_invalidate),
+                            preview);
+
+  gtk_widget_show (dialog);
+
+  run = (gimp_dialog_run (GIMP_DIALOG (dialog)) == GTK_RESPONSE_OK);
+
+  gtk_widget_destroy (dialog);
+
+  return run;
+}
+
+/*
+ * Applies the algorithm
+ */
+static void
+retinex (gint32       drawable_ID,
+         GimpPreview *preview)
+{
+  GeglBuffer *src_buffer;
+  GeglBuffer *dest_buffer;
+  const Babl *format;
+  guchar     *src  = NULL;
+  guchar     *psrc = NULL;
+  gint        x, y, width, height;
+  gint        size, bytes;
+
+  /*
+   * Get the size of the current image or its selection.
+   */
+  if (preview)
+    {
+      src = gimp_zoom_preview_get_source (GIMP_ZOOM_PREVIEW (preview),
+                                          &width, &height, &bytes);
+    }
+  else
+    {
+      if (! gimp_drawable_mask_intersect (drawable_ID,
+                                          &x, &y, &width, &height))
+        return;
+
+      if (gimp_drawable_has_alpha (drawable_ID))
+        format = babl_format ("R'G'B'A u8");
+      else
+        format = babl_format ("R'G'B' u8");
+
+      bytes = babl_format_get_bytes_per_pixel (format);
+
+      /* Allocate memory */
+      size = width * height * bytes;
+      src  = g_try_malloc (sizeof (guchar) * size);
+
+      if (src == NULL)
+        {
+          g_warning ("Failed to allocate memory");
+          return;
+        }
+
+      memset (src, 0, sizeof (guchar) * size);
+
+      /* Fill allocated memory with pixel data */
+      src_buffer = gimp_drawable_get_buffer (drawable_ID);
+
+      gegl_buffer_get (src_buffer, GEGL_RECTANGLE (x, y, width, height), 1.0,
+                       format, src,
+                       GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_NONE);
+    }
+
+  /*
+    Algorithm for Multi-scale Retinex with color Restoration (MSRCR).
+   */
+  psrc = src;
+  MSRCR (psrc, width, height, bytes, preview != NULL);
+
+  if (preview)
+    {
+      gimp_preview_draw_buffer (preview, psrc, width * bytes);
+    }
+  else
+    {
+      dest_buffer = gimp_drawable_get_shadow_buffer (drawable_ID);
+
+      gegl_buffer_set (dest_buffer, GEGL_RECTANGLE (x, y, width, height), 0,
+                       format, psrc,
+                       GEGL_AUTO_ROWSTRIDE);
+
+      g_object_unref (src_buffer);
+      g_object_unref (dest_buffer);
+
+      gimp_progress_update (1.0);
+
+      gimp_drawable_merge_shadow (drawable_ID, TRUE);
+      gimp_drawable_update (drawable_ID, x, y, width, height);
+    }
+
+  g_free (src);
+}
+
+static void
+retinex_preview (gpointer     drawable_ID,
+                 GimpPreview *preview)
+{
+  retinex (GPOINTER_TO_INT (drawable_ID), preview);
+}
+
+/*
+ * calculate scale values for desired distribution.
+ */
+static void
+retinex_scales_distribution (gfloat *scales,
+                             gint    nscales,
+                             gint    mode,
+                             gint    s)
+{
+  if (nscales == 1)
+    { /* For one filter we choose the median scale */
+      scales[0] = (gint) s / 2;
+    }
+  else if (nscales == 2)
+    { /* For two filters we choose the median and maximum scale */
+      scales[0] = (gint) s / 2;
+      scales[1] = (gint) s;
+    }
+  else
+    {
+      gfloat size_step = (gfloat) s / (gfloat) nscales;
+      gint   i;
+
+      switch(mode)
+        {
+        case RETINEX_UNIFORM:
+          for(i = 0; i < nscales; ++i)
+            scales[i] = 2. + (gfloat) i * size_step;
+          break;
+
+        case RETINEX_LOW:
+          size_step = (gfloat) log(s - 2.0) / (gfloat) nscales;
+          for (i = 0; i < nscales; ++i)
+            scales[i] = 2. + pow (10, (i * size_step) / log (10));
+          break;
+
+        case RETINEX_HIGH:
+          size_step = (gfloat) log(s - 2.0) / (gfloat) nscales;
+          for (i = 0; i < nscales; ++i)
+            scales[i] = s - pow (10, (i * size_step) / log (10));
+          break;
+
+        default:
+          break;
+        }
+    }
+}
+
+/*
+ * Calculate the coefficients for the recursive filter algorithm
+ * Fast Computation of gaussian blurring.
+ */
+static void
+compute_coefs3 (gauss3_coefs *c, gfloat sigma)
+{
+  /*
+   * Papers:  "Recursive Implementation of the gaussian filter.",
+   *          Ian T. Young , Lucas J. Van Vliet, Signal Processing 44, Elsevier 1995.
+   * formula: 11b       computation of q
+   *          8c        computation of b0..b1
+   *          10        alpha is normalization constant B
+   */
+  gfloat q, q2, q3;
+
+  if (sigma >= 2.5)
+    {
+      q = 0.98711 * sigma - 0.96330;
+    }
+  else if ((sigma >= 0.5) && (sigma < 2.5))
+    {
+      q = 3.97156 - 4.14554 * (gfloat) sqrt ((double) 1 - 0.26891 * sigma);
+    }
+  else
+    {
+      q = 0.1147705018520355224609375;
+    }
+
+  q2 = q * q;
+  q3 = q * q2;
+  c->b[0] = (1.57825+(2.44413*q)+(1.4281 *q2)+(0.422205*q3));
+  c->b[1] = (        (2.44413*q)+(2.85619*q2)+(1.26661 *q3));
+  c->b[2] = (                   -((1.4281*q2)+(1.26661 *q3)));
+  c->b[3] = (                                 (0.422205*q3));
+  c->B = 1.0-((c->b[1]+c->b[2]+c->b[3])/c->b[0]);
+  c->sigma = sigma;
+  c->N = 3;
+
+/*
+  g_printerr ("q %f\n", q);
+  g_printerr ("q2 %f\n", q2);
+  g_printerr ("q3 %f\n", q3);
+  g_printerr ("c->b[0] %f\n", c->b[0]);
+  g_printerr ("c->b[1] %f\n", c->b[1]);
+  g_printerr ("c->b[2] %f\n", c->b[2]);
+  g_printerr ("c->b[3] %f\n", c->b[3]);
+  g_printerr ("c->B %f\n", c->B);
+  g_printerr ("c->sigma %f\n", c->sigma);
+  g_printerr ("c->N %d\n", c->N);
+*/
+}
+
+static void
+gausssmooth (gfloat *in, gfloat *out, gint size, gint rowstride, gauss3_coefs *c)
+{
+  /*
+   * Papers:  "Recursive Implementation of the gaussian filter.",
+   *          Ian T. Young , Lucas J. Van Vliet, Signal Processing 44, Elsevier 1995.
+   * formula: 9a        forward filter
+   *          9b        backward filter
+   *          fig7      algorithm
+   */
+  gint i,n, bufsize;
+  gfloat *w1,*w2;
+
+  /* forward pass */
+  bufsize = size+3;
+  size -= 1;
+  w1 = (gfloat *) g_try_malloc (bufsize * sizeof (gfloat));
+  w2 = (gfloat *) g_try_malloc (bufsize * sizeof (gfloat));
+  w1[0] = in[0];
+  w1[1] = in[0];
+  w1[2] = in[0];
+  for ( i = 0 , n=3; i <= size ; i++, n++)
+    {
+      w1[n] = (gfloat)(c->B*in[i*rowstride] +
+                       ((c->b[1]*w1[n-1] +
+                         c->b[2]*w1[n-2] +
+                         c->b[3]*w1[n-3] ) / c->b[0]));
+    }
+
+  /* backward pass */
+  w2[size+1]= w1[size+3];
+  w2[size+2]= w1[size+3];
+  w2[size+3]= w1[size+3];
+  for (i = size, n = i; i >= 0; i--, n--)
+    {
+      w2[n]= out[i * rowstride] = (gfloat)(c->B*w1[n+3] +
+                                           ((c->b[1]*w2[n+1] +
+                                             c->b[2]*w2[n+2] +
+                                             c->b[3]*w2[n+3] ) / c->b[0]));
+    }
+
+  g_free (w1);
+  g_free (w2);
+}
+
+/*
+ * This function is the heart of the algo.
+ * (a)  Filterings at several scales and sumarize the results.
+ * (b)  Calculation of the final values.
+ */
+static void
+MSRCR (guchar *src, gint width, gint height, gint bytes, gboolean preview_mode)
+{
+
+  gint          scale,row,col;
+  gint          i,j;
+  gint          size;
+  gint          channel;
+  guchar       *psrc = NULL;            /* backup pointer for src buffer */
+  gfloat       *dst  = NULL;            /* float buffer for algorithm */
+  gfloat       *pdst = NULL;            /* backup pointer for float buffer */
+  gfloat       *in, *out;
+  gint          channelsize;            /* Float memory cache for one channel */
+  gfloat        weight;
+  gauss3_coefs  coef;
+  gfloat        mean, var;
+  gfloat        mini, range, maxi;
+  gfloat        alpha;
+  gfloat        gain;
+  gfloat        offset;
+  gdouble       max_preview = 0.0;
+
+  if (!preview_mode)
+    {
+      gimp_progress_init (_("Retinex: filtering"));
+      max_preview = 3 * rvals.nscales;
+    }
+
+  /* Allocate all the memory needed for algorithm*/
+  size = width * height * bytes;
+  dst = g_try_malloc (size * sizeof (gfloat));
+  if (dst == NULL)
+    {
+      g_warning ("Failed to allocate memory");
+      return;
+    }
+  memset (dst, 0, size * sizeof (gfloat));
+
+  channelsize  = (width * height);
+  in  = (gfloat *) g_try_malloc (channelsize * sizeof (gfloat));
+  if (in == NULL)
+    {
+      g_free (dst);
+      g_warning ("Failed to allocate memory");
+      return; /* do some clever stuff */
+    }
+
+  out  = (gfloat *) g_try_malloc (channelsize * sizeof (gfloat));
+  if (out == NULL)
+    {
+      g_free (in);
+      g_free (dst);
+      g_warning ("Failed to allocate memory");
+      return; /* do some clever stuff */
+    }
+
+
+  /*
+     Calculate the scales of filtering according to the
+     number of filter and their distribution.
+   */
+
+  retinex_scales_distribution (RetinexScales,
+                               rvals.nscales, rvals.scales_mode, rvals.scale);
+
+  /*
+      Filtering according to the various scales.
+      Summerize the results of the various filters according to a
+      specific weight(here equivalent for all).
+  */
+  weight = 1./ (gfloat) rvals.nscales;
+
+  /*
+    The recursive filtering algorithm needs different coefficients according
+    to the selected scale (~ = standard deviation of Gaussian).
+   */
+  for (channel = 0; channel < 3; channel++)
+    {
+      gint pos;
+
+      for (i = 0, pos = channel; i < channelsize ; i++, pos += bytes)
+         {
+            /* 0-255 => 1-256 */
+            in[i] = (gfloat)(src[pos] + 1.0);
+         }
+      for (scale = 0; scale < rvals.nscales; scale++)
+        {
+          compute_coefs3 (&coef, RetinexScales[scale]);
+          /*
+           *  Filtering (smoothing) Gaussian recursive.
+           *
+           *  Filter rows first
+           */
+          for (row=0 ;row < height; row++)
+            {
+              pos =  row * width;
+              gausssmooth (in + pos, out + pos, width, 1, &coef);
+            }
+
+          memcpy(in,  out, channelsize * sizeof(gfloat));
+          memset(out, 0  , channelsize * sizeof(gfloat));
+
+          /*
+           *  Filtering (smoothing) Gaussian recursive.
+           *
+           *  Second columns
+           */
+          for (col=0; col < width; col++)
+            {
+              pos = col;
+              gausssmooth(in + pos, out + pos, height, width, &coef);
+            }
+
+          /*
+             Summarize the filtered values.
+             In fact one calculates a ratio between the original values and the filtered values.
+           */
+          for (i = 0, pos = channel; i < channelsize; i++, pos += bytes)
+            {
+              dst[pos] += weight * (log (src[pos] + 1.) - log (out[i]));
+            }
+
+           if (!preview_mode)
+             gimp_progress_update ((channel * rvals.nscales + scale) /
+                                   max_preview);
+        }
+    }
+  g_free(in);
+  g_free(out);
+
+  /*
+      Final calculation with original value and cumulated filter values.
+      The parameters gain, alpha and offset are constants.
+  */
+  /* Ci(x,y)=log[a Ii(x,y)]-log[ Ei=1-s Ii(x,y)] */
+
+  alpha  = 128.;
+  gain   = 1.;
+  offset = 0.;
+
+  for (i = 0; i < size; i += bytes)
+    {
+      gfloat logl;
+
+      psrc = src+i;
+      pdst = dst+i;
+
+      logl = log((gfloat)psrc[0] + (gfloat)psrc[1] + (gfloat)psrc[2] + 3.);
+
+      pdst[0] = gain * ((log(alpha * (psrc[0]+1.)) - logl) * pdst[0]) + offset;
+      pdst[1] = gain * ((log(alpha * (psrc[1]+1.)) - logl) * pdst[1]) + offset;
+      pdst[2] = gain * ((log(alpha * (psrc[2]+1.)) - logl) * pdst[2]) + offset;
+    }
+
+/*  if (!preview_mode)
+    gimp_progress_update ((2.0 + (rvals.nscales * 3)) /
+                          ((rvals.nscales * 3) + 3));*/
+
+  /*
+      Adapt the dynamics of the colors according to the statistics of the first and second order.
+      The use of the variance makes it possible to control the degree of saturation of the colors.
+  */
+  pdst = dst;
+
+  compute_mean_var (pdst, &mean, &var, size, bytes);
+  mini = mean - rvals.cvar*var;
+  maxi = mean + rvals.cvar*var;
+  range = maxi - mini;
+
+  if (!range)
+    range = 1.0;
+
+  for (i = 0; i < size; i+= bytes)
+    {
+      psrc = src + i;
+      pdst = dst + i;
+
+      for (j = 0 ; j < 3 ; j++)
+        {
+          gfloat c = 255 * ( pdst[j] - mini ) / range;
+
+          psrc[j] = (guchar) CLAMP (c, 0, 255);
+        }
+    }
+
+  g_free (dst);
+}
+
+/*
+ * Calculate the average and variance in one go.
+ */
+static void
+compute_mean_var (gfloat *src, gfloat *mean, gfloat *var, gint size, gint bytes)
+{
+  gfloat vsquared;
+  gint i,j;
+  gfloat *psrc;
+
+  vsquared = 0;
+  *mean = 0;
+  for (i = 0; i < size; i+= bytes)
+    {
+       psrc = src+i;
+       for (j = 0 ; j < 3 ; j++)
+         {
+            *mean += psrc[j];
+            vsquared += psrc[j] * psrc[j];
+         }
+    }
+
+  *mean /= (gfloat) size; /* mean */
+  vsquared /= (gfloat) size; /* mean (x^2) */
+  *var = ( vsquared - (*mean * *mean) );
+  *var = sqrt(*var); /* var */
+}