Python scipy.spatial.qhull.QhullError() Examples

def test_wrong_feasible_point(self):
        halfspaces = np.array([[-1.0, 0.0, 0.0],
                               [0.0, -1.0, 0.0],
                               [1.0, 0.0, -1.0],
                               [0.0, 1.0, -1.0]])
        feasible_point = np.array([0.5, 0.5, 0.5])
        #Feasible point is (ndim,) instead of (ndim-1,)
        assert_raises(ValueError, qhull.HalfspaceIntersection, halfspaces, feasible_point)
        feasible_point = np.array([[0.5], [0.5]])
        #Feasible point is (ndim-1, 1) instead of (ndim-1,)
        assert_raises(ValueError, qhull.HalfspaceIntersection, halfspaces, feasible_point)
        feasible_point = np.array([[0.5, 0.5]])
        #Feasible point is (1, ndim-1) instead of (ndim-1,)
        assert_raises(ValueError, qhull.HalfspaceIntersection, halfspaces, feasible_point)

        feasible_point = np.array([-0.5, -0.5])
        #Feasible point is outside feasible region
        assert_raises(qhull.QhullError, qhull.HalfspaceIntersection, halfspaces, feasible_point) 

def cvHull(pL):
    if len(pL) < 3:
        return pL
    from scipy.spatial import ConvexHull
    from scipy.spatial.qhull import QhullError
    S = [xx.std() for xx in pL]
    try:
        hull = ConvexHull(S)
    except QhullError:
        xL = [xx.x for xx in pL]
        xx, pL = order(xL, pL)
        return [pL[0], pL[-1]]
    P = [pL[i] for i in hull.vertices]
    P.append(pL[hull.vertices[0]])
    return P 

def patch_up_roi(roi):
    """
    After being non-linearly transformed, ROIs tend to have holes in them.
    We perform a couple of computational geometry operations on the ROI to
    fix that up.

    Parameters
    ----------
    roi : 3D binary array
        The ROI after it has been transformed.

    sigma : float
        The sigma for initial Gaussian smoothing.

    truncate : float
        The truncation for the Gaussian

    Returns
    -------
    ROI after dilation and hole-filling
    """

    hole_filled = ndim.binary_fill_holes(roi > 0)
    try:
        return convex_hull_image(hole_filled)
    except QhullError:
        return hole_filled 

def setRegion(self, safe_region):
        debug = DebugData()
        pos = safe_region.point
        try:
            xy_verts = safe_region.xy_polytope()
            if xy_verts.shape[1] == 0:
                raise QhullError("No points returned")
            xyz_verts = np.vstack((xy_verts, pos[2] + 0.02 + np.zeros((1, xy_verts.shape[1]))))
            xyz_verts = np.hstack((xyz_verts, np.vstack((xy_verts, pos[2] + 0.015 + np.zeros((1, xy_verts.shape[1]))))))
            # print xyz_verts.shape
            polyData = vnp.getVtkPolyDataFromNumpyPoints(xyz_verts.T.copy())
            vol_mesh = filterUtils.computeDelaunay3D(polyData)
            for j in range(xy_verts.shape[1]):
                z = pos[2] + 0.005
                p1 = np.hstack((xy_verts[:,j], z))
                if j < xy_verts.shape[1] - 1:
                    p2 = np.hstack((xy_verts[:,j+1], z))
                else:
                    p2 = np.hstack((xy_verts[:,0], z))
                debug.addLine(p1, p2, color=[.7,.7,.7], radius=0.003)
            debug.addPolyData(vol_mesh)
            # self.setPolyData(vol_mesh)
            self.setPolyData(debug.getPolyData())
            self.safe_region = safe_region
        except QhullError:
            print("Could not generate convex hull (polytope is likely unbounded).") 

def convexhull(self, x, y, fill=False, smooth=False):
        try:
            from scipy.spatial import ConvexHull
            from scipy.spatial.qhull import QhullError
        except:
            raise Exception('ConvexHull requires scipy')

        if len(x) < 3:
            raise Exception('convexhull requires at least 3 points')


        points = np.vstack((x,y)).T
        try:
            hull = ConvexHull(points)
            xhull = points[hull.vertices,0]
            yhull = points[hull.vertices,1]

            if smooth:
                xhull, yhull = self.__generate_spline(xhull, yhull, closed=True)

            if fill:
                self.poly(xhull,yhull)
            else:
                self.linestrip(xhull, yhull, 3, closed=True)

        except QhullError as qerr:
            self.linestrip(x, y, 3, closed=False) 

