scikit-image

1
import math
2

3
import re
4
import numpy as np
5
import pytest
6
import scipy.ndimage as ndi
7
from numpy.testing import (
8
    assert_allclose,
9
    assert_almost_equal,
10
    assert_array_almost_equal,
11
    assert_array_equal,
12
    assert_equal,
13
)
14

15
from skimage import data, draw, transform
16
from skimage._shared import testing
17
from skimage.measure._regionprops import (
18
    COL_DTYPES,
19
    OBJECT_COLUMNS,
20
    PROPS,
21
    _inertia_eigvals_to_axes_lengths_3D,
22
    _parse_docs,
23
    _props_to_dict,
24
    _require_intensity_image,
25
    euler_number,
26
    perimeter,
27
    perimeter_crofton,
28
    regionprops,
29
    regionprops_table,
30
)
31
from skimage.segmentation import slic
32

33
SAMPLE = np.array(
34
    [
35
        [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0],
36
        [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
37
        [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0],
38
        [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0],
39
        [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
40
        [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
41
        [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
42
        [1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0],
43
        [0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1],
44
        [0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
45
    ]
46
)
47
INTENSITY_SAMPLE = SAMPLE.copy()
48
INTENSITY_SAMPLE[1, 9:11] = 2
49
INTENSITY_FLOAT_SAMPLE = INTENSITY_SAMPLE.copy().astype(np.float64) / 10.0
50
INTENSITY_FLOAT_SAMPLE_MULTICHANNEL = INTENSITY_FLOAT_SAMPLE[..., np.newaxis] * [
51
    1,
52
    2,
53
    3,
54
]
55

56
SAMPLE_MULTIPLE = np.eye(10, dtype=np.int32)
57
SAMPLE_MULTIPLE[3:5, 7:8] = 2
58
INTENSITY_SAMPLE_MULTIPLE = SAMPLE_MULTIPLE.copy() * 2.0
59

60
SAMPLE_3D = np.zeros((6, 6, 6), dtype=np.uint8)
61
SAMPLE_3D[1:3, 1:3, 1:3] = 1
62
SAMPLE_3D[3, 2, 2] = 1
63
INTENSITY_SAMPLE_3D = SAMPLE_3D.copy()
64

65

66
def get_moment_function(img, spacing=(1, 1)):
67
    rows, cols = img.shape
68
    Y, X = np.meshgrid(
69
        np.linspace(0, rows * spacing[0], rows, endpoint=False),
70
        np.linspace(0, cols * spacing[1], cols, endpoint=False),
71
        indexing='ij',
72
    )
73
    return lambda p, q: np.sum(Y**p * X**q * img)
74

75

76
def get_moment3D_function(img, spacing=(1, 1, 1)):
77
    slices, rows, cols = img.shape
78
    Z, Y, X = np.meshgrid(
79
        np.linspace(0, slices * spacing[0], slices, endpoint=False),
80
        np.linspace(0, rows * spacing[1], rows, endpoint=False),
81
        np.linspace(0, cols * spacing[2], cols, endpoint=False),
82
        indexing='ij',
83
    )
84
    return lambda p, q, r: np.sum(Z**p * Y**q * X**r * img)
85

86

87
def get_central_moment_function(img, spacing=(1, 1)):
88
    rows, cols = img.shape
89
    Y, X = np.meshgrid(
90
        np.linspace(0, rows * spacing[0], rows, endpoint=False),
91
        np.linspace(0, cols * spacing[1], cols, endpoint=False),
92
        indexing='ij',
93
    )
94

95
    Mpq = get_moment_function(img, spacing=spacing)
96
    cY = Mpq(1, 0) / Mpq(0, 0)
97
    cX = Mpq(0, 1) / Mpq(0, 0)
98

99
    return lambda p, q: np.sum((Y - cY) ** p * (X - cX) ** q * img)
100

101

102
def test_all_props():
103
    region = regionprops(SAMPLE, INTENSITY_SAMPLE)[0]
104
    for prop in PROPS:
105
        try:
106
            # access legacy name via dict
107
            assert_almost_equal(region[prop], getattr(region, PROPS[prop]))
108

109
            # skip property access tests for old CamelCase names
110
            # (we intentionally do not provide properties for these)
111
            if prop.lower() == prop:
112
                # access legacy name via attribute
113
                assert_almost_equal(getattr(region, prop), getattr(region, PROPS[prop]))
114

115
        except TypeError:  # the `slice` property causes this
116
            pass
117

118

119
def test_all_props_3d():
120
    region = regionprops(SAMPLE_3D, INTENSITY_SAMPLE_3D)[0]
121
    for prop in PROPS:
122
        try:
123
            assert_almost_equal(region[prop], getattr(region, PROPS[prop]))
124

125
            # skip property access tests for old CamelCase names
126
            # (we intentionally do not provide properties for these)
127
            if prop.lower() == prop:
128
                assert_almost_equal(getattr(region, prop), getattr(region, PROPS[prop]))
129

130
        except (NotImplementedError, TypeError):
131
            pass
132

133

134
def test_num_pixels():
135
    num_pixels = regionprops(SAMPLE)[0].num_pixels
136
    assert num_pixels == 72
137

138
    num_pixels = regionprops(SAMPLE, spacing=(2, 1))[0].num_pixels
139
    assert num_pixels == 72
140

141

142
def test_dtype():
143
    regionprops(np.zeros((10, 10), dtype=int))
144
    regionprops(np.zeros((10, 10), dtype=np.uint))
145
    with pytest.raises(TypeError):
146
        regionprops(np.zeros((10, 10), dtype=float))
147
    with pytest.raises(TypeError):
148
        regionprops(np.zeros((10, 10), dtype=np.float64))
149
    with pytest.raises(TypeError):
150
        regionprops(np.zeros((10, 10), dtype=bool))
151

152

153
def test_ndim():
154
    regionprops(np.zeros((10, 10), dtype=int))
155
    regionprops(np.zeros((10, 10, 1), dtype=int))
156
    regionprops(np.zeros((10, 10, 10), dtype=int))
157
    regionprops(np.zeros((1, 1), dtype=int))
158
    regionprops(np.zeros((1, 1, 1), dtype=int))
159
    with pytest.raises(TypeError):
160
        regionprops(np.zeros((10, 10, 10, 2), dtype=int))
161

162

163
def test_feret_diameter_max():
164
    # comparator result is based on SAMPLE from manually-inspected computations
165
    comparator_result = 18
166
    test_result = regionprops(SAMPLE)[0].feret_diameter_max
167
    assert np.abs(test_result - comparator_result) < 1
168
    comparator_result_spacing = 10
169
    test_result_spacing = regionprops(SAMPLE, spacing=[1, 0.1])[0].feret_diameter_max
170
    assert np.abs(test_result_spacing - comparator_result_spacing) < 1
171
    # square, test that maximum Feret diameter is sqrt(2) * square side
172
    img = np.zeros((20, 20), dtype=np.uint8)
173
    img[2:-2, 2:-2] = 1
174
    feret_diameter_max = regionprops(img)[0].feret_diameter_max
175
    assert np.abs(feret_diameter_max - 16 * np.sqrt(2)) < 1
176
    # Due to marching-squares with a level of .5 the diagonal goes from (0, 0.5) to (16, 15.5).
177
    assert np.abs(feret_diameter_max - np.sqrt(16**2 + (16 - 1) ** 2)) < 1e-6
178

179

180
def test_feret_diameter_max_spacing():
181
    # comparator result is based on SAMPLE from manually-inspected computations
182
    comparator_result = 18
183
    test_result = regionprops(SAMPLE)[0].feret_diameter_max
184
    assert np.abs(test_result - comparator_result) < 1
185
    spacing = (2, 1)
186
    # square, test that maximum Feret diameter is sqrt(2) * square side
187
    img = np.zeros((20, 20), dtype=np.uint8)
188
    img[2:-2, 2:-2] = 1
189
    feret_diameter_max = regionprops(img, spacing=spacing)[0].feret_diameter_max
190
    # For anisotropic spacing the shift is applied to the smaller spacing.
191
    assert (
192
        np.abs(
193
            feret_diameter_max
194
            - np.sqrt(
195
                (spacing[0] * 16 - (spacing[0] <= spacing[1])) ** 2
196
                + (spacing[1] * 16 - (spacing[1] < spacing[0])) ** 2
197
            )
198
        )
199
        < 1e-6
200
    )
201

202

203
def test_feret_diameter_max_3d():
204
    img = np.zeros((20, 20), dtype=np.uint8)
205
    img[2:-2, 2:-2] = 1
206
    img_3d = np.dstack((img,) * 3)
207
    feret_diameter_max = regionprops(img_3d)[0].feret_diameter_max
208
    # Due to marching-cubes with a level of .5 -1=2*0.5 has to be subtracted from two axes.
209
    # There are three combinations (x-1, y-1, z), (x-1, y, z-1), (x, y-1, z-1). The option
210
    # yielding the longest diagonal is the computed max_feret_diameter.
211
    assert (
212
        np.abs(feret_diameter_max - np.sqrt((16 - 1) ** 2 + 16**2 + (3 - 1) ** 2))
213
        < 1e-6
214
    )
215
    spacing = (1, 2, 3)
216
    feret_diameter_max = regionprops(img_3d, spacing=spacing)[0].feret_diameter_max
217
    # The longest of the three options is the max_feret_diameter
218
    assert (
219
        np.abs(
220
            feret_diameter_max
221
            - np.sqrt(
222
                (spacing[0] * (16 - 1)) ** 2
223
                + (spacing[1] * (16 - 0)) ** 2
224
                + (spacing[2] * (3 - 1)) ** 2
225
            )
226
        )
227
        < 1e-6
228
    )
229
    assert (
230
        np.abs(
231
            feret_diameter_max
232
            - np.sqrt(
233
                (spacing[0] * (16 - 1)) ** 2
234
                + (spacing[1] * (16 - 1)) ** 2
235
                + (spacing[2] * (3 - 0)) ** 2
236
            )
237
        )
238
        > 1e-6
239
    )
240
    assert (
241
        np.abs(
242
            feret_diameter_max
243
            - np.sqrt(
244
                (spacing[0] * (16 - 0)) ** 2
245
                + (spacing[1] * (16 - 1)) ** 2
246
                + (spacing[2] * (3 - 1)) ** 2
247
            )
248
        )
249
        > 1e-6
250
    )
251

252

253
@pytest.mark.parametrize(
254
    "sample,spacing",
255
    [
256
        (SAMPLE, None),
257
        (SAMPLE, 1),
258
        (SAMPLE, (1, 1)),
259
        (SAMPLE, (1, 2)),
260
        (SAMPLE_3D, None),
261
        (SAMPLE_3D, 1),
262
        (SAMPLE_3D, (2, 1, 3)),
263
    ],
264
)
265
def test_area(sample, spacing):
266
    area = regionprops(sample, spacing=spacing)[0].area
267
    desired = np.sum(sample * (np.prod(spacing) if spacing else 1))
268
    assert area == desired
269

270

271
def test_bbox():
272
    bbox = regionprops(SAMPLE)[0].bbox
273
    assert_array_almost_equal(bbox, (0, 0, SAMPLE.shape[0], SAMPLE.shape[1]))
274

275
    bbox = regionprops(SAMPLE, spacing=(1, 2))[0].bbox
276
    assert_array_almost_equal(bbox, (0, 0, SAMPLE.shape[0], SAMPLE.shape[1]))
277

278
    SAMPLE_mod = SAMPLE.copy()
279
    SAMPLE_mod[:, -1] = 0
280
    bbox = regionprops(SAMPLE_mod)[0].bbox
281
    assert_array_almost_equal(bbox, (0, 0, SAMPLE.shape[0], SAMPLE.shape[1] - 1))
282
    bbox = regionprops(SAMPLE_mod, spacing=(3, 2))[0].bbox
283
    assert_array_almost_equal(bbox, (0, 0, SAMPLE.shape[0], SAMPLE.shape[1] - 1))
284

285
    bbox = regionprops(SAMPLE_3D)[0].bbox
286
    assert_array_almost_equal(bbox, (1, 1, 1, 4, 3, 3))
287
    bbox = regionprops(SAMPLE_3D, spacing=(0.5, 2, 7))[0].bbox
288
    assert_array_almost_equal(bbox, (1, 1, 1, 4, 3, 3))
289

290

291
def test_area_bbox():
292
    padded = np.pad(SAMPLE, 5, mode='constant')
293
    bbox_area = regionprops(padded)[0].area_bbox
294
    assert_array_almost_equal(bbox_area, SAMPLE.size)
295

296

297
def test_area_bbox_spacing():
298
    spacing = (0.5, 3)
299
    padded = np.pad(SAMPLE, 5, mode='constant')
300
    bbox_area = regionprops(padded, spacing=spacing)[0].area_bbox
301
    assert_array_almost_equal(bbox_area, SAMPLE.size * np.prod(spacing))
302

303

304
def test_moments_central():
305
    mu = regionprops(SAMPLE)[0].moments_central
306
    # determined with OpenCV
307
    assert_almost_equal(mu[2, 0], 436.00000000000045)
308
    # different from OpenCV results, bug in OpenCV
309
    assert_almost_equal(mu[3, 0], -737.333333333333)
310
    assert_almost_equal(mu[1, 1], -87.33333333333303)
311
    assert_almost_equal(mu[2, 1], -127.5555555555593)
312
    assert_almost_equal(mu[0, 2], 1259.7777777777774)
313
    assert_almost_equal(mu[1, 2], 2000.296296296291)
314
    assert_almost_equal(mu[0, 3], -760.0246913580195)
315

316
    # Verify central moment test functions
317
    centralMpq = get_central_moment_function(SAMPLE, spacing=(1, 1))
318
    assert_almost_equal(centralMpq(2, 0), mu[2, 0])
319
    assert_almost_equal(centralMpq(3, 0), mu[3, 0])
320
    assert_almost_equal(centralMpq(1, 1), mu[1, 1])
321
    assert_almost_equal(centralMpq(2, 1), mu[2, 1])
322
    assert_almost_equal(centralMpq(0, 2), mu[0, 2])
323
    assert_almost_equal(centralMpq(1, 2), mu[1, 2])
324
    assert_almost_equal(centralMpq(0, 3), mu[0, 3])
325

326

327
def test_moments_central_spacing():
328
    # Test spacing against verified central moment test function
329
    spacing = (1.8, 0.8)
330
    centralMpq = get_central_moment_function(SAMPLE, spacing=spacing)
331

332
    mu = regionprops(SAMPLE, spacing=spacing)[0].moments_central
333
    assert_almost_equal(mu[2, 0], centralMpq(2, 0))
334
    assert_almost_equal(mu[3, 0], centralMpq(3, 0))
335
    assert_almost_equal(mu[1, 1], centralMpq(1, 1))
336
    assert_almost_equal(mu[2, 1], centralMpq(2, 1))
337
    assert_almost_equal(mu[0, 2], centralMpq(0, 2))
338
    assert_almost_equal(mu[1, 2], centralMpq(1, 2))
339
    assert_almost_equal(mu[0, 3], centralMpq(0, 3))
340

341

342
def test_centroid():
343
    centroid = regionprops(SAMPLE)[0].centroid
344
    # determined with MATLAB
345
    assert_array_almost_equal(centroid, (5.66666666666666, 9.444444444444444))
346

347
    # Verify test moment function with spacing=(1, 1)
348
    Mpq = get_moment_function(SAMPLE, spacing=(1, 1))
349
    cY = Mpq(1, 0) / Mpq(0, 0)
350
    cX = Mpq(0, 1) / Mpq(0, 0)
351

352
    assert_array_almost_equal((cY, cX), centroid)
353

354

355
def test_centroid_spacing():
356
    spacing = (1.8, 0.8)
357
    # Moment
358
    Mpq = get_moment_function(SAMPLE, spacing=spacing)
359
    cY = Mpq(1, 0) / Mpq(0, 0)
360
    cX = Mpq(0, 1) / Mpq(0, 0)
361

362
    centroid = regionprops(SAMPLE, spacing=spacing)[0].centroid
363
    assert_array_almost_equal(centroid, (cY, cX))
364

365

366
def test_centroid_3d():
367
    centroid = regionprops(SAMPLE_3D)[0].centroid
368
    # determined by mean along axis 1 of SAMPLE_3D.nonzero()
369
    assert_array_almost_equal(centroid, (1.66666667, 1.55555556, 1.55555556))
370

371
    # Verify moment 3D test function
372
    Mpqr = get_moment3D_function(SAMPLE_3D, spacing=(1, 1, 1))
373
    cZ = Mpqr(1, 0, 0) / Mpqr(0, 0, 0)
374
    cY = Mpqr(0, 1, 0) / Mpqr(0, 0, 0)
375
    cX = Mpqr(0, 0, 1) / Mpqr(0, 0, 0)
376

377
    assert_array_almost_equal((cZ, cY, cX), centroid)
378

379

380
@pytest.mark.parametrize(
381
    "spacing",
382
    [[2.1, 2.2, 2.3], [2.0, 2.0, 2.0], [2, 2, 2]],
383
)
384
def test_spacing_parameter_3d(spacing):
385
    """Test the _normalize_spacing code."""
386

387
    # Test centroid3d spacing
388
    Mpqr = get_moment3D_function(SAMPLE_3D, spacing=spacing)
389
    cZ = Mpqr(1, 0, 0) / Mpqr(0, 0, 0)
390
    cY = Mpqr(0, 1, 0) / Mpqr(0, 0, 0)
391
    cX = Mpqr(0, 0, 1) / Mpqr(0, 0, 0)
392
    centroid = regionprops(SAMPLE_3D, spacing=spacing)[0].centroid
393
    assert_array_almost_equal(centroid, (cZ, cY, cX))
394

395

396
@pytest.mark.parametrize(
397
    "spacing",
398
    [(1, 1j), 1 + 0j],
399
)
400
def test_spacing_parameter_complex_input(spacing):
401
    """Test the _normalize_spacing code."""
402
    with pytest.raises(
403
        TypeError, match="Element of spacing isn't float or integer type, got"
404
    ):
405
        regionprops(SAMPLE, spacing=spacing)[0].centroid
406

407

408
@pytest.mark.parametrize(
409
    "spacing",
410
    [np.nan, np.inf, -np.inf],
411
)
412
def test_spacing_parameter_nan_inf(spacing):
413
    """Test the _normalize_spacing code."""
414
    with pytest.raises(ValueError):
415
        regionprops(SAMPLE, spacing=spacing)[0].centroid
416

417

418
@pytest.mark.parametrize("spacing", ([1], [[1, 1]], (1, 1, 1)))
419
def test_spacing_mismtaching_shape(spacing):
420
    with pytest.raises(ValueError, match="spacing isn't a scalar nor a sequence"):
421
        regionprops(SAMPLE, spacing=spacing)[0].centroid
422

423

424
@pytest.mark.parametrize("spacing", [[2.1, 2.2], [2.0, 2.0], [2, 2]])
425
def test_spacing_parameter_2d(spacing):
426
    """Test the _normalize_spacing code."""
427
    # Test weight centroid spacing
428
    Mpq = get_moment_function(INTENSITY_SAMPLE, spacing=spacing)
429
    cY = Mpq(0, 1) / Mpq(0, 0)
430
    cX = Mpq(1, 0) / Mpq(0, 0)
431
    centroid = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE, spacing=spacing)[
432
        0
433
    ].centroid_weighted
434
    assert_almost_equal(centroid, (cX, cY))
435

436

437
@pytest.mark.parametrize(
438
    "spacing",
439
    [["bad input"], ["bad input", 1, 2.1]],
440
)
441
def test_spacing_parameter_2d_bad_input(spacing):
442
    """Test the _normalize_spacing code."""
443
    with pytest.raises(ValueError):
444
        regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE, spacing=spacing)[
445
            0
446
        ].centroid_weighted
447

448

449
def test_area_convex():
450
    area = regionprops(SAMPLE)[0].area_convex
451
    assert area == 125
452

453

454
def test_area_convex_spacing():
455
    spacing = (1, 4)
456
    area = regionprops(SAMPLE, spacing=spacing)[0].area_convex
457
    assert area == 125 * np.prod(spacing)
458

459

460
def test_image_convex():
461
    img = regionprops(SAMPLE)[0].image_convex
462
    ref = np.array(
463
        [
464
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0],
465
            [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
466
            [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
467
            [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
468
            [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
469
            [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
470
            [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
471
            [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
472
            [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
473
            [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
474
        ]
475
    )
476
    assert_array_equal(img, ref)
477

478

479
def test_coordinates():
480
    sample = np.zeros((10, 10), dtype=np.int8)
481
    coords = np.array([[3, 2], [3, 3], [3, 4]])
482
    sample[coords[:, 0], coords[:, 1]] = 1
483
    prop_coords = regionprops(sample)[0].coords
484
    assert_array_equal(prop_coords, coords)
485
    prop_coords = regionprops(sample, spacing=(0.5, 1.2))[0].coords
486
    assert_array_equal(prop_coords, coords)
487

488

489
@pytest.mark.parametrize(
490
    "spacing",
491
    [None, 1, 2, (1, 1), (1, 0.5)],
492
)
493
def test_coordinates_scaled(spacing):
494
    sample = np.zeros((10, 10), dtype=np.int8)
495
    coords = np.array([[3, 2], [3, 3], [3, 4]])
496
    sample[coords[:, 0], coords[:, 1]] = 1
497
    prop_coords = regionprops(sample, spacing=spacing)[0].coords_scaled
498
    if spacing is None:
499
        desired_coords = coords
500
    else:
501
        desired_coords = coords * np.array(spacing)
502
    assert_array_equal(prop_coords, desired_coords)
503

504

505
@pytest.mark.parametrize(
506
    "spacing",
507
    [None, 1, 2, (0.2, 3, 2.3)],
508
)
509
def test_coordinates_scaled_3d(spacing):
510
    sample = np.zeros((6, 6, 6), dtype=np.int8)
511
    coords = np.array([[1, 1, 1], [1, 2, 1], [1, 3, 1]])
512
    sample[coords[:, 0], coords[:, 1], coords[:, 2]] = 1
513
    prop_coords = regionprops(sample, spacing=spacing)[0].coords_scaled
514
    if spacing is None:
515
        desired_coords = coords
516
    else:
517
        desired_coords = coords * np.array(spacing)
518
    assert_array_equal(prop_coords, desired_coords)
519

520

521
def test_slice():
522
    padded = np.pad(SAMPLE, ((2, 4), (5, 2)), mode='constant')
523
    nrow, ncol = SAMPLE.shape
524
    result = regionprops(padded)[0].slice
525
    expected = (slice(2, 2 + nrow), slice(5, 5 + ncol))
526
    assert_equal(result, expected)
527

528

529
def test_slice_spacing():
530
    padded = np.pad(SAMPLE, ((2, 4), (5, 2)), mode='constant')
531
    nrow, ncol = SAMPLE.shape
532
    result = regionprops(padded)[0].slice
533
    expected = (slice(2, 2 + nrow), slice(5, 5 + ncol))
534
    spacing = (2, 0.2)
535
    result = regionprops(padded, spacing=spacing)[0].slice
536
    assert_equal(result, expected)
537

538

539
def test_eccentricity():
540
    eps = regionprops(SAMPLE)[0].eccentricity
541
    assert_almost_equal(eps, 0.814629313427)
542

543
    eps = regionprops(SAMPLE, spacing=(1.5, 1.5))[0].eccentricity
544
    assert_almost_equal(eps, 0.814629313427)
545

546
    img = np.zeros((5, 5), dtype=int)
547
    img[2, 2] = 1
548
    eps = regionprops(img)[0].eccentricity
549
    assert_almost_equal(eps, 0)
550

551
    eps = regionprops(img, spacing=(3, 3))[0].eccentricity
552
    assert_almost_equal(eps, 0)
553

554

555
def test_equivalent_diameter_area():
556
    diameter = regionprops(SAMPLE)[0].equivalent_diameter_area
557
    # determined with MATLAB
558
    assert_almost_equal(diameter, 9.57461472963)
559

560
    spacing = (1, 3)
561
    diameter = regionprops(SAMPLE, spacing=spacing)[0].equivalent_diameter_area
562
    equivalent_area = np.pi * (diameter / 2.0) ** 2
563
    assert_almost_equal(equivalent_area, SAMPLE.sum() * np.prod(spacing))
564

565

566
def test_euler_number():
567
    for spacing in [(1, 1), (2.1, 0.9)]:
568
        en = regionprops(SAMPLE, spacing=spacing)[0].euler_number
569
        assert en == 0
570

571
        SAMPLE_mod = SAMPLE.copy()
572
        SAMPLE_mod[7, -3] = 0
573
        en = regionprops(SAMPLE_mod, spacing=spacing)[0].euler_number
574
        assert en == -1
575

576
        en = euler_number(SAMPLE, 1)
577
        assert en == 2
578

579
        en = euler_number(SAMPLE_mod, 1)
580
        assert en == 1
581

582
    en = euler_number(SAMPLE_3D, 1)
583
    assert en == 1
584

585
    en = euler_number(SAMPLE_3D, 3)
586
    assert en == 1
587

588
    # for convex body, Euler number is 1
589
    SAMPLE_3D_2 = np.zeros((100, 100, 100))
590
    SAMPLE_3D_2[40:60, 40:60, 40:60] = 1
591
    en = euler_number(SAMPLE_3D_2, 3)
592
    assert en == 1
593

594
    SAMPLE_3D_2[45:55, 45:55, 45:55] = 0
595
    en = euler_number(SAMPLE_3D_2, 3)
596
    assert en == 2
597

598

599
def test_extent():
600
    extent = regionprops(SAMPLE)[0].extent
601
    assert_almost_equal(extent, 0.4)
602
    extent = regionprops(SAMPLE, spacing=(5, 0.2))[0].extent
603
    assert_almost_equal(extent, 0.4)
604

605

606
def test_moments_hu():
607
    hu = regionprops(SAMPLE)[0].moments_hu
608
    ref = np.array(
609
        [
610
            3.27117627e-01,
611
            2.63869194e-02,
612
            2.35390060e-02,
613
            1.23151193e-03,
614
            1.38882330e-06,
615
            -2.72586158e-05,
616
            -6.48350653e-06,
617
        ]
618
    )
619
    # bug in OpenCV caused in Central Moments calculation?
620
    assert_array_almost_equal(hu, ref)
621

622
    with testing.raises(NotImplementedError):
623
        regionprops(SAMPLE, spacing=(2, 1))[0].moments_hu
624

625

626
def test_image():
627
    img = regionprops(SAMPLE)[0].image
628
    assert_array_equal(img, SAMPLE)
629

630
    img = regionprops(SAMPLE_3D)[0].image
631
    assert_array_equal(img, SAMPLE_3D[1:4, 1:3, 1:3])
632

633

634
def test_label():
635
    label = regionprops(SAMPLE)[0].label
636
    assert_array_equal(label, 1)
637

638
    label = regionprops(SAMPLE_3D)[0].label
639
    assert_array_equal(label, 1)
640

641

642
def test_area_filled():
643
    area = regionprops(SAMPLE)[0].area_filled
644
    assert area == np.sum(SAMPLE)
645

646

647
def test_area_filled_zero():
648
    SAMPLE_mod = SAMPLE.copy()
649
    SAMPLE_mod[7, -3] = 0
650
    area = regionprops(SAMPLE_mod)[0].area_filled
651
    assert area == np.sum(SAMPLE)
652

653

654
def test_area_filled_spacing():
655
    SAMPLE_mod = SAMPLE.copy()
656
    SAMPLE_mod[7, -3] = 0
657

658
    spacing = (2, 1.2)
659
    area = regionprops(SAMPLE, spacing=spacing)[0].area_filled
660
    assert area == np.sum(SAMPLE) * np.prod(spacing)
661

662
    area = regionprops(SAMPLE_mod, spacing=spacing)[0].area_filled
663
    assert area == np.sum(SAMPLE) * np.prod(spacing)
664

665

666
def test_image_filled():
667
    img = regionprops(SAMPLE)[0].image_filled
668
    assert_array_equal(img, SAMPLE)
669
    img = regionprops(SAMPLE, spacing=(1, 4))[0].image_filled
670
    assert_array_equal(img, SAMPLE)
671

672

673
def test_axis_major_length():
674
    length = regionprops(SAMPLE)[0].axis_major_length
675
    # MATLAB has different interpretation of ellipse than found in literature,
676
    # here implemented as found in literature
677
    target_length = 16.7924234999
678
    assert_almost_equal(length, target_length)
679

680
    length = regionprops(SAMPLE, spacing=(2, 2))[0].axis_major_length
681
    assert_almost_equal(length, 2 * target_length)
682

683
    from skimage.draw import ellipse
684

685
    img = np.zeros((20, 24), dtype=np.uint8)
686
    rr, cc = ellipse(11, 11, 7, 9, rotation=np.deg2rad(45))
687
    img[rr, cc] = 1
688

689
    target_length = regionprops(img, spacing=(1, 1))[0].axis_major_length
690
    length_wo_spacing = regionprops(img[::2], spacing=(1, 1))[0].axis_minor_length
691
    assert abs(length_wo_spacing - target_length) > 0.1
692
    length = regionprops(img[:, ::2], spacing=(1, 2))[0].axis_major_length
693
    assert_almost_equal(length, target_length, decimal=0)
694

695

696
def test_intensity_max():
697
    intensity = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE)[0].intensity_max
698
    assert_almost_equal(intensity, 2)
699

700

701
def test_intensity_mean():
702
    intensity = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE)[0].intensity_mean
703
    assert_almost_equal(intensity, 1.02777777777777)
704

705

706
def test_intensity_min():
707
    intensity = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE)[0].intensity_min
708
    assert_almost_equal(intensity, 1)
709

710

711
def test_intensity_std():
712
    intensity = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE)[0].intensity_std
713
    assert_almost_equal(intensity, 0.16433554953054486)
714

715

716
def test_axis_minor_length():
717
    length = regionprops(SAMPLE)[0].axis_minor_length
718
    # MATLAB has different interpretation of ellipse than found in literature,
719
    # here implemented as found in literature
720
    target_length = 9.739302807263
721
    assert_almost_equal(length, target_length)
722

723
    length = regionprops(SAMPLE, spacing=(1.5, 1.5))[0].axis_minor_length
724
    assert_almost_equal(length, 1.5 * target_length)
725

726
    from skimage.draw import ellipse
727

728
    img = np.zeros((10, 12), dtype=np.uint8)
729
    rr, cc = ellipse(5, 6, 3, 5, rotation=np.deg2rad(30))
730
    img[rr, cc] = 1
731

732
    target_length = regionprops(img, spacing=(1, 1))[0].axis_minor_length
733
    length_wo_spacing = regionprops(img[::2], spacing=(1, 1))[0].axis_minor_length
734
    assert abs(length_wo_spacing - target_length) > 0.1
735
    length = regionprops(img[::2], spacing=(2, 1))[0].axis_minor_length
736
    assert_almost_equal(length, target_length, decimal=1)
737

738

739
def test_moments():
740
    m = regionprops(SAMPLE)[0].moments
741
    # determined with OpenCV
742
    assert_almost_equal(m[0, 0], 72.0)
743
    assert_almost_equal(m[0, 1], 680.0)
744
    assert_almost_equal(m[0, 2], 7682.0)
745
    assert_almost_equal(m[0, 3], 95588.0)
746
    assert_almost_equal(m[1, 0], 408.0)
747
    assert_almost_equal(m[1, 1], 3766.0)
748
    assert_almost_equal(m[1, 2], 43882.0)
749
    assert_almost_equal(m[2, 0], 2748.0)
750
    assert_almost_equal(m[2, 1], 24836.0)
751
    assert_almost_equal(m[3, 0], 19776.0)
752

753
    # Verify moment test function
754
    Mpq = get_moment_function(SAMPLE, spacing=(1, 1))
755
    assert_almost_equal(Mpq(0, 0), m[0, 0])
756
    assert_almost_equal(Mpq(0, 1), m[0, 1])
757
    assert_almost_equal(Mpq(0, 2), m[0, 2])
758
    assert_almost_equal(Mpq(0, 3), m[0, 3])
759
    assert_almost_equal(Mpq(1, 0), m[1, 0])
760
    assert_almost_equal(Mpq(1, 1), m[1, 1])
761
    assert_almost_equal(Mpq(1, 2), m[1, 2])
762
    assert_almost_equal(Mpq(2, 0), m[2, 0])
763
    assert_almost_equal(Mpq(2, 1), m[2, 1])
764
    assert_almost_equal(Mpq(3, 0), m[3, 0])
765

766

767
def test_moments_spacing():
768
    # Test moment on spacing
769
    spacing = (2, 0.3)
770
    m = regionprops(SAMPLE, spacing=spacing)[0].moments
771
    Mpq = get_moment_function(SAMPLE, spacing=spacing)
772
    assert_almost_equal(m[0, 0], Mpq(0, 0))
773
    assert_almost_equal(m[0, 1], Mpq(0, 1))
774
    assert_almost_equal(m[0, 2], Mpq(0, 2))
775
    assert_almost_equal(m[0, 3], Mpq(0, 3))
776
    assert_almost_equal(m[1, 0], Mpq(1, 0))
777
    assert_almost_equal(m[1, 1], Mpq(1, 1))
778
    assert_almost_equal(m[1, 2], Mpq(1, 2))
779
    assert_almost_equal(m[2, 0], Mpq(2, 0))
780
    assert_almost_equal(m[2, 1], Mpq(2, 1))
781
    assert_almost_equal(m[3, 0], Mpq(3, 0))
782

783

784
def test_moments_normalized():
785
    nu = regionprops(SAMPLE)[0].moments_normalized
786

787
    # determined with OpenCV
788
    assert_almost_equal(nu[0, 2], 0.24301268861454037)
789
    assert_almost_equal(nu[0, 3], -0.017278118992041805)
790
    assert_almost_equal(nu[1, 1], -0.016846707818929982)
791
    assert_almost_equal(nu[1, 2], 0.045473992910668816)
792
    assert_almost_equal(nu[2, 0], 0.08410493827160502)
793
    assert_almost_equal(nu[2, 1], -0.002899800614433943)
794

795

796
def test_moments_normalized_spacing():
797
    spacing = (3, 3)
798
    nu = regionprops(SAMPLE, spacing=spacing)[0].moments_normalized
799

800
    # Normalized moments are scale invariant.
801
    assert_almost_equal(nu[0, 2], 0.24301268861454037)
802
    assert_almost_equal(nu[0, 3], -0.017278118992041805)
803
    assert_almost_equal(nu[1, 1], -0.016846707818929982)
804
    assert_almost_equal(nu[1, 2], 0.045473992910668816)
805
    assert_almost_equal(nu[2, 0], 0.08410493827160502)
806
    assert_almost_equal(nu[2, 1], -0.002899800614433943)
807

808

809
def test_orientation():
810
    orient = regionprops(SAMPLE)[0].orientation
811
    # determined with MATLAB
812
    target_orient = -1.4663278802756865
813
    assert_almost_equal(orient, target_orient)
814

815
    orient = regionprops(SAMPLE, spacing=(2, 2))[0].orientation
816
    assert_almost_equal(orient, target_orient)
817

818
    # test diagonal regions
819
    diag = np.eye(10, dtype=int)
820
    orient_diag = regionprops(diag)[0].orientation
821
    assert_almost_equal(orient_diag, math.pi / 4)
822
    orient_diag = regionprops(diag, spacing=(1, 2))[0].orientation
823
    assert_almost_equal(orient_diag, np.arccos(0.5 / np.sqrt(1 + 0.5**2)))
824
    orient_diag = regionprops(np.flipud(diag))[0].orientation
825
    assert_almost_equal(orient_diag, -math.pi / 4)
826
    orient_diag = regionprops(np.flipud(diag), spacing=(1, 2))[0].orientation
827
    assert_almost_equal(orient_diag, -np.arccos(0.5 / np.sqrt(1 + 0.5**2)))
828
    orient_diag = regionprops(np.fliplr(diag))[0].orientation
829
    assert_almost_equal(orient_diag, -math.pi / 4)
830
    orient_diag = regionprops(np.fliplr(diag), spacing=(1, 2))[0].orientation
831
    assert_almost_equal(orient_diag, -np.arccos(0.5 / np.sqrt(1 + 0.5**2)))
832
    orient_diag = regionprops(np.fliplr(np.flipud(diag)))[0].orientation
833
    assert_almost_equal(orient_diag, math.pi / 4)
834
    orient_diag = regionprops(np.fliplr(np.flipud(diag)), spacing=(1, 2))[0].orientation
835
    assert_almost_equal(orient_diag, np.arccos(0.5 / np.sqrt(1 + 0.5**2)))
836

837

838
def test_orientation_continuity():
839
    # nearly diagonal array
840
    arr1 = np.array([[0, 0, 1, 1], [0, 0, 1, 0], [0, 1, 0, 0], [1, 0, 0, 0]])
841
    # diagonal array
842
    arr2 = np.array([[0, 0, 0, 2], [0, 0, 2, 0], [0, 2, 0, 0], [2, 0, 0, 0]])
843
    # nearly diagonal array
844
    arr3 = np.array([[0, 0, 0, 3], [0, 0, 3, 3], [0, 3, 0, 0], [3, 0, 0, 0]])
845
    image = np.hstack((arr1, arr2, arr3))
846
    props = regionprops(image)
847
    orientations = [prop.orientation for prop in props]
848
    np.testing.assert_allclose(orientations, orientations[1], rtol=0, atol=0.08)
849
    assert_almost_equal(orientations[0], -0.7144496360953664)
850
    assert_almost_equal(orientations[1], -0.7853981633974483)
851
    assert_almost_equal(orientations[2], -0.8563466906995303)
852

853
    # Test spacing
854
    spacing = (3.2, 1.2)
855
    wmu = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE, spacing=spacing)[
856
        0
857
    ].moments_weighted_central
858
    centralMpq = get_central_moment_function(INTENSITY_SAMPLE, spacing=spacing)
859
    assert_almost_equal(wmu[0, 0], centralMpq(0, 0))
860
    assert_almost_equal(wmu[0, 1], centralMpq(0, 1))
861
    assert_almost_equal(wmu[0, 2], centralMpq(0, 2))
862
    assert_almost_equal(wmu[0, 3], centralMpq(0, 3))
863
    assert_almost_equal(wmu[1, 0], centralMpq(1, 0))
864
    assert_almost_equal(wmu[1, 1], centralMpq(1, 1))
865
    assert_almost_equal(wmu[1, 2], centralMpq(1, 2))
866
    assert_almost_equal(wmu[1, 3], centralMpq(1, 3))
867
    assert_almost_equal(wmu[2, 0], centralMpq(2, 0))
868
    assert_almost_equal(wmu[2, 1], centralMpq(2, 1))
869
    assert_almost_equal(wmu[2, 2], centralMpq(2, 2))
870
    assert_almost_equal(wmu[2, 3], centralMpq(2, 3))
871
    assert_almost_equal(wmu[3, 0], centralMpq(3, 0))
872
    assert_almost_equal(wmu[3, 1], centralMpq(3, 1))
873
    assert_almost_equal(wmu[3, 2], centralMpq(3, 2))
874
    assert_almost_equal(wmu[3, 3], centralMpq(3, 3))
875

876

877
def test_perimeter():
878
    per = regionprops(SAMPLE)[0].perimeter
879
    target_per = 55.2487373415
880
    assert_almost_equal(per, target_per)
881
    per = regionprops(SAMPLE, spacing=(2, 2))[0].perimeter
882
    assert_almost_equal(per, 2 * target_per)
883

884
    per = perimeter(SAMPLE.astype('double'), neighborhood=8)
885
    assert_almost_equal(per, 46.8284271247)
886

887
    with testing.raises(NotImplementedError):
888
        per = regionprops(SAMPLE, spacing=(2, 1))[0].perimeter
889

890

891
def test_perimeter_crofton():
892
    per = regionprops(SAMPLE)[0].perimeter_crofton
893
    target_per_crof = 61.0800637973
894
    assert_almost_equal(per, target_per_crof)
895
    per = regionprops(SAMPLE, spacing=(2, 2))[0].perimeter_crofton
896
    assert_almost_equal(per, 2 * target_per_crof)
897

898
    per = perimeter_crofton(SAMPLE.astype('double'), directions=2)
899
    assert_almost_equal(per, 64.4026493985)
900

901
    with testing.raises(NotImplementedError):
902
        per = regionprops(SAMPLE, spacing=(2, 1))[0].perimeter_crofton
903

904

905
def test_solidity():
906
    solidity = regionprops(SAMPLE)[0].solidity
907
    target_solidity = 0.576
908
    assert_almost_equal(solidity, target_solidity)
909

910
    solidity = regionprops(SAMPLE, spacing=(3, 9))[0].solidity
911
    assert_almost_equal(solidity, target_solidity)
912

913

914
def test_multichannel_centroid_weighted_table():
915
    """Test for https://github.com/scikit-image/scikit-image/issues/6860."""
916
    intensity_image = INTENSITY_FLOAT_SAMPLE_MULTICHANNEL
917
    rp0 = regionprops(SAMPLE, intensity_image=intensity_image[..., 0])[0]
918
    rp1 = regionprops(SAMPLE, intensity_image=intensity_image[..., 0:1])[0]
919
    rpm = regionprops(SAMPLE, intensity_image=intensity_image)[0]
920
    np.testing.assert_almost_equal(
921
        rp0.centroid_weighted, np.squeeze(rp1.centroid_weighted)
922
    )
923
    np.testing.assert_almost_equal(
924
        rp0.centroid_weighted, np.array(rpm.centroid_weighted)[:, 0]
925
    )
926
    assert np.shape(rp0.centroid_weighted) == (SAMPLE.ndim,)
927
    assert np.shape(rp1.centroid_weighted) == (SAMPLE.ndim, 1)
928
    assert np.shape(rpm.centroid_weighted) == (SAMPLE.ndim, intensity_image.shape[-1])
929

930
    table = regionprops_table(
931
        SAMPLE, intensity_image=intensity_image, properties=('centroid_weighted',)
932
    )
933
    # check the number of returned columns is correct
934
    assert len(table) == np.size(rpm.centroid_weighted)
935

936

937
def test_moments_weighted_central():
938
    wmu = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE)[
939
        0
940
    ].moments_weighted_central
941
    ref = np.array(
942
        [
943
            [7.4000000000e01, 3.7303493627e-14, 1.2602837838e03, -7.6561796932e02],
944
            [-2.1316282073e-13, -8.7837837838e01, 2.1571526662e03, -4.2385971907e03],
945
            [4.7837837838e02, -1.4801314828e02, 6.6989799420e03, -9.9501164076e03],
946
            [-7.5943608473e02, -1.2714707125e03, 1.5304076361e04, -3.3156729271e04],
947
        ]
948
    )
949

950
    assert_array_almost_equal(wmu, ref)
951

952
    # Verify test function
953
    centralMpq = get_central_moment_function(INTENSITY_SAMPLE, spacing=(1, 1))
954
    assert_almost_equal(centralMpq(0, 0), ref[0, 0])
955
    assert_almost_equal(centralMpq(0, 1), ref[0, 1])
956
    assert_almost_equal(centralMpq(0, 2), ref[0, 2])
957
    assert_almost_equal(centralMpq(0, 3), ref[0, 3])
958
    assert_almost_equal(centralMpq(1, 0), ref[1, 0])
959
    assert_almost_equal(centralMpq(1, 1), ref[1, 1])
960
    assert_almost_equal(centralMpq(1, 2), ref[1, 2])
961
    assert_almost_equal(centralMpq(1, 3), ref[1, 3])
962
    assert_almost_equal(centralMpq(2, 0), ref[2, 0])
963
    assert_almost_equal(centralMpq(2, 1), ref[2, 1])
964
    assert_almost_equal(centralMpq(2, 2), ref[2, 2])
965
    assert_almost_equal(centralMpq(2, 3), ref[2, 3])
966
    assert_almost_equal(centralMpq(3, 0), ref[3, 0])
967
    assert_almost_equal(centralMpq(3, 1), ref[3, 1])
968
    assert_almost_equal(centralMpq(3, 2), ref[3, 2])
969
    assert_almost_equal(centralMpq(3, 3), ref[3, 3])
970

971

972
def test_centroid_weighted():
973
    sample_for_spacing = np.array(
974
        [
975
            [0, 0, 0, 0, 0, 0],
976
            [0, 0, 0, 0, 0, 0],
977
            [0, 0, 0, 0, 0, 0],
978
            [0, 0, 0, 1, 1, 1],
979
            [0, 0, 0, 1, 1, 1],
980
            [0, 0, 0, 1, 1, 1],
981
        ]
982
    )
983
    target_centroid_wspacing = (4.0, 4.0)
984
    centroid = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE)[
985
        0
986
    ].centroid_weighted
987
    target_centroid = (5.540540540540, 9.445945945945)
988
    assert_array_almost_equal(centroid, target_centroid)
989

990
    # Verify test function
991
    Mpq = get_moment_function(INTENSITY_SAMPLE, spacing=(1, 1))
992
    cY = Mpq(0, 1) / Mpq(0, 0)
993
    cX = Mpq(1, 0) / Mpq(0, 0)
994
    assert_almost_equal((cX, cY), centroid)
995

996
    # Test spacing
997
    spacing = (2, 2)
998
    Mpq = get_moment_function(INTENSITY_SAMPLE, spacing=spacing)
999
    cY = Mpq(0, 1) / Mpq(0, 0)
1000
    cX = Mpq(1, 0) / Mpq(0, 0)
1001
    centroid = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE, spacing=spacing)[
1002
        0
1003
    ].centroid_weighted
1004
    assert_almost_equal(centroid, (cX, cY))
1005
    assert_almost_equal(centroid, 2 * np.array(target_centroid))
1006
    centroid = regionprops(
1007
        sample_for_spacing, intensity_image=sample_for_spacing, spacing=spacing
1008
    )[0].centroid_weighted
1009
    assert_almost_equal(centroid, 2 * np.array(target_centroid_wspacing))
1010

1011
    spacing = (1.3, 0.7)
1012
    Mpq = get_moment_function(INTENSITY_SAMPLE, spacing=spacing)
1013
    cY = Mpq(0, 1) / Mpq(0, 0)
1014
    cX = Mpq(1, 0) / Mpq(0, 0)
1015
    centroid = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE, spacing=spacing)[
1016
        0
1017
    ].centroid_weighted
1018
    assert_almost_equal(centroid, (cX, cY))
1019
    centroid = regionprops(
1020
        sample_for_spacing, intensity_image=sample_for_spacing, spacing=spacing
1021
    )[0].centroid_weighted
1022
    assert_almost_equal(centroid, spacing * np.array(target_centroid_wspacing))
1023

1024

1025
def test_moments_weighted_hu():
1026
    whu = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE)[0].moments_weighted_hu
1027
    ref = np.array(
1028
        [
1029
            3.1750587329e-01,
1030
            2.1417517159e-02,
1031
            2.3609322038e-02,
1032
            1.2565683360e-03,
1033
            8.3014209421e-07,
1034
            -3.5073773473e-05,
1035
            -6.7936409056e-06,
1036
        ]
1037
    )
1038
    assert_array_almost_equal(whu, ref)
1039

1040
    with testing.raises(NotImplementedError):
1041
        regionprops(SAMPLE, spacing=(2, 1))[0].moments_weighted_hu
1042

1043

1044
def test_moments_weighted():
1045
    wm = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE)[0].moments_weighted
1046
    ref = np.array(
1047
        [
1048
            [7.4000000e01, 6.9900000e02, 7.8630000e03, 9.7317000e04],
1049
            [4.1000000e02, 3.7850000e03, 4.4063000e04, 5.7256700e05],
1050
            [2.7500000e03, 2.4855000e04, 2.9347700e05, 3.9007170e06],
1051
            [1.9778000e04, 1.7500100e05, 2.0810510e06, 2.8078871e07],
1052
        ]
1053
    )
1054
    assert_array_almost_equal(wm, ref)
1055

1056
    # Verify test function
1057
    Mpq = get_moment_function(INTENSITY_SAMPLE, spacing=(1, 1))
1058
    assert_almost_equal(Mpq(0, 0), ref[0, 0])
1059
    assert_almost_equal(Mpq(0, 1), ref[0, 1])
1060
    assert_almost_equal(Mpq(0, 2), ref[0, 2])
1061
    assert_almost_equal(Mpq(0, 3), ref[0, 3])
1062
    assert_almost_equal(Mpq(1, 0), ref[1, 0])
1063
    assert_almost_equal(Mpq(1, 1), ref[1, 1])
1064
    assert_almost_equal(Mpq(1, 2), ref[1, 2])
1065
    assert_almost_equal(Mpq(1, 3), ref[1, 3])
