scikit-image

1
"""
2

3
General Description
4
-------------------
5

6
These filters compute the local histogram at each pixel, using a sliding window
7
similar to the method described in [1]_. A histogram is built using a moving
8
window in order to limit redundant computation. The moving window follows a
9
snake-like path:
10

11
...------------------------↘
12
↙--------------------------↙
13
↘--------------------------...
14

15
The local histogram is updated at each pixel as the footprint window
16
moves by, i.e. only those pixels entering and leaving the footprint
17
update the local histogram. The histogram size is 8-bit (256 bins) for 8-bit
18
images and 2- to 16-bit for 16-bit images depending on the maximum value of the
19
image.
20

21
The filter is applied up to the image border, the neighborhood used is
22
adjusted accordingly. The user may provide a mask image (same size as input
23
image) where non zero values are the part of the image participating in the
24
histogram computation. By default the entire image is filtered.
25

26
This implementation outperforms :func:`skimage.morphology.dilation`
27
for large footprints.
28

29
Input images will be cast in unsigned 8-bit integer or unsigned 16-bit integer
30
if necessary. The number of histogram bins is then determined from the maximum
31
value present in the image. Eventually, the output image is cast in the input
32
dtype, or the `output_dtype` if set.
33

34
To do
35
-----
36

37
* add simple examples, adapt documentation on existing examples
38
* add/check existing doc
39
* adapting tests for each type of filter
40

41

42
References
43
----------
44

45
.. [1] Huang, T. ,Yang, G. ;  Tang, G.. "A fast two-dimensional
46
       median filtering algorithm", IEEE Transactions on Acoustics, Speech and
47
       Signal Processing, Feb 1979. Volume: 27 , Issue: 1, Page(s): 13 - 18.
48

49
"""
50

51
import numpy as np
52
from scipy import ndimage as ndi
53

54
from ..._shared.utils import check_nD, warn
55
from ...morphology.footprints import _footprint_is_sequence
56
from ...util import img_as_ubyte
57
from . import generic_cy
58

59

60
__all__ = [
61
    'autolevel',
62
    'equalize',
63
    'gradient',
64
    'maximum',
65
    'mean',
66
    'geometric_mean',
67
    'subtract_mean',
68
    'median',
69
    'minimum',
70
    'modal',
71
    'enhance_contrast',
72
    'pop',
73
    'threshold',
74
    'noise_filter',
75
    'entropy',
76
    'otsu',
77
]
78

79

80
def _preprocess_input(
81
    image,
82
    footprint=None,
83
    out=None,
84
    mask=None,
85
    out_dtype=None,
86
    pixel_size=1,
87
    shift_x=None,
88
    shift_y=None,
89
):
90
    """Preprocess and verify input for filters.rank methods.
91

92
    Parameters
93
    ----------
94
    image : 2-D array (integer or float)
95
        Input image.
96
    footprint : 2-D array (integer or float), optional
97
        The neighborhood expressed as a 2-D array of 1's and 0's.
98
    out : 2-D array (integer or float), optional
99
        If None, a new array is allocated.
100
    mask : ndarray (integer or float), optional
101
        Mask array that defines (>0) area of the image included in the local
102
        neighborhood. If None, the complete image is used (default).
103
    out_dtype : data-type, optional
104
        Desired output data-type. Default is None, which means we cast output
105
        in input dtype.
106
    pixel_size : int, optional
107
        Dimension of each pixel. Default value is 1.
108
    shift_x, shift_y : int, optional
109
        Offset added to the footprint center point. Shift is bounded to the
110
        footprint size (center must be inside of the given footprint).
111

112
    Returns
113
    -------
114
    image : 2-D array (np.uint8 or np.uint16)
115
    footprint : 2-D array (np.uint8)
116
        The neighborhood expressed as a binary 2-D array.
117
    out : 3-D array (same dtype out_dtype or as input)
118
        Output array. The two first dimensions are the spatial ones, the third
119
        one is the pixel vector (length 1 by default).
120
    mask : 2-D array (np.uint8)
121
        Mask array that defines (>0) area of the image included in the local
122
        neighborhood.
123
    n_bins : int
124
        Number of histogram bins.
125

126
    """
127
    check_nD(image, 2)
128
    input_dtype = image.dtype
129
    if input_dtype in (bool, bool) or out_dtype in (bool, bool):
130
        raise ValueError('dtype cannot be bool.')
131
    if input_dtype not in (np.uint8, np.uint16):
132
        message = (
133
            f'Possible precision loss converting image of type '
134
            f'{input_dtype} to uint8 as required by rank filters. '
135
            f'Convert manually using skimage.util.img_as_ubyte to '
136
            f'silence this warning.'
137
        )
138
        warn(message, stacklevel=5)
139
        image = img_as_ubyte(image)
140

141
    if _footprint_is_sequence(footprint):
142
        raise ValueError(
143
            "footprint sequences are not currently supported by rank filters"
144
        )
145

146
    footprint = np.ascontiguousarray(img_as_ubyte(footprint > 0))
147
    if footprint.ndim != image.ndim:
148
        raise ValueError('Image dimensions and neighborhood dimensions' 'do not match')
149

150
    image = np.ascontiguousarray(image)
151

152
    if mask is not None:
153
        mask = img_as_ubyte(mask)
154
        mask = np.ascontiguousarray(mask)
155

156
    if image is out:
157
        raise NotImplementedError("Cannot perform rank operation in place.")
158

159
    if out is None:
160
        if out_dtype is None:
161
            out_dtype = image.dtype
162
        out = np.empty(image.shape + (pixel_size,), dtype=out_dtype)
163
    else:
164
        if len(out.shape) == 2:
165
            out = out.reshape(out.shape + (pixel_size,))
166

167
    if image.dtype in (np.uint8, np.int8):
168
        n_bins = 256
169
    else:
170
        # Convert to a Python int to avoid the potential overflow when we add
171
        # 1 to the maximum of the image.
172
        n_bins = int(max(3, image.max())) + 1
173

174
    if n_bins > 2**10:
175
        warn(
176
            f'Bad rank filter performance is expected due to a '
177
            f'large number of bins ({n_bins}), equivalent to an approximate '
178
            f'bitdepth of {np.log2(n_bins):.1f}.',
179
            stacklevel=2,
180
        )
181

182
    for name, value in zip(("shift_x", "shift_y"), (shift_x, shift_y)):
183
        if np.dtype(type(value)) == bool:
184
            warn(
185
                f"Paramter `{name}` is boolean and will be interpreted as int. "
186
                "This is not officially supported, use int instead.",
187
                category=UserWarning,
188
                stacklevel=4,
189
            )
190

191
    return image, footprint, out, mask, n_bins
192

193

194
def _handle_input_3D(
195
    image,
196
    footprint=None,
197
    out=None,
198
    mask=None,
199
    out_dtype=None,
200
    pixel_size=1,
201
    shift_x=None,
202
    shift_y=None,
203
    shift_z=None,
204
):
205
    """Preprocess and verify input for filters.rank methods.
206

207
    Parameters
208
    ----------
209
    image : 3-D array (integer or float)
210
        Input image.
211
    footprint : 3-D array (integer or float), optional
212
        The neighborhood expressed as a 3-D array of 1's and 0's.
213
    out : 3-D array (integer or float), optional
214
        If None, a new array is allocated.
215
    mask : ndarray (integer or float), optional
216
        Mask array that defines (>0) area of the image included in the local
217
        neighborhood. If None, the complete image is used (default).
218
    out_dtype : data-type, optional
219
        Desired output data-type. Default is None, which means we cast output
220
        in input dtype.
221
    pixel_size : int, optional
222
        Dimension of each pixel. Default value is 1.
223
    shift_x, shift_y, shift_z : int, optional
224
        Offset added to the footprint center point. Shift is bounded to the
225
        footprint size (center must be inside of the given footprint).
226

227
    Returns
228
    -------
229
    image : 3-D array (np.uint8 or np.uint16)
230
    footprint : 3-D array (np.uint8)
231
        The neighborhood expressed as a binary 3-D array.
232
    out : 3-D array (same dtype out_dtype or as input)
233
        Output array. The two first dimensions are the spatial ones, the third
234
        one is the pixel vector (length 1 by default).
235
    mask : 3-D array (np.uint8)
236
        Mask array that defines (>0) area of the image included in the local
237
        neighborhood.
238
    n_bins : int
239
        Number of histogram bins.
240

241
    """
