scikit-image

1
import inspect
2
import itertools
3
import math
4
from collections import OrderedDict
5
from collections.abc import Iterable
6

7
import numpy as np
8
from scipy import ndimage as ndi
9

10
from .._shared.filters import gaussian
11
from .._shared.utils import _supported_float_type, warn
12
from .._shared.version_requirements import require
13
from ..exposure import histogram
14
from ..filters._multiotsu import (
15
    _get_multiotsu_thresh_indices,
16
    _get_multiotsu_thresh_indices_lut,
17
)
18
from ..transform import integral_image
19
from ..util import dtype_limits
20
from ._sparse import _correlate_sparse, _validate_window_size
21

22
__all__ = [
23
    'try_all_threshold',
24
    'threshold_otsu',
25
    'threshold_yen',
26
    'threshold_isodata',
27
    'threshold_li',
28
    'threshold_local',
29
    'threshold_minimum',
30
    'threshold_mean',
31
    'threshold_niblack',
32
    'threshold_sauvola',
33
    'threshold_triangle',
34
    'apply_hysteresis_threshold',
35
    'threshold_multiotsu',
36
]
37

38

39
def _try_all(image, methods=None, figsize=None, num_cols=2, verbose=True):
40
    """Returns a figure comparing the outputs of different methods.
41

42
    Parameters
43
    ----------
44
    image : (M, N) ndarray
45
        Input image.
46
    methods : dict, optional
47
        Names and associated functions.
48
        Functions must take and return an image.
49
    figsize : tuple, optional
50
        Figure size (in inches).
51
    num_cols : int, optional
52
        Number of columns.
53
    verbose : bool, optional
54
        Print function name for each method.
55

56
    Returns
57
    -------
58
    fig, ax : tuple
59
        Matplotlib figure and axes.
60
    """
61
    from matplotlib import pyplot as plt
62

63
    # Compute the image histogram for better performances
64
    nbins = 256  # Default in threshold functions
65
    hist = histogram(image.reshape(-1), nbins, source_range='image')
66

67
    # Handle default value
68
    methods = methods or {}
69

70
    num_rows = math.ceil((len(methods) + 1.0) / num_cols)
71
    fig, ax = plt.subplots(
72
        num_rows, num_cols, figsize=figsize, sharex=True, sharey=True
73
    )
74
    ax = ax.reshape(-1)
75

76
    ax[0].imshow(image, cmap=plt.cm.gray)
77
    ax[0].set_title('Original')
78

79
    i = 1
80
    for name, func in methods.items():
81
        # Use precomputed histogram for supporting functions
82
        sig = inspect.signature(func)
83
        _kwargs = dict(hist=hist) if 'hist' in sig.parameters else {}
84

85
        ax[i].set_title(name)
86
        try:
87
            ax[i].imshow(func(image, **_kwargs), cmap=plt.cm.gray)
88
        except Exception as e:
89
            ax[i].text(
90
                0.5,
91
                0.5,
92
                f"{type(e).__name__}",
93
                ha="center",
94
                va="center",
95
                transform=ax[i].transAxes,
96
            )
97
        i += 1
98
        if verbose:
99
            print(func.__orifunc__)
100

101
    for a in ax:
102
        a.axis('off')
103

104
    fig.tight_layout()
105
    return fig, ax
106

107

108
@require("matplotlib", ">=3.3")
109
def try_all_threshold(image, figsize=(8, 5), verbose=True):
110
    """Returns a figure comparing the outputs of different thresholding methods.
111

112
    Parameters
113
    ----------
114
    image : (M, N) ndarray
115
        Input image.
116
    figsize : tuple, optional
117
        Figure size (in inches).
118
    verbose : bool, optional
119
        Print function name for each method.
120

121
    Returns
122
    -------
123
    fig, ax : tuple
124
        Matplotlib figure and axes.
125

126
    Notes
127
    -----
128
    The following algorithms are used:
129

130
    * isodata
131
    * li
132
    * mean
133
    * minimum
134
    * otsu
135
    * triangle
136
    * yen
137

138
    Examples
139
    --------
140
    .. testsetup::
141
        >>> import pytest; _ = pytest.importorskip('matplotlib')
142

143
    >>> from skimage.data import text
144
    >>> fig, ax = try_all_threshold(text(), figsize=(10, 6), verbose=False)
145
    """
146

147
    def thresh(func):
148
        """
149
        A wrapper function to return a thresholded image.
150
        """
151

152
        def wrapper(im):
153
            return im > func(im)
154

155
        try:
156
            wrapper.__orifunc__ = func.__orifunc__
157
        except AttributeError:
158
            wrapper.__orifunc__ = func.__module__ + '.' + func.__name__
159
        return wrapper
160

161
    # Global algorithms.
162
    methods = OrderedDict(
163
        {
164
            'Isodata': thresh(threshold_isodata),
165
            'Li': thresh(threshold_li),
166
            'Mean': thresh(threshold_mean),
167
            'Minimum': thresh(threshold_minimum),
168
            'Otsu': thresh(threshold_otsu),
169
            'Triangle': thresh(threshold_triangle),
170
            'Yen': thresh(threshold_yen),
171
        }
172
    )
173

174
    return _try_all(image, figsize=figsize, methods=methods, verbose=verbose)
175

176

