Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Find the intersection of two polygons in NumPy using the Sutherland–Hodgman algorithm
# 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment