Python scipy.spatial 模块，ConvexHull() 实例源码

def convex_hull(self):
        """Return a 3D mesh that represents the convex hull of the mesh.
        """
        hull = ss.ConvexHull(self.vertices_)
        hull_tris = hull.simplices
        if self.normals_ is None:
            cvh_mesh = Mesh3D(self.vertices_.copy(), hull_tris.copy(), center_of_mass=self.center_of_mass_)
        else:
            cvh_mesh = Mesh3D(self.vertices_.copy(), hull_tris.copy(), normals=self.normals_.copy(), center_of_mass=self.center_of_mass_)
        cvh_mesh.remove_unreferenced_vertices()
        return cvh_mesh

def determine_convex_hull(formation_en_grid):
    """
    Wraps the scipy.spatial ConvexHull algo for our purposes.
    For now only for 2D phase diagrams
    Adds the points [1.0, 0.0] and [0.0, 1.0], because in material science these
    are always there.

    :params: formation_en_grid: list of points in phase space [[x, formation_energy]]
    :returns: a hul datatype
    """
    import numpy as np
    from scipy.spatial import ConvexHull

    # TODO multi d
    # check if endpoints are in
    if [1.0, 0.0] not in formation_en_grid:
        formation_en_grid.append([1.0, 0.0])
    if [0.0, 0.0] not in formation_en_grid:
        formation_en_grid.append([0.0, 0.0])

    points = np.array(formation_en_grid)
    hull = ConvexHull(points)

    return hull

def compute_vc(self, points):
        """compute the Voronoi cell"""
        #compute S and V
        S = ConvexHull(points).area
        V = ConvexHull(points).volume

        #voronoi cell
        eta = S ** 3 / (36 * np.pi * V ** 2)

        return eta

def convex_hull_area(pts):
    """
    Calculates the surface area from a given point cloud using simplices of
    its convex hull. For the estimation of the synapse contact area, divide by
    a factor of two, in order to get the area of only one face (we assume that
    the contact site is sufficiently thin represented by the points).
    :param pts: np.array of coordinates in nm (scaled)
    :return: Area of the point cloud (nm^2)
    """
    if len(pts) < 4:
        return 0
    area = 0
    try:
        ch = ConvexHull(pts)
        triangles = ch.points[ch.simplices]
        for triangle in triangles:
            area += poly_area(triangle)
    except QhullError as e:
        # warnings.warn("%s encountered during calculation of convex hull "
        #               "area with %d points. Returning 0 nm^2." %
        #               (e, len(pts)), RuntimeWarning)
        pass
    return area

def fit_convex_hull(points):
    """ Creates a feasible set by taking a convex hull of the points given. Returns P = { x : Ax >= b }

    Args:
        points (list): Set of numpy points.

    Returns:
        A (numpy): constraint matrix
        b (numpy): constraint vector
    """
    hull = ConvexHull(points)
    m,n = hull.equations.shape
    A = -1 * hull.equations[:,0:n-1]
    b = hull.equations[:,n-1]
    return np.mat(A), np.mat(b).T

def tri_normals(self, align_to_hull=False):
        """Returns a list of the triangle normals.

        Parameters
        ----------
        align_to_hull : bool
            If true, we re-orient the normals to point outward from
            the mesh by using the convex hull.

        Returns
        -------
        :obj:`numpy.ndarray` of float
            A #triangles by 3 array of floats, where each 3-ndarray
            represents the 3D normal vector of the corresponding triangle.
        """
        # compute normals
        v0 = self.vertices_[self.triangles_[:,0],:]
        v1 = self.vertices_[self.triangles_[:,1],:]
        v2 = self.vertices_[self.triangles_[:,2],:]
        n = np.cross(v1 - v0, v2 - v0)
        normals = n / np.tile(np.linalg.norm(n, axis=1)[:,np.newaxis], [1,3])

        # reverse normal based on alignment with convex hull
        if align_to_hull:
            tri_centers = self.tri_centers()
            hull = ss.ConvexHull(tri_centers)
            hull_tris = hull.simplices
            hull_vertex_ind = hull_tris[0][0]
            hull_vertex = tri_centers[hull_vertex_ind]
            hull_vertex_normal = normals[hull_vertex_ind]
            v = hull_vertex.reshape([1,3])
            n = hull_vertex_normal
            ip = (tri_centers - np.tile(hull_vertex,
                  [tri_centers.shape[0], 1])).dot(n)
            if ip[0] > 0:
                normals = -normals
        return normals

