Find if a point is inside a convex hull for a set of points without computing the hull itself

Question

What is the simplest way to test if a point P is inside a convex hull formed by a set of points X?

I'd like an algorithm that works in a high-dimensional space (say, up to 40 dimensions) that doesn't explicitly compute the convex hull itself. Any ideas?

Answer 1

The problem can be solved by finding a feasible point of a Linear Program. If you're interested in the full details, as opposed to just plugging an LP into an existing solver, I'd recommend reading Chapter 11.4 in Boyd and Vandenberghe's excellent book on convex optimization.

Set A = (X[1] X[2] ... X[n]), that is, the first column is v1, the second v2, etc.

Solve the following LP problem,

minimize (over x): 1
s.t.     Ax = P
         x^T * [1] = 1
         x[i] >= 0  \forall i

where

1. x^T is the transpose of x
2. [1] is the all-1 vector.

The problem has a solution iff the point is in the convex hull.

Answer 2

The point lies outside of the convex hull of the other points if and only if the direction of all the vectors from it to those other points are on less than one half of a circle/sphere/hypersphere around it.

Here is a sketch for the situation of two points, a blue one inside the convex hull (green) and the red one outside:

For the red one, there exist bisections of the circle, such that the vectors from the point to the points on the convex hull intersect only one half of the circle. For the blue point, it is not possible to find such a bisection.

Answer 3

You don't have to compute convex hull itself, as it seems quite troublesome in multidimensional spaces. There's a well-known property of convex hulls:

Any vector (point) v inside convex hull of points [v1, v2, .., vn] can be presented as sum(ki*vi), where 0 <= ki <= 1 and sum(ki) = 1. Correspondingly, no point outside of convex hull will have such representation.

In m-dimensional space, this will give us the set of m linear equations with n unknowns.

edit
I'm not sure about complexity of this new problem in general case, but for m = 2 it seems linear. Perhaps, somebody with more experience in this area will correct me.

Answer 4

I had the same problem with 16 dimensions. Since even qhull didn't work properly as too much faces had to be generated, I developed my own approach by testing, whether a separating hyperplane can be found between the new point and the reference data (I call this "HyperHull" ;) ).

The problem of finding a separating hyperplane can be transformed to a convex quadratic programming problem (see: SVM). I did this in python using cvxopt with less then 170 lines of code (including I/O). The algorithm works without modification in any dimension even if there exists the problem, that as higher the dimension as higher the number of points on the hull (see: On the convex hull of random points in a polytope). Since the hull isn't explicitely constructed but only checked, whether a point is inside or not, the algorithm has very big advantages in higher dimensions compared to e.g. quick hull.

This algorithm can 'naturally' be parallelized and speed up should be equal to number of processors.

Answer 5

Though the original post was three years ago, perhaps this answer will still be of help. The Gilbert-Johnson-Keerthi (GJK) algorithm finds the shortest distance between two convex polytopes, each of which is defined as the convex hull of a set of generators---notably, the convex hull itself does not have to be calculated. In a special case, which is the case being asked about, one of the polytopes is just a point. Why not try using the GJK algorithm to calculate the distance between P and the convex hull of the points X? If that distance is 0, then P is inside X (or at least on its boundary). A GJK implementation in Octave/Matlab, called ClosestPointInConvexPolytopeGJK.m, along with supporting code, is available at http://www.99main.com/~centore/MunsellAndKubelkaMunkToolbox/MunsellAndKubelkaMunkToolbox.html . A simple description of the GJK algorithm is available in Sect. 2 of a paper, at http://www.99main.com/~centore/ColourSciencePapers/GJKinConstrainedLeastSquares.pdf . I've used the GJK algorithm for some very small sets X in 31-dimensional space, and had good results. How the performance of GJK compares to the linear programming methods that others are recommending is uncertain (although any comparisons would be interesting). The GJK method does avoid computing the convex hull, or expressing the hull in terms of linear inequalities, both of which might be time-consuming. Hope this answer helps.

Answer 6

Are you willing to accept a heuristic answer that should usually work but is not guaranteed to? If you are then you could try this random idea.

Let f(x) be the cube of the distance to P times the number of things in X, minus the sum of the cubes of the distance to all of the points in X. Start somewhere random, and use a hill climbing algorithm to maximize f(x) for x in a sphere that is very far away from P. Excepting degenerate cases, if P is not in the convex hull this should have a very good probability of finding the normal to a hyperplane which P is on one side of, and everything in X is on the other side of.

Answer 7

A write-up to test if a point is in a hull space, using scipy.optimize.minimize.

Based on user1071136's answer.

It does go a lot faster if you compute the convex hull, so I added a couple of lines for people who want to do that. I switched from graham scan (2D only) to the scipy qhull algorithm.

scipy.optimize.minimize documentation:
https://docs.scipy.org/doc/scipy/reference/optimize.nonlin.html

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


def hull_test(P, X, use_hull=True, verbose=True, hull_tolerance=1e-5, return_hull=True):
    if use_hull:
        hull = ConvexHull(X)
        X = X[hull.vertices]

    n_points = len(X)

    def F(x, X, P):
        return np.linalg.norm( np.dot( x.T, X ) - P )

    bnds = [[0, None]]*n_points # coefficients for each point must be > 0
    cons = ( {'type': 'eq', 'fun': lambda x: np.sum(x)-1} ) # Sum of coefficients must equal 1
    x0 = np.ones((n_points,1))/n_points # starting coefficients

    result = scipy.optimize.minimize(F, x0, args=(X, P), bounds=bnds, constraints=cons)

    if result.fun < hull_tolerance:
        hull_result = True
    else:
        hull_result = False

    if verbose:
        print( '# boundary points:', n_points)
        print( 'x.T * X - P:', F(result.x,X,P) )
        if hull_result: 
            print( 'Point P is in the hull space of X')
        else: 
            print( 'Point P is NOT in the hull space of X')

    if return_hull:
        return hull_result, X
    else:
        return hull_result

Test on some sample data:

n_dim = 3
n_points = 20
np.random.seed(0)

P = np.random.random(size=(1,n_dim))
X = np.random.random(size=(n_points,n_dim))

_, X_hull = hull_test(P, X, use_hull=True, hull_tolerance=1e-5, return_hull=True)

Output:

# boundary points: 14
x.T * X - P: 2.13984259782e-06
Point P is in the hull space of X

Visualize it:

rows = max(1,n_dim-1)
cols = rows
plt.figure(figsize=(rows*3,cols*3))
for row in range(rows):
    for col in range(row, cols):
        col += 1
        plt.subplot(cols,rows,row*rows+col)
        plt.scatter(P[:,row],P[:,col],label='P',s=300)
        plt.scatter(X[:,row],X[:,col],label='X',alpha=0.5)
        plt.scatter(X_hull[:,row],X_hull[:,col],label='X_hull')
        plt.xlabel('x{}'.format(row))
        plt.ylabel('x{}'.format(col))
plt.tight_layout()