> Van: [email protected] [mailto:programming-
> [email protected]] Namens Bo Jacoby
> Onderwerp: Re: [Jprogramming] Polygon containment
>
> Dan wrote: I need to check whether a point falls within an arbitrary
> polygon.
>
> The winding number formula is: W=.9&o.@(%&0j2p1)@(+/)@:^.@(% _1&|.) . To
> test if (0,0) is within the square ((1,0),(0,1),(-1,0),(0,-1)), type W 1
> 0j1 _1 0j_1 and get the result 1, indicating YES. To test if (4,0) is
> within the square, type W 1 0j1 _1 0j_1 - 4 and get the result 0,
> indicating NO.
>
> (see
> http://en.wikipedia.org/wiki/Wikipedia:Reference_desk/Archives/Mathematics
> /2009_December_17#Easiest_way_to_decide_if_a_point_is_within_a_polygon.3F)
>
> -Bo
>
The basic idea of winding number is brilliant, but the J-implementation
needs some improvement.
First of all W ues log (^.)to determine the angle of a complex number where
J has a primitive to do that.
Furthermore for a point on a vertex you get a problem,
W 1 0j1 _1 0j_1 -1
|NaN error: W
| W 1 0j1 _1 0j_1-1
whereas the edges seems to be inside the polygon
W 1 0j1 _1 0j_1 - 2%~ 1j1
1
At last but certainly not at least, W has rank 0 which, for the " gazillion
points" of Bron, will give a worse performance then necessary.
So I constructed wn which is a dyad, accepting a polygon as x and the (2D)
points as y.
wn=: 4 : 0 NB. x is polygon such that 2[\(,{.)x are edges; y are points
NB. 1 if in polygon, 0 if on edge or vertex, _1 if outside
y=. ,:^:((2 > #...@$)`]) y
ov=. y e. x NB. on vertex
dv=. (,{.) j./"1 x -"1/ y NB. difference vectors
oe=. +./ 1e_14 ((<|)*]) (}. (= -) }:) dv NB. on edge
dv1=. dv #"1~ t=. -. ov +. oe
in=. t #^:_1 [0 ~: 1e_14 ((< |) * ]) (2*pi) %~ +/ {:"1 *. (}. % }:) dv1
<: (+: in) + ov +. oe
)
(1 0, 0 1, _1 0,: 0 _1) wn 2 %~ 1 1*/~i:2
_1 0 1 0 _1
(1 0, 0 1, _1 0,: 0 _1) wn 1 0*/~i:2
_1 0 1 0 _1
I estimate that this solution performs better than others.
Apart from that it works for a large(r) collection of polygons.
R.E. Boss
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
