Which side? Whose side? Left side? Right side?
The points? Both? Just one? Which one?
The segment?
To answer these questions observe the red and green point markers numbered in clockwise order for each of the two polygons.
Our brain can quickly decide which side IF we know the sidedness rule.
If the segments are ordered clockwise, then any point will be inside if it is to the right. If the polygon that we are testing as the home for the points is convex, then an inside point will be to the right of all the line segments.
This is the principle behind the winding number test for point inclusion. It, however, accounts for idiosyncracies that are associated with concave polygons.
Looking at sides
The _side_ code listed below, returns four values.
- r : the side array
- inside : the points inside (right of) a segment
- outside : the points outside (left of) a segment
- equal_ : any points that may intersect (equal to) one another by virtue of their equality or potential intersection on the segment.
# check b4 for inclusion in c1
b4 c1
array([[ 0.00, 10.00], array([[ 0.50, 10.50],
[ 1.00, 13.00], [ 0.50, 15.00],
[ 3.00, 14.00], [ 2.50, 13.50],
[ 2.50, 13.00], [ 6.50, 15.00],
[ 1.00, 11.50], [ 5.50, 10.50],
[ 0.00, 10.00]]) [ 2.00, 11.50],
[ 0.50, 10.50]])The coordinates for b4 are used here to just represent points. They actually form a polygon, yet they need not.
The c1 array represents a polygon, as per previous definition.
The _side_ function returns the set of values shown below. r is the cross-product array, which is used to determine parameters.
Two versions of determining points inside (in_ and inside) are calculated. in_ incorporates winding-number calculations for point
inclusion in polygons. It is interesting to examine under what circumstances these parameters are similar or different.
r, in_, inside, outside, equal_ = _side_(b4, c1)
r
array([[ 2.25, -10.75, -10.25, -24.25, 7.25, 0.25],
[ -2.25, -3.25, 0.25, -22.75, -4.25, -3.25],
[-11.25, 1.75, 1.25, -14.75, -9.75, -2.75],
[ -9.00, -1.00, -2.00, -16.00, -5.75, -1.75],
[ -2.25, -6.25, -5.75, -21.25, 1.00, -1.00],
[ 2.25, -10.75, -10.25, -24.25, 7.25, 0.25]])
(r < 0) r > 0
array([[0, 1, 1, 1, 0, 0], array([[1, 0, 0, 0, 1, 1],
[1, 1, 0, 1, 1, 1], [0, 0, 1, 0, 0, 0],
[1, 0, 0, 1, 1, 1], [0, 1, 1, 0, 0, 0],
[1, 1, 1, 1, 1, 1], [0, 0, 0, 0, 0, 0],
[1, 1, 1, 1, 0, 1], [0, 0, 0, 0, 1, 0],
[0, 1, 1, 1, 0, 0]]) [1, 0, 0, 0, 1, 1]])
inside outside equal_
array([[ 2.50, 13.00]]) array([[ 0.00, 10.00], array([],
[ 1.00, 13.00], shape=(0, 2),
_in_ [ 3.00, 14.00], dtype=float64)
array([[ 1.00, 13.00], [ 1.00, 11.50],
[ 2.50, 13.00], [ 0.00, 10.00]])
[ 1.00, 11.50]])
# -- calculate in_ using winding number
_wn_(b4, c1)
array([[ 1.00, 13.00],
[ 2.50, 13.00],
[ 1.00, 11.50]])
array([ 0, -1, 0, -1, -1, 0]))
#
lt0 = (r < 0).all(axis=-1) # array([0, 0, 0, 1, 0, 0])
gt0 = (r > 0).any(-1) # array([1, 1, 1, 0, 1, 1])
eq0 = (r == 0).any(-1) # array([0, 0, 0, 0, 0, 0])side the code
The number of points in each geometry object need not be the same since array-broadcasting is used to facilitate the calculations.
A term like y0[:, None] and other variants change the shape and dimension of the array as follows:
y0 # -- the base array
# y0.ndim : 1
# y0.shape : (6,)
array([ 10.00, 13.00, 14.00, 13.00, 11.50, 10.00])
# reshape to produce a column array.
y0[:, None] # ndim = 2, shape = (6, 1)
array([[ 10.00],
[ 13.00],
[ 14.00],
[ 13.00],
[ 11.50],
[ 10.00]])
# reshape to change the row array to a 2D array... note the additional [ ] around the coordinates
y0[None, :] # ndim = 2, shape = (1, 6)
array([[ 10.00, 13.00, 14.00, 13.00, 11.50, 10.00]])If you start with 2D arrays, then the transitions move them to 3D arrays.
pnts[:, None]
# ndim = 3
# shape = (6, 1, 2)
#
array([[[ 0.00, 10.00]],
[[ 1.00, 13.00]],
[[ 3.00, 14.00]],
[[ 2.50, 13.00]],
[[ 1.00, 11.50]],
[[ 0.00, 10.00]]])
pnts[None, :]
# ndim = 3
# shape = (1, 6, 2)
#
array([[[ 0.00, 10.00],
[ 1.00, 13.00],
[ 3.00, 14.00],
[ 2.50, 13.00],
[ 1.00, 11.50],
[ 0.00, 10.00]]])The principles of _side_ are used in subsequent code examples for polygon splitting, clipping and other overlay operations.
def _side_(pnts, poly):
"""Return points inside, outside or equal/crossing a convex poly feature."""
if pnts.ndim < 2:
pnts = np.atleast_2d(pnts)
x0, y0 = pnts.T # the x and y points
x2, y2 = poly[:-1].T # poly segment start points
x3, y3 = poly[1:].T # poly segment end points
r = (x3 - x2) * (y0[:, None] - y2) - (y3 - y2) * (x0[:, None] - x2)
# -- from _wn_, winding numbers for concave/convex poly
chk1 = ((y0[:, None] - y2) >= 0.)
chk2 = (y0[:, None] < y3)
chk3 = np.sign(r).astype(np.int)
pos = (chk1 & chk2 & (chk3 > 0)).sum(axis=1, dtype=int)
neg = (~chk1 & ~chk2 & (chk3 < 0)).sum(axis=1, dtype=int)
wn_vals = pos - neg
in_ = pnts[np.nonzero(wn_vals)]
inside = pnts[(r < 0).all(axis=-1)] # all True along row, convex polys only
outside = pnts[(r > 0).any(-1)] # any True along row
equal_ = pnts[(r == 0).any(-1)] # ditto
return r, in_, inside, outside, equal_