scikit-image

1
import itertools
2

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

7
from skimage import draw
8
from skimage._shared import testing
9
from skimage._shared.testing import assert_allclose, assert_almost_equal, assert_equal
10
from skimage._shared.utils import _supported_float_type
11
from skimage.measure import (
12
    centroid,
13
    inertia_tensor,
14
    inertia_tensor_eigvals,
15
    moments,
16
    moments_central,
17
    moments_coords,
18
    moments_coords_central,
19
    moments_hu,
20
    moments_normalized,
21
)
22

23

24
def compare_moments(m1, m2, thresh=1e-8):
25
    """Compare two moments arrays.
26

27
    Compares only values in the upper-left triangle of m1, m2 since
28
    values below the diagonal exceed the specified order and are not computed
29
    when the analytical computation is used.
30

31
    Also, there the first-order central moments will be exactly zero with the
32
    analytical calculation, but will not be zero due to limited floating point
33
    precision when using a numerical computation. Here we just specify the
34
    tolerance as a fraction of the maximum absolute value in the moments array.
35
    """
36
    m1 = m1.copy()
37
    m2 = m2.copy()
38

39
    # make sure location of any NaN values match and then ignore the NaN values
40
    # in the subsequent comparisons
41
    nan_idx1 = np.where(np.isnan(m1.ravel()))[0]
42
    nan_idx2 = np.where(np.isnan(m2.ravel()))[0]
43
    assert len(nan_idx1) == len(nan_idx2)
44
    assert np.all(nan_idx1 == nan_idx2)
45
    m1[np.isnan(m1)] = 0
46
    m2[np.isnan(m2)] = 0
47

48
    max_val = np.abs(m1[m1 != 0]).max()
49
    for orders in itertools.product(*((range(m1.shape[0]),) * m1.ndim)):
50
        if sum(orders) > m1.shape[0] - 1:
51
            m1[orders] = 0
52
            m2[orders] = 0
53
            continue
54
        abs_diff = abs(m1[orders] - m2[orders])
55
        rel_diff = abs_diff / max_val
56
        assert rel_diff < thresh
57

58

59
@pytest.mark.parametrize('anisotropic', [False, True, None])
60
def test_moments(anisotropic):
61
    image = np.zeros((20, 20), dtype=np.float64)
62
    image[14, 14] = 1
63
    image[15, 15] = 1
64
    image[14, 15] = 0.5
65
    image[15, 14] = 0.5
66
    if anisotropic:
67
        spacing = (1.4, 2)
68
    else:
69
        spacing = (1, 1)
70
    if anisotropic is None:
71
        m = moments(image)
72
    else:
73
        m = moments(image, spacing=spacing)
74
    assert_equal(m[0, 0], 3)
75
    assert_almost_equal(m[1, 0] / m[0, 0], 14.5 * spacing[0])
76
    assert_almost_equal(m[0, 1] / m[0, 0], 14.5 * spacing[1])
77

78

79
@pytest.mark.parametrize('anisotropic', [False, True, None])
80
def test_moments_central(anisotropic):
81
    image = np.zeros((20, 20), dtype=np.float64)
82
    image[14, 14] = 1
83
    image[15, 15] = 1
84
    image[14, 15] = 0.5
85
    image[15, 14] = 0.5
86
    if anisotropic:
87
        spacing = (2, 1)
88
    else:
89
        spacing = (1, 1)
90
    if anisotropic is None:
91
        mu = moments_central(image, (14.5, 14.5))
92
        # check for proper centroid computation
93
        mu_calc_centroid = moments_central(image)
94
    else:
95
        mu = moments_central(
96
            image, (14.5 * spacing[0], 14.5 * spacing[1]), spacing=spacing
97
        )
98
        # check for proper centroid computation
99
        mu_calc_centroid = moments_central(image, spacing=spacing)
100

101
    compare_moments(mu, mu_calc_centroid, thresh=1e-14)
102

103
    # shift image by dx=2, dy=2
104
    image2 = np.zeros((20, 20), dtype=np.double)
105
    image2[16, 16] = 1
106
    image2[17, 17] = 1
107
    image2[16, 17] = 0.5
108
    image2[17, 16] = 0.5
109
    if anisotropic is None:
110
        mu2 = moments_central(image2, (14.5 + 2, 14.5 + 2))
111
    else:
112
        mu2 = moments_central(
113
            image2, ((14.5 + 2) * spacing[0], (14.5 + 2) * spacing[1]), spacing=spacing
114
        )
115
    # central moments must be translation invariant
116
    compare_moments(mu, mu2, thresh=1e-14)
117

118

119
def test_moments_coords():
120
    image = np.zeros((20, 20), dtype=np.float64)
121
    image[13:17, 13:17] = 1
122
    mu_image = moments(image)
123

124
    coords = np.array(
125
        [[r, c] for r in range(13, 17) for c in range(13, 17)], dtype=np.float64
126
    )
127
    mu_coords = moments_coords(coords)
128
    assert_almost_equal(mu_coords, mu_image)
129

130

