geometric utilities

geometric utilities

[mikeg2]

Let me preface this by being very clear that I'm not
looking for a free ride on a homework assignment...

I'm trying to find some public domain Ada utilities
for some relatively common geometric problems - e.g.,
polygon area, polygon hull, point in polygon (in
particular), and polygon tesselation (I think that's
the spelling).  I've checked around in PAL and at the
AdaHome site, and not found anything.  Have been
trying to implement this with the aid of various
algorithms books, and, in fact, have pretty much
everything working - except point in polygon, which
doesn't keep the tesselation routine from working,
just keeps it from generating optimal solutions (i.e.,
fewest number of convex polygon partitions for a
non-convex polygon).

If anybody knows of or has some Ada code that does
this, please contact me:

  mike.glasgow@lmco.com
  (301) 640-2834

BTW, I have some C, C++, and Fortran versions - but
something just isn't coming out right in the translation
to Ada - the tricky part for point in polygon has to
do with edge vertices, and that's what's not working
for me.

Thanks.

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

Re: geometric utilities

[Brian Rogoff]

On Wed, 3 Jun 1998 mikeg2@my-dejanews.com wrote:
> I'm trying to find some public domain Ada utilities
> for some relatively common geometric problems - e.g.,
> polygon area, polygon hull, point in polygon (in
> particular), and polygon tesselation (I think that's
> the spelling).  I've checked around in PAL and at the
> AdaHome site, and not found anything.  Have been
> trying to implement this with the aid of various
> algorithms books, and, in fact, have pretty much
> everything working - except point in polygon, which
> doesn't keep the tesselation routine from working,
> just keeps it from generating optimal solutions (i.e.,
> fewest number of convex polygon partitions for a
> non-convex polygon).

Post your current solution, and someone will fix it. That's the beauty of 
the net :-).

Are you trying to implement Point_In_Polygon for non-convex or convex 
polygons? If non-convex, I assume you are using a ray crossing technique, 
and not winding number calculations which are slower, right?

Floating point or integer coordinates? 

I've implemented lots of geometric routines in C and C++ but not Ada.
It should be just the same.

> BTW, I have some C, C++, and Fortran versions - but
> something just isn't coming out right in the translation
> to Ada - the tricky part for point in polygon has to
> do with edge vertices, and that's what's not working
> for me.

If you post a clear description of what you want, I'll try to help.

Your post does bring up one of the problems of using Ada right now, that
is the lack of available source libraries in an easy to find form.

-- Brian

Re: geometric utilities

[Michael F Brenner]

I agree that you should post the non-working solution, because such
geometric utilities are interesting. Perhaps we could all
get together and post an entire geometric suite include both plane
and spherical Voronoi diagrams, and why stop at planes and spheres?

One of the things that is difficult in any language, but is now
possible in Ada-95 is to easily change between geometric spaces.
Such as switching smoothly between plane versus spherical trigonometry.

Actually, we could also switch into discrete topologies, building
tesselations of discrete spaces like braids, knots, etc. 

geometric utilities

[mikeg2]

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

Re: geometric utilities

[mikeg2]

Probably to my detriment, I've continued to try variations rather
than post what I've got so far - but, I'm really starting to
wonder if the algorithm actually works in all cases.

First, I think the lp.p1 := p[i] and lp.p2 := p[i] on line 10 (I
think that's where it was) is intentional per the text - i.e., ignore
*points* that are on the test line.  As it turns out, the intersect
procedure provided in the book doesn't appear to work in this case
(i.e., a point and a line) maybe because my stuff is floating point.
At any rate, I fixed this so that my code does properly detect
intersection of the test line with the ith point.

Using polygon ((1.0, 0.0), (0.0, 1.0), (-1.0, 0.0), (0.0, -1.0))
(a square around 0.0 on a 45 degree angle), and test point (0.0, 0.9),
the algorithm just doesn't seem to work, even when I do it by hand.
Maybe has something to do with j happening to start out at point
3 since point 4 passes the anchor rule (smallest x coordinate among
points with smallest y coordinate).  Point 3 happens to be collinear
with the test line...

Tomorrow I will post my code (not trying to be coy, I just have no way
here at home to get it from my laptop to my PC) - but was hoping in the
meantime someone might have taken a good look at the algorithm.  I did,
by the way, check with Addison Wesley for an errata sheet and sent a
note to Dr. Sedgewick to see if per chance there was a known problem
with the algorithm - but couldn't find anything and haven't yet heard
back from the good Doctor.

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

Re: geometric utilities

[Brian Rogoff]

On Thu, 4 Jun 1998 mikeg2@my-dejanews.com wrote:

> 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.

Mike, 
	The first problem is that you posted Pascal code to an Ada
newsgroup :-). 

	I just wrote down a naive Point_In_Polygon in Ada, untested (I'm 
Adaless for now). It should be easy to modify for whatever you're doing,
and its pretty self explanatory if you understand ray casting. It doesn't
have the same preconditions as what you posted. Try this if you can't
figure out the code you have. This is definitely hacked up on the fly,
and so serves as a counterexample that you can't hack in Ada ;-). 

