sane-backends-1.0.32/sanei/sanei_ir.c

/** @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, see <https://www.gnu.org/licenses/>.
 *
 * 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 <float.h>
#include <limits.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;
}