/*
 * "$Id: despeckle.c,v 1.69 2004/09/15 19:56:01 neo Exp $"
 *
 *   Despeckle (adaptive median) filter for The GIMP -- an image manipulation
 *   program
 *
 *   Copyright 1997-1998 Michael Sweet (mike@easysw.com)
 *
 *   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.
 */

#include "config.h"

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

#include <gtk/gtk.h>

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

#include "libgimp/stdplugins-intl.h"
typedef enum
{
  MEDIAN,
  SMOLKA,
  VECTOR
} FilterMethod;

typedef struct
{
  gint          radius;
  FilterMethod  method;
} FilterValues;

/*
 * Constants...
 */

#define PLUG_IN_NAME     "plug_in_despeckle"
#define PLUG_IN_VERSION  "1.3.2 - 17 May 1998"
#define HELP_ID          "plug-in-despeckle"
#define SCALE_WIDTH      100
#define ENTRY_WIDTH        3
#define MAX_RADIUS        20

#define FILTER_ADAPTIVE  0x01
#define FILTER_RECURSIVE 0x02
#define FILTER_SMOLKA    0x04
#define FILTER_VDF       0x08
#define FILTER_VMF       0x16

#define despeckle_radius (despeckle_vals[0])    /* Radius of filter */
#define filter_type      (despeckle_vals[1])    /* Type of filter */
#define black_level      (despeckle_vals[2])    /* Black level */
#define white_level      (despeckle_vals[3])    /* White level */
#define update_toggle    (despeckle_vals[4])    /* Update the preview? */

#define SWAP(a,b) { guchar temp=(a);(a)=(b);(b)=temp; }
#define SORT(a,b) { if ((a)>(b)) SWAP((a),(b)); }


/*
 * Local functions...
 */

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

static void      despeckle (void);
static void      despeckle_vdf    (guchar *src, 
                                   guchar *dst, 
                                   gint    width, 
                                   gint    height, 
                                   gint    bpp, 
                                   gint    radius);
static void      despeckle_vmf    (guchar *src, 
                                   guchar *dst, 
                                   gint    width, 
                                   gint    height, 
                                   gint    bpp, 
                                   gint    radius);                                   
static void      despeckle_smolka (guchar *src, 
                                   guchar *dst, 
                                   gint    width, 
                                   gint    height, 
                                   gint    bpp, 
                                   gint    radius);

static gboolean  despeckle_dialog              (void);

static void      dialog_iscale_update      (GtkAdjustment *, gint *);
static void      dialog_adaptive_callback  (GtkWidget *, gpointer);
static void      dialog_recursive_callback (GtkWidget *, gpointer);
static void      dialog_smolka_callback    (GtkWidget *, gpointer);
static void      dialog_vdf_callback       (GtkWidget *, gpointer);
static void      my_radio_button_update    (GtkWidget *, gpointer);

static void      preview_update            (GtkWidget *preview);

/*
 * Globals...
 */

GimpPlugInInfo PLUG_IN_INFO =
{
  NULL,  /* init  */
  NULL,  /* quit  */
  query, /* query */
  run    /* run   */
};

static GtkWidget      *preview;                 /* Preview widget */

static GimpDrawable   *drawable = NULL;         /* Current image */
static gint            sel_x1,                  /* Selection bounds */
                       sel_y1,
                       sel_x2,
                       sel_y2;
static gint            sel_width,               /* Selection width */
                       sel_height;              /* Selection height */
static gint            img_bpp;                 /* Bytes-per-pixel in image */

static gint despeckle_vals[5] =
{
  3,                  /* Default value for the radius */
  FILTER_ADAPTIVE,    /* Default value for the filter type */
  7,                  /* Default value for the black level */
  248,                /* Default value for the white level */
  TRUE                /* Default value for the update toggle */
};

static FilterValues fvals =
{
  3.0,  /*  y radius  */
  MEDIAN
};


/*
 * 'main()' - Main entry - just call gimp_main()...
 */

MAIN ()

/*
 * 'query()' - Respond to a plug-in query...
 */

static void
query (void)
{
  static GimpParamDef   args[] =
  {
    { GIMP_PDB_INT32,    "run_mode", "Interactive, non-interactive" },
    { GIMP_PDB_IMAGE,    "image",    "Input image" },
    { GIMP_PDB_DRAWABLE, "drawable", "Input drawable" },
    { GIMP_PDB_INT32,    "radius",   "Filter box radius (default = 3)" },
    { GIMP_PDB_INT32,    "type",     "Filter type (0 = median, 1 = adaptive, 2 = recursive-median, 3 = recursive-adaptive)" },
    { GIMP_PDB_INT32,    "black",    "Black level (-1 to 255)" },
    { GIMP_PDB_INT32,    "white",    "White level (0 to 256)" }
  };

  gimp_install_procedure (PLUG_IN_NAME,
                          "Despeckle filter, typically used to \'despeckle\' "
                          "a photographic image.",
                          "This plug-in selectively performs a median or "
                          "adaptive box filter on an image.",
                          "Michael Sweet <mike@easysw.com>",
                          "Copyright 1997-1998 by Michael Sweet",
                          PLUG_IN_VERSION,
                          N_("Des_peckle..."),
                          "RGB*, GRAY*",
                          GIMP_PLUGIN,
                          G_N_ELEMENTS (args), 0,
                          args, NULL);

  gimp_plugin_menu_register (PLUG_IN_NAME,
                             N_("<Image>/Filters/Enhance"));
}


