Python skimage.measure.marching_cubes_lewiner() Examples

def get_point_cloud(self):
    """Extract a point cloud from the voxel volume.
    """
    tsdf_vol, color_vol = self.get_volume()

    # Marching cubes
    verts = measure.marching_cubes_lewiner(tsdf_vol, level=0)[0]
    verts_ind = np.round(verts).astype(int)
    verts = verts*self._voxel_size + self._vol_origin

    # Get vertex colors
    rgb_vals = color_vol[verts_ind[:, 0], verts_ind[:, 1], verts_ind[:, 2]]
    colors_b = np.floor(rgb_vals / self._color_const)
    colors_g = np.floor((rgb_vals - colors_b*self._color_const) / 256)
    colors_r = rgb_vals - colors_b*self._color_const - colors_g*256
    colors = np.floor(np.asarray([colors_r, colors_g, colors_b])).T
    colors = colors.astype(np.uint8)

    pc = np.hstack([verts, colors])
    return pc 

def get_mesh(self):
    """Compute a mesh from the voxel volume using marching cubes.
    """
    tsdf_vol, color_vol = self.get_volume()

    # Marching cubes
    verts, faces, norms, vals = measure.marching_cubes_lewiner(tsdf_vol, level=0)
    verts_ind = np.round(verts).astype(int)
    verts = verts*self._voxel_size+self._vol_origin  # voxel grid coordinates to world coordinates

    # Get vertex colors
    rgb_vals = color_vol[verts_ind[:,0], verts_ind[:,1], verts_ind[:,2]]
    colors_b = np.floor(rgb_vals/self._color_const)
    colors_g = np.floor((rgb_vals-colors_b*self._color_const)/256)
    colors_r = rgb_vals-colors_b*self._color_const-colors_g*256
    colors = np.floor(np.asarray([colors_r,colors_g,colors_b])).T
    colors = colors.astype(np.uint8)
    return verts, faces, norms, colors 

def extract_surface(voxels, iso_level, dense=False):
    from skimage import measure
    verts, faces, normals, values = measure.marching_cubes_lewiner(voxels, iso_level)
    if dense:
        return augment_mesh(verts, faces)
    else:
        return verts 

def extract_mesh(vox):
    padded = np.pad(vox, 2, mode='constant', constant_values=0.)
    filtered = ndimage.filters.gaussian_filter(padded, 1., mode='constant', cval=0.)
    verts, faces, _, _ = measure.marching_cubes_lewiner(filtered, spacing=(0.1, 0.1, 0.1), gradient_direction='ascent')
    return dict(verts=verts.tolist(), faces=faces.tolist()) 

def write_obj(name,verts,faces,normals,values,affine=None,one=False):
	"""
	Write a .obj file for the output of marching cube algorithm.

	Parameters
	-----------
	name : str
		Ouput file name.
	verts : array
		Spatial coordinates for vertices as returned by skimage.measure.marching_cubes_lewiner().
	faces : array
		List of faces, referencing indices of verts as returned by skimage.measure.marching_cubes_lewiner().
	normals : array
		Normal direction of each vertex as returned by skimage.measure.marching_cubes_lewiner().
	affine : array,optional
		If given, vertices coordinates are affine transformed to create mesh with correct origin and size.
	one : bool
		Specify if faces values should start at 1 or at 0. Different visualization programs use different conventions.

	"""
	if (one) : faces=faces+1
	thefile = open(name,'w')
	if affine is not None:
		for item in verts:
			transformed = f(item[0],item[1],item[2],affine)
			thefile.write("v {0} {1} {2}\n".format(transformed[0],transformed[1],transformed[2]))
	else :
		for item in verts:
			thefile.write("v {0} {1} {2}\n".format(item[0],item[1],item[2]))
	for item in normals:
		thefile.write("vn {0} {1} {2}\n".format(item[0],item[1],item[2]))
	for item in faces:
		thefile.write("f {0}//{0} {1}//{1} {2}//{2}\n".format(item[0],item[1],item[2]))
	thefile.close() 