1066
    assert_almost_equal(Mpq(2, 0), ref[2, 0])
1067
    assert_almost_equal(Mpq(2, 1), ref[2, 1])
1068
    assert_almost_equal(Mpq(2, 2), ref[2, 2])
1069
    assert_almost_equal(Mpq(2, 3), ref[2, 3])
1070
    assert_almost_equal(Mpq(3, 0), ref[3, 0])
1071
    assert_almost_equal(Mpq(3, 1), ref[3, 1])
1072
    assert_almost_equal(Mpq(3, 2), ref[3, 2])
1073
    assert_almost_equal(Mpq(3, 3), ref[3, 3])
1074

1075

1076
def test_moments_weighted_spacing():
1077
    wm = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE)[0].moments_weighted
1078
    ref = np.array(
1079
        [
1080
            [7.4000000e01, 6.9900000e02, 7.8630000e03, 9.7317000e04],
1081
            [4.1000000e02, 3.7850000e03, 4.4063000e04, 5.7256700e05],
1082
            [2.7500000e03, 2.4855000e04, 2.9347700e05, 3.9007170e06],
1083
            [1.9778000e04, 1.7500100e05, 2.0810510e06, 2.8078871e07],
1084
        ]
1085
    )
1086
    assert_array_almost_equal(wm, ref)