/*
 * 'run()' - Run the filter...
 */

static void
run (const gchar      *name,
     gint              nparams,
     const GimpParam  *param,
     gint             *nreturn_vals,
     GimpParam       **return_vals)
{
  GimpRunMode        run_mode;
  GimpPDBStatusType  status;
  GimpParam          *values;

  INIT_I18N ();

  /*
   * Initialize parameter data...
   */

  status   = GIMP_PDB_SUCCESS;
  run_mode = param[0].data.d_int32;

  values = g_new (GimpParam, 1);

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

  *nreturn_vals = 1;
  *return_vals  = values;

  /*
   * Get drawable information...
   */

  drawable = gimp_drawable_get (param[2].data.d_drawable);

  gimp_drawable_mask_bounds (drawable->drawable_id,
                             &sel_x1, &sel_y1, &sel_x2, &sel_y2);

  sel_width  = sel_x2 - sel_x1;
  sel_height = sel_y2 - sel_y1;
  img_bpp    = gimp_drawable_bpp (drawable->drawable_id);

  /*
   * See how we will run
   */

  switch (run_mode)
    {
    case GIMP_RUN_INTERACTIVE :
      /*
       * Possibly retrieve data...
       */

      gimp_get_data (PLUG_IN_NAME, &despeckle_radius);

      /*
       * Get information from the dialog...
       */
      if (gimp_drawable_is_rgb(drawable->drawable_id) ||
          gimp_drawable_is_gray(drawable->drawable_id))
       {
          if (!despeckle_dialog ())
          return;
       }
      break;

    case GIMP_RUN_NONINTERACTIVE:
      /*
       * Make sure all the arguments are present...
       */

      if (nparams < 4 || nparams > 7)
        status = GIMP_PDB_CALLING_ERROR;
      else if (nparams == 4)
        {
          despeckle_radius = param[3].data.d_int32;
          filter_type      = FILTER_ADAPTIVE;
          black_level      = 7;
          white_level      = 248;
        }
      else if (nparams == 5)
        {
          despeckle_radius = param[3].data.d_int32;
          filter_type      = param[4].data.d_int32;
          black_level      = 7;
          white_level      = 248;
        }
      else if (nparams == 6)
        {
          despeckle_radius = param[3].data.d_int32;
          filter_type      = param[4].data.d_int32;
          black_level      = param[5].data.d_int32;
          white_level      = 248;
        }
      else
        {
          despeckle_radius = param[3].data.d_int32;
          filter_type      = param[4].data.d_int32;
          black_level      = param[5].data.d_int32;
          white_level      = param[6].data.d_int32;
        }
      break;

    case GIMP_RUN_WITH_LAST_VALS:
      /*
       * Possibly retrieve data...
       */

      INIT_I18N();
      gimp_get_data (PLUG_IN_NAME, despeckle_vals);
        break;

    default:
      status = GIMP_PDB_CALLING_ERROR;
      break;
    }

  /*
   * Despeckle the image...
   */

  if (status == GIMP_PDB_SUCCESS)
    {
    	if (gimp_drawable_is_rgb(drawable->drawable_id) ||
            gimp_drawable_is_gray(drawable->drawable_id))
        {
          /*
           * Set the tile cache size...
           */

          gimp_tile_cache_ntiles (2 * (drawable->width + gimp_tile_width() - 1) /
                                  gimp_tile_width() + 1);

          /*
           * Run!
           */

          despeckle ();

          /*
           * If run mode is interactive, flush displays...
           */

          if (run_mode != GIMP_RUN_NONINTERACTIVE)
            gimp_displays_flush ();

          /*
           * Store data...
           */

          if (run_mode == GIMP_RUN_INTERACTIVE)
            gimp_set_data (PLUG_IN_NAME,
                           despeckle_vals, sizeof (despeckle_vals));
        }
      else
        status = GIMP_PDB_EXECUTION_ERROR;
    }

  /*
   * Reset the current run status...
   */

  values[0].data.d_status = status;

  /*
   * Detach from the drawable...
   */

  gimp_drawable_detach (drawable);
}


/*
 * 'despeckle()' - Despeckle an image using a median filter.
 *
 * A median filter basically collects pixel values in a region around the
 * target pixel, sorts them, and uses the median value. This code uses a
 * circular row buffer to improve performance.
 *
 * The adaptive filter is based on the median filter but analizes the histogram
 * of the region around the target pixel and adjusts the despeckle radius
 * accordingly.
 */