242
    check_nD(image, 3)
243
    if image.dtype not in (np.uint8, np.uint16):
244
        message = (
245
            f'Possible precision loss converting image of type '
246
            f'{image.dtype} to uint8 as required by rank filters. '
247
            f'Convert manually using skimage.util.img_as_ubyte to '
248
            f'silence this warning.'
249
        )
250
        warn(message, stacklevel=2)
251
        image = img_as_ubyte(image)
252

253
    footprint = np.ascontiguousarray(img_as_ubyte(footprint > 0))
254
    if footprint.ndim != image.ndim:
255
        raise ValueError('Image dimensions and neighborhood dimensions' 'do not match')
256
    image = np.ascontiguousarray(image)
257

258
    if mask is None:
259
        mask = np.ones(image.shape, dtype=np.uint8)
260
    else:
261
        mask = img_as_ubyte(mask)
262
        mask = np.ascontiguousarray(mask)
263

264
    if image is out:
265
        raise NotImplementedError("Cannot perform rank operation in place.")
266

267
    if out is None:
268
        if out_dtype is None:
269
            out_dtype = image.dtype
270
        out = np.empty(image.shape + (pixel_size,), dtype=out_dtype)
271
    else:
272
        out = out.reshape(out.shape + (pixel_size,))
273

274
    is_8bit = image.dtype in (np.uint8, np.int8)
275

276
    if is_8bit:
277
        n_bins = 256
278
    else:
279
        # Convert to a Python int to avoid the potential overflow when we add
280
        # 1 to the maximum of the image.
281
        n_bins = int(max(3, image.max())) + 1
282

283
    if n_bins > 2**10:
284
        warn(
285
            f'Bad rank filter performance is expected due to a '
286
            f'large number of bins ({n_bins}), equivalent to an approximate '
287
            f'bitdepth of {np.log2(n_bins):.1f}.',
288
            stacklevel=2,
289
        )
290

291
    for name, value in zip(
292
        ("shift_x", "shift_y", "shift_z"), (shift_x, shift_y, shift_z)
293
    ):
294
        if np.dtype(type(value)) == bool:
295
            warn(
296
                f"Parameter `{name}` is boolean and will be interpreted as int. "
297
                "This is not officially supported, use int instead.",
298
                category=UserWarning,
299
                stacklevel=4,
300
            )
301

302
    return image, footprint, out, mask, n_bins
303

304

305
def _apply_scalar_per_pixel(
306
    func, image, footprint, out, mask, shift_x, shift_y, out_dtype=None
307
):
308
    """Process the specific cython function to the image.
309

310
    Parameters
311
    ----------
312
    func : function
313
        Cython function to apply.
314
    image : 2-D array (integer or float)
315
        Input image.
316
    footprint : 2-D array (integer or float)
317
        The neighborhood expressed as a 2-D array of 1's and 0's.
318
    out : 2-D array (integer or float)
319
        If None, a new array is allocated.
320
    mask : ndarray (integer or float)
321
        Mask array that defines (>0) area of the image included in the local
322
        neighborhood. If None, the complete image is used (default).
323
    shift_x, shift_y : int
324
        Offset added to the footprint center point. Shift is bounded to the
325
        footprint sizes (center must be inside the given footprint).
326
    out_dtype : data-type, optional
327
        Desired output data-type. Default is None, which means we cast output
328
        in input dtype.
329

330
    """
331
    # preprocess and verify the input
332
    image, footprint, out, mask, n_bins = _preprocess_input(
333
        image, footprint, out, mask, out_dtype, shift_x=shift_x, shift_y=shift_y
334
    )
335

336
    # apply cython function
337
    func(
338
        image,
339
        footprint,
340
        shift_x=shift_x,
341
        shift_y=shift_y,
342
        mask=mask,
343
        out=out,
344
        n_bins=n_bins,
345
    )
346

347
    return np.squeeze(out, axis=-1)
348

349

350
def _apply_scalar_per_pixel_3D(
351
    func, image, footprint, out, mask, shift_x, shift_y, shift_z, out_dtype=None
352
):
353
    image, footprint, out, mask, n_bins = _handle_input_3D(
354
        image,
355
        footprint,
356
        out,
357
        mask,
358
        out_dtype,
359
        shift_x=shift_x,
360
        shift_y=shift_y,
361
        shift_z=shift_z,
362
    )
363

364
    func(
365
        image,
366
        footprint,
367
        shift_x=shift_x,
368
        shift_y=shift_y,
369
        shift_z=shift_z,
370
        mask=mask,
371
        out=out,
372
        n_bins=n_bins,
373
    )
374

375
    return out.reshape(out.shape[:3])
376

377

378
def _apply_vector_per_pixel(
379
    func, image, footprint, out, mask, shift_x, shift_y, out_dtype=None, pixel_size=1
380
):
381
    """
382

383
    Parameters
384
    ----------
385
    func : function
386
        Cython function to apply.
387
    image : 2-D array (integer or float)
388
        Input image.
389
    footprint : 2-D array (integer or float)
390
        The neighborhood expressed as a 2-D array of 1's and 0's.
391
    out : 2-D array (integer or float)
392
        If None, a new array is allocated.
393
    mask : ndarray (integer or float)
394
        Mask array that defines (>0) area of the image included in the local
395
        neighborhood. If None, the complete image is used (default).
396
    shift_x, shift_y : int
397
        Offset added to the footprint center point. Shift is bounded to the
398
        footprint sizes (center must be inside the given footprint).
399
    out_dtype : data-type, optional
400
        Desired output data-type. Default is None, which means we cast output
401
        in input dtype.
402
    pixel_size : int, optional
403
        Dimension of each pixel.
404

405
    Returns
406
    -------
407
    out : 3-D array with float dtype of dimensions (H,W,N), where (H,W) are
408
        the dimensions of the input image and N is n_bins or
409
        ``image.max() + 1`` if no value is provided as a parameter.
410
        Effectively, each pixel is a N-D feature vector that is the histogram.
411
        The sum of the elements in the feature vector will be 1, unless no
412
        pixels in the window were covered by both footprint and mask, in which
413
        case all elements will be 0.
414

415
    """
416
    # preprocess and verify the input
417
    image, footprint, out, mask, n_bins = _preprocess_input(
418
        image,
419
        footprint,
420
        out,
421
        mask,
422
        out_dtype,
423
        pixel_size,
424
        shift_x=shift_x,
425
        shift_y=shift_y,
426
    )
427

428
    # apply cython function
429
    func(
430
        image,
431
        footprint,
432
        shift_x=shift_x,
433
        shift_y=shift_y,
434
        mask=mask,
435
        out=out,
436
        n_bins=n_bins,
437
    )
438

439
    return out
440

441

442
def autolevel(image, footprint, out=None, mask=None, shift_x=0, shift_y=0, shift_z=0):
443
    """Auto-level image using local histogram.
444

445
    This filter locally stretches the histogram of gray values to cover the
446
    entire range of values from "white" to "black".
447

448
    Parameters
449
    ----------
450
    image : ([P,] M, N) ndarray (uint8, uint16)
451
        Input image.
452
    footprint : ndarray
453
        The neighborhood expressed as an ndarray of 1's and 0's.
454
    out : ([P,] M, N) array (same dtype as input)
455
        If None, a new array is allocated.
456
    mask : ndarray (integer or float), optional
457
        Mask array that defines (>0) area of the image included in the local
458
        neighborhood. If None, the complete image is used (default).
459
    shift_x, shift_y, shift_z : int
460
        Offset added to the footprint center point. Shift is bounded to the
461
        footprint sizes (center must be inside the given footprint).
462

463
    Returns
464
    -------
465
    out : ([P,] M, N) ndarray (same dtype as input image)
466
        Output image.
467

468
    Examples
469
    --------
470
    >>> from skimage import data
471
    >>> from skimage.morphology import disk, ball
472
    >>> from skimage.filters.rank import autolevel
473
    >>> import numpy as np
474
    >>> img = data.camera()
475
    >>> rng = np.random.default_rng()
476
    >>> volume = rng.integers(0, 255, size=(10,10,10), dtype=np.uint8)
477
    >>> auto = autolevel(img, disk(5))
478
    >>> auto_vol = autolevel(volume, ball(5))
479

480
    """