1087

1088
    # Test spacing
1089
    spacing = (3.2, 1.2)
1090
    wmu = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE, spacing=spacing)[
1091
        0
1092
    ].moments_weighted
1093
    Mpq = get_moment_function(INTENSITY_SAMPLE, spacing=spacing)
1094
    assert_almost_equal(wmu[0, 0], Mpq(0, 0))
1095
    assert_almost_equal(wmu[0, 1], Mpq(0, 1))
1096
    assert_almost_equal(wmu[0, 2], Mpq(0, 2))
1097
    assert_almost_equal(wmu[0, 3], Mpq(0, 3))
1098
    assert_almost_equal(wmu[1, 0], Mpq(1, 0))
1099
    assert_almost_equal(wmu[1, 1], Mpq(1, 1))
1100
    assert_almost_equal(wmu[1, 2], Mpq(1, 2))
1101
    assert_almost_equal(wmu[1, 3], Mpq(1, 3))
1102
    assert_almost_equal(wmu[2, 0], Mpq(2, 0))
1103
    assert_almost_equal(wmu[2, 1], Mpq(2, 1))
1104
    assert_almost_equal(wmu[2, 2], Mpq(2, 2))
1105
    assert_almost_equal(wmu[2, 3], Mpq(2, 3))
1106
    assert_almost_equal(wmu[3, 0], Mpq(3, 0))
