scikit-image

1
"""Miscellaneous morphology functions."""
2

3
import numpy as np
4
import functools
5
from scipy import ndimage as ndi
6
from scipy.spatial import cKDTree
7

8
from .._shared.utils import warn
9
from ._misc_cy import _remove_objects_by_distance
10

11

12
# Our function names don't exactly correspond to ndimages.
13
# This dictionary translates from our names to scipy's.
14
funcs = ('erosion', 'dilation', 'opening', 'closing')
15
skimage2ndimage = {x: 'grey_' + x for x in funcs}
16

17
# These function names are the same in ndimage.
18
funcs = (
19
    'binary_erosion',
20
    'binary_dilation',
21
    'binary_opening',
22
    'binary_closing',
23
    'black_tophat',
24
    'white_tophat',
25
)
26
skimage2ndimage.update({x: x for x in funcs})
27

28

29
def default_footprint(func):
30
    """Decorator to add a default footprint to morphology functions.
31

32
    Parameters
33
    ----------
34
    func : function
35
        A morphology function such as erosion, dilation, opening, closing,
36
        white_tophat, or black_tophat.
37

38
    Returns
39
    -------
40
    func_out : function
41
        The function, using a default footprint of same dimension
42
        as the input image with connectivity 1.
43

44
    """
45

46
    @functools.wraps(func)
47
    def func_out(image, footprint=None, *args, **kwargs):
48
        if footprint is None:
49
            footprint = ndi.generate_binary_structure(image.ndim, 1)
50
        return func(image, footprint=footprint, *args, **kwargs)
51

52
    return func_out
53

54

55
def _check_dtype_supported(ar):
56
    # Should use `issubdtype` for bool below, but there's a bug in numpy 1.7
57
    if not (ar.dtype == bool or np.issubdtype(ar.dtype, np.integer)):
58
        raise TypeError(
59
            "Only bool or integer image types are supported. " f"Got {ar.dtype}."
60
        )
61

62

63
def remove_small_objects(ar, min_size=64, connectivity=1, *, out=None):
64
    """Remove objects smaller than the specified size.
65

66
    Expects ar to be an array with labeled objects, and removes objects
67
    smaller than min_size. If `ar` is bool, the image is first labeled.
68
    This leads to potentially different behavior for bool and 0-and-1
69
    arrays.
70

71
    Parameters
72
    ----------
73
    ar : ndarray (arbitrary shape, int or bool type)
74
        The array containing the objects of interest. If the array type is
75
        int, the ints must be non-negative.
76
    min_size : int, optional (default: 64)
77
        The smallest allowable object size.
78
    connectivity : int, {1, 2, ..., ar.ndim}, optional (default: 1)
79
        The connectivity defining the neighborhood of a pixel. Used during
80
        labelling if `ar` is bool.
81
    out : ndarray
82
        Array of the same shape as `ar`, into which the output is
83
        placed. By default, a new array is created.
84

85
    Raises
86
    ------
87
    TypeError
88
        If the input array is of an invalid type, such as float or string.
89
    ValueError
90
        If the input array contains negative values.
91

92
    Returns
93
    -------
94
    out : ndarray, same shape and type as input `ar`
95
        The input array with small connected components removed.
96

97
    See Also
98
    --------
99
    skimage.morphology.remove_objects_by_distance
100

101
    Examples
102
    --------
103
    >>> from skimage import morphology
104
    >>> a = np.array([[0, 0, 0, 1, 0],
105
    ...               [1, 1, 1, 0, 0],
106
    ...               [1, 1, 1, 0, 1]], bool)
107
    >>> b = morphology.remove_small_objects(a, 6)
108
    >>> b
109
    array([[False, False, False, False, False],
110
           [ True,  True,  True, False, False],
111
           [ True,  True,  True, False, False]])
112
    >>> c = morphology.remove_small_objects(a, 7, connectivity=2)
113
    >>> c
114
    array([[False, False, False,  True, False],
115
           [ True,  True,  True, False, False],
116
           [ True,  True,  True, False, False]])
117
    >>> d = morphology.remove_small_objects(a, 6, out=a)
118
    >>> d is a
119
    True
120

121
    """
122
    # Raising type error if not int or bool
123
    _check_dtype_supported(ar)
124

125
    if out is None:
126
        out = ar.copy()
127
    else:
128
        out[:] = ar
129

130
    if min_size == 0:  # shortcut for efficiency
131
        return out