177
def threshold_local(
178
    image, block_size=3, method='gaussian', offset=0, mode='reflect', param=None, cval=0
179
):
180
    """Compute a threshold mask image based on local pixel neighborhood.
181

182
    Also known as adaptive or dynamic thresholding. The threshold value is
183
    the weighted mean for the local neighborhood of a pixel subtracted by a
184
    constant. Alternatively the threshold can be determined dynamically by a
185
    given function, using the 'generic' method.
186

187
    Parameters
188
    ----------
189
    image : (M, N[, ...]) ndarray
190
        Grayscale input image.
191
    block_size : int or sequence of int
192
        Odd size of pixel neighborhood which is used to calculate the
193
        threshold value (e.g. 3, 5, 7, ..., 21, ...).
194
    method : {'generic', 'gaussian', 'mean', 'median'}, optional
195
        Method used to determine adaptive threshold for local neighborhood in
196
        weighted mean image.
197

198
        * 'generic': use custom function (see ``param`` parameter)
199
        * 'gaussian': apply gaussian filter (see ``param`` parameter for custom\
200
                      sigma value)
201
        * 'mean': apply arithmetic mean filter
202
        * 'median': apply median rank filter
203

204
        By default, the 'gaussian' method is used.
205
    offset : float, optional
206
        Constant subtracted from weighted mean of neighborhood to calculate
207
        the local threshold value. Default offset is 0.
208
    mode : {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}, optional
209
        The mode parameter determines how the array borders are handled, where
210
        cval is the value when mode is equal to 'constant'.
211
        Default is 'reflect'.
212
    param : {int, function}, optional
213
        Either specify sigma for 'gaussian' method or function object for
214
        'generic' method. This functions takes the flat array of local
215
        neighborhood as a single argument and returns the calculated
216
        threshold for the centre pixel.
217
    cval : float, optional
218
        Value to fill past edges of input if mode is 'constant'.
219

220
    Returns
221
    -------
222
    threshold : (M, N[, ...]) ndarray
223
        Threshold image. All pixels in the input image higher than the
224
        corresponding pixel in the threshold image are considered foreground.
225

226
    References
227
    ----------
228
    .. [1] Gonzalez, R. C. and Wood, R. E. "Digital Image Processing
229
           (2nd Edition)." Prentice-Hall Inc., 2002: 600--612.
230
           ISBN: 0-201-18075-8
231

232
    Examples
233
    --------
234
    >>> from skimage.data import camera
235
    >>> image = camera()[:50, :50]
236
    >>> binary_image1 = image > threshold_local(image, 15, 'mean')
237
    >>> func = lambda arr: arr.mean()
238
    >>> binary_image2 = image > threshold_local(image, 15, 'generic',
239
    ...                                         param=func)
240

241
    """
242

243
    if np.isscalar(block_size):
244
        block_size = (block_size,) * image.ndim
245
    elif len(block_size) != image.ndim:
246
        raise ValueError("len(block_size) must equal image.ndim.")
247
    block_size = tuple(block_size)
248
    if any(b % 2 == 0 for b in block_size):
249
        raise ValueError(
250
            f'block_size must be odd! Given block_size '
251
            f'{block_size} contains even values.'
252
        )
253
    float_dtype = _supported_float_type(image.dtype)
254
    image = image.astype(float_dtype, copy=False)
255
    thresh_image = np.zeros(image.shape, dtype=float_dtype)
256
    if method == 'generic':
257
        ndi.generic_filter(
258
            image, param, block_size, output=thresh_image, mode=mode, cval=cval
259
        )
260
    elif method == 'gaussian':
261
        if param is None:
262
            # automatically determine sigma which covers > 99% of distribution
263
            sigma = tuple([(b - 1) / 6.0 for b in block_size])
264
        else:
265
            sigma = param
266
        gaussian(image, sigma=sigma, out=thresh_image, mode=mode, cval=cval)
267
    elif method == 'mean':
268
        ndi.uniform_filter(image, block_size, output=thresh_image, mode=mode, cval=cval)
269
    elif method == 'median':
270
        ndi.median_filter(image, block_size, output=thresh_image, mode=mode, cval=cval)
271
    else:
272
        raise ValueError(
273
            "Invalid method specified. Please use `generic`, "
274
            "`gaussian`, `mean`, or `median`."
275
        )
276

277
    return thresh_image - offset
278

279

280
def _validate_image_histogram(image, hist, nbins=None, normalize=False):
281
    """Ensure that either image or hist were given, return valid histogram.
282

283
    If hist is given, image is ignored.
284

285
    Parameters
286
    ----------
287
    image : array or None
288
        Grayscale image.
289
    hist : array, 2-tuple of array, or None
290
        Histogram, either a 1D counts array, or an array of counts together
291
        with an array of bin centers.
292
    nbins : int, optional
293
        The number of bins with which to compute the histogram, if `hist` is
294
        None.
295
    normalize : bool
296
        If hist is not given, it will be computed by this function. This
297
        parameter determines whether the computed histogram is normalized
298
        (i.e. entries sum up to 1) or not.
299

300
    Returns
301
    -------
302
    counts : 1D array of float
303
        Each element is the number of pixels falling in each intensity bin.
304
    bin_centers : 1D array
305
        Each element is the value corresponding to the center of each intensity
306
        bin.
307

308
    Raises
309
    ------
310
    ValueError : if image and hist are both None
311
    """
312
    if image is None and hist is None:
313
        raise Exception("Either image or hist must be provided.")
314

315
    if hist is not None:
316
        if isinstance(hist, (tuple, list)):
317
            counts, bin_centers = hist
318
        else:
319
            counts = hist
320
            bin_centers = np.arange(counts.size)
321

322
        if counts[0] == 0 or counts[-1] == 0:
323
            # Trim histogram from both ends by removing starting and
324
            # ending zeroes as in histogram(..., source_range="image")
325
            cond = counts > 0
326
            start = np.argmax(cond)
327
            end = cond.size - np.argmax(cond[::-1])
328
            counts, bin_centers = counts[start:end], bin_centers[start:end]
329
    else:
330
        counts, bin_centers = histogram(
331
            image.reshape(-1), nbins, source_range='image', normalize=normalize
332
        )
333
    return counts.astype('float32', copy=False), bin_centers
334

335

336
def threshold_otsu(image=None, nbins=256, *, hist=None):
337
    """Return threshold value based on Otsu's method.
338

339
    Either image or hist must be provided. If hist is provided, the actual
340
    histogram of the image is ignored.
341

342
    Parameters
343
    ----------
344
    image : (M, N[, ...]) ndarray, optional
345
        Grayscale input image.
346
    nbins : int, optional
347
        Number of bins used to calculate histogram. This value is ignored for
348
        integer arrays.
349
    hist : array, or 2-tuple of arrays, optional
350
        Histogram from which to determine the threshold, and optionally a
351
        corresponding array of bin center intensities. If no hist provided,
352
        this function will compute it from the image.
353

354

355
    Returns
356
    -------
357
    threshold : float
358
        Upper threshold value. All pixels with an intensity higher than
359
        this value are assumed to be foreground.
360

361
    References
362
    ----------
363
    .. [1] Wikipedia, https://en.wikipedia.org/wiki/Otsu's_Method
364

365
    Examples
366
    --------
367
    >>> from skimage.data import camera
368
    >>> image = camera()
369
    >>> thresh = threshold_otsu(image)
370
    >>> binary = image <= thresh
371

372
    Notes
373
    -----
374
    The input image must be grayscale.
375
    """
376
    if image is not None and image.ndim > 2 and image.shape[-1] in (3, 4):
377
        warn(
378
            f'threshold_otsu is expected to work correctly only for '
379
            f'grayscale images; image shape {image.shape} looks like '
380
            f'that of an RGB image.'
381
        )
382

383
    # Check if the image has more than one intensity value; if not, return that
384
    # value
385
    if image is not None:
