#!/usr/bin/env python3
'''Class to compute if a point lies inside/outside/on-side of a polygon.
This is a Python 3 implementation of the Sloan's improved version of the
Nordbeck and Rystedt algorithm, published in the paper:
SLOAN, S.W. (1985): A point-in-polygon program.
Adv. Eng. Software, Vol 7, No. 1, pp 45-47.
This class has 1 public method (is_inside) that returns the minimum distance to
the nearest point of the polygon:
If is_inside < 0 then point is outside the polygon.
If is_inside = 0 then point in on a side of the polygon.
If is_inside > 0 then point is inside the polygon.
'''
import math
import numpy as np
class Polygon:
'''Polygon object.
Input parameters:
x -- A sequence of nodal x-coords.
y -- A sequence of nodal y-coords.
'''
def __init__(self, x, y):
if len(x) != len(y):
raise IndexError('x and y must be equally sized.')
self.x = np.asfarray(x)
self.y = np.asfarray(y)
# Closes the polygon if needed
x1, y1 = x[0], y[0]
xn, yn = x[-1], y[-1]
if x1 != xn or y1 != yn:
self.x = np.concatenate((self.x, [x1]))
self.y = np.concatenate((self.y, [y1]))
# Transform to anti-clockwise if needed
if self._det(self.x, self.y) < 0:
self.x = self.x[::-1]
self.y = self.y[::-1]
@staticmethod
def _det(xvert, yvert):
'''Compute twice the area of the triangle defined by points with using
determinant formula.
Input parameters:
xvert -- A vector of nodal x-coords.
yvert -- A vector of nodal y-coords.
Output parameters:
Twice the area of the triangle defined by the points.
Notes:
_det is positive if points define polygon in anticlockwise order.
_det is negative if points define polygon in clockwise order.
_det is zero if at least two of the points are concident or if
all points are collinear.
'''
xvert = np.asfarray(xvert)
yvert = np.asfarray(yvert)
x_prev = np.concatenate(([xvert[-1]], xvert[:-1]))
y_prev = np.concatenate(([yvert[-1]], yvert[:-1]))
return np.sum(yvert * x_prev - xvert * y_prev)
def is_inside(self, xpoint, ypoint, smalld=1e-12):
'''Check if point is inside a general polygon.
Input parameters:
xpoint -- The x-coord of the point to be tested.
ypoint -- The y-coord of the point to be tested.
smalld -- A small float number.
Output parameters:
mindst -- The distance from the point to the nearest point of the
polygon.
If mindst < 0 then point is outside the polygon.
If mindst = 0 then point in on a side of the polygon.
If mindst > 0 then point is inside the polygon.
Notes:
An improved version of the algorithm of Nordbeck and Rydstedt.
REF: SLOAN, S.W. (1985): A point-in-polygon program. Adv. Eng.
Software, Vol 7, No. 1, pp 45-47.
'''
# If snear = True: Dist to nearest side < nearest vertex
# If snear = False: Dist to nearest vertex < nearest side
x = self.x
y = self.y
n = len(x) - 1 # Number of sides/vertices defining the polygon
mindst = None
# Loop over each side defining polygon
for i in range(n):
# Start of side has coords (x1, y1)
# End of side has coords (x2, y2)
# Point has coords (xpoint, ypoint)
x1 = x[i]
y1 = y[i]
x21 = x[i + 1] - x1
y21 = y[i + 1] - y1
x1p = x1 - xpoint
y1p = y1 - ypoint
# Points on infinite line defined by
# x = x1 + t * (x1 - x2)
# y = y1 + t * (y1 - y2)
# where
# t = 0 at (x1, y1)
# t = 1 at (x2, y2)
# Find where normal passing through (xpoint, ypoint) intersects
# infinite line
t = -(x1p * x21 + y1p * y21) / (x21 ** 2 + y21 ** 2)
if t < 0:
# Normal does not intersects side
# Point is closest to vertex (x1, y1)
# Compute square of distance to this vertex
d = x1p ** 2 + y1p ** 2
if mindst is None or d < mindst:
# Point is closer to (x1, y1) than any other vertex or side
snear = False
mindst = d
j = i
elif t <= 1:
# Normal intersects side
dx = x1p + t * x21
dy = y1p + t * y21
d = dx ** 2 + dy ** 2
if mindst is None or d < mindst:
# Point is closer to this side than to any other side or
# vertex
snear = True
mindst = d
j = i
mindst **= 0.5
if mindst < smalld:
# Point is on side of polygon
mindst = 0
elif snear:
# Point is closer to its nearest side than to its nearest vertex,
# check if point is to left or right of this side.
# If point is to left of side it is inside polygon, else point is
# outside polygon.
area = self._det([x[j], x[j + 1], xpoint],
[y[j], y[j + 1], ypoint])
mindst = math.copysign(mindst, area)
else:
# Point is closer to its nearest vertex than its nearest side,
# check if nearest vertex is concave.
# If the nearest vertex is concave then point is inside the
# polygon, else the point is outside the polygon.
if not j:
x = x[:-1]
y = y[:-1]
area = self._det([x[j + 1], x[j], x[j - 1]],
[y[j + 1], y[j], y[j - 1]])
mindst = math.copysign(mindst, area)
return mindst
# TEST
if __name__ == '__main__':
# Define a triangle
xvert = [1, 10, 10]
yvert = [1, 10, 1]
poly = Polygon(xvert, yvert)
# Test the function in every point of a 12x12 grid
grid = np.zeros((12, 12), dtype=int)
for y in range(12):
for x in range(12):
if poly.is_inside(x, y) >= 0:
grid[y, x] = 1
# Print the result (0=outside)
print(grid)