132

133
    if out.dtype == bool:
134
        footprint = ndi.generate_binary_structure(ar.ndim, connectivity)
135
        ccs = np.zeros_like(ar, dtype=np.int32)
136
        ndi.label(ar, footprint, output=ccs)
137
    else:
138
        ccs = out
139

140
    try:
141
        component_sizes = np.bincount(ccs.ravel())
142
    except ValueError:
143
        raise ValueError(
144
            "Negative value labels are not supported. Try "
145
            "relabeling the input with `scipy.ndimage.label` or "
146
            "`skimage.morphology.label`."
147
        )
148

149
    if len(component_sizes) == 2 and out.dtype != bool:
150
        warn(
151
            "Only one label was provided to `remove_small_objects`. "
152
            "Did you mean to use a boolean array?"
153
        )
154

155
    too_small = component_sizes < min_size
156
    too_small_mask = too_small[ccs]
157
    out[too_small_mask] = 0
158

159
    return out
160

161

162
def remove_small_holes(ar, area_threshold=64, connectivity=1, *, out=None):
163
    """Remove contiguous holes smaller than the specified size.
164

165
    Parameters
166
    ----------
167
    ar : ndarray (arbitrary shape, int or bool type)
168
        The array containing the connected components of interest.
169
    area_threshold : int, optional (default: 64)
170
        The maximum area, in pixels, of a contiguous hole that will be filled.
171
        Replaces `min_size`.
172
    connectivity : int, {1, 2, ..., ar.ndim}, optional (default: 1)
173
        The connectivity defining the neighborhood of a pixel.
174
    out : ndarray
175
        Array of the same shape as `ar` and bool dtype, into which the
176
        output is placed. By default, a new array is created.
177

178
    Raises
179
    ------
180
    TypeError
181
        If the input array is of an invalid type, such as float or string.
182
    ValueError
183
        If the input array contains negative values.
184

185
    Returns
186
    -------
187
    out : ndarray, same shape and type as input `ar`
188
        The input array with small holes within connected components removed.
189

190
    Examples
191
    --------
192
    >>> from skimage import morphology
193
    >>> a = np.array([[1, 1, 1, 1, 1, 0],
194
    ...               [1, 1, 1, 0, 1, 0],
195
    ...               [1, 0, 0, 1, 1, 0],
196
    ...               [1, 1, 1, 1, 1, 0]], bool)
197
    >>> b = morphology.remove_small_holes(a, 2)
198
    >>> b
199
    array([[ True,  True,  True,  True,  True, False],
200
           [ True,  True,  True,  True,  True, False],
201
           [ True, False, False,  True,  True, False],
202
           [ True,  True,  True,  True,  True, False]])
203
    >>> c = morphology.remove_small_holes(a, 2, connectivity=2)
204
    >>> c
205
    array([[ True,  True,  True,  True,  True, False],
206
           [ True,  True,  True, False,  True, False],
207
           [ True, False, False,  True,  True, False],
208
           [ True,  True,  True,  True,  True, False]])
209
    >>> d = morphology.remove_small_holes(a, 2, out=a)
210
    >>> d is a
211
    True
212

213
    Notes
214
    -----
215
    If the array type is int, it is assumed that it contains already-labeled
216
    objects. The labels are not kept in the output image (this function always
217
    outputs a bool image). It is suggested that labeling is completed after
218
    using this function.
219

220
    """
221
    _check_dtype_supported(ar)
222

223
    # Creates warning if image is an integer image
224
    if ar.dtype != bool:
225
        warn(
226
            "Any labeled images will be returned as a boolean array. "
227
            "Did you mean to use a boolean array?",
228
            UserWarning,
229
        )
230

231
    if out is not None:
232
        if out.dtype != bool:
233
            raise TypeError("out dtype must be bool")
234
    else:
235
        out = ar.astype(bool, copy=True)
236

237
    # Creating the inverse of ar
238
    np.logical_not(ar, out=out)
239

240
    # removing small objects from the inverse of ar
241
    out = remove_small_objects(out, area_threshold, connectivity, out=out)
242

243
    np.logical_not(out, out=out)
244

245
    return out
246

247