-- Brian

------------------------------------------------------------------------
--| Geometry_2D.ads
------------------------------------------------------------------------
package Geometry_2D is 
    type Dimension_Type is range 1..2;
    X : constant := 1;
    Y : constant := 2;

    type Point_Type is array (Dimension_Type) of Float;
    type Polygon_Type is array (Positive range <>) of Point_Type;

    function X_Intersect( P1 : Point_Type; P2 : Point_Type ) return Float;

    function Point_In_Polygon ( Point  : Point_Type;
				Polygon : Polygon_Type )
			       return Boolean;
end Geometry_2D;

------------------------------------------------------------------------
--| Geometry_2D.adb
------------------------------------------------------------------------

package body Geometry_2D is
    procedure Increment( N : out Integer ) is
    begin
        N := N + 1;
    end Increment;

    function X_Intersect( P1 : Point_Type; P2 : Point_Type ) return Float
is
    begin
        return (P1(X)*P2(Y) - P2(X)*P1(Y)) / (P2(Y) - P1(Y));
    end X_Intersect;


    function "-" ( P1 : Point_Type; P2 : Point_Type ) return Point_Type is
        Result : Point_Type :=
          (P1(X) - P2(X), (P1(Y) - P2(Y)));
    begin
        return Result;
    end "-";

    function Point_In_Polygon ( Point   : Point_Type;
                                Polygon : Polygon_Type )
                               return Boolean is
        X_Intersection : Float;
        Edge_Point1, Edge_Point2 : Point_Type;
        J, Num_Intersections, Num_Points : Integer;
    begin
        Num_Intersections := 0;
        Num_Points := Polygon'Last;
        for I in 1 .. Num_Points loop
            J := (I + 1) mod Num_Points;

            -- Shift coordinate system so that Point is origin
            Edge_Point1 := Polygon(I) - Point;
            Edge_Point2 := Polygon(J) - Point;

            -- If edge has one point strictly above x-axis and one on or
            -- below, check for intersection. Note that while this
approach
            -- handles most of the degenerate cases, it means that points
on the
            -- boundary get "inconsistent" treatment. The inconsistency is
            -- practically useful though, in a region divided into
polygons each
            -- point will only be in one polygon.

            if (Edge_Point1(Y) >  0.0 and Edge_Point2(Y) <= 0.0) or
               (Edge_Point1(Y) <= 0.0 and Edge_Point2(Y) >  0.0) then
                X_Intersection := X_Intersect(Edge_Point1, Edge_Point2);
                if X_Intersection > 0.0 then
                    Increment(Num_Intersections);
                end if;
            end if;
        end loop;

        if Num_Intersections mod 2 = 1 then
            return True;
        else
            return False;
        end if;
    end Point_In_Polygon;
end Geometry_2D;

Re: geometric utilities

[mikeg2]

In article <Pine.BSF.3.96.980605150858.22568C-100000@shell5.ba.best.com>,
  Brian Rogoff <bpr@shell5.ba.best.com> wrote:

>
> Mike,
> 	The first problem is that you posted Pascal code to an Ada
> newsgroup :-).
>
> 	I just wrote down a naive Point_In_Polygon in Ada, untested (I'm
> Adaless for now). It should be easy to modify for whatever you're doing,
> and its pretty self explanatory if you understand ray casting. It doesn't
> have the same preconditions as what you posted. Try this if you can't
> figure out the code you have. This is definitely hacked up on the fly,
> and so serves as a counterexample that you can't hack in Ada ;-).
>
> -- Brian
>

thanks brian, your code seems to do the trick and was pretty easily integrated
with what i have.  i never could figure out what was wrong with the other
algorithm (sorry about the pascal, that's how it appeared in the textbook ;).

-mike

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