Imported Upstream version 1.0.25upstream/1.0.25

1 files changed, 1205 insertions, 0 deletions

diff --git a/sanei/sanei_ir.c b/sanei/sanei_ir.c
new file mode 100644
index 0000000..42e82ba
--- /dev/null
+++ b/sanei/sanei_ir.c
@@ -0,0 +1,1205 @@
+/** @file sanei_ir.c
+ *
+ * sanei_ir - functions for utilizing the infrared plane
+ *
+ * Copyright (C) 2012 Michael Rickmann <mrickma@gwdg.de>
+ *
+ * This file is part of the SANE package.
+ *
+ * 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.
+ *
+ * The threshold yen, otsu and max_entropy routines have been
+ * adapted from the FOURIER 0.8 library by M. Emre Celebi,
+ * http://sourceforge.net/projects/fourier-ipal/ which is
+ * licensed under the GNU General Public License version 2 or later.
+*/
+
+#include <stdlib.h>
+#include <string.h>
+#include <values.h>
+#include <math.h>
+
+#define BACKEND_NAME sanei_ir	/* name of this module for debugging */
+
+#include "../include/sane/sane.h"
+#include "../include/sane/sanei_debug.h"
+#include "../include/sane/sanei_ir.h"
+#include "../include/sane/sanei_magic.h"
+
+
+double *
+sanei_ir_create_norm_histo (const SANE_Parameters * params, const SANE_Uint *img_data);
+double * sanei_ir_accumulate_norm_histo (double * histo_data);
+
+
+/* Initialize sanei_ir
+ */
+void
+sanei_ir_init (void)
+{
+  DBG_INIT ();
+}
+
+
+/* Create a normalized histogram of a grayscale image, internal
+ */
+double *
+sanei_ir_create_norm_histo (const SANE_Parameters * params,
+                       const SANE_Uint *img_data)
+{
+  int is, i;
+  int num_pixels;
+  int *histo_data;
+  double *histo;
+  double term;
+
+  DBG (10, "sanei_ir_create_norm_histo\n");
+
+  if ((params->format != SANE_FRAME_GRAY)
+      && (params->format != SANE_FRAME_RED)
+      && (params->format != SANE_FRAME_GREEN)
+      && (params->format != SANE_FRAME_BLUE))
+    {
+      DBG (5, "sanei_ir_create_norm_histo: invalid format\n");
+      return NULL;
+    }
+
+  /* Allocate storage for the histogram */
+  histo_data = calloc (HISTOGRAM_SIZE, sizeof (int));
+  histo = malloc (HISTOGRAM_SIZE * sizeof (double));
+  if ((histo == NULL) || (histo_data == NULL))
+    {
+      DBG (5, "sanei_ir_create_norm_histo: no buffers\n");
+      if (histo) free (histo);
+      if (histo_data) free (histo_data);
+      return NULL;
+    }
+
+  num_pixels = params->pixels_per_line * params->lines;
+
+  DBG (1, "sanei_ir_create_norm_histo: %d pixels_per_line, %d lines => %d num_pixels\n", params->pixels_per_line, params->lines, num_pixels);
+  DBG (1, "sanei_ir_create_norm_histo: histo_data[] with %d x %ld bytes\n", HISTOGRAM_SIZE, sizeof(int));
+  /* Populate the histogram */
+  is = 16 - HISTOGRAM_SHIFT; /* Number of data bits to ignore */
+  DBG (1, "sanei_ir_create_norm_histo: depth %d, HISTOGRAM_SHIFT %d => ignore %d bits\n", params->depth, HISTOGRAM_SHIFT, is);
+  for (i = num_pixels; i > 0; i--) {
+      histo_data[*img_data++ >> is]++;
+  }
+
+  /* Calculate the normalized histogram */
+  term = 1.0 / (double) num_pixels;
+  for (i = 0; i < HISTOGRAM_SIZE; i++)
+    histo[i] = term * (double) histo_data[i];
+
+  free (histo_data);
+  return histo;
+}
+
+
+/* Create the normalized histogram of a grayscale image
+ */
+SANE_Status
+sanei_ir_create_norm_histogram (const SANE_Parameters * params,
+                                const SANE_Uint *img_data,
+                                double ** histogram)
+{
+  double *histo;
+
+  DBG (10, "sanei_ir_create_norm_histogram\n");
+
+  histo = sanei_ir_create_norm_histo (params, img_data);
+  if (!histo)
+    return SANE_STATUS_NO_MEM;
+
+  *histogram = histo;
+  return SANE_STATUS_GOOD;
+}
+
+/* Accumulate a normalized histogram, internal
+ */
+double *
+sanei_ir_accumulate_norm_histo (double * histo_data)
+{
+  int i;
+  double *accum_data;
+
+  accum_data = malloc (HISTOGRAM_SIZE * sizeof (double));
+  if (accum_data == NULL)
+    {
+      DBG (5, "sanei_ir_accumulate_norm_histo: Insufficient memory !\n");
+      return NULL;
+    }
+
+  accum_data[0] = histo_data[0];
+  for (i = 1; i < HISTOGRAM_SIZE; i++)
+    accum_data[i] = accum_data[i - 1] + histo_data[i];
+
+  return accum_data;
+}
+
+/* Implements Yen's thresholding method
+ */
+SANE_Status
+sanei_ir_threshold_yen (const SANE_Parameters * params,
+                         double * norm_histo, int *thresh)
+{
+  double *P1 = NULL;            /* cumulative normalized histogram */
+  double *P1_sq = NULL;         /* cumulative normalized histogram */
+  double *P2_sq = NULL;
+  double crit, max_crit;
+  int threshold, i;
+  SANE_Status ret = SANE_STATUS_NO_MEM;
+
+  DBG (10, "sanei_ir_threshold_yen\n");
+
+  P1 = sanei_ir_accumulate_norm_histo (norm_histo);
+  P1_sq = malloc (HISTOGRAM_SIZE * sizeof (double));
+  P2_sq = malloc (HISTOGRAM_SIZE * sizeof (double));
+  if (!P1 || !P1_sq || !P2_sq)
+    {
+      DBG (5, "sanei_ir_threshold_yen: no buffers\n");
+      goto cleanup;
+    }
+
+  /* calculate cumulative squares */
+  P1_sq[0] = norm_histo[0] * norm_histo[0];
+  for (i = 1; i < HISTOGRAM_SIZE; i++)
+      P1_sq[i] = P1_sq[i - 1] + norm_histo[i] * norm_histo[i];
+  P2_sq[HISTOGRAM_SIZE - 1] = 0.0;
+  for (i = HISTOGRAM_SIZE - 2; i >= 0; i--)
+      P2_sq[i] = P2_sq[i + 1] + norm_histo[i + 1] * norm_histo[i + 1];
+
+  /* Find the threshold that maximizes the criterion */
+  threshold = INT_MIN;
+  max_crit = DBL_MIN;
+  for (i = 0; i < HISTOGRAM_SIZE; i++)
+    {
+      crit =
+        -1.0 * SAFE_LOG (P1_sq[i] * P2_sq[i]) +
+        2 * SAFE_LOG (P1[i] * (1.0 - P1[i]));
+      if (crit > max_crit)
+        {
+          max_crit = crit;
+          threshold = i;
+        }
+    }
+
+  if (threshold == INT_MIN)
+    {
+      DBG (5, "sanei_ir_threshold_yen: no threshold found\n");
+      ret = SANE_STATUS_INVAL;
+    }
+  else
+    {
+      ret = SANE_STATUS_GOOD;
+      if (params->depth > 8)
+        {
+          i = 1 << (params->depth - HISTOGRAM_SHIFT);
+          *thresh = threshold * i + i / 2;
+        }
+      else
+        *thresh = threshold;
+      DBG (10, "sanei_ir_threshold_yen: threshold %d\n", *thresh);
+    }
+
+  cleanup:
+    if (P1)
+      free (P1);
+    if (P1_sq)
+      free (P1_sq);
+    if (P2_sq)
+      free (P2_sq);
+    return ret;
+}
+
+
+/* Implements Otsu's thresholding method
+ */
+SANE_Status
+sanei_ir_threshold_otsu (const SANE_Parameters * params,
+                          double * norm_histo, int *thresh)
+{
+  double *cnh = NULL;
+  double *mean = NULL;
+  double total_mean;
+  double bcv, max_bcv;
+  int first_bin, last_bin;
+  int threshold, i;
+  SANE_Status ret = SANE_STATUS_NO_MEM;
+
+  DBG (10, "sanei_ir_threshold_otsu\n");
+
+  cnh = sanei_ir_accumulate_norm_histo (norm_histo);
+  mean = malloc (HISTOGRAM_SIZE * sizeof (double));
+  if (!cnh || !mean)
+    {
+      DBG (5, "sanei_ir_threshold_otsu: no buffers\n");
+      goto cleanup;
+    }
+
+  mean[0] = 0.0;
+  for (i = 1; i < HISTOGRAM_SIZE; i++)
+      mean[i] = mean[i - 1] + i * norm_histo[i];
+  total_mean = mean[HISTOGRAM_SIZE - 1];
+
+  first_bin = 0;
+  for (i = 0; i < HISTOGRAM_SIZE; i++)
+    if (cnh[i] != 0)
+      {
+        first_bin = i;
+        break;
+      }
+  last_bin = HISTOGRAM_SIZE - 1;
+  for (i = HISTOGRAM_SIZE - 1; i >= first_bin; i--)
+    if (1.0 - cnh[i] != 0)
+      {
+        last_bin = i;
+        break;
+      }
+
+  threshold = INT_MIN;
+  max_bcv = 0.0;
+  for (i = first_bin; i <= last_bin; i++)
+    {
+      bcv = total_mean * cnh[i] - mean[i];
+      bcv *= bcv / (cnh[i] * (1.0 - cnh[i]));
+      if (max_bcv < bcv)
+        {
+          max_bcv = bcv;
+          threshold = i;
+        }
+    }
+
+  if (threshold == INT_MIN)
+    {
+      DBG (5, "sanei_ir_threshold_otsu: no threshold found\n");
+      ret = SANE_STATUS_INVAL;
+    }
+  else
+    {
+      ret = SANE_STATUS_GOOD;
+      if (params->depth > 8)
+        {
+          i = 1 << (params->depth - HISTOGRAM_SHIFT);
+          *thresh = threshold * i + i / 2;
+        }
+      else
+        *thresh = threshold;
+      DBG (10, "sanei_ir_threshold_otsu: threshold %d\n", *thresh);
+    }
+  cleanup:
+    if (cnh)
+      free (cnh);
+    if (mean)
+      free (mean);
+    return ret;
+}
+
+
+/* Implements a Maximum Entropy thresholding method
+ */
+SANE_Status
+sanei_ir_threshold_maxentropy (const SANE_Parameters * params,
+                               double * norm_histo, int *thresh)
+{
+ int ih, it;
+ int threshold;
+ int first_bin;
+ int last_bin;
+ double tot_ent, max_ent;       /* entropies */
+ double ent_back, ent_obj;
+ double *P1;                    /* cumulative normalized histogram */
+ double *P2;
+ SANE_Status ret = SANE_STATUS_NO_MEM;
+
+ DBG (10, "sanei_ir_threshold_maxentropy\n");
+
+ /* Calculate the cumulative normalized histogram */
+ P1 = sanei_ir_accumulate_norm_histo (norm_histo);
+ P2 = malloc (HISTOGRAM_SIZE * sizeof (double));
+ if (!P1 || !P2)
+   {
+     DBG (5, "sanei_ir_threshold_maxentropy: no buffers\n");
+     goto cleanup;
+   }
+
+ for ( ih = 0; ih < HISTOGRAM_SIZE; ih++ )
+   P2[ih] = 1.0 - P1[ih];
+
+ first_bin = 0;
+ for ( ih = 0; ih < HISTOGRAM_SIZE; ih++ )
+   if (P1[ih] != 0)
+    {
+     first_bin = ih;
+     break;
+    }
+ last_bin = HISTOGRAM_SIZE - 1;
+ for ( ih = HISTOGRAM_SIZE - 1; ih >= first_bin; ih-- )
+   if (P2[ih] != 0)
+    {
+     last_bin = ih;
+     break;
+    }
+
+ /* Calculate the total entropy each gray-level
+  * and find the threshold that maximizes it
+  */
+ threshold = INT_MIN;
+ max_ent = DBL_MIN;
+ for ( it = first_bin; it <= last_bin; it++ )
+  {
+   /* Entropy of the background pixels */
+   ent_back = 0.0;
+   for ( ih = 0; ih <= it; ih++ )
+     if (norm_histo[ih] != 0)
+       ent_back -= ( norm_histo[ih] / P1[it] ) * log ( norm_histo[ih] / P1[it] );
+
+   /* Entropy of the object pixels */
+   ent_obj = 0.0;
+   for ( ih = it + 1; ih < HISTOGRAM_SIZE; ih++ )
+     if (norm_histo[ih] != 0)
+       ent_obj -= ( norm_histo[ih] / P2[it] ) * log ( norm_histo[ih] / P2[it] );
+
+   /* Total entropy */
+   tot_ent = ent_back + ent_obj;
+
+   if ( max_ent < tot_ent )
+    {
+     max_ent = tot_ent;
+     threshold = it;
+    }
+  }
+
+ if (threshold == INT_MIN)
+   {
+     DBG (5, "sanei_ir_threshold_maxentropy: no threshold found\n");
+     ret = SANE_STATUS_INVAL;
+   }
+ else
+   {
+     ret = SANE_STATUS_GOOD;
+     if (params->depth > 8)
+       {
+         it = 1 << (params->depth - HISTOGRAM_SHIFT);
+         *thresh = threshold * it + it / 2;
+       }
+     else
+       *thresh = threshold;
+     DBG (10, "sanei_ir_threshold_maxentropy: threshold %d\n", *thresh);
+ }
+
+ cleanup:
+   if (P1)
+     free (P1);
+   if (P2)
+     free (P2);
+   return ret;
+}
+
+/* Generate gray scale luminance image from separate R, G, B images
+ */
+SANE_Status
+sanei_ir_RGB_luminance (SANE_Parameters * params, const SANE_Uint **in_img,
+                      SANE_Uint **out_img)
+{
+  SANE_Uint *outi;
+  int itop, i;
+
+  if ((params->depth < 8) || (params->depth > 16) ||
+      (params->format != SANE_FRAME_GRAY))
+    {
+      DBG (5, "sanei_ir_RGB_luminance: invalid format\n");
+      return SANE_STATUS_UNSUPPORTED;
+    }
+
+  itop = params->pixels_per_line * params->lines;
+  outi = malloc (itop * sizeof(SANE_Uint));
+  if (!outi)
+    {
+      DBG (5, "sanei_ir_RGB_luminance: can not allocate out_img\n");
+      return SANE_STATUS_NO_MEM;
+    }
+
+  for (i = itop; i > 0; i--)
+      *(outi++) = (218 * (int) *(in_img[0]++) +
+                   732 * (int) *(in_img[1]++) +
+                   74 * (int) *(in_img[2]++)) >> 10;
+  *out_img = outi;
+  return SANE_STATUS_GOOD;
+}
+
+/* Convert image from >8 bit depth to an 8 bit image
+ */
+SANE_Status
+sanei_ir_to_8bit (SANE_Parameters * params, const SANE_Uint *in_img,
+                 SANE_Parameters * out_params, SANE_Uint **out_img)
+{
+  SANE_Uint *outi;
+  size_t ssize;
+  int i, is;
+
+  if ((params->depth < 8) || (params->depth > 16))
+    {
+      DBG (5, "sanei_ir_to_8bit: invalid format\n");
+      return SANE_STATUS_UNSUPPORTED;
+    }
+  ssize = params->pixels_per_line * params->lines;
+  if (params->format == SANE_FRAME_RGB)
+    ssize *= 3;
+  outi = malloc (ssize * sizeof(SANE_Uint));
+  if (!outi)
+    {
+      DBG (5, "sanei_ir_to_8bit: can not allocate out_img\n");
+      return SANE_STATUS_NO_MEM;
+    }
+
+  if (out_params)
+    {
+      memmove (out_params, params, sizeof(SANE_Parameters));
+      out_params->bytes_per_line = out_params->pixels_per_line;
+      if (params->format == SANE_FRAME_RGB)
+        out_params->bytes_per_line *= 3;
+      out_params->depth = 8;
+    }
+
+  memmove (outi, in_img, ssize * sizeof(SANE_Uint));
+  is = params->depth - 8;
+  for (i = ssize; i > 0; i--) {
+    *outi = *outi >> is, outi += 2;
+  }
+
+  *out_img = outi;
+  return SANE_STATUS_GOOD;
+}
+
+/* allocate and initialize logarithmic lookup table
+ */
+SANE_Status
+sanei_ir_ln_table (int len, double **lut_ln)
+{
+  double *llut;
+  int i;
+
+  DBG (10, "sanei_ir_ln_table\n");
+
+  llut = malloc (len * sizeof (double));
+  if (!llut)
+    {
+      DBG (5, "sanei_ir_ln_table: no table\n");
+      return SANE_STATUS_NO_MEM;
+    }
+  llut[0] = 0;
+  llut[1] = 0;
+  for (i = 2; i < len; i++)
+    llut[i] = log ((double) i);
+
+  *lut_ln = llut;
+  return SANE_STATUS_GOOD;
+}
+
+
+/* Reduce red spectral overlap from an infrared image plane
+ */
+SANE_Status
+sanei_ir_spectral_clean (const SANE_Parameters * params, double *lut_ln,
+			const SANE_Uint *red_data,
+			SANE_Uint *ir_data)
+{
+  const SANE_Uint *rptr;
+  SANE_Uint *iptr;
+  SANE_Int depth;
+  double *llut;
+  double rval, rsum, rrsum;
+  double risum, rfac, radd;
+  double *norm_histo;
+  int64_t isum;
+  int *calc_buf, *calc_ptr;
+  int ival, imin, imax;
+  int itop, len, ssize;
+  int thresh_low, thresh;
+  int irand, i;
+  SANE_Status status;
+
+  DBG (10, "sanei_ir_spectral_clean\n");
+
+  itop = params->pixels_per_line * params->lines;
+  calc_buf = malloc (itop * sizeof (int));		/* could save this */
+  if (!calc_buf)
+    {
+      DBG (5, "sanei_ir_spectral_clean: no buffer\n");
+      return SANE_STATUS_NO_MEM;
+    }
+
+  depth = params->depth;
+  len = 1 << depth;
+  if (lut_ln)
+    llut = lut_ln;
+  else
+    {
+      status = sanei_ir_ln_table (len, &llut);
+      if (status != SANE_STATUS_GOOD) {
+        free (calc_buf);
+        return status;
+      }
+    }
+
+  /* determine not transparent areas to exclude them later
+   * TODO: this has not been tested for negatives
+   */
+  thresh_low = INT_MAX;
+  status =
+      sanei_ir_create_norm_histogram (params, ir_data, &norm_histo);
+  if (status != SANE_STATUS_GOOD)
+    {
+      DBG (5, "sanei_ir_spectral_clean: no buffer\n");
+      free (calc_buf);
+      return SANE_STATUS_NO_MEM;
+    }
+
+  /* TODO: remember only needed if cropping is not ok */
+  status = sanei_ir_threshold_maxentropy (params, norm_histo, &thresh);
+  if (status == SANE_STATUS_GOOD)
+    thresh_low = thresh;
+  status = sanei_ir_threshold_otsu (params, norm_histo, &thresh);
+  if ((status == SANE_STATUS_GOOD) && (thresh < thresh_low))
+    thresh_low = thresh;
+  status = sanei_ir_threshold_yen (params, norm_histo, &thresh);
+  if ((status == SANE_STATUS_GOOD) && (thresh < thresh_low))
+    thresh_low = thresh;
+  if (thresh_low == INT_MAX)
+    thresh_low = 0;
+  else
+    thresh_low /= 2;
+  DBG (10, "sanei_ir_spectral_clean: low threshold %d\n", thresh_low);
+
+  /* calculate linear regression ired (red) from randomly chosen points */
+  ssize = itop / 2;
+  if (SAMPLE_SIZE < ssize)
+    ssize = SAMPLE_SIZE;
+  isum = 0;
+  rsum = rrsum = risum = 0.0;
+  i = ssize;
+  while (i > 0)
+    {
+      irand = rand () % itop;
+      rval = llut[red_data[irand]];
+      ival = ir_data[irand];
+      if (ival > thresh_low)
+        {
+          isum += ival;
+          rsum += rval;
+          rrsum += rval * rval;
+          risum += rval * (double) ival;
+          i--;
+        }
+    }
+
+  /* "a" in ired = b + a * ln (red) */
+  rfac =
+    ((double) ssize * risum -
+    rsum * (double) isum) / ((double) ssize * rrsum - rsum * rsum);
+    radd = ((double) isum - rfac * rsum) / (double) ssize;      /* "b" unused */
+
+  DBG (10, "sanei_ir_spectral_clean: n = %d, ired(red) = %f * ln(red) + %f\n",
+            ssize, rfac, radd);
+
+  /* now calculate ired' = ired - a  * ln (red) */
+  imin = INT_MAX;
+  imax = INT_MIN;
+  rptr = red_data;
+  iptr = ir_data;
+  calc_ptr = calc_buf;
+    for (i = itop; i > 0; i--)
+      {
+	ival = *iptr++ - (int) (rfac * llut[*rptr++] + 0.5);
+	if (ival > imax)
+	  imax = ival;
+	if (ival < imin)
+	  imin = ival;
+	*calc_ptr++ = ival;
+      }
+
+  /* scale the result back into the ired image */
+  calc_ptr = calc_buf;
+  iptr = ir_data;
+  rfac = (double) (len - 1) / (double) (imax - imin);
+    for (i = itop; i > 0; i--)
+      *iptr++ = (double) (*calc_ptr++ - imin) * rfac;
+
+  if (!lut_ln)
+    free (llut);
+  free (calc_buf);
+  free (norm_histo);
+  return SANE_STATUS_GOOD;
+}
+
+
+/* Hopefully fast mean filter
+ * JV: what does this do? Remove local mean?
+ */
+SANE_Status
+sanei_ir_filter_mean (const SANE_Parameters * params,
+		      const SANE_Uint *in_img, SANE_Uint *out_img,
+		      int win_rows, int win_cols)
+{
+  const SANE_Uint *src;
+  SANE_Uint *dest;
+  int num_cols, num_rows;
+  int itop, iadd, isub;
+  int ndiv, the_sum;
+  int nrow, ncol;
+  int hwr, hwc;
+  int *sum;
+  int i, j;
+
+  DBG (10, "sanei_ir_filter_mean, window: %d x%d\n", win_rows, win_cols);
+
+  if (((win_rows & 1) == 0) || ((win_cols & 1) == 0))
+    {
+      DBG (5, "sanei_ir_filter_mean: window even sized\n");
+      return SANE_STATUS_INVAL;
+    }
+
+  num_cols = params->pixels_per_line;
+  num_rows = params->lines;
+
+  sum = malloc (num_cols * sizeof (int));
+  if (!sum)
+    {
+      DBG (5, "sanei_ir_filter_mean: no buffer for sums\n");
+      return SANE_STATUS_NO_MEM;
+    }
+  dest = out_img;
+
+  hwr = win_rows / 2;		/* half window sizes */
+  hwc = win_cols / 2;
+
+  /* pre-pre calculation */
+  for (j = 0; j < num_cols; j++)
+    {
+        sum[j] = 0;
+	src = in_img + j;
+	for (i = 0; i < hwr; i++)
+	  {
+	    sum[j] += *src;
+	    src += num_cols;
+	  }
+    }
+
+  itop = num_rows * num_cols;
+  iadd = hwr * num_cols;
+  isub = (hwr - win_rows) * num_cols;
+  nrow = hwr;
+
+      for (i = 0; i < num_rows; i++)
+	{
+	  /* update row sums if possible */
+	  if (isub >= 0)	/* subtract old row */
+	    {
+	      nrow--;
+	      src = in_img + isub;
+	      for (j = 0; j < num_cols; j++)
+		sum[j] -= *src++;
+	    }
+	  isub += num_cols;
+
+	  if (iadd < itop)	/* add new row */
+	    {
+	      nrow++;
+	      src = in_img + iadd;
+	      for (j = 0; j < num_cols; j++)
+		sum[j] += *src++;
+	    }
+	  iadd += num_cols;
+
+	  /* now we do the image columns using only the precalculated sums */
+
+	  the_sum = 0;		/* precalculation */
+	  for (j = 0; j < hwc; j++)
+	    the_sum += sum[j];
+	  ncol = hwc;
+
+	  /* at the left margin, real index hwc lower */
+	  for (j = hwc; j < win_cols; j++)
+	    {
+	      ncol++;
+	      the_sum += sum[j];
+	      *dest++ = the_sum / (ncol * nrow);
+	    }
+
+	  ndiv = ncol * nrow;
+	  /* in the middle, real index hwc + 1 higher */
+	  for (j = 0; j < num_cols - win_cols; j++)
+	    {
+	      the_sum -= sum[j];
+	      the_sum += sum[j + win_cols];
+	      *dest++ = the_sum / ndiv;
+	    }
+
+	  /* at the right margin, real index hwc + 1 higher */
+	  for (j = num_cols - win_cols; j < num_cols - hwc - 1; j++)
+	    {
+	      ncol--;
+	      the_sum -= sum[j];	/* j - hwc - 1 */
+	      *dest++ = the_sum / (ncol * nrow);
+	    }
+	}
+  free (sum);
+  return SANE_STATUS_GOOD;
+}
+
+
+/* Find noise by adaptive thresholding
+ */
+SANE_Status
+sanei_ir_filter_madmean (const SANE_Parameters * params,
+			 const SANE_Uint *in_img,
+			 SANE_Uint ** out_img, int win_size,
+			 int a_val, int b_val)
+{
+  SANE_Uint *delta_ij, *delta_ptr;
+  SANE_Uint *mad_ij;
+  const SANE_Uint *mad_ptr;
+  SANE_Uint *out_ij, *dest8;
+  double ab_term;
+  int num_rows, num_cols;
+  int threshold, itop;
+  size_t size;
+  int ival, i;
+  int depth;
+  SANE_Status ret = SANE_STATUS_NO_MEM;
+
+  DBG (10, "sanei_ir_filter_madmean\n");
+
+  depth = params->depth;
+  if (depth != 8)
+    {
+      a_val = a_val << (depth - 8);
+      b_val = b_val << (depth - 8);
+    }
+  num_cols = params->pixels_per_line;
+  num_rows = params->lines;
+  itop = num_rows * num_cols;
+  size = itop * sizeof (SANE_Uint);
+  out_ij = malloc (size);
+  delta_ij = malloc (size);
+  mad_ij = malloc (size);
+
+  if (out_ij && delta_ij && mad_ij)
+    {
+      /* get the differences to the local mean */
+      mad_ptr = in_img;
+      if (sanei_ir_filter_mean (params, mad_ptr, delta_ij, win_size, win_size)
+	  == SANE_STATUS_GOOD)
+	{
+	  delta_ptr = delta_ij;
+	    for (i = 0; i < itop; i++)
+	      {
+		ival = *mad_ptr++ - *delta_ptr;
+		*delta_ptr++ = abs (ival);
+	      }
+	  /* make the second filtering window a bit larger */
+	  win_size = MAD_WIN2_SIZE(win_size);
+	  /* and get the local mean differences */
+	  if (sanei_ir_filter_mean
+	      (params, delta_ij, mad_ij, win_size,
+	       win_size) == SANE_STATUS_GOOD)
+	    {
+	      mad_ptr = mad_ij;
+	      delta_ptr = delta_ij;
+	      dest8 = out_ij;
+	      /* construct the noise map */
+	      ab_term = (b_val - a_val) / (double) b_val;
+		for (i = 0; i < itop; i++)
+		  {
+		    /* by calculating the threshold */
+		    ival = *mad_ptr++;
+		    if (ival >= b_val)	/* outlier */
+		      threshold = a_val;
+		    else
+		      threshold = a_val + (double) ival *ab_term;
+		    /* above threshold is noise, indicated by 0 */
+		    if (*delta_ptr++ >= threshold)
+		      *dest8++ = 0;
+		    else
+		      *dest8++ = 255;
+		  }
+	      *out_img = out_ij;
+	      ret = SANE_STATUS_GOOD;
+	    }
+	}
+    }
+  else
+    DBG (5, "sanei_ir_filter_madmean: Cannot allocate buffers\n");
+
+  free (mad_ij);
+  free (delta_ij);
+  return ret;
+}
+
+
+/* Add dark pixels to mask from static threshold
+ */
+void
+sanei_ir_add_threshold (const SANE_Parameters * params,
+			const SANE_Uint *in_img,
+			SANE_Uint * mask_img, int threshold)
+{
+  const SANE_Uint *in_ptr;
+  SANE_Uint *mask_ptr;
+  int itop, i;
+
+  DBG (10, "sanei_ir_add_threshold\n");
+
+  itop = params->pixels_per_line * params->lines;
+  in_ptr = in_img;
+  mask_ptr = mask_img;
+
+    for (i = itop; i > 0; i--)
+      {
+	if (*in_ptr++ <= threshold)
+	  *mask_ptr = 0;
+	mask_ptr++;
+      }
+}
+
+
+/* Calculate minimal Manhattan distances for an image mask
+ */
+void
+sanei_ir_manhattan_dist (const SANE_Parameters * params,
+			const SANE_Uint * mask_img, unsigned int *dist_map,
+			unsigned int *idx_map, unsigned int erode)
+{
+  const SANE_Uint *mask;
+  unsigned int *index, *manhattan;
+  int rows, cols, itop;
+  int i, j;
+
+  DBG (10, "sanei_ir_manhattan_dist\n");
+
+  if (erode != 0)
+    erode = 255;
+
+  /* initialize maps */
+  cols = params->pixels_per_line;
+  rows = params->lines;
+  itop = rows * cols;
+  mask = mask_img;
+  manhattan = dist_map;
+  index = idx_map;
+  for (i = 0; i < itop; i++)
+    {
+      *manhattan++ = *mask++;
+      *index++ = i;
+    }
+
+  /* traverse from top left to bottom right */
+  manhattan = dist_map;
+  index = idx_map;
+  for (i = 0; i < rows; i++)
+    for (j = 0; j < cols; j++)
+      {
+	if (*manhattan == erode)
+	  {
+	    /* take original, distance = 0, index stays the same */
+	    *manhattan = 0;
+	  }
+	else
+	  {
+	    /* assume maximal distance to clean pixel */
+	    *manhattan = cols + rows;
+	    /* or one further away than pixel to the top */
+	    if (i > 0)
+	      if (manhattan[-cols] + 1 < *manhattan)
+		{
+		  *manhattan = manhattan[-cols] + 1;
+		  *index = index[-cols];	/* index follows */
+		}
+	    /* or one further away than pixel to the left */
+	    if (j > 0)
+	      {
+		if (manhattan[-1] + 1 < *manhattan)
+		  {
+		    *manhattan = manhattan[-1] + 1;
+		    *index = index[-1];	/* index follows */
+		  }
+		if (manhattan[-1] + 1 == *manhattan)
+		  if (rand () % 2 == 0)	/* chose index */
+		    *index = index[-1];
+	      }
+	  }
+	manhattan++;
+	index++;
+      }
+
+  /* traverse from bottom right to top left */
+  manhattan = dist_map + itop - 1;
+  index = idx_map + itop - 1;
+  for (i = rows - 1; i >= 0; i--)
+    for (j = cols - 1; j >= 0; j--)
+      {
+	if (i < rows - 1)
+	  {
+	    /* either what we had on the first pass
+	       or one more than the pixel to the bottm */
+	    if (manhattan[+cols] + 1 < *manhattan)
+	      {
+		*manhattan = manhattan[+cols] + 1;
+		*index = index[+cols];	/* index follows */
+	      }
+	    if (manhattan[+cols] + 1 == *manhattan)
+	      if (rand () % 2 == 0)	/* chose index */
+		*index = index[+cols];
+	  }
+	if (j < cols - 1)
+	  {
+	    /* or one more than pixel to the right */
+	    if (manhattan[1] + 1 < *manhattan)
+	      {
+		*manhattan = manhattan[1] + 1;
+		*index = index[1];	/* index follows */
+	      }
+	    if (manhattan[1] + 1 == *manhattan)
+	      if (rand () % 2 == 0)	/* chose index */
+		*index = index[1];
+	  }
+	manhattan--;
+	index--;
+      }
+}
+
+
+/* dilate or erode a mask image */
+
+void
+sanei_ir_dilate (const SANE_Parameters *params, SANE_Uint *mask_img,
+		unsigned int *dist_map, unsigned int *idx_map, int by)
+{
+  SANE_Uint *mask;
+  unsigned int *manhattan;
+  unsigned int erode;
+  unsigned int thresh;
+  int i, itop;
+
+  DBG (10, "sanei_ir_dilate\n");
+
+  if (by == 0)
+    return;
+  if (by > 0)
+    {
+      erode = 0;
+      thresh = by;
+    }
+  else
+    {
+      erode = 1;
+      thresh = -by;
+    }
+
+  itop = params->pixels_per_line * params->lines;
+  mask = mask_img;
+  sanei_ir_manhattan_dist (params, mask_img, dist_map, idx_map, erode);
+
+  manhattan = dist_map;
+  for (i = 0; i < itop; i++)
+    {
+      if (*manhattan++ <= thresh)
+	*mask++ = 0;
+      else
+	*mask++ = 255;
+    }
+
+  return;
+}
+
+
+/* Suggest cropping for dark margins of positive film
+ */
+void
+sanei_ir_find_crop (const SANE_Parameters * params,
+                    unsigned int * dist_map, int inner, int * edges)
+{
+  int width = params->pixels_per_line;
+  int height = params->lines;
+  uint64_t sum_x, sum_y, n;
+  int64_t sum_xx, sum_xy;
+  double a, b, mami;
+  unsigned int *src;
+  int off1, off2, inc, wh, i, j;
+
+  DBG (10, "sanei_ir_find_crop\n");
+
+  /* loop through top, bottom, left, right */
+  for (j = 0; j < 4; j++)
+    {
+      if (j < 2)        /* top, bottom */
+        {
+          off1 = width / 8;     /* only middle 3/4 */
+          off2 = width - off1;
+          n = width - 2 * off1;
+          src = dist_map + off1;        /* first row */
+          inc = 1;
+          wh = width;
+          if (j == 1)                   /* last row */
+            src += (height - 1) * width;
+        }
+      else              /* left, right */
+        {
+          off1 = height / 8;     /* only middle 3/4 */
+          off2 = height - off1;
+          n = height - 2 * off1;
+          src = dist_map + (off1 * width);      /* first column */
+          inc = width;
+          wh = height;
+          if (j == 3)
+            src += width - 1;                   /* last column */
+        }
+
+      /* calculate linear regression */
+      sum_x = 0; sum_y = 0;
+      sum_xx = 0; sum_xy = 0;
+      for (i = off1; i < off2; i++)
+        {
+          sum_x += i;
+          sum_y += *src;
+          sum_xx += i * i;
+          sum_xy += i * (*src);
+          src += inc;
+        }
+      b = ((double) n * (double) sum_xy - (double) sum_x * (double) sum_y)
+          / ((double) n * (double) sum_xx - (double) sum_x * (double) sum_x);
+      a = ((double) sum_y - b * (double) sum_x) / (double) n;
+
+      DBG (10, "sanei_ir_find_crop: y = %f + %f * x\n", a, b);
+
+      /* take maximal/minimal value from either side */
+      mami = a + b * (wh - 1);
+      if (inner)
+        {
+          if (a > mami)
+            mami = a;
+        }
+      else
+        {
+          if (a < mami)
+            mami = a;
+        }
+      edges[j] = mami + 0.5;
+    }
+  edges[1] = height - edges[1];
+  edges[3] = width - edges[3];
+
+  DBG (10, "sanei_ir_find_crop: would crop at top: %d, bot: %d, left %d, right %d\n",
+      edges[0], edges[1], edges[2], edges[3]);
+
+  return;
+}
+
+
+/* Dilate clean image parts into dirty ones and smooth
+ */
+SANE_Status
+sanei_ir_dilate_mean (const SANE_Parameters * params,
+                      SANE_Uint **in_img,
+                      SANE_Uint * mask_img,
+                      int dist_max, int expand, int win_size,
+                      SANE_Bool smooth, int inner,
+                      int *crop)
+{
+  SANE_Uint *color;
+  SANE_Uint *plane;
+  unsigned int *dist_map, *manhattan;
+  unsigned int *idx_map, *index;
+  int dist;
+  int rows, cols;
+  int k, i, itop;
+  SANE_Status ret = SANE_STATUS_NO_MEM;
+
+  DBG (10, "sanei_ir_dilate_mean(): dist max = %d, expand = %d, win size = %d, smooth = %d, inner = %d\n",
+    dist_max, expand, win_size, smooth, inner);
+
+  cols = params->pixels_per_line;
+  rows = params->lines;
+  itop = rows * cols;
+  idx_map = malloc (itop * sizeof (unsigned int));
+  dist_map = malloc (itop * sizeof (unsigned int));
+  plane = malloc (itop * sizeof (SANE_Uint));
+
+  if (!idx_map || !dist_map || !plane)
+    DBG (5, "sanei_ir_dilate_mean: Cannot allocate buffers\n");
+  else
+    {
+      /* expand dirty regions into their half dirty surround*/
+      if (expand > 0)
+	sanei_ir_dilate (params, mask_img, dist_map, idx_map, expand);
+      /* for dirty pixels determine the closest clean ones */
+      sanei_ir_manhattan_dist (params, mask_img, dist_map, idx_map, 1);
+
+      /* use the distance map to find how to crop dark edges */
+      if (crop)
+        sanei_ir_find_crop (params, dist_map, inner, crop);
+
+      /* replace dirty pixels */
+      for (k = 0; k < 3; k++)
+	{
+	  manhattan = dist_map;
+	  index = idx_map;
+	  color = in_img[k];
+	  /* first replacement */
+	    for (i = 0; i < itop; i++)
+	      {
+		dist = *manhattan++;
+		if ((dist != 0) && (dist <= dist_max))
+		  color[i] = color[index[i]];
+	      }
+          /* adapt pixels to their new surround and
+           * smooth the whole image or the replaced pixels only */
+	  ret =
+	    sanei_ir_filter_mean (params, color, plane, win_size, win_size);
+	  if (ret != SANE_STATUS_GOOD)
+	    break;
+	  else
+	    if (smooth)
+              {
+                /* a second mean results in triangular blur */
+                DBG (10, "sanei_ir_dilate_mean(): smoothing whole image\n");
+                ret =
+                  sanei_ir_filter_mean (params, plane, color, win_size,
+                                        win_size);
+                if (ret != SANE_STATUS_GOOD)
+                  break;
+              }
+            else
+              {
+                /* replace with smoothened pixels only */
+                DBG (10, "sanei_ir_dilate_mean(): smoothing replaced pixels only\n");
+                manhattan = dist_map;
+                  for (i = 0; i < itop; i++)
+                    {
+                      dist = *manhattan++;
+                      if ((dist != 0) && (dist <= dist_max))
+                        color[i] = plane[i];
+                    }
+              }
+      }
+    }
+  free (plane);
+  free (dist_map);
+  free (idx_map);
+
+  return ret;
+}