481

482
    np_image = np.asanyarray(image)
483
    if np_image.ndim == 2:
484
        return _apply_scalar_per_pixel(
485
            generic_cy._autolevel,
486
            image,
487
            footprint,
488
            out=out,
489
            mask=mask,
490
            shift_x=shift_x,
491
            shift_y=shift_y,
492
        )
493
    elif np_image.ndim == 3:
494
        return _apply_scalar_per_pixel_3D(
495
            generic_cy._autolevel_3D,
496
            image,
497
            footprint,
498
            out=out,
499
            mask=mask,
500
            shift_x=shift_x,
501
            shift_y=shift_y,
502
            shift_z=shift_z,
503
        )
504
    raise ValueError(f'`image` must have 2 or 3 dimensions, got {np_image.ndim}.')
505

506

507
def equalize(image, footprint, out=None, mask=None, shift_x=0, shift_y=0, shift_z=0):
508
    """Equalize image using local histogram.
509

510
    Parameters
511
    ----------
512
    image : ([P,] M, N) ndarray (uint8, uint16)
513
        Input image.
514
    footprint : ndarray
515
        The neighborhood expressed as an ndarray of 1's and 0's.
516
    out : ([P,] M, N) array (same dtype as input)
517
        If None, a new array is allocated.
518
    mask : ndarray (integer or float), optional
519
        Mask array that defines (>0) area of the image included in the local
520
        neighborhood. If None, the complete image is used (default).
521
    shift_x, shift_y, shift_z : int
522
        Offset added to the footprint center point. Shift is bounded to the
523
        footprint sizes (center must be inside the given footprint).
524

525
    Returns
526
    -------
527
    out : ([P,] M, N) ndarray (same dtype as input image)
528
        Output image.
529

530
    Examples
531
    --------
532
    >>> from skimage import data
533
    >>> from skimage.morphology import disk, ball
534
    >>> from skimage.filters.rank import equalize
535
    >>> import numpy as np
536
    >>> img = data.camera()
537
    >>> rng = np.random.default_rng()
538
    >>> volume = rng.integers(0, 255, size=(10,10,10), dtype=np.uint8)
539
    >>> equ = equalize(img, disk(5))
540
    >>> equ_vol = equalize(volume, ball(5))
541

542
    """
543

544
    np_image = np.asanyarray(image)
545
    if np_image.ndim == 2:
546
        return _apply_scalar_per_pixel(
547
            generic_cy._equalize,
548
            image,
549
            footprint,
550
            out=out,
551
            mask=mask,
552
            shift_x=shift_x,
553
            shift_y=shift_y,
554
        )
555
    elif np_image.ndim == 3:
556
        return _apply_scalar_per_pixel_3D(
557
            generic_cy._equalize_3D,
558
            image,
559
            footprint,
560
            out=out,
561
            mask=mask,
562
            shift_x=shift_x,
563
            shift_y=shift_y,
564
            shift_z=shift_z,
565
        )
566
    raise ValueError(f'`image` must have 2 or 3 dimensions, got {np_image.ndim}.')
567

568

569
def gradient(image, footprint, out=None, mask=None, shift_x=0, shift_y=0, shift_z=0):
570
    """Return local gradient of an image (i.e. local maximum - local minimum).
571

572
    Parameters
573
    ----------
574
    image : ([P,] M, N) ndarray (uint8, uint16)
575
        Input image.
576
    footprint : ndarray
577
        The neighborhood expressed as an ndarray of 1's and 0's.
578
    out : ([P,] M, N) array (same dtype as input)
579
        If None, a new array is allocated.
580
    mask : ndarray (integer or float), optional
581
        Mask array that defines (>0) area of the image included in the local
582
        neighborhood. If None, the complete image is used (default).
583
    shift_x, shift_y, shift_z : int
584
        Offset added to the footprint center point. Shift is bounded to the
585
        footprint sizes (center must be inside the given footprint).
586

587
    Returns
588
    -------
589
    out : ([P,] M, N) ndarray (same dtype as input image)
590
        Output image.
591

592
    Examples
593
    --------
594
    >>> from skimage import data
595
    >>> from skimage.morphology import disk, ball
596
    >>> from skimage.filters.rank import gradient
597
    >>> import numpy as np
598
    >>> img = data.camera()
599
    >>> rng = np.random.default_rng()
600
    >>> volume = rng.integers(0, 255, size=(10,10,10), dtype=np.uint8)
601
    >>> out = gradient(img, disk(5))
602
    >>> out_vol = gradient(volume, ball(5))
603

604
    """
605

606
    np_image = np.asanyarray(image)
607
    if np_image.ndim == 2:
608
        return _apply_scalar_per_pixel(
609
            generic_cy._gradient,
610
            image,
611
            footprint,
612
            out=out,
613
            mask=mask,
614
            shift_x=shift_x,
615
            shift_y=shift_y,
616
        )
617
    elif np_image.ndim == 3:
618
        return _apply_scalar_per_pixel_3D(
619
            generic_cy._gradient_3D,
620
            image,
621
            footprint,
622
            out=out,
623
            mask=mask,
624
            shift_x=shift_x,
625
            shift_y=shift_y,
626
            shift_z=shift_z,
627
        )
628
    raise ValueError(f'`image` must have 2 or 3 dimensions, got {np_image.ndim}.')
629

630

631
def maximum(image, footprint, out=None, mask=None, shift_x=0, shift_y=0, shift_z=0):
632
    """Return local maximum of an image.
633

634
    Parameters
635
    ----------
636
    image : ([P,] M, N) ndarray (uint8, uint16)
637
        Input image.
638
    footprint : ndarray
639
        The neighborhood expressed as an ndarray of 1's and 0's.
640
    out : ([P,] M, N) array (same dtype as input)
641
        If None, a new array is allocated.
642
    mask : ndarray (integer or float), optional
643
        Mask array that defines (>0) area of the image included in the local
644
        neighborhood. If None, the complete image is used (default).
645
    shift_x, shift_y, shift_z : int
646
        Offset added to the footprint center point. Shift is bounded to the
647
        footprint sizes (center must be inside the given footprint).
648

649
    Returns
650
    -------
651
    out : ([P,] M, N) ndarray (same dtype as input image)
652
        Output image.
653

654
    See also
655
    --------
656
    skimage.morphology.dilation
657

658
    Notes
659
    -----
660
    The lower algorithm complexity makes `skimage.filters.rank.maximum`
661
    more efficient for larger images and footprints.
662

663
    Examples
664
    --------
665
    >>> from skimage import data
666
    >>> from skimage.morphology import disk, ball
667
    >>> from skimage.filters.rank import maximum
668
    >>> import numpy as np
669
    >>> img = data.camera()
670
    >>> rng = np.random.default_rng()
671
    >>> volume = rng.integers(0, 255, size=(10,10,10), dtype=np.uint8)
672
    >>> out = maximum(img, disk(5))
673
    >>> out_vol = maximum(volume, ball(5))
674

675
    """
676

677
    np_image = np.asanyarray(image)
678
    if np_image.ndim == 2:
679
        return _apply_scalar_per_pixel(
680
            generic_cy._maximum,
681
            image,
682
            footprint,
683
            out=out,
684
            mask=mask,
685
            shift_x=shift_x,
686
            shift_y=shift_y,
687
        )
688
    elif np_image.ndim == 3:
