42

I have a class describing a Point (has 2 coordinates x and y) and a class describing a Polygon which has a list of Points which correspond to corners (self.corners) I need to check if a Point is in a Polygon

Here is the function that is supposed to check if the Point is in the Polygon. I am using the Ray Casting Method

def in_me(self, point):
        result = False
        n = len(self.corners)
        p1x = int(self.corners[0].x)
        p1y = int(self.corners[0].y)
        for i in range(n+1):
            p2x = int(self.corners[i % n].x)
            p2y = int(self.corners[i % n].y)
            if point.y > min(p1y,p2y):
                if point.x <= max(p1x,p2x):
                    if p1y != p2y:
                        xinters = (point.y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
                        print xinters
                    if p1x == p2x or point.x <= xinters:
                        result = not result
            p1x,p1y = p2x,p2y
         return result

I run a test with following shape and point:

PG1 = (0,0), (0,2), (2,2), (2,0)
point = (1,1)

The script happily returns False even though the point it within the line. I am unable to find the mistake

5
  • 1
    Might be because you're using "/" on integers, which returns an integer (rounded down). You should do all computations with floats instead. Also, if p1y == p2y, xinters might not be defined but still used just afterwards.
    – Armin Rigo
    May 18, 2013 at 15:04
  • Better yet: don't divide at all. Instead of computing xinters, check if (point.x - p1x)*(p2y-p1y) <= (point.y-p1y)*(p2x-p1x). However, casting the vertex coordinates to integers could introduce errors if they aren't already integers to start with.
    – chepner
    May 18, 2013 at 16:19
  • 1
    ...or use Python 3, which doesn't truncate to integers on division. May 18, 2013 at 16:56
  • how would using (point.x - p1x)*(p2y-p1y) <= (point.y-p1y)*(p2x-p1x) make the actual code look like? Since it is a homework assignment, then we have to use Python 2.7 :(
    – Helena
    May 18, 2013 at 17:41
  • @Ulrich & helena: Python 3 division can be enabled in Python 2 using from __future__ import division. Another alternative is to just float() either the numerator or denominator (or a term in one of them in this case).
    – martineau
    Jul 14, 2014 at 13:22

4 Answers 4

65

I would suggest using the Path class from matplotlib

import matplotlib.path as mplPath
import numpy as np

poly = [190, 50, 500, 310]
bbPath = mplPath.Path(np.array([[poly[0], poly[1]],
                     [poly[1], poly[2]],
                     [poly[2], poly[3]],
                     [poly[3], poly[0]]]))

bbPath.contains_point((200, 100))

(There is also a contains_points function if you want to test for multiple points)

7
  • 3
    For this to work, you must first import numpy as np Mar 6, 2015 at 16:54
  • 5
    Anyone checked performance of contains_points against a pure Python implementation ? Jul 11, 2016 at 8:18
  • 1
    Something's wrong, using array = [[100,100],[200,100],[200,200],[100,200],[100,100]] it gives False for point 100,100 and true for point 200,200
    – Maciek
    Mar 10, 2018 at 12:00
  • 1
    Why the variable name 'bbPath'? if (Does 'bb' abbreviate something?): what does 'bb' abbreviate?
    – nda
    Feb 7, 2019 at 2:46
  • 1
    bb means bounding box even though the polygon very like wont be a box :)
    – P.R.
    Feb 11, 2019 at 12:07
3

I'd like to suggest some other changes there:

def contains(self, point):
    if not self.corners:
        return False

    def lines():
        p0 = self.corners[-1]
        for p1 in self.corners:
            yield p0, p1
            p0 = p1

    for p1, p2 in lines():
        ... # perform actual checks here

Notes:

  • A polygon with 5 corners also has 5 bounding lines, not 6, your loop is one off.
  • Using a separate generator expression makes clear that you are checking each line in turn.
  • Checking for an empty number of lines was added. However, how to treat zero-length lines and polygons with a single corner is still open.
  • I'd also consider making the lines() function a normal member instead of a nested utility.
  • Instead of the many nested if structures, you could also check for the inverse and then continue or use and.
1

I was trying to solve the same problem for my project and I got this code from someone in my network.

#!/usr/bin/env python
#
# routine for performing the "point in polygon" inclusion test

# Copyright 2001, softSurfer (www.softsurfer.com)
# This code may be freely used and modified for any purpose
# providing that this copyright notice is included with it.
# SoftSurfer makes no warranty for this code, and cannot be held
# liable for any real or imagined damage resulting from its use.
# Users of this code must verify correctness for their application.

# translated to Python by Maciej Kalisiak <[email protected]>

#   a Point is represented as a tuple: (x,y)

#===================================================================

# is_left(): tests if a point is Left|On|Right of an infinite line.

#   Input: three points P0, P1, and P2
#   Return: >0 for P2 left of the line through P0 and P1
#           =0 for P2 on the line
#           <0 for P2 right of the line
#   See: the January 2001 Algorithm "Area of 2D and 3D Triangles and Polygons"

def is_left(P0, P1, P2):
    return (P1[0] - P0[0]) * (P2[1] - P0[1]) - (P2[0] - P0[0]) * (P1[1] - P0[1])

#===================================================================

# cn_PnPoly(): crossing number test for a point in a polygon
#     Input:  P = a point,
#             V[] = vertex points of a polygon
#     Return: 0 = outside, 1 = inside
# This code is patterned after [Franklin, 2000]

def cn_PnPoly(P, V):
    cn = 0    # the crossing number counter

    # repeat the first vertex at end
    V = tuple(V[:])+(V[0],)

    # loop through all edges of the polygon
    for i in range(len(V)-1):   # edge from V[i] to V[i+1]
        if ((V[i][1] <= P[1] and V[i+1][1] > P[1])   # an upward crossing
            or (V[i][1] > P[1] and V[i+1][1] <= P[1])):  # a downward crossing
            # compute the actual edge-ray intersect x-coordinate
            vt = (P[1] - V[i][1]) / float(V[i+1][1] - V[i][1])
            if P[0] < V[i][0] + vt * (V[i+1][0] - V[i][0]): # P[0] < intersect
                cn += 1  # a valid crossing of y=P[1] right of P[0]

    return cn % 2   # 0 if even (out), and 1 if odd (in)

#===================================================================

# wn_PnPoly(): winding number test for a point in a polygon
#     Input:  P = a point,
#             V[] = vertex points of a polygon
#     Return: wn = the winding number (=0 only if P is outside V[])

def wn_PnPoly(P, V):
    wn = 0   # the winding number counter

    # repeat the first vertex at end
    V = tuple(V[:]) + (V[0],)

    # loop through all edges of the polygon
    for i in range(len(V)-1):     # edge from V[i] to V[i+1]
        if V[i][1] <= P[1]:        # start y <= P[1]
            if V[i+1][1] > P[1]:     # an upward crossing
                if is_left(V[i], V[i+1], P) > 0: # P left of edge
                    wn += 1           # have a valid up intersect
        else:                      # start y > P[1] (no test needed)
            if V[i+1][1] <= P[1]:    # a downward crossing
                if is_left(V[i], V[i+1], P) < 0: # P right of edge
                    wn -= 1           # have a valid down intersect
    return wn
0

Steps:

  • Iterate over all the segments in the polygon
  • Check whether they intersect with a ray going in the increasing-x direction

Using the intersect function from This SO Question

def ccw(A,B,C):
    return (C.y-A.y) * (B.x-A.x) > (B.y-A.y) * (C.x-A.x)

# Return true if line segments AB and CD intersect
def intersect(A,B,C,D):
    return ccw(A,C,D) != ccw(B,C,D) and ccw(A,B,C) != ccw(A,B,D)

def point_in_polygon(pt, poly, inf):
    result = False
    for i in range(len(poly.corners)-1):
        if intersect((poly.corners[i].x, poly.corners[i].y), ( poly.corners[i+1].x, poly.corners[i+1].y), (pt.x, pt.y), (inf, pt.y)):
            result = not result
    if intersect((poly.corners[-1].x, poly.corners[-1].y), (poly.corners[0].x, poly.corners[0].y), (pt.x, pt.y), (inf, pt.y)):
        result = not result
    return result

Please note that the inf parameter should be the maximum point in the x axis in your figure.

1
  • This is incorrect, doesn't work for point [2, 5] with polygon [8, 6], [11, 10], [16, 5], [11, 3] Edit: The issue is probably that the ray goes directly through a point of the polygon, causing two polygon line segments to be toggling result, turning it back to its previous state
    – h345k34cr
    Jul 4, 2017 at 12:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.