scikit-image

1
import itertools
2

3
import numpy as np
4

5
from .._shared.utils import _supported_float_type, check_nD
6
from . import _moments_cy
7
from ._moments_analytical import moments_raw_to_central
8

9

10
def moments_coords(coords, order=3):
11
    """Calculate all raw image moments up to a certain order.
12

13
    The following properties can be calculated from raw image moments:
14
     * Area as: ``M[0, 0]``.
15
     * Centroid as: {``M[1, 0] / M[0, 0]``, ``M[0, 1] / M[0, 0]``}.
16

17
    Note that raw moments are neither translation, scale, nor rotation
18
    invariant.
19

20
    Parameters
21
    ----------
22
    coords : (N, D) double or uint8 array
23
        Array of N points that describe an image of D dimensionality in
24
        Cartesian space.
25
    order : int, optional
26
        Maximum order of moments. Default is 3.
27

28
    Returns
29
    -------
30
    M : (``order + 1``, ``order + 1``, ...) array
31
        Raw image moments. (D dimensions)
32

33
    References
34
    ----------
35
    .. [1] Johannes Kilian. Simple Image Analysis By Moments. Durham
36
           University, version 0.2, Durham, 2001.
37

38
    Examples
39
    --------
40
    >>> coords = np.array([[row, col]
41
    ...                    for row in range(13, 17)
42
    ...                    for col in range(14, 18)], dtype=np.float64)
43
    >>> M = moments_coords(coords)
44
    >>> centroid = (M[1, 0] / M[0, 0], M[0, 1] / M[0, 0])
45
    >>> centroid
46
    (14.5, 15.5)
47
    """
48
    return moments_coords_central(coords, 0, order=order)
49

50

51
def moments_coords_central(coords, center=None, order=3):
52
    """Calculate all central image moments up to a certain order.
53

54
    The following properties can be calculated from raw image moments:
55
     * Area as: ``M[0, 0]``.
56
     * Centroid as: {``M[1, 0] / M[0, 0]``, ``M[0, 1] / M[0, 0]``}.
57

58
    Note that raw moments are neither translation, scale nor rotation
59
    invariant.
60

61
    Parameters
62
    ----------
63
    coords : (N, D) double or uint8 array
64
        Array of N points that describe an image of D dimensionality in
65
        Cartesian space. A tuple of coordinates as returned by
66
        ``np.nonzero`` is also accepted as input.
67
    center : tuple of float, optional
68
        Coordinates of the image centroid. This will be computed if it
69
        is not provided.
70
    order : int, optional
71
        Maximum order of moments. Default is 3.
72

73
    Returns
74
    -------
75
    Mc : (``order + 1``, ``order + 1``, ...) array
76
        Central image moments. (D dimensions)
77

78
    References
79
    ----------
80
    .. [1] Johannes Kilian. Simple Image Analysis By Moments. Durham
81
           University, version 0.2, Durham, 2001.
82

83
    Examples
84
    --------
85
    >>> coords = np.array([[row, col]
86
    ...                    for row in range(13, 17)
87
    ...                    for col in range(14, 18)])
88
    >>> moments_coords_central(coords)
89
    array([[16.,  0., 20.,  0.],
90
           [ 0.,  0.,  0.,  0.],
91
           [20.,  0., 25.,  0.],
92
           [ 0.,  0.,  0.,  0.]])
93

94
    As seen above, for symmetric objects, odd-order moments (columns 1 and 3,
95
    rows 1 and 3) are zero when centered on the centroid, or center of mass,
96
    of the object (the default). If we break the symmetry by adding a new
97
    point, this no longer holds:
98

99
    >>> coords2 = np.concatenate((coords, [[17, 17]]), axis=0)
100
    >>> np.round(moments_coords_central(coords2),
101
    ...          decimals=2)  # doctest: +NORMALIZE_WHITESPACE
102
    array([[17.  ,  0.  , 22.12, -2.49],
103
           [ 0.  ,  3.53,  1.73,  7.4 ],
104
           [25.88,  6.02, 36.63,  8.83],
105
           [ 4.15, 19.17, 14.8 , 39.6 ]])
106

107
    Image moments and central image moments are equivalent (by definition)
108
    when the center is (0, 0):
109

110
    >>> np.allclose(moments_coords(coords),
111
    ...             moments_coords_central(coords, (0, 0)))
112
    True
113
    """