static void
despeckle (void)
{
  GimpPixelRgn  src_rgn,        /* Source image region */
                dst_rgn;        /* Destination image region */
  guchar      **src_rows,       /* Source pixel rows */
               *dst_row,        /* Destination pixel row */
               *src_ptr,        /* Source pixel pointer */
               *sort,           /* Pixel value sort array */
               *sort_ptr;       /* Current sort value */
  gint          sort_count,     /* Number of soft values */
                i, j, t, d,     /* Looping vars */
                x, y,           /* Current location in image */
                row,            /* Current row in src_rows */
                rowcount,       /* Number of rows loaded */
                lasty,          /* Last row loaded in src_rows */
                trow,           /* Looping var */
                startrow,       /* Starting row for loop */
                endrow,         /* Ending row for loop */
                max_row,        /* Maximum number of filled src_rows */
                size,           /* Width/height of the filter box */
                width,          /* Byte width of the image */
                xmin, xmax, tx, /* Looping vars */
                radius,         /* Current radius */
                hist0,          /* Histogram count for 0 values */
                hist255;        /* Histogram count for 255 values */

  /*
   * Let the user know what we're doing...
   */

  gimp_progress_init (_("Despeckling..."));

  /*
   * Setup for filter...
   */

  gimp_pixel_rgn_init (&src_rgn, drawable, sel_x1, sel_y1, sel_width, sel_height,
                       FALSE, FALSE);
  gimp_pixel_rgn_init (&dst_rgn, drawable, sel_x1, sel_y1, sel_width, sel_height,
                       TRUE, TRUE);
  if (fvals.method == SMOLKA || fvals.method == VECTOR )
    {
      guchar *src  = g_new (guchar, sel_width * sel_height * img_bpp);
      guchar *dst  = g_new (guchar, sel_width * sel_height * img_bpp);
      gimp_pixel_rgn_get_rect (&src_rgn, src, sel_x1, sel_y1,sel_width, sel_height);
      if (fvals.method == SMOLKA )
        {
          despeckle_smolka(src, dst, sel_width, sel_height,  img_bpp, despeckle_radius);
        }
      else if (fvals.method == VECTOR )
        {
          despeckle_vdf(src, dst, sel_width, sel_height,  img_bpp, despeckle_radius);
        }

      gimp_pixel_rgn_set_rect (&dst_rgn, dst, sel_x1, sel_y1, sel_width, sel_height);

       /*
        * Update the screen...
        */

       gimp_drawable_flush (drawable);
       gimp_drawable_merge_shadow (drawable->drawable_id, TRUE);
       gimp_drawable_update (drawable->drawable_id,
                             sel_x1, sel_y1, sel_width, sel_height);
       g_free(src);                      
       g_free(dst);
       return;
    }

  size     = despeckle_radius * 2 + 1;
  max_row  = 2 * gimp_tile_height ();
  width    = sel_width * img_bpp;

  src_rows    = g_new (guchar *, max_row);
  src_rows[0] = g_new (guchar, max_row * width);

  for (row = 1; row < max_row; row ++)
    src_rows[row] = src_rows[0] + row * width;

  dst_row = g_new (guchar, width),
  sort    = g_new (guchar, size * size);

  /*
   * Pre-load the first "size" rows for the filter...
   */

  if (sel_height < gimp_tile_height())
    rowcount = sel_height;
  else
    rowcount = gimp_tile_height ();

  gimp_pixel_rgn_get_rect (&src_rgn, src_rows[0], sel_x1, sel_y1, sel_width,
                           rowcount);

  row   = rowcount;
  lasty = sel_y1 + rowcount;

  /*
   * Despeckle...
   */

  for (y = sel_y1 ; y < sel_y2; y ++)
    {
      if ((y + despeckle_radius) >= lasty &&
          lasty < sel_y2)
        {
          /*
           * Load the next block of rows...
           */

          rowcount -= gimp_tile_height();
          if ((i = sel_y2 - lasty) > gimp_tile_height())
            i = gimp_tile_height();

          gimp_pixel_rgn_get_rect (&src_rgn, src_rows[row],
                                   sel_x1, lasty, sel_width, i);

          rowcount += i;
          lasty    += i;
          row      = (row + i) % max_row;
        }

      /*
       * Now find the median pixels and save the results...
       */

      radius = despeckle_radius;

      memcpy (dst_row, src_rows[(row + y - lasty + max_row) % max_row], width);

      if (y >= (sel_y1 + radius) && y < (sel_y2 - radius))
        {
          for (x = 0; x < width; x ++)
            {
              hist0   = 0;
              hist255 = 0;
              xmin    = x - radius * img_bpp;
              xmax    = x + (radius + 1) * img_bpp;

              if (xmin < 0)
                xmin = x % img_bpp;

              if (xmax > width)
                xmax = width;

              startrow = (row + y - lasty - radius + max_row) % max_row;
              endrow   = (row + y - lasty + radius + 1 + max_row) % max_row;

              for (sort_ptr = sort, trow = startrow;
                   trow != endrow;
                   trow = (trow + 1) % max_row)
                for (tx = xmin, src_ptr = src_rows[trow] + xmin;
                     tx < xmax;
                     tx += img_bpp, src_ptr += img_bpp)
                  {
                    if ((*sort_ptr = *src_ptr) <= black_level)
                      hist0 ++;
                    else if (*sort_ptr >= white_level)
                      hist255 ++;

                    if (*sort_ptr < white_level && *sort_ptr > black_level)
                      sort_ptr ++;
                  }

              /*
               * Shell sort the color values...
               */

              sort_count = sort_ptr - sort;

              if (sort_count > 1)
                {
                  for (d = sort_count / 2; d > 0; d = d / 2)
                    for (i = d; i < sort_count; i ++)
                      for (j = i - d, sort_ptr = sort + j;
                           j >= 0 && sort_ptr[0] > sort_ptr[d];
                           j -= d, sort_ptr -= d)
                        {
                          t           = sort_ptr[0];
                          sort_ptr[0] = sort_ptr[d];
                          sort_ptr[d] = t;
                        }

                  /*
                   * Assign the median value...
                   */

                  t = sort_count / 2;

                  if (sort_count & 1)
                    dst_row[x] = (sort[t] + sort[t + 1]) / 2;
                  else
                    dst_row[x] = sort[t];

                  /*
                   * Save the change to the source image too if the user
                   * wants the recursive method...
                   */

                  if (filter_type & FILTER_RECURSIVE)
                    src_rows[(row + y - lasty + max_row) % max_row][x] =
                      dst_row[x];
                }

              /*
               * Check the histogram and adjust the radius accordingly...
               */

              if (filter_type & FILTER_ADAPTIVE)
                {
                  if (hist0 >= radius || hist255 >= radius)
                    {
                      if (radius < despeckle_radius)
                        radius ++;
                    }
                  else if (radius > 1)
                    radius --;
                }
            }
        }

      gimp_pixel_rgn_set_row (&dst_rgn, dst_row, sel_x1, y, sel_width);

      if ((y & 15) == 0)
        gimp_progress_update ((double) (y - sel_y1) / (double) sel_height);
    }

  /*
   * OK, we're done.  Free all memory used...
   */

  g_free (src_rows[0]);
  g_free (src_rows);
  g_free (dst_row);
  g_free (sort);

  /*
   * Update the screen...
   */

  gimp_drawable_flush (drawable);
  gimp_drawable_merge_shadow (drawable->drawable_id, TRUE);
  gimp_drawable_update (drawable->drawable_id,
                        sel_x1, sel_y1, sel_width, sel_height);
}


