Intersections : Polygon overlay operations
- the main polygon in blue,
- the "clipper" in red
In order to perform boolean operations on polygons you need to determine where the overlaying geometries intersect.
Far from simple. Countless algorithms have been proposed over the years and each have their advantages and disadvantages. Most work by cycling through the boundary edges and determining intersection points one at a time. I like NumPy and I was sure that it could be put into action to determine where the intersection points were in a quasi-vectorized fashion while at the same time constructing the associated information as to what crossed what and where and how many times.
We can start. Here are the coordinates, ordered clockwise with the first and last points being the same.
# -- polygon
array([[ 0.00, 0.00],
[ 2.00, 8.00],
[ 8.00, 10.00],
[ 10.00, 10.00],
[ 10.00, 8.00],
[ 9.00, 1.00],
[ 0.00, 0.00]])
# -- clipper
array([[ 0.00, 7.50],
[ 11.00, 10.50],
[ 12.00, 4.50],
[ 8.00, 0.00],
[ 0.00, 7.50]])
I use my 'winding number' incarnation to determine points inside/outside of their respective geometries... that is polygon points within the clipper, and the converse, clipper points within the polygon. This information gives rise to information on which polygon edges/segments are completely inside, or cross or are completely outside. From this intersections are determined.
Details on the equations can be found in Crossings : line segment intersections in 2D
The winding number (_w_) and intersection (_xsect_) functions are included below.
The array of interest is the both array
Rows are clipper segments, columns are polygon segments.
0 1 2 3 4 5
[[0, 1, 1, 0, 0, 0] 0
[0, 0, 0, 0, 0, 0] 1
[0, 0, 0, 0, 1, 1] 2
[1, 0, 0, 0, 0, 1] 3
This gives rise to the whr array which provides the intersection pairs
whr
array([[0, 1],
[0, 2],
[2, 4],
[2, 5],
[3, 0],
[3, 5]])
And ultimately, the intersection points themselves.
x_pnts
array([[ 2.75, 8.25],
[ 9.17, 10.00],
[ 9.02, 1.15],
[ 8.88, 0.99],
[ 1.52, 6.08],
[ 7.15, 0.79]])
def _w_(a, b, all_info):
"""Return winding number and other values."""
x0, y0 = a[:-1].T # point `from` coordinates
x1, y1 = a[1:].T # point `to` coordinates
x1_x0, y1_y0 = (a[1:] - a[:-1]).T
#
x2, y2 = b[:-1].T # polygon `from` coordinates
x3, y3 = b[1:].T # polygon `to` coordinates
x3_x2, y3_y2 = (b[1:] - b[:-1]).T
# reshape poly deltas
x3_x2 = x3_x2[:, None]
y3_y2 = y3_y2[:, None]
# deltas between pnts/poly x and y
x0_x2 = (x0 - x2[:, None])
y0_y2 = y0 - y2[:, None]
#
a_0 = y0_y2 * x3_x2
a_1 = x0_x2 * y3_y2
b_0 = x1_x0 * y0_y2
b_1 = y1_y0 * x0_x2
a_num = a_0 - a_1
b_num = b_0 - b_1
#
# pnts in poly
chk1 = (y0_y2 >= 0.0) # y above poly's first y value, per segment
chk2 = np.less(y0, y3[:, None]) # y above the poly's second point
chk3 = np.sign(a_num).astype(np.int32)
pos = (chk1 & chk2 & (chk3 > 0)).sum(axis=0, dtype=np.int32)
neg = (~chk1 & ~chk2 & (chk3 < 0)).sum(axis=0, dtype=np.int32)
wn_vals = pos - neg
wn_ = np.concatenate((wn_vals, np.array([wn_vals[0]])))
#
if all_info:
denom = (x1_x0 * y3_y2) - (y1_y0 * x3_x2)
return wn_, denom, x0, y0, x1_x0, y1_y0, a_num, b_num
return wn_ def _xsect_(a_num, b_num, denom, x1_x0, y1_y0, x0, y0):
"""Return the intersection."""
with np.errstate(all="ignore"): # ignore all errors
u_a = a_num / denom
u_b = b_num / denom
z0 = np.logical_and(u_a >= 0., u_a <= 1.) # np.isfinite(u_a)`
z1 = np.logical_and(u_b >= 0., u_b <= 1.) # np.isfinite(u_b)
both = z0 & z1
xs = (u_a * x1_x0 + x0)[both]
ys = (u_a * y1_y0 + y0)[both]
x_pnts = []
if xs.size > 0:
x_pnts = np.concatenate((xs[:, None], ys[:, None]), axis=1)
# eq = pnts[eq_]
whr = np.array(np.nonzero(both)).T
return whr, x_pntsOut there for posterity and out of my system 😉