scikit-image

1
import base64
2

3
import numpy as np
4

5
from . import _marching_cubes_lewiner_luts as mcluts
6
from . import _marching_cubes_lewiner_cy
7

8

9
def marching_cubes(
10
    volume,
11
    level=None,
12
    *,
13
    spacing=(1.0, 1.0, 1.0),
14
    gradient_direction='descent',
15
    step_size=1,
16
    allow_degenerate=True,
17
    method='lewiner',
18
    mask=None,
19
):
20
    """Marching cubes algorithm to find surfaces in 3d volumetric data.
21

22
    In contrast with Lorensen et al. approach [2]_, Lewiner et
23
    al. algorithm is faster, resolves ambiguities, and guarantees
24
    topologically correct results. Therefore, this algorithm generally
25
    a better choice.
26

27
    Parameters
28
    ----------
29
    volume : (M, N, P) ndarray
30
        Input data volume to find isosurfaces. Will internally be
31
        converted to float32 if necessary.
32
    level : float, optional
33
        Contour value to search for isosurfaces in `volume`. If not
34
        given or None, the average of the min and max of vol is used.
35
    spacing : length-3 tuple of floats, optional
36
        Voxel spacing in spatial dimensions corresponding to numpy array
37
        indexing dimensions (M, N, P) as in `volume`.
38
    gradient_direction : string, optional
39
        Controls if the mesh was generated from an isosurface with gradient
40
        descent toward objects of interest (the default), or the opposite,
41
        considering the *left-hand* rule.
42
        The two options are:
43
        * descent : Object was greater than exterior
44
        * ascent : Exterior was greater than object
45
    step_size : int, optional
46
        Step size in voxels. Default 1. Larger steps yield faster but
47
        coarser results. The result will always be topologically correct
48
        though.
49
    allow_degenerate : bool, optional
50
        Whether to allow degenerate (i.e. zero-area) triangles in the
51
        end-result. Default True. If False, degenerate triangles are
52
        removed, at the cost of making the algorithm slower.
53
    method: {'lewiner', 'lorensen'}, optional
54
        Whether the method of Lewiner et al. or Lorensen et al. will be used.
55
    mask : (M, N, P) array, optional
56
        Boolean array. The marching cube algorithm will be computed only on
57
        True elements. This will save computational time when interfaces
58
        are located within certain region of the volume M, N, P-e.g. the top
59
        half of the cube-and also allow to compute finite surfaces-i.e. open
60
        surfaces that do not end at the border of the cube.
61

62
    Returns
63
    -------
64
    verts : (V, 3) array
65
        Spatial coordinates for V unique mesh vertices. Coordinate order
66
        matches input `volume` (M, N, P). If ``allow_degenerate`` is set to
67
        True, then the presence of degenerate triangles in the mesh can make
68
        this array have duplicate vertices.
69
    faces : (F, 3) array
70
        Define triangular faces via referencing vertex indices from ``verts``.
71
        This algorithm specifically outputs triangles, so each face has
72
        exactly three indices.
73
    normals : (V, 3) array
74
        The normal direction at each vertex, as calculated from the
75
        data.
76
    values : (V,) array
77
        Gives a measure for the maximum value of the data in the local region
78
        near each vertex. This can be used by visualization tools to apply
79
        a colormap to the mesh.
80

81
    See Also
82
    --------
83
    skimage.measure.mesh_surface_area
84
    skimage.measure.find_contours
85

86
    Notes
87
    -----
88
    The algorithm [1]_ is an improved version of Chernyaev's Marching
89
    Cubes 33 algorithm. It is an efficient algorithm that relies on
90
    heavy use of lookup tables to handle the many different cases,
91
    keeping the algorithm relatively easy. This implementation is
92
    written in Cython, ported from Lewiner's C++ implementation.
93

94
    To quantify the area of an isosurface generated by this algorithm, pass
95
    verts and faces to `skimage.measure.mesh_surface_area`.
96

97
    Regarding visualization of algorithm output, to contour a volume
98
    named `myvolume` about the level 0.0, using the ``mayavi`` package::
99

100
      >>>
101
      >> from mayavi import mlab
102
      >> verts, faces, _, _ = marching_cubes(myvolume, 0.0)
103
      >> mlab.triangular_mesh([vert[0] for vert in verts],
104
                              [vert[1] for vert in verts],
105
                              [vert[2] for vert in verts],
106
                              faces)
107
      >> mlab.show()
108

109
    Similarly using the ``visvis`` package::
110

111
      >>>
112
      >> import visvis as vv
113
      >> verts, faces, normals, values = marching_cubes(myvolume, 0.0)
114
      >> vv.mesh(np.fliplr(verts), faces, normals, values)
115
      >> vv.use().Run()
116

117
    To reduce the number of triangles in the mesh for better performance,
118
    see this `example
119
    <https://docs.enthought.com/mayavi/mayavi/auto/example_julia_set_decimation.html#example-julia-set-decimation>`_
120
    using the ``mayavi`` package.
121

122
    References
123
    ----------
124
    .. [1] Thomas Lewiner, Helio Lopes, Antonio Wilson Vieira and Geovan
125
           Tavares. Efficient implementation of Marching Cubes' cases with
126
           topological guarantees. Journal of Graphics Tools 8(2)
127
           pp. 1-15 (december 2003).
128
           :DOI:`10.1080/10867651.2003.10487582`
129
    .. [2] Lorensen, William and Harvey E. Cline. Marching Cubes: A High
130
           Resolution 3D Surface Construction Algorithm. Computer Graphics
131
           (SIGGRAPH 87 Proceedings) 21(4) July 1987, p. 163-170).
132
           :DOI:`10.1145/37401.37422`
133
    """