def _rubberband(bands, intensities, num_ranges):
    '''Basic rubberband method,
    from p.77 of "IR and Raman Spectroscopy" (OPUS manual)'''
    # create n ranges of equal size in the spectrum
    range_size = len(intensities) // num_ranges
    y = intensities[:range_size * num_ranges].reshape((num_ranges, range_size))
    # find the smallest intensity point in each range
    idx = np.arange(num_ranges) * range_size + np.argmin(y, axis=1)
    # add in the start and end points as well, to avoid weird edge effects
    if idx[0] != 0:
        idx = np.append(0, idx)
    if idx[-1] != len(intensities) - 1:
        idx = np.append(idx, len(intensities) - 1)
    baseline_pts = np.column_stack((bands[idx], intensities[idx]))
    # wrap a rubber band around the baseline points
    hull = ConvexHull(baseline_pts)
    hidx = idx[hull.vertices]
    # take only the bottom side of the hull
    left = np.argmin(bands[hidx])
    right = np.argmax(bands[hidx])
    mask = np.ones(len(hidx), dtype=bool)
    for i in range(len(hidx)):
        if i > right and (i < left or right > left):
            mask[i] = False
        elif i < left and i < right:
            mask[i] = False
    hidx = hidx[mask]
    hidx = hidx[np.argsort(bands[hidx])]
    # interpolate a baseline
    return np.interp(bands, bands[hidx], intensities[hidx])

def compute_polyhedron_shadow(n, light_poly, ground_poly, inf_dist=1000):
    light_hull = ConvexHull(light_poly)
    ground_hull = ConvexHull(ground_poly)
    light_vertices = [light_poly[i] for i in light_hull.vertices]
    ground_vertices = [ground_poly[i] for i in ground_hull.vertices]
    rays = [gv - lv for gv in ground_vertices for lv in light_vertices]
    return ground_vertices + [v + inf_dist * r
                              for v in ground_vertices for r in rays]

def polygons(latitudes, longitudes, clusters, maptype=MAPTYPE):
    """Plot clusters of points on map, including them in a polygon defining their convex hull.

    :param pandas.Series latitudes: series of sample latitudes
    :param pandas.Series longitudes: series of sample longitudes
    :param pandas.Series clusters: marker clusters, as integers
    :param string maptype: type of maps, see GoogleStaticMapsAPI docs for more info

    :return: None
    """
    width = SCALE * MAX_SIZE
    img, pixels = background_and_pixels(latitudes, longitudes, MAX_SIZE, maptype)

    polygons = []
    for c in clusters.unique():
        in_polygon = clusters == c
        if in_polygon.sum() < 3:
            print '[WARN] Cannot draw polygon for cluster {} - only {} samples.'.format(c, in_polygon.sum())
            continue
        cluster_pixels = pixels.loc[clusters == c]
        polygons.append(Polygon(cluster_pixels.iloc[ConvexHull(cluster_pixels).vertices], closed=True))

    plt.figure(figsize=(10, 10))
    ax = plt.subplot(111)
    plt.imshow(np.array(img))                                               # Background map
    p = PatchCollection(polygons, cmap='jet', alpha=0.15)                   # Collection of polygons
    p.set_array(clusters.unique())
    ax.add_collection(p)
    plt.scatter(                                                            # Scatter plot
        pixels['x_pixel'],
        pixels['y_pixel'],
        c=clusters,
        s=width / 40,
        linewidth=0,
        alpha=0.25,
    )
    plt.gca().invert_yaxis()                                                # Origin of map is upper left
    plt.axis([0, width, width, 0])                                          # Remove margin
    plt.axis('off')
    plt.tight_layout()
    plt.show()

def get_edges(a, convex=False):
    a = checkma(a)
    #Need to deal with RGB images here
    #Need to be careful, probably want to take minimum value from masks
    if a.ndim == 3:
        #Assume that the same mask applies to each band
        #Only create edges for one band
        b = a[:,:,0]
        #Could convert to HSL and do edges for L channel
        #b = a[:,:,0].mask + a[:,:,1].mask + a[:,:,2].mask
    else:
        b = a

    #Compute edges along both axes, need both to handle undercuts
    #These are inclusive, indices indicate position of unmasked data
    edges0 = np.ma.notmasked_edges(b, axis=0)
    edges1 = np.ma.notmasked_edges(b, axis=1)
    edges = np.array([np.concatenate([edges0[0][0], edges0[1][0], edges1[1][0], edges1[0][0]]), np.concatenate([edges0[0][1], edges0[1][1], edges1[1][1], edges1[0][1]])])

    #This is a rough outline - needs testing
    if convex:
        from scipy.spatial import ConvexHull
        #hull = ConvexHull(edges.T)
        #edges = edges.T[hull.simplices]
        #This is in scipy v0.14
        #edges0 = edges1 = hull.vertices

    return edges0, edges1, edges

