#!/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)