| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509 |
- """Miscellaneous morphology functions."""
- import numpy as np
- import functools
- import warnings
- from scipy import ndimage as ndi
- from scipy.spatial import cKDTree
- from .._shared.utils import warn, deprecate_parameter, DEPRECATED
- from ._misc_cy import _remove_objects_by_distance
- # Our function names don't exactly correspond to ndimages.
- # This dictionary translates from our names to scipy's.
- funcs = ('erosion', 'dilation', 'opening', 'closing')
- skimage2ndimage = {x: 'grey_' + x for x in funcs}
- # These function names are the same in ndimage.
- funcs = (
- 'black_tophat',
- 'white_tophat',
- )
- skimage2ndimage.update({x: x for x in funcs})
- def default_footprint(func):
- """Decorator to add a default footprint to morphology functions.
- Parameters
- ----------
- func : function
- A morphology function such as erosion, dilation, opening, closing,
- white_tophat, or black_tophat.
- Returns
- -------
- func_out : function
- The function, using a default footprint of same dimension
- as the input image with connectivity 1.
- """
- @functools.wraps(func)
- def func_out(image, footprint=None, *args, **kwargs):
- if footprint is None:
- footprint = ndi.generate_binary_structure(image.ndim, 1)
- return func(image, footprint=footprint, *args, **kwargs)
- return func_out
- def _check_dtype_supported(ar):
- # Should use `issubdtype` for bool below, but there's a bug in numpy 1.7
- if not (ar.dtype == bool or np.issubdtype(ar.dtype, np.integer)):
- raise TypeError(
- "Only bool or integer image types are supported. " f"Got {ar.dtype}."
- )
- @deprecate_parameter(
- deprecated_name="min_size",
- new_name="max_size",
- start_version="0.26.0",
- stop_version="2.0.0",
- template=f"{deprecate_parameter.replace_parameter_template} "
- "Note that the new threshold removes objects smaller than **or equal to** "
- "its value, while the previous parameter only removed smaller ones.",
- )
- def remove_small_objects(
- ar, min_size=DEPRECATED, connectivity=1, *, max_size=64, out=None
- ):
- """Remove objects smaller than the specified size.
- Expects `ar` to be an array with labeled objects, and removes objects
- smaller than or equal to `max_size`. If `ar` is bool, the image is first
- labeled. This leads to potentially different behavior for bool vs. 0-and-1
- arrays.
- Parameters
- ----------
- ar : ndarray (arbitrary shape, int or bool type)
- The array containing the objects of interest. If the array type is
- int, the ints must be non-negative.
- max_size : int, optional (default: 64)
- Remove objects whose contiguous area (or volume, in N-D) contains this
- number of pixels or fewer.
- .. versionadded:: 0.26
- To make the naming clearer, replaces deprecated `min_size`
- which only removed objects strictly smaller than its size.
- connectivity : int, {1, 2, ..., ar.ndim}, optional (default: 1)
- The connectivity defining the neighborhood of a pixel. Used during
- labelling if `ar` is bool.
- out : ndarray
- Array of the same shape as `ar`, into which the output is
- placed. By default, a new array is created.
- Raises
- ------
- TypeError
- If the input array is of an invalid type, such as float or string.
- ValueError
- If the input array contains negative values.
- Returns
- -------
- out : ndarray, same shape and type as input `ar`
- The input array with small connected components removed.
- See Also
- --------
- skimage.morphology.remove_small_holes
- skimage.morphology.remove_objects_by_distance
- Examples
- --------
- >>> from skimage import morphology
- >>> a = np.array([[0, 0, 0, 1, 0],
- ... [1, 1, 1, 0, 0],
- ... [1, 1, 1, 0, 1]], bool)
- >>> b = morphology.remove_small_objects(a, max_size=5)
- >>> b
- array([[False, False, False, False, False],
- [ True, True, True, False, False],
- [ True, True, True, False, False]])
- >>> c = morphology.remove_small_objects(a, max_size=6, connectivity=2)
- >>> c
- array([[False, False, False, True, False],
- [ True, True, True, False, False],
- [ True, True, True, False, False]])
- >>> d = morphology.remove_small_objects(a, max_size=5, out=a)
- >>> d is a
- True
- """
- # Raising type error if not int or bool
- _check_dtype_supported(ar)
- if out is None:
- out = ar.copy()
- else:
- out[:] = ar
- if max_size == 0: # shortcut for efficiency
- return out
- if out.dtype == bool:
- footprint = ndi.generate_binary_structure(ar.ndim, connectivity)
- ccs = np.zeros_like(ar, dtype=np.int32)
- ndi.label(ar, footprint, output=ccs)
- else:
- ccs = out
- try:
- component_sizes = np.bincount(ccs.ravel())
- except ValueError:
- raise ValueError(
- "Negative value labels are not supported. Try "
- "relabeling the input with `scipy.ndimage.label` or "
- "`skimage.morphology.label`."
- )
- if len(component_sizes) == 2 and out.dtype != bool:
- warn(
- "Only one label was provided to `remove_small_objects`. "
- "Did you mean to use a boolean array?"
- )
- if min_size is not DEPRECATED:
- # Exclusive threshold is deprecated behavior
- too_small = component_sizes < min_size
- else:
- # New behavior uses inclusive threshold
- too_small = component_sizes <= max_size
- too_small_mask = too_small[ccs]
- out[too_small_mask] = 0
- return out
- @deprecate_parameter(
- deprecated_name="area_threshold",
- new_name="max_size",
- start_version="0.26.0",
- stop_version="2.0.0",
- template=f"{deprecate_parameter.replace_parameter_template} "
- "Note that the new threshold removes objects smaller than **or equal to** "
- "its value, while the previous parameter only removed smaller ones.",
- )
- def remove_small_holes(
- ar, area_threshold=DEPRECATED, connectivity=1, *, max_size=64, out=None
- ):
- """Remove contiguous holes smaller than the specified size.
- Parameters
- ----------
- ar : ndarray (arbitrary shape, int or bool type)
- The array containing the connected components of interest.
- max_size : int, optional (default: 64)
- Remove holes whose contiguous area (or volume, in N-D) contains this
- number of pixels or fewer.
- .. versionadded:: 0.26
- To make the naming clearer, replaces deprecated `area_threshold`
- which only removed holes strictly smaller than its size.
- connectivity : int, {1, 2, ..., ar.ndim}, optional (default: 1)
- The connectivity defining the neighborhood of a pixel.
- out : ndarray
- Array of the same shape as `ar` and bool dtype, into which the
- output is placed. By default, a new array is created.
- Raises
- ------
- TypeError
- If the input array is of an invalid type, such as float or string.
- ValueError
- If the input array contains negative values.
- Returns
- -------
- out : ndarray, same shape and type as input `ar`
- The input array with small holes within connected components removed.
- See Also
- --------
- skimage.morphology.remove_small_objects
- skimage.morphology.remove_objects_by_distance
- Examples
- --------
- >>> from skimage import morphology
- >>> a = np.array([[1, 1, 1, 1, 1, 0],
- ... [1, 1, 1, 0, 1, 0],
- ... [1, 0, 0, 1, 1, 0],
- ... [1, 1, 1, 1, 1, 0]], bool)
- >>> b = morphology.remove_small_holes(a, max_size=1)
- >>> b
- array([[ True, True, True, True, True, False],
- [ True, True, True, True, True, False],
- [ True, False, False, True, True, False],
- [ True, True, True, True, True, False]])
- >>> c = morphology.remove_small_holes(a, max_size=1, connectivity=2)
- >>> c
- array([[ True, True, True, True, True, False],
- [ True, True, True, False, True, False],
- [ True, False, False, True, True, False],
- [ True, True, True, True, True, False]])
- >>> d = morphology.remove_small_holes(a, max_size=1, out=a)
- >>> d is a
- True
- Notes
- -----
- If the array type is int, it is assumed that it contains already-labeled
- objects. The labels are not kept in the output image (this function always
- outputs a bool image). It is suggested that labeling is completed after
- using this function.
- """
- _check_dtype_supported(ar)
- # Creates warning if image is an integer image
- if ar.dtype != bool:
- warn(
- "Any labeled images will be returned as a boolean array. "
- "Did you mean to use a boolean array?",
- UserWarning,
- )
- if out is not None:
- if out.dtype != bool:
- raise TypeError("out dtype must be bool")
- else:
- out = ar.astype(bool, copy=True)
- # Creating the inverse of ar
- np.logical_not(ar, out=out)
- # removing small objects from the inverse of ar
- with warnings.catch_warnings():
- warnings.filterwarnings(
- "ignore",
- message="Parameter `min_size` is deprecated",
- category=FutureWarning,
- )
- out = remove_small_objects(
- out,
- min_size=area_threshold,
- max_size=max_size,
- connectivity=connectivity,
- out=out,
- )
- np.logical_not(out, out=out)
- return out
- def remove_objects_by_distance(
- label_image,
- min_distance,
- *,
- priority=None,
- p_norm=2,
- spacing=None,
- out=None,
- ):
- """Remove objects, in specified order, until remaining are a minimum distance apart.
- Remove labeled objects from an image until the remaining ones are spaced
- more than a given distance from one another. By default, smaller objects
- are removed first.
- Parameters
- ----------
- label_image : ndarray of integers
- An n-dimensional array containing object labels, e.g. as returned by
- :func:`~.label`. A value of zero is considered background, all other
- object IDs must be positive integers.
- min_distance : int or float
- Remove objects whose distance to other objects is not greater than this
- positive value. Objects with a lower `priority` are removed first.
- priority : ndarray, optional
- Defines the priority with which objects are removed. Expects a
- 1-dimensional array of length
- :func:`np.amax(label_image) + 1 <numpy.amax>` that contains the priority
- for each object's label at the respective index. Objects with a lower value
- are removed first until all remaining objects fulfill the distance
- requirement. If not given, priority is given to objects with a higher
- number of samples and their label value second.
- p_norm : int or float, optional
- The Minkowski distance of order p, used to calculate the distance
- between objects. The default ``2`` corresponds to the Euclidean
- distance, ``1`` to the "Manhattan" distance, and ``np.inf`` to the
- Chebyshev distance.
- spacing : sequence of float, optional
- The pixel spacing along each axis of `label_image`. If not specified,
- a grid spacing of unity (1) is implied.
- out : ndarray, optional
- Array of the same shape and dtype as `image`, into which the output is
- placed. By default, a new array is created.
- Returns
- -------
- out : ndarray
- Array of the same shape as `label_image`, for which objects that violate
- the `min_distance` condition were removed.
- See Also
- --------
- skimage.morphology.remove_small_objects
- Remove objects smaller than the specified size.
- skimage.morphology.remove_small_holes
- Remove holes smaller than the specified size.
- Notes
- -----
- The basic steps of this algorithm work as follows:
- 1. Find the indices for of all given objects and separate them depending on
- if they point to an object's border or not.
- 2. Sort indices by their label value, ensuring that indices which point to
- the same object are next to each other. This optimization allows finding
- all parts of an object, simply by stepping to the neighboring indices.
- 3. Sort boundary indices by `priority`. Use a stable-sort to preserve the
- ordering from the previous sorting step. If `priority` is not given,
- use :func:`numpy.bincount` as a fallback.
- 4. Construct a :class:`scipy.spatial.cKDTree` from the boundary indices.
- 5. Iterate across boundary indices in priority-sorted order, and query the
- kd-tree for objects that are too close. Remove ones that are and don't
- take them into account when evaluating other objects later on.
- The performance of this algorithm depends on the number of samples in
- `label_image` that belong to an object's border.
- Examples
- --------
- >>> import skimage as ski
- >>> ski.morphology.remove_objects_by_distance(np.array([2, 0, 1, 1]), 2)
- array([0, 0, 1, 1])
- >>> ski.morphology.remove_objects_by_distance(
- ... np.array([2, 0, 1, 1]), 2, priority=np.array([0, 1, 9])
- ... )
- array([2, 0, 0, 0])
- >>> label_image = np.array(
- ... [[8, 0, 0, 0, 0, 0, 0, 0, 0, 9, 9],
- ... [8, 8, 8, 0, 0, 0, 0, 0, 0, 9, 9],
- ... [0, 0, 0, 0, 0, 0, 0, 0, 9, 0, 0],
- ... [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
- ... [0, 0, 3, 0, 0, 0, 1, 0, 0, 0, 0],
- ... [2, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
- ... [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
- ... [0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 7]]
- ... )
- >>> ski.morphology.remove_objects_by_distance(
- ... label_image, min_distance=3
- ... )
- array([[8, 0, 0, 0, 0, 0, 0, 0, 0, 9, 9],
- [8, 8, 8, 0, 0, 0, 0, 0, 0, 9, 9],
- [0, 0, 0, 0, 0, 0, 0, 0, 9, 0, 0],
- [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
- [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
- [2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
- [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
- [0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 7]])
- """
- if min_distance < 0:
- raise ValueError(f"min_distance must be >= 0, was {min_distance}")
- if not np.issubdtype(label_image.dtype, np.integer):
- raise ValueError(
- f"`label_image` must be of integer dtype, got {label_image.dtype}"
- )
- if out is None:
- out = label_image.copy(order="C")
- elif out is not label_image:
- out[:] = label_image
- # May create a copy if order is not C, account for that later
- out_raveled = out.ravel(order="C")
- if spacing is not None:
- spacing = np.array(spacing)
- if spacing.shape != (out.ndim,) or spacing.min() <= 0:
- raise ValueError(
- "`spacing` must contain exactly one positive factor "
- "for each dimension of `label_image`"
- )
- indices = np.flatnonzero(out_raveled)
- # Optimization: Split indices into those on the object boundaries and inner
- # ones. The KDTree is built only from the boundary indices, which reduces
- # the size of the critical loop significantly! Remaining indices are only
- # used to remove the inner parts of objects as well.
- if (spacing is None or np.all(spacing[0] == spacing)) and p_norm <= 2:
- # For unity spacing we can make the borders more sparse by using a
- # lower connectivity
- footprint = ndi.generate_binary_structure(out.ndim, 1)
- else:
- footprint = ndi.generate_binary_structure(out.ndim, out.ndim)
- border = (
- ndi.maximum_filter(out, footprint=footprint)
- != ndi.minimum_filter(out, footprint=footprint)
- ).ravel()[indices]
- border_indices = indices[border]
- inner_indices = indices[~border]
- if border_indices.size == 0:
- # Image without any or only one object, return early
- return out
- # Sort by label ID first, so that IDs of the same object are contiguous
- # in the sorted index. This allows fast discovery of the whole object by
- # simple iteration up or down the index!
- border_indices = border_indices[np.argsort(out_raveled[border_indices])]
- inner_indices = inner_indices[np.argsort(out_raveled[inner_indices])]
- if priority is None:
- if not np.can_cast(out.dtype, np.intp, casting="safe"):
- # bincount expects intp (32-bit) on WASM or i386, so down-cast to that
- priority = np.bincount(out_raveled.astype(np.intp, copy=False))
- else:
- priority = np.bincount(out_raveled)
- # `priority` can only be indexed by positive object IDs,
- # `border_indices` contains all unique sorted IDs so check the lowest / first
- smallest_id = out_raveled[border_indices[0]]
- if smallest_id < 0:
- raise ValueError(f"found object with negative ID {smallest_id!r}")
- try:
- # Sort by priority second using a stable sort to preserve the contiguous
- # sorting of objects. Because each pixel in an object has the same
- # priority we don't need to worry about separating objects.
- border_indices = border_indices[
- np.argsort(priority[out_raveled[border_indices]], kind="stable")[::-1]
- ]
- except IndexError as error:
- # Use np.amax only for the exception path to provide a nicer error message
- expected_shape = (np.amax(out_raveled) + 1,)
- if priority.shape != expected_shape:
- raise ValueError(
- "shape of `priority` must be (np.amax(label_image) + 1,), "
- f"expected {expected_shape}, got {priority.shape} instead"
- ) from error
- else:
- raise
- # Construct kd-tree from unraveled border indices (optionally scale by `spacing`)
- unraveled_indices = np.unravel_index(border_indices, out.shape)
- if spacing is not None:
- unraveled_indices = tuple(
- unraveled_indices[dim] * spacing[dim] for dim in range(out.ndim)
- )
- kdtree = cKDTree(data=np.asarray(unraveled_indices, dtype=np.float64).T)
- _remove_objects_by_distance(
- out=out_raveled,
- border_indices=border_indices,
- inner_indices=inner_indices,
- kdtree=kdtree,
- min_distance=min_distance,
- p_norm=p_norm,
- shape=label_image.shape,
- )
- if out_raveled.base is not out:
- # `out_raveled` is a copy, re-assign
- out[:] = out_raveled.reshape(out.shape)
- return out
|