scikit-image

1
"""Utility functions used in the morphology subpackage."""
2

3
import numpy as np
4
from scipy import ndimage as ndi
5

6

7
def _validate_connectivity(image_dim, connectivity, offset):
8
    """Convert any valid connectivity to a footprint and offset.
9

10
    Parameters
11
    ----------
12
    image_dim : int
13
        The number of dimensions of the input image.
14
    connectivity : int, array, or None
15
        The neighborhood connectivity. An integer is interpreted as in
16
        ``scipy.ndimage.generate_binary_structure``, as the maximum number
17
        of orthogonal steps to reach a neighbor. An array is directly
18
        interpreted as a footprint and its shape is validated against
19
        the input image shape. ``None`` is interpreted as a connectivity of 1.
20
    offset : tuple of int, or None
21
        The coordinates of the center of the footprint.
22

23
    Returns
24
    -------
25
    c_connectivity : array of bool
26
        The footprint (structuring element) corresponding to the input
27
        `connectivity`.
28
    offset : array of int
29
        The offset corresponding to the center of the footprint.
30

31
    Raises
32
    ------
33
    ValueError:
34
        If the image dimension and the connectivity or offset dimensions don't
35
        match.
36
    """
37
    if connectivity is None:
38
        connectivity = 1
39

40
    if np.isscalar(connectivity):
41
        c_connectivity = ndi.generate_binary_structure(image_dim, connectivity)
42
    else:
43
        c_connectivity = np.array(connectivity, bool)
44
        if c_connectivity.ndim != image_dim:
45
            raise ValueError("Connectivity dimension must be same as image")
46

47
    if offset is None:
48
        if any([x % 2 == 0 for x in c_connectivity.shape]):
49
            raise ValueError("Connectivity array must have an unambiguous " "center")
50

51
        offset = np.array(c_connectivity.shape) // 2
52

53
    return c_connectivity, offset
54

55

56
def _raveled_offsets_and_distances(
57
    image_shape,
58
    *,
59
    footprint=None,
60
    connectivity=1,
61
    center=None,
62
    spacing=None,
63
    order='C',
64
):
65
    """Compute offsets to neighboring pixels in raveled coordinate space.
66

67
    This function also returns the corresponding distances from the center
68
    pixel given a spacing (assumed to be 1 along each axis by default).
69

70
    Parameters
71
    ----------
72
    image_shape : tuple of int
73
        The shape of the image for which the offsets are being computed.
74
    footprint : array of bool
75
        The footprint of the neighborhood, expressed as an n-dimensional array
76
        of 1s and 0s. If provided, the connectivity argument is ignored.
77
    connectivity : {1, ..., ndim}
78
        The square connectivity of the neighborhood: the number of orthogonal
79
        steps allowed to consider a pixel a neighbor. See
80
        `scipy.ndimage.generate_binary_structure`. Ignored if footprint is
81
        provided.
82
    center : tuple of int
83
        Tuple of indices to the center of the footprint. If not provided, it
84
        is assumed to be the center of the footprint, either provided or
85
        generated by the connectivity argument.
86
    spacing : tuple of float
87
        The spacing between pixels/voxels along each axis.
88
    order : 'C' or 'F'
89
        The ordering of the array, either C or Fortran ordering.
90

91
    Returns
92
    -------
93
    raveled_offsets : ndarray
94
        Linear offsets to a samples neighbors in the raveled image, sorted by
95
        their distance from the center.
96
    distances : ndarray
97
        The pixel distances corresponding to each offset.
98

99
    Notes
100
    -----
101
    This function will return values even if `image_shape` contains a dimension
102
    length that is smaller than `footprint`.
103

104
    Examples
105
    --------
106
    >>> off, d = _raveled_offsets_and_distances(
107
    ...         (4, 5), footprint=np.ones((4, 3)), center=(1, 1)
108
    ...         )
109
    >>> off
110
    array([-5, -1,  1,  5, -6, -4,  4,  6, 10,  9, 11])
111
    >>> d[0]
112
    1.0
113
    >>> d[-1]  # distance from (1, 1) to (3, 2)
114
    2.236...
115
    """
116
    ndim = len(image_shape)
117
    if footprint is None:
118
        footprint = ndi.generate_binary_structure(rank=ndim, connectivity=connectivity)
119
    if center is None:
