numpy: get the indices within a polygon without creating a mask

Question

I have a polygon vector formatted as follows (x1,y1,x2,y2, ...., xn,yn). As an example, consider this polygon array: polyPoints = [3,5,7,8,9,5]

How do I get all the indices (or the coordinates) that are within the polygon generated from these points ?

The answers I looked so far requires you to create the 2D mask before you can get the indices within the polygon.

Straightforward brute force way will probably be to iterate on all integer points from min(x), min(y) to max(x), max(y) and for each point test if inside or outside the polygon. To do the test for each point use raycasting or winding number algorithm. See here for these tests: https://en.wikipedia.org/wiki/Point_in_polygon — Aguy, Feb 07 '18 at 03:09

score 4 · Accepted Answer · answered Feb 07 '18 at 04:37

4

You can use scikit-image:

import numpy as np
from skimage.draw import polygon
points = [3,5,7,8,9,5]
r, c = polygon(points[1::2], points[::2])
print(r, c)

the output is:

[5 5 5 5 5 5 6 6 6 6 7 7] [3 4 5 6 7 8 5 6 7 8 6 7]

answered Feb 07 '18 at 04:37

HYRY

94,853
25
187
187

score 1 · Answer 2 · answered Feb 07 '18 at 04:37

Using a mask is probably as efficient as you can get. This is some algorithm you which is rather inefficient, but probably can be optimized to be close to the mask approach. This essentially does a mask but on lines.

The approach is:

Find equations of lines of all edges
Find bounding box
For each y within bounding box (or x, whichever is smaller), compute the edges which intersect with the horizontal line (y=yi) at that y, and find at which x they intersect.
For each x within the bounding box, find which number of edges to the right of x which the line y=yi intersect. If the number of edges is odd, then the point (x,y) is inside the polygon.

It does work on a simple square geometry.

import numpy as np

# taken from: https://stackoverflow.com/questions/20677795/how-do-i-compute-the-intersection-point-of-two-lines-in-python
def line(p1, p2):
    A = (p1[1] - p2[1])
    B = (p2[0] - p1[0])
    C = (p1[0]*p2[1] - p2[0]*p1[1])
    return A, B, -C


def intersection(L1, L2):
    D  = L1[0] * L2[1] - L1[1] * L2[0]
    Dx = L1[2] * L2[1] - L1[1] * L2[2]
    Dy = L1[0] * L2[2] - L1[2] * L2[0]
    if D != 0:
        x = Dx / D
        y = Dy / D
        return x,y
    else:
        return False


# polyPoints = np.array([0, 0, 4, 0,4, 4, 0, 4])
polyPoints = np.array([[3,5,7,8,9,5]])
polyPoints = polyPoints.reshape(-1, 2)
npoints = polyPoints.shape[0]
polyEgdes = []

for i in range(npoints):
    point1, point2 = polyPoints[i, :], polyPoints[(i+1) % npoints, :]
    polyEgdes.append(line(point1, point2))

# bounding box
boundingBox = np.vstack((polyPoints.min(axis=0), polyPoints.max(axis=0)))


inside_points = []

for y in range(boundingBox[0, 1], boundingBox[1, 1]):
    x_intersect = []
    for l in polyEgdes:
        # y_ins should be same as y
        insect_point = intersection(l, [0, y, 0])
        if insect_point:
            x_intersect.append(insect_point[0])

    x_intersect = np.array(x_intersect)
    for x in range(boundingBox[0, 0]+1, boundingBox[1, 0]-1):
        x_int_points = x_intersect[(x_intersect - x) >= 0]

        if len(x_int_points) % 2 == 1:
            inside_points.append((x, y))

print(inside_points)

numpy: get the indices within a polygon without creating a mask

2 Answers2