def convex_hull(self, model):
        train = np.array(model.data)
        global t1,t2
        t1,t2 = [],[]

        for i,example in enumerate(train):
            if example[2]==1:
                t1.append(example[0:2])
            else:
                t2.append(example[0:2])
        global ch1,ch2,has_margin
        if (has_margin):
            self.ax.lines.pop(len(t1)+len(t2)-1)
        for i in range(ch1+ch2):
            self.ax.lines.pop(len(t1)+len(t2)-1)
        if len(t1) > 2:
            X1 = np.matrix(np.array(t1))
            hull1 = ConvexHull(X1);    
            for simplex in hull1.simplices:
                self.ax.plot(X1[simplex,0],X1[simplex,1], color='grey',linewidth=1)
            ch1 = len(hull1.simplices)

        if len(t2) > 2:
            X2 = np.matrix(np.array(t2))
            hull2 = ConvexHull(X2);
            for simplex in hull2.simplices:
                self.ax.plot(X2[simplex,0],X2[simplex,1], color='0.25')
            #self.ax.plot(cx2,cy2,'+',ms=10, color="blue")
            ch2 = len(hull2.simplices)
        global isFitted
        #print len(train),"\t"
        #if isFitted:
        #    self.draw_distance(model,model.clf.gamma)
        self.set_dim()

def build_from_hull(self):
        """
        Test oriented bounding box algorithm using convex hull points.

        Raises:
            None

        Returns:
            EigenVectors(OpenMaya.MVector)
            CenterPoint(OpenMaya.MVector)
            BoundingExtents(OpenMaya.MVector)
        """
        npPointList = [[self.points[i].x, self.points[i].y, self.points[i].z]
                       for i in xrange(self.points.length())]

        try:
            hull = ConvexHull(npPointList)
        except NameError:
            raise RuntimeError(
                "From hull method unavailable because"
                " scipy cannot be imported."
                "Please install it if you need it.")

        indices = hull.simplices
        vertices = npPointList[indices]
        hullPoints = list(vertices.flatten())
        hullTriPoints = list(indices.flatten())

        hullArray = OpenMaya.MVectorArray()
        for ind in xrange(0, len(hullPoints), 3):
            hullArray.append(
                OpenMaya.MVector(
                    hullPoints[ind], hullPoints[ind+1], hullPoints[ind+2]))

        triPoints = OpenMaya.MIntArray()
        for tri in xrange(len(hullTriPoints)):
            triPoints.append(tri)

        return self.build_from_triangles(points=hullArray, triangles=triPoints)

def generate_one_example(n_nodes):
    points = numpy.random.rand(n_nodes, 2)
    hull = ConvexHull(points)  # scipy.spatial.ConvexHull will generate points in CCW order
    v = hull.vertices
    # v = numpy.roll(v, -list(v).index(numpy.min(v)))  # start from the smallest indice
    return points, v + 1

def oriented_bounds_2D(points):
    '''
    Find an oriented bounding box for a set of 2D points.

    Arguments
    ----------
    points: (n,2) float, 2D points

    Returns
    ----------
    transform: (3,3) float, homogenous 2D transformation matrix to move the input set of 
               points to the FIRST QUADRANT, so no value is negative. 
    rectangle: (2,) float, size of extents once input points are transformed by transform
    '''
    c = ConvexHull(np.asanyarray(points))
    # (n,2,3) line segments
    hull = c.points[c.simplices] 
    # (3,n) points on the hull to check against
    dot_test = c.points[c.vertices].reshape((-1,2)).T
    edge_vectors = unitize(np.diff(hull, axis=1).reshape((-1,2)))
    perp_vectors = np.fliplr(edge_vectors) * [-1.0,1.0]
    bounds = np.zeros((len(edge_vectors), 4))
    for i, edge, perp in zip(range(len(edge_vectors)),
                             edge_vectors, 
                             perp_vectors):
        x = np.dot(edge, dot_test)
        y = np.dot(perp, dot_test)
        bounds[i] = [x.min(), y.min(), x.max(), y.max()]

    extents  = np.diff(bounds.reshape((-1,2,2)), axis=1).reshape((-1,2))
    area     = np.product(extents, axis=1)
    area_min = area.argmin()

    offset = -bounds[area_min][0:2]
    theta  = np.arctan2(*edge_vectors[area_min][::-1])

    transform = transformation_2D(offset, theta)
    rectangle = extents[area_min]

    return transform, rectangle