248
def remove_objects_by_distance(
249
    label_image,
250
    min_distance,
251
    *,
252
    priority=None,
253
    p_norm=2,
254
    spacing=None,
255
    out=None,
256
):
257
    """Remove objects, in specified order, until remaining are a minimum distance apart.
258

259
    Remove labeled objects from an image until the remaining ones are spaced
260
    more than a given distance from one another. By default, smaller objects
261
    are removed first.
262

263
    Parameters
264
    ----------
265
    label_image : ndarray of integers
266
        An n-dimensional array containing object labels, e.g. as returned by
267
        :func:`~.label`. A value of zero is considered background, all other
268
        object IDs must be positive integers.
269
    min_distance : int or float
270
        Remove objects whose distance to other objects is not greater than this
271
        positive value. Objects with a lower `priority` are removed first.
272
    priority : ndarray, optional
273
        Defines the priority with which objects are removed. Expects a
274
        1-dimensional array of length
275
        :func:`np.amax(label_image) + 1 <numpy.amax>` that contains the priority
276
        for each object's label at the respective index. Objects with a lower value
277
        are removed first until all remaining objects fulfill the distance
278
        requirement. If not given, priority is given to objects with a higher
279
        number of samples and their label value second.
280
    p_norm : int or float, optional
281
        The Minkowski distance of order p, used to calculate the distance
282
        between objects. The default ``2`` corresponds to the Euclidean
283
        distance, ``1`` to the "Manhattan" distance, and ``np.inf`` to the
284
        Chebyshev distance.
285
    spacing : sequence of float, optional
286
        The pixel spacing along each axis of `label_image`. If not specified,
287
        a grid spacing of unity (1) is implied.
288
    out : ndarray, optional
289
        Array of the same shape and dtype as `image`, into which the output is
290
        placed. By default, a new array is created.
291

292
    Returns
293
    -------
294
    out : ndarray
295
        Array of the same shape as `label_image`, for which objects that violate
296
        the `min_distance` condition were removed.
297

298
    See Also
299
    --------
300
    skimage.morphology.remove_small_objects
301
        Remove objects smaller than the specified size.
302

303
    Notes
304
    -----
305
    The basic steps of this algorithm work as follows:
306

307
    1. Find the indices for of all given objects and separate them depending on
308
       if they point to an object's border or not.
309
    2. Sort indices by their label value, ensuring that indices which point to
310
       the same object are next to each other. This optimization allows finding
311
       all parts of an object, simply by stepping to the neighboring indices.
312
    3. Sort boundary indices by `priority`. Use a stable-sort to preserve the
313
       ordering from the previous sorting step. If `priority` is not given,
314
       use :func:`numpy.bincount` as a fallback.
315
    4. Construct a :class:`scipy.spatial.cKDTree` from the boundary indices.
316
    5. Iterate across boundary indices in priority-sorted order, and query the
317
       kd-tree for objects that are too close. Remove ones that are and don't
318
       take them into account when evaluating other objects later on.
319

320
    The performance of this algorithm depends on the number of samples in
321
    `label_image` that belong to an object's border.
322

323
    Examples
324
    --------
325
    >>> import skimage as ski
326
    >>> ski.morphology.remove_objects_by_distance(np.array([2, 0, 1, 1]), 2)
327
    array([0, 0, 1, 1])
328
    >>> ski.morphology.remove_objects_by_distance(
329
    ...     np.array([2, 0, 1, 1]), 2, priority=np.array([0, 1, 9])
330
    ... )
331
    array([2, 0, 0, 0])
332
    >>> label_image = np.array(
333
    ...     [[8, 0, 0, 0, 0, 0, 0, 0, 0, 9, 9],
334
    ...      [8, 8, 8, 0, 0, 0, 0, 0, 0, 9, 9],
335
    ...      [0, 0, 0, 0, 0, 0, 0, 0, 9, 0, 0],
336
    ...      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
337
    ...      [0, 0, 3, 0, 0, 0, 1, 0, 0, 0, 0],
338
    ...      [2, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
339
    ...      [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
340
    ...      [0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 7]]
341
    ... )
342
    >>> ski.morphology.remove_objects_by_distance(
343
    ...     label_image, min_distance=3
344
    ... )
345
    array([[8, 0, 0, 0, 0, 0, 0, 0, 0, 9, 9],
346
           [8, 8, 8, 0, 0, 0, 0, 0, 0, 9, 9],
347
           [0, 0, 0, 0, 0, 0, 0, 0, 9, 0, 0],
348
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
349
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
350
           [2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
351
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
352
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 7]])
353
    """
354
    if min_distance < 0:
355
        raise ValueError(f"min_distance must be >= 0, was {min_distance}")
