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

[python] 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.min_x or x > self.max_x or y < self.min_y 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

[/python]

Simple Test:

[python] 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() [/python]

Let's try to optimize it using numpy

[python] 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.min_x or x > self.max_x or y < self.min_y 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

[/python]

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:

[python]

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) [/python]

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 < (vert[j, 0] - vert[i, 0]) * (y - vert[i, 1]) / (vert[j, 1] - vert[i, 1]) + vert[i, 0])] = 1
PNPoly2: 10.531755924224854s
len(vertices): 3585
PNPoly: 127.7502794265747s
/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 < (vert[j, 0] - vert[i, 0]) * (y - vert[i, 1]) / (vert[j, 1] - vert[i, 1]) + vert[i, 0])] = 1
PNPoly2: 34.64983057975769s