Trishock

It’s one of the simplest and most repeated operations in computational geometry: does a point fall within a polygon? Examples of this include assigning a point location to a state or county, potentially for filtering or to then gather information about the area the location falls within. A popular and simple method for solving this problem is the Ray Casting Algorithm which works by counting the number of intersections between a test line that contains the point in question and each edge of the polygon. If the number of intersections to the left of the test line and to the right of the test line are both odd, the point falls within the polygon.

As an exercise, I developed a solution to the point in polygon problem in Python. Most complicated is the test for two segments intersecting one another. I still barely understand the mathematics behind it, but the method I chose to use copy is done using vector cross products. See a thread on StackOverflow here for a detailed explanation of the methodology.

There are many special cases relating to where the point falls (exactly on an edge, exactly on the start or end point of an edge), how the test line is drawn (edge overlaps the test line), and the type of polygon (multi, self-intersecting). While the implementation I cobbled together below isn’t exhaustive, it’s functional for most cases.

def main():
    poly = Polygon()
    poly.AddPoint(Point(0,0))
    poly.AddPoint(Point(2,0))
    poly.AddPoint(Point(2,2))
    poly.AddPoint(Point(0,2))
    print(PointInPolygon(poly, Point(3,1)))
    print(PointInPolygon(poly, Point(1,1)))

def PointInPolygon(polygon, point):
    testline_left = Segment(Point(-999999999,point.y), point)
    testline_right = Segment(point, Point(-999999999,point.y))
    count_left = 0
    count_right = 0
    for e in polygon.GetEdges():
        if EdgesIntersect(testline_left, e):
            count_left += 1
        if EdgesIntersect(testline_right, e):
            count_right += 1
    if count_left % 2 == 0 and count_right % 2 == 0:
        return False
    else:
        return True

def EdgesIntersect(e1, e2):

    a = e1.p1
    b = e1.p2
    c = e2.p1
    d = e2.p2

    cmp = Point(c.x - a.x, c.y - a.y)
    r = Point(b.x - a.x, b.y - a.y)
    s = Point(d.x - c.x, d.y - c.y)

    cmpxr = cmp.x * r.y - cmp.y * r.x
    cmpxs = cmp.x * s.y - cmp.y * s.x
    rxs = r.x * s.y - r.y * s.x

    if cmpxr == 0:
        return (c.x - a.x < 0) != (c.x - b.x < 0)
    if rxs == 0:
        return False

    rxsr = 1 / rxs
    t = cmpxs * rxsr
    u = cmpxr * rxsr

    return t >= 0 and t <= 1 and u >= 0 and u <= 1

class Point:
    x = None
    y = None
    def __init__(self, x, y):
        self.x = x
        self.y = y

class Segment:
    p1 = None
    p2 = None
    def __init__(self, p1, p2):
        self.p1 = p1
        self.p2 = p2

class Polygon:
    points = None
    def __init__(self):
        self.points = []
    def AddPoint(self, p):
        self.points.append(p)
    def GetEdges(self):
        edges = []
        for i in range(len(self.points)):
            if i == len(self.points) - 1:
                i2 = 0
            else:
                i2 = i + 1
            edges.append(Segment(self.points[i], self.points[i2]))
        return edges

if __name__ == '__main__':
    main()

The output of this Python script is seen below.

False
True

For further reading, see:

• http://alienryderflex.com/polygon/
• http://en.wikipedia.org/wiki/Point_in_polygon