134
    use_classic = False
135
    if method == 'lorensen':
136
        use_classic = True
137
    elif method != 'lewiner':
138
        raise ValueError("method should be either 'lewiner' or 'lorensen'")
139
    return _marching_cubes_lewiner(
140
        volume,
141
        level,
142
        spacing,
143
        gradient_direction,
144
        step_size,
145
        allow_degenerate,
146
        use_classic=use_classic,
147
        mask=mask,
148
    )
149

150

151
def _marching_cubes_lewiner(
152
    volume,
153
    level,
154
    spacing,
155
    gradient_direction,
156
    step_size,
157
    allow_degenerate,
158
    use_classic,
159
    mask,
160
):
161
    """Lewiner et al. algorithm for marching cubes. See
162
    marching_cubes_lewiner for documentation.
163

164
    """
165

166
    # Check volume and ensure its in the format that the alg needs
167
    if not isinstance(volume, np.ndarray) or (volume.ndim != 3):
168
        raise ValueError('Input volume should be a 3D numpy array.')
169
    if volume.shape[0] < 2 or volume.shape[1] < 2 or volume.shape[2] < 2:
170
        raise ValueError("Input array must be at least 2x2x2.")
171
    volume = np.ascontiguousarray(volume, np.float32)  # no copy if not necessary
172

173
    # Check/convert other inputs:
174
    # level
175
    if level is None:
176
        level = 0.5 * (volume.min() + volume.max())
177
    else:
178
        level = float(level)
179
        if level < volume.min() or level > volume.max():
180
            raise ValueError("Surface level must be within volume data range.")
181
    # spacing
182
    if len(spacing) != 3:
183
        raise ValueError("`spacing` must consist of three floats.")
184
    # step_size
185
    step_size = int(step_size)
186
    if step_size < 1:
187
        raise ValueError('step_size must be at least one.')
188
    # use_classic
189
    use_classic = bool(use_classic)
190

191
    # Get LutProvider class (reuse if possible)
192
    L = _get_mc_luts()
193

194
    # Check if a mask array is passed
195
    if mask is not None:
196
        if not mask.shape == volume.shape:
197
            raise ValueError('volume and mask must have the same shape.')
198

199
    # Apply algorithm
200
    func = _marching_cubes_lewiner_cy.marching_cubes
201
    vertices, faces, normals, values = func(
202
        volume, level, L, step_size, use_classic, mask
203
    )
204

205
    if not len(vertices):
206
        raise RuntimeError('No surface found at the given iso value.')
207

208
    # Output in z-y-x order, as is common in skimage
209
    vertices = np.fliplr(vertices)
210
    normals = np.fliplr(normals)
211

212
    # Finishing touches to output
213
    faces.shape = -1, 3
214
    if gradient_direction == 'descent':
215
        # MC implementation is right-handed, but gradient_direction is
216
        # left-handed
217
        faces = np.fliplr(faces)
218
    elif not gradient_direction == 'ascent':
219
        raise ValueError(
220
            f"Incorrect input {gradient_direction} in `gradient_direction`, "
221
            "see docstring."
222
        )
223
    if not np.array_equal(spacing, (1, 1, 1)):
224
        vertices = vertices * np.r_[spacing]
225

226
    if allow_degenerate:
227
        return vertices, faces, normals, values
228
    else:
229
        fun = _marching_cubes_lewiner_cy.remove_degenerate_faces
230
        return fun(vertices.astype(np.float32), faces, normals, values)
231

232

233
def _to_array(args):
234
    shape, text = args
235
    byts = base64.decodebytes(text.encode('utf-8'))
236
    ar = np.frombuffer(byts, dtype='int8')
237
    ar.shape = shape
238
    return ar
239

240

241
# Map an edge-index to two relative pixel positions. The edge index
242
# represents a point that lies somewhere in between these pixels.
243
# Linear interpolation should be used to determine where it is exactly.
244
#   0
245
# 3   1   ->  0x
246
#   2         xx
247