131
@pytest.mark.parametrize('dtype', [np.float16, np.float32, np.float64])
132
def test_moments_coords_dtype(dtype):
133
    image = np.zeros((20, 20), dtype=dtype)
134
    image[13:17, 13:17] = 1
135

136
    expected_dtype = _supported_float_type(dtype)
137
    mu_image = moments(image)
138
    assert mu_image.dtype == expected_dtype
139

140
    coords = np.array(
141
        [[r, c] for r in range(13, 17) for c in range(13, 17)], dtype=dtype
142
    )
143
    mu_coords = moments_coords(coords)
144
    assert mu_coords.dtype == expected_dtype
145

146
    assert_almost_equal(mu_coords, mu_image)
147

148

149
def test_moments_central_coords():
150
    image = np.zeros((20, 20), dtype=np.float64)
151
    image[13:17, 13:17] = 1
152
    mu_image = moments_central(image, (14.5, 14.5))
153

154
    coords = np.array(
155
        [[r, c] for r in range(13, 17) for c in range(13, 17)], dtype=np.float64
156
    )
157
    mu_coords = moments_coords_central(coords, (14.5, 14.5))
158
    assert_almost_equal(mu_coords, mu_image)
159

160
    # ensure that center is being calculated normally
161
    mu_coords_calc_centroid = moments_coords_central(coords)
162
    assert_almost_equal(mu_coords_calc_centroid, mu_coords)
163

164
    # shift image by dx=3 dy=3
165
    image = np.zeros((20, 20), dtype=np.float64)
166
    image[16:20, 16:20] = 1
167
    mu_image = moments_central(image, (14.5, 14.5))
168

169
    coords = np.array(
170
        [[r, c] for r in range(16, 20) for c in range(16, 20)], dtype=np.float64
171
    )
172
    mu_coords = moments_coords_central(coords, (14.5, 14.5))
173
    assert_almost_equal(mu_coords, mu_image)
174

175

176
def test_moments_normalized():
177
    image = np.zeros((20, 20), dtype=np.float64)
178
    image[13:17, 13:17] = 1
179
    mu = moments_central(image, (14.5, 14.5))
180
    nu = moments_normalized(mu)
181
    # shift image by dx=-2, dy=-2 and scale non-zero extent by 0.5
182
    image2 = np.zeros((20, 20), dtype=np.float64)
183
    # scale amplitude by 0.7
184
    image2[11:13, 11:13] = 0.7
185
    mu2 = moments_central(image2, (11.5, 11.5))
186
    nu2 = moments_normalized(mu2)
187
    # central moments must be translation and scale invariant
188
    assert_almost_equal(nu, nu2, decimal=1)
189

190

191
@pytest.mark.parametrize('anisotropic', [False, True])
192
def test_moments_normalized_spacing(anisotropic):
193
    image = np.zeros((20, 20), dtype=np.double)
194
    image[13:17, 13:17] = 1
195

196
    if not anisotropic:
197
        spacing1 = (1, 1)
198
        spacing2 = (3, 3)
199
    else:
200
        spacing1 = (1, 2)
201
        spacing2 = (2, 4)
202

203
    mu = moments_central(image, spacing=spacing1)
204
    nu = moments_normalized(mu, spacing=spacing1)
205

206
    mu2 = moments_central(image, spacing=spacing2)
207
    nu2 = moments_normalized(mu2, spacing=spacing2)
208

209
    # result should be invariant to absolute scale of spacing
210
    compare_moments(nu, nu2)
211

212

213
def test_moments_normalized_3d():
214
    image = draw.ellipsoid(1, 1, 10)
215
    mu_image = moments_central(image)
216
    nu = moments_normalized(mu_image)
217
    assert nu[0, 0, 2] > nu[0, 2, 0]
218
    assert_almost_equal(nu[0, 2, 0], nu[2, 0, 0])
219

220
    coords = np.where(image)
221
    mu_coords = moments_coords_central(coords)
222
    assert_almost_equal(mu_image, mu_coords)
223

224

225
@pytest.mark.parametrize('dtype', [np.uint8, np.int32, np.float32, np.float64])
226
@pytest.mark.parametrize('order', [1, 2, 3, 4])
227
@pytest.mark.parametrize('ndim', [2, 3, 4])
228
def test_analytical_moments_calculation(dtype, order, ndim):
229
    if ndim == 2:
230
        shape = (256, 256)
231
    elif ndim == 3:
232
        shape = (64, 64, 64)
233
    else:
234
        shape = (16,) * ndim
235
    rng = np.random.default_rng(1234)
236
    if np.dtype(dtype).kind in 'iu':
237
        x = rng.integers(0, 256, shape, dtype=dtype)
238
    else:
239
        x = rng.standard_normal(shape, dtype=dtype)
240
    # setting center=None will use the analytical expressions
241
    m1 = moments_central(x, center=None, order=order)
242
    # providing explicit centroid will bypass the analytical code path
243
    m2 = moments_central(x, center=centroid(x), order=order)
244
    # ensure numeric and analytical central moments are close
245
    # TODO: np 2 failed w/ thresh = 1e-4
246
    thresh = 1.5e-4 if x.dtype == np.float32 else 1e-9