/*
 * 'despeckle_dialog()' - Popup a dialog window for the filter box size...
 */

static gint
despeckle_dialog (void)
{
  GtkWidget *dialog;
  GtkWidget *main_vbox;
  GtkWidget *vbox;
  GtkWidget *hbox;
  GtkWidget *table;
  GtkWidget *frame;
  GtkWidget *button;
  GtkObject *adj;
  gboolean   run;

  gimp_ui_init ("despeckle", TRUE);

  dialog = gimp_dialog_new (_("Despeckle"), "despeckle",
                            NULL, 0,
                            gimp_standard_help_func, HELP_ID,

                            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, &update_toggle);
  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);

  /*
   * Filter type controls...
   */

  frame = gimp_frame_new (_("Type"));
  gtk_box_pack_start (GTK_BOX (main_vbox), frame, FALSE, FALSE, 0);
  gtk_widget_show (frame);
  
  
   hbox = gtk_hbox_new (FALSE, 12);
  gtk_box_pack_start (GTK_BOX (main_vbox), hbox, FALSE, FALSE, 0);
  gtk_widget_show (hbox);

  /*  parameter settings  */
  frame = gimp_int_radio_group_new (TRUE, _("Filter Method"),
                                    G_CALLBACK (my_radio_button_update),
                                    &fvals.method, fvals.method,
                                    _("_Median"), MEDIAN, NULL,
                                    _("_Smolka"), SMOLKA, NULL,
                                    _("_Vector Median"), VECTOR, NULL,
                                    NULL);

  gtk_box_pack_start (GTK_BOX (hbox), frame, TRUE, TRUE, 0);
  gtk_widget_show (frame);
 
 
  frame = gimp_frame_new (_("Median"));

  vbox = gtk_vbox_new (FALSE, 6);
  gtk_container_add (GTK_CONTAINER (frame), vbox);
  gtk_widget_show (vbox);

  button = gtk_check_button_new_with_mnemonic (_("_Adaptive"));
  gtk_box_pack_start (GTK_BOX (vbox), button, FALSE, FALSE, 0);
  gtk_toggle_button_set_active (GTK_TOGGLE_BUTTON (button),
                                (filter_type & FILTER_ADAPTIVE) ? TRUE : FALSE);
  gtk_widget_show (button);

  g_signal_connect (button, "toggled",
                    G_CALLBACK (dialog_adaptive_callback),
                    NULL);

  button = gtk_check_button_new_with_mnemonic (_("R_ecursive"));
  gtk_box_pack_start (GTK_BOX (vbox), button, FALSE, FALSE, 0);
  gtk_toggle_button_set_active (GTK_TOGGLE_BUTTON (button),
                                (filter_type & FILTER_RECURSIVE) ? TRUE : FALSE);
  gtk_widget_show (button);

  g_signal_connect (button, "toggled",
                    G_CALLBACK (dialog_recursive_callback),
                    NULL);

  gtk_box_pack_start (GTK_BOX (hbox), frame, FALSE, FALSE, 0);
  gtk_widget_show (frame);

  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);

  /*
   * Box size (radius) control...
   */

  adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 0,
                              _("_Radius:"), SCALE_WIDTH, ENTRY_WIDTH,
                              despeckle_radius, 1, MAX_RADIUS, 1, 5, 0,
                              TRUE, 0, 0,
                              NULL, NULL);
  g_signal_connect (adj, "value_changed",
                    G_CALLBACK (dialog_iscale_update),
                    &despeckle_radius);

  /*
   * Black level control...
   */

  adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 1,
                              _("_Black level:"), SCALE_WIDTH, ENTRY_WIDTH,
                              black_level, -1, 255, 1, 8, 0,
                              TRUE, 0, 0,
                              NULL, NULL);
  g_signal_connect (adj, "value_changed",
                    G_CALLBACK (dialog_iscale_update),
                    &black_level);

  /*
   * White level control...
   */

  adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 2,
                              _("_White level:"), SCALE_WIDTH, ENTRY_WIDTH,
                              white_level, 0, 256, 1, 8, 0,
                              TRUE, 0, 0,
                              NULL, NULL);
  g_signal_connect (adj, "value_changed",
                    G_CALLBACK (dialog_iscale_update),
                    &white_level);

  /*
   * Show it and wait for the user to do something...
   */

  gtk_widget_show (dialog);

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

  gtk_widget_destroy (dialog);

  /*
   * Return ok/cancel...
   */

  return run;
}

