/************************************************** * file: nlfilt/nlfilt.c * * Copyright (c) 1997 Eric L. Hernes (erich@rrnet.com) * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. The name of the author may not be used to endorse or promote products * derived from this software without specific prior written permission * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ /* * Algorithm fixes, V2.0 compatibility by David Hodson hodsond@ozemail.com.au */ #include "config.h" #include #include #include #include "libpika/stdplugins-intl.h" #define PLUG_IN_PROC "plug-in-nl-filter" #define PLUG_IN_BINARY "nl-filter" #define PLUG_IN_ROLE "pika-nl-filter" typedef struct { gdouble alpha; gdouble radius; gint filter; } NLFilterValues; typedef enum { filter_alpha_trim, filter_opt_est, filter_edge_enhance } FilterType; typedef struct _Nlfilter Nlfilter; typedef struct _NlfilterClass NlfilterClass; struct _Nlfilter { PikaPlugIn parent_instance; }; struct _NlfilterClass { PikaPlugInClass parent_class; }; #define NLFILTER_TYPE (nlfilter_get_type ()) #define NLFILTER(obj) (G_TYPE_CHECK_INSTANCE_CAST ((obj), NLFILTER_TYPE, Nlfilter)) GType nlfilter_get_type (void) G_GNUC_CONST; static GList * nlfilter_query_procedures (PikaPlugIn *plug_in); static PikaProcedure * nlfilter_create_procedure (PikaPlugIn *plug_in, const gchar *name); static PikaValueArray * nlfilter_run (PikaProcedure *procedure, PikaRunMode run_mode, PikaImage *image, gint n_drawables, PikaDrawable **drawables, PikaProcedureConfig *config, gpointer run_data); static void nlfilter (GObject *config, PikaDrawable *drawable, PikaPreview *preview); static void nlfilter_preview (GtkWidget *widget, GObject *config); static gboolean nlfilter_dialog (PikaProcedure *procedure, GObject *config, PikaDrawable *drawable); static gint nlfiltInit (gdouble alpha, gdouble radius, FilterType filter); static void nlfiltRow (guchar *srclast, guchar *srcthis, guchar *srcnext, guchar *dst, gint width, gint bpp, gint filtno); G_DEFINE_TYPE (Nlfilter, nlfilter, PIKA_TYPE_PLUG_IN) PIKA_MAIN (NLFILTER_TYPE) DEFINE_STD_SET_I18N static void nlfilter_class_init (NlfilterClass *klass) { PikaPlugInClass *plug_in_class = PIKA_PLUG_IN_CLASS (klass); plug_in_class->query_procedures = nlfilter_query_procedures; plug_in_class->create_procedure = nlfilter_create_procedure; plug_in_class->set_i18n = STD_SET_I18N; } static void nlfilter_init (Nlfilter *nlfilter) { } static GList * nlfilter_query_procedures (PikaPlugIn *plug_in) { return g_list_append (NULL, g_strdup (PLUG_IN_PROC)); } static PikaProcedure * nlfilter_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, nlfilter_run, NULL, NULL); pika_procedure_set_image_types (procedure, "RGB, GRAY"); pika_procedure_set_sensitivity_mask (procedure, PIKA_PROCEDURE_SENSITIVE_DRAWABLE); pika_procedure_set_menu_label (procedure, _("_NL Filter...")); pika_procedure_add_menu_path (procedure,"/Filters/Enhance"); pika_procedure_set_documentation (procedure, _("Nonlinear swiss army knife filter"), _("This is the pnmnlfilt, in PIKA's " "clothing. See the pnmnlfilt manpage " "for details."), name); pika_procedure_set_attribution (procedure, "Graeme W. Gill, PIKA 0.99 plug-in " "by Eric L. Hernes", "Graeme W. Gill, Eric L. Hernes", "1997"); PIKA_PROC_ARG_DOUBLE (procedure, "alpha", _("_Alpha"), _("The amount of the filter to apply"), 0, 1, 0.3, G_PARAM_READWRITE); PIKA_PROC_ARG_DOUBLE (procedure, "radius", _("Ra_dius"), _("The filter radius"), 1.0 / 3.0, 1, 1.0 / 3.0, G_PARAM_READWRITE); PIKA_PROC_ARG_INT (procedure, "filter", _("Filter"), _("The Filter to Run, " "0 - alpha trimmed mean; " "1 - optimal estimation " "(alpha controls noise variance); " "2 - edge enhancement"), 0, 2, 0, G_PARAM_READWRITE); } return procedure; } static PikaValueArray * nlfilter_run (PikaProcedure *procedure, PikaRunMode run_mode, PikaImage *image, gint n_drawables, PikaDrawable **drawables, PikaProcedureConfig *config, gpointer run_data) { PikaDrawable *drawable; 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."), pika_procedure_get_name (procedure)); return pika_procedure_new_return_values (procedure, PIKA_PDB_CALLING_ERROR, error); } else { drawable = drawables[0]; } if (run_mode == PIKA_RUN_INTERACTIVE && ! nlfilter_dialog (procedure, G_OBJECT (config), drawable)) return pika_procedure_new_return_values (procedure, PIKA_PDB_CANCEL, NULL); nlfilter (G_OBJECT (config), drawable, NULL); if (run_mode != PIKA_RUN_NONINTERACTIVE) pika_displays_flush (); return pika_procedure_new_return_values (procedure, PIKA_PDB_SUCCESS, NULL); } /* pnmnlfilt.c - 4 in 1 (2 non-linear) filter ** - smooth an anyimage ** - do alpha trimmed mean filtering on an anyimage ** - do optimal estimation smoothing on an anyimage ** - do edge enhancement on an anyimage ** ** Version 1.0 ** ** The implementation of an alpha-trimmed mean filter ** is based on the description in IEEE CG&A May 1990 ** Page 23 by Mark E. Lee and Richard A. Redner. ** ** The paper recommends using a hexagon sampling region around each ** pixel being processed, allowing an effective sub pixel radius to be ** specified. The hexagon values are sythesised by area sampling the ** rectangular pixels with a hexagon grid. The seven hexagon values ** obtained from the 3x3 pixel grid are used to compute the alpha ** trimmed mean. Note that an alpha value of 0.0 gives a conventional ** mean filter (where the radius controls the contribution of ** surrounding pixels), while a value of 0.5 gives a median filter. ** Although there are only seven values to trim from before finding ** the mean, the algorithm has been extended from that described in ** CG&A by using interpolation, to allow a continuous selection of ** alpha value between and including 0.0 to 0.5 The useful values ** for radius are between 0.3333333 (where the filter will have no ** effect because only one pixel is sampled), to 1.0, where all ** pixels in the 3x3 grid are sampled. ** ** The optimal estimation filter is taken from an article "Converting Dithered ** Images Back to Gray Scale" by Allen Stenger, Dr Dobb's Journal, November ** 1992, and this article references "Digital Image Enhancement andNoise Filtering by ** Use of Local Statistics", Jong-Sen Lee, IEEE Transactions on Pattern Analysis and ** Machine Intelligence, March 1980. ** ** Also borrow the technique used in pgmenhance(1) to allow edge ** enhancement if the alpha value is negative. ** ** Author: ** Graeme W. Gill, 30th Jan 1993 ** graeme@labtam.oz.au ** ** Permission is hereby granted, to use, copy, modify, distribute, ** and sell this software and its associated documentation files ** (the "Software") for any purpose without fee, provided ** that: ** ** 1) The above copyright notices and this permission notice ** accompany all source code copies of the Software and ** related documentation. ** and ** ** 2) If executable code based on the Software only is distributed, ** then the accompanying documentation must acknowledge that ** "this software is based in part on the work of Graeme W. Gill". ** and ** ** 3) It is accepted that Graeme W. Gill (the "Author") accepts ** NO LIABILITY for damages of any kind. The Software is ** provided without fee by the Author "AS-IS" and without ** warranty of any kind, express, implied or otherwise, ** including without limitation, any warranty of merchantability ** or fitness for a particular purpose. ** and ** ** 4) These conditions apply to any software derived from or based ** on the Software, not just to the unmodified library. ** */ /* ************************************************** */ /* Hexagon intersecting square area functions */ /* Compute the area of the intersection of a triangle */ /* and a rectangle */ static gdouble triang_area(gdouble, gdouble, gdouble, gdouble, gdouble, gdouble, gdouble, gdouble, gint); static gdouble rectang_area(gdouble, gdouble, gdouble, gdouble, gdouble, gdouble, gdouble, gdouble); static gdouble hex_area(gdouble, gdouble, gdouble, gdouble, gdouble); static gint atfilt0 (gint *p); static gint atfilt1 (gint *p); static gint atfilt2 (gint *p); static gint atfilt3 (gint *p); static gint atfilt4 (gint *p); static gint atfilt5 (gint *p); gint (*atfuncs[6])(gint *) = { atfilt0, atfilt1, atfilt2, atfilt3, atfilt4, atfilt5 }; static gint noisevariance; #define MXIVAL 255 /* maximum input value */ #define NOIVAL (MXIVAL + 1) /* number of possible input values */ #define SCALEB 8 /* scale bits */ #define SCALE (1 << SCALEB) /* scale factor */ #define CSCALEB 2 /* coarse scale bits */ #define CSCALE (1 << CSCALEB) /* coarse scale factor */ #define MXCSVAL (MXIVAL * CSCALE) /* maximum coarse scaled values */ #define NOCSVAL (MXCSVAL + 1) /* number of coarse scaled values */ #define SCTOCSC(x) ((x) >> (SCALEB - CSCALEB)) /* convert from scaled to coarse scaled */ #define CSCTOSC(x) ((x) << (SCALEB - CSCALEB)) /* convert from course scaled to scaled */ /* round and scale floating point to scaled integer */ #define SROUND(x) ((gint)(((x) * (gdouble)SCALE) + 0.5)) /* round and un-scale scaled integer value */ #define RUNSCALE(x) (((x) + (1 << (SCALEB-1))) >> SCALEB) /* rounded un-scale */ #define UNSCALE(x) ((x) >> SCALEB) /* Note: modified by David Hodson, nlfiltRow now accesses * srclast, srcthis, and srcnext from [-bpp] to [width*bpp-1]. * Beware if you use this code anywhere else! */ static void nlfiltRow (guchar *srclast, guchar *srcthis, guchar *srcnext, guchar *dst, gint width, gint bpp, gint filtno) { gint pf[9]; guchar *ip0, *ip1, *ip2, *or, *orend; orend = dst + width * bpp; ip0 = srclast; ip1 = srcthis; ip2 = srcnext; for (or = dst; or < orend; ip0++, ip1++, ip2++, or++) { pf[0] = *ip1; pf[1] = *(ip1 - bpp); pf[2] = *(ip2 - bpp); pf[3] = *(ip2); pf[4] = *(ip2 + bpp); pf[5] = *(ip1 + bpp); pf[6] = *(ip0 + bpp); pf[7] = *(ip0); pf[8] = *(ip0 - bpp); *or=(atfuncs[filtno])(pf); } } /* We restrict radius to the values: 0.333333 <= radius <= 1.0 */ /* so that no fewer and no more than a 3x3 grid of pixels around */ /* the pixel in question needs to be read. Given this, we only */ /* need 3 or 4 weightings per hexagon, as follows: */ /* _ _ */ /* Vertical hex: |_|_| 1 2 */ /* |X|_| 0 3 */ /* _ */ /* _ _|_| 1 */ /* Middle hex: |_| 1 Horizontal hex: |X|_| 0 2 */ /* |X| 0 |_| 3 */ /* |_| 2 */ /* all filters */ gint V0[NOIVAL],V1[NOIVAL],V2[NOIVAL],V3[NOIVAL]; /* vertical hex */ gint M0[NOIVAL],M1[NOIVAL],M2[NOIVAL]; /* middle hex */ gint H0[NOIVAL],H1[NOIVAL],H2[NOIVAL],H3[NOIVAL]; /* horizontal hex */ /* alpha trimmed and edge enhancement only */ gint ALFRAC[NOIVAL * 8]; /* fractional alpha divider table */ /* optimal estimation only */ gint AVEDIV[7 * NOCSVAL]; /* divide by 7 to give average value */ gint SQUARE[2 * NOCSVAL]; /* scaled square lookup table */ /* Table initialisation function - return alpha range */ static gint nlfiltInit (gdouble alpha, gdouble radius, FilterType filter) { gint alpharange; /* alpha range value 0 - 3 */ gdouble meanscale; /* scale for finding mean */ gdouble mmeanscale; /* scale for finding mean - midle hex */ gdouble alphafraction; /* fraction of next largest/smallest * to subtract from sum */ switch (filter) { case filter_alpha_trim: { gdouble noinmean; /* alpha only makes sense in range 0.0 - 0.5 */ alpha /= 2.0; /* number of elements (out of a possible 7) used in the mean */ noinmean = ((0.5 - alpha) * 12.0) + 1.0; mmeanscale = meanscale = 1.0/noinmean; if (alpha == 0.0) { /* mean filter */ alpharange = 0; alphafraction = 0.0; /* not used */ } else if (alpha < (1.0/6.0)) { /* mean of 5 to 7 middle values */ alpharange = 1; alphafraction = (7.0 - noinmean)/2.0; } else if (alpha < (1.0/3.0)) { /* mean of 3 to 5 middle values */ alpharange = 2; alphafraction = (5.0 - noinmean)/2.0; } else { /* mean of 1 to 3 middle values */ /* alpha==0.5 => median filter */ alpharange = 3; alphafraction = (3.0 - noinmean)/2.0; } } break; case filter_opt_est: { gint i; gdouble noinmean = 7.0; /* edge enhancement function */ alpharange = 5; /* compute scaled hex values */ mmeanscale=meanscale=1.0; /* Set up 1:1 division lookup - not used */ alphafraction=1.0/noinmean; /* estimate of noise variance */ noisevariance = alpha * (gdouble)255; noisevariance = noisevariance * noisevariance / 8.0; /* set yp optimal estimation specific stuff */ for (i=0;i<(7*NOCSVAL);i++) { /* divide scaled value by 7 lookup */ AVEDIV[i] = CSCTOSC(i)/7; /* scaled divide by 7 */ } /* compute square and rescale by * (val >> (2 * SCALEB + 2)) table */ for (i=0;i<(2*NOCSVAL);i++) { gint val; /* NOCSVAL offset to cope with -ve input values */ val = CSCTOSC(i - NOCSVAL); SQUARE[i] = (val * val) >> (2 * SCALEB + 2); } } break; case filter_edge_enhance: { if (alpha == 1.0) alpha = 0.99; alpharange = 4; /* mean of 7 and scaled by -alpha/(1-alpha) */ meanscale = 1.0 * (-alpha/((1.0 - alpha) * 7.0)); /* middle pixel has 1/(1-alpha) as well */ mmeanscale = 1.0 * (1.0/(1.0 - alpha) + meanscale); alphafraction = 0.0; /* not used */ } break; default: g_printerr ("unknown filter %d\n", filter); return -1; } /* * Setup pixel weighting tables - * note we pre-compute mean division here too. */ { gint i; gdouble hexhoff,hexvoff; gdouble tabscale,mtabscale; gdouble v0,v1,v2,v3,m0,m1,m2,h0,h1,h2,h3; /* horizontal offset of vertical hex centers */ hexhoff = radius/2; /* vertical offset of vertical hex centers */ hexvoff = 3.0 * radius/sqrt(12.0); /* * scale tables to normalise by hexagon * area, and number of hexes used in mean */ tabscale = meanscale / (radius * hexvoff); mtabscale = mmeanscale / (radius * hexvoff); v0 = hex_area(0.0,0.0,hexhoff,hexvoff,radius) * tabscale; v1 = hex_area(0.0,1.0,hexhoff,hexvoff,radius) * tabscale; v2 = hex_area(1.0,1.0,hexhoff,hexvoff,radius) * tabscale; v3 = hex_area(1.0,0.0,hexhoff,hexvoff,radius) * tabscale; m0 = hex_area(0.0,0.0,0.0,0.0,radius) * mtabscale; m1 = hex_area(0.0,1.0,0.0,0.0,radius) * mtabscale; m2 = hex_area(0.0,-1.0,0.0,0.0,radius) * mtabscale; h0 = hex_area(0.0,0.0,radius,0.0,radius) * tabscale; h1 = hex_area(1.0,1.0,radius,0.0,radius) * tabscale; h2 = hex_area(1.0,0.0,radius,0.0,radius) * tabscale; h3 = hex_area(1.0,-1.0,radius,0.0,radius) * tabscale; for (i=0; i <= MXIVAL; i++) { gdouble fi; fi = (gdouble)i; V0[i] = SROUND(fi * v0); V1[i] = SROUND(fi * v1); V2[i] = SROUND(fi * v2); V3[i] = SROUND(fi * v3); M0[i] = SROUND(fi * m0); M1[i] = SROUND(fi * m1); M2[i] = SROUND(fi * m2); H0[i] = SROUND(fi * h0); H1[i] = SROUND(fi * h1); H2[i] = SROUND(fi * h2); H3[i] = SROUND(fi * h3); } /* set up alpha fraction lookup table used on big/small */ for (i=0; i < (NOIVAL * 8); i++) { ALFRAC[i] = SROUND((gdouble)i * alphafraction); } } return alpharange; } /* Core pixel processing function - hand it 3x3 pixels and return result. */ /* Mean filter */ static gint atfilt0(gint32 *p) { gint retv; /* map to scaled hexagon values */ retv = M0[p[0]] + M1[p[3]] + M2[p[7]]; retv += H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]]; retv += V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]]; retv += V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]]; retv += H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]]; retv += V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]]; retv += V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]]; return UNSCALE(retv); } /* Mean of 5 - 7 middle values */ static gint atfilt1 (gint32 *p) { gint h0,h1,h2,h3,h4,h5,h6; /* hexagon values 2 3 */ /* 1 0 4 */ /* 6 5 */ gint big,small; /* map to scaled hexagon values */ h0 = M0[p[0]] + M1[p[3]] + M2[p[7]]; h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]]; h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]]; h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]]; h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]]; h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]]; h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]]; /* sum values and also discover the largest and smallest */ big = small = h0; #define CHECK(xx) \ h0 += xx; \ if (xx > big) \ big = xx; \ else if (xx < small) \ small = xx; CHECK(h1) CHECK(h2) CHECK(h3) CHECK(h4) CHECK(h5) CHECK(h6) #undef CHECK /* Compute mean of middle 5-7 values */ return UNSCALE(h0 -ALFRAC[(big + small)>>SCALEB]); } /* Mean of 3 - 5 middle values */ static gint atfilt2 (gint32 *p) { gint h0,h1,h2,h3,h4,h5,h6; /* hexagon values 2 3 */ /* 1 0 4 */ /* 6 5 */ gint big0,big1,small0,small1; /* map to scaled hexagon values */ h0 = M0[p[0]] + M1[p[3]] + M2[p[7]]; h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]]; h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]]; h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]]; h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]]; h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]]; h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]]; /* sum values and also discover the 2 largest and 2 smallest */ big0 = small0 = h0; small1 = G_MAXINT; big1 = 0; #define CHECK(xx) \ h0 += xx; \ if (xx > big1) \ { \ if (xx > big0) \ { \ big1 = big0; \ big0 = xx; \ } \ else \ big1 = xx; \ } \ if (xx < small1) \ { \ if (xx < small0) \ { \ small1 = small0; \ small0 = xx; \ } \ else \ small1 = xx; \ } CHECK(h1) CHECK(h2) CHECK(h3) CHECK(h4) CHECK(h5) CHECK(h6) #undef CHECK /* Compute mean of middle 3-5 values */ return UNSCALE(h0 -big0 -small0 -ALFRAC[(big1 + small1)>>SCALEB]); } /* * Mean of 1 - 3 middle values. * If only 1 value, then this is a median filter. */ static gint32 atfilt3(gint32 *p) { gint h0,h1,h2,h3,h4,h5,h6; /* hexagon values 2 3 */ /* 1 0 4 */ /* 6 5 */ gint big0,big1,big2,small0,small1,small2; /* map to scaled hexagon values */ h0 = M0[p[0]] + M1[p[3]] + M2[p[7]]; h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]]; h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]]; h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]]; h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]]; h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]]; h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]]; /* sum values and also discover the 3 largest and 3 smallest */ big0 = small0 = h0; small1 = small2 = G_MAXINT; big1 = big2 = 0; #define CHECK(xx) \ h0 += xx; \ if (xx > big2) \ { \ if (xx > big1) \ { \ if (xx > big0) \ { \ big2 = big1; \ big1 = big0; \ big0 = xx; \ } \ else \ { \ big2 = big1; \ big1 = xx; \ } \ } \ else \ big2 = xx; \ } \ if (xx < small2) \ { \ if (xx < small1) \ { \ if (xx < small0) \ { \ small2 = small1; \ small1 = small0; \ small0 = xx; \ } \ else \ { \ small2 = small1; \ small1 = xx; \ } \ } \ else \ small2 = xx; \ } CHECK(h1) CHECK(h2) CHECK(h3) CHECK(h4) CHECK(h5) CHECK(h6) #undef CHECK /* Compute mean of middle 1-3 values */ return UNSCALE(h0-big0-big1-small0-small1-ALFRAC[(big2+small2)>>SCALEB]); } /* Edge enhancement */ static gint atfilt4 (gint *p) { gint hav; /* map to scaled hexagon values and compute enhance value */ hav = M0[p[0]] + M1[p[3]] + M2[p[7]]; hav += H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]]; hav += V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]]; hav += V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]]; hav += H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]]; hav += V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]]; hav += V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]]; if (hav < 0) hav = 0; hav = UNSCALE(hav); if (hav > (gdouble)255) hav = (gdouble)255; return hav; } /* Optimal estimation - do smoothing in inverse proportion */ /* to the local variance. */ /* notice we use the globals noisevariance */ gint atfilt5(gint *p) { gint mean,variance,temp; gint h0,h1,h2,h3,h4,h5,h6; /* hexagon values 2 3 */ /* 1 0 4 */ /* 6 5 */ /* map to scaled hexagon values */ h0 = M0[p[0]] + M1[p[3]] + M2[p[7]]; h1 = H0[p[0]] + H1[p[2]] + H2[p[1]] + H3[p[8]]; h2 = V0[p[0]] + V1[p[3]] + V2[p[2]] + V3[p[1]]; h3 = V0[p[0]] + V1[p[3]] + V2[p[4]] + V3[p[5]]; h4 = H0[p[0]] + H1[p[4]] + H2[p[5]] + H3[p[6]]; h5 = V0[p[0]] + V1[p[7]] + V2[p[6]] + V3[p[5]]; h6 = V0[p[0]] + V1[p[7]] + V2[p[8]] + V3[p[1]]; mean = h0 + h1 + h2 + h3 + h4 + h5 + h6; /* compute scaled mean by dividing by 7 */ mean = AVEDIV[SCTOCSC(mean)]; /* compute scaled variance */ temp = (h1 - mean); variance = SQUARE[NOCSVAL + SCTOCSC(temp)]; /* and rescale to keep */ temp = (h2 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)]; /* within 32 bit limits */ temp = (h3 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)]; temp = (h4 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)]; temp = (h5 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)]; temp = (h6 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)]; /* (temp = h0 - mean) */ temp = (h0 - mean); variance += SQUARE[NOCSVAL + SCTOCSC(temp)]; if (variance != 0) /* avoid possible divide by 0 */ /* optimal estimate */ temp = mean + (variance * temp) / (variance + noisevariance); else temp = h0; if (temp < 0) temp = 0; temp = RUNSCALE(temp); if (temp > (gdouble)255) temp = (gdouble)255; return temp; } /* Triangle orientation is per geometric axes (not graphical axies) */ #define NW 0 /* North west triangle /| */ #define NE 1 /* North east triangle |\ */ #define SW 2 /* South west triangle \| */ #define SE 3 /* South east triangle |/ */ #define STH 2 #define EST 1 #define SWAPI(a,b) (t = a, a = -b, b = -t) /* compute the area of overlap of a hexagon diameter d, */ /* centered at hx,hy, with a unit square of center sx,sy. */ static gdouble hex_area (gdouble sx, gdouble sy, gdouble hx, gdouble hy, gdouble d) { gdouble hx0,hx1,hx2,hy0,hy1,hy2,hy3; gdouble sx0,sx1,sy0,sy1; /* compute square coordinates */ sx0 = sx - 0.5; sy0 = sy - 0.5; sx1 = sx + 0.5; sy1 = sy + 0.5; /* compute hexagon coordinates */ hx0 = hx - d/2.0; hx1 = hx; hx2 = hx + d/2.0; hy0 = hy - 0.5773502692 * d; /* d / sqrt(3) */ hy1 = hy - 0.2886751346 * d; /* d / sqrt(12) */ hy2 = hy + 0.2886751346 * d; /* d / sqrt(12) */ hy3 = hy + 0.5773502692 * d; /* d / sqrt(3) */ return triang_area(sx0,sy0,sx1,sy1,hx0,hy2,hx1,hy3,NW) + triang_area(sx0,sy0,sx1,sy1,hx1,hy2,hx2,hy3,NE) + rectang_area(sx0,sy0,sx1,sy1,hx0,hy1,hx2,hy2) + triang_area(sx0,sy0,sx1,sy1,hx0,hy0,hx1,hy1,SW) + triang_area(sx0,sy0,sx1,sy1,hx1,hy0,hx2,hy1,SE); } static gdouble triang_area (gdouble rx0, gdouble ry0, gdouble rx1, gdouble ry1, gdouble tx0, gdouble ty0, gdouble tx1, gdouble ty1, gint tt) { gdouble a,b,c,d; gdouble lx0,ly0,lx1,ly1; /* Convert everything to a NW triangle */ if (tt & STH) { gdouble t; SWAPI(ry0,ry1); SWAPI(ty0,ty1); } if (tt & EST) { gdouble t; SWAPI(rx0,rx1); SWAPI(tx0,tx1); } /* Compute overlapping box */ if (tx0 > rx0) rx0 = tx0; if (ty0 > ry0) ry0 = ty0; if (tx1 < rx1) rx1 = tx1; if (ty1 < ry1) ry1 = ty1; if (rx1 <= rx0 || ry1 <= ry0) return 0.0; /* Need to compute diagonal line intersection with the box */ /* First compute co-efficients to formulas x = a + by and y = c + dx */ b = (tx1 - tx0)/(ty1 - ty0); a = tx0 - b * ty0; d = (ty1 - ty0)/(tx1 - tx0); c = ty0 - d * tx0; /* compute top or right intersection */ tt = 0; ly1 = ry1; lx1 = a + b * ly1; if (lx1 <= rx0) return (rx1 - rx0) * (ry1 - ry0); else if (lx1 > rx1) { /* could be right hand side */ lx1 = rx1; ly1 = c + d * lx1; if (ly1 <= ry0) return (rx1 - rx0) * (ry1 - ry0); tt = 1; /* right hand side intersection */ } /* compute left or bottom intersection */ lx0 = rx0; ly0 = c + d * lx0; if (ly0 >= ry1) return (rx1 - rx0) * (ry1 - ry0); else if (ly0 < ry0) { /* could be right hand side */ ly0 = ry0; lx0 = a + b * ly0; if (lx0 >= rx1) return (rx1 - rx0) * (ry1 - ry0); tt |= 2; /* bottom intersection */ } if (tt == 0) { /* top and left intersection */ /* rectangle minus triangle */ return ((rx1 - rx0) * (ry1 - ry0)) - (0.5 * (lx1 - rx0) * (ry1 - ly0)); } else if (tt == 1) { /* right and left intersection */ return ((rx1 - rx0) * (ly0 - ry0)) + (0.5 * (rx1 - rx0) * (ly1 - ly0)); } else if (tt == 2) { /* top and bottom intersection */ return ((rx1 - lx1) * (ry1 - ry0)) + (0.5 * (lx1 - lx0) * (ry1 - ry0)); } else { /* tt == 3 */ /* right and bottom intersection */ /* triangle */ return (0.5 * (rx1 - lx0) * (ly1 - ry0)); } } /* Compute rectangle area */ static gdouble rectang_area (gdouble rx0, gdouble ry0, gdouble rx1, gdouble ry1, gdouble tx0, gdouble ty0, gdouble tx1, gdouble ty1) { /* Compute overlapping box */ if (tx0 > rx0) rx0 = tx0; if (ty0 > ry0) ry0 = ty0; if (tx1 < rx1) rx1 = tx1; if (ty1 < ry1) ry1 = ty1; if (rx1 <= rx0 || ry1 <= ry0) return 0.0; return (rx1 - rx0) * (ry1 - ry0); } static void nlfilter (GObject *config, PikaDrawable *drawable, PikaPreview *preview) { GeglBuffer *src_buffer; GeglBuffer *dest_buffer; const Babl *format; guchar *srcbuf, *dstbuf; guchar *lastrow, *thisrow, *nextrow, *temprow; gint x1, y1, y2; gint width, height, bpp; gint filtno, y, rowsize, exrowsize, p_update; gdouble alpha; gdouble radius; gint filter; g_object_get (config, "alpha", &alpha, "radius", &radius, "filter", &filter, NULL); if (preview) { pika_preview_get_position (preview, &x1, &y1); pika_preview_get_size (preview, &width, &height); y2 = y1 + height; } else { if (! pika_drawable_mask_intersect (drawable, &x1, &y1, &width, &height)) return; y2 = y1 + height; } if (pika_drawable_has_alpha (drawable)) format = babl_format ("R'G'B'A u8"); else format = babl_format ("R'G'B' u8"); src_buffer = pika_drawable_get_buffer (drawable); if (preview) dest_buffer = gegl_buffer_new (gegl_buffer_get_extent (src_buffer), format); else dest_buffer = pika_drawable_get_shadow_buffer (drawable); bpp = babl_format_get_bytes_per_pixel (format); rowsize = width * bpp; exrowsize = (width + 2) * bpp; p_update = width / 20 + 1; /* source buffer gives one pixel margin all around destination buffer */ srcbuf = g_new0 (guchar, exrowsize * 3); dstbuf = g_new0 (guchar, rowsize); /* pointers to second pixel in each source row */ lastrow = srcbuf + bpp; thisrow = lastrow + exrowsize; nextrow = thisrow + exrowsize; filtno = nlfiltInit (alpha, radius, filter); if (!preview) pika_progress_init (_("NL Filter")); /* first row */ gegl_buffer_get (src_buffer, GEGL_RECTANGLE (x1, y1, width, 1), 1.0, format, thisrow, GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_NONE); /* copy thisrow[0] to thisrow[-1], thisrow[width-1] to thisrow[width] */ memcpy (thisrow - bpp, thisrow, bpp); memcpy (thisrow + rowsize, thisrow + rowsize - bpp, bpp); /* copy whole thisrow to lastrow */ memcpy (lastrow - bpp, thisrow - bpp, exrowsize); for (y = y1; y < y2 - 1; y++) { if (((y % p_update) == 0) && !preview) pika_progress_update ((gdouble) y / (gdouble) height); gegl_buffer_get (src_buffer, GEGL_RECTANGLE (x1, y + 1, width, 1), 1.0, format, nextrow, GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_NONE); memcpy (nextrow - bpp, nextrow, bpp); memcpy (nextrow + rowsize, nextrow + rowsize - bpp, bpp); nlfiltRow (lastrow, thisrow, nextrow, dstbuf, width, bpp, filtno); gegl_buffer_set (dest_buffer, GEGL_RECTANGLE (x1, y, width, 1), 0, format, dstbuf, GEGL_AUTO_ROWSTRIDE); /* rotate row buffers */ temprow = lastrow; lastrow = thisrow; thisrow = nextrow; nextrow = temprow; } /* last row */ memcpy (nextrow - bpp, thisrow - bpp, exrowsize); nlfiltRow (lastrow, thisrow, nextrow, dstbuf, width, bpp, filtno); gegl_buffer_set (dest_buffer, GEGL_RECTANGLE (x1, y2 - 1, width, 1), 0, format, dstbuf, GEGL_AUTO_ROWSTRIDE); g_free (srcbuf); g_free (dstbuf); if (preview) { guchar *buf = g_new (guchar, width * height * bpp); gegl_buffer_get (dest_buffer, GEGL_RECTANGLE (x1, y1, width, height), 1.0, format, buf, GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_NONE); pika_preview_draw_buffer (PIKA_PREVIEW (preview), buf, width * bpp); g_free (buf); } g_object_unref (src_buffer); g_object_unref (dest_buffer); if (! preview) { pika_progress_update (1.0); pika_drawable_merge_shadow (drawable, TRUE); pika_drawable_update (drawable, x1, y1, width, height); pika_displays_flush (); } } static void nlfilter_preview (GtkWidget *widget, GObject *config) { PikaPreview *preview = PIKA_PREVIEW (widget); PikaDrawable *drawable = g_object_get_data (config, "drawable"); nlfilter (config, drawable, preview); } static gboolean nlfilter_dialog (PikaProcedure *procedure, GObject *config, PikaDrawable *drawable) { GtkWidget *dialog; GtkWidget *preview; GtkListStore *store; GtkWidget *frame; GtkWidget *scale; gboolean run; pika_ui_init (PLUG_IN_BINARY); dialog = pika_procedure_dialog_new (procedure, PIKA_PROCEDURE_CONFIG (config), _("NL Filter")); 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_drawable_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 (_("Alpha trimmed mean"), 0, _("Optimal estimation"), 1, _("Edge enhancement"), 2, NULL); frame = pika_procedure_dialog_get_int_radio (PIKA_PROCEDURE_DIALOG (dialog), "filter", PIKA_INT_STORE (store)); gtk_widget_set_margin_bottom (frame, 12); scale = pika_procedure_dialog_get_scale_entry (PIKA_PROCEDURE_DIALOG (dialog), "alpha", 1.0); gtk_widget_set_margin_bottom (scale, 12); scale = pika_procedure_dialog_get_scale_entry (PIKA_PROCEDURE_DIALOG (dialog), "radius", 1.0); gtk_widget_set_margin_bottom (scale, 12); pika_procedure_dialog_fill_box (PIKA_PROCEDURE_DIALOG (dialog), "nlfilter-vbox", "filter", "alpha", "radius", NULL); g_object_set_data (config, "drawable", drawable); g_signal_connect (preview, "invalidated", G_CALLBACK (nlfilter_preview), config); g_signal_connect_swapped (config, "notify", G_CALLBACK (pika_preview_invalidate), preview); pika_procedure_dialog_fill (PIKA_PROCEDURE_DIALOG (dialog), "nlfilter-vbox", NULL); gtk_widget_show (dialog); run = pika_procedure_dialog_run (PIKA_PROCEDURE_DIALOG (dialog)); gtk_widget_destroy (dialog); return run; }