summaryrefslogtreecommitdiff
path: root/sanei/sanei_ir.c
diff options
context:
space:
mode:
authorJörg Frings-Fürst <debian@jff-webhosting.net>2015-10-08 12:32:49 +0200
committerJörg Frings-Fürst <debian@jff-webhosting.net>2015-10-08 12:32:49 +0200
commitcfd27ef2ad8b005fd47ab41ef29b71d9e3d48201 (patch)
tree56cb8849a5df91b19b72fd4f48c79c343cd9bc5e /sanei/sanei_ir.c
parentba8d09abe681600aad991f4a75e904615b7ed29f (diff)
Imported Upstream version 1.0.25upstream/1.0.25
Diffstat (limited to 'sanei/sanei_ir.c')
-rw-r--r--sanei/sanei_ir.c1205
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;
+}