implicit_surfaces.py 文件源码

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