Python spatial.Delaunay类代码示例

def show_hull(cname, ccol):
    ccoords = coords[names==cname]

    hull = ConvexHull(ccoords)

    xs = ccoords.ix[:,0]
    ys = ccoords.ix[:,1]
    zs = ccoords.ix[:,2]
    #ax.scatter(xs, ys, zs, c=ccol, marker='o')

    for simplex in hull.simplices:
        s = ccoords.irow(simplex)
        #print s
        sx = list(s.ix[:,0])
        sy = list(s.ix[:,1])
        sz = list(s.ix[:,2])
        sx.append(sx[0])
        sy.append(sy[0])
        sz.append(sz[0])
        ax.plot(sx, sy, sz, ccol, alpha=0.2)

    hulld = Delaunay(ccoords.irow(hull.vertices))
    hulld.find_simplex(coords)
    hcol = ['grey' if x<0 else 'green' for x in hulld.find_simplex(coords)]
    hxs = coords.ix[:,0]
    hys = coords.ix[:,1]
    hzs = coords.ix[:,2]
    ax.scatter(hxs, hys, hzs, c=hcol, marker='o', alpha=0.2)

def zonotope_quadrature_rule(avmap, N, NX=10000):
    """
    Description of zonotope_quadrature_rule

    Arguments:
        vert:
        W1:
        N:
        NX: (default=10000)
    Outputs:
        points:
        weights:
    """
    
    vert = avmap.domain.vertY
    W1 = avmap.domain.subspaces.W1    
    
    # number of dimensions
    m, n = W1.shape

    # points
    y = np.vstack((vert, maximin_design(vert, N)))
    T = Delaunay(y)
    c = []
    for t in T.simplices:
        c.append(np.mean(T.points[t], axis=0))
    points = np.array(c)

    # approximate weights
    Y_samples = np.dot(np.random.uniform(-1.0, 1.0, size=(NX,m)), W1)
    I = T.find_simplex(Y_samples)
    weights = np.zeros((T.nsimplex, 1))
    for i in range(T.nsimplex):
        weights[i] = np.sum(I==i) / float(NX)
    return points.reshape((T.nsimplex,n)), weights.reshape((T.nsimplex,1))

    def _simplex_interp(self, teff, logg, fe, model_indecies):
        """
        Perform barycentric interpolation on simplecies

        Args:
            teff (float): effective temperature
            logg (float): surface gravity
            fe (float): metalicity
            model_indecies (array): models to use in the interpolation

        Returns:
            array: spectrum at teff, logg, fe

        """

        ndim = 3
        model_table = np.array(self.model_table['teff logg fe'.split()])
        points = model_table[model_indecies]
        tri = Delaunay(points) # Delaunay triangulation
        p = np.array([teff, logg, fe]) # cartesian coordinates

        simplex = tri.find_simplex(p)
        r = tri.transform[simplex,ndim,:]
        Tinv = tri.transform[simplex,:ndim,:ndim]
        c = Tinv.dot(p - r)
        c = np.hstack([c,1.0-c.sum()]) # barycentric coordinates

        # (3,many array)
        model_indecies_interp = model_indecies[tri.simplices[simplex]]
        model_spectra = self.model_spectra[model_indecies_interp,:]
        flux = np.dot(c,model_spectra)
        if simplex==-1:
            return np.ones_like(flux)
        return flux

    def __init__(self, target_image_dimensions, target_vertices=None):
        self.target_image_dimensions = target_image_dimensions
        # if target vertices are not set, assume we want to fill the entire targetImage
        if not target_vertices:
            rows, columns = target_image_dimensions[:2]
            self.target_vertices = np.array([[0, 0],
                                             [0, columns],
                                             [rows, 0],
                                             [rows, columns]])
        else:
            self.target_vertices = target_vertices

        r0, c0 = np.round(np.mean(self.target_vertices, axis=0)).astype(np.int32)
        self.target_vertices = np.array(sorted(self.target_vertices, key=lambda (r, c): np.arctan2(r0 - r, c0 - c)))
        # get barycentric coordinates and corresponding points in target (new) image
        nys = np.arange(self.target_image_dimensions[0])
        nxs = np.arange(self.target_image_dimensions[1])
        image_pixels = np.transpose([np.repeat(nys, len(nxs)), np.tile(nxs, len(nys))])

        triangles = Delaunay(self.target_vertices)
        memberships = triangles.find_simplex(image_pixels)  # returns the triangle that each pixel is a member of
        Ts = triangles.transform[memberships, :2]  # transformation matrices
        prs = image_pixels - triangles.transform[memberships, 2]  # intermediate transformation

        barycentric_coordinates = np.array([Ts[i].dot(pr) for i, pr in enumerate(prs)])
        barycentric_coordinates = np.hstack((barycentric_coordinates,
                                            1 - np.sum(barycentric_coordinates, axis=1, keepdims=True)))

        target_vertices_indices = triangles.simplices[memberships]

        self.targetRow = image_pixels[:, 0]
        self.targetCol = image_pixels[:, 1]
        self.barycentric = barycentric_coordinates
        self.indices = target_vertices_indices

    def draw_dimension(self, dimNumber):
        '''
        Method to draw the points on the axes using the current dimension number
        '''

        self.ax.cla()

        dim0 = self.combinations[self.dimNumber][0]
        dim1 = self.combinations[self.dimNumber][1]

        self.ax.plot(self.points[:,dim0][self.inCluster], self.points[:, dim1][self.inCluster], 'g.')
        
        self.ax.plot(self.points[:, dim0][self.outsideCluster], self.points[:,dim1][self.outsideCluster], marker='.', color='0.8', linestyle='None')

        self.ax.set_xlabel('Dimension {}'.format(dim0))
        self.ax.set_ylabel('Dimension {}'.format(dim1))

        
        plt.title('press c to cut, < or > to switch dimensions')
        self.fig.canvas.draw()
        
    
        ''' Method to take the current points from mouse input, 
        convert them to a convex hull, and then update the inCluster and 
outsideCluster attributes based on the points that fall within the hull'''

        if not isinstance(hull, Delaunay):
            hull=Delaunay(hull)

        self.inCluster = hull.find_simplex(points)>=0
        self.outsideCluster=np.logical_not(self.inCluster)

	def triangulate(self, size):
		shape = self.mImg.shape
		spacing = max(shape) / size
		sigma = spacing / 4
		coords = self._get_distributed_points(spacing, sigma, include_corners = True)
		
		tri = Delaunay(coords)
		im_pts = self.get_xy_features()
		# pt_tri_membership becomes a map which is the same size as the
		# original image (first two dimensions only). each element contains
		# the triangle ID of that point in the source image
		pt_tri_membership = tri.find_simplex(im_pts.astype(dtype = np.double))
		pt_tri_membership.resize(shape[0], shape[1])
		num_tri = np.max(pt_tri_membership)

		tri_map = np.copy(self.mImg) 
		# replace elements of each triangle with the mean value of the color
		# channels from the original image
		for tri in range(num_tri + 1):
			this_tri = pt_tri_membership == tri
			if not np.any(this_tri):
				continue
			for col in range(shape[2]):
				tri_map[this_tri, col] = np.mean(self.mImg[this_tri, col])
		return tri_map

 def compute(self):
     """ABC API"""
     Delaunay.__init__(self,
                       self._points,
                       self._furthest_site,
                       self._incremental,
                       self._qhull_options)

 def compute(self):
     """ABC API"""
     self.id = "D({},{})".format(self._furthest_site, self._qhull_options)
     Delaunay.__init__(self,
                       self._points,
                       self._furthest_site,
                       self._incremental,
                       self._qhull_options)