386
        first_pixel = image.reshape(-1)[0]
387
        if np.all(image == first_pixel):
388
            return first_pixel
389

390
    counts, bin_centers = _validate_image_histogram(image, hist, nbins)
391

392
    # class probabilities for all possible thresholds
393
    weight1 = np.cumsum(counts)
394
    weight2 = np.cumsum(counts[::-1])[::-1]
395
    # class means for all possible thresholds
396
    mean1 = np.cumsum(counts * bin_centers) / weight1
397
    mean2 = (np.cumsum((counts * bin_centers)[::-1]) / weight2[::-1])[::-1]
398

399
    # Clip ends to align class 1 and class 2 variables:
400
    # The last value of ``weight1``/``mean1`` should pair with zero values in
401
    # ``weight2``/``mean2``, which do not exist.
402
    variance12 = weight1[:-1] * weight2[1:] * (mean1[:-1] - mean2[1:]) ** 2
403

404
    idx = np.argmax(variance12)
405
    threshold = bin_centers[idx]
406

407
    return threshold
408

409

410
def threshold_yen(image=None, nbins=256, *, hist=None):
411
    """Return threshold value based on Yen's method.
412
    Either image or hist must be provided. In case hist is given, the actual
413
    histogram of the image is ignored.
414

415
    Parameters
416
    ----------
417
    image : (M, N[, ...]) ndarray
418
        Grayscale input image.
419
    nbins : int, optional
420
        Number of bins used to calculate histogram. This value is ignored for
421
        integer arrays.
422
    hist : array, or 2-tuple of arrays, optional
423
        Histogram from which to determine the threshold, and optionally a
424
        corresponding array of bin center intensities.
425
        An alternative use of this function is to pass it only hist.
426

427
    Returns
428
    -------
429
    threshold : float
430
        Upper threshold value. All pixels with an intensity higher than
431
        this value are assumed to be foreground.
432

433
    References
434
    ----------
435
    .. [1] Yen J.C., Chang F.J., and Chang S. (1995) "A New Criterion
436
           for Automatic Multilevel Thresholding" IEEE Trans. on Image
437
           Processing, 4(3): 370-378. :DOI:`10.1109/83.366472`
438
    .. [2] Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding
439
           Techniques and Quantitative Performance Evaluation" Journal of
440
           Electronic Imaging, 13(1): 146-165, :DOI:`10.1117/1.1631315`
441
           http://www.busim.ee.boun.edu.tr/~sankur/SankurFolder/Threshold_survey.pdf
442
    .. [3] ImageJ AutoThresholder code, http://fiji.sc/wiki/index.php/Auto_Threshold
443

444
    Examples
445
    --------
446
    >>> from skimage.data import camera
447
    >>> image = camera()
448
    >>> thresh = threshold_yen(image)
449
    >>> binary = image <= thresh
450
    """
451
    counts, bin_centers = _validate_image_histogram(image, hist, nbins)
452

453
    # On blank images (e.g. filled with 0) with int dtype, `histogram()`
454
    # returns ``bin_centers`` containing only one value. Speed up with it.
455
    if bin_centers.size == 1:
456
        return bin_centers[0]
457

458
    # Calculate probability mass function
459
    pmf = counts.astype('float32', copy=False) / counts.sum()
460
    P1 = np.cumsum(pmf)  # Cumulative normalized histogram
461
    P1_sq = np.cumsum(pmf**2)
462
    # Get cumsum calculated from end of squared array:
463
    P2_sq = np.cumsum(pmf[::-1] ** 2)[::-1]
464
    # P2_sq indexes is shifted +1. I assume, with P1[:-1] it's help avoid
465
    # '-inf' in crit. ImageJ Yen implementation replaces those values by zero.
466
    crit = np.log(((P1_sq[:-1] * P2_sq[1:]) ** -1) * (P1[:-1] * (1.0 - P1[:-1])) ** 2)
467
    return bin_centers[crit.argmax()]
468

469

470
def threshold_isodata(image=None, nbins=256, return_all=False, *, hist=None):
471
    """Return threshold value(s) based on ISODATA method.
472

473
    Histogram-based threshold, known as Ridler-Calvard method or inter-means.
474
    Threshold values returned satisfy the following equality::
475

476
        threshold = (image[image <= threshold].mean() +
477
                     image[image > threshold].mean()) / 2.0
478

479
    That is, returned thresholds are intensities that separate the image into
480
    two groups of pixels, where the threshold intensity is midway between the
481
    mean intensities of these groups.
482

483
    For integer images, the above equality holds to within one; for floating-
484
    point images, the equality holds to within the histogram bin-width.
485

486
    Either image or hist must be provided. In case hist is given, the actual
487
    histogram of the image is ignored.
488

489
    Parameters
490
    ----------
491
    image : (M, N[, ...]) ndarray
492
        Grayscale input image.
493
    nbins : int, optional
494
        Number of bins used to calculate histogram. This value is ignored for
495
        integer arrays.
496
    return_all : bool, optional
497
        If False (default), return only the lowest threshold that satisfies
498
        the above equality. If True, return all valid thresholds.
499
    hist : array, or 2-tuple of arrays, optional
500
        Histogram to determine the threshold from and a corresponding array
501
        of bin center intensities. Alternatively, only the histogram can be
502
        passed.
503

504
    Returns
505
    -------
506
    threshold : float or int or array
507
        Threshold value(s).
508

509
    References
510
    ----------
511
    .. [1] Ridler, TW & Calvard, S (1978), "Picture thresholding using an
512
           iterative selection method"
513
           IEEE Transactions on Systems, Man and Cybernetics 8: 630-632,
514
           :DOI:`10.1109/TSMC.1978.4310039`
515
    .. [2] Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding
516
           Techniques and Quantitative Performance Evaluation" Journal of
517
           Electronic Imaging, 13(1): 146-165,
518
           http://www.busim.ee.boun.edu.tr/~sankur/SankurFolder/Threshold_survey.pdf
519
           :DOI:`10.1117/1.1631315`
520
    .. [3] ImageJ AutoThresholder code,
521
           http://fiji.sc/wiki/index.php/Auto_Threshold
522

523
    Examples
524
    --------
525
    >>> from skimage.data import coins
526
    >>> image = coins()
527
    >>> thresh = threshold_isodata(image)
528
    >>> binary = image > thresh
529
    """
530
    counts, bin_centers = _validate_image_histogram(image, hist, nbins)
531

532
    # image only contains one unique value
533
    if len(bin_centers) == 1:
534
        if return_all:
535
            return bin_centers
536
        else:
537
            return bin_centers[0]
538