/*
 * 'preview_update()' - Update the preview window.
 */

static void
preview_update (GtkWidget *widget)
{
  GimpPixelRgn  src_rgn;        /* Source image region */
  guchar       *sort_ptr,       /* Current preview_sort value */
               *src_ptr,        /* Current source pixel */
               *dst_ptr;        /* Current destination pixel */
  gint          sort_count,     /* Number of soft values */
                i, j, t, d,     /* Looping vars */
                x, y,           /* Current location in image */
                x1, y1,         /* Offset of the preview from the drawable */
                size,           /* Width/height of the filter box */
                width,          /* Byte width of the image */
                xmin, xmax, tx, /* Looping vars */
                radius,         /* Current radius */
                hist0,          /* Histogram count for 0 values */
                hist255;        /* Histogram count for 255 values */
  guchar       *rgba;           /* Output image */
  GimpPreview  *preview;        /* The preview widget */

  guchar       *preview_src;    /* Source pixel rows */
  guchar       *preview_dst;    /* Destination pixel row */
  guchar       *preview_sort;   /* Pixel value sort array */

  preview = GIMP_PREVIEW (widget);
  rgba = g_new (guchar, preview->width * preview->height * img_bpp);
  /*
   * Setup for filter...
   */
  gimp_preview_get_position (preview, &x1, &y1);

  gimp_pixel_rgn_init (&src_rgn, drawable,
                       x1, y1, preview->width, preview->height,
                       FALSE, FALSE);

  /*
   * Pre-load the preview rectangle...
   */


  width = preview->width * img_bpp;

  preview_src  = g_new (guchar, width * preview->height);
  gimp_pixel_rgn_get_rect (&src_rgn, preview_src, x1, y1,
                           preview->width, preview->height);

  if (fvals.method == SMOLKA || fvals.method == VECTOR )
    {
    
      if (fvals.method == SMOLKA )
        {
          despeckle_smolka(preview_src, rgba, preview->width, preview->height,  img_bpp, despeckle_radius);
        }
      else if (fvals.method == VECTOR )
        {
          despeckle_vdf(preview_src, rgba, preview->width, preview->height,  img_bpp, despeckle_radius);
        }

      /*
       * Update the screen...
       */

      gimp_drawable_preview_draw_buffer (GIMP_DRAWABLE_PREVIEW (preview),
                                         rgba, preview->width * img_bpp);
      g_free (rgba);
      g_free (preview_src);
      return;      
    }
  /*
   * Despeckle...
   */
  size  = despeckle_radius * 2 + 1;
  preview_dst  = g_new (guchar, width);
  preview_sort = g_new (guchar, size * size);
  for (y = 0; y < preview->height; y ++)
    {
      /*
       * Now find the median pixels and save the results...
       */

      radius = despeckle_radius;

      memcpy (preview_dst, preview_src + y * width, width);

      if (y >= radius && y < (preview->height - radius))
        {
          for (x = 0, dst_ptr = preview_dst; x < width; x ++, dst_ptr ++)
            {
              hist0   = 0;
              hist255 = 0;
              xmin    = x - radius * img_bpp;
              xmax    = x + (radius + 1) * img_bpp;

              if (xmin < 0)
                xmin = x % img_bpp;

              if (xmax > width)
                xmax = width;

              for (i = -radius, sort_ptr = preview_sort,
                     src_ptr = preview_src + width * (y - radius);
                   i <= radius;
                   i ++, src_ptr += width)
                for (tx = xmin; tx < xmax; tx += img_bpp)
                  {
                    if ((*sort_ptr = src_ptr[tx]) <= black_level)
                      hist0 ++;
                    else if (*sort_ptr >= white_level)
                      hist255 ++;

                    if (*sort_ptr < white_level && *sort_ptr > black_level)
                      sort_ptr ++;
                  }

              /*
               * Shell preview_sort the color values...
               */

              sort_count = sort_ptr - preview_sort;

              if (sort_count > 1)
                {
                  for (d = sort_count / 2; d > 0; d = d / 2)
                    for (i = d; i < sort_count; i ++)
                      for (j = i - d, sort_ptr = preview_sort + j;
                           j >= 0 && sort_ptr[0] > sort_ptr[d];
                           j -= d, sort_ptr -= d)
                        {
                          t           = sort_ptr[0];
                          sort_ptr[0] = sort_ptr[d];
                          sort_ptr[d] = t;
                        }

                  /*
                   * Assign the median value...
                   */

                  t = sort_count / 2;

                  if (sort_count & 1)
                    *dst_ptr = (preview_sort[t] + preview_sort[t + 1]) / 2;
                  else
                    *dst_ptr = preview_sort[t];

                  /*
                   * Save the change to the source image too if the user
                   * wants the recursive method...
                   */

                  if (filter_type & FILTER_RECURSIVE)
                    preview_src[y * width + x] = *dst_ptr;
                }

              /*
               * Check the histogram and adjust the radius accordingly...
               */

              if (filter_type & FILTER_ADAPTIVE)
                {
                  if (hist0 >= radius || hist255 >= radius)
                    {
                      if (radius < despeckle_radius)
                        radius ++;
                    }
                  else if (radius > 1)
                    radius --;
                }
            }
        }

      /*
       * Draw this row...
       */

      memcpy (rgba + y * width, preview_dst, width);
    }
  
  /*
   * Update the screen...
   */

  gimp_drawable_preview_draw_buffer (GIMP_DRAWABLE_PREVIEW (preview),
                                     rgba, preview->width * img_bpp);
  g_free (rgba);

  g_free (preview_src);
  g_free (preview_dst);
  g_free (preview_sort);
}

