853 lines
26 KiB
C
853 lines
26 KiB
C
/* PIKA - Photo and Image Kooker Application
|
|
* a rebranding of The GNU Image Manipulation Program (created with heckimp)
|
|
* A derived work which may be trivial. However, any changes may be (C)2023 by Aldercone Studio
|
|
*
|
|
* Original copyright, applying to most contents (license remains unchanged):
|
|
* 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 "libpika/pika.h"
|
|
#include "libpika/pikaui.h"
|
|
|
|
#include "libpika/stdplugins-intl.h"
|
|
|
|
|
|
#define PLUG_IN_PROC "plug-in-retinex"
|
|
#define PLUG_IN_BINARY "contrast-retinex"
|
|
#define PLUG_IN_ROLE "pika-contrast-retinex"
|
|
#define MAX_RETINEX_SCALES 8
|
|
#define MIN_GAUSSIAN_SCALE 16
|
|
#define MAX_GAUSSIAN_SCALE 250
|
|
|
|
|
|
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;
|
|
|
|
|
|
typedef struct _Retinex Retinex;
|
|
typedef struct _RetinexClass RetinexClass;
|
|
|
|
struct _Retinex
|
|
{
|
|
PikaPlugIn parent_instance;
|
|
};
|
|
|
|
struct _RetinexClass
|
|
{
|
|
PikaPlugInClass parent_class;
|
|
};
|
|
|
|
|
|
#define RETINEX_TYPE (retinex_get_type ())
|
|
#define RETINEX(obj) (G_TYPE_CHECK_INSTANCE_CAST ((obj), RETINEX_TYPE, Retinex))
|
|
|
|
GType retinex_get_type (void) G_GNUC_CONST;
|
|
|
|
static GList * retinex_query_procedures (PikaPlugIn *plug_in);
|
|
static PikaProcedure * retinex_create_procedure (PikaPlugIn *plug_in,
|
|
const gchar *name);
|
|
|
|
static PikaValueArray * retinex_run (PikaProcedure *procedure,
|
|
PikaRunMode run_mode,
|
|
PikaImage *image,
|
|
gint n_drawables,
|
|
PikaDrawable **drawables,
|
|
PikaProcedureConfig *config,
|
|
gpointer run_data);
|
|
|
|
static gboolean retinex_dialog (PikaProcedure *procedure,
|
|
GObject *config,
|
|
PikaDrawable *drawable);
|
|
static void retinex (GObject *config,
|
|
PikaDrawable *drawable,
|
|
PikaPreview *preview);
|
|
static void retinex_preview (GtkWidget *widget,
|
|
GObject *config);
|
|
|
|
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);
|
|
|
|
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 (GObject *config,
|
|
guchar *src,
|
|
gint width,
|
|
gint height,
|
|
gint bytes,
|
|
gboolean preview_mode);
|
|
|
|
|
|
G_DEFINE_TYPE (Retinex, retinex, PIKA_TYPE_PLUG_IN)
|
|
|
|
PIKA_MAIN (RETINEX_TYPE)
|
|
DEFINE_STD_SET_I18N
|
|
|
|
|
|
static void
|
|
retinex_class_init (RetinexClass *klass)
|
|
{
|
|
PikaPlugInClass *plug_in_class = PIKA_PLUG_IN_CLASS (klass);
|
|
|
|
plug_in_class->query_procedures = retinex_query_procedures;
|
|
plug_in_class->create_procedure = retinex_create_procedure;
|
|
plug_in_class->set_i18n = STD_SET_I18N;
|
|
}
|
|
|
|
static void
|
|
retinex_init (Retinex *retinex)
|
|
{
|
|
}
|
|
|
|
static GList *
|
|
retinex_query_procedures (PikaPlugIn *plug_in)
|
|
{
|
|
return g_list_append (NULL, g_strdup (PLUG_IN_PROC));
|
|
}
|
|
|
|
static PikaProcedure *
|
|
retinex_create_procedure (PikaPlugIn *plug_in,
|
|
const gchar *name)
|
|
{
|
|
PikaProcedure *procedure = NULL;
|
|
|
|
if (! strcmp (name, PLUG_IN_PROC))
|
|
{
|
|
procedure = pika_image_procedure_new (plug_in, name,
|
|
PIKA_PDB_PROC_TYPE_PLUGIN,
|
|
retinex_run, NULL, NULL);
|
|
|
|
pika_procedure_set_image_types (procedure, "RGB*");
|
|
pika_procedure_set_sensitivity_mask (procedure,
|
|
PIKA_PROCEDURE_SENSITIVE_DRAWABLE);
|
|
|
|
pika_procedure_set_menu_label (procedure, _("Retine_x..."));
|
|
pika_procedure_add_menu_path (procedure, "<Image>/Colors/Tone Mapping");
|
|
|
|
pika_procedure_set_documentation (procedure,
|
|
_("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."),
|
|
name);
|
|
pika_procedure_set_attribution (procedure,
|
|
"Fabien Pelisson",
|
|
"Fabien Pelisson",
|
|
"2003");
|
|
|
|
PIKA_PROC_ARG_INT (procedure, "scale",
|
|
_("Scal_e"),
|
|
_("Biggest scale value"),
|
|
MIN_GAUSSIAN_SCALE, MAX_GAUSSIAN_SCALE, 240,
|
|
G_PARAM_READWRITE);
|
|
|
|
PIKA_PROC_ARG_INT (procedure, "nscales",
|
|
_("Scale _division"),
|
|
_("Number of scales"),
|
|
0, MAX_RETINEX_SCALES, 3,
|
|
G_PARAM_READWRITE);
|
|
|
|
PIKA_PROC_ARG_INT (procedure, "scales-mode",
|
|
_("Le_vel"),
|
|
_("Retinex distribution through scales "
|
|
"{ Uniform (0), Low (1), High (2) }"),
|
|
RETINEX_UNIFORM, RETINEX_HIGH, RETINEX_UNIFORM,
|
|
G_PARAM_READWRITE);
|
|
|
|
PIKA_PROC_ARG_DOUBLE (procedure, "cvar",
|
|
_("Dy_namic"),
|
|
_("Variance value"),
|
|
0, 4, 1.2,
|
|
G_PARAM_READWRITE);
|
|
}
|
|
|
|
return procedure;
|
|
}
|
|
|
|
static PikaValueArray *
|
|
retinex_run (PikaProcedure *procedure,
|
|
PikaRunMode run_mode,
|
|
PikaImage *image,
|
|
gint n_drawables,
|
|
PikaDrawable **drawables,
|
|
PikaProcedureConfig *config,
|
|
gpointer run_data)
|
|
{
|
|
PikaDrawable *drawable;
|
|
gint x, y, width, height;
|
|
|
|
gegl_init (NULL, NULL);
|
|
|
|
if (n_drawables != 1)
|
|
{
|
|
GError *error = NULL;
|
|
|
|
g_set_error (&error, PIKA_PLUG_IN_ERROR, 0,
|
|
_("Procedure '%s' only works with one drawable."),
|
|
PLUG_IN_PROC);
|
|
|
|
return pika_procedure_new_return_values (procedure,
|
|
PIKA_PDB_CALLING_ERROR,
|
|
error);
|
|
}
|
|
else
|
|
{
|
|
drawable = drawables[0];
|
|
}
|
|
|
|
if (! pika_drawable_mask_intersect (drawable, &x, &y, &width, &height) ||
|
|
width < MIN_GAUSSIAN_SCALE ||
|
|
height < MIN_GAUSSIAN_SCALE)
|
|
return pika_procedure_new_return_values (procedure,
|
|
PIKA_PDB_EXECUTION_ERROR,
|
|
NULL);
|
|
|
|
if (run_mode == PIKA_RUN_INTERACTIVE && ! retinex_dialog (procedure, G_OBJECT (config), drawable))
|
|
return pika_procedure_new_return_values (procedure,
|
|
PIKA_PDB_CANCEL,
|
|
NULL);
|
|
|
|
if (pika_drawable_is_rgb (drawable))
|
|
{
|
|
pika_progress_init (_("Retinex"));
|
|
|
|
retinex (G_OBJECT (config), drawable, NULL);
|
|
|
|
if (run_mode != PIKA_RUN_NONINTERACTIVE)
|
|
pika_displays_flush ();
|
|
}
|
|
else
|
|
{
|
|
return pika_procedure_new_return_values (procedure,
|
|
PIKA_PDB_EXECUTION_ERROR,
|
|
NULL);
|
|
}
|
|
|
|
return pika_procedure_new_return_values (procedure, PIKA_PDB_SUCCESS, NULL);
|
|
}
|
|
|
|
|
|
static gboolean
|
|
retinex_dialog (PikaProcedure *procedure,
|
|
GObject *config,
|
|
PikaDrawable *drawable)
|
|
{
|
|
GtkWidget *dialog;
|
|
GtkWidget *preview;
|
|
GtkWidget *combo;
|
|
GtkWidget *scale;
|
|
GtkListStore *store;
|
|
gboolean run;
|
|
|
|
pika_ui_init (PLUG_IN_BINARY);
|
|
|
|
dialog = pika_procedure_dialog_new (procedure,
|
|
PIKA_PROCEDURE_CONFIG (config),
|
|
_("Retinex Image Enhancement"));
|
|
|
|
pika_dialog_set_alternative_button_order (GTK_DIALOG (dialog),
|
|
GTK_RESPONSE_OK,
|
|
GTK_RESPONSE_CANCEL,
|
|
-1);
|
|
|
|
pika_window_set_transient (GTK_WINDOW (dialog));
|
|
|
|
preview = pika_zoom_preview_new_from_drawable (drawable);
|
|
gtk_box_pack_start (GTK_BOX (gtk_dialog_get_content_area (GTK_DIALOG (dialog))),
|
|
preview, TRUE, TRUE, 0);
|
|
gtk_widget_set_margin_bottom (preview, 12);
|
|
gtk_widget_show (preview);
|
|
|
|
store = pika_int_store_new (_("Uniform"), filter_uniform,
|
|
_("Low"), filter_low,
|
|
_("High"), filter_high,
|
|
NULL);
|
|
combo = pika_procedure_dialog_get_int_combo (PIKA_PROCEDURE_DIALOG (dialog),
|
|
"scales-mode",
|
|
PIKA_INT_STORE (store));
|
|
gtk_widget_set_margin_bottom (combo, 12);
|
|
|
|
scale = pika_procedure_dialog_get_scale_entry (PIKA_PROCEDURE_DIALOG (dialog),
|
|
"scale", 1.0);
|
|
gtk_widget_set_margin_bottom (scale, 12);
|
|
|
|
scale = pika_procedure_dialog_get_scale_entry (PIKA_PROCEDURE_DIALOG (dialog),
|
|
"nscales", 1.0);
|
|
gtk_widget_set_margin_bottom (scale, 12);
|
|
|
|
scale = pika_procedure_dialog_get_scale_entry (PIKA_PROCEDURE_DIALOG (dialog),
|
|
"cvar", 1.0);
|
|
gtk_widget_set_margin_bottom (scale, 12);
|
|
|
|
pika_procedure_dialog_fill_box (PIKA_PROCEDURE_DIALOG (dialog),
|
|
"retinex-vbox", "scales-mode",
|
|
"scale", "nscales", "cvar", NULL);
|
|
|
|
g_object_set_data (config, "drawable", drawable);
|
|
|
|
g_signal_connect (preview, "invalidated",
|
|
G_CALLBACK (retinex_preview),
|
|
config);
|
|
|
|
g_signal_connect_swapped (config, "notify",
|
|
G_CALLBACK (pika_preview_invalidate),
|
|
preview);
|
|
|
|
pika_procedure_dialog_fill (PIKA_PROCEDURE_DIALOG (dialog),
|
|
"retinex-vbox",
|
|
NULL);
|
|
|
|
gtk_widget_show (dialog);
|
|
|
|
run = pika_procedure_dialog_run (PIKA_PROCEDURE_DIALOG (dialog));
|
|
|
|
gtk_widget_destroy (dialog);
|
|
|
|
return run;
|
|
}
|
|
|
|
/*
|
|
* Applies the algorithm
|
|
*/
|
|
static void
|
|
retinex (GObject *config,
|
|
PikaDrawable *drawable,
|
|
PikaPreview *preview)
|
|
{
|
|
GeglBuffer *src_buffer = NULL;
|
|
GeglBuffer *dest_buffer;
|
|
const Babl *format = NULL;
|
|
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 = pika_zoom_preview_get_source (PIKA_ZOOM_PREVIEW (preview),
|
|
&width, &height, &bytes);
|
|
}
|
|
else
|
|
{
|
|
if (! pika_drawable_mask_intersect (drawable,
|
|
&x, &y, &width, &height))
|
|
return;
|
|
|
|
if (pika_drawable_has_alpha (drawable))
|
|
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 = pika_drawable_get_buffer (drawable);
|
|
|
|
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 (config, psrc, width, height, bytes, preview != NULL);
|
|
|
|
if (preview)
|
|
{
|
|
pika_preview_draw_buffer (preview, psrc, width * bytes);
|
|
}
|
|
else
|
|
{
|
|
dest_buffer = pika_drawable_get_shadow_buffer (drawable);
|
|
|
|
gegl_buffer_set (dest_buffer, GEGL_RECTANGLE (x, y, width, height), 0,
|
|
format, psrc,
|
|
GEGL_AUTO_ROWSTRIDE);
|
|
|
|
if (src_buffer)
|
|
g_object_unref (src_buffer);
|
|
g_object_unref (dest_buffer);
|
|
|
|
pika_progress_update (1.0);
|
|
|
|
pika_drawable_merge_shadow (drawable, TRUE);
|
|
pika_drawable_update (drawable, x, y, width, height);
|
|
}
|
|
|
|
g_free (src);
|
|
}
|
|
|
|
static void
|
|
retinex_preview (GtkWidget *widget,
|
|
GObject *config)
|
|
{
|
|
PikaPreview *preview = PIKA_PREVIEW (widget);
|
|
PikaDrawable *drawable = g_object_get_data (config, "drawable");
|
|
|
|
retinex (config, drawable, 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 (GObject *config,
|
|
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;
|
|
gint scales_mode;
|
|
gint config_scale;
|
|
gint nscales;
|
|
gdouble cvar;
|
|
|
|
g_object_get (config,
|
|
"scales-mode", &scales_mode,
|
|
"scale", &config_scale,
|
|
"nscales", &nscales,
|
|
"cvar", &cvar,
|
|
NULL);
|
|
|
|
if (!preview_mode)
|
|
{
|
|
pika_progress_init (_("Retinex: filtering"));
|
|
max_preview = 3 * 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,
|
|
nscales, scales_mode, config_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) 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 < 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)
|
|
pika_progress_update ((channel * 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)
|
|
pika_progress_update ((2.0 + (nscales * 3)) /
|
|
((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 - cvar * var;
|
|
maxi = mean + 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 */
|
|
}
|