689
        return _apply_scalar_per_pixel_3D(
690
            generic_cy._maximum_3D,
691
            image,
692
            footprint,
693
            out=out,
694
            mask=mask,
695
            shift_x=shift_x,
696
            shift_y=shift_y,
697
            shift_z=shift_z,
698
        )
699
    raise ValueError(f'`image` must have 2 or 3 dimensions, got {np_image.ndim}.')
700

701

702
def mean(image, footprint, out=None, mask=None, shift_x=0, shift_y=0, shift_z=0):
703
    """Return local mean of an image.
704

705
    Parameters
706
    ----------
707
    image : ([P,] M, N) ndarray (uint8, uint16)
708
        Input image.
709
    footprint : ndarray
710
        The neighborhood expressed as an ndarray of 1's and 0's.
711
    out : ([P,] M, N) array (same dtype as input)
712
        If None, a new array is allocated.
713
    mask : ndarray (integer or float), optional
714
        Mask array that defines (>0) area of the image included in the local
715
        neighborhood. If None, the complete image is used (default).
716
    shift_x, shift_y, shift_z : int
717
        Offset added to the footprint center point. Shift is bounded to the
718
        footprint sizes (center must be inside the given footprint).
719

720
    Returns
721
    -------
722
    out : ([P,] M, N) ndarray (same dtype as input image)
723
        Output image.
724

725
    Examples
726
    --------
727
    >>> from skimage import data
728
    >>> from skimage.morphology import disk, ball
729
    >>> from skimage.filters.rank import mean
730
    >>> import numpy as np
731
    >>> img = data.camera()
732
    >>> rng = np.random.default_rng()
733
    >>> volume = rng.integers(0, 255, size=(10,10,10), dtype=np.uint8)
734
    >>> avg = mean(img, disk(5))
735
    >>> avg_vol = mean(volume, ball(5))
736

737
    """
738

739
    np_image = np.asanyarray(image)
740
    if np_image.ndim == 2:
741
        return _apply_scalar_per_pixel(
742
            generic_cy._mean,
743
            image,
744
            footprint,
745
            out=out,
746
            mask=mask,
747
            shift_x=shift_x,
748
            shift_y=shift_y,
749
        )
750
    elif np_image.ndim == 3:
751
        return _apply_scalar_per_pixel_3D(
752
            generic_cy._mean_3D,
753
            image,
754
            footprint,
755
            out=out,
756
            mask=mask,
757
            shift_x=shift_x,
758
            shift_y=shift_y,
759
            shift_z=shift_z,
760
        )
761
    raise ValueError(f'`image` must have 2 or 3 dimensions, got {np_image.ndim}.')
762

763

764
def geometric_mean(
765
    image, footprint, out=None, mask=None, shift_x=0, shift_y=0, shift_z=0
766
):
767
    """Return local geometric mean of an image.
768

769
    Parameters
770
    ----------
771
    image : ([P,] M, N) ndarray (uint8, uint16)
772
        Input image.
773
    footprint : ndarray
774
        The neighborhood expressed as an ndarray of 1's and 0's.
775
    out : ([P,] M, N) array (same dtype as input)
776
        If None, a new array is allocated.
777
    mask : ndarray (integer or float), optional
778
        Mask array that defines (>0) area of the image included in the local
779
        neighborhood. If None, the complete image is used (default).
780
    shift_x, shift_y, shift_z : int
781
        Offset added to the footprint center point. Shift is bounded to the
782
        footprint sizes (center must be inside the given footprint).
783

784
    Returns
785
    -------
786
    out : ([P,] M, N) ndarray (same dtype as input image)
787
        Output image.
788

789
    Examples
790
    --------
791
    >>> from skimage import data
792
    >>> from skimage.morphology import disk, ball
793
    >>> from skimage.filters.rank import mean
794
    >>> import numpy as np
795
    >>> img = data.camera()
796
    >>> rng = np.random.default_rng()
797
    >>> volume = rng.integers(0, 255, size=(10,10,10), dtype=np.uint8)
798
    >>> avg = geometric_mean(img, disk(5))
799
    >>> avg_vol = geometric_mean(volume, ball(5))
800

801
    References
802
    ----------
803
    .. [1] Gonzalez, R. C. and Woods, R. E. "Digital Image Processing
804
           (3rd Edition)." Prentice-Hall Inc, 2006.
805

806
    """
807

808
    np_image = np.asanyarray(image)
809
    if np_image.ndim == 2:
810
        return _apply_scalar_per_pixel(
811
            generic_cy._geometric_mean,
812
            image,
813
            footprint,
814
            out=out,
815
            mask=mask,
816
            shift_x=shift_x,
817
            shift_y=shift_y,
818
        )
819
    elif np_image.ndim == 3:
820
        return _apply_scalar_per_pixel_3D(
821
            generic_cy._geometric_mean_3D,
822
            image,
823
            footprint,
824
            out=out,
825
            mask=mask,
826
            shift_x=shift_x,
827
            shift_y=shift_y,
828
            shift_z=shift_z,
829
        )
830
    raise ValueError(f'`image` must have 2 or 3 dimensions, got {np_image.ndim}.')
831

832

833
def subtract_mean(
834
    image, footprint, out=None, mask=None, shift_x=0, shift_y=0, shift_z=0
835
):
836
    """Return image subtracted from its local mean.
837

838
    Parameters
839
    ----------
840
    image : ([P,] M, N) ndarray (uint8, uint16)
841
        Input image.
842
    footprint : ndarray
843
        The neighborhood expressed as an ndarray of 1's and 0's.
844
    out : ([P,] M, N) array (same dtype as input)
845
        If None, a new array is allocated.
846
    mask : ndarray (integer or float), optional
847
        Mask array that defines (>0) area of the image included in the local
848
        neighborhood. If None, the complete image is used (default).
849
    shift_x, shift_y, shift_z : int
850
        Offset added to the footprint center point. Shift is bounded to the
851
        footprint sizes (center must be inside the given footprint).
852

853
    Returns
854
    -------
855
    out : ([P,] M, N) ndarray (same dtype as input image)
856
        Output image.
857

858
    Notes
859
    -----
860
    Subtracting the mean value may introduce underflow. To compensate
861
    this potential underflow, the obtained difference is downscaled by
862
    a factor of 2 and shifted by `n_bins / 2 - 1`, the median value of
863
    the local histogram (`n_bins = max(3, image.max()) +1` for 16-bits
864
    images and 256 otherwise).
865

866
    Examples
867
    --------
868
    >>> from skimage import data
869
    >>> from skimage.morphology import disk, ball
870
    >>> from skimage.filters.rank import subtract_mean
871
    >>> import numpy as np
872
    >>> img = data.camera()
873
    >>> rng = np.random.default_rng()
874
    >>> volume = rng.integers(0, 255, size=(10,10,10), dtype=np.uint8)
875
    >>> out = subtract_mean(img, disk(5))
876
    >>> out_vol = subtract_mean(volume, ball(5))
877

878
    """
879

880
    np_image = np.asanyarray(image)
881
    if np_image.ndim == 2:
882
        return _apply_scalar_per_pixel(
883
            generic_cy._subtract_mean,
884
            image,
885
            footprint,
886
            out=out,
887
            mask=mask,
888
            shift_x=shift_x,
889
            shift_y=shift_y,
890
        )
891
    elif np_image.ndim == 3:
892
        return _apply_scalar_per_pixel_3D(
893
            generic_cy._subtract_mean_3D,
894
            image,
895
            footprint,
896
            out=out,
897
            mask=mask,
898
            shift_x=shift_x,
899
            shift_y=shift_y,
900
            shift_z=shift_z,
901
        )
902
    raise ValueError(f'`image` must have 2 or 3 dimensions, got {np_image.ndim}.')
903

904

