CODEWITHSUNDEEP
pythonnumpy2dpolygon
IssamLaradji
IssamLaradji

Reputation: 6865

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

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.

Upvotes: 2

Views: 1636

Answers (2)

Gerges
Gerges

Reputation: 6509

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:

  1. Find equations of lines of all edges

  2. Find bounding box

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

  4. 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)

Upvotes: 1

HYRY
HYRY

Reputation: 97301

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]

Upvotes: 4

Related Questions