1import inspect 2import itertools 3import math 4from collections import OrderedDict 5from collections.abc import Iterable 6 7import numpy as np 8from scipy import ndimage as ndi 9 10from .._shared.filters import gaussian 11from .._shared.utils import _supported_float_type, deprecate_kwarg, warn 12from .._shared.version_requirements import require 13from ..exposure import histogram 14from ..filters._multiotsu import (_get_multiotsu_thresh_indices, 15 _get_multiotsu_thresh_indices_lut) 16from ..transform import integral_image 17from ..util import dtype_limits 18from ._sparse import _correlate_sparse, _validate_window_size 19 20__all__ = ['try_all_threshold', 21 'threshold_otsu', 22 'threshold_yen', 23 'threshold_isodata', 24 'threshold_li', 25 'threshold_local', 26 'threshold_minimum', 27 'threshold_mean', 28 'threshold_niblack', 29 'threshold_sauvola', 30 'threshold_triangle', 31 'apply_hysteresis_threshold', 32 'threshold_multiotsu'] 33 34 35def _try_all(image, methods=None, figsize=None, num_cols=2, verbose=True): 36 """Returns a figure comparing the outputs of different methods. 37 38 Parameters 39 ---------- 40 image : (N, M) ndarray 41 Input image. 42 methods : dict, optional 43 Names and associated functions. 44 Functions must take and return an image. 45 figsize : tuple, optional 46 Figure size (in inches). 47 num_cols : int, optional 48 Number of columns. 49 verbose : bool, optional 50 Print function name for each method. 51 52 Returns 53 ------- 54 fig, ax : tuple 55 Matplotlib figure and axes. 56 """ 57 from matplotlib import pyplot as plt 58 59 # Compute the image histogram for better performances 60 nbins = 256 # Default in threshold functions 61 hist = histogram(image.ravel(), nbins, source_range='image') 62 63 # Handle default value 64 methods = methods or {} 65 66 num_rows = math.ceil((len(methods) + 1.) / num_cols) 67 fig, ax = plt.subplots(num_rows, num_cols, figsize=figsize, 68 sharex=True, sharey=True) 69 ax = ax.ravel() 70 71 ax[0].imshow(image, cmap=plt.cm.gray) 72 ax[0].set_title('Original') 73 74 i = 1 75 for name, func in methods.items(): 76 # Use precomputed histogram for supporting functions 77 sig = inspect.signature(func) 78 _kwargs = dict(hist=hist) if 'hist' in sig.parameters else {} 79 80 ax[i].set_title(name) 81 try: 82 ax[i].imshow(func(image, **_kwargs), cmap=plt.cm.gray) 83 except Exception as e: 84 ax[i].text(0.5, 0.5, "%s" % type(e).__name__, 85 ha="center", va="center", transform=ax[i].transAxes) 86 i += 1 87 if verbose: 88 print(func.__orifunc__) 89 90 for a in ax: 91 a.axis('off') 92 93 fig.tight_layout() 94 return fig, ax 95 96 97@require("matplotlib", ">=3.0.3") 98def try_all_threshold(image, figsize=(8, 5), verbose=True): 99 """Returns a figure comparing the outputs of different thresholding methods. 100 101 Parameters 102 ---------- 103 image : (N, M) ndarray 104 Input image. 105 figsize : tuple, optional 106 Figure size (in inches). 107 verbose : bool, optional 108 Print function name for each method. 109 110 Returns 111 ------- 112 fig, ax : tuple 113 Matplotlib figure and axes. 114 115 Notes 116 ----- 117 The following algorithms are used: 118 119 * isodata 120 * li 121 * mean 122 * minimum 123 * otsu 124 * triangle 125 * yen 126 127 Examples 128 -------- 129 >>> from skimage.data import text 130 >>> fig, ax = try_all_threshold(text(), figsize=(10, 6), verbose=False) 131 """ 132 def thresh(func): 133 """ 134 A wrapper function to return a thresholded image. 135 """ 136 def wrapper(im): 137 return im > func(im) 138 try: 139 wrapper.__orifunc__ = func.__orifunc__ 140 except AttributeError: 141 wrapper.__orifunc__ = func.__module__ + '.' + func.__name__ 142 return wrapper 143 144 # Global algorithms. 145 methods = OrderedDict({'Isodata': thresh(threshold_isodata), 146 'Li': thresh(threshold_li), 147 'Mean': thresh(threshold_mean), 148 'Minimum': thresh(threshold_minimum), 149 'Otsu': thresh(threshold_otsu), 150 'Triangle': thresh(threshold_triangle), 151 'Yen': thresh(threshold_yen)}) 152 153 return _try_all(image, figsize=figsize, 154 methods=methods, verbose=verbose) 155 156 157def threshold_local(image, block_size=3, method='gaussian', offset=0, 158 mode='reflect', param=None, cval=0): 159 """Compute a threshold mask image based on local pixel neighborhood. 160 161 Also known as adaptive or dynamic thresholding. The threshold value is 162 the weighted mean for the local neighborhood of a pixel subtracted by a 163 constant. Alternatively the threshold can be determined dynamically by a 164 given function, using the 'generic' method. 165 166 Parameters 167 ---------- 168 image : (N, M[, ..., P]) ndarray 169 Grayscale input image. 170 block_size : int or sequence of int 171 Odd size of pixel neighborhood which is used to calculate the 172 threshold value (e.g. 3, 5, 7, ..., 21, ...). 173 method : {'generic', 'gaussian', 'mean', 'median'}, optional 174 Method used to determine adaptive threshold for local neighbourhood in 175 weighted mean image. 176 177 * 'generic': use custom function (see ``param`` parameter) 178 * 'gaussian': apply gaussian filter (see ``param`` parameter for custom\ 179 sigma value) 180 * 'mean': apply arithmetic mean filter 181 * 'median': apply median rank filter 182 183 By default the 'gaussian' method is used. 184 offset : float, optional 185 Constant subtracted from weighted mean of neighborhood to calculate 186 the local threshold value. Default offset is 0. 187 mode : {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}, optional 188 The mode parameter determines how the array borders are handled, where 189 cval is the value when mode is equal to 'constant'. 190 Default is 'reflect'. 191 param : {int, function}, optional 192 Either specify sigma for 'gaussian' method or function object for 193 'generic' method. This functions takes the flat array of local 194 neighbourhood as a single argument and returns the calculated 195 threshold for the centre pixel. 196 cval : float, optional 197 Value to fill past edges of input if mode is 'constant'. 198 199 Returns 200 ------- 201 threshold : (N, M[, ..., P]) ndarray 202 Threshold image. All pixels in the input image higher than the 203 corresponding pixel in the threshold image are considered foreground. 204 205 References 206 ---------- 207 .. [1] Gonzalez, R. C. and Wood, R. E. "Digital Image Processing 208 (2nd Edition)." Prentice-Hall Inc., 2002: 600--612. 209 ISBN: 0-201-18075-8 210 211 Examples 212 -------- 213 >>> from skimage.data import camera 214 >>> image = camera()[:50, :50] 215 >>> binary_image1 = image > threshold_local(image, 15, 'mean') 216 >>> func = lambda arr: arr.mean() 217 >>> binary_image2 = image > threshold_local(image, 15, 'generic', 218 ... param=func) 219 220 """ 221 222 if np.isscalar(block_size): 223 block_size = (block_size,) * image.ndim 224 elif len(block_size) != image.ndim: 225 raise ValueError("len(block_size) must equal image.ndim.") 226 block_size = tuple(block_size) 227 if any(b % 2 == 0 for b in block_size): 228 raise ValueError(f'block_size must be odd! Given block_size ' 229 f'{block_size} contains even values.') 230 float_dtype = _supported_float_type(image) 231 image = image.astype(float_dtype, copy=False) 232 thresh_image = np.zeros(image.shape, dtype=float_dtype) 233 if method == 'generic': 234 ndi.generic_filter(image, param, block_size, 235 output=thresh_image, mode=mode, cval=cval) 236 elif method == 'gaussian': 237 if param is None: 238 # automatically determine sigma which covers > 99% of distribution 239 sigma = tuple([(b - 1) / 6.0 for b in block_size]) 240 else: 241 sigma = param 242 gaussian(image, sigma, output=thresh_image, mode=mode, cval=cval) 243 elif method == 'mean': 244 ndi.uniform_filter(image, block_size, output=thresh_image, mode=mode, 245 cval=cval) 246 elif method == 'median': 247 ndi.median_filter(image, block_size, output=thresh_image, mode=mode, 248 cval=cval) 249 else: 250 raise ValueError("Invalid method specified. Please use `generic`, " 251 "`gaussian`, `mean`, or `median`.") 252 253 return thresh_image - offset 254 255 256def _validate_image_histogram(image, hist, nbins=None, normalize=False): 257 """Ensure that either image or hist were given, return valid histogram. 258 259 If hist is given, image is ignored. 260 261 Parameters 262 ---------- 263 image : array or None 264 Grayscale image. 265 hist : array, 2-tuple of array, or None 266 Histogram, either a 1D counts array, or an array of counts together 267 with an array of bin centers. 268 nbins : int, optional 269 The number of bins with which to compute the histogram, if `hist` is 270 None. 271 normalize : bool 272 If hist is not given, it will be computed by this function. This 273 parameter determines whether the computed histogram is normalized 274 (i.e. entries sum up to 1) or not. 275 276 Returns 277 ------- 278 counts : 1D array of float 279 Each element is the number of pixels falling in each intensity bin. 280 bin_centers : 1D array 281 Each element is the value corresponding to the center of each intensity 282 bin. 283 284 Raises 285 ------ 286 ValueError : if image and hist are both None 287 """ 288 if image is None and hist is None: 289 raise Exception("Either image or hist must be provided.") 290 291 if hist is not None: 292 if isinstance(hist, (tuple, list)): 293 counts, bin_centers = hist 294 else: 295 counts = hist 296 bin_centers = np.arange(counts.size) 297 298 if counts[0] == 0 or counts[-1] == 0: 299 # Trim histogram from both ends by removing starting and 300 # ending zeroes as in histogram(..., source_range="image") 301 cond = counts > 0 302 start = np.argmax(cond) 303 end = cond.size - np.argmax(cond[::-1]) 304 counts, bin_centers = counts[start:end], bin_centers[start:end] 305 else: 306 counts, bin_centers = histogram( 307 image.ravel(), nbins, source_range='image', normalize=normalize 308 ) 309 return counts.astype(float), bin_centers 310 311 312def threshold_otsu(image=None, nbins=256, *, hist=None): 313 """Return threshold value based on Otsu's method. 314 315 Either image or hist must be provided. If hist is provided, the actual 316 histogram of the image is ignored. 317 318 Parameters 319 ---------- 320 image : (N, M[, ..., P]) ndarray, optional 321 Grayscale input image. 322 nbins : int, optional 323 Number of bins used to calculate histogram. This value is ignored for 324 integer arrays. 325 hist : array, or 2-tuple of arrays, optional 326 Histogram from which to determine the threshold, and optionally a 327 corresponding array of bin center intensities. If no hist provided, 328 this function will compute it from the image. 329 330 331 Returns 332 ------- 333 threshold : float 334 Upper threshold value. All pixels with an intensity higher than 335 this value are assumed to be foreground. 336 337 References 338 ---------- 339 .. [1] Wikipedia, https://en.wikipedia.org/wiki/Otsu's_Method 340 341 Examples 342 -------- 343 >>> from skimage.data import camera 344 >>> image = camera() 345 >>> thresh = threshold_otsu(image) 346 >>> binary = image <= thresh 347 348 Notes 349 ----- 350 The input image must be grayscale. 351 """ 352 if image is not None and image.ndim > 2 and image.shape[-1] in (3, 4): 353 warn(f'threshold_otsu is expected to work correctly only for ' 354 f'grayscale images; image shape {image.shape} looks like ' 355 f'that of an RGB image.') 356 357 # Check if the image has more than one intensity value; if not, return that 358 # value 359 if image is not None: 360 first_pixel = image.ravel()[0] 361 if np.all(image == first_pixel): 362 return first_pixel 363 364 counts, bin_centers = _validate_image_histogram(image, hist, nbins) 365 366 # class probabilities for all possible thresholds 367 weight1 = np.cumsum(counts) 368 weight2 = np.cumsum(counts[::-1])[::-1] 369 # class means for all possible thresholds 370 mean1 = np.cumsum(counts * bin_centers) / weight1 371 mean2 = (np.cumsum((counts * bin_centers)[::-1]) / weight2[::-1])[::-1] 372 373 # Clip ends to align class 1 and class 2 variables: 374 # The last value of ``weight1``/``mean1`` should pair with zero values in 375 # ``weight2``/``mean2``, which do not exist. 376 variance12 = weight1[:-1] * weight2[1:] * (mean1[:-1] - mean2[1:]) ** 2 377 378 idx = np.argmax(variance12) 379 threshold = bin_centers[idx] 380 381 return threshold 382 383 384def threshold_yen(image=None, nbins=256, *, hist=None): 385 """Return threshold value based on Yen's method. 386 Either image or hist must be provided. In case hist is given, the actual 387 histogram of the image is ignored. 388 389 Parameters 390 ---------- 391 image : (N, M[, ..., P]) ndarray 392 Grayscale input image. 393 nbins : int, optional 394 Number of bins used to calculate histogram. This value is ignored for 395 integer arrays. 396 hist : array, or 2-tuple of arrays, optional 397 Histogram from which to determine the threshold, and optionally a 398 corresponding array of bin center intensities. 399 An alternative use of this function is to pass it only hist. 400 401 Returns 402 ------- 403 threshold : float 404 Upper threshold value. All pixels with an intensity higher than 405 this value are assumed to be foreground. 406 407 References 408 ---------- 409 .. [1] Yen J.C., Chang F.J., and Chang S. (1995) "A New Criterion 410 for Automatic Multilevel Thresholding" IEEE Trans. on Image 411 Processing, 4(3): 370-378. :DOI:`10.1109/83.366472` 412 .. [2] Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding 413 Techniques and Quantitative Performance Evaluation" Journal of 414 Electronic Imaging, 13(1): 146-165, :DOI:`10.1117/1.1631315` 415 http://www.busim.ee.boun.edu.tr/~sankur/SankurFolder/Threshold_survey.pdf 416 .. [3] ImageJ AutoThresholder code, http://fiji.sc/wiki/index.php/Auto_Threshold 417 418 Examples 419 -------- 420 >>> from skimage.data import camera 421 >>> image = camera() 422 >>> thresh = threshold_yen(image) 423 >>> binary = image <= thresh 424 """ 425 counts, bin_centers = _validate_image_histogram(image, hist, nbins) 426 427 # On blank images (e.g. filled with 0) with int dtype, `histogram()` 428 # returns ``bin_centers`` containing only one value. Speed up with it. 429 if bin_centers.size == 1: 430 return bin_centers[0] 431 432 # Calculate probability mass function 433 pmf = counts.astype(np.float32) / counts.sum() 434 P1 = np.cumsum(pmf) # Cumulative normalized histogram 435 P1_sq = np.cumsum(pmf ** 2) 436 # Get cumsum calculated from end of squared array: 437 P2_sq = np.cumsum(pmf[::-1] ** 2)[::-1] 438 # P2_sq indexes is shifted +1. I assume, with P1[:-1] it's help avoid 439 # '-inf' in crit. ImageJ Yen implementation replaces those values by zero. 440 crit = np.log(((P1_sq[:-1] * P2_sq[1:]) ** -1) * 441 (P1[:-1] * (1.0 - P1[:-1])) ** 2) 442 return bin_centers[crit.argmax()] 443 444 445def threshold_isodata(image=None, nbins=256, return_all=False, *, hist=None): 446 """Return threshold value(s) based on ISODATA method. 447 448 Histogram-based threshold, known as Ridler-Calvard method or inter-means. 449 Threshold values returned satisfy the following equality:: 450 451 threshold = (image[image <= threshold].mean() + 452 image[image > threshold].mean()) / 2.0 453 454 That is, returned thresholds are intensities that separate the image into 455 two groups of pixels, where the threshold intensity is midway between the 456 mean intensities of these groups. 457 458 For integer images, the above equality holds to within one; for floating- 459 point images, the equality holds to within the histogram bin-width. 460 461 Either image or hist must be provided. In case hist is given, the actual 462 histogram of the image is ignored. 463 464 Parameters 465 ---------- 466 image : (N, M[, ..., P]) ndarray 467 Grayscale input image. 468 nbins : int, optional 469 Number of bins used to calculate histogram. This value is ignored for 470 integer arrays. 471 return_all : bool, optional 472 If False (default), return only the lowest threshold that satisfies 473 the above equality. If True, return all valid thresholds. 474 hist : array, or 2-tuple of arrays, optional 475 Histogram to determine the threshold from and a corresponding array 476 of bin center intensities. Alternatively, only the histogram can be 477 passed. 478 479 Returns 480 ------- 481 threshold : float or int or array 482 Threshold value(s). 483 484 References 485 ---------- 486 .. [1] Ridler, TW & Calvard, S (1978), "Picture thresholding using an 487 iterative selection method" 488 IEEE Transactions on Systems, Man and Cybernetics 8: 630-632, 489 :DOI:`10.1109/TSMC.1978.4310039` 490 .. [2] Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding 491 Techniques and Quantitative Performance Evaluation" Journal of 492 Electronic Imaging, 13(1): 146-165, 493 http://www.busim.ee.boun.edu.tr/~sankur/SankurFolder/Threshold_survey.pdf 494 :DOI:`10.1117/1.1631315` 495 .. [3] ImageJ AutoThresholder code, 496 http://fiji.sc/wiki/index.php/Auto_Threshold 497 498 Examples 499 -------- 500 >>> from skimage.data import coins 501 >>> image = coins() 502 >>> thresh = threshold_isodata(image) 503 >>> binary = image > thresh 504 """ 505 counts, bin_centers = _validate_image_histogram(image, hist, nbins) 506 507 # image only contains one unique value 508 if len(bin_centers) == 1: 509 if return_all: 510 return bin_centers 511 else: 512 return bin_centers[0] 513 514 counts = counts.astype(np.float32) 515 516 # csuml and csumh contain the count of pixels in that bin or lower, and 517 # in all bins strictly higher than that bin, respectively 518 csuml = np.cumsum(counts) 519 csumh = csuml[-1] - csuml 520 521 # intensity_sum contains the total pixel intensity from each bin 522 intensity_sum = counts * bin_centers 523 524 # l and h contain average value of all pixels in that bin or lower, and 525 # in all bins strictly higher than that bin, respectively. 526 # Note that since exp.histogram does not include empty bins at the low or 527 # high end of the range, csuml and csumh are strictly > 0, except in the 528 # last bin of csumh, which is zero by construction. 529 # So no worries about division by zero in the following lines, except 530 # for the last bin, but we can ignore that because no valid threshold 531 # can be in the top bin. 532 # To avoid the division by zero, we simply skip over the last element in 533 # all future computation. 534 csum_intensity = np.cumsum(intensity_sum) 535 lower = csum_intensity[:-1] / csuml[:-1] 536 higher = (csum_intensity[-1] - csum_intensity[:-1]) / csumh[:-1] 537 538 # isodata finds threshold values that meet the criterion t = (l + m)/2 539 # where l is the mean of all pixels <= t and h is the mean of all pixels 540 # > t, as calculated above. So we are looking for places where 541 # (l + m) / 2 equals the intensity value for which those l and m figures 542 # were calculated -- which is, of course, the histogram bin centers. 543 # We only require this equality to be within the precision of the bin 544 # width, of course. 545 all_mean = (lower + higher) / 2.0 546 bin_width = bin_centers[1] - bin_centers[0] 547 548 # Look only at thresholds that are below the actual all_mean value, 549 # for consistency with the threshold being included in the lower pixel 550 # group. Otherwise can get thresholds that are not actually fixed-points 551 # of the isodata algorithm. For float images, this matters less, since 552 # there really can't be any guarantees anymore anyway. 553 distances = all_mean - bin_centers[:-1] 554 thresholds = bin_centers[:-1][(distances >= 0) & (distances < bin_width)] 555 556 if return_all: 557 return thresholds 558 else: 559 return thresholds[0] 560 561 562# Computing a histogram using np.histogram on a uint8 image with bins=256 563# doesn't work and results in aliasing problems. We use a fully specified set 564# of bins to ensure that each uint8 value false into its own bin. 565_DEFAULT_ENTROPY_BINS = tuple(np.arange(-0.5, 255.51, 1)) 566 567 568def _cross_entropy(image, threshold, bins=_DEFAULT_ENTROPY_BINS): 569 """Compute cross-entropy between distributions above and below a threshold. 570 571 Parameters 572 ---------- 573 image : array 574 The input array of values. 575 threshold : float 576 The value dividing the foreground and background in ``image``. 577 bins : int or array of float, optional 578 The number of bins or the bin edges. (Any valid value to the ``bins`` 579 argument of ``np.histogram`` will work here.) For an exact calculation, 580 each unique value should have its own bin. The default value for bins 581 ensures exact handling of uint8 images: ``bins=256`` results in 582 aliasing problems due to bin width not being equal to 1. 583 584 Returns 585 ------- 586 nu : float 587 The cross-entropy target value as defined in [1]_. 588 589 Notes 590 ----- 591 See Li and Lee, 1993 [1]_; this is the objective function ``threshold_li`` 592 minimizes. This function can be improved but this implementation most 593 closely matches equation 8 in [1]_ and equations 1-3 in [2]_. 594 595 References 596 ---------- 597 .. [1] Li C.H. and Lee C.K. (1993) "Minimum Cross Entropy Thresholding" 598 Pattern Recognition, 26(4): 617-625 599 :DOI:`10.1016/0031-3203(93)90115-D` 600 .. [2] Li C.H. and Tam P.K.S. (1998) "An Iterative Algorithm for Minimum 601 Cross Entropy Thresholding" Pattern Recognition Letters, 18(8): 771-776 602 :DOI:`10.1016/S0167-8655(98)00057-9` 603 """ 604 histogram, bin_edges = np.histogram(image, bins=bins, density=True) 605 bin_centers = np.convolve(bin_edges, [0.5, 0.5], mode='valid') 606 t = np.flatnonzero(bin_centers > threshold)[0] 607 m0a = np.sum(histogram[:t]) # 0th moment, background 608 m0b = np.sum(histogram[t:]) 609 m1a = np.sum(histogram[:t] * bin_centers[:t]) # 1st moment, background 610 m1b = np.sum(histogram[t:] * bin_centers[t:]) 611 mua = m1a / m0a # mean value, background 612 mub = m1b / m0b 613 nu = -m1a * np.log(mua) - m1b * np.log(mub) 614 return nu 615 616 617def threshold_li(image, *, tolerance=None, initial_guess=None, 618 iter_callback=None): 619 """Compute threshold value by Li's iterative Minimum Cross Entropy method. 620 621 Parameters 622 ---------- 623 image : (N, M[, ..., P]) ndarray 624 Grayscale input image. 625 626 tolerance : float, optional 627 Finish the computation when the change in the threshold in an iteration 628 is less than this value. By default, this is half the smallest 629 difference between intensity values in ``image``. 630 631 initial_guess : float or Callable[[array[float]], float], optional 632 Li's iterative method uses gradient descent to find the optimal 633 threshold. If the image intensity histogram contains more than two 634 modes (peaks), the gradient descent could get stuck in a local optimum. 635 An initial guess for the iteration can help the algorithm find the 636 globally-optimal threshold. A float value defines a specific start 637 point, while a callable should take in an array of image intensities 638 and return a float value. Example valid callables include 639 ``numpy.mean`` (default), ``lambda arr: numpy.quantile(arr, 0.95)``, 640 or even :func:`skimage.filters.threshold_otsu`. 641 642 iter_callback : Callable[[float], Any], optional 643 A function that will be called on the threshold at every iteration of 644 the algorithm. 645 646 Returns 647 ------- 648 threshold : float 649 Upper threshold value. All pixels with an intensity higher than 650 this value are assumed to be foreground. 651 652 References 653 ---------- 654 .. [1] Li C.H. and Lee C.K. (1993) "Minimum Cross Entropy Thresholding" 655 Pattern Recognition, 26(4): 617-625 656 :DOI:`10.1016/0031-3203(93)90115-D` 657 .. [2] Li C.H. and Tam P.K.S. (1998) "An Iterative Algorithm for Minimum 658 Cross Entropy Thresholding" Pattern Recognition Letters, 18(8): 771-776 659 :DOI:`10.1016/S0167-8655(98)00057-9` 660 .. [3] Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding 661 Techniques and Quantitative Performance Evaluation" Journal of 662 Electronic Imaging, 13(1): 146-165 663 :DOI:`10.1117/1.1631315` 664 .. [4] ImageJ AutoThresholder code, http://fiji.sc/wiki/index.php/Auto_Threshold 665 666 Examples 667 -------- 668 >>> from skimage.data import camera 669 >>> image = camera() 670 >>> thresh = threshold_li(image) 671 >>> binary = image > thresh 672 """ 673 # Remove nan: 674 image = image[~np.isnan(image)] 675 if image.size == 0: 676 return np.nan 677 678 # Make sure image has more than one value; otherwise, return that value 679 # This works even for np.inf 680 if np.all(image == image.flat[0]): 681 return image.flat[0] 682 683 # At this point, the image only contains np.inf, -np.inf, or valid numbers 684 image = image[np.isfinite(image)] 685 # if there are no finite values in the image, return 0. This is because 686 # at this point we *know* that there are *both* inf and -inf values, 687 # because inf == inf evaluates to True. We might as well separate them. 688 if image.size == 0: 689 return 0. 690 691 # Li's algorithm requires positive image (because of log(mean)) 692 image_min = np.min(image) 693 image -= image_min 694 if image.dtype.kind in 'iu': 695 tolerance = tolerance or 0.5 696 else: 697 tolerance = tolerance or np.min(np.diff(np.unique(image))) / 2 698 699 # Initial estimate for iteration. See "initial_guess" in the parameter list 700 if initial_guess is None: 701 t_next = np.mean(image) 702 elif callable(initial_guess): 703 t_next = initial_guess(image) 704 elif np.isscalar(initial_guess): # convert to new, positive image range 705 t_next = initial_guess - image_min 706 image_max = np.max(image) + image_min 707 if not 0 < t_next < np.max(image): 708 msg = (f'The initial guess for threshold_li must be within the ' 709 f'range of the image. Got {initial_guess} for image min ' 710 f'{image_min} and max {image_max}.') 711 raise ValueError(msg) 712 else: 713 raise TypeError('Incorrect type for `initial_guess`; should be ' 714 'a floating point value, or a function mapping an ' 715 'array to a floating point value.') 716 717 # initial value for t_curr must be different from t_next by at 718 # least the tolerance. Since the image is positive, we ensure this 719 # by setting to a large-enough negative number 720 t_curr = -2 * tolerance 721 722 # Callback on initial iterations 723 if iter_callback is not None: 724 iter_callback(t_next + image_min) 725 726 # Stop the iterations when the difference between the 727 # new and old threshold values is less than the tolerance 728 729 if image.dtype.kind in 'iu': 730 hist, bin_centers = histogram(image.reshape(-1), 731 source_range='image') 732 hist = hist.astype(float) 733 while abs(t_next - t_curr) > tolerance: 734 t_curr = t_next 735 foreground = bin_centers > t_curr 736 background = ~foreground 737 738 mean_fore = np.average(bin_centers[foreground], 739 weights=hist[foreground]) 740 mean_back = np.average(bin_centers[background], 741 weights=hist[background]) 742 743 t_next = ((mean_back - mean_fore) 744 / (np.log(mean_back) - np.log(mean_fore))) 745 746 if iter_callback is not None: 747 iter_callback(t_next + image_min) 748 749 else: 750 while abs(t_next - t_curr) > tolerance: 751 t_curr = t_next 752 foreground = (image > t_curr) 753 mean_fore = np.mean(image[foreground]) 754 mean_back = np.mean(image[~foreground]) 755 756 t_next = ((mean_back - mean_fore) 757 / (np.log(mean_back) - np.log(mean_fore))) 758 759 if iter_callback is not None: 760 iter_callback(t_next + image_min) 761 762 threshold = t_next + image_min 763 return threshold 764 765 766@deprecate_kwarg({'max_iter': 'max_num_iter'}, removed_version="1.0", 767 deprecated_version="0.19") 768def threshold_minimum(image=None, nbins=256, max_num_iter=10000, *, hist=None): 769 """Return threshold value based on minimum method. 770 771 The histogram of the input ``image`` is computed if not provided and 772 smoothed until there are only two maxima. Then the minimum in between is 773 the threshold value. 774 775 Either image or hist must be provided. In case hist is given, the actual 776 histogram of the image is ignored. 777 778 Parameters 779 ---------- 780 image : (N, M[, ..., P]) ndarray, optional 781 Grayscale input image. 782 nbins : int, optional 783 Number of bins used to calculate histogram. This value is ignored for 784 integer arrays. 785 max_num_iter : int, optional 786 Maximum number of iterations to smooth the histogram. 787 hist : array, or 2-tuple of arrays, optional 788 Histogram to determine the threshold from and a corresponding array 789 of bin center intensities. Alternatively, only the histogram can be 790 passed. 791 792 Returns 793 ------- 794 threshold : float 795 Upper threshold value. All pixels with an intensity higher than 796 this value are assumed to be foreground. 797 798 Raises 799 ------ 800 RuntimeError 801 If unable to find two local maxima in the histogram or if the 802 smoothing takes more than 1e4 iterations. 803 804 References 805 ---------- 806 .. [1] C. A. Glasbey, "An analysis of histogram-based thresholding 807 algorithms," CVGIP: Graphical Models and Image Processing, 808 vol. 55, pp. 532-537, 1993. 809 .. [2] Prewitt, JMS & Mendelsohn, ML (1966), "The analysis of cell 810 images", Annals of the New York Academy of Sciences 128: 1035-1053 811 :DOI:`10.1111/j.1749-6632.1965.tb11715.x` 812 813 Examples 814 -------- 815 >>> from skimage.data import camera 816 >>> image = camera() 817 >>> thresh = threshold_minimum(image) 818 >>> binary = image > thresh 819 """ 820 821 def find_local_maxima_idx(hist): 822 # We can't use scipy.signal.argrelmax 823 # as it fails on plateaus 824 maximum_idxs = list() 825 direction = 1 826 827 for i in range(hist.shape[0] - 1): 828 if direction > 0: 829 if hist[i + 1] < hist[i]: 830 direction = -1 831 maximum_idxs.append(i) 832 else: 833 if hist[i + 1] > hist[i]: 834 direction = 1 835 836 return maximum_idxs 837 838 counts, bin_centers = _validate_image_histogram(image, hist, nbins) 839 840 smooth_hist = counts.astype(np.float64, copy=False) 841 842 for counter in range(max_num_iter): 843 smooth_hist = ndi.uniform_filter1d(smooth_hist, 3) 844 maximum_idxs = find_local_maxima_idx(smooth_hist) 845 if len(maximum_idxs) < 3: 846 break 847 848 if len(maximum_idxs) != 2: 849 raise RuntimeError('Unable to find two maxima in histogram') 850 elif counter == max_num_iter - 1: 851 raise RuntimeError('Maximum iteration reached for histogram' 852 'smoothing') 853 854 # Find lowest point between the maxima 855 threshold_idx = np.argmin(smooth_hist[maximum_idxs[0]:maximum_idxs[1] + 1]) 856 857 return bin_centers[maximum_idxs[0] + threshold_idx] 858 859 860def threshold_mean(image): 861 """Return threshold value based on the mean of grayscale values. 862 863 Parameters 864 ---------- 865 image : (N, M[, ..., P]) ndarray 866 Grayscale input image. 867 868 Returns 869 ------- 870 threshold : float 871 Upper threshold value. All pixels with an intensity higher than 872 this value are assumed to be foreground. 873 874 References 875 ---------- 876 .. [1] C. A. Glasbey, "An analysis of histogram-based thresholding 877 algorithms," CVGIP: Graphical Models and Image Processing, 878 vol. 55, pp. 532-537, 1993. 879 :DOI:`10.1006/cgip.1993.1040` 880 881 Examples 882 -------- 883 >>> from skimage.data import camera 884 >>> image = camera() 885 >>> thresh = threshold_mean(image) 886 >>> binary = image > thresh 887 """ 888 return np.mean(image) 889 890 891def threshold_triangle(image, nbins=256): 892 """Return threshold value based on the triangle algorithm. 893 894 Parameters 895 ---------- 896 image : (N, M[, ..., P]) ndarray 897 Grayscale input image. 898 nbins : int, optional 899 Number of bins used to calculate histogram. This value is ignored for 900 integer arrays. 901 902 Returns 903 ------- 904 threshold : float 905 Upper threshold value. All pixels with an intensity higher than 906 this value are assumed to be foreground. 907 908 References 909 ---------- 910 .. [1] Zack, G. W., Rogers, W. E. and Latt, S. A., 1977, 911 Automatic Measurement of Sister Chromatid Exchange Frequency, 912 Journal of Histochemistry and Cytochemistry 25 (7), pp. 741-753 913 :DOI:`10.1177/25.7.70454` 914 .. [2] ImageJ AutoThresholder code, 915 http://fiji.sc/wiki/index.php/Auto_Threshold 916 917 Examples 918 -------- 919 >>> from skimage.data import camera 920 >>> image = camera() 921 >>> thresh = threshold_triangle(image) 922 >>> binary = image > thresh 923 """ 924 # nbins is ignored for integer arrays 925 # so, we recalculate the effective nbins. 926 hist, bin_centers = histogram(image.ravel(), nbins, source_range='image') 927 nbins = len(hist) 928 929 # Find peak, lowest and highest gray levels. 930 arg_peak_height = np.argmax(hist) 931 peak_height = hist[arg_peak_height] 932 arg_low_level, arg_high_level = np.where(hist > 0)[0][[0, -1]] 933 934 # Flip is True if left tail is shorter. 935 flip = arg_peak_height - arg_low_level < arg_high_level - arg_peak_height 936 if flip: 937 hist = hist[::-1] 938 arg_low_level = nbins - arg_high_level - 1 939 arg_peak_height = nbins - arg_peak_height - 1 940 941 # If flip == True, arg_high_level becomes incorrect 942 # but we don't need it anymore. 943 del(arg_high_level) 944 945 # Set up the coordinate system. 946 width = arg_peak_height - arg_low_level 947 x1 = np.arange(width) 948 y1 = hist[x1 + arg_low_level] 949 950 # Normalize. 951 norm = np.sqrt(peak_height**2 + width**2) 952 peak_height /= norm 953 width /= norm 954 955 # Maximize the length. 956 # The ImageJ implementation includes an additional constant when calculating 957 # the length, but here we omit it as it does not affect the location of the 958 # minimum. 959 length = peak_height * x1 - width * y1 960 arg_level = np.argmax(length) + arg_low_level 961 962 if flip: 963 arg_level = nbins - arg_level - 1 964 965 return bin_centers[arg_level] 966 967 968def _mean_std(image, w): 969 """Return local mean and standard deviation of each pixel using a 970 neighborhood defined by a rectangular window size ``w``. 971 The algorithm uses integral images to speedup computation. This is 972 used by :func:`threshold_niblack` and :func:`threshold_sauvola`. 973 974 Parameters 975 ---------- 976 image : (N, M[, ..., P]) ndarray 977 Grayscale input image. 978 w : int, or iterable of int 979 Window size specified as a single odd integer (3, 5, 7, …), 980 or an iterable of length ``image.ndim`` containing only odd 981 integers (e.g. ``(1, 5, 5)``). 982 983 Returns 984 ------- 985 m : ndarray of float, same shape as ``image`` 986 Local mean of the image. 987 s : ndarray of float, same shape as ``image`` 988 Local standard deviation of the image. 989 990 References 991 ---------- 992 .. [1] F. Shafait, D. Keysers, and T. M. Breuel, "Efficient 993 implementation of local adaptive thresholding techniques 994 using integral images." in Document Recognition and 995 Retrieval XV, (San Jose, USA), Jan. 2008. 996 :DOI:`10.1117/12.767755` 997 """ 998 999 if not isinstance(w, Iterable): 1000 w = (w,) * image.ndim 1001 _validate_window_size(w) 1002 1003 float_dtype = _supported_float_type(image.dtype) 1004 pad_width = tuple((k // 2 + 1, k // 2) for k in w) 1005 padded = np.pad(image.astype(float_dtype, copy=False), pad_width, 1006 mode='reflect') 1007 1008 # Note: keep float64 integral images for accuracy. Outputs of 1009 # _correlate_sparse can later be safely cast to float_dtype 1010 integral = integral_image(padded, dtype=np.float64) 1011 padded *= padded 1012 integral_sq = integral_image(padded, dtype=np.float64) 1013 1014 1015 # Create lists of non-zero kernel indices and values 1016 kernel_indices = list(itertools.product(*tuple([(0, _w) for _w in w]))) 1017 kernel_values = [(-1) ** (image.ndim % 2 != np.sum(indices) % 2) 1018 for indices in kernel_indices] 1019 1020 total_window_size = np.prod(w) 1021 kernel_shape = tuple(_w + 1 for _w in w) 1022 m = _correlate_sparse(integral, kernel_shape, kernel_indices, 1023 kernel_values) 1024 m = m.astype(float_dtype, copy=False) 1025 m /= total_window_size 1026 g2 = _correlate_sparse(integral_sq, kernel_shape, kernel_indices, 1027 kernel_values) 1028 g2 = g2.astype(float_dtype, copy=False) 1029 g2 /= total_window_size 1030 # Note: we use np.clip because g2 is not guaranteed to be greater than 1031 # m*m when floating point error is considered 1032 s = np.sqrt(np.clip(g2 - m * m, 0, None)) 1033 return m, s 1034 1035 1036def threshold_niblack(image, window_size=15, k=0.2): 1037 """Applies Niblack local threshold to an array. 1038 1039 A threshold T is calculated for every pixel in the image using the 1040 following formula:: 1041 1042 T = m(x,y) - k * s(x,y) 1043 1044 where m(x,y) and s(x,y) are the mean and standard deviation of 1045 pixel (x,y) neighborhood defined by a rectangular window with size w 1046 times w centered around the pixel. k is a configurable parameter 1047 that weights the effect of standard deviation. 1048 1049 Parameters 1050 ---------- 1051 image : (N, M[, ..., P]) ndarray 1052 Grayscale input image. 1053 window_size : int, or iterable of int, optional 1054 Window size specified as a single odd integer (3, 5, 7, …), 1055 or an iterable of length ``image.ndim`` containing only odd 1056 integers (e.g. ``(1, 5, 5)``). 1057 k : float, optional 1058 Value of parameter k in threshold formula. 1059 1060 Returns 1061 ------- 1062 threshold : (N, M) ndarray 1063 Threshold mask. All pixels with an intensity higher than 1064 this value are assumed to be foreground. 1065 1066 Notes 1067 ----- 1068 This algorithm is originally designed for text recognition. 1069 1070 The Bradley threshold is a particular case of the Niblack 1071 one, being equivalent to 1072 1073 >>> from skimage import data 1074 >>> image = data.page() 1075 >>> q = 1 1076 >>> threshold_image = threshold_niblack(image, k=0) * q 1077 1078 for some value ``q``. By default, Bradley and Roth use ``q=1``. 1079 1080 1081 References 1082 ---------- 1083 .. [1] W. Niblack, An introduction to Digital Image Processing, 1084 Prentice-Hall, 1986. 1085 .. [2] D. Bradley and G. Roth, "Adaptive thresholding using Integral 1086 Image", Journal of Graphics Tools 12(2), pp. 13-21, 2007. 1087 :DOI:`10.1080/2151237X.2007.10129236` 1088 1089 Examples 1090 -------- 1091 >>> from skimage import data 1092 >>> image = data.page() 1093 >>> threshold_image = threshold_niblack(image, window_size=7, k=0.1) 1094 """ 1095 m, s = _mean_std(image, window_size) 1096 return m - k * s 1097 1098 1099def threshold_sauvola(image, window_size=15, k=0.2, r=None): 1100 """Applies Sauvola local threshold to an array. Sauvola is a 1101 modification of Niblack technique. 1102 1103 In the original method a threshold T is calculated for every pixel 1104 in the image using the following formula:: 1105 1106 T = m(x,y) * (1 + k * ((s(x,y) / R) - 1)) 1107 1108 where m(x,y) and s(x,y) are the mean and standard deviation of 1109 pixel (x,y) neighborhood defined by a rectangular window with size w 1110 times w centered around the pixel. k is a configurable parameter 1111 that weights the effect of standard deviation. 1112 R is the maximum standard deviation of a grayscale image. 1113 1114 Parameters 1115 ---------- 1116 image : (N, M[, ..., P]) ndarray 1117 Grayscale input image. 1118 window_size : int, or iterable of int, optional 1119 Window size specified as a single odd integer (3, 5, 7, …), 1120 or an iterable of length ``image.ndim`` containing only odd 1121 integers (e.g. ``(1, 5, 5)``). 1122 k : float, optional 1123 Value of the positive parameter k. 1124 r : float, optional 1125 Value of R, the dynamic range of standard deviation. 1126 If None, set to the half of the image dtype range. 1127 1128 Returns 1129 ------- 1130 threshold : (N, M) ndarray 1131 Threshold mask. All pixels with an intensity higher than 1132 this value are assumed to be foreground. 1133 1134 Notes 1135 ----- 1136 This algorithm is originally designed for text recognition. 1137 1138 References 1139 ---------- 1140 .. [1] J. Sauvola and M. Pietikainen, "Adaptive document image 1141 binarization," Pattern Recognition 33(2), 1142 pp. 225-236, 2000. 1143 :DOI:`10.1016/S0031-3203(99)00055-2` 1144 1145 Examples 1146 -------- 1147 >>> from skimage import data 1148 >>> image = data.page() 1149 >>> t_sauvola = threshold_sauvola(image, window_size=15, k=0.2) 1150 >>> binary_image = image > t_sauvola 1151 """ 1152 if r is None: 1153 imin, imax = dtype_limits(image, clip_negative=False) 1154 r = 0.5 * (imax - imin) 1155 m, s = _mean_std(image, window_size) 1156 return m * (1 + k * ((s / r) - 1)) 1157 1158 1159def apply_hysteresis_threshold(image, low, high): 1160 """Apply hysteresis thresholding to ``image``. 1161 1162 This algorithm finds regions where ``image`` is greater than ``high`` 1163 OR ``image`` is greater than ``low`` *and* that region is connected to 1164 a region greater than ``high``. 1165 1166 Parameters 1167 ---------- 1168 image : array, shape (M,[ N, ..., P]) 1169 Grayscale input image. 1170 low : float, or array of same shape as ``image`` 1171 Lower threshold. 1172 high : float, or array of same shape as ``image`` 1173 Higher threshold. 1174 1175 Returns 1176 ------- 1177 thresholded : array of bool, same shape as ``image`` 1178 Array in which ``True`` indicates the locations where ``image`` 1179 was above the hysteresis threshold. 1180 1181 Examples 1182 -------- 1183 >>> image = np.array([1, 2, 3, 2, 1, 2, 1, 3, 2]) 1184 >>> apply_hysteresis_threshold(image, 1.5, 2.5).astype(int) 1185 array([0, 1, 1, 1, 0, 0, 0, 1, 1]) 1186 1187 References 1188 ---------- 1189 .. [1] J. Canny. A computational approach to edge detection. 1190 IEEE Transactions on Pattern Analysis and Machine Intelligence. 1191 1986; vol. 8, pp.679-698. 1192 :DOI:`10.1109/TPAMI.1986.4767851` 1193 """ 1194 low = np.clip(low, a_min=None, a_max=high) # ensure low always below high 1195 mask_low = image > low 1196 mask_high = image > high 1197 # Connected components of mask_low 1198 labels_low, num_labels = ndi.label(mask_low) 1199 # Check which connected components contain pixels from mask_high 1200 sums = ndi.sum(mask_high, labels_low, np.arange(num_labels + 1)) 1201 connected_to_high = sums > 0 1202 thresholded = connected_to_high[labels_low] 1203 return thresholded 1204 1205 1206def threshold_multiotsu(image=None, classes=3, nbins=256, *, hist=None): 1207 r"""Generate `classes`-1 threshold values to divide gray levels in `image`, 1208 following Otsu's method for multiple classes. 1209 1210 The threshold values are chosen to maximize the total sum of pairwise 1211 variances between the thresholded graylevel classes. See Notes and [1]_ 1212 for more details. 1213 1214 Either image or hist must be provided. If hist is provided, the actual 1215 histogram of the image is ignored. 1216 1217 Parameters 1218 ---------- 1219 image : (N, M[, ..., P]) ndarray, optional 1220 Grayscale input image. 1221 classes : int, optional 1222 Number of classes to be thresholded, i.e. the number of resulting 1223 regions. 1224 nbins : int, optional 1225 Number of bins used to calculate the histogram. This value is ignored 1226 for integer arrays. 1227 hist : array, or 2-tuple of arrays, optional 1228 Histogram from which to determine the threshold, and optionally a 1229 corresponding array of bin center intensities. If no hist provided, 1230 this function will compute it from the image (see notes). 1231 1232 Returns 1233 ------- 1234 thresh : array 1235 Array containing the threshold values for the desired classes. 1236 1237 Raises 1238 ------ 1239 ValueError 1240 If ``image`` contains less grayscale value then the desired 1241 number of classes. 1242 1243 Notes 1244 ----- 1245 This implementation relies on a Cython function whose complexity 1246 is :math:`O\left(\frac{Ch^{C-1}}{(C-1)!}\right)`, where :math:`h` 1247 is the number of histogram bins and :math:`C` is the number of 1248 classes desired. 1249 1250 If no hist is given, this function will make use of 1251 `skimage.exposure.histogram`, which behaves differently than 1252 `np.histogram`. While both allowed, use the former for consistent 1253 behaviour. 1254 1255 The input image must be grayscale. 1256 1257 References 1258 ---------- 1259 .. [1] Liao, P-S., Chen, T-S. and Chung, P-C., "A fast algorithm for 1260 multilevel thresholding", Journal of Information Science and 1261 Engineering 17 (5): 713-727, 2001. Available at: 1262 <https://ftp.iis.sinica.edu.tw/JISE/2001/200109_01.pdf> 1263 :DOI:`10.6688/JISE.2001.17.5.1` 1264 .. [2] Tosa, Y., "Multi-Otsu Threshold", a java plugin for ImageJ. 1265 Available at: 1266 <http://imagej.net/plugins/download/Multi_OtsuThreshold.java> 1267 1268 Examples 1269 -------- 1270 >>> from skimage.color import label2rgb 1271 >>> from skimage import data 1272 >>> image = data.camera() 1273 >>> thresholds = threshold_multiotsu(image) 1274 >>> regions = np.digitize(image, bins=thresholds) 1275 >>> regions_colorized = label2rgb(regions) 1276 """ 1277 if image is not None and image.ndim > 2 and image.shape[-1] in (3, 4): 1278 warn(f'threshold_multiotsu is expected to work correctly only for ' 1279 f'grayscale images; image shape {image.shape} looks like ' 1280 f'that of an RGB image.') 1281 1282 # calculating the histogram and the probability of each gray level. 1283 prob, bin_centers = _validate_image_histogram(image, hist, nbins, 1284 normalize=True) 1285 prob = prob.astype('float32') 1286 1287 nvalues = np.count_nonzero(prob) 1288 if nvalues < classes: 1289 msg = (f'The input image has only {nvalues} different values. ' 1290 f'It cannot be thresholded in {classes} classes.') 1291 raise ValueError(msg) 1292 elif nvalues == classes: 1293 thresh_idx = np.where(prob > 0)[0][:-1] 1294 else: 1295 # Get threshold indices 1296 try: 1297 thresh_idx = _get_multiotsu_thresh_indices_lut(prob, classes - 1) 1298 except MemoryError: 1299 # Don't use LUT if the number of bins is too large (if the 1300 # image is uint16 for example): in this case, the 1301 # allocated memory is too large. 1302 thresh_idx = _get_multiotsu_thresh_indices(prob, classes - 1) 1303 1304 thresh = bin_centers[thresh_idx] 1305 1306 return thresh 1307