Find the intersection of two polygons in NumPy using the Sutherland–Hodgman algorithm

polygon-intersection.py

# Line intersection utility
# Works with broadcasting
def intersect_lines(line1, line2):
    """Computes the intersection of line1 and line2,
       both described by tuple of origin and direction vector"""
    p, r = line1
    q, s = line2
    rs = np.cross(r, s)[..., None]
    d = (q - p) / (rs + 1e-24)
    t = np.cross(d, s)[..., None]
    u = np.cross(d, r)[..., None]
    out = p + r * t
    hit = (np.abs(rs) > 0.) & (t >= 0) & (t <= 1.) & (u >= 0) & (u <= 1.)
    return out, hit[..., 0]

def norm_polygon(poly):
    """Returns a normalised polygon which is oriented clockwise"""
    # Check if polygon is ccw
    P0, P1 = poly, np.roll(poly, 1, 0)
    ccw = np.sum((P1[:, 0] - P0[:, 0]) * (P1[:, 1] + P0[:, 1])) < 0.
    if ccw:
      return np.roll(poly[::-1], 1, 0)
    return poly

# Optimized Sutherland–Hodgman algorithm
# https://en.wikipedia.org/wiki/Sutherland%E2%80%93Hodgman_algorithm
def clip_polygon(poly, clip_poly):
    """Computes the intersection of both polygons
         (requires normalized polygons, i.e. call
            norm_polygon beforehand)"""
    out = []
    clip_edges = np.roll(clip_poly, -1, 0) - clip_poly
    for clip_point, clip_direction in zip(clip_poly, clip_edges):
        out.clear()
        cur, prev = poly, np.roll(poly, 1, 0)
        # Find orientation of line
        orientation = np.sign(np.cross(clip_direction[None], cur - clip_point))
        inter, hit = intersect_lines(
            (prev, (cur - prev)),
            (clip_point[None], clip_direction[None]))
        # Apply algorithm
        for i in range(len(poly)):
            if orientation[i] > 0.:
                if orientation[i-1] <= 0.:
                    out.append(inter[i])
                out.append(cur[i])
            elif orientation[i-1] > 0.:
                out.append(inter[i])
        if len(out) <= 2:
            return None
        poly = np.stack(out, 0)
    return poly