def inside_convex_poly(pts):
    """Returns a function that checks if inputs are inside the convex hull of polyhedron defined by pts

    Alternative method to check is to get faces of the convex hull, then check if each normal is pointed away from each point.
    As it turns out, this is vastly slower than using qhull's find_simplex, even though the simplex is not needed.
    """
    tri = Delaunay(pts)
    return lambda x: tri.find_simplex(x) != -1

def chromaticity_diagram_visual(
        samples=256,
        cmfs='CIE 1931 2 Degree Standard Observer',
        transformation='CIE 1931',
        parent=None):
    """
    Creates a chromaticity diagram visual based on
    :class:`colour_analysis.visuals.Primitive` class.

    Parameters
    ----------
    samples : int, optional
        Inner samples count used to construct the chromaticity diagram
        triangulation.
    cmfs : unicode, optional
        Standard observer colour matching functions used for the chromaticity
        diagram boundaries.
    transformation : unicode, optional
        {'CIE 1931', 'CIE 1960 UCS', 'CIE 1976 UCS'}

        Chromaticity diagram transformation.
    parent : Node, optional
        Parent of the chromaticity diagram in the `SceneGraph`.

    Returns
    -------
    Primitive
        Chromaticity diagram visual.
    """

    cmfs = get_cmfs(cmfs)

    illuminant = DEFAULT_PLOTTING_ILLUMINANT

    XYZ_to_ij = (
        CHROMATICITY_DIAGRAM_TRANSFORMATIONS[transformation]['XYZ_to_ij'])
    ij_to_XYZ = (
        CHROMATICITY_DIAGRAM_TRANSFORMATIONS[transformation]['ij_to_XYZ'])

    ij_c = XYZ_to_ij(cmfs.values, illuminant)

    triangulation = Delaunay(ij_c, qhull_options='QJ')
    samples = np.linspace(0, 1, samples)
    ii, jj = np.meshgrid(samples, samples)
    ij = tstack((ii, jj))
    ij = np.vstack((ij_c, ij[triangulation.find_simplex(ij) > 0]))

    ij_p = np.hstack((ij, np.full((ij.shape[0], 1), 0)))
    triangulation = Delaunay(ij, qhull_options='QJ')
    RGB = normalise(XYZ_to_sRGB(ij_to_XYZ(ij, illuminant), illuminant),
                    axis=-1)

    diagram = Primitive(vertices=ij_p,
                        faces=triangulation.simplices,
                        vertex_colours=RGB,
                        parent=parent)

    return diagram