356
    if not np.issubdtype(label_image.dtype, np.integer):
357
        raise ValueError(
358
            f"`label_image` must be of integer dtype, got {label_image.dtype}"
359
        )
360
    if out is None:
361
        out = label_image.copy(order="C")
362
    elif out is not label_image:
363
        out[:] = label_image
364
    # May create a copy if order is not C, account for that later
365
    out_raveled = out.ravel(order="C")
366

367
    if spacing is not None:
368
        spacing = np.array(spacing)
369
        if spacing.shape != (out.ndim,) or spacing.min() <= 0:
370
            raise ValueError(
371
                "`spacing` must contain exactly one positive factor "
372
                "for each dimension of `label_image`"
373
            )
374

375
    indices = np.flatnonzero(out_raveled)
376
    # Optimization: Split indices into those on the object boundaries and inner
377
    # ones. The KDTree is built only from the boundary indices, which reduces
378
    # the size of the critical loop significantly! Remaining indices are only
379
    # used to remove the inner parts of objects as well.
380
    if (spacing is None or np.all(spacing[0] == spacing)) and p_norm <= 2:
381
        # For unity spacing we can make the borders more sparse by using a
382
        # lower connectivity
383
        footprint = ndi.generate_binary_structure(out.ndim, 1)
384
    else:
385
        footprint = ndi.generate_binary_structure(out.ndim, out.ndim)
386
    border = (
387
        ndi.maximum_filter(out, footprint=footprint)
388
        != ndi.minimum_filter(out, footprint=footprint)
389
    ).ravel()[indices]
390
    border_indices = indices[border]
391
    inner_indices = indices[~border]
392

393
    if border_indices.size == 0:
394
        # Image without any or only one object, return early
395
        return out
396

397
    # Sort by label ID first, so that IDs of the same object are contiguous
398
    # in the sorted index. This allows fast discovery of the whole object by
399
    # simple iteration up or down the index!
400
    border_indices = border_indices[np.argsort(out_raveled[border_indices])]
401
    inner_indices = inner_indices[np.argsort(out_raveled[inner_indices])]
402

403
    if priority is None:
404
        if not np.can_cast(out.dtype, np.intp, casting="safe"):
405
            # bincount expects intp (32-bit) on WASM or i386, so down-cast to that
406
            priority = np.bincount(out_raveled.astype(np.intp, copy=False))
407
        else:
408
            priority = np.bincount(out_raveled)
409
    # `priority` can only be indexed by positive object IDs,
410
    # `border_indices` contains all unique sorted IDs so check the lowest / first
411
    smallest_id = out_raveled[border_indices[0]]
412
    if smallest_id < 0:
413
        raise ValueError(f"found object with negative ID {smallest_id!r}")
414

415
    try:
416
        # Sort by priority second using a stable sort to preserve the contiguous
417
        # sorting of objects. Because each pixel in an object has the same
418
        # priority we don't need to worry about separating objects.
419
        border_indices = border_indices[
420
            np.argsort(priority[out_raveled[border_indices]], kind="stable")[::-1]
421
        ]
422
    except IndexError as error:
423
        # Use np.amax only for the exception path to provide a nicer error message
424
        expected_shape = (np.amax(out_raveled) + 1,)
425
        if priority.shape != expected_shape:
426
            raise ValueError(
427
                "shape of `priority` must be (np.amax(label_image) + 1,), "
428
                f"expected {expected_shape}, got {priority.shape} instead"
429
            ) from error
430
        else:
431
            raise
432

433
    # Construct kd-tree from unraveled border indices (optionally scale by `spacing`)
434
    unraveled_indices = np.unravel_index(border_indices, out.shape)
435
    if spacing is not None:
436
        unraveled_indices = tuple(
437
            unraveled_indices[dim] * spacing[dim] for dim in range(out.ndim)
438
        )
439
    kdtree = cKDTree(data=np.asarray(unraveled_indices, dtype=np.float64).T)
440

441
    _remove_objects_by_distance(
442
        out=out_raveled,
443
        border_indices=border_indices,
444
        inner_indices=inner_indices,
445
        kdtree=kdtree,
446
        min_distance=min_distance,
447
        p_norm=p_norm,
448
        shape=label_image.shape,
449
    )
450

451
    if out_raveled.base is not out:
452
        # `out_raveled` is a copy, re-assign
453
        out[:] = out_raveled.reshape(out.shape)
454
    return out
455