static void
dialog_iscale_update (GtkAdjustment *adjustment,
                      gint          *value)
{
  *value = adjustment->value;

  gimp_preview_invalidate (GIMP_PREVIEW (preview));
}

static void
dialog_adaptive_callback (GtkWidget *widget,
                          gpointer   data)
{
  if (GTK_TOGGLE_BUTTON (widget)->active)
    filter_type |= FILTER_ADAPTIVE;
  else
    filter_type &= ~FILTER_ADAPTIVE;

  gimp_preview_invalidate (GIMP_PREVIEW (preview));
}

static void
dialog_recursive_callback (GtkWidget *widget,
                           gpointer  data)
{
  if (GTK_TOGGLE_BUTTON (widget)->active)
    filter_type |= FILTER_RECURSIVE;
  else
    filter_type &= ~FILTER_RECURSIVE;

  gimp_preview_invalidate (GIMP_PREVIEW (preview));
}

static void
my_radio_button_update (GtkWidget *widget,
			  gpointer   data)
{
  gint *toggle_val;

  if (gtk_toggle_button_get_active (GTK_TOGGLE_BUTTON (widget)))
    {
      toggle_val = (gint *) data;

      *toggle_val = (gint) gtk_object_get_user_data (GTK_OBJECT (widget));
    }

  gimp_toggle_button_sensitive_update (GTK_TOGGLE_BUTTON (widget));
  gimp_preview_invalidate (GIMP_PREVIEW (preview));
}

/*
 * Filter algorithm according to :
 * NEW FILTERING TECHNIQUE FOR THE IMPULSIVE NOISE REDUCTION IN COLOR IMAGES 
 * B. SMOLKA, A. CHYDZINSKI, K. N. PLATANIOTIS, AND A. N. VENETSANOPOULOS
 * URL: http://dx.doi.org/10.1155/S1024123X04110016
 */
static void 
despeckle_smolka(guchar *src, guchar *dst, gint width, gint height, gint bpp, gint radius)
{
  gint          pos1, pos2, x, y, i,j,g,h,r,s,w,z,delta,c;
  gdouble       sim,max,beta, *expo;  
  guchar        *pixel;

  /* create odd radius */
  radius+=radius%2;
  radius=radius-1;
  delta = radius /2;


  pixel = g_new (guchar, bpp);

  /* beta - similarity function coefficient */
  beta = 0.00504;
  expo = g_new(gdouble, 256*256);
  for (i=0 ; i < 256 ; i++ )
    {
      for (j=0 ; j < 256 ; j++ )
        {
          expo[i+j*256]=exp(-beta*abs(i-j));
        }
    }
  for ( i=0 ; i < width ;i++)
    {
       for ( j=0 ; j < height ;j++)
          {
             max=-1;
             pos2 = (i+j*width)*bpp;
             for ( c = 0 ; c < bpp ; c++)
               {
                  pixel[c]=src[pos2+c];
               }
             for (g=-delta ; g <= delta ; g++)
               {
                 for (h=-delta ; h <= delta ; h++)
                   {                                  
                     w=i+g;
                     z=j+h;
                     sim=0;
                     for (r=-delta ; r <= delta ; r++)
                       {
                         for (s=-delta ; s <= delta ; s++)
                           {
                              x=i+r;
                              y=j+s;
                              if (     (w!=x || z!=y ) 
                                    && (r!=0 || s!=0 ) 
                                    && (x>-1 && y >-1)
                                    && (w>-1 && z >-1)
                                  )
                                {                             
                                   pos1 = (x+y*width)*bpp;
                                   pos2 = (w+z*width)*bpp;
                                   /* sum precalculated simularity value*/                                     
                                   for ( c = 0 ; c < bpp ; c++)
                                     {
                                       sim+=expo[src[pos1+c]+src[pos2+c]*256];
                                     }
                                }
                           }
                       }                     
                     if ( sim>max &&  w>-1 && z >-1 )
                       {
                          max=sim;
                          pos2 = (w+z*width)*bpp;
                          for ( c = 0 ; c < bpp ; c++)
                            {
                               pixel[c]=src[pos2+c];
                            }
                       }
                   }
               }
               pos2 = (i+j*width)*bpp;
               for ( c = 0 ; c < bpp ; c++)
                 {
                    dst[pos2+c]=pixel[c];
                 }
          }
    }
  g_free(pixel);                       
  g_free(expo);              
}

