1/** @file sanei_ir.c 2 * 3 * sanei_ir - functions for utilizing the infrared plane 4 * 5 * Copyright (C) 2012 Michael Rickmann <mrickma@gwdg.de> 6 * 7 * This file is part of the SANE package. 8 * 9 * This program is free software; you can redistribute it and/or 10 * modify it under the terms of the GNU General Public License as 11 * published by the Free Software Foundation; either version 2 of the 12 * License, or (at your option) any later version. 13 * 14 * This program is distributed in the hope that it will be useful, but 15 * WITHOUT ANY WARRANTY; without even the implied warranty of 16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 17 * General Public License for more details. 18 * 19 * You should have received a copy of the GNU General Public License 20 * along with this program. If not, see <https://www.gnu.org/licenses/>. 21 * 22 * The threshold yen, otsu and max_entropy routines have been 23 * adapted from the FOURIER 0.8 library by M. Emre Celebi, 24 * http://sourceforge.net/projects/fourier-ipal/ which is 25 * licensed under the GNU General Public License version 2 or later. 26*/ 27 28#include <stdlib.h> 29#include <string.h> 30#include <float.h> 31#include <limits.h> 32#include <math.h> 33 34#define BACKEND_NAME sanei_ir /* name of this module for debugging */ 35 36#include "../include/sane/sane.h" 37#include "../include/sane/sanei_debug.h" 38#include "../include/sane/sanei_ir.h" 39#include "../include/sane/sanei_magic.h" 40 41 42double * 43sanei_ir_create_norm_histo (const SANE_Parameters * params, const SANE_Uint *img_data); 44double * sanei_ir_accumulate_norm_histo (double * histo_data); 45 46 47/* Initialize sanei_ir 48 */ 49void 50sanei_ir_init (void) 51{ 52 DBG_INIT (); 53} 54 55 56/* Create a normalized histogram of a grayscale image, internal 57 */ 58double * 59sanei_ir_create_norm_histo (const SANE_Parameters * params, 60 const SANE_Uint *img_data) 61{ 62 int is, i; 63 int num_pixels; 64 int *histo_data; 65 double *histo; 66 double term; 67 68 DBG (10, "sanei_ir_create_norm_histo\n"); 69 70 if ((params->format != SANE_FRAME_GRAY) 71 && (params->format != SANE_FRAME_RED) 72 && (params->format != SANE_FRAME_GREEN) 73 && (params->format != SANE_FRAME_BLUE)) 74 { 75 DBG (5, "sanei_ir_create_norm_histo: invalid format\n"); 76 return NULL; 77 } 78 79 /* Allocate storage for the histogram */ 80 histo_data = calloc (HISTOGRAM_SIZE, sizeof (int)); 81 histo = malloc (HISTOGRAM_SIZE * sizeof (double)); 82 if ((histo == NULL) || (histo_data == NULL)) 83 { 84 DBG (5, "sanei_ir_create_norm_histo: no buffers\n"); 85 if (histo) free (histo); 86 if (histo_data) free (histo_data); 87 return NULL; 88 } 89 90 num_pixels = params->pixels_per_line * params->lines; 91 92 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); 93 DBG (1, "sanei_ir_create_norm_histo: histo_data[] with %d x %ld bytes\n", HISTOGRAM_SIZE, sizeof(int)); 94 /* Populate the histogram */ 95 is = 16 - HISTOGRAM_SHIFT; /* Number of data bits to ignore */ 96 DBG (1, "sanei_ir_create_norm_histo: depth %d, HISTOGRAM_SHIFT %d => ignore %d bits\n", params->depth, HISTOGRAM_SHIFT, is); 97 for (i = num_pixels; i > 0; i--) { 98 histo_data[*img_data++ >> is]++; 99 } 100 101 /* Calculate the normalized histogram */ 102 term = 1.0 / (double) num_pixels; 103 for (i = 0; i < HISTOGRAM_SIZE; i++) 104 histo[i] = term * (double) histo_data[i]; 105 106 free (histo_data); 107 return histo; 108} 109 110 111/* Create the normalized histogram of a grayscale image 112 */ 113SANE_Status 114sanei_ir_create_norm_histogram (const SANE_Parameters * params, 115 const SANE_Uint *img_data, 116 double ** histogram) 117{ 118 double *histo; 119 120 DBG (10, "sanei_ir_create_norm_histogram\n"); 121 122 histo = sanei_ir_create_norm_histo (params, img_data); 123 if (!histo) 124 return SANE_STATUS_NO_MEM; 125 126 *histogram = histo; 127 return SANE_STATUS_GOOD; 128} 129 130/* Accumulate a normalized histogram, internal 131 */ 132double * 133sanei_ir_accumulate_norm_histo (double * histo_data) 134{ 135 int i; 136 double *accum_data; 137 138 accum_data = malloc (HISTOGRAM_SIZE * sizeof (double)); 139 if (accum_data == NULL) 140 { 141 DBG (5, "sanei_ir_accumulate_norm_histo: Insufficient memory !\n"); 142 return NULL; 143 } 144 145 accum_data[0] = histo_data[0]; 146 for (i = 1; i < HISTOGRAM_SIZE; i++) 147 accum_data[i] = accum_data[i - 1] + histo_data[i]; 148 149 return accum_data; 150} 151 152/* Implements Yen's thresholding method 153 */ 154SANE_Status 155sanei_ir_threshold_yen (const SANE_Parameters * params, 156 double * norm_histo, int *thresh) 157{ 158 double *P1 = NULL; /* cumulative normalized histogram */ 159 double *P1_sq = NULL; /* cumulative normalized histogram */ 160 double *P2_sq = NULL; 161 double crit, max_crit; 162 int threshold, i; 163 SANE_Status ret = SANE_STATUS_NO_MEM; 164 165 DBG (10, "sanei_ir_threshold_yen\n"); 166 167 P1 = sanei_ir_accumulate_norm_histo (norm_histo); 168 P1_sq = malloc (HISTOGRAM_SIZE * sizeof (double)); 169 P2_sq = malloc (HISTOGRAM_SIZE * sizeof (double)); 170 if (!P1 || !P1_sq || !P2_sq) 171 { 172 DBG (5, "sanei_ir_threshold_yen: no buffers\n"); 173 goto cleanup; 174 } 175 176 /* calculate cumulative squares */ 177 P1_sq[0] = norm_histo[0] * norm_histo[0]; 178 for (i = 1; i < HISTOGRAM_SIZE; i++) 179 P1_sq[i] = P1_sq[i - 1] + norm_histo[i] * norm_histo[i]; 180 P2_sq[HISTOGRAM_SIZE - 1] = 0.0; 181 for (i = HISTOGRAM_SIZE - 2; i >= 0; i--) 182 P2_sq[i] = P2_sq[i + 1] + norm_histo[i + 1] * norm_histo[i + 1]; 183 184 /* Find the threshold that maximizes the criterion */ 185 threshold = INT_MIN; 186 max_crit = DBL_MIN; 187 for (i = 0; i < HISTOGRAM_SIZE; i++) 188 { 189 crit = 190 -1.0 * SAFE_LOG (P1_sq[i] * P2_sq[i]) + 191 2 * SAFE_LOG (P1[i] * (1.0 - P1[i])); 192 if (crit > max_crit) 193 { 194 max_crit = crit; 195 threshold = i; 196 } 197 } 198 199 if (threshold == INT_MIN) 200 { 201 DBG (5, "sanei_ir_threshold_yen: no threshold found\n"); 202 ret = SANE_STATUS_INVAL; 203 } 204 else 205 { 206 ret = SANE_STATUS_GOOD; 207 if (params->depth > 8) 208 { 209 i = 1 << (params->depth - HISTOGRAM_SHIFT); 210 *thresh = threshold * i + i / 2; 211 } 212 else 213 *thresh = threshold; 214 DBG (10, "sanei_ir_threshold_yen: threshold %d\n", *thresh); 215 } 216 217 cleanup: 218 if (P1) 219 free (P1); 220 if (P1_sq) 221 free (P1_sq); 222 if (P2_sq) 223 free (P2_sq); 224 return ret; 225} 226 227 228/* Implements Otsu's thresholding method 229 */ 230SANE_Status 231sanei_ir_threshold_otsu (const SANE_Parameters * params, 232 double * norm_histo, int *thresh) 233{ 234 double *cnh = NULL; 235 double *mean = NULL; 236 double total_mean; 237 double bcv, max_bcv; 238 int first_bin, last_bin; 239 int threshold, i; 240 SANE_Status ret = SANE_STATUS_NO_MEM; 241 242 DBG (10, "sanei_ir_threshold_otsu\n"); 243 244 cnh = sanei_ir_accumulate_norm_histo (norm_histo); 245 mean = malloc (HISTOGRAM_SIZE * sizeof (double)); 246 if (!cnh || !mean) 247 { 248 DBG (5, "sanei_ir_threshold_otsu: no buffers\n"); 249 goto cleanup; 250 } 251 252 mean[0] = 0.0; 253 for (i = 1; i < HISTOGRAM_SIZE; i++) 254 mean[i] = mean[i - 1] + i * norm_histo[i]; 255 total_mean = mean[HISTOGRAM_SIZE - 1]; 256 257 first_bin = 0; 258 for (i = 0; i < HISTOGRAM_SIZE; i++) 259 if (cnh[i] != 0) 260 { 261 first_bin = i; 262 break; 263 } 264 last_bin = HISTOGRAM_SIZE - 1; 265 for (i = HISTOGRAM_SIZE - 1; i >= first_bin; i--) 266 if (1.0 - cnh[i] != 0) 267 { 268 last_bin = i; 269 break; 270 } 271 272 threshold = INT_MIN; 273 max_bcv = 0.0; 274 for (i = first_bin; i <= last_bin; i++) 275 { 276 bcv = total_mean * cnh[i] - mean[i]; 277 bcv *= bcv / (cnh[i] * (1.0 - cnh[i])); 278 if (max_bcv < bcv) 279 { 280 max_bcv = bcv; 281 threshold = i; 282 } 283 } 284 285 if (threshold == INT_MIN) 286 { 287 DBG (5, "sanei_ir_threshold_otsu: no threshold found\n"); 288 ret = SANE_STATUS_INVAL; 289 } 290 else 291 { 292 ret = SANE_STATUS_GOOD; 293 if (params->depth > 8) 294 { 295 i = 1 << (params->depth - HISTOGRAM_SHIFT); 296 *thresh = threshold * i + i / 2; 297 } 298 else 299 *thresh = threshold; 300 DBG (10, "sanei_ir_threshold_otsu: threshold %d\n", *thresh); 301 } 302 cleanup: 303 if (cnh) 304 free (cnh); 305 if (mean) 306 free (mean); 307 return ret; 308} 309 310 311/* Implements a Maximum Entropy thresholding method 312 */ 313SANE_Status 314sanei_ir_threshold_maxentropy (const SANE_Parameters * params, 315 double * norm_histo, int *thresh) 316{ 317 int ih, it; 318 int threshold; 319 int first_bin; 320 int last_bin; 321 double tot_ent, max_ent; /* entropies */ 322 double ent_back, ent_obj; 323 double *P1; /* cumulative normalized histogram */ 324 double *P2; 325 SANE_Status ret = SANE_STATUS_NO_MEM; 326 327 DBG (10, "sanei_ir_threshold_maxentropy\n"); 328 329 /* Calculate the cumulative normalized histogram */ 330 P1 = sanei_ir_accumulate_norm_histo (norm_histo); 331 P2 = malloc (HISTOGRAM_SIZE * sizeof (double)); 332 if (!P1 || !P2) 333 { 334 DBG (5, "sanei_ir_threshold_maxentropy: no buffers\n"); 335 goto cleanup; 336 } 337 338 for ( ih = 0; ih < HISTOGRAM_SIZE; ih++ ) 339 P2[ih] = 1.0 - P1[ih]; 340 341 first_bin = 0; 342 for ( ih = 0; ih < HISTOGRAM_SIZE; ih++ ) 343 if (P1[ih] != 0) 344 { 345 first_bin = ih; 346 break; 347 } 348 last_bin = HISTOGRAM_SIZE - 1; 349 for ( ih = HISTOGRAM_SIZE - 1; ih >= first_bin; ih-- ) 350 if (P2[ih] != 0) 351 { 352 last_bin = ih; 353 break; 354 } 355 356 /* Calculate the total entropy each gray-level 357 * and find the threshold that maximizes it 358 */ 359 threshold = INT_MIN; 360 max_ent = DBL_MIN; 361 for ( it = first_bin; it <= last_bin; it++ ) 362 { 363 /* Entropy of the background pixels */ 364 ent_back = 0.0; 365 for ( ih = 0; ih <= it; ih++ ) 366 if (norm_histo[ih] != 0) 367 ent_back -= ( norm_histo[ih] / P1[it] ) * log ( norm_histo[ih] / P1[it] ); 368 369 /* Entropy of the object pixels */ 370 ent_obj = 0.0; 371 for ( ih = it + 1; ih < HISTOGRAM_SIZE; ih++ ) 372 if (norm_histo[ih] != 0) 373 ent_obj -= ( norm_histo[ih] / P2[it] ) * log ( norm_histo[ih] / P2[it] ); 374 375 /* Total entropy */ 376 tot_ent = ent_back + ent_obj; 377 378 if ( max_ent < tot_ent ) 379 { 380 max_ent = tot_ent; 381 threshold = it; 382 } 383 } 384 385 if (threshold == INT_MIN) 386 { 387 DBG (5, "sanei_ir_threshold_maxentropy: no threshold found\n"); 388 ret = SANE_STATUS_INVAL; 389 } 390 else 391 { 392 ret = SANE_STATUS_GOOD; 393 if (params->depth > 8) 394 { 395 it = 1 << (params->depth - HISTOGRAM_SHIFT); 396 *thresh = threshold * it + it / 2; 397 } 398 else 399 *thresh = threshold; 400 DBG (10, "sanei_ir_threshold_maxentropy: threshold %d\n", *thresh); 401 } 402 403 cleanup: 404 if (P1) 405 free (P1); 406 if (P2) 407 free (P2); 408 return ret; 409} 410 411/* Generate gray scale luminance image from separate R, G, B images 412 */ 413SANE_Status 414sanei_ir_RGB_luminance (SANE_Parameters * params, const SANE_Uint **in_img, 415 SANE_Uint **out_img) 416{ 417 SANE_Uint *outi; 418 int itop, i; 419 420 if ((params->depth < 8) || (params->depth > 16) || 421 (params->format != SANE_FRAME_GRAY)) 422 { 423 DBG (5, "sanei_ir_RGB_luminance: invalid format\n"); 424 return SANE_STATUS_UNSUPPORTED; 425 } 426 427 itop = params->pixels_per_line * params->lines; 428 outi = malloc (itop * sizeof(SANE_Uint)); 429 if (!outi) 430 { 431 DBG (5, "sanei_ir_RGB_luminance: can not allocate out_img\n"); 432 return SANE_STATUS_NO_MEM; 433 } 434 435 for (i = itop; i > 0; i--) 436 *(outi++) = (218 * (int) *(in_img[0]++) + 437 732 * (int) *(in_img[1]++) + 438 74 * (int) *(in_img[2]++)) >> 10; 439 *out_img = outi; 440 return SANE_STATUS_GOOD; 441} 442 443/* Convert image from >8 bit depth to an 8 bit image 444 */ 445SANE_Status 446sanei_ir_to_8bit (SANE_Parameters * params, const SANE_Uint *in_img, 447 SANE_Parameters * out_params, SANE_Uint **out_img) 448{ 449 SANE_Uint *outi; 450 size_t ssize; 451 int i, is; 452 453 if ((params->depth < 8) || (params->depth > 16)) 454 { 455 DBG (5, "sanei_ir_to_8bit: invalid format\n"); 456 return SANE_STATUS_UNSUPPORTED; 457 } 458 ssize = params->pixels_per_line * params->lines; 459 if (params->format == SANE_FRAME_RGB) 460 ssize *= 3; 461 outi = malloc (ssize * sizeof(SANE_Uint)); 462 if (!outi) 463 { 464 DBG (5, "sanei_ir_to_8bit: can not allocate out_img\n"); 465 return SANE_STATUS_NO_MEM; 466 } 467 468 if (out_params) 469 { 470 memmove (out_params, params, sizeof(SANE_Parameters)); 471 out_params->bytes_per_line = out_params->pixels_per_line; 472 if (params->format == SANE_FRAME_RGB) 473 out_params->bytes_per_line *= 3; 474 out_params->depth = 8; 475 } 476 477 memmove (outi, in_img, ssize * sizeof(SANE_Uint)); 478 is = params->depth - 8; 479 for (i = ssize; i > 0; i--) { 480 *outi = *outi >> is, outi += 2; 481 } 482 483 *out_img = outi; 484 return SANE_STATUS_GOOD; 485} 486 487/* allocate and initialize logarithmic lookup table 488 */ 489SANE_Status 490sanei_ir_ln_table (int len, double **lut_ln) 491{ 492 double *llut; 493 int i; 494 495 DBG (10, "sanei_ir_ln_table\n"); 496 497 llut = malloc (len * sizeof (double)); 498 if (!llut) 499 { 500 DBG (5, "sanei_ir_ln_table: no table\n"); 501 return SANE_STATUS_NO_MEM; 502 } 503 llut[0] = 0; 504 llut[1] = 0; 505 for (i = 2; i < len; i++) 506 llut[i] = log ((double) i); 507 508 *lut_ln = llut; 509 return SANE_STATUS_GOOD; 510} 511 512 513/* Reduce red spectral overlap from an infrared image plane 514 */ 515SANE_Status 516sanei_ir_spectral_clean (const SANE_Parameters * params, double *lut_ln, 517 const SANE_Uint *red_data, 518 SANE_Uint *ir_data) 519{ 520 const SANE_Uint *rptr; 521 SANE_Uint *iptr; 522 SANE_Int depth; 523 double *llut; 524 double rval, rsum, rrsum; 525 double risum, rfac, radd; 526 double *norm_histo; 527 int64_t isum; 528 int *calc_buf, *calc_ptr; 529 int ival, imin, imax; 530 int itop, len, ssize; 531 int thresh_low, thresh; 532 int irand, i; 533 SANE_Status status; 534 535 DBG (10, "sanei_ir_spectral_clean\n"); 536 537 itop = params->pixels_per_line * params->lines; 538 calc_buf = malloc (itop * sizeof (int)); /* could save this */ 539 if (!calc_buf) 540 { 541 DBG (5, "sanei_ir_spectral_clean: no buffer\n"); 542 return SANE_STATUS_NO_MEM; 543 } 544 545 depth = params->depth; 546 len = 1 << depth; 547 if (lut_ln) 548 llut = lut_ln; 549 else 550 { 551 status = sanei_ir_ln_table (len, &llut); 552 if (status != SANE_STATUS_GOOD) { 553 free (calc_buf); 554 return status; 555 } 556 } 557 558 /* determine not transparent areas to exclude them later 559 * TODO: this has not been tested for negatives 560 */ 561 thresh_low = INT_MAX; 562 status = 563 sanei_ir_create_norm_histogram (params, ir_data, &norm_histo); 564 if (status != SANE_STATUS_GOOD) 565 { 566 DBG (5, "sanei_ir_spectral_clean: no buffer\n"); 567 free (calc_buf); 568 return SANE_STATUS_NO_MEM; 569 } 570 571 /* TODO: remember only needed if cropping is not ok */ 572 status = sanei_ir_threshold_maxentropy (params, norm_histo, &thresh); 573 if (status == SANE_STATUS_GOOD) 574 thresh_low = thresh; 575 status = sanei_ir_threshold_otsu (params, norm_histo, &thresh); 576 if ((status == SANE_STATUS_GOOD) && (thresh < thresh_low)) 577 thresh_low = thresh; 578 status = sanei_ir_threshold_yen (params, norm_histo, &thresh); 579 if ((status == SANE_STATUS_GOOD) && (thresh < thresh_low)) 580 thresh_low = thresh; 581 if (thresh_low == INT_MAX) 582 thresh_low = 0; 583 else 584 thresh_low /= 2; 585 DBG (10, "sanei_ir_spectral_clean: low threshold %d\n", thresh_low); 586 587 /* calculate linear regression ired (red) from randomly chosen points */ 588 ssize = itop / 2; 589 if (SAMPLE_SIZE < ssize) 590 ssize = SAMPLE_SIZE; 591 isum = 0; 592 rsum = rrsum = risum = 0.0; 593 i = ssize; 594 while (i > 0) 595 { 596 irand = rand () % itop; 597 rval = llut[red_data[irand]]; 598 ival = ir_data[irand]; 599 if (ival > thresh_low) 600 { 601 isum += ival; 602 rsum += rval; 603 rrsum += rval * rval; 604 risum += rval * (double) ival; 605 i--; 606 } 607 } 608 609 /* "a" in ired = b + a * ln (red) */ 610 rfac = 611 ((double) ssize * risum - 612 rsum * (double) isum) / ((double) ssize * rrsum - rsum * rsum); 613 radd = ((double) isum - rfac * rsum) / (double) ssize; /* "b" unused */ 614 615 DBG (10, "sanei_ir_spectral_clean: n = %d, ired(red) = %f * ln(red) + %f\n", 616 ssize, rfac, radd); 617 618 /* now calculate ired' = ired - a * ln (red) */ 619 imin = INT_MAX; 620 imax = INT_MIN; 621 rptr = red_data; 622 iptr = ir_data; 623 calc_ptr = calc_buf; 624 for (i = itop; i > 0; i--) 625 { 626 ival = *iptr++ - (int) (rfac * llut[*rptr++] + 0.5); 627 if (ival > imax) 628 imax = ival; 629 if (ival < imin) 630 imin = ival; 631 *calc_ptr++ = ival; 632 } 633 634 /* scale the result back into the ired image */ 635 calc_ptr = calc_buf; 636 iptr = ir_data; 637 rfac = (double) (len - 1) / (double) (imax - imin); 638 for (i = itop; i > 0; i--) 639 *iptr++ = (double) (*calc_ptr++ - imin) * rfac; 640 641 if (!lut_ln) 642 free (llut); 643 free (calc_buf); 644 free (norm_histo); 645 return SANE_STATUS_GOOD; 646} 647 648 649/* Hopefully fast mean filter 650 * JV: what does this do? Remove local mean? 651 */ 652SANE_Status 653sanei_ir_filter_mean (const SANE_Parameters * params, 654 const SANE_Uint *in_img, SANE_Uint *out_img, 655 int win_rows, int win_cols) 656{ 657 const SANE_Uint *src; 658 SANE_Uint *dest; 659 int num_cols, num_rows; 660 int itop, iadd, isub; 661 int ndiv, the_sum; 662 int nrow, ncol; 663 int hwr, hwc; 664 int *sum; 665 int i, j; 666 667 DBG (10, "sanei_ir_filter_mean, window: %d x%d\n", win_rows, win_cols); 668 669 if (((win_rows & 1) == 0) || ((win_cols & 1) == 0)) 670 { 671 DBG (5, "sanei_ir_filter_mean: window even sized\n"); 672 return SANE_STATUS_INVAL; 673 } 674 675 num_cols = params->pixels_per_line; 676 num_rows = params->lines; 677 678 sum = malloc (num_cols * sizeof (int)); 679 if (!sum) 680 { 681 DBG (5, "sanei_ir_filter_mean: no buffer for sums\n"); 682 return SANE_STATUS_NO_MEM; 683 } 684 dest = out_img; 685 686 hwr = win_rows / 2; /* half window sizes */ 687 hwc = win_cols / 2; 688 689 /* pre-pre calculation */ 690 for (j = 0; j < num_cols; j++) 691 { 692 sum[j] = 0; 693 src = in_img + j; 694 for (i = 0; i < hwr; i++) 695 { 696 sum[j] += *src; 697 src += num_cols; 698 } 699 } 700 701 itop = num_rows * num_cols; 702 iadd = hwr * num_cols; 703 isub = (hwr - win_rows) * num_cols; 704 nrow = hwr; 705 706 for (i = 0; i < num_rows; i++) 707 { 708 /* update row sums if possible */ 709 if (isub >= 0) /* subtract old row */ 710 { 711 nrow--; 712 src = in_img + isub; 713 for (j = 0; j < num_cols; j++) 714 sum[j] -= *src++; 715 } 716 isub += num_cols; 717 718 if (iadd < itop) /* add new row */ 719 { 720 nrow++; 721 src = in_img + iadd; 722 for (j = 0; j < num_cols; j++) 723 sum[j] += *src++; 724 } 725 iadd += num_cols; 726 727 /* now we do the image columns using only the precalculated sums */ 728 729 the_sum = 0; /* precalculation */ 730 for (j = 0; j < hwc; j++) 731 the_sum += sum[j]; 732 ncol = hwc; 733 734 /* at the left margin, real index hwc lower */ 735 for (j = hwc; j < win_cols; j++) 736 { 737 ncol++; 738 the_sum += sum[j]; 739 *dest++ = the_sum / (ncol * nrow); 740 } 741 742 ndiv = ncol * nrow; 743 /* in the middle, real index hwc + 1 higher */ 744 for (j = 0; j < num_cols - win_cols; j++) 745 { 746 the_sum -= sum[j]; 747 the_sum += sum[j + win_cols]; 748 *dest++ = the_sum / ndiv; 749 } 750 751 /* at the right margin, real index hwc + 1 higher */ 752 for (j = num_cols - win_cols; j < num_cols - hwc - 1; j++) 753 { 754 ncol--; 755 the_sum -= sum[j]; /* j - hwc - 1 */ 756 *dest++ = the_sum / (ncol * nrow); 757 } 758 } 759 free (sum); 760 return SANE_STATUS_GOOD; 761} 762 763 764/* Find noise by adaptive thresholding 765 */ 766SANE_Status 767sanei_ir_filter_madmean (const SANE_Parameters * params, 768 const SANE_Uint *in_img, 769 SANE_Uint ** out_img, int win_size, 770 int a_val, int b_val) 771{ 772 SANE_Uint *delta_ij, *delta_ptr; 773 SANE_Uint *mad_ij; 774 const SANE_Uint *mad_ptr; 775 SANE_Uint *out_ij, *dest8; 776 double ab_term; 777 int num_rows, num_cols; 778 int threshold, itop; 779 size_t size; 780 int ival, i; 781 int depth; 782 SANE_Status ret = SANE_STATUS_NO_MEM; 783 784 DBG (10, "sanei_ir_filter_madmean\n"); 785 786 depth = params->depth; 787 if (depth != 8) 788 { 789 a_val = a_val << (depth - 8); 790 b_val = b_val << (depth - 8); 791 } 792 num_cols = params->pixels_per_line; 793 num_rows = params->lines; 794 itop = num_rows * num_cols; 795 size = itop * sizeof (SANE_Uint); 796 out_ij = malloc (size); 797 delta_ij = malloc (size); 798 mad_ij = malloc (size); 799 800 if (out_ij && delta_ij && mad_ij) 801 { 802 /* get the differences to the local mean */ 803 mad_ptr = in_img; 804 if (sanei_ir_filter_mean (params, mad_ptr, delta_ij, win_size, win_size) 805 == SANE_STATUS_GOOD) 806 { 807 delta_ptr = delta_ij; 808 for (i = 0; i < itop; i++) 809 { 810 ival = *mad_ptr++ - *delta_ptr; 811 *delta_ptr++ = abs (ival); 812 } 813 /* make the second filtering window a bit larger */ 814 win_size = MAD_WIN2_SIZE(win_size); 815 /* and get the local mean differences */ 816 if (sanei_ir_filter_mean 817 (params, delta_ij, mad_ij, win_size, 818 win_size) == SANE_STATUS_GOOD) 819 { 820 mad_ptr = mad_ij; 821 delta_ptr = delta_ij; 822 dest8 = out_ij; 823 /* construct the noise map */ 824 ab_term = (b_val - a_val) / (double) b_val; 825 for (i = 0; i < itop; i++) 826 { 827 /* by calculating the threshold */ 828 ival = *mad_ptr++; 829 if (ival >= b_val) /* outlier */ 830 threshold = a_val; 831 else 832 threshold = a_val + (double) ival *ab_term; 833 /* above threshold is noise, indicated by 0 */ 834 if (*delta_ptr++ >= threshold) 835 *dest8++ = 0; 836 else 837 *dest8++ = 255; 838 } 839 *out_img = out_ij; 840 ret = SANE_STATUS_GOOD; 841 } 842 } 843 } 844 else 845 DBG (5, "sanei_ir_filter_madmean: Cannot allocate buffers\n"); 846 847 free (mad_ij); 848 free (delta_ij); 849 return ret; 850} 851 852 853/* Add dark pixels to mask from static threshold 854 */ 855void 856sanei_ir_add_threshold (const SANE_Parameters * params, 857 const SANE_Uint *in_img, 858 SANE_Uint * mask_img, int threshold) 859{ 860 const SANE_Uint *in_ptr; 861 SANE_Uint *mask_ptr; 862 int itop, i; 863 864 DBG (10, "sanei_ir_add_threshold\n"); 865 866 itop = params->pixels_per_line * params->lines; 867 in_ptr = in_img; 868 mask_ptr = mask_img; 869 870 for (i = itop; i > 0; i--) 871 { 872 if (*in_ptr++ <= threshold) 873 *mask_ptr = 0; 874 mask_ptr++; 875 } 876} 877 878 879/* Calculate minimal Manhattan distances for an image mask 880 */ 881void 882sanei_ir_manhattan_dist (const SANE_Parameters * params, 883 const SANE_Uint * mask_img, unsigned int *dist_map, 884 unsigned int *idx_map, unsigned int erode) 885{ 886 const SANE_Uint *mask; 887 unsigned int *index, *manhattan; 888 int rows, cols, itop; 889 int i, j; 890 891 DBG (10, "sanei_ir_manhattan_dist\n"); 892 893 if (erode != 0) 894 erode = 255; 895 896 /* initialize maps */ 897 cols = params->pixels_per_line; 898 rows = params->lines; 899 itop = rows * cols; 900 mask = mask_img; 901 manhattan = dist_map; 902 index = idx_map; 903 for (i = 0; i < itop; i++) 904 { 905 *manhattan++ = *mask++; 906 *index++ = i; 907 } 908 909 /* traverse from top left to bottom right */ 910 manhattan = dist_map; 911 index = idx_map; 912 for (i = 0; i < rows; i++) 913 for (j = 0; j < cols; j++) 914 { 915 if (*manhattan == erode) 916 { 917 /* take original, distance = 0, index stays the same */ 918 *manhattan = 0; 919 } 920 else 921 { 922 /* assume maximal distance to clean pixel */ 923 *manhattan = cols + rows; 924 /* or one further away than pixel to the top */ 925 if (i > 0) 926 if (manhattan[-cols] + 1 < *manhattan) 927 { 928 *manhattan = manhattan[-cols] + 1; 929 *index = index[-cols]; /* index follows */ 930 } 931 /* or one further away than pixel to the left */ 932 if (j > 0) 933 { 934 if (manhattan[-1] + 1 < *manhattan) 935 { 936 *manhattan = manhattan[-1] + 1; 937 *index = index[-1]; /* index follows */ 938 } 939 if (manhattan[-1] + 1 == *manhattan) 940 if (rand () % 2 == 0) /* chose index */ 941 *index = index[-1]; 942 } 943 } 944 manhattan++; 945 index++; 946 } 947 948 /* traverse from bottom right to top left */ 949 manhattan = dist_map + itop - 1; 950 index = idx_map + itop - 1; 951 for (i = rows - 1; i >= 0; i--) 952 for (j = cols - 1; j >= 0; j--) 953 { 954 if (i < rows - 1) 955 { 956 /* either what we had on the first pass 957 or one more than the pixel to the bottm */ 958 if (manhattan[+cols] + 1 < *manhattan) 959 { 960 *manhattan = manhattan[+cols] + 1; 961 *index = index[+cols]; /* index follows */ 962 } 963 if (manhattan[+cols] + 1 == *manhattan) 964 if (rand () % 2 == 0) /* chose index */ 965 *index = index[+cols]; 966 } 967 if (j < cols - 1) 968 { 969 /* or one more than pixel to the right */ 970 if (manhattan[1] + 1 < *manhattan) 971 { 972 *manhattan = manhattan[1] + 1; 973 *index = index[1]; /* index follows */ 974 } 975 if (manhattan[1] + 1 == *manhattan) 976 if (rand () % 2 == 0) /* chose index */ 977 *index = index[1]; 978 } 979 manhattan--; 980 index--; 981 } 982} 983 984 985/* dilate or erode a mask image */ 986 987void 988sanei_ir_dilate (const SANE_Parameters *params, SANE_Uint *mask_img, 989 unsigned int *dist_map, unsigned int *idx_map, int by) 990{ 991 SANE_Uint *mask; 992 unsigned int *manhattan; 993 unsigned int erode; 994 unsigned int thresh; 995 int i, itop; 996 997 DBG (10, "sanei_ir_dilate\n"); 998 999 if (by == 0) 1000 return; 1001 if (by > 0) 1002 { 1003 erode = 0; 1004 thresh = by; 1005 } 1006 else 1007 { 1008 erode = 1; 1009 thresh = -by; 1010 } 1011 1012 itop = params->pixels_per_line * params->lines; 1013 mask = mask_img; 1014 sanei_ir_manhattan_dist (params, mask_img, dist_map, idx_map, erode); 1015 1016 manhattan = dist_map; 1017 for (i = 0; i < itop; i++) 1018 { 1019 if (*manhattan++ <= thresh) 1020 *mask++ = 0; 1021 else 1022 *mask++ = 255; 1023 } 1024 1025 return; 1026} 1027 1028 1029/* Suggest cropping for dark margins of positive film 1030 */ 1031void 1032sanei_ir_find_crop (const SANE_Parameters * params, 1033 unsigned int * dist_map, int inner, int * edges) 1034{ 1035 int width = params->pixels_per_line; 1036 int height = params->lines; 1037 uint64_t sum_x, sum_y, n; 1038 int64_t sum_xx, sum_xy; 1039 double a, b, mami; 1040 unsigned int *src; 1041 int off1, off2, inc, wh, i, j; 1042 1043 DBG (10, "sanei_ir_find_crop\n"); 1044 1045 /* loop through top, bottom, left, right */ 1046 for (j = 0; j < 4; j++) 1047 { 1048 if (j < 2) /* top, bottom */ 1049 { 1050 off1 = width / 8; /* only middle 3/4 */ 1051 off2 = width - off1; 1052 n = width - 2 * off1; 1053 src = dist_map + off1; /* first row */ 1054 inc = 1; 1055 wh = width; 1056 if (j == 1) /* last row */ 1057 src += (height - 1) * width; 1058 } 1059 else /* left, right */ 1060 { 1061 off1 = height / 8; /* only middle 3/4 */ 1062 off2 = height - off1; 1063 n = height - 2 * off1; 1064 src = dist_map + (off1 * width); /* first column */ 1065 inc = width; 1066 wh = height; 1067 if (j == 3) 1068 src += width - 1; /* last column */ 1069 } 1070 1071 /* calculate linear regression */ 1072 sum_x = 0; sum_y = 0; 1073 sum_xx = 0; sum_xy = 0; 1074 for (i = off1; i < off2; i++) 1075 { 1076 sum_x += i; 1077 sum_y += *src; 1078 sum_xx += i * i; 1079 sum_xy += i * (*src); 1080 src += inc; 1081 } 1082 b = ((double) n * (double) sum_xy - (double) sum_x * (double) sum_y) 1083 / ((double) n * (double) sum_xx - (double) sum_x * (double) sum_x); 1084 a = ((double) sum_y - b * (double) sum_x) / (double) n; 1085 1086 DBG (10, "sanei_ir_find_crop: y = %f + %f * x\n", a, b); 1087 1088 /* take maximal/minimal value from either side */ 1089 mami = a + b * (wh - 1); 1090 if (inner) 1091 { 1092 if (a > mami) 1093 mami = a; 1094 } 1095 else 1096 { 1097 if (a < mami) 1098 mami = a; 1099 } 1100 edges[j] = mami + 0.5; 1101 } 1102 edges[1] = height - edges[1]; 1103 edges[3] = width - edges[3]; 1104 1105 DBG (10, "sanei_ir_find_crop: would crop at top: %d, bot: %d, left %d, right %d\n", 1106 edges[0], edges[1], edges[2], edges[3]); 1107 1108 return; 1109} 1110 1111 1112/* Dilate clean image parts into dirty ones and smooth 1113 */ 1114SANE_Status 1115sanei_ir_dilate_mean (const SANE_Parameters * params, 1116 SANE_Uint **in_img, 1117 SANE_Uint * mask_img, 1118 int dist_max, int expand, int win_size, 1119 SANE_Bool smooth, int inner, 1120 int *crop) 1121{ 1122 SANE_Uint *color; 1123 SANE_Uint *plane; 1124 unsigned int *dist_map, *manhattan; 1125 unsigned int *idx_map, *index; 1126 int dist; 1127 int rows, cols; 1128 int k, i, itop; 1129 SANE_Status ret = SANE_STATUS_NO_MEM; 1130 1131 DBG (10, "sanei_ir_dilate_mean(): dist max = %d, expand = %d, win size = %d, smooth = %d, inner = %d\n", 1132 dist_max, expand, win_size, smooth, inner); 1133 1134 cols = params->pixels_per_line; 1135 rows = params->lines; 1136 itop = rows * cols; 1137 idx_map = malloc (itop * sizeof (unsigned int)); 1138 dist_map = malloc (itop * sizeof (unsigned int)); 1139 plane = malloc (itop * sizeof (SANE_Uint)); 1140 1141 if (!idx_map || !dist_map || !plane) 1142 DBG (5, "sanei_ir_dilate_mean: Cannot allocate buffers\n"); 1143 else 1144 { 1145 /* expand dirty regions into their half dirty surround*/ 1146 if (expand > 0) 1147 sanei_ir_dilate (params, mask_img, dist_map, idx_map, expand); 1148 /* for dirty pixels determine the closest clean ones */ 1149 sanei_ir_manhattan_dist (params, mask_img, dist_map, idx_map, 1); 1150 1151 /* use the distance map to find how to crop dark edges */ 1152 if (crop) 1153 sanei_ir_find_crop (params, dist_map, inner, crop); 1154 1155 /* replace dirty pixels */ 1156 for (k = 0; k < 3; k++) 1157 { 1158 manhattan = dist_map; 1159 index = idx_map; 1160 color = in_img[k]; 1161 /* first replacement */ 1162 for (i = 0; i < itop; i++) 1163 { 1164 dist = *manhattan++; 1165 if ((dist != 0) && (dist <= dist_max)) 1166 color[i] = color[index[i]]; 1167 } 1168 /* adapt pixels to their new surround and 1169 * smooth the whole image or the replaced pixels only */ 1170 ret = 1171 sanei_ir_filter_mean (params, color, plane, win_size, win_size); 1172 if (ret != SANE_STATUS_GOOD) 1173 break; 1174 else 1175 if (smooth) 1176 { 1177 /* a second mean results in triangular blur */ 1178 DBG (10, "sanei_ir_dilate_mean(): smoothing whole image\n"); 1179 ret = 1180 sanei_ir_filter_mean (params, plane, color, win_size, 1181 win_size); 1182 if (ret != SANE_STATUS_GOOD) 1183 break; 1184 } 1185 else 1186 { 1187 /* replace with smoothened pixels only */ 1188 DBG (10, "sanei_ir_dilate_mean(): smoothing replaced pixels only\n"); 1189 manhattan = dist_map; 1190 for (i = 0; i < itop; i++) 1191 { 1192 dist = *manhattan++; 1193 if ((dist != 0) && (dist <= dist_max)) 1194 color[i] = plane[i]; 1195 } 1196 } 1197 } 1198 } 1199 free (plane); 1200 free (dist_map); 1201 free (idx_map); 1202 1203 return ret; 1204} 1205