905
def median(
906
    image,
907
    footprint=None,
908
    out=None,
909
    mask=None,
910
    shift_x=0,
911
    shift_y=0,
912
    shift_z=0,
913
):
914
    """Return local median of an image.
915

916
    Parameters
917
    ----------
918
    image : ([P,] M, N) ndarray (uint8, uint16)
919
        Input image.
920
    footprint : ndarray
921
        The neighborhood expressed as an ndarray of 1's and 0's. If None, a
922
        full square of size 3 is used.
923
    out : ([P,] M, N) array (same dtype as input)
924
        If None, a new array is allocated.
925
    mask : ndarray (integer or float), optional
926
        Mask array that defines (>0) area of the image included in the local
927
        neighborhood. If None, the complete image is used (default).
928
    shift_x, shift_y, shift_z : int
929
        Offset added to the footprint center point. Shift is bounded to the
930
        footprint sizes (center must be inside the given footprint).
931

932
    Returns
933
    -------
934
    out : ([P,] M, N) ndarray (same dtype as input image)
935
        Output image.
936

937
    See also
938
    --------
939
    skimage.filters.median : Implementation of a median filtering which handles
940
        images with floating precision.
941

942
    Examples
943
    --------
944
    >>> from skimage import data
945
    >>> from skimage.morphology import disk, ball
946
    >>> from skimage.filters.rank import median
947
    >>> import numpy as np
948
    >>> img = data.camera()
949
    >>> rng = np.random.default_rng()
950
    >>> volume = rng.integers(0, 255, size=(10,10,10), dtype=np.uint8)
951
    >>> med = median(img, disk(5))
952
    >>> med_vol = median(volume, ball(5))
953

954
    """
955

956
    np_image = np.asanyarray(image)
957
    if footprint is None:
958
        footprint = ndi.generate_binary_structure(image.ndim, image.ndim)
959
    if np_image.ndim == 2:
960
        return _apply_scalar_per_pixel(
961
            generic_cy._median,
962
            image,
963
            footprint,
964
            out=out,
965
            mask=mask,
966
            shift_x=shift_x,
967
            shift_y=shift_y,
968
        )
969
    elif np_image.ndim == 3:
970
        return _apply_scalar_per_pixel_3D(
971
            generic_cy._median_3D,
972
            image,
973
            footprint,
974
            out=out,
975
            mask=mask,
976
            shift_x=shift_x,
977
            shift_y=shift_y,
978
            shift_z=shift_z,
979
        )
980
    raise ValueError(f'`image` must have 2 or 3 dimensions, got {np_image.ndim}.')
981

982

983
def minimum(image, footprint, out=None, mask=None, shift_x=0, shift_y=0, shift_z=0):
984
    """Return local minimum of an image.
985

986
    Parameters
987
    ----------
988
    image : ([P,] M, N) ndarray (uint8, uint16)
989
        Input image.
990
    footprint : ndarray
991
        The neighborhood expressed as an ndarray of 1's and 0's.
992
    out : ([P,] M, N) array (same dtype as input)
993
        If None, a new array is allocated.
994
    mask : ndarray (integer or float), optional
995
        Mask array that defines (>0) area of the image included in the local
996
        neighborhood. If None, the complete image is used (default).
997
    shift_x, shift_y, shift_z : int
998
        Offset added to the footprint center point. Shift is bounded to the
999
        footprint sizes (center must be inside the given footprint).
1000

1001
    Returns
1002
    -------
1003
    out : ([P,] M, N) ndarray (same dtype as input image)
1004
        Output image.
1005

1006
    See also
1007
    --------
1008
    skimage.morphology.erosion
1009

1010
    Notes
1011
    -----
1012
    The lower algorithm complexity makes `skimage.filters.rank.minimum` more
1013
    efficient for larger images and footprints.
1014

1015
    Examples
1016
    --------
1017
    >>> from skimage import data
1018
    >>> from skimage.morphology import disk, ball
1019
    >>> from skimage.filters.rank import minimum
1020
    >>> import numpy as np
1021
    >>> img = data.camera()
1022
    >>> rng = np.random.default_rng()
1023
    >>> volume = rng.integers(0, 255, size=(10,10,10), dtype=np.uint8)
1024
    >>> out = minimum(img, disk(5))
1025
    >>> out_vol = minimum(volume, ball(5))
1026

1027
    """
1028

1029
    np_image = np.asanyarray(image)
1030
    if np_image.ndim == 2:
1031
        return _apply_scalar_per_pixel(
1032
            generic_cy._minimum,
1033
            image,
1034
            footprint,
1035
            out=out,
1036
            mask=mask,
1037
            shift_x=shift_x,
1038
            shift_y=shift_y,
1039
        )
1040
    elif np_image.ndim == 3:
1041
        return _apply_scalar_per_pixel_3D(
1042
            generic_cy._minimum_3D,
1043
            image,
1044
            footprint,
1045
            out=out,
1046
            mask=mask,
1047
            shift_x=shift_x,
1048
            shift_y=shift_y,
1049
            shift_z=shift_z,
1050
        )
1051
    raise ValueError(f'`image` must have 2 or 3 dimensions, got {np_image.ndim}.')
1052

1053

1054
def modal(image, footprint, out=None, mask=None, shift_x=0, shift_y=0, shift_z=0):
1055
    """Return local mode of an image.
1056

1057
    The mode is the value that appears most often in the local histogram.
1058

1059
    Parameters
1060
    ----------
1061
    image : ([P,] M, N) ndarray (uint8, uint16)
1062
        Input image.
1063
    footprint : ndarray
1064
        The neighborhood expressed as an ndarray of 1's and 0's.
1065
    out : ([P,] M, N) array (same dtype as input)
1066
        If None, a new array is allocated.
1067
    mask : ndarray (integer or float), optional
1068
        Mask array that defines (>0) area of the image included in the local
1069
        neighborhood. If None, the complete image is used (default).
1070
    shift_x, shift_y, shift_z : int
1071
        Offset added to the footprint center point. Shift is bounded to the
1072
        footprint sizes (center must be inside the given footprint).
1073

1074
    Returns
1075
    -------
1076
    out : ([P,] M, N) ndarray (same dtype as input image)
1077
        Output image.
1078

1079
    Examples
1080
    --------
1081
    >>> from skimage import data
1082
    >>> from skimage.morphology import disk, ball
1083
    >>> from skimage.filters.rank import modal
1084
    >>> import numpy as np
1085
    >>> img = data.camera()
1086
    >>> rng = np.random.default_rng()
1087
    >>> volume = rng.integers(0, 255, size=(10,10,10), dtype=np.uint8)
1088
    >>> out = modal(img, disk(5))
1089
    >>> out_vol = modal(volume, ball(5))
1090

1091
    """
1092

1093
    np_image = np.asanyarray(image)
1094
    if np_image.ndim == 2:
1095
        return _apply_scalar_per_pixel(
1096
            generic_cy._modal,
1097
            image,
1098
            footprint,
1099
            out=out,
1100
            mask=mask,
1101
            shift_x=shift_x,
1102
            shift_y=shift_y,
1103
        )
1104
    elif np_image.ndim == 3:
1105
        return _apply_scalar_per_pixel_3D(
1106
            generic_cy._modal_3D,
1107
            image,
1108
            footprint,
1109
            out=out,
1110
            mask=mask,
1111
            shift_x=shift_x,
1112
            shift_y=shift_y,
1113
            shift_z=shift_z,
1114
        )
1115
    raise ValueError(f'`image` must have 2 or 3 dimensions, got {np_image.ndim}.')
1116

1117

1118
def enhance_contrast(
1119
    image, footprint, out=None, mask=None, shift_x=0, shift_y=0, shift_z=0
1120
):
1121
    """Enhance contrast of an image.
1122

1123
    This replaces each pixel by the local maximum if the pixel gray value is
1124
    closer to the local maximum than the local minimum. Otherwise it is
1125
    replaced by the local minimum.
1126

1127
    Parameters
1128
    ----------
1129
    image : ([P,] M, N) ndarray (uint8, uint16)
1130
        Input image.
1131
    footprint : ndarray
1132
        The neighborhood expressed as an ndarray of 1's and 0's.
1133
    out : ([P,] M, N) array (same dtype as input)
1134
        If None, a new array is allocated.
1135
    mask : ndarray (integer or float), optional
1136
        Mask array that defines (>0) area of the image included in the local
1137
        neighborhood. If None, the complete image is used (default).
1138
    shift_x, shift_y, shift_z : int
1139
        Offset added to the footprint center point. Shift is bounded to the
1140
        footprint sizes (center must be inside the given footprint).
1141

1142
    Returns
1143
    -------
1144
    out : ([P,] M, N) ndarray (same dtype as input image)
1145
        Output image
1146

1147
    Examples
1148
    --------
1149
    >>> from skimage import data
1150
    >>> from skimage.morphology import disk, ball
1151
    >>> from skimage.filters.rank import enhance_contrast
1152
    >>> import numpy as np
1153
    >>> img = data.camera()
1154
    >>> rng = np.random.default_rng()
1155
    >>> volume = rng.integers(0, 255, size=(10,10,10), dtype=np.uint8)
1156
    >>> out = enhance_contrast(img, disk(5))
1157
    >>> out_vol = enhance_contrast(volume, ball(5))
1158

1159
    """