114
    if isinstance(coords, tuple):
115
        # This format corresponds to coordinate tuples as returned by
116
        # e.g. np.nonzero: (row_coords, column_coords).
117
        # We represent them as an npoints x ndim array.
118
        coords = np.stack(coords, axis=-1)
119
    check_nD(coords, 2)
120
    ndim = coords.shape[1]
121

122
    float_type = _supported_float_type(coords.dtype)
123
    if center is None:
124
        center = np.mean(coords, axis=0, dtype=float)
125

126
    # center the coordinates
127
    coords = coords.astype(float_type, copy=False) - center
128

129
    # generate all possible exponents for each axis in the given set of points
130
    # produces a matrix of shape (N, D, order + 1)
131
    coords = np.stack([coords**c for c in range(order + 1)], axis=-1)
132

133
    # add extra dimensions for proper broadcasting
134
    coords = coords.reshape(coords.shape + (1,) * (ndim - 1))
135

136
    calc = 1
137

138
    for axis in range(ndim):
139
        # isolate each point's axis
140
        isolated_axis = coords[:, axis]
141

142
        # rotate orientation of matrix for proper broadcasting
143
        isolated_axis = np.moveaxis(isolated_axis, 1, 1 + axis)
144

145
        # calculate the moments for each point, one axis at a time
146
        calc = calc * isolated_axis
147

148
    # sum all individual point moments to get our final answer
149
    Mc = np.sum(calc, axis=0)
150

151
    return Mc
152

153

154
def moments(image, order=3, *, spacing=None):
155
    """Calculate all raw image moments up to a certain order.
156

157
    The following properties can be calculated from raw image moments:
158
     * Area as: ``M[0, 0]``.
159
     * Centroid as: {``M[1, 0] / M[0, 0]``, ``M[0, 1] / M[0, 0]``}.
160

161
    Note that raw moments are neither translation, scale nor rotation
162
    invariant.
163

164
    Parameters
165
    ----------
166
    image : (N[, ...]) double or uint8 array
167
        Rasterized shape as image.
168
    order : int, optional
169
        Maximum order of moments. Default is 3.
170
    spacing: tuple of float, shape (ndim,)
171
        The pixel spacing along each axis of the image.
172

173
    Returns
174
    -------
175
    m : (``order + 1``, ``order + 1``) array
176
        Raw image moments.
177

178
    References
179
    ----------
180
    .. [1] Wilhelm Burger, Mark Burge. Principles of Digital Image Processing:
181
           Core Algorithms. Springer-Verlag, London, 2009.
182
    .. [2] B. Jähne. Digital Image Processing. Springer-Verlag,
183
           Berlin-Heidelberg, 6. edition, 2005.
184
    .. [3] T. H. Reiss. Recognizing Planar Objects Using Invariant Image
185
           Features, from Lecture notes in computer science, p. 676. Springer,
186
           Berlin, 1993.
187
    .. [4] https://en.wikipedia.org/wiki/Image_moment
188

189
    Examples
190
    --------
191
    >>> image = np.zeros((20, 20), dtype=np.float64)
192
    >>> image[13:17, 13:17] = 1
193
    >>> M = moments(image)
194
    >>> centroid = (M[1, 0] / M[0, 0], M[0, 1] / M[0, 0])
195
    >>> centroid
196
    (14.5, 14.5)
197
    """
198
    return moments_central(image, (0,) * image.ndim, order=order, spacing=spacing)
199

200