def image_to_surface(image_path, obj_file_path, voxel_size=1.0,
                     threshold=0, invert_axes=None, orientation="saggital",
                     step_size=1):
    """
    Saves the surface of an image as an .obj file

    :param image_path: str
    :param output_file: obj_file_path
    :param voxel_size: float (Default value = 1.0)
    :param threshold: float (Default value = 0)
    :param invert_axes: tuple (Default value = None)
    :param obj_file_path: 
    :param orientation:  (Default value = "saggital")
    :param step_size:  (Default value = 1)

    """

    image = brainio.load_any(image_path)

    image = reorient_image(image, invert_axes=invert_axes,
                           orientation=orientation)
    verts, faces, normals, values = \
        measure.marching_cubes_lewiner(image, threshold, step_size=step_size)

    # Scale to atlas spacing
    if voxel_size != 1.:
        verts = verts * voxel_size

    faces = faces + 1

    marching_cubes_to_obj((verts, faces, normals, values), obj_file_path) 

def surface_boundaries(self):
        """
        This function tries to find the isosurface using no interpolation to find the 
        correct positions of the surface to be able to shift to the interpolated one
        to the correct position

        Returns
        -------
        list of tuples 
            DESCRIPTION. [(mins[0],maxs[0]),(mins[1],maxs[1]),(mins[2],maxs[2])]

        """
        
        padding_x = self.padding[0]
        padding_y = self.padding[1]
        padding_z = self.padding[2]

        eigen_matrix = np.pad(self.V_matrix,
                              ((padding_x, padding_x), (padding_y, padding_y), 
                                (padding_z, padding_z)), "wrap")
        try:
            verts, faces, normals, values = measure.marching_cubes_lewiner(
                eigen_matrix, self.fermi)
            for ix in range(3):
                verts[:, ix] -= verts[:, ix].min()
                verts[:, ix] -= (verts[:, ix].max() - verts[:, ix].min()) / 2 #+self.origin[ix]
                verts[:, ix] *= self.dxyz[ix] 
            mins = verts.min(axis=0)
            maxs = verts.max(axis=0)
            
            return [(mins[0],maxs[0]),(mins[1],maxs[1]),(mins[2],maxs[2])]
        except:
            return None 

def mkoutersurf(image, radius, outfile):
    #radius information is currently ignored
    #it is a little tougher to deal with the morphology in python

    fill = nib.load( image )
    filld = fill.get_data()
    filld[filld==1] = 255

    gaussian = np.ones((2,2))*.25

    image_f = np.zeros((256,256,256))

    for slice in range(256):
        temp = filld[:,:,slice]
        image_f[:,:,slice] = convolve(temp, gaussian, 'same')

    image2 = np.zeros((256,256,256))
    image2[np.where(image_f <= 25)] = 0
    image2[np.where(image_f > 25)] = 255

    strel15 = generate_binary_structure(3, 1)

    BW2 = grey_closing(image2, structure=strel15)
    thresh = np.max(BW2)/2
    BW2[np.where(BW2 <= thresh)] = 0
    BW2[np.where(BW2 > thresh)] = 255

    v, f, _, _ = measure.marching_cubes_lewiner(BW2, 100)

    v2 = np.transpose(
             np.vstack( ( 128 - v[:,0],
                          v[:,2] - 128,
                          128 - v[:,1], )))
    
    write_surface(outfile, v2, f) 

def create_mesh_from_volume(volume, gradient_direction="ascent"):
    verts, faces, normals, values = measure.marching_cubes_lewiner(volume,
                                                                   0.0,
                                                                   spacing=(
                                                                       1.0, -1.0, 1.0),
                                                                   gradient_direction=gradient_direction)
    mesh = o3d.geometry.TriangleMesh()
    mesh.vertices = o3d.utility.Vector3dVector(verts)
    mesh.triangles = o3d.utility.Vector3iVector(faces)
    mesh.triangle_normals = o3d.utility.Vector3dVector(normals)
    return mesh 