1160

1161
    np_image = np.asanyarray(image)
1162
    if np_image.ndim == 2:
1163
        return _apply_scalar_per_pixel(
1164
            generic_cy._enhance_contrast,
1165
            image,
1166
            footprint,
1167
            out=out,
1168
            mask=mask,
1169
            shift_x=shift_x,
1170
            shift_y=shift_y,
1171
        )
1172
    elif np_image.ndim == 3:
1173
        return _apply_scalar_per_pixel_3D(
1174
            generic_cy._enhance_contrast_3D,
1175
            image,
1176
            footprint,
1177
            out=out,
1178
            mask=mask,
1179
            shift_x=shift_x,
1180
            shift_y=shift_y,
1181
            shift_z=shift_z,
1182
        )
1183
    raise ValueError(f'`image` must have 2 or 3 dimensions, got {np_image.ndim}.')
1184

1185

1186
def pop(image, footprint, out=None, mask=None, shift_x=0, shift_y=0, shift_z=0):
1187
    """Return the local number (population) of pixels.
1188

1189
    The number of pixels is defined as the number of pixels which are included
1190
    in the footprint and the mask.
1191

1192
    Parameters
1193
    ----------
1194
    image : ([P,] M, N) ndarray (uint8, uint16)
1195
        Input image.
1196
    footprint : ndarray
1197
        The neighborhood expressed as an ndarray of 1's and 0's.
1198
    out : ([P,] M, N) array (same dtype as input)
1199
        If None, a new array is allocated.
1200
    mask : ndarray (integer or float), optional
1201
        Mask array that defines (>0) area of the image included in the local
1202
        neighborhood. If None, the complete image is used (default).
1203
    shift_x, shift_y, shift_z : int
1204
        Offset added to the footprint center point. Shift is bounded to the
1205
        footprint sizes (center must be inside the given footprint).
1206

1207
    Returns
1208
    -------
1209
    out : ([P,] M, N) ndarray (same dtype as input image)
1210
        Output image.
1211

1212
    Examples
1213
    --------
1214
    >>> from skimage.morphology import square, cube # Need to add 3D example
1215
    >>> import skimage.filters.rank as rank
1216
    >>> img = 255 * np.array([[0, 0, 0, 0, 0],
1217
    ...                       [0, 1, 1, 1, 0],
1218
    ...                       [0, 1, 1, 1, 0],
1219
    ...                       [0, 1, 1, 1, 0],
1220
    ...                       [0, 0, 0, 0, 0]], dtype=np.uint8)
1221
    >>> rank.pop(img, square(3))
1222
    array([[4, 6, 6, 6, 4],
1223
           [6, 9, 9, 9, 6],
1224
           [6, 9, 9, 9, 6],
1225
           [6, 9, 9, 9, 6],
1226
           [4, 6, 6, 6, 4]], dtype=uint8)
1227

1228
    """
1229

1230
    np_image = np.asanyarray(image)
1231
    if np_image.ndim == 2:
1232
        return _apply_scalar_per_pixel(
1233
            generic_cy._pop,
1234
            image,
1235
            footprint,
1236
            out=out,
1237
            mask=mask,
1238
            shift_x=shift_x,
1239
            shift_y=shift_y,
1240
        )
1241
    elif np_image.ndim == 3:
1242
        return _apply_scalar_per_pixel_3D(
1243
            generic_cy._pop_3D,
1244
            image,
1245
            footprint,
1246
            out=out,
1247
            mask=mask,
1248
            shift_x=shift_x,
1249
            shift_y=shift_y,
1250
            shift_z=shift_z,
1251
        )
1252
    raise ValueError(f'`image` must have 2 or 3 dimensions, got {np_image.ndim}.')
1253

1254

1255
def sum(image, footprint, out=None, mask=None, shift_x=0, shift_y=0, shift_z=0):
1256
    """Return the local sum of pixels.
1257

1258
    Note that the sum may overflow depending on the data type of the input
1259
    array.
1260

1261
    Parameters
1262
    ----------
1263
    image : ([P,] M, N) ndarray (uint8, uint16)
1264
        Input image.
1265
    footprint : ndarray
1266
        The neighborhood expressed as an ndarray of 1's and 0's.
1267
    out : ([P,] M, N) array (same dtype as input)
1268
        If None, a new array is allocated.
1269
    mask : ndarray (integer or float), optional
1270
        Mask array that defines (>0) area of the image included in the local
1271
        neighborhood. If None, the complete image is used (default).
1272
    shift_x, shift_y, shift_z : int
1273
        Offset added to the footprint center point. Shift is bounded to the
1274
        footprint sizes (center must be inside the given footprint).
1275

1276
    Returns
1277
    -------
1278
    out : ([P,] M, N) ndarray (same dtype as input image)
1279
        Output image.
1280

1281
    Examples
1282
    --------
1283
    >>> from skimage.morphology import square, cube # Need to add 3D example
1284
    >>> import skimage.filters.rank as rank         # Cube seems to fail but
1285
    >>> img = np.array([[0, 0, 0, 0, 0],            # Ball can pass
1286
    ...                 [0, 1, 1, 1, 0],
1287
    ...                 [0, 1, 1, 1, 0],
1288
    ...                 [0, 1, 1, 1, 0],
1289
    ...                 [0, 0, 0, 0, 0]], dtype=np.uint8)
1290
    >>> rank.sum(img, square(3))
1291
    array([[1, 2, 3, 2, 1],
1292
           [2, 4, 6, 4, 2],
1293
           [3, 6, 9, 6, 3],
1294
           [2, 4, 6, 4, 2],
1295
           [1, 2, 3, 2, 1]], dtype=uint8)
1296

1297
    """
1298

1299
    np_image = np.asanyarray(image)
1300
    if np_image.ndim == 2:
1301
        return _apply_scalar_per_pixel(
1302
            generic_cy._sum,
1303
            image,
1304
            footprint,
1305
            out=out,
1306
            mask=mask,
1307
            shift_x=shift_x,
1308
            shift_y=shift_y,
1309
        )
1310
    elif np_image.ndim == 3:
1311
        return _apply_scalar_per_pixel_3D(
1312
            generic_cy._sum_3D,
1313
            image,
1314
            footprint,
1315
            out=out,
1316
            mask=mask,
1317
            shift_x=shift_x,
1318
            shift_y=shift_y,
1319
            shift_z=shift_z,
1320
        )
1321
    raise ValueError(f'`image` must have 2 or 3 dimensions, got {np_image.ndim}.')
1322

1323

