PNPOLY: How do I find if a point lies within a polygon?

This is a python version to translate W. Randolph Franklin’s work. And recorded for further usage.

class PNPoly(object):
    """
    PNPOLY - Point Inclusion in Polygon Test
    W. Randolph Franklin (WRF)

    See: https://wrfranklin.org/Research/Short_Notes/pnpoly.html
    """

    def __init__(self, vertices) -> None:
        self.vertices = vertices
        assert len(vertices) >= 3, "vertices must be >= 3"
        self.min_x, self.max_x, self.min_y, self.max_y = self._make_border()

    def _make_border(self):
        vertices = self.vertices
        min_x, min_y = vertices[0]
        max_x, max_y = vertices[0]
        for i in range(1, len(vertices)):
            x, y = vertices[i]
            min_x = min(min_x, x)
            min_y = min(min_y, y)
            max_x = max(max_x, x)
            max_y = max(max_y, y)
        return min_x, max_x, min_y, max_y

    def __call__(self, x, y):
        return self.is_inside(x, y)

    def is_inside(self, x, y):
        """given a polygon and a point p, returns True if p is inside the polygon, False otherwise

        Args:
            x (int): x coordinate of the point
            y (int): y coordinate of the point
        """
        if x  self.max_x or y  self.max_y:
            return False
        vert = self.vertices
        nvert = len(vert)
        i = 0
        j = nvert - 1
        c = False
        for i in range(nvert):
            if ((vert[i][1] > y) != (vert[j][1] > y)) and (x < (vert[j][0] - vert[i][0]) * (y - vert[i][1]) / (vert[j][1] - vert[i][1]) + vert[i][0]):
                c = not c
            j = i
        return c

Simple Test:

def test():
    vertices = [(0,0), (2,0), (2,1), (1,1), (1,2), (2,2), (2,3), (0,3)]
    pnpoly = PNPoly(vertices)
    cases = [
        (0,0, True),
        (1.5, 0.5, True),
        (1.5, 1.5, False),
        (1.5, 2.5, True),
        (1.5, 3.5, False),
    ]
    for x, y, expected in cases:
        assert pnpoly(x, y) == expected, "x: {}, y: {}, expected: {}".format(x, y, expected)
        print("OK: x: {}, y: {}, expected: {}".format(x, y, expected))

test()

Let’s try to optimize it using numpy

import numpy as np

class PNPoly2(object):
    def __init__(self, vertices) -> None:
        assert len(vertices) >= 3, "vertices must be >= 3"
        self.vert = np.array(vertices)
        self.min_x, self.max_x, self.min_y, self.max_y = self._make_border()

    def _make_border(self):
        min_x, min_y = np.min(self.vert, axis=0)
        max_x, max_y = np.max(self.vert, axis=0)
        return min_x, max_x, min_y, max_y

    def __call__(self, x, y):
        return self.is_inside(x, y)

    def is_inside(self, x, y):
        """given a polygon and a point p, returns True if p is inside the polygon, False otherwise

        Args:
            x (int): x coordinate of the point
            y (int): y coordinate of the point
        """
        if x  self.max_x or y  self.max_y:
            return False
        # vert = np.array(self.vertices)
        vert = self.vert
        nvert = len(vert)
        i = np.arange(nvert)
        j = np.concatenate((i[1:], [0]))
        c = np.zeros(nvert, dtype=int)
        c[((vert[i, 1] > y) != (vert[j, 1] > y)) & (x < (vert[j, 0] - vert[i, 0]) * (y - vert[i, 1]) / (vert[j, 1] - vert[i, 1]) + vert[i, 0])] = 1
        return c.sum() % 2 == 1

According to the test, in the small size of the vertices, the overhead is too much ( np.arrange … blah operation). I have to generate much more vertices to let the average performance in numpy get better.

here is the test code:

def make_random_cases(seed):
    import random

    p = PNPoly(seed)
    min_x, max_x, min_y, max_y = p._make_border()
    rand_min_x = min_x - (max_x - min_x) * 0.5
    rand_max_x = max_x + (max_x - min_x) * 0.5
    rand_min_y = min_y - (max_y - min_y) * 0.5
    rand_max_y = max_y + (max_y - min_y) * 0.5
    for i in range(1000):
        x = random.uniform(rand_min_x, rand_max_x)
        y = random.uniform(rand_min_y, rand_max_y)
        expected = p.is_inside(x, y)
        yield x, y, expected

def add_mid_point(vertices):
    new_vertices = []
    for i in range(len(vertices)):
        vi = vertices[i]
        new_vertices.append(vi)
        if i == len(vertices) - 1:
            break
        vj = vertices[i + 1]
        if vi[0] == vj[0]:
            new_vertices.append((vi[0], (vi[1] + vj[1]) / 2))
        elif vi[1] == vj[1]:
            new_vertices.append(((vi[0] + vj[0]) / 2, vi[1]))
        else:
            new_vertices.append(((vi[0] + vj[0]) / 2, (vi[1] + vj[1]) / 2))
    return new_vertices

def test(n_more_verts):
    vertices = [(0, 0), (2, 0), (2, 1), (1, 1), (1, 2), (2, 2), (2, 3), (0, 3)]
    for i in range(n_more_verts):
        vertices = add_mid_point(vertices)
    print("len(vertices): {}".format(len(vertices)))
    test_cases = []
    for x, y, expected in make_random_cases(vertices):
        test_cases.append((x, y, expected))

    def make_test(pnpoly):
        for x, y, expected in test_cases:
            result = pnpoly(x, y)
            assert result == expected, "x: {}, y: {}, expected: {}, actually: {}".format(x, y, expected, result)

    time_in = time.time()
    for i in range(1000):
        pnpoly = PNPoly(vertices)
        make_test(pnpoly)
    time_out = time.time()
    print("PNPoly: {}s".format(time_out - time_in))

    time_in = time.time()
    for i in range(1000):
        pnpoly = PNPoly2(vertices)
        make_test(pnpoly)
    time_out = time.time()
    print("PNPoly2: {}s".format(time_out - time_in))

test(5)
test(9)

here is the test output

len(vertices): 225
PNPoly: 7.810054063796997s
/data/yu/aigc/pipe/monet/monet/mask_process/pnploy.py:85: RuntimeWarning: divide by zero encountered in divide
  c[((vert[i, 1] > y) != (vert[j, 1] > y)) & (x  y) != (vert[j, 1] > y)) & (x < (vert[j, 0] - vert[i, 0]) * (y - vert[i, 1]) / (vert[j, 1] - vert[i, 1]) + vert[i, 0])] = 1
PNPoly2: 34.64983057975769s