summaryrefslogtreecommitdiffstats
path: root/plug-ins/common/contrast-retinex.c
diff options
context:
space:
mode:
authorDaniel Baumann <daniel.baumann@progress-linux.org>2024-04-07 18:30:19 +0000
committerDaniel Baumann <daniel.baumann@progress-linux.org>2024-04-07 18:30:19 +0000
commit5c1676dfe6d2f3c837a5e074117b45613fd29a72 (patch)
treecbffb45144febf451e54061db2b21395faf94bfe /plug-ins/common/contrast-retinex.c
parentInitial commit. (diff)
downloadgimp-5c1676dfe6d2f3c837a5e074117b45613fd29a72.tar.xz
gimp-5c1676dfe6d2f3c837a5e074117b45613fd29a72.zip
Adding upstream version 2.10.34.upstream/2.10.34upstream
Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>
Diffstat (limited to 'plug-ins/common/contrast-retinex.c')
-rw-r--r--plug-ins/common/contrast-retinex.c840
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 */
+}