/*
 * Copyright (C) 2005 Peter Heckert
 *
 * peter /dot/ heckert /at/ arcor /dot/ de
 *
 *
 * 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 2 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, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 */


/*
 * This is the common unsharp-plugin for gimp-2.2, enhanced by me.
 * Gimp >= 2.2 is required.
 *
 *
 * INSTALL with: $ gimptool-2.2  --install unsharp2-#.##.c (#.## is the version number)
 * (gimp-devel must be installed)
 *
 *
 * It will install in your user-plugin-folder 
 * To uninstall, simply delete it from  ~/.gimp-2.2/plugins
 * Or type gimptool --help and read and learn ;-)
 *
 * It will show up as "Unsharp mask 2 V #.##" 
 * The original Unsharp mask is not overridden; it is still available. 
 *
 *===== Version History ==============================================================================
 *
 * 29 May 2005 Version 0.10
 *
 * This version had an error. I uploaded wrong file. It doesnt work. Sorry.
 *
 * 29 May 2005 Version 0.11
 *
 * This should work as advertised ;-)
 * 
 * ==== Usage ========================================================================================
 *
 * Radius is the radius as known from standard USM
 * Because here a fast iir-filter is used, very high radius is possible without slowdown.
 * This comes in handy, when USM is used to enhance local contrast.
 * I do this often, and this was very very slow with gimp's original USM.
 * ....................................................................................................
 *
 * Amount+ and Amount- are the amount factors, as known from USM, however, "+" means amount of
 * brightening and "-" means amount of darkening 
 * Threshold is the treshold as known from gimps USM
 * ....................................................................................................
 *
 * Gamma changes the gamma factor for internal processing. High gamma values can reduce sharpening halos.
 * The gamma of the resulting image is not changed.
 * (Try gamma of 5.0 and +/- amount upto 5.0. This sounds unreasonable, but it works very good )
 *
 * ....................................................................................................
 *
 * Highlight-,Light-,Midtone-,Shadow Amount 
 * 
 * This can be used, to constrain the filter amount to darker or lighter areas of the image. 
 * Especially when USM is used to enhance local contrast it is useful.
 * For normal sharpening these sliders can be left at 1.0
 *
 * Preserve currently does nothing in sharpening mode.
 * In a later version this will protect color, that means it will prevent color artifacts 
 *
 * .......... Special feature, Blur and denoise ......................................................
 *
 * When both, Amout+ and Amount- are set to zero, then USM2 goes into blur mode.
 * The image is blurred with given radius.
 * The blurred image is copied into the Output image, where the given Threshold is used as
 * criteria. As long Threshold is 0.0 no blur happens. When Threshold is larger, then we get an effect
 * similar to selective gaussian blur. When Threshold is 1.0, then we get pure gaussian blur.
 * Gamma value can improve denoising results, it works with best selectivity at gamma values <= 0.5
 * Highlight Amount .... Shadow Amount can be used to constrain the filtereffect to brightness ranges.
 *
 * Because the same treshold can be used for sharpening and for denoising, this comes in handy when
 * USM has amplified some noise in sharpening mode, or when denoised pictures should be sharpened.
 *
 * Preserve protects original luminace while denoising/ blurring.
 * This comes in handy when we want to remove color noise in shadows before or after sharpening.
 *
 * I have also written  another denoising tool (dcamnoise). However, sometimes, when I see the results
 * which I get with this modified USM version here, then I am astonished what it can do, when the parameters are
 * carefully adjusted. Especially with low frequency noise (that is already filtered by camera firmware) it
 * is often better.
 *
 */ 


#define STANDALONE


#ifndef STANDALONE
#include "config.h"
#endif


#include <stdlib.h>
#include <string.h>
#include <stdio.h>

#include <gtk/gtk.h>

#include <libgimp/gimp.h>
#include <libgimp/gimpui.h>


#ifdef STANDALONE

#define INIT_I18N voidproc
void voidproc(void){};
#define _(x) (x)
#define N_(x)  (x)

#else

#include "libgimp/stdplugins-intl.h"

#endif


#define PLUG_IN_VERSION "0.11"

#define SCALE_WIDTH   150
#define ENTRY_WIDTH     4

// This is for testing only!
#define float double
#define gfloat gdouble


/* Uncomment this line to get a rough estimate of how long the plug-in
 * takes to run.
 */

/*  #define TIMER  */


typedef struct
{
  gdouble  radius;
  gdouble  amount;
  gdouble  threshold;
  gdouble  threshold2;
  gdouble  gamma;
  gdouble  high,light,midtone,shadow; 
  gdouble  preserve;
  gboolean update_preview;
} UnsharpMaskParams;

typedef struct
{
  gboolean  run;
} UnsharpMaskInterface;

/* local function prototypes */
static void      query (void);
static void      run   (const gchar      *name,
                        gint              nparams,
                        const GimpParam  *param,
                        gint             *nreturn_vals,
                        GimpParam       **return_vals);