static void
despeckle_vmf (guchar *src, guchar *dst, gint width, gint height, gint bpp, gint radius)
{
  gint          iv,ih,jv,jh,delta, x, y,pos,count,size;
  guchar        Ri,Rj,Gi,Gj,Bi,Bj;
  guchar        *Rvm, *Gvm, *Bvm;
  gdouble       *Dvm, *Svm;
  gdouble       normi, normj, sum, val, dvm;

  size = radius*radius;
  Dvm = g_new(gdouble, size*sizeof(gdouble));
  Svm = g_new(gdouble, size*sizeof(gdouble));
  Rvm = g_new(guchar, size);
  Gvm = g_new(guchar, size);
  Bvm = g_new(guchar, size);
  delta =  radius/2;
  for ( x = 0 ;  x < width; x++ )
    {
       for ( y = 0 ;  y < height; y++ )
         {
            pos = (x+(y*width))*bpp;
            Dvm[0] = sqrt(radius*radius*2+1);
            Rvm[0] = 255;
            Gvm[0] = 255;
            Bvm[0] = 255;  
            /* Make sure Svm is ininialized to a sufficient large value */
            Svm[0] = (4*radius*radius)+1; 
            count=0;
            for ( jh = -delta ; jh <= delta ; jh++)
	      {        
	         for ( jv = -delta ; jv <= delta ; jv++)
                   {  
                      if (     x+jh>=0 
                            && y+jv>=0
                            && x+jh<width
                            && y+jv<height
                         )
                        {
                           pos = (x+jh+((y+jv)*width))*bpp;
                           Rj = *(src+pos+0);
                           Gj = *(src+pos+1);
                           Bj = *(src+pos+2);
                           normj = sqrt(Rj*Rj+Gj*Gj+Bj*Bj);
                           sum = 0;                                                                  
                           for ( ih = -delta ; ih <= delta ; ih++)
                             {         
	                        for ( iv = -delta ; iv <= delta ; iv++)
                                  {
                                     if (     x+ih>=0 
                                           && y+iv>=0
                                           && x+ih<width
                                           && y+iv<height
                                         )
                                       {
                                          pos = (x+ih+((y+iv)*width))*bpp;
                                          Ri = *(src+pos+0);
                                          Gi = *(src+pos+1);
                                          Bi = *(src+pos+2);
                                          normi = sqrt(Ri*Ri+Gi*Gi+Bi*Bi);                                        
                                          /* Vector dot product           */               
                                          /* A . B = |A| * |B| * cos 0    */
                                          /* cos 0  =  A . B / |A| * |B|  */
                                          val = (Ri*Rj+Gi*Gj*Bi*Bj)/(normi*normj);
                                          /* Rounding error */
                                          val = ( val > 1 ) ? 1 : val;
                                          val = ( val < -1 ) ? -1 : val;
                                          sum+=fabs(acos(val));                                
                                       }
                                  }
                             }
                           if ( sum <= Svm[0] )
                             {                                           
                                dvm=sqrt(jh*jh+jv*jv);
                                count++;
                                Rvm[count] = Rj;
                                Gvm[count] = Gj;
                                Bvm[count] = Bj;
                                Dvm[count] = dvm;
                                Svm[count] = sum;
                              }
                        }
                      }
                   }
                 if ( count == 0 )
                   {
                      pos = (x+(y*width))*bpp;
                      Rvm[0] = *(src+pos+0);
                      Gvm[0] = *(src+pos+1);
                      Bvm[0] = *(src+pos+2);
                   }  
                 while ( count >= 1 && Svm[count] <= Svm[0])
                   {
                      if ( Dvm[count-1] <= Dvm[0] )
                        { 
                           Rvm[0] = Rvm[count];
                           Gvm[0] = Gvm[count];
                           Bvm[0] = Bvm[count];
                           Dvm[0] = Dvm[count];
                           Svm[0] = Svm[count];
                        }
                      count--;
                    }
                  pos = (x+(y*width))*bpp;
                  *(dst+pos+0) = Rvm[0];
                  *(dst+pos+1) = Gvm[0];
                  *(dst+pos+2) = Bvm[0];
         }             	           
    }

  g_free(Dvm);
  g_free(Svm);
  g_free(Rvm);
  g_free(Gvm);
  g_free(Bvm);
}