539
    counts = counts.astype('float32', copy=False)
540

541
    # csuml and csumh contain the count of pixels in that bin or lower, and
542
    # in all bins strictly higher than that bin, respectively
543
    csuml = np.cumsum(counts)
544
    csumh = csuml[-1] - csuml
545

546
    # intensity_sum contains the total pixel intensity from each bin
547
    intensity_sum = counts * bin_centers
548

549
    # l and h contain average value of all pixels in that bin or lower, and
550
    # in all bins strictly higher than that bin, respectively.
551
    # Note that since exp.histogram does not include empty bins at the low or
552
    # high end of the range, csuml and csumh are strictly > 0, except in the
553
    # last bin of csumh, which is zero by construction.
554
    # So no worries about division by zero in the following lines, except
555
    # for the last bin, but we can ignore that because no valid threshold
556
    # can be in the top bin.
557
    # To avoid the division by zero, we simply skip over the last element in
558
    # all future computation.
559
    csum_intensity = np.cumsum(intensity_sum)
560
    lower = csum_intensity[:-1] / csuml[:-1]
561
    higher = (csum_intensity[-1] - csum_intensity[:-1]) / csumh[:-1]
562

563
    # isodata finds threshold values that meet the criterion t = (l + m)/2
564
    # where l is the mean of all pixels <= t and h is the mean of all pixels
565
    # > t, as calculated above. So we are looking for places where
566
    # (l + m) / 2 equals the intensity value for which those l and m figures
567
    # were calculated -- which is, of course, the histogram bin centers.
568
    # We only require this equality to be within the precision of the bin
569
    # width, of course.
570
    all_mean = (lower + higher) / 2.0
571
    bin_width = bin_centers[1] - bin_centers[0]
572

573
    # Look only at thresholds that are below the actual all_mean value,
574
    # for consistency with the threshold being included in the lower pixel
575
    # group. Otherwise, can get thresholds that are not actually fixed-points
576
    # of the isodata algorithm. For float images, this matters less, since
577
    # there really can't be any guarantees anymore anyway.
578
    distances = all_mean - bin_centers[:-1]
579
    thresholds = bin_centers[:-1][(distances >= 0) & (distances < bin_width)]
580

581
    if return_all:
582
        return thresholds
583
    else:
584
        return thresholds[0]
585

586

587
# Computing a histogram using np.histogram on a uint8 image with bins=256
588
# doesn't work and results in aliasing problems. We use a fully specified set
589
# of bins to ensure that each uint8 value false into its own bin.
590
_DEFAULT_ENTROPY_BINS = tuple(np.arange(-0.5, 255.51, 1))
591

592

593
def _cross_entropy(image, threshold, bins=_DEFAULT_ENTROPY_BINS):
594
    """Compute cross-entropy between distributions above and below a threshold.
595

596
    Parameters
597
    ----------
598
    image : array
599
        The input array of values.
600
    threshold : float
601
        The value dividing the foreground and background in ``image``.
602
    bins : int or array of float, optional
603
        The number of bins or the bin edges. (Any valid value to the ``bins``
604
        argument of ``np.histogram`` will work here.) For an exact calculation,
605
        each unique value should have its own bin. The default value for bins
606
        ensures exact handling of uint8 images: ``bins=256`` results in
607
        aliasing problems due to bin width not being equal to 1.
608

609
    Returns
610
    -------
611
    nu : float
612
        The cross-entropy target value as defined in [1]_.
613

614
    Notes
615
    -----
616
    See Li and Lee, 1993 [1]_; this is the objective function ``threshold_li``
617
    minimizes. This function can be improved but this implementation most
618
    closely matches equation 8 in [1]_ and equations 1-3 in [2]_.
619

620
    References
621
    ----------
622
    .. [1] Li C.H. and Lee C.K. (1993) "Minimum Cross Entropy Thresholding"
623
           Pattern Recognition, 26(4): 617-625
624
           :DOI:`10.1016/0031-3203(93)90115-D`
625
    .. [2] Li C.H. and Tam P.K.S. (1998) "An Iterative Algorithm for Minimum
626
           Cross Entropy Thresholding" Pattern Recognition Letters, 18(8): 771-776
627
           :DOI:`10.1016/S0167-8655(98)00057-9`
628
    """
629
    histogram, bin_edges = np.histogram(image, bins=bins, density=True)
630
    bin_centers = np.convolve(bin_edges, [0.5, 0.5], mode='valid')
631
    t = np.flatnonzero(bin_centers > threshold)[0]
632
    m0a = np.sum(histogram[:t])  # 0th moment, background
633
    m0b = np.sum(histogram[t:])
634
    m1a = np.sum(histogram[:t] * bin_centers[:t])  # 1st moment, background
635
    m1b = np.sum(histogram[t:] * bin_centers[t:])
636
    mua = m1a / m0a  # mean value, background
637
    mub = m1b / m0b
638
    nu = -m1a * np.log(mua) - m1b * np.log(mub)
639
    return nu
640

641

642
def threshold_li(image, *, tolerance=None, initial_guess=None, iter_callback=None):
643
    """Compute threshold value by Li's iterative Minimum Cross Entropy method.
644

645
    Parameters
646
    ----------
647
    image : (M, N[, ...]) ndarray
648
        Grayscale input image.
649
    tolerance : float, optional
650
        Finish the computation when the change in the threshold in an iteration
651
        is less than this value. By default, this is half the smallest
652
        difference between intensity values in ``image``.
653
    initial_guess : float or Callable[[array[float]], float], optional
654
        Li's iterative method uses gradient descent to find the optimal
655
        threshold. If the image intensity histogram contains more than two
656
        modes (peaks), the gradient descent could get stuck in a local optimum.
657
        An initial guess for the iteration can help the algorithm find the
658
        globally-optimal threshold. A float value defines a specific start
659
        point, while a callable should take in an array of image intensities
660
        and return a float value. Example valid callables include
661
        ``numpy.mean`` (default), ``lambda arr: numpy.quantile(arr, 0.95)``,
662
        or even :func:`skimage.filters.threshold_otsu`.
663
    iter_callback : Callable[[float], Any], optional
664
        A function that will be called on the threshold at every iteration of
665
        the algorithm.
666

667
    Returns
668
    -------
669
    threshold : float
670
        Upper threshold value. All pixels with an intensity higher than
671
        this value are assumed to be foreground.
672

673
    References
674
    ----------
675
    .. [1] Li C.H. and Lee C.K. (1993) "Minimum Cross Entropy Thresholding"
676
           Pattern Recognition, 26(4): 617-625
677
           :DOI:`10.1016/0031-3203(93)90115-D`
678
    .. [2] Li C.H. and Tam P.K.S. (1998) "An Iterative Algorithm for Minimum
679
           Cross Entropy Thresholding" Pattern Recognition Letters, 18(8): 771-776
680
           :DOI:`10.1016/S0167-8655(98)00057-9`
681
    .. [3] Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding
682
           Techniques and Quantitative Performance Evaluation" Journal of
683
           Electronic Imaging, 13(1): 146-165
684
           :DOI:`10.1117/1.1631315`
685
    .. [4] ImageJ AutoThresholder code, http://fiji.sc/wiki/index.php/Auto_Threshold
686

687
    Examples
688
    --------
689
    >>> from skimage.data import camera
690
    >>> image = camera()
691
    >>> thresh = threshold_li(image)
692
    >>> binary = image > thresh
693
    """
