2

I am having trouble understanding how to compute whether an n-dimensional point is within an n-dimensional ConvexHull.

A very similar question (the same) was asked here: What's an efficient way to find if a point lies in the convex hull of a point cloud?

However, the answers are confusing me or don't seem to work for me, and I have no idea why.

def in_hull(p, hull):
    """ Copied and from the Top Original answer """
    from scipy.spatial import Delaunay
    if not isinstance(hull,Delaunay):
        hull = Delaunay(hull)

    return hull.find_simplex(p)>=0

This function was giving me a lot of false or unwanted results with real data, that I am using it on. However, when debugging I wrote a simple script to test some obvious expectation I would have:

If I construct a ConvexHull out of a group of points, when I check that group of points for "membership", they should all be "members".

results_all = []
for _ in range(5000):
    cloud = np.random.rand(5000, 2)
    result = in_hull(cloud, cloud)
    results_all.append(np.all(result))

arr = np.array(results_all)
print(np.sum(np.logical_not(arr)))

While rare, this seems to fail on randomly generated data (3 out of 5000), the problem is bigger on real data. What I mean by fail, is that I actually get some situations, where not all points are considered members.

Is there something I am doing wrong? Or perhaps completely misunderstanding? I am fairly confused at this point, so would love an explanation of what is going on.

In the end, I want; given a ConvexHull, computed at some previous stage; be able to determine if points lie within the hull.

SirSteel
  • 123
  • 1
  • 10

1 Answers1

5

It seems to be an edge case problem with the find_simplex method of the Delaunay object for almost flat simplex (triangle).

Here is a code to find and plot a faulty case with only 3 points:

import matplotlib.pylab as plt
from scipy.spatial import Delaunay
from scipy.spatial import delaunay_plot_2d

for _ in range(5000):
    cloud = np.random.rand(3, 2)

    tri = Delaunay(cloud)

    if np.any( tri.find_simplex(cloud)<0 ):
        print('break at', _)

        delaunay_plot_2d(tri);
        id_break = np.where(tri.find_simplex(cloud)<0)
        plt.plot( *cloud[id_break].ravel(), 'or' );
        break

faulty example

The other method proposed here seems to work well:

hull = ConvexHull(cloud)

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

[ point_in_hull(point, hull) for point in cloud ]
# [True, True, True]
xdze2
  • 3,986
  • 2
  • 12
  • 29