1
/* GIMP - The GNU Image Manipulation Program
2
* Copyright (C) 1995 Spencer Kimball and Peter Mattis
4
* This program is free software; you can redistribute it and/or modify
5
* it under the terms of the GNU General Public License as published by
6
* the Free Software Foundation; either version 2 of the License, or
7
* (at your option) any later version.
9
* This program is distributed in the hope that it will be useful,
10
* but WITHOUT ANY WARRANTY; without even the implied warranty of
11
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12
* GNU General Public License for more details.
14
* You should have received a copy of the GNU General Public License
15
* along with this program; if not, write to the Free Software
16
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
23
#include "libgimp/gimp.h"
24
#include "libgimp/gimpui.h"
26
#include "libgimp/stdplugins-intl.h"
29
#define PLUG_IN_PROC "plug-in-retinex"
30
#define PLUG_IN_BINARY "contrast-retinex"
31
#define MAX_RETINEX_SCALES 8
32
#define MIN_GAUSSIAN_SCALE 16
33
#define MAX_GAUSSIAN_SCALE 250
34
#define SCALE_WIDTH 150
54
Definit comment sont repartis les
55
differents filtres en fonction de
56
l'echelle (~= ecart type de la gaussienne)
58
#define RETINEX_UNIFORM 0
60
#define RETINEX_HIGH 2
62
static gfloat RetinexScales[MAX_RETINEX_SCALES];
74
* Declare local functions.
76
static void query (void);
77
static void run (const gchar *name,
79
const GimpParam *param,
81
GimpParam **return_vals);
84
static gboolean retinex_dialog (GimpDrawable *drawable);
85
static void retinex (GimpDrawable *drawable,
86
GimpPreview *preview);
88
static void retinex_scales_distribution (gfloat *scales,
93
static void compute_mean_var (gfloat *src,
101
static void compute_coefs3 (gauss3_coefs *c,
104
static void gausssmooth (gfloat *in,
111
* MSRCR = MultiScale Retinex with Color Restoration
113
static void MSRCR (guchar *src,
117
gboolean preview_mode);
123
static RetinexParams rvals =
127
RETINEX_UNIFORM, /* Echelles reparties uniformement */
131
static GimpPlugInInfo PLUG_IN_INFO =
133
NULL, /* init_proc */
134
NULL, /* quit_proc */
135
query, /* query_proc */
144
static const GimpParamDef args[] =
146
{ GIMP_PDB_INT32, "run-mode", "Interactive, non-interactive" },
147
{ GIMP_PDB_IMAGE, "image", "Input image (unused)" },
148
{ GIMP_PDB_DRAWABLE, "drawable", "Input drawable" },
149
{ GIMP_PDB_INT32, "scale", "Biggest scale value" },
150
{ GIMP_PDB_INT32, "nscales", "Number of scales" },
151
{ GIMP_PDB_INT32, "scales-mode", "Retinex distribution through scales" },
152
{ GIMP_PDB_FLOAT, "cvar", "Variance value" }
155
gimp_install_procedure (PLUG_IN_PROC,
156
N_("Enhance contrast using the Retinex method"),
157
"The Retinex Image Enhancement Algorithm is an "
158
"automatic image enhancement method that enhances "
159
"a digital image in terms of dynamic range "
160
"compression, color independence from the spectral "
161
"distribution of the scene illuminant, and "
162
"color/lightness rendition.",
169
G_N_ELEMENTS (args), 0,
172
gimp_plugin_menu_register (PLUG_IN_PROC, "<Image>/Colors/Modify");
176
run (const gchar *name,
178
const GimpParam *param,
180
GimpParam **return_vals)
182
static GimpParam values[1];
183
GimpDrawable *drawable;
184
GimpRunMode run_mode;
185
GimpPDBStatusType status = GIMP_PDB_SUCCESS;
186
gint x, y, width, height;
188
run_mode = param[0].data.d_int32;
193
*return_vals = values;
195
values[0].type = GIMP_PDB_STATUS;
196
values[0].data.d_status = status;
198
drawable = gimp_drawable_get (param[2].data.d_drawable);
200
if (! gimp_drawable_mask_intersect (drawable->drawable_id,
201
&x, &y, &width, &height) ||
202
width < MIN_GAUSSIAN_SCALE ||
203
height < MIN_GAUSSIAN_SCALE)
205
status = GIMP_PDB_EXECUTION_ERROR;
206
gimp_drawable_detach (drawable);
207
values[0].data.d_status = status;
211
gimp_tile_cache_ntiles (2 * (drawable->width / gimp_tile_width () + 1));
215
case GIMP_RUN_INTERACTIVE:
216
/* Possibly retrieve data */
217
gimp_get_data (PLUG_IN_PROC, &rvals);
219
/* First acquire information with a dialog */
220
if (! retinex_dialog (drawable))
224
case GIMP_RUN_NONINTERACTIVE:
225
/* Make sure all the arguments are there! */
228
status = GIMP_PDB_CALLING_ERROR;
232
rvals.scale = (param[3].data.d_int32);
233
rvals.nscales = (param[4].data.d_int32);
234
rvals.scales_mode = (param[5].data.d_int32);
235
rvals.cvar = (param[6].data.d_float);
239
case GIMP_RUN_WITH_LAST_VALS:
240
gimp_get_data (PLUG_IN_PROC, &rvals);
247
if ((status == GIMP_PDB_SUCCESS) &&
248
(gimp_drawable_is_rgb (drawable->drawable_id)))
250
gimp_progress_init (_("Retinex"));
252
retinex (drawable, NULL);
254
if (run_mode != GIMP_RUN_NONINTERACTIVE)
255
gimp_displays_flush ();
258
if (run_mode == GIMP_RUN_INTERACTIVE)
259
gimp_set_data (PLUG_IN_PROC, &rvals, sizeof (RetinexParams));
263
status = GIMP_PDB_EXECUTION_ERROR;
266
gimp_drawable_detach (drawable);
268
values[0].data.d_status = status;
273
retinex_dialog (GimpDrawable *drawable)
276
GtkWidget *main_vbox;
284
gimp_ui_init (PLUG_IN_BINARY, FALSE);
286
dialog = gimp_dialog_new (_("Retinex Image Enhancement"), PLUG_IN_BINARY,
288
gimp_standard_help_func, PLUG_IN_PROC,
290
GTK_STOCK_CANCEL, GTK_RESPONSE_CANCEL,
291
GTK_STOCK_OK, GTK_RESPONSE_OK,
295
gtk_dialog_set_alternative_button_order (GTK_DIALOG (dialog),
300
gimp_window_set_transient (GTK_WINDOW (dialog));
302
main_vbox = gtk_vbox_new (FALSE, 12);
303
gtk_container_set_border_width (GTK_CONTAINER (main_vbox), 12);
304
gtk_container_add (GTK_CONTAINER (GTK_DIALOG (dialog)->vbox), main_vbox);
305
gtk_widget_show (main_vbox);
307
preview = gimp_zoom_preview_new (drawable);
308
gtk_box_pack_start (GTK_BOX (main_vbox), preview, TRUE, TRUE, 0);
309
gtk_widget_show (preview);
311
g_signal_connect_swapped (preview, "invalidated",
312
G_CALLBACK (retinex),
315
table = gtk_table_new (4, 3, FALSE);
316
gtk_table_set_col_spacings (GTK_TABLE (table), 6);
317
gtk_table_set_row_spacings (GTK_TABLE (table), 6);
318
gtk_box_pack_start (GTK_BOX (main_vbox), table, FALSE, FALSE, 0);
319
gtk_widget_show (table);
321
combo = gimp_int_combo_box_new (_("Uniform"), filter_uniform,
322
_("Low"), filter_low,
323
_("High"), filter_high,
326
gimp_int_combo_box_connect (GIMP_INT_COMBO_BOX (combo), rvals.scales_mode,
327
G_CALLBACK (gimp_int_combo_box_get_active),
329
g_signal_connect_swapped (combo, "changed",
330
G_CALLBACK (gimp_preview_invalidate),
333
gimp_table_attach_aligned (GTK_TABLE (table), 0, 0,
334
_("_Level:"), 0.0, 0.5,
336
gtk_widget_show (combo);
338
adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 1,
339
_("_Scale:"), SCALE_WIDTH, ENTRY_WIDTH,
341
MIN_GAUSSIAN_SCALE, MAX_GAUSSIAN_SCALE, 1, 1, 0,
342
TRUE, 0, 0, NULL, NULL);
343
scale = GIMP_SCALE_ENTRY_SCALE (adj);
344
gtk_range_set_update_policy (GTK_RANGE (scale), GTK_UPDATE_DISCONTINUOUS);
346
g_signal_connect (adj, "value-changed",
347
G_CALLBACK (gimp_int_adjustment_update),
349
g_signal_connect_swapped (adj, "value-changed",
350
G_CALLBACK (gimp_preview_invalidate),
353
adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 2,
354
_("Scale _division:"), SCALE_WIDTH, ENTRY_WIDTH,
356
0, MAX_RETINEX_SCALES, 1, 1, 0,
357
TRUE, 0, 0, NULL, NULL);
358
scale = GIMP_SCALE_ENTRY_SCALE (adj);
359
gtk_range_set_update_policy (GTK_RANGE (scale), GTK_UPDATE_DISCONTINUOUS);
361
g_signal_connect (adj, "value-changed",
362
G_CALLBACK (gimp_int_adjustment_update),
364
g_signal_connect_swapped (adj, "value-changed",
365
G_CALLBACK (gimp_preview_invalidate),
368
adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 3,
369
_("Dy_namic:"), SCALE_WIDTH, ENTRY_WIDTH,
370
rvals.cvar, 0, 4, 0.1, 0.1, 1,
371
TRUE, 0, 0, NULL, NULL);
372
scale = GIMP_SCALE_ENTRY_SCALE (adj);
373
gtk_range_set_update_policy (GTK_RANGE (scale), GTK_UPDATE_DISCONTINUOUS);
375
g_signal_connect (adj, "value-changed",
376
G_CALLBACK (gimp_float_adjustment_update),
378
g_signal_connect_swapped (adj, "value-changed",
379
G_CALLBACK (gimp_preview_invalidate),
382
gtk_widget_show (dialog);
384
run = (gimp_dialog_run (GIMP_DIALOG (dialog)) == GTK_RESPONSE_OK);
386
gtk_widget_destroy (dialog);
392
* Applies the algorithm
395
retinex (GimpDrawable *drawable,
396
GimpPreview *preview)
398
gint x, y, width, height;
402
GimpPixelRgn dst_rgn, src_rgn;
404
bytes = drawable->bpp;
407
* Get the size of the current image or its selection.
411
src = gimp_zoom_preview_get_source (GIMP_ZOOM_PREVIEW (preview),
412
&width, &height, &bytes);
416
if (! gimp_drawable_mask_intersect (drawable->drawable_id,
417
&x, &y, &width, &height))
420
/* Allocate memory */
421
size = width * height * bytes;
422
src = g_try_malloc (sizeof (guchar) * size);
426
g_warning ("Failed to allocate memory");
430
memset (src, 0, sizeof (guchar) * size);
432
/* Fill allocated memory with pixel data */
433
gimp_pixel_rgn_init (&src_rgn, drawable,
436
gimp_pixel_rgn_get_rect (&src_rgn, src, x, y, width, height);
440
Algorithm for Multi-scale Retinex with color Restoration (MSRCR).
443
MSRCR (psrc, width, height, bytes, preview != NULL);
447
gimp_preview_draw_buffer (preview, psrc, width * bytes);
451
gimp_pixel_rgn_init (&dst_rgn, drawable,
454
gimp_pixel_rgn_set_rect (&dst_rgn, psrc, x, y, width, height);
456
gimp_progress_update (1.0);
458
gimp_drawable_flush (drawable);
459
gimp_drawable_merge_shadow (drawable->drawable_id, TRUE);
460
gimp_drawable_update (drawable->drawable_id, x, y, width, height);
468
* calculate scale values for desired distribution.
471
retinex_scales_distribution(gfloat* scales, gint nscales, gint mode, gint s)
474
{ /* For one filter we choose the median scale */
475
scales[0] = (gint) s / 2;
477
else if (nscales == 2)
478
{ /* For two filters whe choose the median and maximum scale */
479
scales[0] = (gint) s / 2;
480
scales[1] = (gint) s;
484
gfloat size_step = (gfloat) s / (gfloat) nscales;
489
case RETINEX_UNIFORM:
490
for(i = 0; i < nscales; ++i)
491
scales[i] = 2. + (gfloat) i * size_step;
495
size_step = (gfloat) log(s - 2.0) / (gfloat) nscales;
496
for (i = 0; i < nscales; ++i)
497
scales[i] = 2. + pow (10, (i * size_step) / log (10));
501
size_step = (gfloat) log(s - 2.0) / (gfloat) nscales;
502
for (i = 0; i < nscales; ++i)
503
scales[i] = s - pow (10, (i * size_step) / log (10));
513
* Calculate the coefficients for the recursive filter algorithm
514
* Fast Computation of gaussian blurring.
517
compute_coefs3 (gauss3_coefs *c, gfloat sigma)
520
* Papers: "Recursive Implementation of the gaussian filter.",
521
* Ian T. Young , Lucas J. Van Vliet, Signal Processing 44, Elsevier 1995.
522
* formula: 11b computation of q
523
* 8c computation of b0..b1
524
* 10 alpha is normalization constant B
532
q = 0.98711 * sigma - 0.96330;
534
else if ((sigma >= 0.5) && (sigma < 2.5))
536
q = 3.97156 - 4.14554 * (gfloat) sqrt ((double) 1 - 0.26891 * sigma);
540
q = 0.1147705018520355224609375;
545
c->b[0] = (1.57825+(2.44413*q)+(1.4281 *q2)+(0.422205*q3));
546
c->b[1] = ( (2.44413*q)+(2.85619*q2)+(1.26661 *q3));
547
c->b[2] = ( -((1.4281*q2)+(1.26661 *q3)));
548
c->b[3] = ( (0.422205*q3));
549
c->B = 1.0-((c->b[1]+c->b[2]+c->b[3])/c->b[0]);
554
g_printerr ("q %f\n", q);
555
g_printerr ("q2 %f\n", q2);
556
g_printerr ("q3 %f\n", q3);
557
g_printerr ("c->b[0] %f\n", c->b[0]);
558
g_printerr ("c->b[1] %f\n", c->b[1]);
559
g_printerr ("c->b[2] %f\n", c->b[2]);
560
g_printerr ("c->b[3] %f\n", c->b[3]);
561
g_printerr ("c->B %f\n", c->B);
562
g_printerr ("c->sigma %f\n", c->sigma);
563
g_printerr ("c->N %d\n", c->N);
568
gausssmooth (gfloat *in, gfloat *out, gint size, gint rowstride, gauss3_coefs *c)
571
* Papers: "Recursive Implementation of the gaussian filter.",
572
* Ian T. Young , Lucas J. Van Vliet, Signal Processing 44, Elsevier 1995.
573
* formula: 9a forward filter
583
w1 = (gfloat *) g_try_malloc (bufsize * sizeof (gfloat));
584
w2 = (gfloat *) g_try_malloc (bufsize * sizeof (gfloat));
588
for ( i = 0 , n=3; i <= size ; i++, n++)
590
w1[n] = (gfloat)(c->B*in[i*rowstride] +
593
c->b[3]*w1[n-3] ) / c->b[0]));
597
w2[size+1]= w1[size+3];
598
w2[size+2]= w1[size+3];
599
w2[size+3]= w1[size+3];
600
for (i = size, n = i; i >= 0; i--, n--)
602
w2[n]= out[i * rowstride] = (gfloat)(c->B*w1[n] +
605
c->b[3]*w2[n+3] ) / c->b[0]));
613
* This function is the heart of the algo.
614
* (a) Filterings at several scales and sumarize the results.
615
* (b) Calculation of the final values.
618
MSRCR (guchar *src, gint width, gint height, gint bytes, gboolean preview_mode)
626
guchar *psrc = NULL; /* backup pointer for src buffer */
627
gfloat *dst = NULL; /* float buffer for algorithm */
628
gfloat *pdst = NULL; /* backup pointer for float buffer */
630
gint channelsize; /* Float memory cache for one channel */
634
gfloat mini, range, maxi;
638
gdouble max_preview = 0.0;
642
gimp_progress_init (_("Retinex: filtering"));
643
max_preview = 3 * rvals.nscales;
646
/* Allocate all the memory needed for algorithm*/
647
size = width * height * bytes;
648
dst = g_try_malloc (size * sizeof (gfloat));
651
g_warning ("Failed to allocate memory");
654
memset (dst, 0, size * sizeof (gfloat));
656
channelsize = (width * height);
657
in = (gfloat *) g_try_malloc (channelsize * sizeof (gfloat));
661
g_warning ("Failed to allocate memory");
662
return; /* do some clever stuff */
665
out = (gfloat *) g_try_malloc (channelsize * sizeof (gfloat));
670
g_warning ("Failed to allocate memory");
671
return; /* do some clever stuff */
676
Calculate the scales of filtering according to the
677
number of filter and their distribution.
680
retinex_scales_distribution (RetinexScales,
681
rvals.nscales, rvals.scales_mode, rvals.scale);
684
Filtering according to the various scales.
685
Summerize the results of the various filters according to a
686
specific weight(here equivalent for all).
688
weight = 1./ (gfloat) rvals.nscales;
691
The recursive filtering algorithm needs different coefficients according
692
to the selected scale (~ = standard deviation of Gaussian).
695
for (channel = 0; channel < 3; channel++)
697
for (i = 0, pos = channel; i < channelsize ; i++, pos += bytes)
700
in[i] = (gfloat)(src[pos] + 1.0);
702
for (scale = 0; scale < rvals.nscales; scale++)
704
compute_coefs3 (&coef, RetinexScales[scale]);
706
* Filtering (smoothing) Gaussian recursive.
710
for (row=0 ;row < height; row++)
713
gausssmooth (in + pos, out + pos, width, 1, &coef);
716
memcpy(in, out, channelsize * sizeof(gfloat));
717
memset(out, 0 , channelsize * sizeof(gfloat));
720
* Filtering (smoothing) Gaussian recursive.
724
for (col=0; col < width; col++)
727
gausssmooth(in + pos, out + pos, height, width, &coef);
731
Summarize the filtered values.
732
In fact one calculates a ratio between the original values and the filtered values.
734
for (i = 0, pos = channel; i < channelsize; i++, pos += bytes)
736
dst[pos] += weight * (log (src[pos] + 1.) - log (out[i]));
740
gimp_progress_update ((channel * rvals.nscales + scale) /
748
Final calculation with original value and cumulated filter values.
749
The parameters gain, alpha and offset are constants.
751
/* Ci(x,y)=log[a Ii(x,y)]-log[ Ei=1-s Ii(x,y)] */
757
for (i = 0; i < size; i += bytes)
764
logl = log((gfloat)psrc[0] + (gfloat)psrc[1] + (gfloat)psrc[2] + 3.);
766
pdst[0] = gain * ((log(alpha * (psrc[0]+1.)) - logl) * pdst[0]) + offset;
767
pdst[1] = gain * ((log(alpha * (psrc[1]+1.)) - logl) * pdst[1]) + offset;
768
pdst[2] = gain * ((log(alpha * (psrc[2]+1.)) - logl) * pdst[2]) + offset;
771
/* if (!preview_mode)
772
gimp_progress_update ((2.0 + (rvals.nscales * 3)) /
773
((rvals.nscales * 3) + 3));*/
776
Adapt the dynamics of the colors according to the statistics of the first and second order.
777
The use of the variance makes it possible to control the degree of saturation of the colors.
781
compute_mean_var (pdst, &mean, &var, size, bytes);
782
mini = mean - rvals.cvar*var;
783
maxi = mean + rvals.cvar*var;
789
for (i = 0; i < size; i+= bytes)
794
for (j = 0 ; j < 3 ; j++)
796
gfloat c = 255 * ( pdst[j] - mini ) / range;
798
psrc[j] = (guchar) CLAMP (c, 0, 255);
806
* Calculate the average and variance in one go.
809
compute_mean_var (gfloat *src, gfloat *mean, gfloat *var, gint size, gint bytes)
817
for (i = 0; i < size; i+= bytes)
820
for (j = 0 ; j < 3 ; j++)
823
vsquared += psrc[j] * psrc[j];
827
*mean /= (gfloat) size; /* mean */
828
vsquared /= (gfloat) size; /* mean (x^2) */
829
*var = ( vsquared - (*mean * *mean) );
830
*var = sqrt(*var); /* var */
added
removed
plug-ins/common/contrast-retinex.c