static void      blur_line           (const gdouble  *ctable,
                                      const gdouble  *cmatrix,
                                      gint            cmatrix_length,
                                      const guchar   *src,
                                      guchar         *dest,
                                      gint            len,
                                      glong           bytes);
static gint      gen_convolve_matrix (gdouble         std_dev,
                                      gdouble       **cmatrix);
static gdouble * gen_lookup_table    (const gdouble  *cmatrix,
                                      gint            cmatrix_length);
static void      unsharp_region      (GimpPixelRgn   *srcPTR,
                                      GimpPixelRgn   *dstPTR,
                                      gint            bytes,
                                      gdouble         radius,
                                      gdouble         amount,
                                      gint            x1,
                                      gint            x2,
                                      gint            y1,
                                      gint            y2,
                                      gboolean        show_progress);

static void      unsharp_mask        (GimpDrawable   *drawable,
                                      gdouble         radius,
                                      gdouble         amount);

static gboolean  unsharp_mask_dialog (GimpDrawable   *drawable);
static void      preview_update      (GimpPreview    *preview);


/* create a few globals, set default values */
static UnsharpMaskParams unsharp_params =
  {
    5.0, /* default radius = 5 */
    0.5, /* default amount = .5 */
    0.0,   /* default threshold = 0 */
    0.0,   /* PH: threshold2 */
    1.8, /* PH: gamma */
    1.0,1.0,1.0,1.0, /* Highlight,Light,Midtone,Shadow Noise */ 
    0.0,   
    TRUE /* default is to update the preview */
  };

/* Setting PLUG_IN_INFO */
GimpPlugInInfo PLUG_IN_INFO =
  {
    NULL,  /* init_proc  */
    NULL,  /* quit_proc  */
    query, /* query_proc */
    run,   /* run_proc   */
  };


gfloat lut_gamma, inv_gamma;  
  
gfloat lut[256];  
  
static void 
lut_init(gfloat g)
{
 int i;
 
 for (i = 1; i< 255; i++) lut[i] = pow((gfloat) i/255.0, g);
 lut[0]= 0.0;
 lut[255] = 1.0;

 lut_gamma = g;
 inv_gamma = 1.0/g;

}




MAIN ()

static void
query (void)
{
  static GimpParamDef args[] =
    {
      { GIMP_PDB_INT32,    "run_mode",  "Interactive, non-interactive" },
      { GIMP_PDB_IMAGE,    "image",     "(unused)" },
      { GIMP_PDB_DRAWABLE, "drawable",  "Drawable to draw on" },
      { GIMP_PDB_FLOAT,    "radius",    "Radius of gaussian blur (in pixels > 1.0)" },
      { GIMP_PDB_FLOAT,    "amount",    "Strength of effect" },
      { GIMP_PDB_FLOAT,    "threshold", "Threshold" },
      { GIMP_PDB_FLOAT,    "threshold2", "Threshold2" },      
      { GIMP_PDB_FLOAT,    "gamma", "Gamma" } ,      
      { GIMP_PDB_FLOAT,    "high", "Highlight Noise" },      
      { GIMP_PDB_FLOAT,    "light", "Light Noise" },      
      { GIMP_PDB_FLOAT,    "midtone", "Midtone Noise" },      
      { GIMP_PDB_FLOAT,    "shadow", "Shadow Noise" },      
      
      { GIMP_PDB_FLOAT,    "preserve", "Preserve Luminance or Color" }      
    
    
    };

  gimp_install_procedure ("plug_in_unsharp_mask2-"PLUG_IN_VERSION,
                          "An unsharp mask filter",
                          "The unsharp mask is a sharpening filter that works "
                          "by comparing using the difference of the image and "
                          "a blurred version of the image.  It is commonly "
                          "used on photographic images, and is provides a much "
                          "more pleasing result than the standard sharpen "
                          "filter.",

                          "This is an experimental version with hacked features ",
			  "",
			  "",
			    

			  N_("_Unsharp Mask2-"PLUG_IN_VERSION"..."),
                          "GRAY*, RGB*",
                          GIMP_PLUGIN,
                          G_N_ELEMENTS (args), 0,
                          args, NULL);

  gimp_plugin_menu_register ("plug_in_unsharp_mask2-"PLUG_IN_VERSION,
                             "<Image>/Filters/Enhance");
}