def hadrian_min(vectorized_f, xbnds, ybnds, xtol, ytol, swarm=8, mx_iters=5,
        inc=False):
    """
    hadrian_min is a stochastic, hill climbing minimization algorithm.  It
    uses a stratified sampling technique (Latin Hypercube) to get good
    coverage of potential new points.  It also uses vectorized function
    evaluations to drive concurrent function evaluations.

    It is named after the Roman Emperor Hadrian, the most famous Latin hill
    mountain climber of ancient times.
    """
    assert xbnds[1] > xbnds[0]
    assert ybnds[1] > ybnds[0]
    assert xtol > 0
    assert ytol > 0
    # simplexes are simplex indexes
    # vertexes are vertex indexes
    # points are spatial coordinates

    bnds = np.vstack((xbnds, ybnds))
    points = latin_sample(np.vstack((xbnds, ybnds)), swarm)
    z = vectorized_f(points)

    # exclude corners from possibilities, but add them to the triangulation
    # this bounds the domain, but ensures they don't get picked
    points = np.append(points, list(product(xbnds, ybnds)), axis=0)
    z = np.append(z, 4*[z.max()])

    tri = Delaunay(points, incremental=inc)
    del points

    for step in range(1, mx_iters+1):
        i, vertexes = get_minimum_neighbors(tri, z)
        disp = tri.points[i] - tri.points[np.unique(vertexes)]
        disp /= np.array([xtol, ytol])
        err = err_mean(disp)

        if err < 1.:
            return tri.points[i], z[i], tri.points, z, 1

        tri_points = tri.points[vertexes]
        bnds = np.cumsum(area(tri_points))
        bnds /= bnds[-1]
        indx = np.searchsorted(bnds, np.random.rand(swarm))

        new_pts = sample_triangle(tri_points[indx])
        new_z = vectorized_f(new_pts)

        if inc:
            tri.add_points(new_pts)
        else:
            # make a new triangulation if I can't append points
            points = np.append(tri.points, new_pts, axis=0)
            tri = Delaunay(points)
        z = np.append(z, new_z)

    return None, None, tri.points, z, step

def interp_weights(xyz, uvw,d=3):

    tri = Delaunay(xyz)
    simplex = tri.find_simplex(uvw)
    vertices = np.take(tri.simplices, simplex, axis=0)
    temp = np.take(tri.transform, simplex, axis=0)
    delta = uvw - temp[:, d]
    bary = np.einsum('njk,nk->nj', temp[:, :d, :], delta)
    return vertices, np.hstack((bary, 1 - bary.sum(axis=1, keepdims=True)))

def check_polygon_membership(cluster_boundary, points_to_check):
    cluster_boundary_points = [[x[0], x[1]] for x in cluster_boundary]
    hull = cluster_boundary_points
    if not isinstance(hull, Delaunay):
        hull = Delaunay(hull)
    gps_coords_to_check = [[x[0], x[1]] for x in points_to_check]
    point_in_cluster = []
    for coord in gps_coords_to_check:
        point_in_cluster.append(1 if hull.find_simplex(coord) >= 0 else 0)
    return point_in_cluster

    def cut_cluster(self, points, hull):
        ''' Method to take the current points from mouse input, 
        convert them to a convex hull, and then update the inCluster and 
outsideCluster attributes based on the points that fall within the hull'''

        if not isinstance(hull, Delaunay):
            hull=Delaunay(hull)

        self.inCluster = hull.find_simplex(points)>=0
        self.outsideCluster=np.logical_not(self.inCluster)

 def inside_hull(self, labels):
     """This method checks whether a requested set of labels is inside the
     convex hull of the training data.  Returns a bool.  This method relies
     on Delauynay Triangulation and is thus exceedlingly slow for large
     dimensionality.
     """
     L = flatten_struct(self.training_labels, use_labels=self.used_labels)
     l = flatten_struct(labels, use_labels=self.used_labels)
     hull = Delaunay(L.T)
     return hull.find_simplex(l.T) >= 0