def convex_hull(mesh, clean=True):
    '''
    Get a new Trimesh object representing the convex hull of the 
    current mesh. Requires scipy >.12.

    Argments
    --------
    clean: boolean, if True will fix normals and winding
           to be coherent (as qhull/scipy outputs are not)

    Returns
    --------
    convex: Trimesh object of convex hull of current mesh
    '''

    type_trimesh = type_named(mesh, 'Trimesh')
    c = ConvexHull(mesh.vertices.view(np.ndarray).reshape((-1,3)))

    vid = np.sort(c.vertices)
    mask = np.zeros(len(c.points), dtype=np.int64)
    mask[vid] = np.arange(len(vid))

    faces    = mask[c.simplices]
    vertices = c.points[vid].copy()

    convex = type_trimesh(vertices = vertices, 
                          faces    = faces,
                          process  = True)
    if clean:
        # the normals and triangle winding returned by scipy/qhull's
        # ConvexHull are apparently random, so we need to completely fix them
        convex.fix_normals()
    return convex

def planar_hull(points, normal, origin=None, input_convex=False):
    '''
    Find the convex outline of a set of points projected to a plane.

    Arguments
    -----------
    points: (n,3) float, input points
    normal: (3) float vector, normal vector of plane
    origin: (3) float, location of plane origin
    input_convex: bool, if True we assume the input points are already from
                  a convex hull which provides a speedup. 

    Returns
    -----------
    hull_lines: (n,2,2) set of unordered line segments
    T:          (4,4) float, transformation matrix 
    '''
    if origin is None:
        origin = np.zeros(3)
    if not input_convex:
        pass
    planar, T = project_to_plane(points,
                                 plane_normal     = normal,
                                 plane_origin     = origin,
                                 return_planar    = False,
                                 return_transform = True)
    hull_edges = ConvexHull(planar[:,0:2]).simplices
    hull_lines = planar[hull_edges]
    planar_z = planar[:,2]
    height = np.array([planar_z.min(),
                       planar_z.max()])
    return hull_lines, T, height

def calculate_boundary(vertices):

    # find the points that actually comprise a convex hull
    hull = ConvexHull(list(vertices))

    # keep only vertices that actually comprise a convex hull and arrange in CCW order
    vertices = vertices[hull.vertices]

    # get the real number of vertices
    nVertices = vertices.shape[0]

    # initialize normals array
    unit_normals = np.zeros([nVertices, 2])

    # determine if point is inside or outside of each face, and distance from each face
    for j in range(0, nVertices):

        # calculate the unit normal vector of the current face (taking points CCW)
        if j < nVertices - 1:  # all but the set of point that close the shape
            normal = np.array([vertices[j+1, 1]-vertices[j, 1],
                               -(vertices[j+1, 0]-vertices[j, 0])])
            unit_normals[j] = normal/np.linalg.norm(normal)
        else:   # the set of points that close the shape
            normal = np.array([vertices[0, 1]-vertices[j, 1],
                               -(vertices[0, 0]-vertices[j, 0])])
            unit_normals[j] = normal/np.linalg.norm(normal)

    return vertices, unit_normals

def get_normals(hull, number_fitnodes=12):
    """
    Calculate normals from given hull points using local convex hull fitting.
    Orientation of normals is found using local center of mass.
    :param hull: 3D coordinates of points representing cell hull
    :type hull: np.array
    :return: normals for each hull point
    """
    normals = np.zeros_like(hull, dtype=np.float)
    hull_tree = spatial.cKDTree(hull)
    dists, nearest_nodes_ixs = hull_tree.query(hull, k=number_fitnodes,
                                                distance_upper_bound=1000)
    for ii, nearest_ixs in enumerate(nearest_nodes_ixs):
        nearest_nodes = hull[nearest_ixs[dists[ii] != np.inf]]
        ch = ConvexHull(nearest_nodes, qhull_options='QJ Pp')
        triangles = ch.points[ch.simplices]
        normal = np.zeros((3), dtype=np.float)
        # average normal
        for triangle in triangles:
            cnt = 0
            n_help = unit_normal(triangle[0], triangle[1], triangle[2])
            if not np.any(np.isnan(n_help)):
                normal += np.abs(n_help)
        normal /= np.linalg.norm(normal)
        normal_sign = (hull[ii] - np.mean(nearest_nodes, axis=0))/\
                      np.abs(hull[ii] - np.mean(nearest_nodes, axis=0))
        normals[ii] = normal * normal_sign
    return normals

