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

For those that may be interested, I wrote my own implementation that does not rely on scipy.spatial at all (although scipy.optimize is used to compute the Chebyshev center if requested), which is both faster and provides extra functionalities. However, it only supports 2D space. The module is publicly available on Github.

Answer 2

Vectorized, with plotting (based on other great answers here: 1, 2 ):

import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial import ConvexHull


def points_in_hull(points: np.ndarray, hull: ConvexHull, tolerance: float = 1e-12):
    return np.all(np.add(points @ hull.equations[:, :-1].T, hull.equations[:, -1]) <= tolerance, axis=1)


points = np.array([(1, 2), (3, 4), (3, 6), (2, 4.5), (2.5, 5)])
hull = ConvexHull(points)

np.random.seed(1)
random_points = np.random.uniform(0, 6, (1000, 2))

in_hull = points_in_hull(random_points, hull)
random_points_in_hull = random_points[in_hull]

# plot
fig, ax = plt.subplots(figsize=(20, 10))
plt.scatter(random_points[:, 0], random_points[:, 1], color='blue')
plt.scatter(points[:, 0], points[:, 1], color='green')
plt.scatter(random_points_in_hull[:, 0], random_points_in_hull[:, 1], color='red')
for simplex in hull.simplices:
    plt.plot(points[simplex, 0], points[simplex, 1], 'k-')
plt.show()

Answer 3

For those interested, I made a vectorized version of @charlie-brummit answer:

def points_in_hull(p, hull, tol=1e-12):
    return np.all(hull.equations[:,:-1] @ p.T + np.repeat(hull.equations[:,-1][None,:], len(p), axis=0).T <= tol, 0)

where p is now a [N,2] array. It is ~6 times faster than the recommended solution (@Sildoreth answer), and ~10 times faster than the original.

Adapted demonstration without the for loop: (pasted from below to avoid searching in the thread)

import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial import ConvexHull

points = np.array([(1, 2), (3, 4), (3, 6), (2, 4.5), (2.5, 5)])
hull = ConvexHull(points)

np.random.seed(1)
random_points = np.random.uniform(0, 6, (100, 2))

for simplex in hull.simplices:
    plt.plot(points[simplex, 0], points[simplex, 1])

plt.scatter(*points.T, alpha=.5, color='k', s=200, marker='v')

in_hull = points_in_hull(random_points, hull)
plt.scatter(random_points[in_hull,0], random_points[in_hull,1], marker='x', color='g')
plt.scatter(random_points[~in_hull,0], random_points[~in_hull,1], marker='d', color='m')

Answer 4

Building on the work of @Charlie Brummitt, I implemented a more efficient version enabling to check if multiple points are in the convex hull at the same time and replacing any loop by faster linear algebra.

import numpy as np
from scipy.spatial.qhull import _Qhull

def in_hull(points, queries):
    hull = _Qhull(b"i", points,
                  options=b"",
                  furthest_site=False,
                  incremental=False, 
                  interior_point=None)
    equations = hull.get_simplex_facet_array()[2].T
    return np.all(queries @ equations[:-1] < - equations[-1], axis=1)

# ============== Demonstration ================

points = np.random.rand(8, 2)
queries = np.random.rand(3, 2)
print(in_hull(points, queries))

Note that I'm using the lower-level _Qhull class for efficiency.

Answer 5

To piggy-back off of this answer, to check all points in a numpy array at once, this worked for me:

import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial import ConvexHull

points = np.array([(1, 2), (3, 4), (3, 6), (2, 4.5), (2.5, 5)])
hull = ConvexHull(points)

np.random.seed(1)
random_points = np.random.uniform(0, 6, (100, 2))

# get array of boolean values indicating in hull if True
in_hull = np.all(np.add(np.dot(random_points, hull.equations[:,:-1].T),
                        hull.equations[:,-1]) <= tolerance, axis=1)

random_points_in_hull = random_points[in_hull]

Answer 6

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 7

Based on this post, here is my quick-and-dirty solution for a convex regions with 4 sides (you can easily extend it to more)

def same_sign(arr): return np.all(arr > 0) if arr[0] > 0 else np.all(arr < 0)

def inside_quad(pts, pt):
    a =  pts - pt
    d = np.zeros((4,2))
    d[0,:] = pts[1,:]-pts[0,:]
    d[1,:] = pts[2,:]-pts[1,:]
    d[2,:] = pts[3,:]-pts[2,:]
    d[3,:] = pts[0,:]-pts[3,:]
    res = np.cross(a,d)
    return same_sign(res), res

points = np.array([(1, 2), (3, 4), (3, 6), (2.5, 5)])

np.random.seed(1)
random_points = np.random.uniform(0, 6, (1000, 2))

print wlk1.inside_quad(points, random_points[0])
res = np.array([inside_quad(points, p)[0] for p in random_points])
print res[:4]
plt.plot(random_points[:,0], random_points[:,1], 'b.')
plt.plot(random_points[res][:,0], random_points[res][:,1], 'r.')

Answer 8

Use the equations attribute of ConvexHull:

def point_in_hull(point, hull, tolerance=1e-12):
    return all(
        (np.dot(eq[:-1], point) + eq[-1] <= tolerance)
        for eq in hull.equations)

