python类marching_cubes()的实例源码

def plot_3D(img, threshold=-400):
    verts, faces = measure.marching_cubes(img, threshold)

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

    mesh = Poly3DCollection(verts[faces], alpha=0.1)
    face_color = [0.5, 0.5, 1]
    mesh.set_facecolor(face_color)
    ax.add_collection3d(mesh)

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

    plt.show()

def generate_mesh(image, isovalue=0, channel=0):
    """
    Creates and returns a Mesh object
    :param image: an AICSImage object
    :param isovalue: The value that is used to pick the isosurface returned by the marching cubes algorithm
                     For more info: https://www.youtube.com/watch?v=5fNbCFjqWao @ 40:00 mins
    :param channel: The channel in the image that is used to extract the isosurface
    :return: A Mesh object
    """
    if not isinstance(image, AICSImage):
        raise ValueError("Meshes can only be generated with AICSImage objects!")
    if channel >= image.size_c:
        raise IndexError("Channel provided for mesh generation is out of bounds for image data!")
    image_stack = image.get_image_data("ZYX", C=channel)
    # Use marching cubes to obtain the surface mesh of the membrane wall
    verts, faces, normals, values = measure.marching_cubes(image_stack, isovalue, allow_degenerate=False)
    return Mesh(verts, faces, normals, values)

def marching_cubes(field,iso=0.5):
    try:
        from skimage.measure import marching_cubes
        surface_points, surface_triangles = marching_cubes(density_field,iso)

    except ImportError:
        print "Please try to install SciKit-Image!"

        from mayavi import mlab
        from mayavi.mlab import contour3d

        mlab.clf()
        surface = mlab.contour3d(field,contours=[iso])

        my_actor=surface.actor.actors[0] 
        poly_data_object=my_actor.mapper.input 
        surface_points = (np.array(poly_data_object.points) - np.array([abs(grid_points/2.),abs(grid_points/2.),abs(grid_points/2.)])[np.newaxis,:])*(grid_max/abs(grid_points/2.))
        surface_triangles = poly_data_object.polys.data.to_array().reshape([-1,4]) 
        surface_triangles = surface_triangles[:,1:]

    return surface_points, surface_triangles

def plot_3d_cubic(image):
    '''
        plot the 3D cubic
    :param image:   image saved as npy file path
    :return:
    '''
    from skimage import measure, morphology
    from mpl_toolkits.mplot3d.art3d import Poly3DCollection
    image = np.load(image)
    verts, faces = measure.marching_cubes(image,0)
    fig = plt.figure(figsize=(40, 40))
    ax = fig.add_subplot(111, projection='3d')
    # Fancy indexing: `verts[faces]` to generate a collection of triangles
    mesh = Poly3DCollection(verts[faces], alpha=0.1)
    face_color = [0.5, 0.5, 1]
    mesh.set_facecolor(face_color)
    ax.add_collection3d(mesh)
    ax.set_xlim(0, image.shape[0])
    ax.set_ylim(0, image.shape[1])
    ax.set_zlim(0, image.shape[2])
    plt.show()

# LUNA2016 data prepare ,first step: truncate HU to -1000 to 400

def plot_3d_cubic(image):
    '''
        plot the 3D cubic
    :param image:   image saved as npy file path
    :return:
    '''
    from skimage import measure, morphology
    from mpl_toolkits.mplot3d.art3d import Poly3DCollection
    image = np.load(image)
    verts, faces = measure.marching_cubes(image,0)
    fig = plt.figure(figsize=(40, 40))
    ax = fig.add_subplot(111, projection='3d')
    # Fancy indexing: `verts[faces]` to generate a collection of triangles
    mesh = Poly3DCollection(verts[faces], alpha=0.1)
    face_color = [0.5, 0.5, 1]
    mesh.set_facecolor(face_color)
    ax.add_collection3d(mesh)
    ax.set_xlim(0, image.shape[0])
    ax.set_ylim(0, image.shape[1])
    ax.set_zlim(0, image.shape[2])
    plt.show()

