What's an efficient way to find if a point lies in the convex hull of a point cloud?

Question

I have a point cloud of coordinates in numpy. For a high number of points, I want to find out if the points lie in the convex hull of the point cloud.

I tried pyhull but I cant figure out how to check if a point is in the ConvexHull:

hull = ConvexHull(np.array([(1, 2), (3, 4), (3, 6)]))
for s in hull.simplices:
    s.in_simplex(np.array([2, 3]))

raises LinAlgError: Array must be square.

Answer 1

Here is an easy solution that requires only scipy:

def in_hull(p, hull):
    """
    Test if points in `p` are in `hull`

    `p` should be a `NxK` coordinates of `N` points in `K` dimensions
    `hull` is either a scipy.spatial.Delaunay object or the `MxK` array of the 
    coordinates of `M` points in `K`dimensions for which Delaunay triangulation
    will be computed
    """
    from scipy.spatial import Delaunay
    if not isinstance(hull,Delaunay):
        hull = Delaunay(hull)

    return hull.find_simplex(p)>=0

It returns a boolean array where True values indicate points that lie in the given convex hull. It can be used like this:

tested = np.random.rand(20,3)
cloud  = np.random.rand(50,3)

print in_hull(tested,cloud)

If you have matplotlib installed, you can also use the following function that calls the first one and plots the results. For 2D data only, given by Nx2 arrays:

def plot_in_hull(p, hull):
    """
    plot relative to `in_hull` for 2d data
    """
    import matplotlib.pyplot as plt
    from matplotlib.collections import PolyCollection, LineCollection

    from scipy.spatial import Delaunay
    if not isinstance(hull,Delaunay):
        hull = Delaunay(hull)

    # plot triangulation
    poly = PolyCollection(hull.points[hull.vertices], facecolors='w', edgecolors='b')
    plt.clf()
    plt.title('in hull')
    plt.gca().add_collection(poly)
    plt.plot(hull.points[:,0], hull.points[:,1], 'o', hold=1)


    # plot the convex hull
    edges = set()
    edge_points = []

    def add_edge(i, j):
        """Add a line between the i-th and j-th points, if not in the list already"""
        if (i, j) in edges or (j, i) in edges:
            # already added
            return
        edges.add( (i, j) )
        edge_points.append(hull.points[ [i, j] ])

    for ia, ib in hull.convex_hull:
        add_edge(ia, ib)

    lines = LineCollection(edge_points, color='g')
    plt.gca().add_collection(lines)
    plt.show()    

    # plot tested points `p` - black are inside hull, red outside
    inside = in_hull(p,hull)
    plt.plot(p[ inside,0],p[ inside,1],'.k')
    plt.plot(p[-inside,0],p[-inside,1],'.r')

Answer 2

I would not use a convex hull algorithm, because you do not need to compute the convex hull, you just want to check whether your point can be expressed as a convex combination of the set of points of whom a subset defines a convex hull. Moreover, finding the convex hull is computationally expensive, especially in higher dimensions.

In fact, the mere problem of finding out whether a point can be expressed as a convex combination of another set of points can be formulated as a linear programming problem.

import numpy as np
from scipy.optimize import linprog

def in_hull(points, x):
    n_points = len(points)
    n_dim = len(x)
    c = np.zeros(n_points)
    A = np.r_[points.T,np.ones((1,n_points))]
    b = np.r_[x, np.ones(1)]
    lp = linprog(c, A_eq=A, b_eq=b)
    return lp.success

n_points = 10000
n_dim = 10
Z = np.random.rand(n_points,n_dim)
x = np.random.rand(n_dim)
print(in_hull(Z, x))

For the example, I solved the problem for 10000 points in 10 dimensions. The executions time is in the ms range. Would not want to know how long this would take with QHull.

Answer 3

Hi I am not sure about how to use your program library to achieve this. But there is a simple algorithm to achieve this described in words:

1. create a point that is definitely outside your hull. Call it Y
2. produce a line segment connecting your point in question (X) to the new point Y.
3. loop around all edge segments of your convex hull. check for each of them if the segment intersects with XY.
4. If the number of intersection you counted is even (including 0), X is outside the hull. Otherwise X is inside the hull.
5. if so occurs XY pass through one of your vertexes on the hull, or directly overlap with one of your hull's edge, move Y a little bit.
6. the above works for concave hull as well. You can see in below illustration (Green dot is the X point you are trying to determine. Yellow marks the intersection points.