Searching the internet, I found an interesting and (IMO) fundamental
solution to the problem of point in polygon, see "Fast Winding Number
Inclusion of a Point in a Polygon" on
http://www.softsurfer.com/Archive/algorithm_0103/algorithm_0103.htm .
( http://www.softsurfer.com/algorithms.htm is a real treasure trove for
geometric algorithms.)
Some quotes:
- So, for both geometric correctness and efficiency reasons, the wn
algorithm should always be preferred for determining the inclusion of a
point in a polygon.
- Further, one must decide whether a point on the polygon's boundary is
inside or outside. A standard convention is to say that a point on a left
or bottom edge is inside, and a point on a right or top edge is outside.
This way, if two distinct polygons share a common boundary segment, then a
point on that segment will be in one polygon or the other, but not both at
the same time. This avoids a number of problems that might occur especially
in computer graphics displays.
- It is efficient to first test that a point is inside the bounding box or
ball (....)
I translated the wn_PnPoly() function to
wn=: 4 : 0 NB. x is polygon, y is points
XX=. |:(, {.)x [ Y=. |: y
in=. *./ Y ((>: {.) *. (<: {:))"1 (<./,>./)"1 XX
'Xy Yy'=. XX ;&{: Yin=. in #"1 Y
T =. Xy <:/ Yy
'A B'=. (}: ,: }.) T
'p n'=. (A * -.B) ,: (-.A) * B
S=. |:(t0=. $A) #: s=. I., p +. n
R=. * (({.S) {"1 (}. - }:)"1 XX) ([: -/ (*|.)) (({:S) {"1 Yin) - ({.S) {"1
XX
t=. t0 $ R s} t=.0 #~ */t0
in=. (+/ (p * 0 < t) - n * 0 > t) (I.1=in)} in
)
The second line of wn ( in=: ...) uses the bounding box of minimal and
maximal coordinates of the polygon. Only for points within that box the
algorithm is applied.
>From the second quote it is hard to compare the outcome of wn with other
solutions.
PLG=: +. sprl=: +/\(k...@-.,1,k=.(*0j1&^)@>:&i.)50
NB. see for sprl
http://jsoftware.com/pipermail/programming/2010-April/019376.html
PTS=: 50 * (3 :'0 ?...@$~ 2 ,~ 10^y') 4
NB. 50 since
NB. (<./,>./)"1 |: +.sprl
NB. _1 49
NB. 1 50
('wn';'W2a') dspl rnkng scores 'PLG wn PTS';'PLG W2a PTS'
+----+----+-----+----+----+
|verb|rank|et*sz|time|size|
+----+----+-----+----+----+
|wn | 0 |1.00 |1.00|1.00|
+----+----+-----+----+----+
|W2a | 1 |9.82 |5.54|1.77|
+----+----+-----+----+----+
(PLG wn PTS)-:PLG W2a PTS
1
If the majority of points are outside the polygon, the improvement is even
larger:
('wn';'W2a') dspl rnkng scores 'PLG wn 3* PTS';'PLG W2a 3* PTS'
+----+----+-----+-----+----+
|verb|rank|et*sz|time |size|
+----+----+-----+-----+----+
|wn | 0 | 1.00| 1.00|1.00|
+----+----+-----+-----+----+
|W2a | 1 |46.96|10.68|4.40|
+----+----+-----+-----+----+
PLG (wn -: W2a) 3* PTS
1
I don't think Bron will get a much more efficient solution.
Apart from that, it works for all polygons.
R.E. Boss
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
