scikit-image

1
import numpy as np
2
import pytest
3
from numpy.testing import assert_allclose
4

5
from skimage.draw import ellipsoid, ellipsoid_stats
6
from skimage.measure import marching_cubes, mesh_surface_area
7

8

9
def test_marching_cubes_isotropic():
10
    ellipsoid_isotropic = ellipsoid(6, 10, 16, levelset=True)
11
    _, surf = ellipsoid_stats(6, 10, 16)
12

13
    # Classic
14
    verts, faces = marching_cubes(ellipsoid_isotropic, 0.0, method='lorensen')[:2]
15
    surf_calc = mesh_surface_area(verts, faces)
16
    # Test within 1% tolerance for isotropic. Will always underestimate.
17
    assert surf > surf_calc and surf_calc > surf * 0.99
18

19
    # Lewiner
20
    verts, faces = marching_cubes(ellipsoid_isotropic, 0.0)[:2]
21
    surf_calc = mesh_surface_area(verts, faces)
22
    # Test within 1% tolerance for isotropic. Will always underestimate.
23
    assert surf > surf_calc and surf_calc > surf * 0.99
24

25

26
def test_marching_cubes_anisotropic():
27
    # test spacing as numpy array (and not just tuple)
28
    spacing = np.array([1.0, 10 / 6.0, 16 / 6.0])
29
    ellipsoid_anisotropic = ellipsoid(6, 10, 16, spacing=spacing, levelset=True)
30
    _, surf = ellipsoid_stats(6, 10, 16)
31

32
    # Classic
33
    verts, faces = marching_cubes(
34
        ellipsoid_anisotropic, 0.0, spacing=spacing, method='lorensen'
35
    )[:2]
36
    surf_calc = mesh_surface_area(verts, faces)
37
    # Test within 1.5% tolerance for anisotropic. Will always underestimate.
38
    assert surf > surf_calc and surf_calc > surf * 0.985
39

40
    # Lewiner
41
    verts, faces = marching_cubes(ellipsoid_anisotropic, 0.0, spacing=spacing)[:2]
42
    surf_calc = mesh_surface_area(verts, faces)
43
    # Test within 1.5% tolerance for anisotropic. Will always underestimate.
44
    assert surf > surf_calc and surf_calc > surf * 0.985
45

46
    # Test marching cube with mask
47
    with pytest.raises(ValueError):
48
        verts, faces = marching_cubes(
49
            ellipsoid_anisotropic, 0.0, spacing=spacing, mask=np.array([])
50
        )[:2]
51

52
    # Test spacing together with allow_degenerate=False
53
    marching_cubes(ellipsoid_anisotropic, 0, spacing=spacing, allow_degenerate=False)
54

55

56
def test_invalid_input():
57
    # Classic
58
    with pytest.raises(ValueError):
59
        marching_cubes(np.zeros((2, 2, 1)), 0, method='lorensen')
60
    with pytest.raises(ValueError):
61
        marching_cubes(np.zeros((2, 2, 1)), 1, method='lorensen')
62
    with pytest.raises(ValueError):
63
        marching_cubes(np.ones((3, 3, 3)), 1, spacing=(1, 2), method='lorensen')
64
    with pytest.raises(ValueError):
65
        marching_cubes(np.zeros((20, 20)), 0, method='lorensen')
66

67
    # Lewiner
68
    with pytest.raises(ValueError):
69
        marching_cubes(np.zeros((2, 2, 1)), 0)
70
    with pytest.raises(ValueError):
71
        marching_cubes(np.zeros((2, 2, 1)), 1)
72
    with pytest.raises(ValueError):
73
        marching_cubes(np.ones((3, 3, 3)), 1, spacing=(1, 2))
74
    with pytest.raises(ValueError):
75
        marching_cubes(np.zeros((20, 20)), 0)
76

77
    # invalid method name
78
    ellipsoid_isotropic = ellipsoid(6, 10, 16, levelset=True)
79
    with pytest.raises(ValueError):
80
        marching_cubes(ellipsoid_isotropic, 0.0, method='abcd')
81

82

83
def test_both_algs_same_result_ellipse():
84
    # Performing this test on data that does not have ambiguities
85

86
    sphere_small = ellipsoid(1, 1, 1, levelset=True)
87

88
    vertices1, faces1 = marching_cubes(sphere_small, 0, allow_degenerate=False)[:2]
89
    vertices2, faces2 = marching_cubes(
90
        sphere_small, 0, allow_degenerate=False, method='lorensen'
91
    )[:2]
92