1107
    assert_almost_equal(wmu[3, 1], Mpq(3, 1))
1108
    assert_almost_equal(wmu[3, 2], Mpq(3, 2))
1109
    assert_almost_equal(wmu[3, 3], Mpq(3, 3), decimal=6)
1110

1111

1112
def test_moments_weighted_normalized():
1113
    wnu = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE)[
1114
        0
1115
    ].moments_weighted_normalized
1116
    ref = np.array(
1117
        [
1118
            [np.nan, np.nan, 0.2301467830, -0.0162529732],
1119
            [np.nan, -0.0160405109, 0.0457932622, -0.0104598869],
1120
            [0.0873590903, -0.0031421072, 0.0165315478, -0.0028544152],
1121
            [-0.0161217406, -0.0031376984, 0.0043903193, -0.0011057191],
1122
        ]
1123
    )
1124
    assert_array_almost_equal(wnu, ref)
1125

1126

1127
def test_moments_weighted_normalized_spacing():
1128
    spacing = (3, 3)
1129
    wnu = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE, spacing=spacing)[
1130
        0
1131
    ].moments_weighted_normalized
1132

1133
    np.array(
1134
        [
1135
            [np.nan, np.nan, 0.2301467830, -0.0162529732],
1136
            [np.nan, -0.0160405109, 0.0457932622, -0.0104598869],
1137
            [0.0873590903, -0.0031421072, 0.0165315478, -0.0028544152],
1138
            [-0.0161217406, -0.0031376984, 0.0043903193, -0.0011057191],
1139
        ]
1140
    )