120
        center = tuple(s // 2 for s in footprint.shape)
121

122
    if not footprint.ndim == ndim == len(center):
123
        raise ValueError(
124
            "number of dimensions in image shape, footprint and its"
125
            "center index does not match"
126
        )
127

128
    offsets = np.stack(
129
        [(idx - c) for idx, c in zip(np.nonzero(footprint), center)], axis=-1
130
    )
131

132
    if order == 'F':
133
        offsets = offsets[:, ::-1]
134
        image_shape = image_shape[::-1]
135
    elif order != 'C':
136
        raise ValueError("order must be 'C' or 'F'")
137

138
    # Scale offsets in each dimension and sum
139
    ravel_factors = image_shape[1:] + (1,)
140
    ravel_factors = np.cumprod(ravel_factors[::-1])[::-1]
141
    raveled_offsets = (offsets * ravel_factors).sum(axis=1)
142

143
    # Sort by distance
144
    if spacing is None:
145
        spacing = np.ones(ndim)
146
    weighted_offsets = offsets * spacing
147
    distances = np.sqrt(np.sum(weighted_offsets**2, axis=1))
148
    sorted_raveled_offsets = raveled_offsets[np.argsort(distances, kind="stable")]
149
    sorted_distances = np.sort(distances, kind="stable")
150

151
    # If any dimension in image_shape is smaller than footprint.shape
152
    # duplicates might occur, remove them
153
    if any(x < y for x, y in zip(image_shape, footprint.shape)):
154
        # np.unique reorders, which we don't want
155
        _, indices = np.unique(sorted_raveled_offsets, return_index=True)
156
        indices = np.sort(indices, kind="stable")
157
        sorted_raveled_offsets = sorted_raveled_offsets[indices]
158
        sorted_distances = sorted_distances[indices]
159

160
    # Remove "offset to center"
161
    sorted_raveled_offsets = sorted_raveled_offsets[1:]
162
    sorted_distances = sorted_distances[1:]
163

164
    return sorted_raveled_offsets, sorted_distances
165

166

167
def _offsets_to_raveled_neighbors(image_shape, footprint, center, order='C'):
168
    """Compute offsets to a samples neighbors if the image would be raveled.
169

170
    Parameters
171
    ----------
172
    image_shape : tuple
173
        The shape of the image for which the offsets are computed.
174
    footprint : ndarray
175
        The footprint (structuring element) determining the neighborhood
176
        expressed as an n-D array of 1's and 0's.
177
    center : tuple
178
        Tuple of indices to the center of `footprint`.
179
    order : {"C", "F"}, optional
180
        Whether the image described by `image_shape` is in row-major (C-style)
181
        or column-major (Fortran-style) order.
182

183
    Returns
184
    -------
185
    raveled_offsets : ndarray
186
        Linear offsets to a samples neighbors in the raveled image, sorted by
187
        their distance from the center.
188

189
    Notes
190
    -----
191
    This function will return values even if `image_shape` contains a dimension
192
    length that is smaller than `footprint`.
193

194
    Examples
195
    --------
196
    >>> _offsets_to_raveled_neighbors((4, 5), np.ones((4, 3)), (1, 1))
197
    array([-5, -1,  1,  5, -6, -4,  4,  6, 10,  9, 11])
198
    >>> _offsets_to_raveled_neighbors((2, 3, 2), np.ones((3, 3, 3)), (1, 1, 1))
199
    array([-6, -2, -1,  1,  2,  6, -8, -7, -5, -4, -3,  3,  4,  5,  7,  8, -9,
200
            9])
201
    """
202
    raveled_offsets = _raveled_offsets_and_distances(
203
        image_shape, footprint=footprint, center=center, order=order
204
    )[0]
205

206
    return raveled_offsets
207

208

209
def _resolve_neighborhood(footprint, connectivity, ndim, enforce_adjacency=True):
210
    """Validate or create a footprint (structuring element).
211

212
    Depending on the values of `connectivity` and `footprint` this function
213
    either creates a new footprint (`footprint` is None) using `connectivity`
214
    or validates the given footprint (`footprint` is not None).
215

216
    Parameters
217
    ----------
218
    footprint : ndarray
219
        The footprint (structuring) element used to determine the neighborhood
220
        of each evaluated pixel (``True`` denotes a connected pixel). It must
221
        be a boolean array and have the same number of dimensions as `image`.
222
        If neither `footprint` nor `connectivity` are given, all adjacent
223
        pixels are considered as part of the neighborhood.
224
    connectivity : int
225
        A number used to determine the neighborhood of each evaluated pixel.
226
        Adjacent pixels whose squared distance from the center is less than or
227
        equal to `connectivity` are considered neighbors. Ignored if
228
        `footprint` is not None.
229
    ndim : int
230
        Number of dimensions `footprint` ought to have.
231
    enforce_adjacency : bool
232
        A boolean that determines whether footprint must only specify direct
233
        neighbors.
234

235
    Returns
236
    -------
237
    footprint : ndarray
238
        Validated or new footprint specifying the neighborhood.
239

240
    Examples
241
    --------
242
    >>> _resolve_neighborhood(None, 1, 2)
243
    array([[False,  True, False],
244
           [ True,  True,  True],
245
           [False,  True, False]])
246
    >>> _resolve_neighborhood(None, None, 3).shape
247
    (3, 3, 3)
248
    """
249
    if footprint is None:
250
        if connectivity is None:
251
            connectivity = ndim
252
        footprint = ndi.generate_binary_structure(ndim, connectivity)
253
    else:
254
        # Validate custom structured element
255
        footprint = np.asarray(footprint, dtype=bool)
256
        # Must specify neighbors for all dimensions
257
        if footprint.ndim != ndim:
258
            raise ValueError(
259
                "number of dimensions in image and footprint do not" "match"
260
            )
261
        # Must only specify direct neighbors
262
        if enforce_adjacency and any(s != 3 for s in footprint.shape):
263
            raise ValueError("dimension size in footprint is not 3")
264
        elif any((s % 2 != 1) for s in footprint.shape):
265
            raise ValueError("footprint size must be odd along all dimensions")
266

267
    return footprint
268

269

270
def _set_border_values(image, value, border_width=1):
271
    """Set edge values along all axes to a constant value.
272

273
    Parameters
274
    ----------
275
    image : ndarray
276
        The array to modify inplace.
277
    value : scalar
278
        The value to use. Should be compatible with `image`'s dtype.
279
    border_width : int or sequence of tuples
280
        A sequence with one 2-tuple per axis where the first and second values
281
        are the width of the border at the start and end of the axis,
282
        respectively. If an int is provided, a uniform border width along all
283
        axes is used.
284

285
    Examples
286
    --------
287
    >>> image = np.zeros((4, 5), dtype=int)
288
    >>> _set_border_values(image, 1)
289
    >>> image
290
    array([[1, 1, 1, 1, 1],
291
           [1, 0, 0, 0, 1],
292
           [1, 0, 0, 0, 1],
293
           [1, 1, 1, 1, 1]])
294
    >>> image = np.zeros((8, 8), dtype=int)
295
    >>> _set_border_values(image, 1, border_width=((1, 1), (2, 3)))
296
    >>> image
297
    array([[1, 1, 1, 1, 1, 1, 1, 1],
298
           [1, 1, 0, 0, 0, 1, 1, 1],
299
           [1, 1, 0, 0, 0, 1, 1, 1],
300
           [1, 1, 0, 0, 0, 1, 1, 1],
301
           [1, 1, 0, 0, 0, 1, 1, 1],
302
           [1, 1, 0, 0, 0, 1, 1, 1],
303
           [1, 1, 0, 0, 0, 1, 1, 1],
304
           [1, 1, 1, 1, 1, 1, 1, 1]])
305
    """
306
    if np.isscalar(border_width):
307
        border_width = ((border_width, border_width),) * image.ndim
308
    elif len(border_width) != image.ndim:
309
        raise ValueError('length of `border_width` must match image.ndim')
310
    for axis, npad in enumerate(border_width):
311
        if len(npad) != 2:
312
            raise ValueError('each sequence in `border_width` must have ' 'length 2')
313
        w_start, w_end = npad
314
        if w_start == w_end == 0:
315
            continue
316
        elif w_start == w_end == 1:
317
            # Index first and last element in the current dimension
318
            sl = (slice(None),) * axis + ((0, -1),) + (...,)
319
            image[sl] = value
320
            continue
321
        if w_start > 0:
322
            # set first w_start entries along axis to value
323
            sl = (slice(None),) * axis + (slice(0, w_start),) + (...,)
324
            image[sl] = value
325
        if w_end > 0:
326
            # set last w_end entries along axis to value
327
            sl = (slice(None),) * axis + (slice(-w_end, None),) + (...,)
328
            image[sl] = value
329