201
def moments_central(image, center=None, order=3, *, spacing=None, **kwargs):
202
    """Calculate all central image moments up to a certain order.
203

204
    The center coordinates (cr, cc) can be calculated from the raw moments as:
205
    {``M[1, 0] / M[0, 0]``, ``M[0, 1] / M[0, 0]``}.
206

207
    Note that central moments are translation invariant but not scale and
208
    rotation invariant.
209

210
    Parameters
211
    ----------
212
    image : (N[, ...]) double or uint8 array
213
        Rasterized shape as image.
214
    center : tuple of float, optional
215
        Coordinates of the image centroid. This will be computed if it
216
        is not provided.
217
    order : int, optional
218
        The maximum order of moments computed.
219
    spacing: tuple of float, shape (ndim,)
220
        The pixel spacing along each axis of the image.
221

222
    Returns
223
    -------
224
    mu : (``order + 1``, ``order + 1``) array
225
        Central image moments.
226

227
    References
228
    ----------
229
    .. [1] Wilhelm Burger, Mark Burge. Principles of Digital Image Processing:
230
           Core Algorithms. Springer-Verlag, London, 2009.
231
    .. [2] B. Jähne. Digital Image Processing. Springer-Verlag,
232
           Berlin-Heidelberg, 6. edition, 2005.
233
    .. [3] T. H. Reiss. Recognizing Planar Objects Using Invariant Image
234
           Features, from Lecture notes in computer science, p. 676. Springer,
235
           Berlin, 1993.
236
    .. [4] https://en.wikipedia.org/wiki/Image_moment
237

238
    Examples
239
    --------
240
    >>> image = np.zeros((20, 20), dtype=np.float64)
241
    >>> image[13:17, 13:17] = 1
242
    >>> M = moments(image)
243
    >>> centroid = (M[1, 0] / M[0, 0], M[0, 1] / M[0, 0])
244
    >>> moments_central(image, centroid)
245
    array([[16.,  0., 20.,  0.],
246
           [ 0.,  0.,  0.,  0.],
247
           [20.,  0., 25.,  0.],
248
           [ 0.,  0.,  0.,  0.]])
249
    """
250
    if center is None:
251
        # Note: No need for an explicit call to centroid.
252
        #       The centroid will be obtained from the raw moments.
253
        moments_raw = moments(image, order=order, spacing=spacing)
254
        return moments_raw_to_central(moments_raw)
255
    float_dtype = _supported_float_type(image.dtype)
256
    if spacing is None:
257
        spacing = np.ones(image.ndim, dtype=float_dtype)
258
    calc = image.astype(float_dtype, copy=False)
259
    for dim, dim_length in enumerate(image.shape):
260
        delta = np.arange(dim_length, dtype=float_dtype) * spacing[dim] - center[dim]
261
        powers_of_delta = delta[:, np.newaxis] ** np.arange(
262
            order + 1, dtype=float_dtype
263
        )
264
        calc = np.rollaxis(calc, dim, image.ndim)
265
        calc = np.dot(calc, powers_of_delta)
266
        calc = np.rollaxis(calc, -1, dim)
267
    return calc
268

269

270
def moments_normalized(mu, order=3, spacing=None):
271
    """Calculate all normalized central image moments up to a certain order.
272

273
    Note that normalized central moments are translation and scale invariant
274
    but not rotation invariant.
275

276
    Parameters
277
    ----------
278
    mu : (M[, ...], M) array
279
        Central image moments, where M must be greater than or equal
280
        to ``order``.
281
    order : int, optional
282
        Maximum order of moments. Default is 3.
283
    spacing: tuple of float, shape (ndim,)
284
        The pixel spacing along each axis of the image.
285

286
    Returns
287
    -------
288
    nu : (``order + 1``[, ...], ``order + 1``) array
289
        Normalized central image moments.
290

291
    References
292
    ----------
293
    .. [1] Wilhelm Burger, Mark Burge. Principles of Digital Image Processing:
294
           Core Algorithms. Springer-Verlag, London, 2009.
295
    .. [2] B. Jähne. Digital Image Processing. Springer-Verlag,
296
           Berlin-Heidelberg, 6. edition, 2005.
297
    .. [3] T. H. Reiss. Recognizing Planar Objects Using Invariant Image
298
           Features, from Lecture notes in computer science, p. 676. Springer,
299
           Berlin, 1993.
300
    .. [4] https://en.wikipedia.org/wiki/Image_moment
301

302
    Examples
303
    --------
304
    >>> image = np.zeros((20, 20), dtype=np.float64)
305
    >>> image[13:17, 13:17] = 1
306
    >>> m = moments(image)
307
    >>> centroid = (m[0, 1] / m[0, 0], m[1, 0] / m[0, 0])
308
    >>> mu = moments_central(image, centroid)
309
    >>> moments_normalized(mu)
310
    array([[       nan,        nan, 0.078125  , 0.        ],
311
           [       nan, 0.        , 0.        , 0.        ],
312
           [0.078125  , 0.        , 0.00610352, 0.        ],
313
           [0.        , 0.        , 0.        , 0.        ]])
314
    """