def convexhull_example(self, length, scales):
        points = np.random.uniform(0, 1, [length, self.input_size])
        target = -1 * np.ones([length])
        ch = ConvexHull(points).vertices
        argmin = np.argsort(ch)[0]
        ch = list(ch[argmin:]) + list(ch[:argmin])
        target[:len(ch)] = np.array(ch)
        target += 1
        return points, target

def _in_volume_convex(points, volume, remote_instance=None, approximate=False, ignore_axis=[]):
    """ Uses scipy to test if points are within a given CATMAID volume.
    The idea is to test if adding the point to the cloud would change the
    convex hull. 
    """

    if remote_instance is None:
        if 'remote_instance' in sys.modules:
            remote_instance = sys.modules['remote_instance']
        elif 'remote_instance' in globals():
            remote_instance = globals()['remote_instance']

    if type(volume) == type(str()):
        volume = fetch.get_volume(volume, remote_instance)

    verts = volume['vertices']

    if not approximate:
        intact_hull = ConvexHull(verts)
        intact_verts = list(intact_hull.vertices)

        if isinstance(points, list):
            points = np.array(points)
        elif isinstance(points, pd.DataFrame):
            points = points.to_matrix()

        return [list(ConvexHull(np.append(verts, list([p]), axis=0)).vertices) == intact_verts for p in points]
    else:
        bbox = [(min([v[0] for v in verts]), max([v[0] for v in verts])),
                (min([v[1] for v in verts]), max([v[1] for v in verts])),
                (min([v[2] for v in verts]), max([v[2] for v in verts]))
                ]

        for a in ignore_axis:
            bbox[a] = (float('-inf'), float('inf'))

        return [False not in [bbox[0][0] < p.x < bbox[0][1], bbox[1][0] < p.y < bbox[1][1], bbox[2][0] < p.z < bbox[2][1], ] for p in points]

def get_hull(self):
        points = self.get_points_hull()
        return ConvexHull(points)

def run_hpr(grid, viewpoint, gamma, boundary):
    """Runs the HPR operator for a single viewpoint and its neighborhood

    Runs the HPR operator for a single viewpoint and its neighborhood.

    Args:
        grid: The support grid containing the data
        viewpoint: The viewpoint to be used by the HPR operator.
        gamma: The exponent used to flip the points.
        boundary: An array of flags used as output.
    """
    candidates = grid.get_candidates(viewpoint)

    if (candidates.shape[0] == 0):
        return

    if (candidates.shape[0] <= grid.dimension):
        boundary[candidates] = True
        return

    flipped = exponential_flip(grid.points[candidates], viewpoint, gamma)

    # add the viewpoint to the end of the list
    flipped = np.vstack([flipped, np.zeros([1, grid.dimension])])

    hull = ConvexHull(flipped)
    visible_idx = hull.vertices

    # remove the index corresponding to the viewpoint
    visible_idx.sort()
    visible_idx = np.delete(visible_idx, -1)

    visible_idx = candidates[visible_idx]

    boundary[visible_idx] = True