# LUNA2016 data prepare ,first step: truncate HU to -1000 to 400

def plot_3d_cubic(image):
    '''
        plot the 3D cubic
    :param image:   image saved as npy file path
    :return:
    '''
    from skimage import measure, morphology
    from mpl_toolkits.mplot3d.art3d import Poly3DCollection
    image = np.load(image)
    verts, faces = measure.marching_cubes(image,0)
    fig = plt.figure(figsize=(40, 40))
    ax = fig.add_subplot(111, projection='3d')
    # Fancy indexing: `verts[faces]` to generate a collection of triangles
    mesh = Poly3DCollection(verts[faces], alpha=0.1)
    face_color = [0.5, 0.5, 1]
    mesh.set_facecolor(face_color)
    ax.add_collection3d(mesh)
    ax.set_xlim(0, image.shape[0])
    ax.set_ylim(0, image.shape[1])
    ax.set_zlim(0, image.shape[2])
    plt.show()

# LUNA2016 data prepare ,first step: truncate HU to -1000 to 400

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)

    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 three_d_print(self):
        """This will produce a 3d printable stl based on self.volume_data. It is to be used for the final "print" button, and needs to be fed high quality data."""

        name = raw_input('What should the filename be?') + '.stl'


        verts, faces = measure.marching_cubes(self.volume_data, 0)   #Marching Cubes algorithm

        solid = mesh.Mesh(np.zeros(faces.shape[0], dtype=mesh.Mesh.dtype))
        for i, f in enumerate(faces):
            for j in range(3):
                solid.vectors[i][j] = verts[f[j],:]


        solid.save(name)

#Here be UI

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 estimate_surface_area(self):
        """
        Estimate the surface area by summing the areas of a trianglation
        of the nodules surface in 3d. Returned units are mm^2.
        """
        mask = self.get_boolean_mask()
        mask = np.pad(mask, [(1,1), (1,1), (1,1)], 'constant') # Cap the ends.
        dist = dtrans(mask) - dtrans(~mask)

        rxy  = self.scan.pixel_spacing
        rz   = self.scan.slice_thickness
        verts, faces, _, _ = marching_cubes(dist, 0, spacing=(rxy, rxy, rz))
        return mesh_surface_area(verts, faces)

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

def implicit_surface(density_field,size,resolution,iso=0.5):
    import numpy as np
    from scipy.cluster.vq                       import kmeans, vq
    from openalea.container import array_dict

    from skimage.measure import marching_cubes
    surface_points, surface_triangles = marching_cubes(density_field,iso)

    surface_points = (np.array(surface_points))*(size*resolution/np.array(density_field.shape)) - size*resolution/2.

    points_ids = np.arange(len(surface_points))
    points_to_delete = []
    for p,point in enumerate(surface_points):
        matching_points = np.sort(np.where(vq(surface_points,np.array([point]))[1] == 0)[0])
        if len(matching_points) > 1:
            points_to_fuse = matching_points[1:]
            for m_p in points_to_fuse:
                surface_triangles[np.where(surface_triangles==m_p)] = matching_points[0]
                points_to_delete.append(m_p)

    points_to_delete = np.unique(points_to_delete)
    print len(points_to_delete),"points deleted"
    surface_points = np.delete(surface_points,points_to_delete,0)
    points_ids = np.delete(points_ids,points_to_delete,0)
    surface_triangles = array_dict(np.arange(len(surface_points)),points_ids).values(surface_triangles)

    for p,point in enumerate(surface_points):
        matching_points = np.where(vq(surface_points,np.array([point]))[1] == 0)[0]
        if len(matching_points) > 1:
            print p,point
            raw_input()

    triangles_to_delete = []
    for t,triangle in enumerate(surface_triangles):
        if len(np.unique(triangle)) < 3:
            triangles_to_delete.append(t)
        # elif triangle.max() >= len(surface_points):
        #     triangles_to_delete.append(t)
    surface_triangles = np.delete(surface_triangles,triangles_to_delete,0)

    return surface_points, surface_triangles