315
    if np.any(np.array(mu.shape) <= order):
316
        raise ValueError("Shape of image moments must be >= `order`")
317
    if spacing is None:
318
        spacing = np.ones(mu.ndim)
319
    nu = np.zeros_like(mu)
320
    mu0 = mu.ravel()[0]
321
    scale = min(spacing)
322
    for powers in itertools.product(range(order + 1), repeat=mu.ndim):
323
        if sum(powers) < 2:
324
            nu[powers] = np.nan
325
        else:
326
            nu[powers] = (mu[powers] / scale ** sum(powers)) / (
327
                mu0 ** (sum(powers) / nu.ndim + 1)
328
            )
329
    return nu
330

331

332
def moments_hu(nu):
333
    """Calculate Hu's set of image moments (2D-only).
334

335
    Note that this set of moments is proved to be translation, scale and
336
    rotation invariant.
337

338
    Parameters
339
    ----------
340
    nu : (M, M) array
341
        Normalized central image moments, where M must be >= 4.
342

343
    Returns
344
    -------
345
    nu : (7,) array
346
        Hu's set of image moments.
347

348
    References
349
    ----------
350
    .. [1] M. K. Hu, "Visual Pattern Recognition by Moment Invariants",
351
           IRE Trans. Info. Theory, vol. IT-8, pp. 179-187, 1962
352
    .. [2] Wilhelm Burger, Mark Burge. Principles of Digital Image Processing:
353
           Core Algorithms. Springer-Verlag, London, 2009.
354
    .. [3] B. Jähne. Digital Image Processing. Springer-Verlag,
355
           Berlin-Heidelberg, 6. edition, 2005.
356
    .. [4] T. H. Reiss. Recognizing Planar Objects Using Invariant Image
357
           Features, from Lecture notes in computer science, p. 676. Springer,
358
           Berlin, 1993.
359
    .. [5] https://en.wikipedia.org/wiki/Image_moment
360

361
    Examples
362
    --------
363
    >>> image = np.zeros((20, 20), dtype=np.float64)
364
    >>> image[13:17, 13:17] = 0.5
365
    >>> image[10:12, 10:12] = 1
366
    >>> mu = moments_central(image)
367
    >>> nu = moments_normalized(mu)
368
    >>> moments_hu(nu)
369
    array([0.74537037, 0.35116598, 0.10404918, 0.04064421, 0.00264312,
370
           0.02408546, 0.        ])
371
    """
372
    dtype = np.float32 if nu.dtype == 'float32' else np.float64
373
    return _moments_cy.moments_hu(nu.astype(dtype, copy=False))
374

375

376
def centroid(image, *, spacing=None):
377
    """Return the (weighted) centroid of an image.
378

379
    Parameters
380
    ----------
381
    image : array
382
        The input image.
383
    spacing: tuple of float, shape (ndim,)
384
        The pixel spacing along each axis of the image.
385

386
    Returns
387
    -------
388
    center : tuple of float, length ``image.ndim``
389
        The centroid of the (nonzero) pixels in ``image``.
390

391
    Examples
392
    --------
393
    >>> image = np.zeros((20, 20), dtype=np.float64)
394
    >>> image[13:17, 13:17] = 0.5
395
    >>> image[10:12, 10:12] = 1
396
    >>> centroid(image)
397
    array([13.16666667, 13.16666667])
398
    """
399
    M = moments_central(image, center=(0,) * image.ndim, order=1, spacing=spacing)
400
    center = (
401
        M[tuple(np.eye(image.ndim, dtype=int))]  # array of weighted sums
402
        # for each axis
403
        / M[(0,) * image.ndim]
404
    )  # weighted sum of all points
405
    return center
406

407