1141

1142
    # Normalized moments are scale invariant
1143
    assert_almost_equal(wnu[0, 2], 0.2301467830)
1144
    assert_almost_equal(wnu[0, 3], -0.0162529732)
1145
    assert_almost_equal(wnu[1, 1], -0.0160405109)
1146
    assert_almost_equal(wnu[1, 2], 0.0457932622)
1147
    assert_almost_equal(wnu[1, 3], -0.0104598869)
1148
    assert_almost_equal(wnu[2, 0], 0.0873590903)
1149
    assert_almost_equal(wnu[2, 1], -0.0031421072)
1150
    assert_almost_equal(wnu[2, 2], 0.0165315478)
1151
    assert_almost_equal(wnu[2, 3], -0.0028544152)
1152
    assert_almost_equal(wnu[3, 0], -0.0161217406)
1153
    assert_almost_equal(wnu[3, 1], -0.0031376984)
1154
    assert_almost_equal(wnu[3, 2], 0.0043903193)
1155
    assert_almost_equal(wnu[3, 3], -0.0011057191)
1156

1157

1158
def test_offset_features():
1159
    props = regionprops(SAMPLE)[0]
1160
    offset = np.array([1024, 2048])
1161
    props_offset = regionprops(SAMPLE, offset=offset)[0]
