Finding a partial or complete inclusion of a set of points in convex hull of other set of points in 2D in python

Question

I have a set of two sets of 2D points. I want to see if set B is included completely or partially in the convex hull of set A according to euclidean coordinates.

To explain inclusion the following example might help

Lets consider the following sets

   A={(5,5),(10,10),(5,10),(0,5)}

   B={(3,3),(5,8)} partially included in convex hull of A

   C={(1,5),(5,8)} fully included in convex hull of A

   D={(1,1),(3,3)} is not included in convex hull of A

Thanks a lot

A subproblem of this problem is the "point location" problem. You would first compute a polygon from the base set, compute a point location data structure and then locate individual points. This is usually done with a sweep algorithm. There are solutions in O(n log n), but none of them is trivial to implement. — Helmut Grohne, Jul 24 '12 at 10:44
Maybe see http://stackoverflow.com/questions/1076778/good-geometry-library-in-python. From a first look I'd expect CGAL to be able to do this. — Helmut Grohne, Jul 24 '12 at 11:30

score 2 · Answer 1 · answered Jul 24 '12 at 13:46

Matplotlib has a point_in_poly function that is pretty fast. This is taken straight from the matplotlib documentation: nxutils

In [25]: import numpy as np
In [26]: import matplotlib.nxutils as nx
In [27]: verts = np.array([ [0,0], [0, 1], [1, 1], [1,0]], float)
In [28]: nx.pnpoly( 0.5, 0.5, verts)
Out[28]: 1
In [29]: nx.pnpoly( 0.5, 1.5, verts)
Out[29]: 0
In [30]: points = np.random.rand(10,2)*2
In [31]: points
Out[31]:
array([[ 1.03597426,  0.61029911],
       [ 1.94061056,  0.65233947],
       [ 1.08593748,  1.16010789],
       [ 0.9255139 ,  1.79098751],
       [ 1.54564936,  1.15604046],
       [ 1.71514397,  1.26147554],
       [ 1.19133536,  0.56787764],
       [ 0.40939549,  0.35190339],
       [ 1.8944715 ,  0.61785408],
       [ 0.03128518,  0.48144145]])
In [32]: nx.points_inside_poly(points, verts)
Out[32]: array([False, False, False, False, False, False, False,  True, False, True], dtype=bool)

After that, its just a matter of testing each point in the set and figuring out if both, one, or neither are inside the vertices.

thanks a lot for your answer, the other answer contains complete steps. Thanks again — Shan, Aug 01 '12 at 13:40

score 2 · Accepted Answer · answered Jul 25 '12 at 08:44

One way to find the convex hull of a set of points in Python is to use the Delaunay triangulation function in scipy.spatial. Given a set of points, it returns an object which has a convex_hull attribute - this is an array consisting of pairs of indices into the original set of points, which correspond to edges on the polygon. Annoyingly these are not ordered, so the polygon containing these points needs to be reconstructed (eg. as follows):

import numpy as np
import matplotlib.nxutils
import scipy.spatial

def find_convex_hull(points):
  triangulation = scipy.spatial.Delaunay(points)
  unordered = list(triangulation.convex_hull)
  ordered = list(unordered.pop(0))
  while len(unordered) > 0:
    next = (i for i, seg in enumerate(unordered) if ordered[-1] in seg).next()
    ordered += [point for point in unordered.pop(next) if point != ordered[-1]]
  return points[ordered]

As suggested by @user1443118, the points_inside_poly function in matplotlib.nxutils can then be used to test if points lie in the resulting polygon, which corresponds to the convex hull. This leads to the following function for calculating the degree of intersection.

def inclusion(points_a, points_b):
  ch_a = find_convex_hull(points_a)
  return (1.0 * matplotlib.nxutils.points_inside_poly(points_b, ch_a)).mean()

So given some sets of points (with properties as in your original example), illustrated below:

A = np.random.randn(100, 2)
B = np.array([2,0]) + 0.5 * np.random.randn(100, 2)
C = 0.5 * np.random.randn(100, 2)
D = np.array([5,0]) + 0.5 * np.random.randn(100, 2)

enter image description here

The degree of inclusion can be calculated as follows:

>>> inclusion(A, B)
0.44

>>> inclusion(A, C)
1.0

>>> inclusion(A, D)
0.0

Finally, however, it's worth noting that the the points_in_poly function does not always signify points on the polygon boundary as being inside (see here for an explanation why the underlying function behaves this way). For this reason set C in your original example would only be partially included as point (1,5) lies on the convex hull of A and is not counted.

Finding a partial or complete inclusion of a set of points in convex hull of other set of points in 2D in python

2 Answers2