408
def inertia_tensor(image, mu=None, *, spacing=None):
409
    """Compute the inertia tensor of the input image.
410

411
    Parameters
412
    ----------
413
    image : array
414
        The input image.
415
    mu : array, optional
416
        The pre-computed central moments of ``image``. The inertia tensor
417
        computation requires the central moments of the image. If an
418
        application requires both the central moments and the inertia tensor
419
        (for example, `skimage.measure.regionprops`), then it is more
420
        efficient to pre-compute them and pass them to the inertia tensor
421
        call.
422
    spacing: tuple of float, shape (ndim,)
423
        The pixel spacing along each axis of the image.
424

425
    Returns
426
    -------
427
    T : array, shape ``(image.ndim, image.ndim)``
428
        The inertia tensor of the input image. :math:`T_{i, j}` contains
429
        the covariance of image intensity along axes :math:`i` and :math:`j`.
430

431
    References
432
    ----------
433
    .. [1] https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor
434
    .. [2] Bernd Jähne. Spatio-Temporal Image Processing: Theory and
435
           Scientific Applications. (Chapter 8: Tensor Methods) Springer, 1993.
436
    """
437
    if mu is None:
438
        mu = moments_central(
439
            image, order=2, spacing=spacing
440
        )  # don't need higher-order moments
441
    mu0 = mu[(0,) * image.ndim]
442
    result = np.zeros((image.ndim, image.ndim), dtype=mu.dtype)
443

444
    # nD expression to get coordinates ([2, 0], [0, 2]) (2D),
445
    # ([2, 0, 0], [0, 2, 0], [0, 0, 2]) (3D), etc.
446
    corners2 = tuple(2 * np.eye(image.ndim, dtype=int))
447
    d = np.diag(result)
448
    d.flags.writeable = True
449
    # See https://ocw.mit.edu/courses/aeronautics-and-astronautics/
450
    #             16-07-dynamics-fall-2009/lecture-notes/MIT16_07F09_Lec26.pdf
451
    # Iii is the sum of second-order moments of every axis *except* i, not the
452
    # second order moment of axis i.
453
    # See also https://github.com/scikit-image/scikit-image/issues/3229
454
    d[:] = (np.sum(mu[corners2]) - mu[corners2]) / mu0
455

456
    for dims in itertools.combinations(range(image.ndim), 2):
457
        mu_index = np.zeros(image.ndim, dtype=int)
458
        mu_index[list(dims)] = 1
459
        result[dims] = -mu[tuple(mu_index)] / mu0
460
        result.T[dims] = -mu[tuple(mu_index)] / mu0
461
    return result
462

463

464
def inertia_tensor_eigvals(image, mu=None, T=None, *, spacing=None):
465
    """Compute the eigenvalues of the inertia tensor of the image.
466

467
    The inertia tensor measures covariance of the image intensity along
468
    the image axes. (See `inertia_tensor`.) The relative magnitude of the
469
    eigenvalues of the tensor is thus a measure of the elongation of a
470
    (bright) object in the image.
471

472
    Parameters
473
    ----------
474
    image : array
475
        The input image.
476
    mu : array, optional
477
        The pre-computed central moments of ``image``.
478
    T : array, shape ``(image.ndim, image.ndim)``
479
        The pre-computed inertia tensor. If ``T`` is given, ``mu`` and
480
        ``image`` are ignored.
481
    spacing: tuple of float, shape (ndim,)
482
        The pixel spacing along each axis of the image.
483

484
    Returns
485
    -------
486
    eigvals : list of float, length ``image.ndim``
487
        The eigenvalues of the inertia tensor of ``image``, in descending
488
        order.
489

490
    Notes
491
    -----
492
    Computing the eigenvalues requires the inertia tensor of the input image.
493
    This is much faster if the central moments (``mu``) are provided, or,
494
    alternatively, one can provide the inertia tensor (``T``) directly.
495
    """
496
    if T is None:
497
        T = inertia_tensor(image, mu, spacing=spacing)
498
    eigvals = np.linalg.eigvalsh(T)
499
    # Floating point precision problems could make a positive
500
    # semidefinite matrix have an eigenvalue that is very slightly
501
    # negative. This can cause problems down the line, so set values
502
    # very near zero to zero.
503
    eigvals = np.clip(eigvals, 0, None, out=eigvals)
504
    return sorted(eigvals, reverse=True)
505