1162

1163
    assert_allclose(props.centroid, props_offset.centroid - offset)
1164

1165

1166
def test_label_sequence():
1167
    a = np.empty((2, 2), dtype=int)
1168
    a[:, :] = 2
1169
    ps = regionprops(a)
1170
    assert len(ps) == 1
1171
    assert ps[0].label == 2
1172

1173

1174
def test_pure_background():
1175
    a = np.zeros((2, 2), dtype=int)
1176
    ps = regionprops(a)
1177
    assert len(ps) == 0
1178

1179

1180
def test_invalid():
1181
    ps = regionprops(SAMPLE)
1182

1183
    def get_intensity_image():
1184
        ps[0].image_intensity
1185

1186
    with pytest.raises(AttributeError):
1187
        get_intensity_image()
1188

1189

1190
def test_invalid_size():
1191
    wrong_intensity_sample = np.array([[1], [1]])
1192
    with pytest.raises(ValueError):
1193
        regionprops(SAMPLE, wrong_intensity_sample)
1194

1195

1196
def test_equals():
1197
    arr = np.zeros((100, 100), dtype=int)
1198
    arr[0:25, 0:25] = 1
1199
    arr[50:99, 50:99] = 2
1200

1201
    regions = regionprops(arr)
1202
    r1 = regions[0]
1203

1204
    regions = regionprops(arr)
1205
    r2 = regions[0]
1206
    r3 = regions[1]
1207

1208
    assert_equal(r1 == r2, True, "Same regionprops are not equal")
1209
    assert_equal(r1 != r3, True, "Different regionprops are equal")
1210

1211

1212
def test_iterate_all_props():
1213
    region = regionprops(SAMPLE)[0]
1214
    p0 = {p: region[p] for p in region}
1215

1216
    region = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE)[0]
1217
    p1 = {p: region[p] for p in region}
1218

1219
    assert len(p0) < len(p1)
1220

1221

1222
def test_cache():
1223
    SAMPLE_mod = SAMPLE.copy()
1224
    region = regionprops(SAMPLE_mod)[0]
1225
    f0 = region.image_filled
1226
    region._label_image[:10] = 1
1227
    f1 = region.image_filled
1228

1229
    # Changed underlying image, but cache keeps result the same
1230
    assert_array_equal(f0, f1)
1231

1232
    # Now invalidate cache
1233
    region._cache_active = False
1234
    f1 = region.image_filled
1235

1236
    assert np.any(f0 != f1)
1237

1238

1239
def test_disabled_cache_is_empty():
1240
    SAMPLE_mod = SAMPLE.copy()
1241
    region = regionprops(SAMPLE_mod, cache=False)[0]
1242
    # Access one property to trigger cache
1243
    _ = region.image_filled
1244

1245
    # Cache should be empty
1246
    assert region._cache == dict()
1247

1248

1249
def test_docstrings_and_props():
1250
    def foo():
1251
        """foo"""
1252

1253
    has_docstrings = bool(foo.__doc__)
1254

1255
    region = regionprops(SAMPLE)[0]
1256

1257
    docs = _parse_docs()
1258
    props = [m for m in dir(region) if not m.startswith('_')]
1259

1260
    nr_docs_parsed = len(docs)
1261
    nr_props = len(props)
1262
    if has_docstrings:
1263
        assert_equal(nr_docs_parsed, nr_props)
1264
        ds = docs['moments_weighted_normalized']
1265
        assert 'iteration' not in ds
1266
        assert len(ds.split('\n')) > 3
1267
    else:
1268
        assert_equal(nr_docs_parsed, 0)
1269

1270

1271
def test_props_to_dict():
1272
    regions = regionprops(SAMPLE)
1273
    out = _props_to_dict(regions)
1274
    assert out == {
1275
        'label': np.array([1]),
1276
        'bbox-0': np.array([0]),
1277
        'bbox-1': np.array([0]),
1278
        'bbox-2': np.array([10]),
1279
        'bbox-3': np.array([18]),
1280
    }
1281

1282
    regions = regionprops(SAMPLE)
1283
    out = _props_to_dict(regions, properties=('label', 'area', 'bbox'), separator='+')
1284
    assert out == {
1285
        'label': np.array([1]),
1286
        'area': np.array([72]),
1287
        'bbox+0': np.array([0]),
1288
        'bbox+1': np.array([0]),
1289
        'bbox+2': np.array([10]),
1290
        'bbox+3': np.array([18]),
1291
    }
1292

1293

1294
def test_regionprops_table():
1295
    out = regionprops_table(SAMPLE)
1296
    assert out == {
1297
        'label': np.array([1]),
1298
        'bbox-0': np.array([0]),
1299
        'bbox-1': np.array([0]),
1300
        'bbox-2': np.array([10]),
1301
        'bbox-3': np.array([18]),
1302
    }
1303

1304
    out = regionprops_table(SAMPLE, properties=('label', 'area', 'bbox'), separator='+')
1305
    assert out == {
1306
        'label': np.array([1]),
1307
        'area': np.array([72]),
1308
        'bbox+0': np.array([0]),
1309
        'bbox+1': np.array([0]),
1310
        'bbox+2': np.array([10]),
1311
        'bbox+3': np.array([18]),
1312
    }
1313

1314

1315
def test_regionprops_table_deprecated_vector_property():
1316
    out = regionprops_table(SAMPLE, properties=('local_centroid',))
1317
    for key in out.keys():
1318
        # key reflects the deprecated name, not its new (centroid_local) value
1319
        assert key.startswith('local_centroid')
1320

1321

1322
def test_regionprops_table_deprecated_scalar_property():
1323
    out = regionprops_table(SAMPLE, properties=('bbox_area',))
1324
    assert list(out.keys()) == ['bbox_area']
1325

1326

1327
def test_regionprops_table_equal_to_original():
1328
    regions = regionprops(SAMPLE, INTENSITY_FLOAT_SAMPLE)
1329
    out_table = regionprops_table(
1330
        SAMPLE, INTENSITY_FLOAT_SAMPLE, properties=COL_DTYPES.keys()
1331
    )
1332

1333
    for prop, dtype in COL_DTYPES.items():
1334
        for i, reg in enumerate(regions):
1335
            rp = reg[prop]
1336
            if np.isscalar(rp) or prop in OBJECT_COLUMNS or dtype is np.object_:
1337
                assert_array_equal(rp, out_table[prop][i])
1338
            else:
