only-pie/sanei/sanei_ir.c

Go to the documentation of this file.

#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 SANEI_IR_bufptr 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 SANEI_IR_bufptr img_data)
{
  uint8_t *img_data8;
  uint16_t *img_data16;
  int is, i;
  int num_pixels;
  int *histo_data = NULL;
  double *histo = NULL;
  double term;

  DBG (10, "sanei_ir_create_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");
      free (histo);
      free (histo_data);
      return NULL;
    }

  num_pixels = params->pixels_per_line * params->lines;

  /* Populate the histogram */
  is = params->depth - HISTOGRAM_SHIFT;
  if (params->depth > 8)
    {
      img_data16 = img_data.b16;
      if (is == 8)
        for (i = num_pixels; i > 0; i--)
          histo_data[*img_data16++ >> 8]++;
      else
        for (i = num_pixels; i > 0; i--)
          histo_data[*img_data16++ >> is]++;
    }
  else
    {
      img_data8 = img_data.b8;
      if (is == 0)
        for (i = num_pixels; i > 0; i--)
          histo_data[*img_data8++]++;
      else
        for (i = num_pixels; i > 0; i--)
          histo_data[*img_data8++ >> 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 SANEI_IR_bufptr img_data,
                                double ** histogram)
{
  double *histo;

  DBG (10, "sanei_ir_create_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
sane_ir_RGB_luminance (SANE_Parameters * params, const SANEI_IR_bufptr * in_img,
                      SANEI_IR_bufptr * out_img, SANE_Bool scale8)
{
  SANEI_IR_bufptr cplane[3];    /* R, G, B gray scale planes */
  SANEI_IR_bufptr outi, dest;
  size_t ssize;
  int itop, is, i;

  if ((params->depth < 8) || (params->depth > 16) ||
      (params->format != SANE_FRAME_GRAY))
    {
      DBG (5, "sane_ir_RGB_luminance: invalid format\n");
      return SANE_STATUS_UNSUPPORTED;
    }

  for (i = 0; i < 3; i++)
    cplane[i] = in_img[i];
  itop = params->pixels_per_line * params->lines;
  ssize = itop;
  if ((params->depth > 8) && !scale8)
    ssize *= 2;
  outi.b8 = malloc (ssize);
  if (!outi.b8)
    {
      DBG (5, "sane_ir_to_8bit: can not allocate out_img\n");
      return SANE_STATUS_NO_MEM;
    }
  dest = outi;

  if (params->depth > 8)
    {
      if (scale8)
        {
          is = params->depth + 2;
          for (i = ssize; i > 0; i--)
            *(dest.b8)++ = (218 * (int) *(cplane[0].b16++) +
                            732 * (int) *(cplane[1].b16++) +
                            74 * (int) *(cplane[2].b16++)) >> is;
          params->depth = 8;
          params->bytes_per_line = params->pixels_per_line;
        }
      else
        for (i = ssize; i > 0; i--)
          *(dest.b16++) = (218 * (int) *(cplane[0].b16++) +
                           732 * (int) *(cplane[1].b16++) +
                           74 * (int) *(cplane[2].b16++)) >> 10;
    }
  else
    {
      for (i = ssize; i > 0; i--)
        *(dest.b8++) = (218 * (int) *(cplane[0].b8++) +
                        732 * (int) *(cplane[1].b8++) +
                        74 * (int) *(cplane[2].b8++)) >> 10;
    }

  *out_img = outi;
  return SANE_STATUS_GOOD;
}

/* Convert image from >8 bit depth to an 8 bit image
 */
SANE_Status
sane_ir_to_8bit (SANE_Parameters * params, const SANEI_IR_bufptr in_img,
                 SANE_Parameters * out_params, SANEI_IR_bufptr * out_img)
{
  SANEI_IR_bufptr outi, src, dest;
  size_t ssize;
  int i, is;

  if ((params->depth < 8) || (params->depth > 16))
    {
      DBG (5, "sane_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.b8 = malloc (ssize);
  if (!outi.b8)
    {
      DBG (5, "sane_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;
    }
  if (params->depth == 8)
    memmove (outi.b8, in_img.b8, ssize);
  else
    {
      dest.b8 = outi.b8;
      src.b16 = in_img.b16;
      if (params->depth == 16)
        {
          for (i = ssize; i > 0; i--)
            *dest.b8++ = *src.b16++ >> 8;
        }
      else
        {
          is = params->depth - 8;
          for (i = ssize; i > 0; i--)
            *dest.b8++ = *src.b16++ >> is;
        }
    }

  *out_img = outi;
  return SANE_STATUS_GOOD;
}

/* allocate and initialize logarithmic lookup table
 */
SANE_Status
sane_ir_ln_table (int len, double **lut_ln)
{
  double *llut;
  int i;

  DBG (10, "sane_ir_ln_table\n");

  llut = malloc (len * sizeof (double));
  if (!llut)
    {
      DBG (5, "sane_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
sane_ir_spectral_clean (const SANE_Parameters * params, double *lut_ln,
                        const SANEI_IR_bufptr red_data,
                        SANEI_IR_bufptr ir_data)
{
  SANEI_IR_bufptr rptr, 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, "sane_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, "sane_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 = sane_ir_ln_table (len, &llut);
      if (status != SANE_STATUS_GOOD)
        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, "sane_ir_spectral_clean: no buffer\n");
      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, "sane_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;
      if (depth > 8)
        {
          rval = llut[red_data.b16[irand]];
          ival = ir_data.b16[irand];
        }
      else
        {
          rval = llut[red_data.b8[irand]];
          ival = ir_data.b8[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, "sane_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;
  if (depth > 8)
    for (i = itop; i > 0; i--)
      {
        ival = *iptr.b16++ - (int) (rfac * llut[*rptr.b16++] + 0.5);
        if (ival > imax)
          imax = ival;
        if (ival < imin)
          imin = ival;
        *calc_ptr++ = ival;
      }
  else
    for (i = itop; i > 0; i--)
      {
        ival = *iptr.b8++ - (int) (rfac * llut[*rptr.b8++] + 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);
  if (depth > 8)
    for (i = itop; i > 0; i--)
      *iptr.b16++ = (double) (*calc_ptr++ - imin) * rfac;
  else
    for (i = itop; i > 0; i--)
      *iptr.b8++ = (double) (*calc_ptr++ - imin) * rfac;

  if (!lut_ln)
    free (llut);
  free (calc_buf);
  free (norm_histo);
  return SANE_STATUS_GOOD;
}


/* Hopefully fast mean filter
 */
SANE_Status
sanei_ir_filter_mean (const SANE_Parameters * params,
                      const SANEI_IR_bufptr in_img, SANEI_IR_bufptr out_img,
                      int win_rows, int win_cols)
{
  SANEI_IR_bufptr src, dest;
  int num_cols, num_rows;
  int itop, iadd, isub;
  int ndiv, the_sum;
  int nrow, ncol;
  int hwr, hwc;
  int depth;
  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;
  depth = params->depth;

  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;

  if (depth > 8)                /* pre-pre calculation */
    for (j = 0; j < num_cols; j++)
      {
        sum[j] = 0;
        src.b16 = in_img.b16 + j;
        for (i = 0; i < hwr; i++)
          {
            sum[j] += *src.b16;
            src.b16 += num_cols;
          }
      }
  else
    for (j = 0; j < num_cols; j++)
      {
        sum[j] = 0;
        src.b8 = in_img.b8 + j;
        for (i = 0; i < hwr; i++)
          {
            sum[j] += *src.b8;
            src.b8 += num_cols;
          }
      }

  itop = num_rows * num_cols;
  iadd = hwr * num_cols;
  isub = (hwr - win_rows) * num_cols;
  nrow = hwr;

  if (depth > 8)
    {
      for (i = 0; i < num_rows; i++)
        {
          /* update row sums if possible */
          if (isub >= 0)        /* subtract old row */
            {
              nrow--;
              src.b16 = in_img.b16 + isub;
              for (j = 0; j < num_cols; j++)
                sum[j] -= *src.b16++;
            }
          isub += num_cols;

          if (iadd < itop)      /* add new row */
            {
              nrow++;
              src.b16 = in_img.b16 + iadd;
              for (j = 0; j < num_cols; j++)
                sum[j] += *src.b16++;
            }
          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.b16++ = 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.b16++ = 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.b16++ = the_sum / (ncol * nrow);
            }
        }
    }
  else
    {
      for (i = 0; i < num_rows; i++)
        {
          /* update row sums if possible */
          if (isub >= 0)        /* subtract old row */
            {
              nrow--;
              src.b8 = in_img.b8 + isub;
              for (j = 0; j < num_cols; j++)
                sum[j] -= *src.b8++;
            }
          isub += num_cols;

          if (iadd < itop)      /* add new row */
            {
              nrow++;
              src.b8 = in_img.b8 + iadd;
              for (j = 0; j < num_cols; j++)
                sum[j] += *src.b8++;
            }
          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.b8++ = 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.b8++ = 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.b8++ = the_sum / (ncol * nrow);
            }
        }
    }
  return SANE_STATUS_GOOD;
}


/* Find noise by adaptive thresholding
 */
SANE_Status
sanei_ir_filter_madmean (const SANE_Parameters * params,
                         const SANEI_IR_bufptr in_img,
                         SANE_Byte ** out_img, int win_size,
                         int a_val, int b_val)
{
  SANEI_IR_bufptr delta_ij, delta_ptr;
  SANEI_IR_bufptr mad_ij, mad_ptr;
  SANE_Byte *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 (uint8_t);
  out_ij = malloc (size);
  if (depth > 8)
    size = itop * sizeof (uint16_t);
  delta_ij.b8 = malloc (size);
  mad_ij.b8 = malloc (size);

  if (out_ij && delta_ij.b8 && mad_ij.b8)
    {
      /* 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;
          if (depth > 8)
            for (i = 0; i < itop; i++)
              {
                ival = *mad_ptr.b16++ - *delta_ptr.b16;
                *delta_ptr.b16++ = abs (ival);
              }
          else
            for (i = 0; i < itop; i++)
              {
                ival = *mad_ptr.b8++ - *delta_ptr.b8;
                *delta_ptr.b8++ = 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;
              if (depth > 8)
                for (i = 0; i < itop; i++)
                  {
                    /* by calculating the threshold */
                    ival = *mad_ptr.b16++;
                    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.b16++ >= threshold)
                      *dest8++ = 0;
                    else
                      *dest8++ = 255;
                  }
              else
                for (i = 0; i < itop; i++)
                  {
                    /* by calculating the threshold */
                    ival = *mad_ptr.b8++;
                    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.b8++ >= 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.b8);
  free (delta_ij.b8);
  return ret;
}


/* Add dark pixels to mask from static threshold
 */
void
sanei_ir_add_threshold (const SANE_Parameters * params,
                        const SANEI_IR_bufptr in_img,
                        SANE_Byte * mask_img, int threshold)
{
  SANEI_IR_bufptr in_ptr;
  SANE_Byte *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;

  if (params->depth > 8)
    for (i = itop; i > 0; i--)
      {
        if (*in_ptr.b16++ <= threshold)
          *mask_ptr = 0;
        mask_ptr++;
      }
  else
    for (i = itop; i > 0; i--)
      {
        if (*in_ptr.b8++ <= threshold)
          *mask_ptr = 0;
        mask_ptr++;
      }
}


/* Calculate minimal Manhattan distances for an image mask
 */
void
sane_ir_manhattan_dist (const SANE_Parameters * params,
                        const SANE_Byte * mask_img, unsigned int *dist_map,
                        unsigned int *idx_map, unsigned int erode)
{
  const SANE_Byte *mask;
  unsigned int *index, *manhattan;
  int rows, cols, itop;
  int i, j;

  DBG (10, "sane_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
sane_ir_dilate (const SANE_Parameters * params, SANE_Byte * mask_img,
                unsigned int *dist_map, unsigned int *idx_map, int by)
{
  SANE_Byte *mask;
  unsigned int *manhattan;
  unsigned int erode;
  unsigned int thresh;
  int i, itop;

  DBG (10, "sane_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;
  sane_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,
                      SANEI_IR_bufptr * in_img,
                      SANE_Byte * mask_img,
                      int dist_max, int expand, int win_size,
                      SANE_Bool smooth, int inner,
                      int *crop)
{
  SANEI_IR_bufptr color;
  SANEI_IR_bufptr plane;
  unsigned int *dist_map, *manhattan;
  unsigned int *idx_map, *index;
  int depth, dist;
  int rows, cols;
  int k, i, itop;
  SANE_Status ret = SANE_STATUS_NO_MEM;

  DBG (10, "sanei_ir_dilate_mean\n");

  depth = params->depth;
  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));
  if (depth > 8)
    plane.b16 = malloc (itop * sizeof (uint16_t));
  else
    plane.b8 = malloc (itop * sizeof (uint8_t));

  if (!idx_map || !dist_map || !plane.b8)
    DBG (5, "sanei_ir_dilate_mean: Cannot allocate buffers\n");
  else
    {
      /* expand dirty regions into their half dirty surround*/
      if (expand > 0)
        sane_ir_dilate (params, mask_img, dist_map, idx_map, expand);
      /* for dirty pixels determine the closest clean ones */
      sane_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 */
          if (depth > 8)
            for (i = 0; i < itop; i++)
              {
                dist = *manhattan++;
                if ((dist != 0) && (dist <= dist_max))
                  color.b16[i] = color.b16[index[i]];
              }
          else
            for (i = 0; i < itop; i++)
              {
                dist = *manhattan++;
                if ((dist != 0) && (dist <= dist_max))
                  color.b8[i] = color.b8[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 */
                ret =
                  sanei_ir_filter_mean (params, plane, color, win_size,
                                        win_size);
                if (ret != SANE_STATUS_GOOD)
                  break;
              }
            else
              {
                /* replace with smoothened pixels only */
                manhattan = dist_map;
                if (depth > 8)
                  for (i = 0; i < itop; i++)
                    {
                      dist = *manhattan++;
                      if ((dist != 0) && (dist <= dist_max))
                        color.b16[i] = plane.b16[i];
                    }
                else
                  for (i = 0; i < itop; i++)
                    {
                      dist = *manhattan++;
                      if ((dist != 0) && (dist <= dist_max))
                        color.b8[i] = plane.b8[i];
                    }
              }
      }
    }
  free (plane.b8);
  free (dist_map);
  free (idx_map);

  return ret;
}