def compute_polygon_shadow(n, light_poly, ground_poly):
    """Polygons are described by their convex hulls.

    n -- normal vector used in ZMP computations
    light_poly -- polygon acting as light source
    ground_poly -- polygon whose shadow is projected under the light one

    """
    t = [n[2] - n[1], n[0] - n[2], n[1] - n[0]]
    n = array(n) / pymanoid.toolbox.norm(n)
    t = array(t) / pymanoid.toolbox.norm(t)
    b = cross(n, t)

    light_proj = [array([dot(t, p), dot(b, p)]) for p in light_poly]
    light_hull = ConvexHull(light_proj)
    ground_proj = [array([dot(t, p), dot(b, p)]) for p in ground_poly]
    ground_hull = ConvexHull(ground_proj)
    vertex2poly = {i: j for (i, j) in enumerate(ground_hull.vertices)}
    light_vertices = [light_proj[i] for i in light_hull.vertices]
    ground_vertices = [ground_proj[i] for i in ground_hull.vertices]
    mink_diff = [gv - lv for gv in ground_vertices for lv in light_vertices]

    try:
        u_low, u_high = pymanoid.draw.pick_2d_extreme_rays(mink_diff)
    except pymanoid.exceptions.UnboundedPolyhedron:
        big_dist = 1000  # [m]
        vertices = [
            array([-big_dist, -big_dist, light_poly[0][2]]),
            array([-big_dist, +big_dist, light_poly[0][2]]),
            array([+big_dist, -big_dist, light_poly[0][2]]),
            array([+big_dist, +big_dist, light_poly[0][2]])]
        return vertices, []

    nb_vertices = len(ground_vertices)
    vertex_indices = range(len(ground_vertices))

    def f_low(i):
        return cross(u_low, ground_vertices[i])

    def f_high(i):
        return cross(u_high, ground_vertices[i])

    v_low = min(vertex_indices, key=f_low)
    v_high = max(vertex_indices, key=f_high)
    vertices = [ground_poly[vertex2poly[vertex_index % nb_vertices]]
                for vertex_index in xrange(v_high, (v_low + nb_vertices + 1))]
    rays = [u_low[0] * t + u_low[1] * b,
            u_high[0] * t + u_high[1] * b]
    return vertices, rays

def sortCorners(corners):
    '''
    sort the corners of a given quadrilateral of the type
    corners : [ [xi,yi],... ]

    to an anti-clockwise order starting with the bottom left corner

    or (if plotted as image where y increases to the bottom):
    clockwise, starting top left
    '''
    corners = np.asarray(corners)
    # bring edges in order:
    corners2 = corners[ConvexHull(corners).vertices]

    if len(corners2) == 3:
        # sometimes ConvexHull one point is missing because it is
        # within the hull triangle
        # find the right position of set corner as the minimum perimeter
        # built with that point as different indices
        for c in corners:
            if c not in corners2:
                break
        perimeter = []
        for n in range(0, 4):
            corners3 = np.insert(corners2, n, c, axis=0)
            perimeter.append(
                np.linalg.norm(
                    np.diff(
                        corners3,
                        axis=0),
                    axis=1).sum())
        n = np.argmin(perimeter)
        corners2 = np.insert(corners2, n, c, axis=0)

    # find the edge with the right angle to the quad middle:
    mn = corners2.mean(axis=0)
    d = (corners2 - mn)
    ascent = np.arctan2(d[:, 1], d[:, 0])
    bl = np.abs(BL_ANGLE + ascent).argmin()
    # build a index list starting with bl:
    i = list(range(bl, 4))
    i.extend(list(range(0, bl)))
    return corners2[i]

def lambda_cone(accuracy, runtime, ax, c, conesize, lines=1, aspect_equal=1):
    #ax.scatter(runtime, accuracy, c=tradeoff, lw=0)

    P = zip(runtime, accuracy)
    #P.extend([(0,0), (runtime.max()+0.5e7, accuracy.max())])
    P.extend([(runtime.min()-.1*runtime.ptp(), 0),
              (runtime.max()+.1*runtime.ptp(), accuracy.max())])

#    tradeoff = np.array(list(tradeoff) + [np.inf, 0.0])

    hull = ConvexHull(P)
    v = hull.points[hull.vertices]

    ddd = []

    V = len(hull.vertices)
    for i in range(V):
        x1,y1 = v[(i-1) % V]
        x,y = v[i]
        x3,y3 = v[(i+1) % V]

        # slopes give a range for lambda values to try.
        m12 = (y-y1)/(x-x1)
        m23 = (y3-y)/(x3-x)

        # Note: This filter skips cases where the convex hull contains a
        # segments (edge) that is beneath the fronter (we know this because of
        # the orientation of the edge)..
        if m12 >= m23:
            continue

        ddd.append([m12, m23, x, y])

        ax.add_artist(Polygon([[x,y],
                               point(x, y, m12, -conesize),
                               point(x, y, m23, -conesize),
                               [x,y]],
                              linewidth=0,
                              alpha=0.5,
                              color=c,
                          ))

        if lines:
            ax.plot([x1,x,x3], [y1,y,y3], c=c, alpha=0.5)

    if aspect_equal:
        ax.set_aspect('equal')
    ax.grid(True)
    return ddd