1339
                shape = rp.shape if isinstance(rp, np.ndarray) else (len(rp),)
1340
                for ind in np.ndindex(shape):
1341
                    modified_prop = "-".join(map(str, (prop,) + ind))
1342
                    loc = ind if len(ind) > 1 else ind[0]
1343
                    assert_equal(rp[loc], out_table[modified_prop][i])
1344

1345

1346
def test_regionprops_table_no_regions():
1347
    out = regionprops_table(
1348
        np.zeros((2, 2), dtype=int), properties=('label', 'area', 'bbox'), separator='+'
1349
    )
1350
    assert len(out) == 6
1351
    assert len(out['label']) == 0
1352
    assert len(out['area']) == 0
1353
    assert len(out['bbox+0']) == 0
1354
    assert len(out['bbox+1']) == 0
1355
    assert len(out['bbox+2']) == 0
1356
    assert len(out['bbox+3']) == 0
1357

1358

1359
def test_column_dtypes_correct():
1360
    msg = 'mismatch with expected type,'
1361
    region = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE)[0]
1362
    for col in COL_DTYPES:
1363
        r = region[col]
1364

1365
        if col in OBJECT_COLUMNS:
1366
            assert COL_DTYPES[col] == object
1367
            continue
1368

1369
        t = type(np.ravel(r)[0])
1370

1371
        if np.issubdtype(t, np.floating):
1372
            assert COL_DTYPES[col] == float, f'{col} dtype {t} {msg} {COL_DTYPES[col]}'
1373
        elif np.issubdtype(t, np.integer):
1374
            assert COL_DTYPES[col] == int, f'{col} dtype {t} {msg} {COL_DTYPES[col]}'
1375
        else:
1376
            assert False, f'{col} dtype {t} {msg} {COL_DTYPES[col]}'
1377

1378

1379
def test_all_documented_items_in_col_dtypes():
1380
    numpydoc_docscrape = pytest.importorskip("numpydoc.docscrape")
1381
    docstring = numpydoc_docscrape.FunctionDoc(regionprops)
1382
    notes_lines = docstring['Notes']
1383
    property_lines = filter(lambda line: line.startswith('**'), notes_lines)
1384
    pattern = r'\*\*(?P<property_name>[a-z_]+)\*\*.*'
1385
    property_names = {
1386
        re.search(pattern, property_line).group('property_name')
1387
        for property_line in property_lines
1388
    }
1389
    column_keys = set(COL_DTYPES.keys())
1390
    assert column_keys == property_names
1391

1392

1393
def pixelcount(regionmask):
1394
    """a short test for an extra property"""
1395
    return np.sum(regionmask)
1396

1397

1398
def intensity_median(regionmask, image_intensity):
1399
    return np.median(image_intensity[regionmask])
1400

1401

1402
def bbox_list(regionmask):
1403
    """Extra property whose output shape is dependent on mask shape."""
1404
    return [1] * regionmask.shape[1]
1405

1406

1407
def too_many_args(regionmask, image_intensity, superfluous):
1408
    return 1
1409

1410

1411
def too_few_args():
1412
    return 1
1413

1414

1415
def test_extra_properties():
1416
    region = regionprops(SAMPLE, extra_properties=(pixelcount,))[0]
1417
    assert region.pixelcount == np.sum(SAMPLE == 1)
1418

1419

1420
def test_extra_properties_intensity():
1421
    region = regionprops(
1422
        SAMPLE, intensity_image=INTENSITY_SAMPLE, extra_properties=(intensity_median,)
1423
    )[0]
1424
    assert region.intensity_median == np.median(INTENSITY_SAMPLE[SAMPLE == 1])
1425

1426

1427
@pytest.mark.parametrize('intensity_prop', _require_intensity_image)
1428
def test_intensity_image_required(intensity_prop):
1429
    region = regionprops(SAMPLE)[0]
1430
    with pytest.raises(AttributeError) as e:
1431
        getattr(region, intensity_prop)
1432
    expected_error = (
1433
        f"Attribute '{intensity_prop}' unavailable when `intensity_image` has "
1434
        f"not been specified."
1435
    )
1436
    assert expected_error == str(e.value)
1437

1438

1439
def test_extra_properties_no_intensity_provided():
1440
    with pytest.raises(AttributeError):
1441
        region = regionprops(SAMPLE, extra_properties=(intensity_median,))[0]
1442
        _ = region.intensity_median
1443

1444

1445
def test_extra_properties_nr_args():
1446
    with pytest.raises(AttributeError):
1447
        region = regionprops(SAMPLE, extra_properties=(too_few_args,))[0]
1448
        _ = region.too_few_args
1449
    with pytest.raises(AttributeError):
1450
        region = regionprops(SAMPLE, extra_properties=(too_many_args,))[0]
1451
        _ = region.too_many_args
1452

1453

1454
def test_extra_properties_mixed():
1455
    # mixed properties, with and without intensity
1456
    region = regionprops(
1457
        SAMPLE,
1458
        intensity_image=INTENSITY_SAMPLE,
1459
        extra_properties=(intensity_median, pixelcount),
1460
    )[0]
1461
    assert region.intensity_median == np.median(INTENSITY_SAMPLE[SAMPLE == 1])
1462
    assert region.pixelcount == np.sum(SAMPLE == 1)
1463

1464

1465
def test_extra_properties_table():
1466
    out = regionprops_table(
1467
        SAMPLE_MULTIPLE,
1468
        intensity_image=INTENSITY_SAMPLE_MULTIPLE,
1469
        properties=('label',),
1470
        extra_properties=(intensity_median, pixelcount, bbox_list),
1471
    )
1472
    assert_array_almost_equal(out['intensity_median'], np.array([2.0, 4.0]))
1473
    assert_array_equal(out['pixelcount'], np.array([10, 2]))
1474

1475
    assert out['bbox_list'].dtype == np.object_
1476
    assert out["bbox_list"][0] == [1] * 10
1477
    assert out["bbox_list"][1] == [1] * 1
1478

1479

1480
def test_multichannel():
1481
    """Test that computing multichannel properties works."""
1482
    astro = data.astronaut()[::4, ::4]
1483
    astro_green = astro[..., 1]
1484
    labels = slic(astro.astype(float), start_label=1)
1485

1486
    segment_idx = np.max(labels) // 2
1487
    region = regionprops(labels, astro_green, extra_properties=[intensity_median])[
1488
        segment_idx
1489
    ]
1490
    region_multi = regionprops(labels, astro, extra_properties=[intensity_median])[
1491
        segment_idx
1492
    ]
1493

1494
    for prop in list(PROPS.keys()) + ["intensity_median"]:
1495
        p = region[prop]
1496
        p_multi = region_multi[prop]
1497
        if np.shape(p) == np.shape(p_multi):
1498
            # property does not depend on multiple channels
1499
            assert_array_equal(p, p_multi)
1500
        else:
1501
            # property uses multiple channels, returns props stacked along
1502
            # final axis
1503
            assert_allclose(p, np.asarray(p_multi)[..., 1], rtol=1e-12, atol=1e-12)
1504

1505

1506
def test_3d_ellipsoid_axis_lengths():
1507
    """Verify that estimated axis lengths are correct.
1508

1509
    Uses an ellipsoid at an arbitrary position and orientation.
1510
    """
1511
    # generate a centered ellipsoid with non-uniform half-lengths (radii)
1512
    half_lengths = (20, 10, 50)
1513
    e = draw.ellipsoid(*half_lengths).astype(int)
1514

1515
    # Pad by asymmetric amounts so the ellipse isn't centered. Also, pad enough
1516
    # that the rotated ellipse will still be within the original volume.
1517
    e = np.pad(e, pad_width=[(30, 18), (30, 12), (40, 20)], mode='constant')
1518

1519
    # apply rotations to the ellipsoid
1520
    R = transform.EuclideanTransform(rotation=[0.2, 0.3, 0.4], dimensionality=3)
1521
    e = ndi.affine_transform(e, R.params)
1522

1523
    # Compute regionprops
1524
    rp = regionprops(e)[0]
1525

1526
    # estimate principal axis lengths via the inertia tensor eigenvalues
1527
    evs = rp.inertia_tensor_eigvals
1528
    axis_lengths = _inertia_eigvals_to_axes_lengths_3D(evs)
1529
    expected_lengths = sorted([2 * h for h in half_lengths], reverse=True)
1530
    for ax_len_expected, ax_len in zip(expected_lengths, axis_lengths):
1531
        # verify accuracy to within 1%
1532
        assert abs(ax_len - ax_len_expected) < 0.01 * ax_len_expected
1533

1534
    # verify that the axis length regionprops also agree
1535
    assert abs(rp.axis_major_length - axis_lengths[0]) < 1e-7
1536
    assert abs(rp.axis_minor_length - axis_lengths[-1]) < 1e-7
1537