Point in a polygon algorithm

Question

I am implementing a Point in polygon algorithm.

Inputs:

• M, N: size of the matrix
• poly: a list of tuples that represent the points of polygon.

Output:

A mask matrix which is ones everywhere, except the point in the polygon should be 0.

import numpy as np
import cv2

def getABC(x1, y1, x2, y2):
    A = y2 - y1
    B = x1 - x2
    C = A*x1 + B*y1
    return (A, B, C)

def polygon(M, N, poly):
    out = np.ones((M, N))*255

    n = len(poly)
    for i in range(M):
        intersection_x = i
        intersection_y = []

        # check through all edges
        for edge in range(n + 1):
            v1_x, v1_y = poly[edge % n]
            v2_x, v2_y = poly[(edge + 1) % n]

            A1, B1, C1 = getABC(v1_x, v1_y, v2_x, v2_y)
            A2 = 1
            B2 = 0
            C2 = i

            # find intersection
            det = A1*B2 - A2*B1
            if (det != 0):
                tmp = (A1 * C2 - A2 * C1)/det
                if tmp >= min(v1_y, v2_y) and tmp <= max(v1_y, v2_y):
                    intersection_y.append(tmp)

    intersection_y = list(set(intersection_y))
    print intersection_y

    if len(intersection_y) == 1:
        intersection_y.append(intersection_y[0])

    for k in range(1, len(intersection_y), 2):
        out[intersection_x, intersection_y[k - 1]:intersection_y[k]] = 0

    return out



poly = [(10,20), (10,40), (30,20), (30,40)]

out = polygon(100,100, poly)

cv2.imwrite("out.png", out)

Is this the optimal algorithm? I use a scanning line, find the intersection with edges and make all the point between the intersections to be zero.

Answer 1

When you say "points of polygon", I am assuming you are referring to vertices.

I don't think what you have is the fastest. A simple improvement to this could be to divide your matrix in to a grid of p x p cells, where p is a parameter, and classify each grid-cell as completely inside or completely outside of the polygon. This can be done by checking if all the vertices of the grid-cell are inside the polygon. This way, large regions inside grid cells can be directly labeled as inside or outside, without needing to test every point inside the cell against all the edges of the polygon. For the grid-cells that are not completely inside the polygon, you just go by what you are doing. The optimal choice of p depends on the size of polygons you may encounter.

More details about the grid method here (with comparisons to other strategies): http://erich.realtimerendering.com/ptinpoly/

There could be better ones. I just suggested a simple improvement, which works well in practice.