247
    compare_moments(m1, m2, thresh=thresh)
248

249

250
def test_moments_normalized_invalid():
251
    with testing.raises(ValueError):
252
        moments_normalized(np.zeros((3, 3)), 3)
253
    with testing.raises(ValueError):
254
        moments_normalized(np.zeros((3, 3)), 4)
255

256

257
def test_moments_hu():
258
    image = np.zeros((20, 20), dtype=np.float64)
259
    image[13:15, 13:17] = 1
260
    mu = moments_central(image, (13.5, 14.5))
261
    nu = moments_normalized(mu)
262
    hu = moments_hu(nu)
263
    # shift image by dx=2, dy=3, scale by 0.5 and rotate by 90deg
264
    image2 = np.zeros((20, 20), dtype=np.float64)
265
    image2[11, 11:13] = 1
266
    image2 = image2.T
267
    mu2 = moments_central(image2, (11.5, 11))
268
    nu2 = moments_normalized(mu2)
269
    hu2 = moments_hu(nu2)
270
    # central moments must be translation and scale invariant
271
    assert_almost_equal(hu, hu2, decimal=1)
272

273

274
@pytest.mark.parametrize('dtype', [np.float16, np.float32, np.float64])
275
def test_moments_dtype(dtype):
276
    image = np.zeros((20, 20), dtype=dtype)
277
    image[13:15, 13:17] = 1
278

279
    expected_dtype = _supported_float_type(dtype)
280
    mu = moments_central(image, (13.5, 14.5))
281
    assert mu.dtype == expected_dtype
282

283
    nu = moments_normalized(mu)
284
    assert nu.dtype == expected_dtype
285

286
    hu = moments_hu(nu)
287
    assert hu.dtype == expected_dtype
288

289

290
@pytest.mark.parametrize('dtype', [np.float16, np.float32, np.float64])
291
def test_centroid(dtype):
292
    image = np.zeros((20, 20), dtype=dtype)
293
    image[14, 14:16] = 1
294
    image[15, 14:16] = 1 / 3
295
    image_centroid = centroid(image)
296
    if dtype == np.float16:
297
        rtol = 1e-3
298
    elif dtype == np.float32:
299
        rtol = 1e-5
300
    else:
301
        rtol = 1e-7
302
    assert_allclose(image_centroid, (14.25, 14.5), rtol=rtol)
303

304

305
@pytest.mark.parametrize('dtype', [np.float16, np.float32, np.float64])
306
def test_inertia_tensor_2d(dtype):
307
    image = np.zeros((40, 40), dtype=dtype)
308
    image[15:25, 5:35] = 1  # big horizontal rectangle (aligned with axis 1)
309
    expected_dtype = _supported_float_type(image.dtype)
310

311
    T = inertia_tensor(image)
312
    assert T.dtype == expected_dtype
313
    assert T[0, 0] > T[1, 1]
314
    np.testing.assert_allclose(T[0, 1], 0)
315

316
    v0, v1 = inertia_tensor_eigvals(image, T=T)
317
    assert v0.dtype == expected_dtype
318
    assert v1.dtype == expected_dtype
319
    np.testing.assert_allclose(np.sqrt(v0 / v1), 3, rtol=0.01, atol=0.05)
320

321

322
def test_inertia_tensor_3d():
323
    image = draw.ellipsoid(10, 5, 3)
324
    T0 = inertia_tensor(image)
325
    eig0, V0 = np.linalg.eig(T0)
326
    # principal axis of ellipse = eigenvector of smallest eigenvalue
327
    v0 = V0[:, np.argmin(eig0)]
328

329
    assert np.allclose(v0, [1, 0, 0]) or np.allclose(-v0, [1, 0, 0])
330

331
    imrot = ndi.rotate(image.astype(float), 30, axes=(0, 1), order=1)
332
    Tr = inertia_tensor(imrot)
333
    eigr, Vr = np.linalg.eig(Tr)
334
    vr = Vr[:, np.argmin(eigr)]
335

336
    # Check that axis has rotated by expected amount
337
    pi, cos, sin = np.pi, np.cos, np.sin
338
    R = np.array(
339
        [[cos(pi / 6), -sin(pi / 6), 0], [sin(pi / 6), cos(pi / 6), 0], [0, 0, 1]]
340
    )
341
    expected_vr = R @ v0
342
    assert np.allclose(vr, expected_vr, atol=1e-3, rtol=0.01) or np.allclose(
343
        -vr, expected_vr, atol=1e-3, rtol=0.01
344
    )
345

346

347
def test_inertia_tensor_eigvals():
348
    # Floating point precision problems could make a positive
349
    # semidefinite matrix have an eigenvalue that is very slightly
350
    # negative.  Check that we have caught and fixed this problem.
351
    image = np.array(
352
        [
353
            [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
354
            [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
355
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
356
        ]
357
    )
358
    # mu = np.array([[3, 0, 98], [0, 14, 0], [2, 0, 98]])
359
    eigvals = inertia_tensor_eigvals(image=image)
360
    assert min(eigvals) >= 0
361