Re: geometric utilities

[Brian Rogoff <bpr@shell5.ba.best.com>]

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;