scikit-image

1
"""Analytical transformations from raw image moments to central moments.
2

3
The expressions for the 2D central moments of order <=2 are often given in
4
textbooks. Expressions for higher orders and dimensions were generated in SymPy
5
using ``tools/precompute/moments_sympy.py`` in the GitHub repository.
6

7
"""
8

9
import itertools
10
import math
11

12
import numpy as np
13

14

15
def _moments_raw_to_central_fast(moments_raw):
16
    """Analytical formulae for 2D and 3D central moments of order < 4.
17

18
    `moments_raw_to_central` will automatically call this function when
19
    ndim < 4 and order < 4.
20

21
    Parameters
22
    ----------
23
    moments_raw : ndarray
24
        The raw moments.
25

26
    Returns
27
    -------
28
    moments_central : ndarray
29
        The central moments.
30
    """
31
    ndim = moments_raw.ndim
32
    order = moments_raw.shape[0] - 1
33
    float_dtype = moments_raw.dtype
34
    # convert to float64 during the computation for better accuracy
35
    moments_raw = moments_raw.astype(np.float64, copy=False)
36
    moments_central = np.zeros_like(moments_raw)
37
    if order >= 4 or ndim not in [2, 3]:
38
        raise ValueError("This function only supports 2D or 3D moments of order < 4.")
39
    m = moments_raw
40
    if ndim == 2:
41
        cx = m[1, 0] / m[0, 0]
42
        cy = m[0, 1] / m[0, 0]
43
        moments_central[0, 0] = m[0, 0]
44
        # Note: 1st order moments are both 0
45
        if order > 1:
46
            # 2nd order moments
47
            moments_central[1, 1] = m[1, 1] - cx * m[0, 1]
48
            moments_central[2, 0] = m[2, 0] - cx * m[1, 0]
49
            moments_central[0, 2] = m[0, 2] - cy * m[0, 1]
50
        if order > 2:
51
            # 3rd order moments
52
            moments_central[2, 1] = (
53
                m[2, 1]
54
                - 2 * cx * m[1, 1]
55
                - cy * m[2, 0]
56
                + cx**2 * m[0, 1]
57
                + cy * cx * m[1, 0]
58
            )
59
            moments_central[1, 2] = (
60
                m[1, 2] - 2 * cy * m[1, 1] - cx * m[0, 2] + 2 * cy * cx * m[0, 1]
61
            )
62
            moments_central[3, 0] = m[3, 0] - 3 * cx * m[2, 0] + 2 * cx**2 * m[1, 0]
63
            moments_central[0, 3] = m[0, 3] - 3 * cy * m[0, 2] + 2 * cy**2 * m[0, 1]
64
    else:
65
        # 3D case
66
        cx = m[1, 0, 0] / m[0, 0, 0]
67
        cy = m[0, 1, 0] / m[0, 0, 0]
68
        cz = m[0, 0, 1] / m[0, 0, 0]
69
        moments_central[0, 0, 0] = m[0, 0, 0]
70
        # Note: all first order moments are 0
71
        if order > 1:
72
            # 2nd order moments
73
            moments_central[0, 0, 2] = -cz * m[0, 0, 1] + m[0, 0, 2]
74
            moments_central[0, 1, 1] = -cy * m[0, 0, 1] + m[0, 1, 1]
75
            moments_central[0, 2, 0] = -cy * m[0, 1, 0] + m[0, 2, 0]
76
            moments_central[1, 0, 1] = -cx * m[0, 0, 1] + m[1, 0, 1]
77
            moments_central[1, 1, 0] = -cx * m[0, 1, 0] + m[1, 1, 0]
78
            moments_central[2, 0, 0] = -cx * m[1, 0, 0] + m[2, 0, 0]
79
        if order > 2:
80
            # 3rd order moments
81
            moments_central[0, 0, 3] = (
82
                2 * cz**2 * m[0, 0, 1] - 3 * cz * m[0, 0, 2] + m[0, 0, 3]
83
            )
84
            moments_central[0, 1, 2] = (
85
                -cy * m[0, 0, 2] + 2 * cz * (cy * m[0, 0, 1] - m[0, 1, 1]) + m[0, 1, 2]
86
            )