93
    # Order is different, best we can do is test equal shape and same
94
    # vertices present
95
    assert _same_mesh(vertices1, faces1, vertices2, faces2)
96

97

98
def _same_mesh(vertices1, faces1, vertices2, faces2, tol=1e-10):
99
    """Compare two meshes, using a certain tolerance and invariant to
100
    the order of the faces.
101
    """
102
    # Unwind vertices
103
    triangles1 = vertices1[np.array(faces1)]
104
    triangles2 = vertices2[np.array(faces2)]
105
    # Sort vertices within each triangle
106
    triang1 = [np.concatenate(sorted(t, key=lambda x: tuple(x))) for t in triangles1]
107
    triang2 = [np.concatenate(sorted(t, key=lambda x: tuple(x))) for t in triangles2]
108
    # Sort the resulting 9-element "tuples"
109
    triang1 = np.array(sorted([tuple(x) for x in triang1]))
110
    triang2 = np.array(sorted([tuple(x) for x in triang2]))
111
    return triang1.shape == triang2.shape and np.allclose(triang1, triang2, 0, tol)
112

113

114
def test_both_algs_same_result_donut():
115
    # Performing this test on data that does not have ambiguities
116
    n = 48
117
    a, b = 2.5 / n, -1.25
118

119
    vol = np.empty((n, n, n), 'float32')
120
    for iz in range(vol.shape[0]):
121
        for iy in range(vol.shape[1]):
122
            for ix in range(vol.shape[2]):
123
                # Double-torii formula by Thomas Lewiner
124
                z, y, x = float(iz) * a + b, float(iy) * a + b, float(ix) * a + b
125
                vol[iz, iy, ix] = (
126
                    ((8 * x) ** 2 + (8 * y - 2) ** 2 + (8 * z) ** 2 + 16 - 1.85 * 1.85)
127
                    * (
128
                        (8 * x) ** 2
129
                        + (8 * y - 2) ** 2
130
                        + (8 * z) ** 2
131
                        + 16
132
                        - 1.85 * 1.85
133
                    )
134
                    - 64 * ((8 * x) ** 2 + (8 * y - 2) ** 2)
135
                ) * (
136
                    (
137
                        (8 * x) ** 2
138
                        + ((8 * y - 2) + 4) * ((8 * y - 2) + 4)
139
                        + (8 * z) ** 2
140
                        + 16
141
                        - 1.85 * 1.85
142
                    )
143
                    * (
144
                        (8 * x) ** 2
145
                        + ((8 * y - 2) + 4) * ((8 * y - 2) + 4)
146
                        + (8 * z) ** 2
147
                        + 16
148
                        - 1.85 * 1.85
149
                    )
150
                    - 64 * (((8 * y - 2) + 4) * ((8 * y - 2) + 4) + (8 * z) ** 2)
151
                ) + 1025
152

153
    vertices1, faces1 = marching_cubes(vol, 0, method='lorensen')[:2]
154
    vertices2, faces2 = marching_cubes(vol, 0)[:2]
155

156
    # Old and new alg are different
157
    assert not _same_mesh(vertices1, faces1, vertices2, faces2)
158

159

160
def test_masked_marching_cubes():
161
    ellipsoid_scalar = ellipsoid(6, 10, 16, levelset=True)
162
    mask = np.ones_like(ellipsoid_scalar, dtype=bool)
163
    mask[:10, :, :] = False
164
    mask[:, :, 20:] = False
165
    ver, faces, _, _ = marching_cubes(ellipsoid_scalar, 0, mask=mask)
166
    area = mesh_surface_area(ver, faces)
167

168
    assert_allclose(area, 299.56878662109375, rtol=0.01)
169

170

171
def test_masked_marching_cubes_empty():
172
    ellipsoid_scalar = ellipsoid(6, 10, 16, levelset=True)
173
    mask = np.array([])
174
    with pytest.raises(ValueError):
175
        _ = marching_cubes(ellipsoid_scalar, 0, mask=mask)
176

177

178
def test_masked_marching_cubes_all_true():
179
    ellipsoid_scalar = ellipsoid(6, 10, 16, levelset=True)
180
    mask = np.ones_like(ellipsoid_scalar, dtype=bool)
181
    ver_m, faces_m, _, _ = marching_cubes(ellipsoid_scalar, 0, mask=mask)
182
    ver, faces, _, _ = marching_cubes(ellipsoid_scalar, 0, mask=mask)
183
    assert_allclose(ver_m, ver, rtol=0.00001)
184
    assert_allclose(faces_m, faces, rtol=0.00001)
185