static void
run (const gchar      *name,
     gint              nparams,
     const GimpParam  *param,
     gint             *nreturn_vals,
     GimpParam       **return_vals)
{
  static GimpParam   values[1];
  GimpPDBStatusType  status = GIMP_PDB_SUCCESS;
  GimpDrawable      *drawable;
  GimpRunMode        run_mode;
#ifdef TIMER
  GTimer            *timer = g_timer_new ();
#endif

  run_mode = param[0].data.d_int32;

  *return_vals  = values;
  *nreturn_vals = 1;

  values[0].type          = GIMP_PDB_STATUS;
  values[0].data.d_status = status;

  INIT_I18N ();

  /*
   * Get drawable information...
   */
  drawable = gimp_drawable_get (param[2].data.d_drawable);
  gimp_tile_cache_ntiles (2 * (drawable->width / gimp_tile_width () + 1));

  switch (run_mode)
    {
    case GIMP_RUN_INTERACTIVE:
      gimp_get_data ("plug_in_unsharp_mask2-"PLUG_IN_VERSION, &unsharp_params);
      /* Reset default values show preview unmodified */

      /* initialize pixel regions and buffer */
      if (! unsharp_mask_dialog (drawable))
        return;

      break;

    case GIMP_RUN_NONINTERACTIVE:
      if (nparams != 13)
        {
          status = GIMP_PDB_CALLING_ERROR;
        }
      else
        {
          unsharp_params.radius = param[3].data.d_float;
          unsharp_params.amount = param[4].data.d_float;
          unsharp_params.threshold = param[5].data.d_float;
          unsharp_params.threshold2 = param[6].data.d_float;
          unsharp_params.gamma = param[7].data.d_float;
          
          unsharp_params.high = param[8].data.d_float; 
          unsharp_params.light = param[9].data.d_float; 
          unsharp_params.midtone = param[10].data.d_float; 
          unsharp_params.shadow = param[11].data.d_float; 
          
          unsharp_params.preserve = param[12].data.d_float; 

          
          /* make sure there are legal values */
          if ((unsharp_params.radius < 0.0) ||
              (unsharp_params.amount < 0.0))
            status = GIMP_PDB_CALLING_ERROR;
        }
      break;

    case GIMP_RUN_WITH_LAST_VALS:
      gimp_get_data ("plug_in_unsharp_mask2-"PLUG_IN_VERSION, &unsharp_params);
      break;

    default:
      break;
    }

  if (status == GIMP_PDB_SUCCESS)
    {
      drawable = gimp_drawable_get (param[2].data.d_drawable);

      /* here we go */
      unsharp_mask (drawable, unsharp_params.radius, unsharp_params.amount);

      gimp_displays_flush ();

      /* set data for next use of filter */
      gimp_set_data ("plug_in_unsharp_mask2-"PLUG_IN_VERSION, &unsharp_params,
                     sizeof (UnsharpMaskParams));

      gimp_drawable_detach(drawable);
      values[0].data.d_status = status;
    }

#ifdef TIMER
  g_printerr ("%f seconds\n", g_timer_elapsed (timer, NULL));
  g_timer_destroy (timer);
#endif
}



static struct iir_param
{
 gdouble B,b1,b2,b3,b0,r,q;
 gdouble *p;
} iir; 

static void iir_init(double r)
{
  
  iir.r = r;
  gdouble q;
  
  if ( r >= 2.5) q = 0.98711 * r - 0.96330;
  else q = 3.97156 - 4.14554 * sqrt(1.0-0.26891 * r);
  
  iir.q = q;
  iir.b0 = 1.57825 + ((0.422205 * q  + 1.4281) * q + 2.44413) *  q;
  iir.b1 = ((1.26661 * q +2.85619) * q + 2.44413) * q / iir.b0;
  iir.b2 = - ((1.26661*q +1.4281) * q * q ) / iir.b0;
  iir.b3 = 0.422205 * q * q * q / iir.b0;
  iir.B = 1.0 - (iir.b1 + iir.b2 + iir.b3);

}


static void iir_filter(gdouble *data, gint width)
/* 
 * Very fast gaussian blur with infinite impulse response filter
 * The row is blurred in forward direction and then in backward direction
 * So we achieve zero phase errors and symmetric impulse response
 * and good isotropy
 *
 * Theory for this filter can be found at:
 * <http://www.ph.tn.tudelft.nl/~lucas/publications/1995/SP95TYLV/SP95TYLV.pdf>
 * It is usable for radius downto 0.5. Lower radius must be done with the old
 * method. The old method also is very fast at low radius, so this doesnt matter 
 *
 * Double floating point precision is necessary for radius > 50, as my experiments
 * have shown. On my system (Duron, 1,2 GHz) the speed difference between double
 * and float is neglectable. 
 */