1324
def threshold(image, footprint, out=None, mask=None, shift_x=0, shift_y=0, shift_z=0):
1325
    """Local threshold of an image.
1326

1327
    The resulting binary mask is True if the gray value of the center pixel is
1328
    greater than the local mean.
1329

1330
    Parameters
1331
    ----------
1332
    image : ([P,] M, N) ndarray (uint8, uint16)
1333
        Input image.
1334
    footprint : ndarray
1335
        The neighborhood expressed as an ndarray of 1's and 0's.
1336
    out : ([P,] M, N) array (same dtype as input)
1337
        If None, a new array is allocated.
1338
    mask : ndarray (integer or float), optional
1339
        Mask array that defines (>0) area of the image included in the local
1340
        neighborhood. If None, the complete image is used (default).
1341
    shift_x, shift_y, shift_z : int
1342
        Offset added to the footprint center point. Shift is bounded to the
1343
        footprint sizes (center must be inside the given footprint).
1344

1345
    Returns
1346
    -------
1347
    out : ([P,] M, N) ndarray (same dtype as input image)
1348
        Output image.
1349

1350
    Examples
1351
    --------
1352
    >>> from skimage.morphology import square, cube # Need to add 3D example
1353
    >>> from skimage.filters.rank import threshold
1354
    >>> img = 255 * np.array([[0, 0, 0, 0, 0],
1355
    ...                       [0, 1, 1, 1, 0],
1356
    ...                       [0, 1, 1, 1, 0],
1357
    ...                       [0, 1, 1, 1, 0],
1358
    ...                       [0, 0, 0, 0, 0]], dtype=np.uint8)
1359
    >>> threshold(img, square(3))
1360
    array([[0, 0, 0, 0, 0],
1361
           [0, 1, 1, 1, 0],
1362
           [0, 1, 0, 1, 0],
1363
           [0, 1, 1, 1, 0],
1364
           [0, 0, 0, 0, 0]], dtype=uint8)
1365

1366
    """
1367

1368
    np_image = np.asanyarray(image)
1369
    if np_image.ndim == 2:
1370
        return _apply_scalar_per_pixel(
1371
            generic_cy._threshold,
1372
            image,
1373
            footprint,
1374
            out=out,
1375
            mask=mask,
1376
            shift_x=shift_x,
1377
            shift_y=shift_y,
1378
        )
1379
    elif np_image.ndim == 3:
1380
        return _apply_scalar_per_pixel_3D(
1381
            generic_cy._threshold_3D,
1382
            image,
1383
            footprint,
1384
            out=out,
1385
            mask=mask,
1386
            shift_x=shift_x,
1387
            shift_y=shift_y,
1388
            shift_z=shift_z,
1389
        )
1390
    raise ValueError(f'`image` must have 2 or 3 dimensions, got {np_image.ndim}.')
1391

1392

1393
def noise_filter(
1394
    image, footprint, out=None, mask=None, shift_x=0, shift_y=0, shift_z=0
1395
):
1396
    """Noise feature.
1397

1398
    Parameters
1399
    ----------
1400
    image : ([P,] M, N) ndarray (uint8, uint16)
1401
        Input image.
1402
    footprint : ndarray
1403
        The neighborhood expressed as an ndarray of 1's and 0's.
1404
    out : ([P,] M, N) array (same dtype as input)
1405
        If None, a new array is allocated.
1406
    mask : ndarray (integer or float), optional
1407
        Mask array that defines (>0) area of the image included in the local
1408
        neighborhood. If None, the complete image is used (default).
1409
    shift_x, shift_y, shift_z : int
1410
        Offset added to the footprint center point. Shift is bounded to the
1411
        footprint sizes (center must be inside the given footprint).
1412

1413
    References
1414
    ----------
1415
    .. [1] N. Hashimoto et al. Referenceless image quality evaluation
1416
                     for whole slide imaging. J Pathol Inform 2012;3:9.
1417

1418
    Returns
1419
    -------
1420
    out : ([P,] M, N) ndarray (same dtype as input image)
1421
        Output image.
1422

1423
    Examples
1424
    --------
1425
    >>> from skimage import data
1426
    >>> from skimage.morphology import disk, ball
1427
    >>> from skimage.filters.rank import noise_filter
1428
    >>> import numpy as np
1429
    >>> img = data.camera()
1430
    >>> rng = np.random.default_rng()
1431
    >>> volume = rng.integers(0, 255, size=(10,10,10), dtype=np.uint8)
1432
    >>> out = noise_filter(img, disk(5))
1433
    >>> out_vol = noise_filter(volume, ball(5))
1434

1435
    """
1436

1437
    np_image = np.asanyarray(image)
1438
    if _footprint_is_sequence(footprint):
1439
        raise ValueError(
1440
            "footprint sequences are not currently supported by rank filters"
1441
        )
1442
    if np_image.ndim == 2:
1443
        # ensure that the central pixel in the footprint is empty
1444
        centre_r = int(footprint.shape[0] / 2) + shift_y
1445
        centre_c = int(footprint.shape[1] / 2) + shift_x
1446
        # make a local copy
1447
        footprint_cpy = footprint.copy()
1448
        footprint_cpy[centre_r, centre_c] = 0
1449

1450
        return _apply_scalar_per_pixel(
1451
            generic_cy._noise_filter,
1452
            image,
1453
            footprint_cpy,
1454
            out=out,
1455
            mask=mask,
1456
            shift_x=shift_x,
1457
            shift_y=shift_y,
1458
        )
1459
    elif np_image.ndim == 3:
1460
        # ensure that the central pixel in the footprint is empty
1461
        centre_r = int(footprint.shape[0] / 2) + shift_y
1462
        centre_c = int(footprint.shape[1] / 2) + shift_x
1463
        centre_z = int(footprint.shape[2] / 2) + shift_z
1464
        # make a local copy
1465
        footprint_cpy = footprint.copy()
1466
        footprint_cpy[centre_r, centre_c, centre_z] = 0
1467

1468
        return _apply_scalar_per_pixel_3D(
1469
            generic_cy._noise_filter_3D,
1470
            image,
1471
            footprint_cpy,
1472
            out=out,
1473
            mask=mask,
1474
            shift_x=shift_x,
1475
            shift_y=shift_y,
1476
            shift_z=shift_z,
1477
        )
1478

1479
    raise ValueError(f'`image` must have 2 or 3 dimensions, got {np_image.ndim}.')
1480

1481

1482
def entropy(image, footprint, out=None, mask=None, shift_x=0, shift_y=0, shift_z=0):
1483
    """Local entropy.
1484

1485
    The entropy is computed using base 2 logarithm i.e. the filter returns the
1486
    minimum number of bits needed to encode the local gray level
1487
    distribution.
1488

1489
    Parameters
1490
    ----------
1491
    image : ([P,] M, N) ndarray (uint8, uint16)
1492
        Input image.
1493
    footprint : ndarray
1494
        The neighborhood expressed as an ndarray of 1's and 0's.
1495
    out : ([P,] M, N) array (same dtype as input)
1496
        If None, a new array is allocated.
1497
    mask : ndarray (integer or float), optional
1498
        Mask array that defines (>0) area of the image included in the local
1499
        neighborhood. If None, the complete image is used (default).
1500
    shift_x, shift_y, shift_z : int
1501
        Offset added to the footprint center point. Shift is bounded to the
1502
        footprint sizes (center must be inside the given footprint).
1503

1504
    Returns
1505
    -------
1506
    out : ([P,] M, N) ndarray (float)
1507
        Output image.
1508

1509
    References
1510
    ----------
1511
    .. [1] `https://en.wikipedia.org/wiki/Entropy_(information_theory) <https://en.wikipedia.org/wiki/Entropy_(information_theory)>`_
1512

1513
    Examples
1514
    --------
1515
    >>> from skimage import data
1516
    >>> from skimage.filters.rank import entropy
1517
    >>> from skimage.morphology import disk, ball
1518
    >>> import numpy as np
1519
    >>> img = data.camera()
1520
    >>> rng = np.random.default_rng()
1521
    >>> volume = rng.integers(0, 255, size=(10,10,10), dtype=np.uint8)
1522
    >>> ent = entropy(img, disk(5))
1523
    >>> ent_vol = entropy(volume, ball(5))
1524

1525
    """
1526

1527
    np_image = np.asanyarray(image)
1528
    if np_image.ndim == 2:
1529
        return _apply_scalar_per_pixel(
1530
            generic_cy._entropy,
1531
            image,
1532
            footprint,
1533
            out=out,
1534
            mask=mask,
1535
            shift_x=shift_x,
1536
            shift_y=shift_y,
1537
            out_dtype=np.float64,
1538
        )
1539
    elif np_image.ndim == 3:
