1import numpy as np
2import scipy.ndimage as ndi
3from . import _marching_cubes_classic_cy
4
5
6def _marching_cubes_classic(volume, level, spacing, gradient_direction):
7    """Lorensen et al. algorithm for marching cubes. See
8    marching_cubes_classic for documentation.
9
10    """
11    # Check inputs and ensure `volume` is C-contiguous for memoryviews
12    if volume.ndim != 3:
13        raise ValueError("Input volume must have 3 dimensions.")
14    if level is None:
15        level = 0.5 * (volume.min() + volume.max())
16    else:
17        level = float(level)
18        if level < volume.min() or level > volume.max():
19            raise ValueError("Surface level must be within volume data range.")
20    if len(spacing) != 3:
21        raise ValueError("`spacing` must consist of three floats.")
22
23    volume = np.array(volume, dtype=np.float64, order="C")
24
25    # Extract raw triangles using marching cubes in Cython
26    #   Returns a list of length-3 lists, each sub-list containing three
27    #   tuples. The tuples hold (x, y, z) coordinates for triangle vertices.
28    # Note: this algorithm is fast, but returns degenerate "triangles" which
29    #   have repeated vertices - and equivalent vertices are redundantly
30    #   placed in every triangle they connect with.
31    raw_faces = _marching_cubes_classic_cy.iterate_and_store_3d(volume,
32                                                                float(level))
33
34    # Find and collect unique vertices, storing triangle verts as indices.
35    # Returns a true mesh with no degenerate faces.
36    verts, faces = _marching_cubes_classic_cy.unpack_unique_verts(raw_faces)
37
38    verts = np.asarray(verts)
39    faces = np.asarray(faces)
40
41    # Fancy indexing to define two vector arrays from triangle vertices
42    faces = _correct_mesh_orientation(volume, verts[faces], faces, spacing,
43                                      gradient_direction)
44
45    # Adjust for non-isotropic spacing in `verts` at time of return
46    return verts * np.r_[spacing], faces
47
48
49def mesh_surface_area(verts, faces):
50    """
51    Compute surface area, given vertices & triangular faces
52
53    Parameters
54    ----------
55    verts : (V, 3) array of floats
56        Array containing (x, y, z) coordinates for V unique mesh vertices.
57    faces : (F, 3) array of ints
58        List of length-3 lists of integers, referencing vertex coordinates as
59        provided in `verts`
60
61    Returns
62    -------
63    area : float
64        Surface area of mesh. Units now [coordinate units] ** 2.
65
66    Notes
67    -----
68    The arguments expected by this function are the first two outputs from
69    `skimage.measure.marching_cubes`. For unit correct output, ensure correct
70    `spacing` was passed to `skimage.measure.marching_cubes`.
71
72    This algorithm works properly only if the ``faces`` provided are all
73    triangles.
74
75    See Also
76    --------
77    skimage.measure.marching_cubes
78    skimage.measure.marching_cubes_classic
79
80    """
81    # Fancy indexing to define two vector arrays from triangle vertices
82    actual_verts = verts[faces]
83    a = actual_verts[:, 0, :] - actual_verts[:, 1, :]
84    b = actual_verts[:, 0, :] - actual_verts[:, 2, :]
85    del actual_verts
86
87    # Area of triangle in 3D = 1/2 * Euclidean norm of cross product
88    return ((np.cross(a, b) ** 2).sum(axis=1) ** 0.5).sum() / 2.
89
90
91def _correct_mesh_orientation(volume, actual_verts, faces,
92                              spacing=(1., 1., 1.),
93                              gradient_direction='descent'):
94    """
95    Correct orientations of mesh faces.
96
97    Parameters
98    ----------
99    volume : (M, N, P) array of doubles
100        Input data volume to find isosurfaces. Will be cast to `np.float64`.
101    actual_verts : (F, 3, 3) array of floats
102        Array with (face, vertex, coords) index coordinates.
103    faces : (F, 3) array of ints
104        List of length-3 lists of integers, referencing vertex coordinates as
105        provided in `verts`.
106    spacing : length-3 tuple of floats
107        Voxel spacing in spatial dimensions corresponding to numpy array
108        indexing dimensions (M, N, P) as in `volume`.
109    gradient_direction : string
110        Controls if the mesh was generated from an isosurface with gradient
111        descent toward objects of interest (the default), or the opposite.
112        The two options are:
113        * descent : Object was greater than exterior
114        * ascent : Exterior was greater than object
115
116    Returns
117    -------
118    faces_corrected (F, 3) array of ints
119        Corrected list of faces referencing vertex coordinates in `verts`.
120
121    Notes
122    -----
123    Certain applications and mesh processing algorithms require all faces
124    to be oriented in a consistent way. Generally, this means a normal vector
125    points "out" of the meshed shapes. This algorithm corrects the output from
126    `skimage.measure.marching_cubes_classic` by flipping the orientation of
127    mis-oriented faces.
128
129    Because marching cubes could be used to find isosurfaces either on
130    gradient descent (where the desired object has greater values than the
131    exterior) or ascent (where the desired object has lower values than the
132    exterior), the ``gradient_direction`` kwarg allows the user to inform this
133    algorithm which is correct. If the resulting mesh appears to be oriented
134    completely incorrectly, try changing this option.
135
136    The arguments expected by this function are the exact outputs from
137    `skimage.measure.marching_cubes_classic` except `actual_verts`, which is an
138    uncorrected version of the fancy indexing operation `verts[faces]`.
139    Only `faces` is corrected and returned as the vertices do not change,
140    only the order in which they are referenced.
141
142    This algorithm assumes ``faces`` provided are exclusively triangles.
143
144    See Also
145    --------
146    skimage.measure.marching_cubes_classic
147    skimage.measure.mesh_surface_area
148
149    """
150    # Calculate gradient of `volume`, then interpolate to vertices in `verts`
151    grad_x, grad_y, grad_z = np.gradient(volume)
152
153    a = actual_verts[:, 0, :] - actual_verts[:, 1, :]
154    b = actual_verts[:, 0, :] - actual_verts[:, 2, :]
155
156    # Find triangle centroids
157    centroids = (actual_verts.sum(axis=1) / 3.).T
158
159    del actual_verts
160
161    # Interpolate face centroids into each gradient axis
162    grad_centroids_x = ndi.map_coordinates(grad_x, centroids)
163    grad_centroids_y = ndi.map_coordinates(grad_y, centroids)
164    grad_centroids_z = ndi.map_coordinates(grad_z, centroids)
165
166    # Combine and normalize interpolated gradients
167    grad_centroids = np.c_[grad_centroids_x, grad_centroids_y,
168                           grad_centroids_z]
169    grad_centroids = (grad_centroids /
170                      (np.sum(grad_centroids ** 2,
171                              axis=1) ** 0.5)[:, np.newaxis])
172
173    # Find normal vectors for each face via cross product
174    crosses = np.cross(a, b)
175    crosses = crosses / (np.sum(crosses ** 2, axis=1) ** (0.5))[:, np.newaxis]
176
177    # Take dot product
178    dotproducts = (grad_centroids * crosses).sum(axis=1)
179
180    # Find mis-oriented faces
181    if 'descent' in gradient_direction:
182        # Faces with incorrect orientations have dot product < 0
183        indices = (dotproducts < 0).nonzero()[0]
184    elif 'ascent' in gradient_direction:
185        # Faces with incorrection orientation have dot product > 0
186        indices = (dotproducts > 0).nonzero()[0]
187    else:
188        raise ValueError("Incorrect input %s in `gradient_direction`, see "
189                         "docstring." % (gradient_direction))
190
191    # Correct orientation and return, without modifying original data
192    faces_corrected = faces.copy()
193    faces_corrected[indices] = faces_corrected[indices, ::-1]
194
195    return faces_corrected
196