694
    # Remove nan:
695
    image = image[~np.isnan(image)]
696
    if image.size == 0:
697
        return np.nan
698

699
    # Make sure image has more than one value; otherwise, return that value
700
    # This works even for np.inf
701
    if np.all(image == image.flat[0]):
702
        return image.flat[0]
703

704
    # At this point, the image only contains np.inf, -np.inf, or valid numbers
705
    image = image[np.isfinite(image)]
706
    # if there are no finite values in the image, return 0. This is because
707
    # at this point we *know* that there are *both* inf and -inf values,
708
    # because inf == inf evaluates to True. We might as well separate them.
709
    if image.size == 0:
710
        return 0.0
711

712
    # Li's algorithm requires positive image (because of log(mean))
713
    image_min = np.min(image)
714
    image -= image_min
715
    if image.dtype.kind in 'iu':
716
        tolerance = tolerance or 0.5
717
    else:
718
        tolerance = tolerance or np.min(np.diff(np.unique(image))) / 2
719

720
    # Initial estimate for iteration. See "initial_guess" in the parameter list
721
    if initial_guess is None:
722
        t_next = np.mean(image)
723
    elif callable(initial_guess):
724
        t_next = initial_guess(image)
725
    elif np.isscalar(initial_guess):  # convert to new, positive image range
726
        t_next = initial_guess - float(image_min)
727
        image_max = np.max(image) + image_min
728
        if not 0 < t_next < np.max(image):
729
            msg = (
730
                f'The initial guess for threshold_li must be within the '
731
                f'range of the image. Got {initial_guess} for image min '
732
                f'{image_min} and max {image_max}.'
733
            )
734
            raise ValueError(msg)
735
        t_next = image.dtype.type(t_next)
736
    else:
737
        raise TypeError(
738
            'Incorrect type for `initial_guess`; should be '
739
            'a floating point value, or a function mapping an '
740
            'array to a floating point value.'
741
        )
742

743
    # initial value for t_curr must be different from t_next by at
744
    # least the tolerance. Since the image is positive, we ensure this
745
    # by setting to a large-enough negative number
746
    t_curr = -2 * tolerance
747

748
    # Callback on initial iterations
749
    if iter_callback is not None:
750
        iter_callback(t_next + image_min)
751

752
    # Stop the iterations when the difference between the
753
    # new and old threshold values is less than the tolerance
754
    # or if the background mode has only one value left,
755
    # since log(0) is not defined.
756

757
    if image.dtype.kind in 'iu':
758
        hist, bin_centers = histogram(image.reshape(-1), source_range='image')
759
        hist = hist.astype('float32', copy=False)
760
        while abs(t_next - t_curr) > tolerance:
761
            t_curr = t_next
762
            foreground = bin_centers > t_curr
763
            background = ~foreground
764

765
            mean_fore = np.average(bin_centers[foreground], weights=hist[foreground])
766
            mean_back = np.average(bin_centers[background], weights=hist[background])
767

768
            if mean_back == 0:
769
                break
770

771
            t_next = (mean_back - mean_fore) / (np.log(mean_back) - np.log(mean_fore))
772

773
            if iter_callback is not None:
774
                iter_callback(t_next + image_min)
775

776
    else:
777
        while abs(t_next - t_curr) > tolerance:
778
            t_curr = t_next
779
            foreground = image > t_curr
780
            mean_fore = np.mean(image[foreground])
781
            mean_back = np.mean(image[~foreground])
782

783
            if mean_back == 0.0:
784
                break
785

786
            t_next = (mean_back - mean_fore) / (np.log(mean_back) - np.log(mean_fore))
787

788
            if iter_callback is not None:
789
                iter_callback(t_next + image_min)
790

791
    threshold = t_next + image_min
792
    return threshold
793

794

795
def threshold_minimum(image=None, nbins=256, max_num_iter=10000, *, hist=None):
796
    """Return threshold value based on minimum method.
797

798
    The histogram of the input ``image`` is computed if not provided and
799
    smoothed until there are only two maxima. Then the minimum in between is
800
    the threshold value.
801

802
    Either image or hist must be provided. In case hist is given, the actual
803
    histogram of the image is ignored.
804

805
    Parameters
806
    ----------
807
    image : (M, N[, ...]) ndarray, optional
808
        Grayscale input image.
809
    nbins : int, optional
810
        Number of bins used to calculate histogram. This value is ignored for
811
        integer arrays.
812
    max_num_iter : int, optional
813
        Maximum number of iterations to smooth the histogram.
814
    hist : array, or 2-tuple of arrays, optional
815
        Histogram to determine the threshold from and a corresponding array
816
        of bin center intensities. Alternatively, only the histogram can be
817
        passed.
818

819
    Returns
820
    -------
821
    threshold : float
822
        Upper threshold value. All pixels with an intensity higher than
823
        this value are assumed to be foreground.
824

825
    Raises
826
    ------
827
    RuntimeError
828
        If unable to find two local maxima in the histogram or if the
829
        smoothing takes more than 1e4 iterations.
830

831
    References
832
    ----------
833
    .. [1] C. A. Glasbey, "An analysis of histogram-based thresholding
834
           algorithms," CVGIP: Graphical Models and Image Processing,
835
           vol. 55, pp. 532-537, 1993.
836
    .. [2] Prewitt, JMS & Mendelsohn, ML (1966), "The analysis of cell
837
           images", Annals of the New York Academy of Sciences 128: 1035-1053
838
           :DOI:`10.1111/j.1749-6632.1965.tb11715.x`
839

840
    Examples
841
    --------
842
    >>> from skimage.data import camera
843
    >>> image = camera()
844
    >>> thresh = threshold_minimum(image)
845
    >>> binary = image > thresh
846
    """
847

848
    def find_local_maxima_idx(hist):