{
  gint r = ceil(iir.r/2.0);
  if (!r) return; 
  if (width < 3) return; /* We are lazy and dont filter this */
  
  gdouble lval = 0.0, rval = 0.0;
  gdouble *const lp = data, *const rp = data + width-1;
  
  gint i, sum;
    
  /* estimate *filtered* startup values for left and right edge */
  /* It's only rough estimation, we cannot really know, what's outside image bounds */
  /* Probably we could insert more exact FIR-algorithm for radius < 2.5 here and use this */
  /* to calculate startup-values also. Note that radius for a single pass = radius/sqrt(2.0) */
  
  
  for (i = 1, sum = 0; i <= r; i++)
    {
    if ( r-i < width)
       {
       sum += i;
       lval += lp[r-i] * (gfloat) i;
       }
     }
  
  lval /=   (gdouble) sum;
  
  gint start = width - 1 - r;   
    
  for (i = 1, sum=0 ; i <= r  ; i++)
    {
    if (start + i < width ) 
       {
       sum += i;
       rval += lp[start + i ] * (gfloat) i;
       }
    }
  
  rval /=  (gdouble) sum;
        
  /* Forward convolution pass */  
  data[0] = iir.b3 * lval + iir.b2*lval + iir.b1*lval + iir.B*data[0];
  data++;
  data[0] = iir.b3 * lval + iir.b2*lval + iir.b1*data[-1] + iir.B*data[0];
  data++;
  data[0] = iir.b3 * lval + iir.b2*data[-2] + iir.b1*data[-1] + iir.B*data[0];
  data++;
  
  
  {
  /* Hoping compiler will use optimal alternative, if not enough registers */
  register double d1 = data[-1], d2 = data[-2],d3 = data[-3];
  
    while (data <=  rp){
      
      *data *= iir.B;
      *data +=  iir.b3 * d3;      
      *data +=  iir.b2 * (d3 = d2);    
      *data +=  iir.b1 * (d2 = d1); 
      d1 = *data++;
    } 
  
  
  }
  
  data--;   
  /* Backward convolution pass */
  data[0] = iir.B * data[0] + iir.b1*rval + iir.b2*rval + iir.b3*rval;
  data--;
  data[0] = iir.B * data[0] + iir.b1*data[1] + iir.b2*rval + iir.b3*rval;
  data--;
  data[0] = iir.B * data[0] + iir.b1*data[1] + iir.b2*data[2] + iir.b3*rval;
  data--;
  
  {
  /* Hoping compiler will use optimal alternative, if not enough registers */
  register double d1 = data[1], d2 = data[2], d3 = data[3];
  
    while (data >=  lp){
      *data *= iir.B;
      *data +=  iir.b3 * d3;      
      *data +=  iir.b2 * (d3 = d2);    
      *data +=  iir.b1 * (d2 = d1); 
      d1 = *data--;
    } 
  
  }


}

/* This function is written as if it is blurring a column at a time,
 * even though it can operate on rows, too.  There is no difference
 * in the processing of the lines, at least to the blur_line function.
 */
static void
blur_line (const gdouble *ctable,
           const gdouble *cmatrix,
           gint           cmatrix_length,
           const guchar  *src,
           guchar        *dest,
           gint           len,   /* length of src and dest */
           glong          bytes) /* Bits per plane */  
{
  gfloat scale;
  gfloat sum;
  gint    b,i = 0;
  gint    j = 0;
  gint    row;

  for (b = 0; b<bytes; b++){
     gint idx;
     for (row =b, idx=0 ; idx < len; row +=bytes, idx++) iir.p[idx] = lut[src[row]]; 
        
     iir_filter(iir.p, len);
     
     for (row =b, idx=0; idx < len; row +=bytes, idx++){
        gfloat value = (pow( iir.p[idx],inv_gamma)*255.0+0.5);       
        dest[row] = CLAMP( (gint) value, 0, 255); 
     }
  }
  
}

static void
unsharp_mask (GimpDrawable *drawable,
              gdouble       radius,
              gdouble       amount)
{
  GimpPixelRgn srcPR, destPR;
  gint         x1, y1, x2, y2;

  /* initialize pixel regions */
  gimp_pixel_rgn_init (&srcPR, drawable,
                       0, 0, drawable->width, drawable->height, FALSE, FALSE);
  gimp_pixel_rgn_init (&destPR, drawable,
                       0, 0, drawable->width, drawable->height, TRUE, TRUE);

  /* Get the input */
  gimp_drawable_mask_bounds (drawable->drawable_id, &x1, &y1, &x2, &y2);

  unsharp_region (&srcPR, &destPR, drawable->bpp,
                  radius, amount,
                  x1, x2, y1, y2,
                  TRUE);

  gimp_drawable_flush (drawable);
  gimp_drawable_merge_shadow (drawable->drawable_id, TRUE);
  gimp_drawable_update (drawable->drawable_id, x1, y1, x2 - x1, y2 - y1);
}


static gfloat noise_factor(gfloat value)
{
 
 gfloat v;
 
 if (value <= 0.25) return unsharp_params.shadow; 
 
 if (value <= 0.5){
   v = (0.5-value)/0.25;
   return (unsharp_params.shadow * v + unsharp_params.midtone * (1.0-v)); 
 }
 
 if (value  <= 0.75){ 
   v = (0.75-value)/0.25;
   return (unsharp_params.midtone * v + unsharp_params.light * (1.0-v)); 
 }
 
 v = (1.0-value)/0.25;
 return (unsharp_params.light *v + unsharp_params.high * (1.0-v)); 

}


/* Perform an unsharp mask on the region, given a source region, dest.
 * region, width and height of the regions, and corner coordinates of
 * a subregion to act upon.  Everything outside the subregion is unaffected.
 */
