1 /**************************************************
2 * file: nlfilt/nlfilt.c
3 *
4 * Copyright (c) 1997 Eric L. Hernes (erich@rrnet.com)
5 * All rights reserved.
6 *
7 * Redistribution and use in source and binary forms, with or without
8 * modification, are permitted provided that the following conditions
9 * are met:
10 * 1. Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * 2. The name of the author may not be used to endorse or promote products
13 * derived from this software without specific prior written permission
14 *
15 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
16 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
17 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
18 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
19 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
20 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
24 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 */
26
27 /*
28 * Algorithm fixes, V2.0 compatibility by David Hodson hodsond@ozemail.com.au
29 */
30
31 #include "config.h"
32
33 #include <string.h>
34
35 #include <libgimp/gimp.h>
36 #include <libgimp/gimpui.h>
37
38 #include "libgimp/stdplugins-intl.h"
39
40
41 #define PLUG_IN_PROC "plug-in-nlfilt"
42 #define PLUG_IN_BINARY "nl-filter"
43 #define PLUG_IN_ROLE "gimp-nl-filter"
44
45
46 typedef struct
47 {
48 gdouble alpha;
49 gdouble radius;
50 gint filter;
51 } NLFilterValues;
52
53 typedef enum
54 {
55 filter_alpha_trim,
56 filter_opt_est,
57 filter_edge_enhance
58 } FilterType;
59
60 static NLFilterValues nlfvals =
61 {
62 0.3,
63 0.3,
64 0
65 };
66
67 /* function protos */
68
69 static void query (void);
70 static void run (const gchar *name,
71 gint nparam,
72 const GimpParam *param,
73 gint *nretvals,
74 GimpParam **retvals);
75
76 static void nlfilter (gint32 drawable_id,
77 GimpPreview *preview);
78 static void nlfilter_preview (gpointer drawable_id,
79 GimpPreview *preview);
80
81 static gboolean nlfilter_dialog (gint32 drawable_id);
82
83 static gint nlfiltInit (gdouble alpha,
84 gdouble radius,
85 FilterType filter);
86
87 static void nlfiltRow (guchar *srclast,
88 guchar *srcthis,
89 guchar *srcnext,
90 guchar *dst,
91 gint width,
92 gint bpp,
93 gint filtno);
94
95
96 const GimpPlugInInfo PLUG_IN_INFO =
97 {
98 NULL, /* init_proc */
99 NULL, /* quit_proc */
100 query, /* query_proc */
101 run, /* run_proc */
102 };
103
MAIN()104 MAIN ()
105
106 static void
107 query (void)
108 {
109 static const GimpParamDef args[] =
110 {
111 { GIMP_PDB_INT32, "run-mode", "The run mode { RUN-INTERACTIVE (0), RUN-NONINTERACTIVE (1) }" },
112 { GIMP_PDB_IMAGE, "img", "The Image to Filter" },
113 { GIMP_PDB_DRAWABLE, "drw", "The Drawable" },
114 { GIMP_PDB_FLOAT, "alpha", "The amount of the filter to apply" },
115 { GIMP_PDB_FLOAT, "radius", "The filter radius" },
116 { GIMP_PDB_INT32, "filter", "The Filter to Run, "
117 "0 - alpha trimmed mean; "
118 "1 - optimal estimation (alpha controls noise variance); "
119 "2 - edge enhancement" }
120 };
121
122 gimp_install_procedure (PLUG_IN_PROC,
123 N_("Nonlinear swiss army knife filter"),
124 "This is the pnmnlfilt, in gimp's clothing. "
125 "See the pnmnlfilt manpage for details.",
126 "Graeme W. Gill, gimp 0.99 plug-in by Eric L. Hernes",
127 "Graeme W. Gill, Eric L. Hernes",
128 "1997",
129 N_("_NL Filter..."),
130 "RGB,GRAY",
131 GIMP_PLUGIN,
132 G_N_ELEMENTS (args), 0,
133 args, NULL);
134
135 gimp_plugin_menu_register (PLUG_IN_PROC, "<Image>/Filters/Enhance");
136 }
137
138 static void
run(const gchar * name,gint nparams,const GimpParam * param,gint * nreturn_vals,GimpParam ** return_vals)139 run (const gchar *name,
140 gint nparams,
141 const GimpParam *param,
142 gint *nreturn_vals,
143 GimpParam **return_vals)
144 {
145 static GimpParam values[1];
146 GimpRunMode run_mode;
147 gint32 drawable_id;
148 GimpPDBStatusType status = GIMP_PDB_SUCCESS;
149
150 INIT_I18N ();
151 gegl_init (NULL, NULL);
152
153 run_mode = param[0].data.d_int32;
154 drawable_id = param[2].data.d_drawable;
155
156 *nreturn_vals = 1;
157 *return_vals = values;
158
159 values[0].type = GIMP_PDB_STATUS;
160 values[0].data.d_status = status;
161
162 switch (run_mode)
163 {
164 case GIMP_RUN_INTERACTIVE:
165 gimp_get_data (PLUG_IN_PROC, &nlfvals);
166
167 if (! nlfilter_dialog (drawable_id))
168 return;
169 break;
170
171 case GIMP_RUN_NONINTERACTIVE:
172 if (nparams != 6)
173 {
174 status = GIMP_PDB_CALLING_ERROR;
175 }
176 else
177 {
178 nlfvals.alpha = param[3].data.d_float;
179 nlfvals.radius = param[4].data.d_float;
180 nlfvals.filter = param[5].data.d_int32;
181 }
182
183 break;
184
185 case GIMP_RUN_WITH_LAST_VALS:
186 gimp_get_data (PLUG_IN_PROC, &nlfvals);
187 break;
188
189 default:
190 break;
191 }
192
193 if (status == GIMP_PDB_SUCCESS)
194 {
195 nlfilter (drawable_id, NULL);
196
197 /* Store data */
198 if (run_mode == GIMP_RUN_INTERACTIVE)
199 gimp_set_data (PLUG_IN_PROC, &nlfvals, sizeof (NLFilterValues));
200 }
201
202 values[0].data.d_status = status;
203 }
204
205 /* pnmnlfilt.c - 4 in 1 (2 non-linear) filter
206 ** - smooth an anyimage
207 ** - do alpha trimmed mean filtering on an anyimage
208 ** - do optimal estimation smoothing on an anyimage
209 ** - do edge enhancement on an anyimage
210 **
211 ** Version 1.0
212 **
213 ** The implementation of an alpha-trimmed mean filter
214 ** is based on the description in IEEE CG&A May 1990
215 ** Page 23 by Mark E. Lee and Richard A. Redner.
216 **
217 ** The paper recommends using a hexagon sampling region around each
218 ** pixel being processed, allowing an effective sub pixel radius to be
219 ** specified. The hexagon values are sythesised by area sampling the
220 ** rectangular pixels with a hexagon grid. The seven hexagon values
221 ** obtained from the 3x3 pixel grid are used to compute the alpha
222 ** trimmed mean. Note that an alpha value of 0.0 gives a conventional
223 ** mean filter (where the radius controls the contribution of
224 ** surrounding pixels), while a value of 0.5 gives a median filter.
225 ** Although there are only seven values to trim from before finding
226 ** the mean, the algorithm has been extended from that described in
227 ** CG&A by using interpolation, to allow a continuous selection of
228 ** alpha value between and including 0.0 to 0.5 The useful values
229 ** for radius are between 0.3333333 (where the filter will have no
230 ** effect because only one pixel is sampled), to 1.0, where all
231 ** pixels in the 3x3 grid are sampled.
232 **
233 ** The optimal estimation filter is taken from an article "Converting Dithered
234 ** Images Back to Gray Scale" by Allen Stenger, Dr Dobb's Journal, November
235 ** 1992, and this article references "Digital Image Enhancement andNoise Filtering by
236 ** Use of Local Statistics", Jong-Sen Lee, IEEE Transactions on Pattern Analysis and
237 ** Machine Intelligence, March 1980.
238 **
239 ** Also borrow the technique used in pgmenhance(1) to allow edge
240 ** enhancement if the alpha value is negative.
241 **
242 ** Author:
243 ** Graeme W. Gill, 30th Jan 1993
244 ** graeme@labtam.oz.au
245 **
246 ** Permission is hereby granted, to use, copy, modify, distribute,
247 ** and sell this software and its associated documentation files
248 ** (the "Software") for any purpose without fee, provided
249 ** that:
250 **
251 ** 1) The above copyright notices and this permission notice
252 ** accompany all source code copies of the Software and
253 ** related documentation.
254 ** and
255 **
256 ** 2) If executable code based on the Software only is distributed,
257 ** then the accompanying documentation must acknowledge that
258 ** "this software is based in part on the work of Graeme W. Gill".
259 ** and
260 **
261 ** 3) It is accepted that Graeme W. Gill (the "Author") accepts
262 ** NO LIABILITY for damages of any kind. The Software is
263 ** provided without fee by the Author "AS-IS" and without
264 ** warranty of any kind, express, implied or otherwise,
265 ** including without limitation, any warranty of merchantability
266 ** or fitness for a particular purpose.
267 ** and
268 **
269 ** 4) These conditions apply to any software derived from or based
270 ** on the Software, not just to the unmodified library.
271 **
272 */
273
274 /* ************************************************** */
275 /* Hexagon intersecting square area functions */
276 /* Compute the area of the intersection of a triangle */
277 /* and a rectangle */
278
279 static gdouble triang_area(gdouble, gdouble, gdouble, gdouble, gdouble,
280 gdouble, gdouble, gdouble, gint);
281 static gdouble rectang_area(gdouble, gdouble, gdouble, gdouble,
282 gdouble, gdouble, gdouble, gdouble);
283 static gdouble hex_area(gdouble, gdouble, gdouble, gdouble, gdouble);
284
285 static gint atfilt0 (gint *p);
286 static gint atfilt1 (gint *p);
287 static gint atfilt2 (gint *p);
288 static gint atfilt3 (gint *p);
289 static gint atfilt4 (gint *p);
290 static gint atfilt5 (gint *p);
291
292 gint (*atfuncs[6])(gint *) =
293 {
294 atfilt0,
295 atfilt1,
296 atfilt2,
297 atfilt3,
298 atfilt4,
299 atfilt5
300 };
301
302 static gint noisevariance;
303
304 #define MXIVAL 255 /* maximum input value */
305 #define NOIVAL (MXIVAL + 1) /* number of possible input values */
306
307 #define SCALEB 8 /* scale bits */
308 #define SCALE (1 << SCALEB) /* scale factor */
309
310 #define CSCALEB 2 /* coarse scale bits */
311 #define CSCALE (1 << CSCALEB) /* coarse scale factor */
312 #define MXCSVAL (MXIVAL * CSCALE) /* maximum coarse scaled values */
313 #define NOCSVAL (MXCSVAL + 1) /* number of coarse scaled values */
314 #define SCTOCSC(x) ((x) >> (SCALEB - CSCALEB)) /* convert from scaled to coarse scaled */
315 #define CSCTOSC(x) ((x) << (SCALEB - CSCALEB)) /* convert from course scaled to scaled */
316
317 /* round and scale floating point to scaled integer */
318 #define SROUND(x) ((gint)(((x) * (gdouble)SCALE) + 0.5))
319 /* round and un-scale scaled integer value */
320 #define RUNSCALE(x) (((x) + (1 << (SCALEB-1))) >> SCALEB) /* rounded un-scale */
321 #define UNSCALE(x) ((x) >> SCALEB)
322
323 /* Note: modified by David Hodson, nlfiltRow now accesses
324 * srclast, srcthis, and srcnext from [-bpp] to [width*bpp-1].
325 * Beware if you use this code anywhere else!
326 */
327 static void
nlfiltRow(guchar * srclast,guchar * srcthis,guchar * srcnext,guchar * dst,gint width,gint bpp,gint filtno)328 nlfiltRow (guchar *srclast, guchar *srcthis, guchar *srcnext, guchar *dst,
329 gint width, gint bpp, gint filtno)
330 {
331 gint pf[9];
332 guchar *ip0, *ip1, *ip2, *or, *orend;
333
334 orend = dst + width * bpp;
335 ip0 = srclast;
336 ip1 = srcthis;
337 ip2 = srcnext;
338
339 for (or = dst; or < orend; ip0++, ip1++, ip2++, or++)
340 {
341 pf[0] = *ip1;
342 pf[1] = *(ip1 - bpp);
343 pf[2] = *(ip2 - bpp);
344 pf[3] = *(ip2);
345 pf[4] = *(ip2 + bpp);
346 pf[5] = *(ip1 + bpp);
347 pf[6] = *(ip0 + bpp);
348 pf[7] = *(ip0);
349 pf[8] = *(ip0 - bpp);
350 *or=(atfuncs[filtno])(pf);
351 }
352 }
353
354 /* We restrict radius to the values: 0.333333 <= radius <= 1.0 */
355 /* so that no fewer and no more than a 3x3 grid of pixels around */
356 /* the pixel in question needs to be read. Given this, we only */
357 /* need 3 or 4 weightings per hexagon, as follows: */
358 /* _ _ */
359 /* Virtical hex: |_|_| 1 2 */
360 /* |X|_| 0 3 */
361 /* _ */
362 /* _ _|_| 1 */
363 /* Middle hex: |_| 1 Horizontal hex: |X|_| 0 2 */
364 /* |X| 0 |_| 3 */
365 /* |_| 2 */
366
367 /* all filters */
368 gint V0[NOIVAL],V1[NOIVAL],V2[NOIVAL],V3[NOIVAL]; /* vertical hex */
369 gint M0[NOIVAL],M1[NOIVAL],M2[NOIVAL]; /* middle hex */
370 gint H0[NOIVAL],H1[NOIVAL],H2[NOIVAL],H3[NOIVAL]; /* horizontal hex */
371
372 /* alpha trimmed and edge enhancement only */
373 gint ALFRAC[NOIVAL * 8]; /* fractional alpha divider table */
374
375 /* optimal estimation only */
376 gint AVEDIV[7 * NOCSVAL]; /* divide by 7 to give average value */
377 gint SQUARE[2 * NOCSVAL]; /* scaled square lookup table */
378
379 /* Table initialisation function - return alpha range */
380 static gint
nlfiltInit(gdouble alpha,gdouble radius,FilterType filter)381 nlfiltInit (gdouble alpha, gdouble radius, FilterType filter)
382 {
383 gint alpharange; /* alpha range value 0 - 3 */
384 gdouble meanscale; /* scale for finding mean */
385 gdouble mmeanscale; /* scale for finding mean - midle hex */
386 gdouble alphafraction; /* fraction of next largest/smallest
387 * to subtract from sum
388 */
389 switch (filter)
390 {
391 case filter_alpha_trim:
392 {
393 gdouble noinmean;
394 /* alpha only makes sense in range 0.0 - 0.5 */
395 alpha /= 2.0;
396 /* number of elements (out of a possible 7) used in the mean */
397 noinmean = ((0.5 - alpha) * 12.0) + 1.0;
398 mmeanscale = meanscale = 1.0/noinmean;
399 if (alpha == 0.0) { /* mean filter */
400 alpharange = 0;
401 alphafraction = 0.0; /* not used */
402 } else if (alpha < (1.0/6.0)) { /* mean of 5 to 7 middle values */
403 alpharange = 1;
404 alphafraction = (7.0 - noinmean)/2.0;
405 } else if (alpha < (1.0/3.0)) { /* mean of 3 to 5 middle values */
406 alpharange = 2;
407 alphafraction = (5.0 - noinmean)/2.0;
408 } else { /* mean of 1 to 3 middle values */
409 /* alpha==0.5 => median filter */
410 alpharange = 3;
411 alphafraction = (3.0 - noinmean)/2.0;
412 }
413 }
414 break;
415 case filter_opt_est: {
416 gint i;
417 gdouble noinmean = 7.0;
418
419 /* edge enhancement function */
420 alpharange = 5;
421
422 /* compute scaled hex values */
423 mmeanscale=meanscale=1.0;
424
425 /* Set up 1:1 division lookup - not used */
426 alphafraction=1.0/noinmean;
427
428 /* estimate of noise variance */
429 noisevariance = alpha * (gdouble)255;
430 noisevariance = noisevariance * noisevariance / 8.0;
431
432 /* set yp optimal estimation specific stuff */
433
434 for (i=0;i<(7*NOCSVAL);i++) { /* divide scaled value by 7 lookup */
435 AVEDIV[i] = CSCTOSC(i)/7; /* scaled divide by 7 */
436 }
437 /* compute square and rescale by
438 * (val >> (2 * SCALEB + 2)) table
439 */
440 for (i=0;i<(2*NOCSVAL);i++) {
441 gint val;
442 /* NOCSVAL offset to cope with -ve input values */
443 val = CSCTOSC(i - NOCSVAL);
444 SQUARE[i] = (val * val) >> (2 * SCALEB + 2);
445 }
446 }
447 break;
448 case filter_edge_enhance: {
449 if (alpha == 1.0) alpha = 0.99;
450 alpharange = 4;
451 /* mean of 7 and scaled by -alpha/(1-alpha) */
452 meanscale = 1.0 * (-alpha/((1.0 - alpha) * 7.0));
453
454 /* middle pixel has 1/(1-alpha) as well */
455 mmeanscale = 1.0 * (1.0/(1.0 - alpha) + meanscale);
456 alphafraction = 0.0; /* not used */
457 }
458 break;
459 default:
460 g_printerr ("unknown filter %d\n", filter);
461 return -1;
462 }
463 /*
464 * Setup pixel weighting tables -
465 * note we pre-compute mean division here too.
466 */
467 {
468 gint i;
469 gdouble hexhoff,hexvoff;
470 gdouble tabscale,mtabscale;
471 gdouble v0,v1,v2,v3,m0,m1,m2,h0,h1,h2,h3;
472
473 /* horizontal offset of virtical hex centers */
474 hexhoff = radius/2;
475
476 /* vertical offset of virtical hex centers */
477 hexvoff = 3.0 * radius/sqrt(12.0);
478
479 /*
480 * scale tables to normalise by hexagon
481 * area, and number of hexes used in mean
482 */
483 tabscale = meanscale / (radius * hexvoff);
484 mtabscale = mmeanscale / (radius * hexvoff);
485 v0 = hex_area(0.0,0.0,hexhoff,hexvoff,radius) * tabscale;
486 v1 = hex_area(0.0,1.0,hexhoff,hexvoff,radius) * tabscale;
487 v2 = hex_area(1.0,1.0,hexhoff,hexvoff,radius) * tabscale;
488 v3 = hex_area(1.0,0.0,hexhoff,hexvoff,radius) * tabscale;
489 m0 = hex_area(0.0,0.0,0.0,0.0,radius) * mtabscale;
490 m1 = hex_area(0.0,1.0,0.0,0.0,radius) * mtabscale;
491 m2 = hex_area(0.0,-1.0,0.0,0.0,radius) * mtabscale;
492 h0 = hex_area(0.0,0.0,radius,0.0,radius) * tabscale;
493 h1 = hex_area(1.0,1.0,radius,0.0,radius) * tabscale;
494 h2 = hex_area(1.0,0.0,radius,0.0,radius) * tabscale;
495 h3 = hex_area(1.0,-1.0,radius,0.0,radius) * tabscale;
496
497 for (i=0; i <= MXIVAL; i++) {
498 gdouble fi;
499 fi = (gdouble)i;
500 V0[i] = SROUND(fi * v0);
501 V1[i] = SROUND(fi * v1);
502 V2[i] = SROUND(fi * v2);
503 V3[i] = SROUND(fi * v3);
504 M0[i] = SROUND(fi * m0);
505 M1[i] = SROUND(fi * m1);
506 M2[i] = SROUND(fi * m2);
507 H0[i] = SROUND(fi * h0);
508 H1[i] = SROUND(fi * h1);
509 H2[i] = SROUND(fi * h2);
510 H3[i] = SROUND(fi * h3);
511 }
512 /* set up alpha fraction lookup table used on big/small */
513 for (i=0; i < (NOIVAL * 8); i++) {
514 ALFRAC[i] = SROUND((gdouble)i * alphafraction);
515 }
516 }
517 return alpharange;
518 }
519
520 /* Core pixel processing function - hand it 3x3 pixels and return result. */
521 /* Mean filter */
522 static gint
atfilt0(gint32 * p)523 atfilt0(gint32 *p)
524 {
525 gint retv;
526 /* map to scaled hexagon values */
527 retv = M0[p[0]] + M1[p[3]] + M2[p[7]];
528 retv += H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
529 retv += V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
530 retv += V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
531 retv += H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
532 retv += V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
533 retv += V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
534 return UNSCALE(retv);
535 }
536
537 /* Mean of 5 - 7 middle values */
538 static gint
atfilt1(gint32 * p)539 atfilt1 (gint32 *p)
540 {
541 gint h0,h1,h2,h3,h4,h5,h6; /* hexagon values 2 3 */
542 /* 1 0 4 */
543 /* 6 5 */
544 gint big,small;
545 /* map to scaled hexagon values */
546 h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
547 h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
548 h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
549 h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
550 h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
551 h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
552 h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
553 /* sum values and also discover the largest and smallest */
554 big = small = h0;
555 #define CHECK(xx) \
556 h0 += xx; \
557 if (xx > big) \
558 big = xx; \
559 else if (xx < small) \
560 small = xx;
561 CHECK(h1)
562 CHECK(h2)
563 CHECK(h3)
564 CHECK(h4)
565 CHECK(h5)
566 CHECK(h6)
567 #undef CHECK
568 /* Compute mean of middle 5-7 values */
569 return UNSCALE(h0 -ALFRAC[(big + small)>>SCALEB]);
570 }
571
572 /* Mean of 3 - 5 middle values */
573 static gint
atfilt2(gint32 * p)574 atfilt2 (gint32 *p)
575 {
576 gint h0,h1,h2,h3,h4,h5,h6; /* hexagon values 2 3 */
577 /* 1 0 4 */
578 /* 6 5 */
579 gint big0,big1,small0,small1;
580 /* map to scaled hexagon values */
581 h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
582 h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
583 h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
584 h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
585 h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
586 h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
587 h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
588 /* sum values and also discover the 2 largest and 2 smallest */
589 big0 = small0 = h0;
590 small1 = G_MAXINT;
591 big1 = 0;
592 #define CHECK(xx) \
593 h0 += xx; \
594 if (xx > big1) \
595 { \
596 if (xx > big0) \
597 { \
598 big1 = big0; \
599 big0 = xx; \
600 } \
601 else \
602 big1 = xx; \
603 } \
604 if (xx < small1) \
605 { \
606 if (xx < small0) \
607 { \
608 small1 = small0; \
609 small0 = xx; \
610 } \
611 else \
612 small1 = xx; \
613 }
614 CHECK(h1)
615 CHECK(h2)
616 CHECK(h3)
617 CHECK(h4)
618 CHECK(h5)
619 CHECK(h6)
620 #undef CHECK
621 /* Compute mean of middle 3-5 values */
622 return UNSCALE(h0 -big0 -small0 -ALFRAC[(big1 + small1)>>SCALEB]);
623 }
624
625 /*
626 * Mean of 1 - 3 middle values.
627 * If only 1 value, then this is a median filter.
628 */
629 static gint32
atfilt3(gint32 * p)630 atfilt3(gint32 *p)
631 {
632 gint h0,h1,h2,h3,h4,h5,h6; /* hexagon values 2 3 */
633 /* 1 0 4 */
634 /* 6 5 */
635 gint big0,big1,big2,small0,small1,small2;
636 /* map to scaled hexagon values */
637 h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
638 h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
639 h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
640 h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
641 h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
642 h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
643 h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
644 /* sum values and also discover the 3 largest and 3 smallest */
645 big0 = small0 = h0;
646 small1 = small2 = G_MAXINT;
647 big1 = big2 = 0;
648 #define CHECK(xx) \
649 h0 += xx; \
650 if (xx > big2) \
651 { \
652 if (xx > big1) \
653 { \
654 if (xx > big0) \
655 { \
656 big2 = big1; \
657 big1 = big0; \
658 big0 = xx; \
659 } \
660 else \
661 { \
662 big2 = big1; \
663 big1 = xx; \
664 } \
665 } \
666 else \
667 big2 = xx; \
668 } \
669 if (xx < small2) \
670 { \
671 if (xx < small1) \
672 { \
673 if (xx < small0) \
674 { \
675 small2 = small1; \
676 small1 = small0; \
677 small0 = xx; \
678 } \
679 else \
680 { \
681 small2 = small1; \
682 small1 = xx; \
683 } \
684 } \
685 else \
686 small2 = xx; \
687 }
688 CHECK(h1)
689 CHECK(h2)
690 CHECK(h3)
691 CHECK(h4)
692 CHECK(h5)
693 CHECK(h6)
694 #undef CHECK
695 /* Compute mean of middle 1-3 values */
696 return UNSCALE(h0-big0-big1-small0-small1-ALFRAC[(big2+small2)>>SCALEB]);
697 }
698
699 /* Edge enhancement */
700 static gint
atfilt4(gint * p)701 atfilt4 (gint *p)
702 {
703 gint hav;
704 /* map to scaled hexagon values and compute enhance value */
705 hav = M0[p[0]] + M1[p[3]] + M2[p[7]];
706 hav += H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
707 hav += V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
708 hav += V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
709 hav += H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
710 hav += V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
711 hav += V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
712 if (hav < 0)
713 hav = 0;
714 hav = UNSCALE(hav);
715 if (hav > (gdouble)255)
716 hav = (gdouble)255;
717 return hav;
718 }
719
720 /* Optimal estimation - do smoothing in inverse proportion */
721 /* to the local variance. */
722 /* notice we use the globals noisevariance */
723 gint
atfilt5(gint * p)724 atfilt5(gint *p) {
725 gint mean,variance,temp;
726 gint h0,h1,h2,h3,h4,h5,h6; /* hexagon values 2 3 */
727 /* 1 0 4 */
728 /* 6 5 */
729 /* map to scaled hexagon values */
730 h0 = M0[p[0]] + M1[p[3]] + M2[p[7]];
731 h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]];
732 h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]];
733 h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]];
734 h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]];
735 h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]];
736 h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]];
737 mean = h0 + h1 + h2 + h3 + h4 + h5 + h6;
738 /* compute scaled mean by dividing by 7 */
739 mean = AVEDIV[SCTOCSC(mean)];
740
741 /* compute scaled variance */
742 temp = (h1 - mean); variance = SQUARE[NOCSVAL + SCTOCSC(temp)];
743
744 /* and rescale to keep */
745 temp = (h2 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
746
747 /* within 32 bit limits */
748 temp = (h3 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
749 temp = (h4 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
750 temp = (h5 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
751 temp = (h6 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
752 /* (temp = h0 - mean) */
753 temp = (h0 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)];
754 if (variance != 0) /* avoid possible divide by 0 */
755 /* optimal estimate */
756 temp = mean + (variance * temp) / (variance + noisevariance);
757 else temp = h0;
758 if (temp < 0)
759 temp = 0;
760 temp = RUNSCALE(temp);
761 if (temp > (gdouble)255) temp = (gdouble)255;
762 return temp;
763 }
764
765
766 /* Triangle orientation is per geometric axes (not graphical axies) */
767
768 #define NW 0 /* North west triangle /| */
769 #define NE 1 /* North east triangle |\ */
770 #define SW 2 /* South west triangle \| */
771 #define SE 3 /* South east triangle |/ */
772 #define STH 2
773 #define EST 1
774
775 #define SWAPI(a,b) (t = a, a = -b, b = -t)
776
777 /* compute the area of overlap of a hexagon diameter d, */
778 /* centered at hx,hy, with a unit square of center sx,sy. */
779 static gdouble
hex_area(gdouble sx,gdouble sy,gdouble hx,gdouble hy,gdouble d)780 hex_area (gdouble sx, gdouble sy, gdouble hx, gdouble hy, gdouble d)
781 {
782 gdouble hx0,hx1,hx2,hy0,hy1,hy2,hy3;
783 gdouble sx0,sx1,sy0,sy1;
784
785 /* compute square co-ordinates */
786 sx0 = sx - 0.5;
787 sy0 = sy - 0.5;
788 sx1 = sx + 0.5;
789 sy1 = sy + 0.5;
790 /* compute hexagon co-ordinates */
791 hx0 = hx - d/2.0;
792 hx1 = hx;
793 hx2 = hx + d/2.0;
794 hy0 = hy - 0.5773502692 * d; /* d / sqrt(3) */
795 hy1 = hy - 0.2886751346 * d; /* d / sqrt(12) */
796 hy2 = hy + 0.2886751346 * d; /* d / sqrt(12) */
797 hy3 = hy + 0.5773502692 * d; /* d / sqrt(3) */
798
799 return
800 triang_area(sx0,sy0,sx1,sy1,hx0,hy2,hx1,hy3,NW) +
801 triang_area(sx0,sy0,sx1,sy1,hx1,hy2,hx2,hy3,NE) +
802 rectang_area(sx0,sy0,sx1,sy1,hx0,hy1,hx2,hy2) +
803 triang_area(sx0,sy0,sx1,sy1,hx0,hy0,hx1,hy1,SW) +
804 triang_area(sx0,sy0,sx1,sy1,hx1,hy0,hx2,hy1,SE);
805 }
806
807 static gdouble
triang_area(gdouble rx0,gdouble ry0,gdouble rx1,gdouble ry1,gdouble tx0,gdouble ty0,gdouble tx1,gdouble ty1,gint tt)808 triang_area (gdouble rx0, gdouble ry0, gdouble rx1, gdouble ry1, gdouble tx0,
809 gdouble ty0, gdouble tx1, gdouble ty1, gint tt)
810 {
811 gdouble a,b,c,d;
812 gdouble lx0,ly0,lx1,ly1;
813 /* Convert everything to a NW triangle */
814 if (tt & STH) {
815 gdouble t;
816 SWAPI(ry0,ry1);
817 SWAPI(ty0,ty1);
818 } if (tt & EST) {
819 gdouble t;
820 SWAPI(rx0,rx1);
821 SWAPI(tx0,tx1);
822 }
823 /* Compute overlapping box */
824 if (tx0 > rx0)
825 rx0 = tx0;
826 if (ty0 > ry0)
827 ry0 = ty0;
828 if (tx1 < rx1)
829 rx1 = tx1;
830 if (ty1 < ry1)
831 ry1 = ty1;
832 if (rx1 <= rx0 || ry1 <= ry0)
833 return 0.0;
834 /* Need to compute diagonal line intersection with the box */
835 /* First compute co-efficients to formulas x = a + by and y = c + dx */
836 b = (tx1 - tx0)/(ty1 - ty0);
837 a = tx0 - b * ty0;
838 d = (ty1 - ty0)/(tx1 - tx0);
839 c = ty0 - d * tx0;
840
841 /* compute top or right intersection */
842 tt = 0;
843 ly1 = ry1;
844 lx1 = a + b * ly1;
845 if (lx1 <= rx0)
846 return (rx1 - rx0) * (ry1 - ry0);
847 else if (lx1 > rx1) { /* could be right hand side */
848 lx1 = rx1;
849 ly1 = c + d * lx1;
850 if (ly1 <= ry0)
851 return (rx1 - rx0) * (ry1 - ry0);
852 tt = 1; /* right hand side intersection */
853 }
854 /* compute left or bottom intersection */
855 lx0 = rx0;
856 ly0 = c + d * lx0;
857 if (ly0 >= ry1)
858 return (rx1 - rx0) * (ry1 - ry0);
859 else if (ly0 < ry0) { /* could be right hand side */
860 ly0 = ry0;
861 lx0 = a + b * ly0;
862 if (lx0 >= rx1)
863 return (rx1 - rx0) * (ry1 - ry0);
864 tt |= 2; /* bottom intersection */
865 }
866
867 if (tt == 0) { /* top and left intersection */
868 /* rectangle minus triangle */
869 return ((rx1 - rx0) * (ry1 - ry0))
870 - (0.5 * (lx1 - rx0) * (ry1 - ly0));
871 }
872 else if (tt == 1) { /* right and left intersection */
873 return ((rx1 - rx0) * (ly0 - ry0))
874 + (0.5 * (rx1 - rx0) * (ly1 - ly0));
875 } else if (tt == 2) { /* top and bottom intersection */
876 return ((rx1 - lx1) * (ry1 - ry0))
877 + (0.5 * (lx1 - lx0) * (ry1 - ry0));
878 } else { /* tt == 3 */ /* right and bottom intersection */
879 /* triangle */
880 return (0.5 * (rx1 - lx0) * (ly1 - ry0));
881 }
882 }
883
884 /* Compute rectangle area */
885 static gdouble
rectang_area(gdouble rx0,gdouble ry0,gdouble rx1,gdouble ry1,gdouble tx0,gdouble ty0,gdouble tx1,gdouble ty1)886 rectang_area (gdouble rx0, gdouble ry0, gdouble rx1, gdouble ry1, gdouble tx0,
887 gdouble ty0, gdouble tx1, gdouble ty1)
888 {
889 /* Compute overlapping box */
890 if (tx0 > rx0)
891 rx0 = tx0;
892 if (ty0 > ry0)
893 ry0 = ty0;
894 if (tx1 < rx1)
895 rx1 = tx1;
896 if (ty1 < ry1)
897 ry1 = ty1;
898 if (rx1 <= rx0 || ry1 <= ry0)
899 return 0.0;
900 return (rx1 - rx0) * (ry1 - ry0);
901 }
902
903 static void
nlfilter(gint32 drawable_id,GimpPreview * preview)904 nlfilter (gint32 drawable_id,
905 GimpPreview *preview)
906 {
907 GeglBuffer *src_buffer;
908 GeglBuffer *dest_buffer;
909 const Babl *format;
910 guchar *srcbuf, *dstbuf;
911 guchar *lastrow, *thisrow, *nextrow, *temprow;
912 gint x1, y1, y2;
913 gint width, height, bpp;
914 gint filtno, y, rowsize, exrowsize, p_update;
915
916 if (preview)
917 {
918 gimp_preview_get_position (preview, &x1, &y1);
919 gimp_preview_get_size (preview, &width, &height);
920 y2 = y1 + height;
921 }
922 else
923 {
924 if (! gimp_drawable_mask_intersect (drawable_id,
925 &x1, &y1, &width, &height))
926 return;
927
928 y2 = y1 + height;
929 }
930
931 if (gimp_drawable_has_alpha (drawable_id))
932 format = babl_format ("R'G'B'A u8");
933 else
934 format = babl_format ("R'G'B' u8");
935
936 src_buffer = gimp_drawable_get_buffer (drawable_id);
937
938 if (preview)
939 dest_buffer = gegl_buffer_new (gegl_buffer_get_extent (src_buffer), format);
940 else
941 dest_buffer = gimp_drawable_get_shadow_buffer (drawable_id);
942
943 bpp = babl_format_get_bytes_per_pixel (format);
944
945 rowsize = width * bpp;
946 exrowsize = (width + 2) * bpp;
947 p_update = width / 20 + 1;
948
949 /* source buffer gives one pixel margin all around destination buffer */
950 srcbuf = g_new0 (guchar, exrowsize * 3);
951 dstbuf = g_new0 (guchar, rowsize);
952
953 /* pointers to second pixel in each source row */
954 lastrow = srcbuf + bpp;
955 thisrow = lastrow + exrowsize;
956 nextrow = thisrow + exrowsize;
957
958 filtno = nlfiltInit (nlfvals.alpha, nlfvals.radius, nlfvals.filter);
959
960 if (!preview)
961 gimp_progress_init (_("NL Filter"));
962
963 /* first row */
964 gegl_buffer_get (src_buffer, GEGL_RECTANGLE (x1, y1, width, 1), 1.0,
965 format, thisrow,
966 GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_NONE);
967
968 /* copy thisrow[0] to thisrow[-1], thisrow[width-1] to thisrow[width] */
969 memcpy (thisrow - bpp, thisrow, bpp);
970 memcpy (thisrow + rowsize, thisrow + rowsize - bpp, bpp);
971 /* copy whole thisrow to lastrow */
972 memcpy (lastrow - bpp, thisrow - bpp, exrowsize);
973
974 for (y = y1; y < y2 - 1; y++)
975 {
976 if (((y % p_update) == 0) && !preview)
977 gimp_progress_update ((gdouble) y / (gdouble) height);
978
979 gegl_buffer_get (src_buffer, GEGL_RECTANGLE (x1, y + 1, width, 1), 1.0,
980 format, nextrow,
981 GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_NONE);
982
983 memcpy (nextrow - bpp, nextrow, bpp);
984 memcpy (nextrow + rowsize, nextrow + rowsize - bpp, bpp);
985 nlfiltRow (lastrow, thisrow, nextrow, dstbuf, width, bpp, filtno);
986
987 gegl_buffer_set (dest_buffer, GEGL_RECTANGLE (x1, y, width, 1), 0,
988 format, dstbuf,
989 GEGL_AUTO_ROWSTRIDE);
990
991 /* rotate row buffers */
992 temprow = lastrow; lastrow = thisrow;
993 thisrow = nextrow; nextrow = temprow;
994 }
995
996 /* last row */
997 memcpy (nextrow - bpp, thisrow - bpp, exrowsize);
998 nlfiltRow (lastrow, thisrow, nextrow, dstbuf, width, bpp, filtno);
999
1000 gegl_buffer_set (dest_buffer, GEGL_RECTANGLE (x1, y2 - 1, width, 1), 0,
1001 format, dstbuf,
1002 GEGL_AUTO_ROWSTRIDE);
1003
1004 g_free (srcbuf);
1005 g_free (dstbuf);
1006
1007 if (preview)
1008 {
1009 guchar *buf = g_new (guchar, width * height * bpp);
1010
1011 gegl_buffer_get (dest_buffer, GEGL_RECTANGLE (x1, y1, width, height), 1.0,
1012 format, buf,
1013 GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_NONE);
1014
1015 gimp_preview_draw_buffer (GIMP_PREVIEW (preview), buf, width * bpp);
1016
1017 g_free (buf);
1018 }
1019
1020 g_object_unref (src_buffer);
1021 g_object_unref (dest_buffer);
1022
1023 if (! preview)
1024 {
1025 gimp_progress_update (1.0);
1026
1027 gimp_drawable_merge_shadow (drawable_id, TRUE);
1028 gimp_drawable_update (drawable_id, x1, y1, width, height);
1029 gimp_displays_flush ();
1030 }
1031 }
1032
1033 static void
nlfilter_preview(gpointer drawable_id,GimpPreview * preview)1034 nlfilter_preview (gpointer drawable_id,
1035 GimpPreview *preview)
1036 {
1037 nlfilter (GPOINTER_TO_INT (drawable_id), preview);
1038 }
1039
1040 static gboolean
nlfilter_dialog(gint32 drawable_id)1041 nlfilter_dialog (gint32 drawable_id)
1042 {
1043 GtkWidget *dialog;
1044 GtkWidget *main_vbox;
1045 GtkWidget *preview;
1046 GtkWidget *frame;
1047 GtkWidget *alpha_trim;
1048 GtkWidget *opt_est;
1049 GtkWidget *edge_enhance;
1050 GtkWidget *table;
1051 GtkObject *adj;
1052 gboolean run;
1053
1054 gimp_ui_init (PLUG_IN_BINARY, TRUE);
1055
1056 dialog = gimp_dialog_new (_("NL Filter"), PLUG_IN_ROLE,
1057 NULL, 0,
1058 gimp_standard_help_func, PLUG_IN_PROC,
1059
1060 _("_Cancel"), GTK_RESPONSE_CANCEL,
1061 _("_OK"), GTK_RESPONSE_OK,
1062
1063 NULL);
1064
1065 gtk_dialog_set_alternative_button_order (GTK_DIALOG (dialog),
1066 GTK_RESPONSE_OK,
1067 GTK_RESPONSE_CANCEL,
1068 -1);
1069
1070 gimp_window_set_transient (GTK_WINDOW (dialog));
1071
1072 main_vbox = gtk_box_new (GTK_ORIENTATION_VERTICAL, 12);
1073 gtk_container_set_border_width (GTK_CONTAINER (main_vbox), 12);
1074 gtk_box_pack_start (GTK_BOX (gtk_dialog_get_content_area (GTK_DIALOG (dialog))),
1075 main_vbox, TRUE, TRUE, 0);
1076 gtk_widget_show (main_vbox);
1077
1078 preview = gimp_drawable_preview_new_from_drawable_id (drawable_id);
1079 gtk_box_pack_start (GTK_BOX (main_vbox), preview, TRUE, TRUE, 0);
1080 gtk_widget_show (preview);
1081
1082 g_signal_connect_swapped (preview, "invalidated",
1083 G_CALLBACK (nlfilter_preview),
1084 GINT_TO_POINTER (drawable_id));
1085
1086 frame = gimp_int_radio_group_new (TRUE, _("Filter"),
1087 G_CALLBACK (gimp_radio_button_update),
1088 &nlfvals.filter, nlfvals.filter,
1089
1090 _("_Alpha trimmed mean"),
1091 filter_alpha_trim, &alpha_trim,
1092 _("Op_timal estimation"),
1093 filter_opt_est, &opt_est,
1094 _("_Edge enhancement"),
1095 filter_edge_enhance, &edge_enhance,
1096
1097 NULL);
1098
1099 gtk_box_pack_start (GTK_BOX (main_vbox), frame, FALSE, FALSE, 0);
1100 gtk_widget_show (frame);
1101
1102 g_signal_connect_swapped (alpha_trim, "toggled",
1103 G_CALLBACK (gimp_preview_invalidate),
1104 preview);
1105 g_signal_connect_swapped (opt_est, "toggled",
1106 G_CALLBACK (gimp_preview_invalidate),
1107 preview);
1108 g_signal_connect_swapped (edge_enhance, "toggled",
1109 G_CALLBACK (gimp_preview_invalidate),
1110 preview);
1111
1112 table = gtk_table_new (2, 3, FALSE);
1113 gtk_table_set_col_spacings (GTK_TABLE (table), 6);
1114 gtk_table_set_row_spacings (GTK_TABLE (table), 6);
1115 gtk_box_pack_start (GTK_BOX (main_vbox), table, FALSE, FALSE, 0);
1116 gtk_widget_show (table);
1117
1118 adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 0,
1119 _("A_lpha:"), 0, 0,
1120 nlfvals.alpha, 0.0, 1.0, 0.05, 0.1, 2,
1121 TRUE, 0, 0,
1122 NULL, NULL);
1123 g_signal_connect (adj, "value-changed",
1124 G_CALLBACK (gimp_double_adjustment_update),
1125 &nlfvals.alpha);
1126 g_signal_connect_swapped (adj, "value-changed",
1127 G_CALLBACK (gimp_preview_invalidate),
1128 preview);
1129
1130 adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 1,
1131 _("_Radius:"), 0, 0,
1132 nlfvals.radius, 1.0 / 3.0, 1.0, 0.05, 0.1, 2,
1133 TRUE, 0, 0,
1134 NULL, NULL);
1135 g_signal_connect (adj, "value-changed",
1136 G_CALLBACK (gimp_double_adjustment_update),
1137 &nlfvals.radius);
1138 g_signal_connect_swapped (adj, "value-changed",
1139 G_CALLBACK (gimp_preview_invalidate),
1140 preview);
1141
1142 gtk_widget_show (dialog);
1143
1144 run = (gimp_dialog_run (GIMP_DIALOG (dialog)) == GTK_RESPONSE_OK);
1145
1146 gtk_widget_destroy (dialog);
1147
1148 return run;
1149 }
1150