geometric utilities

[mikeg2@my-dejanews.com]

OK, I'll do this in a couple of pieces.  First, below is an excerpt from
Robert Sedgewick's "Algorithms", 2nd Edition, Addison Wesley, 1988, describing
the point-in-polygon algorithm that I'm trying to use.  BTW, in answer to
Brian's questions - floating point (the application this is for deals with
regions defined via lat/long points, which I convert to float), must work on a
non-convex polygon, and the algorithm does use a ray-crossing technique.

Quoting Sedgewick, starting on page 353 (this is a little long, but i'm
including it all so you'll have what i have):

... A straightforard solution to this problem immediately suggests itself:
draw a long line segment from the point in any direction (long enough so that
its other endpoint is guaranteed to be outside the polygon) and count the
number of lines from the polygon that it crosses.  If the number is odd, the
point must inside; if it is even, the point is outside.  ...

The situation is not quite so simple, however, because some intersections
might occur right at the vertices of the input polygon. ...

The need to handle cases where polygon vertices fall on the test lines forces
us to do more than just count the line segments in the polygon intersecting
the test line.  Essentially, we want to travel around the polygon,
incrementing the intersection counter whenever we go from one side of the test
line to another.  One way to implement this is to simply ignore points that
fall on the test line, as in the following program:

01  function inside(t: point): boolean;
02    var count,i,j: integer;
03        lt,lp: line;
04    begin
05    count:=0;  j:=0;
06    p[0]:=p[N];p[N+1]:=p[1];
07    lt.p1:=t; lt.p2:=t; lt.p2.x:=maxint;
08    for i := 1 to N do
09      begin
10      lp.p1:=p[i]; lp.p2:=p[i];
11      if not intersect (lp,lt) then
12        begin
13        lp.p2:=p[j]; j:=i;
14        if intersect (lp,lt) then count:=count+1;
15        end;
16      end;
17    inside := ((count mode 2)=1);
18    end;

This program uses a horizontal test line for ease of calculation. ... The
variable j is maintained as the index of the last point on the polygon known
not to lie on the test line.  The program assumes that p[1] is the point with
the smallest x coordinate among all the points with the smallest y coordinate,
so that if p[1] is on the test line, then p[0] cannot be.  The same polygon
can be represented by N different p arrays, but as this illustrates it is
sometimes convenient to fix a standard rule for p[1].  ... If the next point
on the polygon that is not on the test line is on the same side of the test
line as the jth point, then we need not increment the intersection counter
(count); otherwise, we have an intersection. ...

--------------------- end of quoted material ---------------------------

I left the code formatted exactly as it appears in the text - hey, it was
written in 1988.

p (the polygon) and N (the 'length of the polygon) are global.  Big
assumption about the ordering of points in p - I had to write extra code
to deal with the notion of p[1] is the point with the smallest x coordinate
among the points with the smallest y coordinate.

Note on line 10 lp is set to essentially a point - maybe that's a bug?
lp.p2 := p[i+1]?  I'll go test that.  My code to follow as soon as I have a
chance.


-----== Posted via Deja News, The Leader in Internet Discussion ==-----
http://www.dejanews.com/   Now offering spam-free web-based newsreading