@@ -51,11 +51,15 @@ function _inpolygon(::GI.AbstractPointTrait, points, ::GI.AbstractGeometryTrait,
5151 return first (inpoly2 ([(GI. x (points), GI. y (points))], nodes, edges; kw... ))
5252end
5353
54+
55+ using PolygonInbounds: vertex, edgecount, flipio!, edgeindex, searchfirst
56+
5457# Copied from PolygonInbounds, to add extra keyword arguments.
5558# PR to include these when this has solidified
5659function inpoly2 (vert, node, edge= zeros (Int);
5760 atol:: T = 0.0 , rtol:: T = NaN , iyperm= nothing ,
5861 vmin= nothing , vmax= nothing , pmin= nothing , pmax= nothing
62+
5963) where T<: AbstractFloat
6064 rtol = ! isnan (rtol) ? rtol : iszero (atol) ? eps (T)^ 0.85 : zero (T)
6165 poly = PolygonInbounds. PolygonMesh (node, edge)
@@ -70,7 +74,8 @@ function inpoly2(vert, node, edge=zeros(Int);
7074 tol = max (abs (rtol * lbar), abs (atol))
7175
7276 ac = PolygonInbounds. areacount (poly)
73- stat = ac > 1 ? falses (length (points), 2 , ac) : falses (length (points), 2 )
77+ # stat = ac > 1 ? fill(false, length(points), 1, ac) :
78+ stat = fill (false , length (points), 2 )
7479 # flip coordinates so expected effort is minimal
7580 dvert = vmax .- vmin
7681 if isnothing (iyperm)
@@ -83,3 +88,80 @@ function inpoly2(vert, node, edge=zeros(Int);
8388 PolygonInbounds. inpoly2! (points, iyperm, poly, ix, tol, stat)
8489 return stat
8590end
91+
92+ function inpoly2! (points, iyperm, poly, ix:: Integer , veps:: T , stat:: S ) where {N,T<: AbstractFloat ,S<: AbstractArray{Bool,N} }
93+ nvrt = length (points) # number of points to be checked
94+ nedg = edgecount (poly) # number of edges of the polygon mesh
95+ vepsx = vepsy = veps
96+ iy = 3 - ix
97+
98+ # ----------------------------------- loop over polygon edges
99+ for epos = 1 : nedg
100+
101+ inod = edgeindex (poly, epos, 1 ) # from
102+ jnod = edgeindex (poly, epos, 2 ) # to
103+ # swap order of vertices
104+ if vertex (poly, inod, iy) > vertex (poly, jnod, iy)
105+ inod, jnod = jnod, inod
106+ end
107+
108+ # ------------------------------- calc. edge bounding-box
109+ xone = vertex (poly, inod, ix)
110+ yone = vertex (poly, inod, iy)
111+ xtwo = vertex (poly, jnod, ix)
112+ ytwo = vertex (poly, jnod, iy)
113+
114+ xmin0 = min (xone, xtwo)
115+ xmax0 = max (xone, xtwo)
116+ xmin = xmin0 - vepsx
117+ xmax = xmax0 + vepsx
118+ ymin = yone - vepsy
119+ ymax = ytwo + vepsy
120+
121+ ydel = ytwo - yone
122+ xdel = xtwo - xone
123+ xysq = xdel^ 2 + ydel^ 2
124+ feps = sqrt (xysq) * veps
125+
126+ # find top points[:,iy] < ymin by binary search
127+ ilow = searchfirst (points, iy, iyperm, ymin)
128+ # ------------------------------- calc. edge-intersection
129+ # loop over all points with y ∈ [ymin,ymax)
130+ for jpos = ilow: nvrt
131+ jorig = iyperm[jpos]
132+ ypos = vertex (points, jorig, iy)
133+ ypos > ymax && break
134+ xpos = vertex (points, jorig, ix)
135+
136+ if xpos >= xmin
137+ if xpos <= xmax
138+ # --------- inside extended bounding box of edge
139+ mul1 = ydel * (xpos - xone)
140+ mul2 = xdel * (ypos - yone)
141+ # if abs(mul2 - mul1) <= feps
142+ # ------- distance from line through edge less veps
143+ # mul3 = xdel * (2xpos-xone-xtwo) + ydel * (2ypos-yone-ytwo)
144+ # if abs(mul3) <= xysq ||
145+ # hypot(xpos- xone, ypos - yone) <= veps ||
146+ # hypot(xpos- xtwo, ypos - ytwo) <= veps
147+ # # ---- round boundaries around endpoints of edge
148+ # setonbounds!(poly, stat, jorig, epos)
149+ # end
150+ # ----- left of line && ypos exact to avoid multiple counting
151+ # end
152+ if mul1 < mul2 && yone <= ypos < ytwo
153+ elseif mul1 < mul2 && yone <= ypos < ytwo
154+ # ----- left of line && ypos exact to avoid multiple counting
155+ flipio! (poly, stat, jorig, epos)
156+ end
157+ end
158+ else # xpos < xmin - left of bounding box
159+ if yone <= ypos < ytwo
160+ # ----- ypos exact to avoid multiple counting
161+ flipio! (poly, stat, jorig, epos)
162+ end
163+ end
164+ end
165+ end
166+ stat
167+ end
0 commit comments