Python skimage.measure.marching_cubes() Examples

def visualize_voxel_spectral(points, vis_size=128):
  """Function to visualize voxel (spectral)."""
  points = np.rint(points)
  points = np.swapaxes(points, 0, 2)
  fig = p.figure(figsize=(1, 1), dpi=vis_size)
  verts, faces = measure.marching_cubes(points, 0, spacing=(0.1, 0.1, 0.1))
  ax = fig.add_subplot(111, projection='3d')
  ax.plot_trisurf(
      verts[:, 0], verts[:, 1], faces, verts[:, 2], cmap='Spectral_r', lw=0.1)
  ax.set_axis_off()
  fig.tight_layout(pad=0)
  fig.canvas.draw()
  data = np.fromstring(
      fig.canvas.tostring_rgb(), dtype=np.uint8, sep='').reshape(
          vis_size, vis_size, 3)
  p.close('all')
  return data 

def plot_3d(image, threshold=-300):
    # Position the scan upright, 
    # so the head of the patient would be at the top facing the camera
    p = image.transpose(2,1,0)

    #p = image
    
    verts, faces = measure.marching_cubes(p, threshold)

    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111, projection='3d')

    # Fancy indexing: `verts[faces]` to generate a collection of triangles
    mesh = Poly3DCollection(verts[faces], alpha=0.70)
    face_color = [0.45, 0.45, 0.75]
    mesh.set_facecolor(face_color)
    ax.add_collection3d(mesh)

    ax.set_xlim(0, p.shape[0])
    ax.set_ylim(0, p.shape[1])
    ax.set_zlim(0, p.shape[2])

    plt.show() 

def save_mesh (binary, path):
    binary = mesh.pad(binary, dtype=np.float) 
    binary = gaussian_filter(binary, 2, mode='constant')
    verts, faces = measure.marching_cubes(binary, 0.5)
    Three(path, verts, faces) 

def getVFByMarchingCubes(voxels, threshold=0.5):
    v, f =  sk.marching_cubes(voxels, level=threshold)
    return v, f 

def create_surf(subject, subcortical_code, subcortical_name):
    aseg = nib.load(op.join(SUBJECTS_DIR, subject, 'mri', 'aseg.mgz')).get_data()
    t1_header = nib.load(op.join(SUBJECTS_DIR, subject, 'mri', 'T1.mgz')).header
    aseg[aseg != subcortical_code] = 255
    aseg[aseg == subcortical_code] = 10
    aseg = np.array(aseg, dtype=np.float)
    aseg_smooth = scipy.ndimage.gaussian_filter(aseg, sigma=1)
    verts_vox, faces, _, _ = measure.marching_cubes(aseg_smooth, 100)
    # Doesn't seem to fix the normals directions that should be out...
    faces = measure.correct_mesh_orientation(aseg_smooth, verts_vox, faces)
    verts = utils.apply_trans(t1_header.get_vox2ras_tkr(), verts_vox)
    fol = utils.make_dir(op.join(MMVT_DIR, subject, 'subcortical_test'))
    ply_file_name = op.join(fol, '{}.ply'.format(subcortical_name))
    utils.write_ply_file(verts, faces, ply_file_name, False) 

def _assign_layers(self):
        """ There are no layers in the Willard-Chandler method.

            This function identifies the dividing surface and stores the
            triangulated isosurface, the density and the particles.

        """
        self.reset_labels()
        # we assign an empty group for consistency
        self._layers, self.normal = self.universe.atoms[:0], None

        self.prepare_box()

        self._define_cluster_group()

        self.centered_positions = None
        if self.do_center is True:
            self.center()

        pos = self.cluster_group.positions
        box = self.universe.dimensions[:3]

        ngrid, spacing = utilities.compute_compatible_mesh_params(
            self.mesh, box)
        self.spacing, self.ngrid = spacing, ngrid
        grid = utilities.generate_grid_in_box(box, ngrid, order='xyz')
        kernel, _ = utilities.density_map(pos, grid, self.alpha, box)

        kernel.pos = pos.copy()
        self.density_field = kernel.evaluate_pbc_fast(grid)

        # Thomas Lewiner, Helio Lopes, Antonio Wilson Vieira and Geovan
        # Tavares. Efficient implementation of Marching Cubes’ cases with
        # topological guarantees. Journal of Graphics Tools 8(2) pp. 1-15
        # (december 2003). DOI: 10.1080/10867651.2003.10487582
        volume = self.density_field.reshape(
            tuple(np.array(ngrid[::-1]).astype(int)))
        verts, faces, normals, values = marching_cubes(
            volume, None, spacing=tuple(spacing))
        # note that len(normals) == len(verts): they are normals
        # at the vertices, and not normals of the faces
        # verts and normals have x and z flipped because skimage uses zyx ordering
        self.triangulated_surface = [np.fliplr(verts), faces, np.fliplr(normals)]
        self.surface_area = measure.mesh_surface_area(verts, faces)
        verts += spacing[::-1] / 2. 

def compute_marching_cubes_regular_grid(self, level: float, scalar_field,
                                            mask_array=None,
                                            rescale=False, **kwargs):
        """Compute the surface (vertices and edges) of a given surface by computing
         marching cubes (by skimage)

        Args:
            level (float): value of the scalar field at the surface
            scalar_field (np.array): scalar_field vector objects
            mask_array (np.array): mask vector with trues where marching cubes has to be performed
            rescale (bool): if True surfaces will be located between 0 and 1
            **kwargs: skimage.measure.marching_cubes_lewiner args (see below)

        Returns:
            list: vertices, simplices, normals, values

        See Also:

            :func:`skimage.measure.marching_cubes`

        """

        vertices, simplices, normals, values = measure.marching_cubes(
            scalar_field.reshape(self.grid.regular_grid.resolution[0],
                                 self.grid.regular_grid.resolution[1],
                                 self.grid.regular_grid.resolution[2]),
            level,
            spacing=self.grid.regular_grid.get_dx_dy_dz(rescale=rescale),
            mask=mask_array,
            **kwargs)

        if rescale is True:
            loc_0 = self.grid.regular_grid.extent_r[[0, 2, 4]] + \
                    np.array(self.grid.regular_grid.get_dx_dy_dz(rescale=True)) / 2

            vertices += np.array(loc_0).reshape(1, 3)

        else:
            loc_0 = self.grid.regular_grid.extent[[0, 2, 4]] + \
                    np.array(self.grid.regular_grid.get_dx_dy_dz(rescale=False)) / 2
            vertices += np.array(loc_0).reshape(1, 3)

        return [vertices, simplices, normals, values]

    # def mask_topo(self, mask_matrix):
    #     """Add the masked elements of the topography to the masking matrix"""
    #     x = ~self.grid.regular_grid.mask_topo
    #     a = (np.swapaxes(x, 0, 1) * mask_matrix)
    #     return a