1540
        return _apply_scalar_per_pixel_3D(
1541
            generic_cy._entropy_3D,
1542
            image,
1543
            footprint,
1544
            out=out,
1545
            mask=mask,
1546
            shift_x=shift_x,
1547
            shift_y=shift_y,
1548
            shift_z=shift_z,
1549
            out_dtype=np.float64,
1550
        )
1551
    raise ValueError(f'`image` must have 2 or 3 dimensions, got {np_image.ndim}.')
1552

1553

1554
def otsu(image, footprint, out=None, mask=None, shift_x=0, shift_y=0, shift_z=0):
1555
    """Local Otsu's threshold value for each pixel.
1556

1557
    Parameters
1558
    ----------
1559
    image : ([P,] M, N) ndarray (uint8, uint16)
1560
        Input image.
1561
    footprint : ndarray
1562
        The neighborhood expressed as an ndarray of 1's and 0's.
1563
    out : ([P,] M, N) array (same dtype as input)
1564
        If None, a new array is allocated.
1565
    mask : ndarray (integer or float), optional
1566
        Mask array that defines (>0) area of the image included in the local
1567
        neighborhood. If None, the complete image is used (default).
1568
    shift_x, shift_y, shift_z : int
1569
        Offset added to the footprint center point. Shift is bounded to the
1570
        footprint sizes (center must be inside the given footprint).
1571

1572
    Returns
1573
    -------
1574
    out : ([P,] M, N) ndarray (same dtype as input image)
1575
        Output image.
1576

1577
    References
1578
    ----------
1579
    .. [1] https://en.wikipedia.org/wiki/Otsu's_method
1580

1581
    Examples
1582
    --------
1583
    >>> from skimage import data
1584
    >>> from skimage.filters.rank import otsu
1585
    >>> from skimage.morphology import disk, ball
1586
    >>> import numpy as np
1587
    >>> img = data.camera()
1588
    >>> rng = np.random.default_rng()
1589
    >>> volume = rng.integers(0, 255, size=(10,10,10), dtype=np.uint8)
1590
    >>> local_otsu = otsu(img, disk(5))
1591
    >>> thresh_image = img >= local_otsu
1592
    >>> local_otsu_vol = otsu(volume, ball(5))
1593
    >>> thresh_image_vol = volume >= local_otsu_vol
1594

1595
    """
1596

1597
    np_image = np.asanyarray(image)
1598
    if np_image.ndim == 2:
1599
        return _apply_scalar_per_pixel(
1600
            generic_cy._otsu,
1601
            image,
1602
            footprint,
1603
            out=out,
1604
            mask=mask,
1605
            shift_x=shift_x,
1606
            shift_y=shift_y,
1607
        )
1608
    elif np_image.ndim == 3:
1609
        return _apply_scalar_per_pixel_3D(
1610
            generic_cy._otsu_3D,
1611
            image,
1612
            footprint,
1613
            out=out,
1614
            mask=mask,
1615
            shift_x=shift_x,
1616
            shift_y=shift_y,
1617
            shift_z=shift_z,
1618
        )
1619
    raise ValueError(f'`image` must have 2 or 3 dimensions, got {np_image.ndim}.')
1620

1621

1622
def windowed_histogram(
1623
    image, footprint, out=None, mask=None, shift_x=0, shift_y=0, n_bins=None
1624
):
1625
    """Normalized sliding window histogram
1626

1627
    Parameters
1628
    ----------
1629
    image : 2-D array (integer or float)
1630
        Input image.
1631
    footprint : 2-D array (integer or float)
1632
        The neighborhood expressed as a 2-D array of 1's and 0's.
1633
    out : 2-D array (integer or float), optional
1634
        If None, a new array is allocated.
1635
    mask : ndarray (integer or float), optional
1636
        Mask array that defines (>0) area of the image included in the local
1637
        neighborhood. If None, the complete image is used (default).
1638
    shift_x, shift_y : int, optional
1639
        Offset added to the footprint center point. Shift is bounded to the
1640
        footprint sizes (center must be inside the given footprint).
1641
    n_bins : int or None
1642
        The number of histogram bins. Will default to ``image.max() + 1``
1643
        if None is passed.
1644

1645
    Returns
1646
    -------
1647
    out : 3-D array (float)
1648
        Array of dimensions (H,W,N), where (H,W) are the dimensions of the
1649
        input image and N is n_bins or ``image.max() + 1`` if no value is
1650
        provided as a parameter. Effectively, each pixel is a N-D feature
1651
        vector that is the histogram. The sum of the elements in the feature
1652
        vector will be 1, unless no pixels in the window were covered by both
1653
        footprint and mask, in which case all elements will be 0.
1654

1655
    Examples
1656
    --------
1657
    >>> from skimage import data
1658
    >>> from skimage.filters.rank import windowed_histogram
1659
    >>> from skimage.morphology import disk, ball
1660
    >>> import numpy as np
1661
    >>> img = data.camera()
1662
    >>> rng = np.random.default_rng()
1663
    >>> volume = rng.integers(0, 255, size=(10,10,10), dtype=np.uint8)
1664
    >>> hist_img = windowed_histogram(img, disk(5))
1665

1666
    """
1667

1668
    if n_bins is None:
1669
        n_bins = int(image.max()) + 1
1670

1671
    return _apply_vector_per_pixel(
1672
        generic_cy._windowed_hist,
1673
        image,
1674
        footprint,
1675
        out=out,
1676
        mask=mask,
1677
        shift_x=shift_x,
1678
        shift_y=shift_y,
1679
        out_dtype=np.float64,
1680
        pixel_size=n_bins,
1681
    )
1682

1683

1684
def majority(
1685
    image,
1686
    footprint,
1687
    *,
1688
    out=None,
1689
    mask=None,
1690
    shift_x=0,
1691
    shift_y=0,
1692
    shift_z=0,
1693
):
1694
    """Assign to each pixel the most common value within its neighborhood.
1695

1696
    Parameters
1697
    ----------
1698
    image : ndarray
1699
        Image array (uint8, uint16 array).
1700
    footprint : 2-D array (integer or float)
1701
        The neighborhood expressed as a 2-D array of 1's and 0's.
1702
    out : ndarray (integer or float), optional
1703
        If None, a new array will be allocated.
1704
    mask : ndarray (integer or float), optional
1705
        Mask array that defines (>0) area of the image included in the local
1706
        neighborhood. If None, the complete image is used (default).
1707
    shift_x, shift_y : int, optional
1708
        Offset added to the footprint center point. Shift is bounded to the
1709
        footprint sizes (center must be inside the given footprint).
1710

1711
    Returns
1712
    -------
1713
    out : 2-D array (same dtype as input image)
1714
        Output image.
1715

1716
    Examples
1717
    --------
1718
    >>> from skimage import data
1719
    >>> from skimage.filters.rank import majority
1720
    >>> from skimage.morphology import disk, ball
1721
    >>> import numpy as np
1722
    >>> img = data.camera()
1723
    >>> rng = np.random.default_rng()
1724
    >>> volume = rng.integers(0, 255, size=(10,10,10), dtype=np.uint8)
1725
    >>> maj_img = majority(img, disk(5))
1726
    >>> maj_img_vol = majority(volume, ball(5))
1727

1728
    """
1729

1730
    np_image = np.asanyarray(image)
1731
    if np_image.ndim == 2:
1732
        return _apply_scalar_per_pixel(
1733
            generic_cy._majority,
1734
            image,
1735
            footprint,
1736
            out=out,
1737
            mask=mask,
1738
            shift_x=shift_x,
1739
            shift_y=shift_y,
1740
        )
1741
    elif np_image.ndim == 3:
1742
        return _apply_scalar_per_pixel_3D(
1743
            generic_cy._majority_3D,
1744
            image,
1745
            footprint,
1746
            out=out,
1747
            mask=mask,
1748
            shift_x=shift_x,
1749
            shift_y=shift_y,
1750
            shift_z=shift_z,
1751
        )
1752
    raise ValueError(f'`image` must have 2 or 3 dimensions, got {np_image.ndim}.')
1753