Point in convex polygon

Question

I'm trying to make a function that will return True if the given (x,y) point is inside a convex polygon. I'm trying to make it without numpy or any similar imports, just pure python code.

I've already found a sample solution, which seems OK at first sight, but it's not working correctly, and I can't figure out why. The code is as follows:

def point_in_poly(x,y,poly):

  n = len(poly)
  inside = False

  p1x,p1y = poly[0]
  for i in range(n+1):
      p2x,p2y = poly[i % n]
      if y > min(p1y,p2y):
          if y <= max(p1y,p2y):
              if x <= max(p1x,p2x):
                  if p1y != p2y:
                      xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
                  if p1x == p2x or x <= xints:
                      inside = not inside
      p1x,p1y = p2x,p2y

  return inside

If I test it for (9,9), for the following polygon, it gives me True:

polygon = [(0,10),(10,10),(10,0),(0,0)]
point_x = 9
point_y = 9
print point_in_poly(point_x,point_y,polygon)

But when I change the order of the points of the polygon, for the same point, it gives me False:

polygon = [(0,0), (0,10), (10,0), (10,10)]
point_x = 9
point_y = 9
print point_in_poly(point_x,point_y,polygon)

Anybody knows the reason? Thanks!

Answer 1

In the particular case you are having problems with is special: polygon = [(0,0), (0,10), (10,0), (10,10)]

Changing the order of points in a polygon can have significant impact on algorithms.

If you draw your polygon on a graph you'll see you have a horizontal hourglass shape. The polygon border overlaps itself. In geospatial analysis this overlap is not allowed because visually and logically you now have two closed polygons with a common intersection point. By the way most geospatial software doesn't deal well with triangles either.

In this case the point at 9,9 will trick the ray casting algorithm used in your method above because it can easily cross the doubled-over polygon boundary twice.

Please run the following code to see what is going on. (9,9) is on the line and this algorithm doesn't account for it. (5,8) is way outside:

import turtle as t
polygon = [(0,0), (0,100), (100,0), (100,100)]
t.goto(0,0)
fp = None
for p in polygon:
  t.goto(p)
  if not fp: fp=p
t.goto(fp)
t.up()
t.goto(90,90)
t.write("90,90")
t.dot(10)
t.goto(50,80)
t.write("50,80")
t.dot(10)
t.done()

This code handles the (9,9) edge case: http://geospatialpython.com/2011/08/point-in-polygon-2-on-line.html

Answer 2

the point 9,0 is not inside the polygon [(0,10),(10,10),(10,0),(0,0)] its on the edge. Points exactly on the edge can be considered in or out depending on the specifics of your algorithm.

Answer 3

For polygons with lots of points it often pays to check if the point falls outside of the bounding box first:

def point_in_poly(x, y, poly):
    # Check bounding box first
    xl = [p[0] for p in poly]
    yl = [p[1] for p in poly]
    if x < min(xl) or x > max(xl) or y < min(yl) or y > max(yl):
       return False
   # Now check all points.

You can find a lot of algoriths dealing with polygons and meshes at one of Paul Bourke's webpages.

For a lot of these algorithms it pays to use numpy, though, since a lot of the steps have to be done for each point in the poly array.