849
        # We can't use scipy.signal.argrelmax
850
        # as it fails on plateaus
851
        maximum_idxs = list()
852
        direction = 1
853

854
        for i in range(hist.shape[0] - 1):
855
            if direction > 0:
856
                if hist[i + 1] < hist[i]:
857
                    direction = -1
858
                    maximum_idxs.append(i)
859
            else:
860
                if hist[i + 1] > hist[i]:
861
                    direction = 1
862

863
        return maximum_idxs
864

865
    counts, bin_centers = _validate_image_histogram(image, hist, nbins)
866

867
    smooth_hist = counts.astype('float32', copy=False)
868

869
    for counter in range(max_num_iter):
870
        smooth_hist = ndi.uniform_filter1d(smooth_hist, 3)
871
        maximum_idxs = find_local_maxima_idx(smooth_hist)
872
        if len(maximum_idxs) < 3:
873
            break
874

875
    if len(maximum_idxs) != 2:
876
        raise RuntimeError('Unable to find two maxima in histogram')
877
    elif counter == max_num_iter - 1:
878
        raise RuntimeError('Maximum iteration reached for histogram' 'smoothing')
879

880
    # Find the lowest point between the maxima
881
    threshold_idx = np.argmin(smooth_hist[maximum_idxs[0] : maximum_idxs[1] + 1])
882

883
    return bin_centers[maximum_idxs[0] + threshold_idx]
884

885

886
def threshold_mean(image):
887
    """Return threshold value based on the mean of grayscale values.
888

889
    Parameters
890
    ----------
891
    image : (M, N[, ...]) ndarray
892
        Grayscale input image.
893

894
    Returns
895
    -------
896
    threshold : float
897
        Upper threshold value. All pixels with an intensity higher than
898
        this value are assumed to be foreground.
899

900
    References
901
    ----------
902
    .. [1] C. A. Glasbey, "An analysis of histogram-based thresholding
903
        algorithms," CVGIP: Graphical Models and Image Processing,
904
        vol. 55, pp. 532-537, 1993.
905
        :DOI:`10.1006/cgip.1993.1040`
906

907
    Examples
908
    --------
909
    >>> from skimage.data import camera
910
    >>> image = camera()
911
    >>> thresh = threshold_mean(image)
912
    >>> binary = image > thresh
913
    """
914
    return np.mean(image)
915

916

917
def threshold_triangle(image, nbins=256):
918
    """Return threshold value based on the triangle algorithm.
919

920
    Parameters
921
    ----------
922
    image : (M, N[, ...]) ndarray
923
        Grayscale input image.
924
    nbins : int, optional
925
        Number of bins used to calculate histogram. This value is ignored for
926
        integer arrays.
927

928
    Returns
929
    -------
930
    threshold : float
931
        Upper threshold value. All pixels with an intensity higher than
932
        this value are assumed to be foreground.
933

934
    References
935
    ----------
936
    .. [1] Zack, G. W., Rogers, W. E. and Latt, S. A., 1977,
937
       Automatic Measurement of Sister Chromatid Exchange Frequency,
938
       Journal of Histochemistry and Cytochemistry 25 (7), pp. 741-753
939
       :DOI:`10.1177/25.7.70454`
940
    .. [2] ImageJ AutoThresholder code,
941
       http://fiji.sc/wiki/index.php/Auto_Threshold
942

943
    Examples
944
    --------
945
    >>> from skimage.data import camera
946
    >>> image = camera()
947
    >>> thresh = threshold_triangle(image)
948
    >>> binary = image > thresh
949
    """
950
    # nbins is ignored for integer arrays
951
    # so, we recalculate the effective nbins.
952
    hist, bin_centers = histogram(image.reshape(-1), nbins, source_range='image')
953
    nbins = len(hist)
954

955
    # Find peak, lowest and highest gray levels.
956
    arg_peak_height = np.argmax(hist)
957
    peak_height = hist[arg_peak_height]
958
    arg_low_level, arg_high_level = np.flatnonzero(hist)[[0, -1]]
959

960
    if arg_low_level == arg_high_level:
961
        # Image has constant intensity.
962
        return image.ravel()[0]
963

964
    # Flip is True if left tail is shorter.
965
    flip = arg_peak_height - arg_low_level < arg_high_level - arg_peak_height
966
    if flip:
967
        hist = hist[::-1]
968
        arg_low_level = nbins - arg_high_level - 1
969
        arg_peak_height = nbins - arg_peak_height - 1
970

971
    # If flip == True, arg_high_level becomes incorrect
972
    # but we don't need it anymore.
973
    del arg_high_level
974

975
    # Set up the coordinate system.
976
    width = arg_peak_height - arg_low_level
977
    x1 = np.arange(width)
978
    y1 = hist[x1 + arg_low_level]
979

980
    # Normalize.
981
    norm = np.sqrt(peak_height**2 + width**2)
982
    peak_height /= norm
983
    width /= norm
984

985
    # Maximize the length.
986
    # The ImageJ implementation includes an additional constant when calculating
987
    # the length, but here we omit it as it does not affect the location of the
988
    # minimum.
989
    length = peak_height * x1 - width * y1
990
    arg_level = np.argmax(length) + arg_low_level
991

992
    if flip:
993
        arg_level = nbins - arg_level - 1
994

995
    return bin_centers[arg_level]
996

997

998
def _mean_std(image, w):
999
    """Return local mean and standard deviation of each pixel using a
1000
    neighborhood defined by a rectangular window size ``w``.
1001
    The algorithm uses integral images to speedup computation. This is
1002
    used by :func:`threshold_niblack` and :func:`threshold_sauvola`.
1003

1004
    Parameters
1005
    ----------
1006
    image : (M, N[, ...]) ndarray
1007
        Grayscale input image.
1008
    w : int, or iterable of int
1009
        Window size specified as a single odd integer (3, 5, 7, …),
1010
        or an iterable of length ``image.ndim`` containing only odd
1011
        integers (e.g. ``(1, 5, 5)``).
1012

1013
    Returns
1014
    -------
1015
    m : ndarray of float, same shape as ``image``
1016
        Local mean of the image.
1017
    s : ndarray of float, same shape as ``image``
1018
        Local standard deviation of the image.
1019

1020
    References
1021
    ----------
1022
    .. [1] F. Shafait, D. Keysers, and T. M. Breuel, "Efficient
1023
           implementation of local adaptive thresholding techniques
1024
           using integral images." in Document Recognition and
1025
           Retrieval XV, (San Jose, USA), Jan. 2008.
1026
           :DOI:`10.1117/12.767755`
1027
    """
1028

