mirror
/
sane-project-backends
kopia lustrzana https://gitlab.com/sane-project/backends
1
0
Forkuj 0
Kod Zgłoszenia Projekty Wydania Wiki Aktywność
sane-project-backends/sanei/sanei_ir.c

1205 wiersze
30 KiB
C
Czysty Bezpośredni odnośnik Wina Historia

/** @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;
}
Reference in New Issue View Git Blame Kopiuj bezpośredni odnośnik