static void
unsharp_region (GimpPixelRgn *srcPR,
                GimpPixelRgn *destPR,
                gint          bytes, /* Bytes per pixel */
                gdouble       radius,
                gdouble       amount,
                gint          x1,    /* Corners of subregion */
                gint          x2,
                gint          y1,
                gint          y2,
                gboolean      show_progress)
{
  guchar  *src;
  guchar  *dest;
  gint     width   = x2 - x1;
  gint     height  = y2 - y1;
  gdouble *cmatrix = NULL;
  gint     cmatrix_length;
  gdouble *ctable;
  gint     row, col;
  
    
  lut_init(unsharp_params.gamma);
  
  
  
  gfloat     amountm =  unsharp_params.threshold;
  gfloat     threshold = unsharp_params.threshold2;
  gboolean sharpen = amount != 0.0 || amountm != 0.0;
  if (!sharpen) amount = amountm = 1.0; 
  
  if (show_progress)
    gimp_progress_init (_("Blurring..."));

  /* generate convolution matrix
     and make sure it's smaller than each dimension */
  cmatrix_length = gen_convolve_matrix (radius, &cmatrix);
  iir_init(unsharp_params.radius);
  /* generate lookup table */
  ctable = gen_lookup_table (cmatrix, cmatrix_length);
  
  /* allocate buffers */
  src  = g_new (guchar, MAX (width, height) * bytes);
  dest = g_new (guchar, MAX (width, height) * bytes);
  iir.p = g_new(gdouble, MAX (width, height));
  
  /* blur the rows */
  for (row = 0; row < height; row++)
    {
      gimp_pixel_rgn_get_row (srcPR, src, x1, y1 + row, width);
      blur_line (ctable, cmatrix, cmatrix_length, src, dest, width, bytes);
      gimp_pixel_rgn_set_row (destPR, dest, x1, y1 + row, width);

      if (show_progress && row % 8 == 0)
        gimp_progress_update ((gdouble) row / (3 * height));
    }

  /* blur the cols */
  for (col = 0; col < width; col++)
    {
      gimp_pixel_rgn_get_col (destPR, src, x1 + col, y1, height);
      blur_line (ctable, cmatrix, cmatrix_length, src, dest, height, bytes);
      gimp_pixel_rgn_set_col (destPR, dest, x1 + col, y1, height);

      if (show_progress && col % 8 == 0)
        gimp_progress_update ((gdouble) col / (3 * width) + 0.33);
    }

  if (show_progress)
    gimp_progress_init (_("Merging..."));

  /* merge the source and destination (which currently contains
     the blurred version) images */
  for (row = 0; row < height; row++){
      const guchar *s = src;
      guchar       *d = dest;
      gfloat        value;
      gint          u, v;

      /* get source row */
      gimp_pixel_rgn_get_row (srcPR, src, x1, y1 + row, width);

      /* get dest row */
      gimp_pixel_rgn_get_row (destPR, dest, x1, y1 + row, width);

      /* combine the two */
      
      //gfloat t2=threshold2;
           
      for (u = 0; u < width; u++){
          
          gfloat lum,red,green,blue;
          
          red = (gfloat) s[0]/255.0;
          if (bytes > 2) 
            {
            green = (gfloat) s[1]/255.0;
            blue =  (gfloat) s[2]/255.0;
            }
          else  green = blue = red;

//define FR 0.3
//define FG 0.59
//#define FB 0.11
          
#define FR 0.212671
#define FG 0.715160
#define FB 0.072169
                            
          lum = (FR*red + FG*green + FB*blue);
          //gfloat a = noise_factor(lum);
  
          
          for (v = 0; v < bytes; v++){
       
              gfloat value = (float) s[v]/255.0;
              gfloat diff = value - (float) d[v]/255.0;
       
            
            gfloat a = noise_factor((gfloat) d[v]/255.0);
            
            if (sharpen){
                
                if ( diff > threshold) diff = diff - threshold; 
                else if ( diff < -threshold) diff = diff +threshold;
                else diff = 0.0; 
                
               //diff *= 1.0 - exp(-diff*diff/(2.0*threshold*threshold));
            } else { /* blur the image */
                 
                    
                diff = - diff;
                
                if (diff > threshold) diff = threshold;
                else if (diff < -threshold) diff = -threshold;
                
               // if (threshold < 1.0)  diff *= -exp(-diff*diff/(2.0*threshold*threshold));
               // else diff *= -1;
            } 
            if (diff >= 0.0) value += diff * a *amount;
            else value += diff * a *amountm;
            value = value *255.0+0.5;
//            value = pow(value,inv_gamma)*255.0+0.5;
            d[v] = CLAMP(value,0,255);
            
        } /* next color[v] */
                       
        if (!sharpen){    
          gfloat lum2,red2,green2,blue2; 
           
          red2 = (gfloat) d[0]/255.0;
          if (bytes > 2) 
              {
              green2 = (gfloat) d[1]/255.0;
              blue2 =  (gfloat) d[2]/255.0;
              }
          else  green2 = blue2 = red2;
                  
          lum2 = (FR*red2 + FG*green2 + FB*blue2);
           // gfloat dl = (0.01+lum)/(0.01+lum2);
          gfloat dl = lum - lum2;            
  
          gfloat p = unsharp_params.preserve;      
          gfloat pm = 1.0 - p;
              
          red2 = p*(red2+dl) + pm * red2;
          green2 = p*(green2+dl) + pm * green2;
          blue2 = p*(blue2+dl) + pm * blue2;

              
          red2 = red2 *255.0 + 0.5;         
          green2 = green2 *255.0 + 0.5;         
          blue2 = blue2 *255.0 + 0.5;         
              
          d[0] = CLAMP(red2,0,255);
          if (bytes > 2) d[1] = CLAMP(green2,0,255), d[2] = CLAMP(blue2,0,255);    
        }     
        d += bytes;
        s +=bytes;
   } /* next pixel[u] */

   if (show_progress && row % 8 == 0)
        gimp_progress_update ((gdouble) row / (3 * height) + 0.67);

   gimp_pixel_rgn_set_row (destPR, dest, x1, y1 + row, width);
 } /* next row */

  if (show_progress)
    gimp_progress_update (0.0);
  g_free (iir.p);
  g_free (dest);
  g_free (src);
  g_free (ctable);
  g_free (cmatrix);
}

