Point in polygon ray casting algorithm

Question

Given a matrix of size M, N and a list of points form a polygon, I need to return a binary matrix where it marks all points inside the polygon. Here is the code: using a horizontal line scanning, intersect with polygon and mark all the points between intersection to be true.

import cv2
import math
import numpy as np

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):
    """
    Return a mask matrix
    Where points inside the polygon is 1, outside is 0
    """
    out = np.zeros((M, N)). astype(bool)

    n = len(poly)
    for i in range(M):

        # horizontal scanning
        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]

        v1_x = int(v1_x)
        v1_y = int(v1_y)
        v2_x = int(v2_x)
        v2_y = int(v2_y)

            assert (v1_x <= M and v2_x <= M and v1_y <= N and v2_y <= N)

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

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

        intersection_y = sorted(list(set(intersection_y)))



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

    return out

if __name__ == "__main__":
out = polygon(512, 512,         
    [[
      321.929203539823, 
      414.4070796460177
    ], 
    [
      164.74414597854175, 
      512.0
    ], 
    [
      0.0, 
      509.9846825337252
    ], 
    [
      0.0, 
      221.4867256637168
    ], 
    [
      170.60176991150445, 
      2.902654867256672
    ], 
    [
      320.1592920353982, 
      91.39823008849561
    ], 
    [
      271.4867256637168, 
      201.1327433628319
    ], 
    [
      348.4778761061947, 
      228.56637168141594
    ], 
    [
      359.98230088495575, 
      302.9026548672567
    ], 
    [
      220.15929203539824, 
      329.4513274336283
    ]])

    cv2.imwrite("polygon.png", out * 255)

It performs wrong in the case of a corner, for example the input above will generate:

I haven't been able to fix this problem. Depending on the number of intersections, I will assign l0-l1 l2-l3 l4-l5 or l1-l2 l3-l4 l5-l6 to be true.

Answer 1

You can fix by applying the following algorithm:
1. If a line intersects an edge inside it(not in one of the end points), always keep the point of intersection.
2. If a line intersects an edge in an end point, keep it if and only if it the left end of the edge(with smaller x-coordinate).
3. If an edge is vertical, just ignore it.

After that, just sort the points without removing duplicates.