11

I am making a convex hull using the scipy.spatial package http://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.ConvexHull.html#scipy.spatial.ConvexHull

I've tried searching, but have yet to know if there is an easy way to find if any point (latitude/longitude) lies inside the convex hull. Any suggestions?

bjornasm
  • 2,211
  • 7
  • 37
  • 62

2 Answers2

18

The method that I've used before is using the Path class matplotlib. This has a contains_point method which does exactly that. As well as a contains_points which allows you to query an array of points.

To use this you would do

from scipy.spatial import ConvexHull
from matplotlib.path import Path

hull = ConvexHull( points )

hull_path = Path( points[hull.vertices] )

print hull_path.contains_point((1,2)) # Is (1,2) in the convex hull?
Simon Gibbons
  • 6,969
  • 1
  • 21
  • 34
6

The matplotlib Path method works, but only for 2-dimensional data. A general method that works in arbitrary dimensions is to use the halfspace equations for the facets of the hull.

import numpy as np
from scipy.spatial import ConvexHull

# Assume points are shape (n, d), and that hull has f facets.
hull = ConvexHull(points)

# A is shape (f, d) and b is shape (f, 1).
A, b = hull.equations[:, :-1], hull.equations[:, -1:]

eps = np.finfo(np.float32).eps

def contained(x):
  # The hull is defined as all points x for which Ax + b <= 0.
  # We compare to a small positive value to account for floating
  # point issues.
  #
  # Assuming x is shape (m, d), output is boolean shape (m,).
  return np.all(np.asarray(x) @ A.T + b.T < eps, axis=1)

# To test one point:
print('Point (2,1) is in the hull?', contained([[2, 1]]))
lmjohns3
  • 7,422
  • 5
  • 36
  • 56
  • I'm not sure this works in general. I fed the same points I used to construct the hull into `contained` and found that some of them were not classified as contained. I thought that was because I was using `<` rather than `<=`, but upon changing to the non-strict inequality I still found that points were not contained. – Galen Nov 15 '21 at 19:02
  • 2
    Maybe it is just a floating point precision issue. Subtracting `np.finfo(np.float32).eps` before evaluating the inequality gave the desired result for the `<= 0` comparison. – Galen Nov 15 '21 at 19:32
  • 2
    Yes, thank you! I'd found the same thing in code I have been using but forgot to update this answer. Added the eps to the comparison now. – lmjohns3 Nov 16 '21 at 04:22