def extract_mesh(input_val, params, indicators, input_holder, params_holder,
                 points_holder, sess, args):
  """Extracting meshes from an indicator function.

  Args:
    input_val: np.array, [1, height, width, channel], input image.
    params: tf.Operation, hyperplane parameter hook.
    indicators: tf.Operation, indicator hook.
    input_holder: tf.Placeholder, input image placeholder.
    params_holder: tf.Placeholder, hyperplane parameter placeholder.
    points_holder: tf.Placeholder, query point placeholder.
    sess: tf.Session, running sess.
    args: tf.app.flags.FLAGS, configurations.

  Returns:
    mesh: trimesh.Trimesh, the extracted mesh.
  """
  mesh_extractor = mise.MISE(64, 1, args.level_set)
  points = mesh_extractor.query()
  params_val = sess.run(params, {input_holder: input_val})

  while points.shape[0] != 0:
    orig_points = points
    points = points.astype(np.float32)
    points = (
        (np.expand_dims(points, axis=0) / mesh_extractor.resolution - 0.5) *
        args.vis_scale)
    n_points = points.shape[1]
    values = []
    for i in range(0, n_points, 100000):  # Add this to prevent OOM.
      value = sess.run(indicators, {
          params_holder: params_val,
          points_holder: points[:, i:i + 100000]
      })
      values.append(value)
    values = np.concatenate(values, axis=1)
    values = values[0, :, 0].astype(np.float64)
    mesh_extractor.update(orig_points, values)
    points = mesh_extractor.query()

  value_grid = mesh_extractor.to_dense()
  value_grid = np.pad(value_grid, 1, "constant", constant_values=-1e6)
  verts, faces, normals, unused_var = measure.marching_cubes_lewiner(
      value_grid, min(args.level_set,
                      value_grid.max() * 0.75))
  del normals
  verts -= 1
  verts /= np.array([
      value_grid.shape[0] - 3, value_grid.shape[1] - 3, value_grid.shape[2] - 3
  ],
                    dtype=np.float32)
  verts = args.vis_scale * (verts - 0.5)
  faces = np.stack([faces[..., 1], faces[..., 0], faces[..., 2]], axis=-1)
  return trimesh.Trimesh(vertices=verts, faces=faces) 

def mesh_region(region: bool, strel=None):
    r"""
    Creates a tri-mesh of the provided region using the marching cubes
    algorithm

    Parameters
    ----------
    im : ND-array
        A boolean image with ``True`` values indicating the region of interest

    strel : ND-array
        The structuring element to use when blurring the region.  The blur is
        perfomed using a simple convolution filter.  The point is to create a
        greyscale region to allow the marching cubes algorithm some freedom
        to conform the mesh to the surface.  As the size of ``strel`` increases
        the region will become increasingly blurred and inaccurate. The default
        is a spherical element with a radius of 1.

    Returns
    -------
    mesh : tuple
        A named-tuple containing ``faces``, ``verts``, ``norm``, and ``val``
        as returned by ``scikit-image.measure.marching_cubes`` function.

    """
    im = region
    if im.ndim != im.squeeze().ndim:
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
                      ' unexpected behavior.')
    if strel is None:
        if region.ndim == 3:
            strel = ball(1)
        if region.ndim == 2:
            strel = disk(1)
    pad_width = np.amax(strel.shape)
    if im.ndim == 3:
        padded_mask = np.pad(im, pad_width=pad_width, mode='constant')
        padded_mask = spim.convolve(padded_mask * 1.0,
                                    weights=strel) / np.sum(strel)
    else:
        padded_mask = np.reshape(im, (1,) + im.shape)
        padded_mask = np.pad(padded_mask, pad_width=pad_width, mode='constant')
    verts, faces, norm, val = marching_cubes_lewiner(padded_mask)
    result = namedtuple('mesh', ('verts', 'faces', 'norm', 'val'))
    result.verts = verts - pad_width
    result.faces = faces
    result.norm = norm
    result.val = val
    return result 

def convert_voxel(img_data, affine = None, threshold = None, data_mask = None, absthreshold = None):
	"""
	Converts a voxel image to a surface including outputs voxel values to paint vertex surface.
	
	Parameters
	----------
	img_data : array
		image array
	affine : array
		 affine [4x4] to convert vertices values to native space (Default = None)
	data_mask : array
		use a mask to create a surface backbone (Default = None)
	threshold : float
		threshold for output of voxels (Default = None)
	absthreshold : float
		threshold for output of abs(voxels) (Default = None)
		
	Returns
	-------
		v : array
			vertices
		f : array
			faces
		values : array
			scalar values
	
	"""
	try:
		from skimage import measure
	except:
		print("Error skimage is required")
		quit()

	if threshold is not None:
		print("Zeroing data less than threshold = %1.2f" % threshold)
		img_data[img_data<threshold] = 0
	if absthreshold is not None:
		print("Zeroing absolute values less than threshold = %1.2f" % absthreshold)
		img_data[np.abs(img_data)<absthreshold] = 0
	if data_mask is not None:
		print("Including mask")
		data_mask *= .1
		data_mask[img_data!=0] = img_data[img_data!=0]
		del img_data
		img_data = np.copy(data_mask)
	try:
		v, f, _, values = measure.marching_cubes_lewiner(img_data)
		if affine is not None:
			print("Applying affine transformation")
			v = nib.affines.apply_affine(affine,v)
	except:
		print("No voxels above threshold")
		v = f = values = []
	return v, f, values

# Check if okay to remove