def compute_pendular_accel_cone(self, com_vertices=None, zdd_max=None,
                                    reduced=False):
        """
        Compute the pendular COM acceleration cone of the stance.

        The pendular cone is the reduction of the Contact Wrench Cone when the
        angular momentum at the COM is zero.

        Parameters
        ----------
        com_vertices : list of (3,) arrays, optional
            Vertices of a COM bounding polytope.
        zdd_max : scalar, optional
            Maximum vertical acceleration in the output cone.
        reduced : bool, optional
            If ``True``, returns the reduced 2D form rather than a 3D cone.

        Returns
        -------
        vertices : list of (3,) arrays
            List of 3D vertices of the (truncated) COM acceleration cone, or of
            the 2D vertices of the reduced form if ``reduced`` is ``True``.

        Notes
        -----
        The method is based on a rewriting of the CWC formula, followed by a 2D
        convex hull on dual vertices. The algorithm is described in [Caron16]_.

        When ``com`` is a list of vertices, the returned cone corresponds to
        COM accelerations that are feasible from *all* COM located inside the
        polytope. See [Caron16]_ for details on this conservative criterion.
        """
        def expand_reduced_pendular_cone(reduced_hull, zdd_max=None):
            g = -gravity[2]  # gravity constant (positive)
            zdd = +g if zdd_max is None else zdd_max
            vertices_at_zdd = [
                array([a * (g + zdd), b * (g + zdd), zdd])
                for (a, b) in reduced_hull]
            return [gravity] + vertices_at_zdd

        if com_vertices is None:
            com_vertices = [self.com.p]
        CWC_O = self.compute_wrench_inequalities([0., 0., 0.])
        B_list, c_list = [], []
        for (i, v) in enumerate(com_vertices):
            B = CWC_O[:, :3] + cross(CWC_O[:, 3:], v)
            c = dot(B, gravity)
            B_list.append(B)
            c_list.append(c)
        B = vstack(B_list)
        c = hstack(c_list)
        try:
            g = -gravity[2]  # gravity constant (positive)
            B_2d = hstack([B[:, j].reshape((B.shape[0], 1)) for j in [0, 1]])
            sigma = c / g  # see Equation (30) in [CK16]
            reduced_hull = compute_polygon_hull(B_2d, sigma)
            if reduced:
                return reduced_hull
            return expand_reduced_pendular_cone(reduced_hull, zdd_max)
        except QhullError:
            raise Exception("Cannot compute 2D polar for acceleration cone") 

def draw_polygon(points, normal, combined='g-#', color=None, faces=None,
                 linewidth=1., pointsize=0.01):
    """
    Draw a polygon defined as the convex hull of a set of points. The normal
    vector n of the plane containing the polygon must also be supplied.

    Parameters
    ----------
    points : list of arrays
        List of coplanar 3D points.
    normal : array, shape=(3,)
        Unit vector normal to the drawing plane.
    combined : string, optional
        Drawing spec in matplotlib fashion. Default: 'g-#'.
    color : char or triplet, optional
        Color letter or RGB values, default is 'g' for green.
    faces : string
        Faces of the polyhedron to draw. Use '.' for vertices, '-' for edges
        and '#' for facets.
    linewidth : scalar
        Thickness of drawn line.
    pointsize : scalar
        Vertex size.

    Returns
    -------
    handles : list of openravepy.GraphHandle
        OpenRAVE graphical handles. Must be stored in some variable, otherwise
        the drawn object will vanish instantly.
    """
    assert abs(1. - norm(normal)) < 1e-10
    n = normal
    t = array([n[2] - n[1], n[0] - n[2], n[1] - n[0]], dtype=float)
    t /= norm(t)
    b = cross(n, t)
    points2d = [[dot(t, x), dot(b, x)] for x in points]
    try:
        hull = ConvexHull(points2d)
    except QhullError:
        warn("QhullError: maybe polygon is empty?")
        return []
    except IndexError:
        warn("Qhull raised an IndexError for points2d=%s" % repr(points2d))
        return []
    return draw_polytope(
        points, combined, color, faces, linewidth, pointsize, hull=hull) 

def get_onion_layers(structure):
    """
    Returns the different layers of a finite structure

    :param structure:
    :return:
    """
    assert (not structure.is_periodic)
    layers = []
    cur_st = structure.copy()

    morelayers = True
    while morelayers:
        pos = cur_st.positions
        if len(pos) <= 4:
            core = range(cur_st.natom)
            layers.append(core)
            break

        st = pychemia.Structure(positions=pos, symbols=len(pos)*['H'], periodicity=False)
        st.canonical_form()
        # print('The current volume is %7.3f' % st.volume)
        if st.volume < 0.1:
            core = range(cur_st.natom)
            layers.append(core)
            break

        try:
            voro = scipy.spatial.Voronoi(pos)
            surface = [i for i in range(cur_st.natom) if -1 in voro.regions[voro.point_region[i]]]
        except qhull.QhullError:
            surface = range(cur_st.natom)
            morelayers = False
        layers.append(surface)
        if not morelayers:
            break
        core = [i for i in range(cur_st.natom) if i not in surface]
        if len(core) == 0:
            break
        symbols = list(np.array(cur_st.symbols)[np.array(core)])
        positions = cur_st.positions[np.array(core)]
        cur_st = pychemia.Structure(symbols=symbols, positions=positions, periodicity=False)
    new_layers = [layers[0]]
    included = list(layers[0])
    acumulator = 0
    for i in range(1, len(layers)):
        noincluded = [j for j in range(structure.natom) if j not in included]
        # print 'Layer: %3d   Atoms on Surface: %3d     Internal Atoms: %3d' % (i,
        #                                                                      len(included) - acumulator,
        #                                                                      len(noincluded))
        acumulator = len(included)
        relabel = [noincluded[j] for j in layers[i]]
        included += relabel
        new_layers.append(relabel)

    return new_layers