87
            moments_central[0, 2, 1] = (
88
                cy**2 * m[0, 0, 1]
89
                - 2 * cy * m[0, 1, 1]
90
                + cz * (cy * m[0, 1, 0] - m[0, 2, 0])
91
                + m[0, 2, 1]
92
            )
93
            moments_central[0, 3, 0] = (
94
                2 * cy**2 * m[0, 1, 0] - 3 * cy * m[0, 2, 0] + m[0, 3, 0]
95
            )
96
            moments_central[1, 0, 2] = (
97
                -cx * m[0, 0, 2] + 2 * cz * (cx * m[0, 0, 1] - m[1, 0, 1]) + m[1, 0, 2]
98
            )
99
            moments_central[1, 1, 1] = (
100
                -cx * m[0, 1, 1]
101
                + cy * (cx * m[0, 0, 1] - m[1, 0, 1])
102
                + cz * (cx * m[0, 1, 0] - m[1, 1, 0])
103
                + m[1, 1, 1]
104
            )
105
            moments_central[1, 2, 0] = (
106
                -cx * m[0, 2, 0] - 2 * cy * (-cx * m[0, 1, 0] + m[1, 1, 0]) + m[1, 2, 0]
107
            )
108
            moments_central[2, 0, 1] = (
109
                cx**2 * m[0, 0, 1]
110
                - 2 * cx * m[1, 0, 1]
111
                + cz * (cx * m[1, 0, 0] - m[2, 0, 0])
112
                + m[2, 0, 1]
113
            )
114
            moments_central[2, 1, 0] = (
115
                cx**2 * m[0, 1, 0]
116
                - 2 * cx * m[1, 1, 0]
117
                + cy * (cx * m[1, 0, 0] - m[2, 0, 0])
118
                + m[2, 1, 0]
119
            )
120
            moments_central[3, 0, 0] = (
121
                2 * cx**2 * m[1, 0, 0] - 3 * cx * m[2, 0, 0] + m[3, 0, 0]
122
            )
123

124
    return moments_central.astype(float_dtype, copy=False)
125

126

127
def moments_raw_to_central(moments_raw):
128
    ndim = moments_raw.ndim
129
    order = moments_raw.shape[0] - 1
130
    if ndim in [2, 3] and order < 4:
131
        return _moments_raw_to_central_fast(moments_raw)
132

133
    moments_central = np.zeros_like(moments_raw)
134
    m = moments_raw
135
    # centers as computed in centroid above
136
    centers = tuple(m[tuple(np.eye(ndim, dtype=int))] / m[(0,) * ndim])
137

138
    if ndim == 2:
139
        # This is the general 2D formula from
140
        # https://en.wikipedia.org/wiki/Image_moment#Central_moments
141
        for p in range(order + 1):
142
            for q in range(order + 1):
143
                if p + q > order:
144
                    continue
145
                for i in range(p + 1):
146
                    term1 = math.comb(p, i)
147
                    term1 *= (-centers[0]) ** (p - i)
148
                    for j in range(q + 1):
149
                        term2 = math.comb(q, j)
150
                        term2 *= (-centers[1]) ** (q - j)
151
                        moments_central[p, q] += term1 * term2 * m[i, j]
152
        return moments_central
153

154
    # The nested loops below are an n-dimensional extension of the 2D formula
155
    # given at https://en.wikipedia.org/wiki/Image_moment#Central_moments
156

157
    # iterate over all [0, order] (inclusive) on each axis
158
    for orders in itertools.product(*((range(order + 1),) * ndim)):
159
        # `orders` here is the index into the `moments_central` output array
160
        if sum(orders) > order:
161
            # skip any moment that is higher than the requested order
162
            continue
163
        # loop over terms from `m` contributing to `moments_central[orders]`
164
        for idxs in itertools.product(*[range(o + 1) for o in orders]):
165
            val = m[idxs]
166
            for i_order, c, idx in zip(orders, centers, idxs):
167
                val *= math.comb(i_order, idx)
168
                val *= (-c) ** (i_order - idx)
169
            moments_central[orders] += val
170

171
    return moments_central
172