/* generates a 1-D convolution matrix to be used for each pass of
 * a two-pass gaussian blur.  Returns the length of the matrix.
 */
static gint
gen_convolve_matrix (gdouble   radius,
                     gdouble **cmatrix_p)
{
  gdouble *cmatrix;
  gdouble  std_dev;
  gdouble  sum;
  gint     matrix_length;
  gint     i, j;

  /* we want to generate a matrix that goes out a certain radius
   * from the center, so we have to go out ceil(rad-0.5) pixels,
   * inlcuding the center pixel.  Of course, that's only in one direction,
   * so we have to go the same amount in the other direction, but not count
   * the center pixel again.  So we double the previous result and subtract
   * one.
   * The radius parameter that is passed to this function is used as
   * the standard deviation, and the radius of effect is the
   * standard deviation * 2.  It's a little confusing.
   */
  radius = fabs (radius) + 1.0;

  std_dev = radius - 1.0;
  radius = std_dev * 2;

  /* go out 'radius' in each direction */
  matrix_length = 2 * ceil (radius - 0.5) + 1;
  if (matrix_length <= 0)
    matrix_length = 1;

  *cmatrix_p = g_new (gdouble, matrix_length);
  cmatrix = *cmatrix_p;

  /*  Now we fill the matrix by doing a numeric integration approximation
   * from -2*std_dev to 2*std_dev, sampling 50 points per pixel.
   * We do the bottom half, mirror it to the top half, then compute the
   * center point.  Otherwise asymmetric quantization errors will occur.
   *  The formula to integrate is e^-(x^2/2s^2).
   */

  /* first we do the top (right) half of matrix */
  for (i = matrix_length / 2 + 1; i < matrix_length; i++)
    {
      gdouble base_x = i - (matrix_length / 2) - 0.5;

      sum = 0;
      for (j = 1; j <= 50; j++)
        {
          if (base_x + 0.02 * j <= radius)
            sum += exp (- SQR (base_x + 0.02 * j) / (2 * SQR (std_dev)));
        }

      cmatrix[i] = sum / 50;
    }

  /* mirror the thing to the bottom half */
  for (i = 0; i <= matrix_length / 2; i++)
    cmatrix[i] = cmatrix[matrix_length - 1 - i];

  /* find center val -- calculate an odd number of quanta to make it symmetric,
   * even if the center point is weighted slightly higher than others. */
  sum = 0;
  for (j = 0; j <= 50; j++)
    sum += exp (- SQR (0.5 + 0.02 * j) / (2 * SQR (std_dev)));

  cmatrix[matrix_length / 2] = sum / 51;

  /* normalize the distribution by scaling the total sum to one */
  sum = 0;
  for (i = 0; i < matrix_length; i++)
    sum += cmatrix[i];

  for (i = 0; i < matrix_length; i++)
    cmatrix[i] = cmatrix[i] / sum;

  return matrix_length;
}

/* ----------------------- gen_lookup_table ----------------------- */
/* generates a lookup table for every possible product of 0-255 and
   each value in the convolution matrix.  The returned array is
   indexed first by matrix position, then by input multiplicand (?)
   value.

This is not used anymore.
At large radius the LUT gets larger than processor cache, and will
slowdown overall operation.
Also modern processors do Multiplication very fast and therefore there
is no noticable win at small radius.

*/
static gdouble *
gen_lookup_table (const gdouble *cmatrix,
                  gint           cmatrix_length)
{
  gdouble       *lookup_table   = g_new (gdouble, cmatrix_length * 256);
  gdouble       *lookup_table_p = lookup_table;
  const gdouble *cmatrix_p      = cmatrix;
  gint           i, j;

  for (i = 0; i < cmatrix_length; i++)
    {
      for (j = 0; j < 256; j++)
        *(lookup_table_p++) = *cmatrix_p * (gdouble)j;

      cmatrix_p++;
    }

  return lookup_table;
}