1029
    if not isinstance(w, Iterable):
1030
        w = (w,) * image.ndim
1031
    _validate_window_size(w)
1032

1033
    float_dtype = _supported_float_type(image.dtype)
1034
    pad_width = tuple((k // 2 + 1, k // 2) for k in w)
1035
    padded = np.pad(image.astype(float_dtype, copy=False), pad_width, mode='reflect')
1036

1037
    # Note: keep float64 integral images for accuracy. Outputs of
1038
    # _correlate_sparse can later be safely cast to float_dtype
1039
    integral = integral_image(padded, dtype=np.float64)
1040
    padded *= padded
1041
    integral_sq = integral_image(padded, dtype=np.float64)
1042

1043
    # Create lists of non-zero kernel indices and values
1044
    kernel_indices = list(itertools.product(*tuple([(0, _w) for _w in w])))
1045
    kernel_values = [
1046
        (-1) ** (image.ndim % 2 != np.sum(indices) % 2) for indices in kernel_indices
1047
    ]
1048

1049
    total_window_size = math.prod(w)
1050
    kernel_shape = tuple(_w + 1 for _w in w)
1051
    m = _correlate_sparse(integral, kernel_shape, kernel_indices, kernel_values)
1052
    m = m.astype(float_dtype, copy=False)
1053
    m /= total_window_size
1054
    g2 = _correlate_sparse(integral_sq, kernel_shape, kernel_indices, kernel_values)
1055
    g2 = g2.astype(float_dtype, copy=False)
1056
    g2 /= total_window_size
1057
    # Note: we use np.clip because g2 is not guaranteed to be greater than
1058
    # m*m when floating point error is considered
1059
    s = np.sqrt(np.clip(g2 - m * m, 0, None))
1060
    return m, s
1061

1062

1063
def threshold_niblack(image, window_size=15, k=0.2):
1064
    """Applies Niblack local threshold to an array.
1065

1066
    A threshold T is calculated for every pixel in the image using the
1067
    following formula::
1068

1069
        T = m(x,y) - k * s(x,y)
1070

1071
    where m(x,y) and s(x,y) are the mean and standard deviation of
1072
    pixel (x,y) neighborhood defined by a rectangular window with size w
1073
    times w centered around the pixel. k is a configurable parameter
1074
    that weights the effect of standard deviation.
1075

1076
    Parameters
1077
    ----------
1078
    image : (M, N[, ...]) ndarray
1079
        Grayscale input image.
1080
    window_size : int, or iterable of int, optional
1081
        Window size specified as a single odd integer (3, 5, 7, …),
1082
        or an iterable of length ``image.ndim`` containing only odd
1083
        integers (e.g. ``(1, 5, 5)``).
1084
    k : float, optional
1085
        Value of parameter k in threshold formula.
1086

1087
    Returns
1088
    -------
1089
    threshold : (M, N[, ...]) ndarray
1090
        Threshold mask. All pixels with an intensity higher than
1091
        this value are assumed to be foreground.
1092

1093
    Notes
1094
    -----
1095
    This algorithm is originally designed for text recognition.
1096

1097
    The Bradley threshold is a particular case of the Niblack
1098
    one, being equivalent to
1099

1100
    >>> from skimage import data
1101
    >>> image = data.page()
1102
    >>> q = 1
1103
    >>> threshold_image = threshold_niblack(image, k=0) * q
1104

1105
    for some value ``q``. By default, Bradley and Roth use ``q=1``.
1106

1107

1108
    References
1109
    ----------
1110
    .. [1] W. Niblack, An introduction to Digital Image Processing,
1111
           Prentice-Hall, 1986.
1112
    .. [2] D. Bradley and G. Roth, "Adaptive thresholding using Integral
1113
           Image", Journal of Graphics Tools 12(2), pp. 13-21, 2007.
1114
           :DOI:`10.1080/2151237X.2007.10129236`
1115

1116
    Examples
1117
    --------
1118
    >>> from skimage import data
1119
    >>> image = data.page()
1120
    >>> threshold_image = threshold_niblack(image, window_size=7, k=0.1)
1121
    """
1122
    m, s = _mean_std(image, window_size)
1123
    return m - k * s
1124

1125

1126
def threshold_sauvola(image, window_size=15, k=0.2, r=None):
1127
    """Applies Sauvola local threshold to an array. Sauvola is a
1128
    modification of Niblack technique.
1129

1130
    In the original method a threshold T is calculated for every pixel
1131
    in the image using the following formula::
1132

1133
        T = m(x,y) * (1 + k * ((s(x,y) / R) - 1))
1134

1135
    where m(x,y) and s(x,y) are the mean and standard deviation of
1136
    pixel (x,y) neighborhood defined by a rectangular window with size w
1137
    times w centered around the pixel. k is a configurable parameter
1138
    that weights the effect of standard deviation.
1139
    R is the maximum standard deviation of a grayscale image.
1140

1141
    Parameters
1142
    ----------
1143
    image : (M, N[, ...]) ndarray
1144
        Grayscale input image.
1145
    window_size : int, or iterable of int, optional
1146
        Window size specified as a single odd integer (3, 5, 7, …),
1147
        or an iterable of length ``image.ndim`` containing only odd
1148
        integers (e.g. ``(1, 5, 5)``).
1149
    k : float, optional
1150
        Value of the positive parameter k.
1151
    r : float, optional
1152
        Value of R, the dynamic range of standard deviation.
1153
        If None, set to the half of the image dtype range.
1154

1155
    Returns
1156
    -------
1157
    threshold : (M, N[, ...]) ndarray
1158
        Threshold mask. All pixels with an intensity higher than
1159
        this value are assumed to be foreground.
1160

1161
    Notes
1162
    -----
1163
    This algorithm is originally designed for text recognition.
1164

1165
    References
1166
    ----------
1167
    .. [1] J. Sauvola and M. Pietikainen, "Adaptive document image
1168
           binarization," Pattern Recognition 33(2),
1169
           pp. 225-236, 2000.
1170
           :DOI:`10.1016/S0031-3203(99)00055-2`
1171

1172
    Examples
1173
    --------
1174
    >>> from skimage import data
1175
    >>> image = data.page()
1176
    >>> t_sauvola = threshold_sauvola(image, window_size=15, k=0.2)
1177
    >>> binary_image = image > t_sauvola
1178
    """
1179
    if r is None:
1180
        imin, imax = dtype_limits(image, clip_negative=False)
1181
        r = 0.5 * (imax - imin)
1182
    m, s = _mean_std(image, window_size)
1183
    return m * (1 + k * ((s / r) - 1))
1184

1185

1186
def apply_hysteresis_threshold(image, low, high):
1187
    """Apply hysteresis thresholding to ``image``.
1188

1189
    This algorithm finds regions where ``image`` is greater than ``high``
1190
    OR ``image`` is greater than ``low`` *and* that region is connected to
1191
    a region greater than ``high``.
1192

1193
    Parameters
1194
    ----------
1195
    image : (M[, ...]) ndarray
1196
        Grayscale input image.
1197
    low : float, or array of same shape as ``image``
1198
        Lower threshold.
1199
    high : float, or array of same shape as ``image``
1200
        Higher threshold.
1201

1202
    Returns
1203
    -------
1204
    thresholded : (M[, ...]) array of bool
1205
        Array in which ``True`` indicates the locations where ``image``
1206
        was above the hysteresis threshold.
1207

1208
    Examples
1209
    --------
1210
    >>> image = np.array([1, 2, 3, 2, 1, 2, 1, 3, 2])
1211
    >>> apply_hysteresis_threshold(image, 1.5, 2.5).astype(int)
1212
    array([0, 1, 1, 1, 0, 0, 0, 1, 1])
1213

1214
    References
1215
    ----------
1216
    .. [1] J. Canny. A computational approach to edge detection.
1217
           IEEE Transactions on Pattern Analysis and Machine Intelligence.
1218
           1986; vol. 8, pp.679-698.
1219
           :DOI:`10.1109/TPAMI.1986.4767851`
1220
    """
1221
    low = np.clip(low, a_min=None, a_max=high)  # ensure low always below high
1222
    mask_low = image > low
1223
    mask_high = image > high
1224
    # Connected components of mask_low
1225
    labels_low, num_labels = ndi.label(mask_low)
1226
    # Check which connected components contain pixels from mask_high
1227
    sums = ndi.sum(mask_high, labels_low, np.arange(num_labels + 1))
1228
    connected_to_high = sums > 0
1229
    thresholded = connected_to_high[labels_low]
1230
    return thresholded
1231

1232

1233
def threshold_multiotsu(image=None, classes=3, nbins=256, *, hist=None):
1234
    r"""Generate `classes`-1 threshold values to divide gray levels in `image`,
1235
    following Otsu's method for multiple classes.
1236

1237
    The threshold values are chosen to maximize the total sum of pairwise
1238
    variances between the thresholded graylevel classes. See Notes and [1]_
1239
    for more details.
1240

1241
    Either image or hist must be provided. If hist is provided, the actual
1242
    histogram of the image is ignored.
1243

1244
    Parameters
1245
    ----------
1246
    image : (M, N[, ...]) ndarray, optional
1247
        Grayscale input image.
1248
    classes : int, optional
1249
        Number of classes to be thresholded, i.e. the number of resulting
1250
        regions.
1251
    nbins : int, optional
1252
        Number of bins used to calculate the histogram. This value is ignored
1253
        for integer arrays.
1254
    hist : array, or 2-tuple of arrays, optional
1255
        Histogram from which to determine the threshold, and optionally a
1256
        corresponding array of bin center intensities. If no hist provided,
1257
        this function will compute it from the image (see notes).
1258

1259
    Returns
1260
    -------
1261
    thresh : array
1262
        Array containing the threshold values for the desired classes.
1263

1264
    Raises
1265
    ------
1266
    ValueError
1267
         If ``image`` contains less grayscale value then the desired
1268
         number of classes.
1269

1270
    Notes
1271
    -----
1272
    This implementation relies on a Cython function whose complexity
1273
    is :math:`O\left(\frac{Ch^{C-1}}{(C-1)!}\right)`, where :math:`h`
1274
    is the number of histogram bins and :math:`C` is the number of
1275
    classes desired.
1276

1277
    If no hist is given, this function will make use of
1278
    `skimage.exposure.histogram`, which behaves differently than
1279
    `np.histogram`. While both allowed, use the former for consistent
1280
    behaviour.
1281

1282
    The input image must be grayscale.
1283

1284
    References
1285
    ----------
1286
    .. [1] Liao, P-S., Chen, T-S. and Chung, P-C., "A fast algorithm for
1287
           multilevel thresholding", Journal of Information Science and
1288
           Engineering 17 (5): 713-727, 2001. Available at:
1289
           <https://ftp.iis.sinica.edu.tw/JISE/2001/200109_01.pdf>
1290
           :DOI:`10.6688/JISE.2001.17.5.1`
1291
    .. [2] Tosa, Y., "Multi-Otsu Threshold", a java plugin for ImageJ.
1292
           Available at:
1293
           <http://imagej.net/plugins/download/Multi_OtsuThreshold.java>
1294

1295
    Examples
1296
    --------
1297
    >>> from skimage.color import label2rgb
1298
    >>> from skimage import data
1299
    >>> image = data.camera()
1300
    >>> thresholds = threshold_multiotsu(image)
1301
    >>> regions = np.digitize(image, bins=thresholds)
1302
    >>> regions_colorized = label2rgb(regions)
1303
    """
1304
    if image is not None and image.ndim > 2 and image.shape[-1] in (3, 4):
1305
        warn(
1306
            f'threshold_multiotsu is expected to work correctly only for '
1307
            f'grayscale images; image shape {image.shape} looks like '
1308
            f'that of an RGB image.'
1309
        )
1310

1311
    # calculating the histogram and the probability of each gray level.
1312
    prob, bin_centers = _validate_image_histogram(image, hist, nbins, normalize=True)
1313
    prob = prob.astype('float32', copy=False)
1314

1315
    nvalues = np.count_nonzero(prob)
1316
    if nvalues < classes:
1317
        msg = (
1318
            f'After discretization into bins, the input image has '
1319
            f'only {nvalues} different values. It cannot be thresholded '
1320
            f'in {classes} classes. If there are more unique values '
1321
            f'before discretization, try increasing the number of bins '
1322
            f'(`nbins`).'
1323
        )
1324
        raise ValueError(msg)
1325
    elif nvalues == classes:
1326
        thresh_idx = np.flatnonzero(prob)[:-1]
1327
    else:
1328
        # Get threshold indices
1329
        try:
1330
            thresh_idx = _get_multiotsu_thresh_indices_lut(prob, classes - 1)
1331
        except MemoryError:
1332
            # Don't use LUT if the number of bins is too large (if the
1333
            # image is uint16 for example): in this case, the
1334
            # allocated memory is too large.
1335
            thresh_idx = _get_multiotsu_thresh_indices(prob, classes - 1)
1336

1337
    thresh = bin_centers[thresh_idx]
1338

1339
    return thresh
1340