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