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.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
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.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
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 < (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