/*
 * 'despeckle()' - Despeckle an image using a vector median filter.
 *
 * A median filter basically collects pixel values in a region around the
 * target pixel, sorts them, and uses the median value. This code uses a
 * circular row buffer to improve performance.
 *
 * The adaptive filter is based on the median filter but analizes the histogram
 * of the region around the target pixel and adjusts the despeckle radius
 * accordingly.
 */

static void
despeckle_vdf (guchar *src, guchar *dst, gint width, gint height, gint bpp, gint radius)
{
  gint          iv,ih,jv,jh,delta, x, y,pos,count,size,c;
  gdouble       *Dvm, *Svm;
  gdouble       normi, normj, sum, val, dvm;
  guchar        *pixels, *pixeli, *pixelj, *pixel;
  
  /* create odd radius */
  radius+=radius%2;
  radius=radius-1;
  delta = radius /2;
  size = radius*radius;
  /*
   * Despeckle...
   */   
  Dvm = g_new(gdouble, size*sizeof(gdouble));
  Svm = g_new(gdouble, size*sizeof(gdouble));
  pixel = g_new(guchar, bpp);
  pixeli = g_new(guchar, bpp);
  pixelj = g_new(guchar, bpp);
  pixels  = g_new(guchar, size*bpp+bpp);
  for ( x = 0 ;  x < width; x++ )
    {
       for ( y = 0 ;  y < height; y++ )
         {
            pos = (x+(y*width))*bpp;
            Dvm[0] = sqrt(size*2+1);
            for ( c = 0 ; c < bpp ; c++ )
              {
                 pixels[c] = src[pos+c];
              }
            
            /* Make sure Svm is ininialized to a sufficient large value */
            Svm[0] = (4*size)+1; 
            count=0;
            for ( jh = -delta ; jh <= delta ; jh++)
	      {        
	         for ( jv = -delta ; jv <= delta ; jv++)
                   {  
                      if (     x+jh>=0 
                            && y+jv>=0
                            && x+jh<width
                            && y+jv<height
                         )
                        {
                           pos = (x+jh+((y+jv)*width))*bpp;
                           normj=0;                                                      
                           for ( c = 0 ; c < bpp ; c++ )
                             {
                               pixelj[c] = src[pos+c];
                               normj+= pixelj[c] * pixelj[c];
                             }
                           normj = sqrt(normj);
                           sum = 0;                                                                  
                           for ( ih = -delta ; ih <= delta ; ih++)
                             {         
	                        for ( iv = -delta ; iv <= delta ; iv++)
                                  {
                                     if (     x+ih>=0 
                                           && y+iv>=0
                                           && x+ih<width
                                           && y+iv<height
                                         )
                                       {
                                          pos = (x+ih+((y+iv)*width))*bpp;
                                          normi=0;                                                      
                                          for ( c = 0 ; c < bpp ; c++ )
                                            {                                               
                                               pixeli[c] = src[pos+c];
                                               normi+= pixeli[c] * pixeli[c];
                                            }
                                          normi = sqrt(normi);
                                          /* Vector dot product           */               
                                          /* A . B = |A| * |B| * cos 0    */
                                          /* cos 0  =  A . B / |A| * |B|  */
                                          val = 0;
                                          for ( c = 0 ; c < bpp ; c++ )
                                            {
                                               val+=pixeli[c]* pixelj[c];
                                            }
                                          val /= normi*normj;
                                          /* Rounding error */
                                          val = ( val > 1 ) ? 1 : val;
                                          val = ( val < -1 ) ? -1 : val;
                                          sum+=fabs(acos(val));                                
                                       }
                                  }
                             }
                           if ( sum <= Svm[0] )
                             {                                           
                                dvm=sqrt(jh*jh+jv*jv);
                                count++;
                                for ( c = 0 ; c < bpp ; c++ )
                                  {
                                     pos=count*bpp;
                                     pixels[pos+c] = pixelj[c];
                                  } 
                                Dvm[count] = dvm;
                                Svm[count] = sum;
                              }
                        }
                      }
                   }
                 if ( count == 0 )
                   {
                      pos = (x+(y*width))*bpp;
                      for ( c = 0 ; c < bpp ; c++ )
                        {
                          pixels[c] = src[pos+c];
                        } 
                   }  
                 while ( count >= 1 && Svm[count] <= Svm[0])
                   {
                      if ( Dvm[count-1] <= Dvm[0] )
                        { 
                          for ( c = 0 ; c < bpp ; c++ )
                            {
                              pos = count*bpp;
                              pixels[c] = pixels[pos+c];
                            }
                           Dvm[0] = Dvm[count];
                           Svm[0] = Svm[count];
                        }
                      count--;
                    }
                  pos = (x+(y*width))*bpp;
                  for ( c = 0 ; c < bpp ; c++ )
                    {
                       dst[pos+c] = pixels[c];
                    }
         }        
    }

  /*
   * OK, we're done.  Free all memory used...
   */
  g_free(Dvm);
  g_free(Svm);
  g_free(pixel);
  g_free(pixeli);
  g_free(pixelj);
  g_free(pixels);
}