248
# fmt: off
249
EDGETORELATIVEPOSX = np.array([ [0,1],[1,1],[1,0],[0,0], [0,1],[1,1],[1,0],[0,0], [0,0],[1,1],[1,1],[0,0] ], 'int8')
250
EDGETORELATIVEPOSY = np.array([ [0,0],[0,1],[1,1],[1,0], [0,0],[0,1],[1,1],[1,0], [0,0],[0,0],[1,1],[1,1] ], 'int8')
251
EDGETORELATIVEPOSZ = np.array([ [0,0],[0,0],[0,0],[0,0], [1,1],[1,1],[1,1],[1,1], [0,1],[0,1],[0,1],[0,1] ], 'int8')
252
# fmt: on
253

254

255
def _get_mc_luts():
256
    """Kind of lazy obtaining of the luts."""
257
    if not hasattr(mcluts, 'THE_LUTS'):
258
        mcluts.THE_LUTS = _marching_cubes_lewiner_cy.LutProvider(
259
            EDGETORELATIVEPOSX,
260
            EDGETORELATIVEPOSY,
261
            EDGETORELATIVEPOSZ,
262
            _to_array(mcluts.CASESCLASSIC),
263
            _to_array(mcluts.CASES),
264
            _to_array(mcluts.TILING1),
265
            _to_array(mcluts.TILING2),
266
            _to_array(mcluts.TILING3_1),
267
            _to_array(mcluts.TILING3_2),
268
            _to_array(mcluts.TILING4_1),
269
            _to_array(mcluts.TILING4_2),
270
            _to_array(mcluts.TILING5),
271
            _to_array(mcluts.TILING6_1_1),
272
            _to_array(mcluts.TILING6_1_2),
273
            _to_array(mcluts.TILING6_2),
274
            _to_array(mcluts.TILING7_1),
275
            _to_array(mcluts.TILING7_2),
276
            _to_array(mcluts.TILING7_3),
277
            _to_array(mcluts.TILING7_4_1),
278
            _to_array(mcluts.TILING7_4_2),
279
            _to_array(mcluts.TILING8),
280
            _to_array(mcluts.TILING9),
281
            _to_array(mcluts.TILING10_1_1),
282
            _to_array(mcluts.TILING10_1_1_),
283
            _to_array(mcluts.TILING10_1_2),
284
            _to_array(mcluts.TILING10_2),
285
            _to_array(mcluts.TILING10_2_),
286
            _to_array(mcluts.TILING11),
287
            _to_array(mcluts.TILING12_1_1),
288
            _to_array(mcluts.TILING12_1_1_),
289
            _to_array(mcluts.TILING12_1_2),
290
            _to_array(mcluts.TILING12_2),
291
            _to_array(mcluts.TILING12_2_),
292
            _to_array(mcluts.TILING13_1),
293
            _to_array(mcluts.TILING13_1_),
294
            _to_array(mcluts.TILING13_2),
295
            _to_array(mcluts.TILING13_2_),
296
            _to_array(mcluts.TILING13_3),
297
            _to_array(mcluts.TILING13_3_),
298
            _to_array(mcluts.TILING13_4),
299
            _to_array(mcluts.TILING13_5_1),
300
            _to_array(mcluts.TILING13_5_2),
301
            _to_array(mcluts.TILING14),
302
            _to_array(mcluts.TEST3),
303
            _to_array(mcluts.TEST4),
304
            _to_array(mcluts.TEST6),
305
            _to_array(mcluts.TEST7),
306
            _to_array(mcluts.TEST10),
307
            _to_array(mcluts.TEST12),
308
            _to_array(mcluts.TEST13),
309
            _to_array(mcluts.SUBCONFIG13),
310
        )
311

312
    return mcluts.THE_LUTS
313

314

315
def mesh_surface_area(verts, faces):
316
    """Compute surface area, given vertices and triangular faces.
317

318
    Parameters
319
    ----------
320
    verts : (V, 3) array of floats
321
        Array containing coordinates for V unique mesh vertices.
322
    faces : (F, 3) array of ints
323
        List of length-3 lists of integers, referencing vertex coordinates as
324
        provided in `verts`.
325

326
    Returns
327
    -------
328
    area : float
329
        Surface area of mesh. Units now [coordinate units] ** 2.
330

331
    Notes
332
    -----
333
    The arguments expected by this function are the first two outputs from
334
    `skimage.measure.marching_cubes`. For unit correct output, ensure correct
335
    `spacing` was passed to `skimage.measure.marching_cubes`.
336

337
    This algorithm works properly only if the ``faces`` provided are all
338
    triangles.
339

340
    See Also
341
    --------
342
    skimage.measure.marching_cubes
343

344
    """
345
    # Fancy indexing to define two vector arrays from triangle vertices
346
    actual_verts = verts[faces]
347
    a = actual_verts[:, 0, :] - actual_verts[:, 1, :]
348
    b = actual_verts[:, 0, :] - actual_verts[:, 2, :]
349
    del actual_verts
350

351
    # Area of triangle in 3D = 1/2 * Euclidean norm of cross product
352
    return ((np.cross(a, b) ** 2).sum(axis=1) ** 0.5).sum() / 2.0
353