static gboolean
unsharp_mask_dialog (GimpDrawable *drawable)
{
  GtkWidget *dialog;
  GtkWidget *main_vbox;
  GtkWidget *preview;
  GtkWidget *table;
  GtkObject *adj;
  gboolean   run;

  gimp_ui_init ("unsharp-"PLUG_IN_VERSION, TRUE);

  dialog = gimp_dialog_new (_("Unsharp Mask 2 V "PLUG_IN_VERSION), "unsharp2-"PLUG_IN_VERSION,
                            NULL, 0,
                            gimp_standard_help_func, "plug-in-unsharp-mask2-"PLUG_IN_VERSION,

                            GTK_STOCK_CANCEL, GTK_RESPONSE_CANCEL,
                            GTK_STOCK_OK,     GTK_RESPONSE_OK,

                            NULL);

  main_vbox = gtk_vbox_new (FALSE, 12);
  gtk_container_set_border_width (GTK_CONTAINER (main_vbox), 12);
  gtk_container_add (GTK_CONTAINER (GTK_DIALOG (dialog)->vbox), main_vbox);
  gtk_widget_show (main_vbox);

  preview = gimp_drawable_preview_new (drawable,
                                       &unsharp_params.update_preview);
  gtk_box_pack_start (GTK_BOX (main_vbox), preview, TRUE, TRUE, 0);
  gtk_widget_show (preview);

  g_signal_connect (preview, "invalidated",
                    G_CALLBACK (preview_update),
                    NULL);

  table = gtk_table_new (3, 3, FALSE);
  gtk_table_set_col_spacings (GTK_TABLE (table), 6);
  gtk_table_set_row_spacings (GTK_TABLE (table), 6);
  gtk_box_pack_start (GTK_BOX (main_vbox), table, FALSE, FALSE, 0);
  gtk_widget_show (table);

  adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 0,
                              _("_Radius:"), SCALE_WIDTH, ENTRY_WIDTH,
                              unsharp_params.radius, 0.1, 200.0, 0.1, 1.0, 1,
                              TRUE, 0, 0,
                              NULL, NULL);
  gimp_scale_entry_set_logarithmic(adj,TRUE);

  g_signal_connect (adj, "value_changed",
                    G_CALLBACK (gimp_double_adjustment_update),
                    &unsharp_params.radius);
  g_signal_connect_swapped (adj, "value_changed",
                            G_CALLBACK (gimp_preview_invalidate),
                            preview);

  adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 1,
                              _("_Amount +:"), SCALE_WIDTH, ENTRY_WIDTH,
                              unsharp_params.amount, 0.0, 5.0, 0.01, 0.1, 2,
                              TRUE, 0, 0,
                              NULL, NULL);

  g_signal_connect (adj, "value_changed",
                    G_CALLBACK (gimp_double_adjustment_update),
                    &unsharp_params.amount);
  g_signal_connect_swapped (adj, "value_changed",
                            G_CALLBACK (gimp_preview_invalidate),
                            preview);

  adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 2,
                              _("_Amount -:"), SCALE_WIDTH, ENTRY_WIDTH,
                              unsharp_params.threshold,
                              0.0, 5.0, 0.01, 0.1, 2,
                              TRUE, 0, 0,
                              NULL, NULL);

  //gimp_scale_entry_set_logarithmic(adj,TRUE);
  
  g_signal_connect (adj, "value_changed",
                    G_CALLBACK (gimp_double_adjustment_update),
                    &unsharp_params.threshold);
  g_signal_connect_swapped (adj, "value_changed",
                            G_CALLBACK (gimp_preview_invalidate),
                            preview);


  adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 3,
                              "_Threshold :", SCALE_WIDTH, ENTRY_WIDTH,
                              unsharp_params.threshold2,
                              0.0, 1.0, 0.01, 0.01, 2,
                              TRUE, 0, 0,
                              NULL, NULL);
  gimp_scale_entry_set_logarithmic(adj,TRUE);

  g_signal_connect (adj, "value_changed",
                    G_CALLBACK (gimp_double_adjustment_update),
                    &unsharp_params.threshold2);
  g_signal_connect_swapped (adj, "value_changed",
                            G_CALLBACK (gimp_preview_invalidate),
                            preview);



  /* PH: Added Gamma interface */

  adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 4,
                              "_Gamma:", SCALE_WIDTH, ENTRY_WIDTH,
                              unsharp_params.gamma,
                              0.01, 5.0, 0.01, 0.1, 2,
                              TRUE, 0, 0,
                              NULL, NULL);

  g_signal_connect (adj, "value_changed",
                    G_CALLBACK (gimp_double_adjustment_update),
                    &unsharp_params.gamma);
  g_signal_connect_swapped (adj, "value_changed",
                            G_CALLBACK (gimp_preview_invalidate),
                            preview);