In words, a point is in the hull if and only if for every equation (describing the facets) the dot product between the point and the normal vector (eq[:-1]) plus the offset (eq[-1]) is less than or equal to zero. You may want to compare to a small, positive constant tolerance = 1e-12 rather than to zero because of issues of numerical precision (otherwise, you may find that a vertex of the convex hull is not in the convex hull).

Demonstration:

import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial import ConvexHull

points = np.array([(1, 2), (3, 4), (3, 6), (2, 4.5), (2.5, 5)])
hull = ConvexHull(points)

np.random.seed(1)
random_points = np.random.uniform(0, 6, (100, 2))

for simplex in hull.simplices:
    plt.plot(points[simplex, 0], points[simplex, 1])

plt.scatter(*points.T, alpha=.5, color='k', s=200, marker='v')

for p in random_points:
    point_is_in_hull = point_in_hull(p, hull)
    marker = 'x' if point_is_in_hull else 'd'
    color = 'g' if point_is_in_hull else 'm'
    plt.scatter(p[0], p[1], marker=marker, color=color)

Answer 9

First, obtain the convex hull for your point cloud.

Then loop over all of the edges of the convex hull in counter-clockwise order. For each of the edges, check whether your target point lies to the "left" of that edge. When doing this, treat the edges as vectors pointing counter-clockwise around the convex hull. If the target point is to the "left" of all of the vectors, then it is contained by the polygon; otherwise, it lies outside the polygon.

This other Stack Overflow topic includes a solution to finding which "side" of a line a point is on: Determine Which Side of a Line a Point Lies

---

The runtime complexity of this approach (once you already have the convex hull) is O(n) where n is the number of edges that the convex hull has.

Note that this will work only for convex polygons. But you're dealing with a convex hull, so it should suit your needs.

It looks like you already have a way to get the convex hull for your point cloud. But if you find that you have to implement your own, Wikipedia has a nice list of convex hull algorithms here: Convex Hull Algorithms

Answer 10

Just for completness, here is a poor's man solution:

import pylab
import numpy
from scipy.spatial import ConvexHull

def is_p_inside_points_hull(points, p):
    global hull, new_points # Remove this line! Just for plotting!
    hull = ConvexHull(points)
    new_points = numpy.append(points, p, axis=0)
    new_hull = ConvexHull(new_points)
    if list(hull.vertices) == list(new_hull.vertices):
        return True
    else:
        return False

# Test:
points = numpy.random.rand(10, 2)   # 30 random points in 2-D
# Note: the number of points must be greater than the dimention.
p = numpy.random.rand(1, 2) # 1 random point in 2-D
print is_p_inside_points_hull(points, p)

# Plot:
pylab.plot(points[:,0], points[:,1], 'o')
for simplex in hull.simplices:
    pylab.plot(points[simplex,0], points[simplex,1], 'k-')
pylab.plot(p[:,0], p[:,1], '^r')
pylab.show()

The idea is simple: the vertices of the convex hull of a set of points P won't change if you add a point p that falls "inside" the hull; the vertices of the convex hull for [P1, P2, ..., Pn] and [P1, P2, ..., Pn, p] are the same. But if p falls "outside", then the vertices must change. This works for n-dimensions, but you have to compute the ConvexHull twice.

Two example plots in 2-D:

False:

True:

Answer 11

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 12

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.

Answer 13

It looks like you are using a 2D point cloud, so I'd like to direct you to the inclusion test for point-in-polygon testing of convex polygons.

Scipy's convex hull algorithm allows for finding convex hulls in 2 or more dimensions which is more complicated than it needs to be for a 2D point cloud. Therefore, I recommend using a different algorithm, such as this one. This is because as you really need for point-in-polygon testing of a convex hull is the list of convex hull points in clockwise order, and a point that is inside of the polygon.

The time performance of this approach is as followed:

• O(N log N) to construct the convex hull
• O(h) in preprocessing to calculate (and store) the wedge angles from the interior point
• O(log h) per point-in-polygon query.

Where N is the number of points in the point cloud and h is the number of points in the point clouds convex hull.

Answer 14

If you want to keep with scipy, you have to convex hull (you did so)

>>> from scipy.spatial import ConvexHull
>>> points = np.random.rand(30, 2)   # 30 random points in 2-D
>>> hull = ConvexHull(points)

then build the list of points on the hull. Here is the code from doc to plot the hull

>>> import matplotlib.pyplot as plt
>>> plt.plot(points[:,0], points[:,1], 'o')
>>> for simplex in hull.simplices:
>>>     plt.plot(points[simplex,0], points[simplex,1], 'k-')

So starting from that, I would propose for computing list of points on the hull

pts_hull = [(points[simplex,0], points[simplex,1]) 
                            for simplex in hull.simplices] 

(although i did not try)

And you can also come with your own code for computing the hull, returning the x,y points.

If you want to know if a point from your original dataset is on the hull, then you are done.

I what you want is to know if a any point is inside the hull or outside, you must do a bit of work more. What you will have to do could be

• for all edges joining two simplices of your hull: decide whether your point is above or under

• if point is below all lines, or above all lines, it is outside the hull

As a speed up, as soon as a point has been above one line and below one another, it is inside the hull.