/*

  adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 5,
                              "_Threshold - Gamma:", SCALE_WIDTH, ENTRY_WIDTH,
                              unsharp_params.dgamma,
                              -5.0, 5.0, 0.01, 0.1, 2,
                              TRUE, 0, 0,
                              NULL, NULL);

  g_signal_connect (adj, "value_changed",
                    G_CALLBACK (gimp_double_adjustment_update),
                    &unsharp_params.dgamma);
  g_signal_connect_swapped (adj, "value_changed",
                            G_CALLBACK (gimp_preview_invalidate),
                            preview);
*/
//-----------------------------------------------------------------------------  
  
  adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 5,
                              "_Highlight amount:", SCALE_WIDTH, ENTRY_WIDTH,
                              unsharp_params.high,
                              0.0, 1.0, 0.01, 0.1, 2,
                              TRUE, 0, 0,
                              NULL, NULL);

  g_signal_connect (adj, "value_changed",
                    G_CALLBACK (gimp_double_adjustment_update),
                    &unsharp_params.high);
  g_signal_connect_swapped (adj, "value_changed",
                            G_CALLBACK (gimp_preview_invalidate),
                            preview);

  adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 6,
                              "Light amount:", SCALE_WIDTH, ENTRY_WIDTH,
                              unsharp_params.light,
                              0.0, 1.0, 0.01, 0.1, 2,
                              TRUE, 0, 0,
                              NULL, NULL);

  g_signal_connect (adj, "value_changed",
                    G_CALLBACK (gimp_double_adjustment_update),
                    &unsharp_params.light);
  g_signal_connect_swapped (adj, "value_changed",
                            G_CALLBACK (gimp_preview_invalidate),
                            preview);

  adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 7,
                              "_Midtone amount:", SCALE_WIDTH, ENTRY_WIDTH,
                              unsharp_params.midtone,
                              0.0, 1.0, 0.01, 0.1, 2,
                              TRUE, 0, 0,
                              NULL, NULL);

  g_signal_connect (adj, "value_changed",
                    G_CALLBACK (gimp_double_adjustment_update),
                    &unsharp_params.midtone);
  g_signal_connect_swapped (adj, "value_changed",
                            G_CALLBACK (gimp_preview_invalidate),
                            preview);

  
  adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 8,
                              "_Shadow amount:", SCALE_WIDTH, ENTRY_WIDTH,
                              unsharp_params.shadow,
                              0.0, 1.0, 0.01, 0.1, 2,
                              TRUE, 0, 0,
                              NULL, NULL);

  g_signal_connect (adj, "value_changed",
                    G_CALLBACK (gimp_double_adjustment_update),
                    &unsharp_params.shadow);
  g_signal_connect_swapped (adj, "value_changed",
                            G_CALLBACK (gimp_preview_invalidate),
                            preview);

  
  
  adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 9,
                              "_Preserve:", SCALE_WIDTH, ENTRY_WIDTH,
                              unsharp_params.preserve,
                              0.0, 1.0, 0.01, 0.1, 2,
                              TRUE, 0, 0,
                              NULL, NULL);

  g_signal_connect (adj, "value_changed",
                    G_CALLBACK (gimp_double_adjustment_update),
                    &unsharp_params.preserve);
  g_signal_connect_swapped (adj, "value_changed",
                            G_CALLBACK (gimp_preview_invalidate),
                            preview);



  gtk_widget_show (dialog);

  run = (gimp_dialog_run (GIMP_DIALOG (dialog)) == GTK_RESPONSE_OK);

  gtk_widget_destroy (dialog);

  return run;
}

static void
preview_update (GimpPreview *preview)
{
  GimpDrawable *drawable;
  gint          x1, x2;
  gint          y1, y2;
  gint          x, y;
  gint          width, height;
  gint          border;
  GimpPixelRgn  srcPR;
  GimpPixelRgn  destPR;

  drawable =
    gimp_drawable_preview_get_drawable (GIMP_DRAWABLE_PREVIEW (preview));

  gimp_pixel_rgn_init (&srcPR, drawable,
                       0, 0, drawable->width, drawable->height, FALSE, FALSE);
  gimp_pixel_rgn_init (&destPR, drawable,
                       0, 0, drawable->width, drawable->height, TRUE, TRUE);

  gimp_preview_get_position (preview, &x, &y);
  gimp_preview_get_size (preview, &width, &height);

  /* enlarge the region to avoid artefacts at the edges of the preview */
  border = 2.0 * unsharp_params.radius + 0.5;
  
  if (border > width/2) border = width/2; //Speed up preview 
  x1 = MAX (0, x - border);
  y1 = MAX (0, y - border);
  x2 = MIN (x + width  + border, drawable->width);
  y2 = MIN (y + height + border, drawable->height);

  unsharp_region (&srcPR, &destPR, drawable->bpp,
                  unsharp_params.radius, unsharp_params.amount,
                  x1, x2, y1, y2,
                  FALSE);

  gimp_pixel_rgn_init (&destPR, drawable, x, y, width, height, FALSE, TRUE);
  gimp_drawable_preview_draw_region (GIMP_DRAWABLE_PREVIEW (preview), &destPR);
}

