From: "Geir Wikran" <gwikran@online.no>
Subject: Re: Again: PoinInPolygon
Date: 09 Oct 1999 00:00:00 GMT
Message-ID: <7tnube$7du10@forums.borland.com>
References: <7tkvu6$erb15@forums.borland.com>
Organization: Another Netscape Collabra Server User
X-MimeOLE: Produced By Microsoft MimeOLE V5.00.2615.200
Newsgroups: borland.public.delphi.graphics


> I have a pretty standard PointInPolygon-routine (got it from MSDN)

Gues you refere to
http://support.microsoft.com/support/kb/articles/Q121/9/60.asp

Here is my implementation of the function. I use a differen method for
calculating intersection points (see function LineInLine belowe for the
original line intersection calculations and testing). The function expects
an array of points (TPolygon) where the lines are going from point 0 to
point 1, point 1 to point 2, 2 to 3 ..., n-1 to n, n to 0.

The PointInPolygon function was done in a hurry, and has not been tested.
Feel free to send me a note if there are any questions.

type
  TPolygon = array[0..0] of TPoint;

function PointInPolygon(P: TPoint; Polygon: TPolygon; Points: Integer):
Integer;
  { Returns:  1 = point P inside polygon }
  {           0 = point P on edge of polygon }
  {          -1 = point P outside polygon }
var
  Index      : Integer;
  Count      : Integer;

  { Vars for calculating intersection between partition and division line. }
  { Use real types (Extended) to avoid aritmetic overflow.                 }
  Denominator: Extended;
  Numerator  : Extended;
  Fraction   : Extended;
  XPar,YPar  : Extended; {starting point of partition line}
  XDiv,YDiv  : Extended; {starting point of division line}
  dXPar,dYPar: Extended; {delta (slope) of partition line}
  dXDiv,dYDiv: Extended; {delta (slope) of division line}

begin
  Result := -1;
  if Points < 2 then Exit;

  { Compute the delta X and delta Y for the partition line. }
  { As partition line for the calculations we use a line    }
  { starting in point P and extending 10 values rightwards. }
  XPar := P.X;
  dXPar := XPar+10; {extend orignal point to make a line}
  dXPar := dXPar-XPar;
  YPar := P.Y;
  dYPar := YPar;
  dYPar := dYPar-YPar;

  Count := 0;
  Index := 0;
  while (Index < Points) and (Result <> 0) do begin
    { Compute the delta X and delta Y for the current division (edge)
line. }
    XDiv := Polygon[Index].X;
    dXDiv := Polygon[(Index+1) mod Points].X;
    dXDiv := dXDiv-XDiv;
    YDiv := Polygon[Index].Y;
    dYDiv := Polygon[(Index+1) mod Points].Y;
    dYDiv := dYDiv-YDiv;

    { Compute denominator and numerator to find the relation }
    { of the division line to the partition line.}
    Denominator := (dYPar*dXDiv)-(dXPar*dYDiv);
    Numerator := ((XPar-XDiv)*dYPar)+((YDiv-YPar)*dXPar);

    if not (Denominator = 0.0) then begin
    { The lines are not parallel. }
      { Calculat the fraction of the division line for the intersection
point. }
      Fraction := Numerator/Denominator;

      if (Fraction >= 0.0) and (Fraction < 1.0) then begin
      { The intersection point is between the starting and ending }
      { point of the division line, or on it's starting point. }
      { For each line we include the starting point but not the }
      { ending point as part of the line. The ending point of a }
      { given line will be the starting point of the next line. }

        { Recompute denominator, numerator, and fraction to find }
        { the relation of the partition line to the division line.}
        Denominator := (dYDiv*dXPar)-(dXDiv*dYPar);
        Numerator := ((XDiv-XPar)*dYDiv)+((YPar-YDiv)*dXDiv);
        Fraction := Numerator/Denominator;

        if Fraction = 0.0 then Result := 0 { point P is on the edge }
        else if Fraction > 0.0 then Inc(Count);
      end;
    end;

    Inc(Index);
  end;

  if Result <> 0 then Result := Count and 1;
end;

{-----------------------------------------------------------------------}
type
  TLinesInLine = set of (
    liDivToRight,  { division line is to the right of partition line }
    liDivToLeft,   { division line is to the left of partition line  }
    liAtDivStart,  { intersection at start of division line          }
    liAtDivEnd,    { intersection at end of division line            }
    liOffDivStart, { intersection beyond start of division line      }
    liOffDivEnd,   { intersection beyond end of division line        }
    liIntInRange,  { intersection point is in range of TPoint type   }
    liParallel     { division line is parallel to partition line     }
  );

function LineInLine(P1,P2,D1,D2: TPoint; var P: TPoint): TLinesInLine;
  { Returnes the point where the partition line (from P1 to P2)
rsects     }
  { the division line (from D1 to D2) or the point where the lines'
finite   }
  { extentions will intersect. The intersection point gives meaning only
when  }
  { the lines are not
          }
  { (see type TLinesIntersection for the returned
       }
var
  Denominator: Extended;
  Numerator  : Extended;
  Fraction   : Extended;
  XInt,YInt  : Extended;
  XPar,YPar  : Extended;
  XDiv,YDiv  : Extended;
  dXPar,dYPar: Extended;
  dXDiv,dYDiv: Extended;
begin
  { Compute the delta X and delta Y for both lines, and move all the  }
  { parameters into real types to avoid arithmetic overflow error. }
  XPar := P1.X;
  YPar := P1.Y;
  XDiv := D1.X;
  YDiv := D1.Y;
  dXPar := P2.X; dXPar := dXPar-P1.X;
  dYPar := P2.Y; dYPar := dYPar-P1.Y;
  dXDiv := D2.X; dXDiv := dXDiv-D1.X;
  dYDiv := D2.Y; dYDiv := dYDiv-D1.Y;
  Result := [];
  { First we compute the denominator and numerator to find  }
  { the relation of the division line to the partition line.}
  Denominator := (dYPar*dXDiv)-(dXPar*dYDiv);
  Numerator := ((XPar-XDiv)*dYPar)+((YDiv-YPar)*dXPar);
  { Find the intersection in relation to the division line. }
  if Denominator = 0.0 then begin { The two lines are parallel }
    Result := Result + [liParallel];
    { Find which side of the partition line the division line is on. }
    case CoordinateMode of
      cmWindows:
        if Numerator > 0.0 then Result := Result + [liDivToRight]
        else if Numerator < 0.0 then Result := Result + [liDivToLeft];
      cmCartesian:
        if Numerator > 0.0 then Result := Result + [liDivToLeft]
        else if Numerator < 0.0 then Result := Result + [liDivToRight];
    end;
  end
  else begin { The lines are NOT parallel }
    Fraction := Numerator/Denominator;
    { Compute the intersection point }
    XInt := Round(XDiv+dXDiv*Fraction);
    YInt := Round(YDiv+dYDiv*Fraction);
    if ((XInt >= -2147483648.0) and (XInt <= 2147483647.0))
    and ((YInt >= -2147483648.0) and (YInt <= 2147483647.0)) then begin
      { Intersection point is within range of the TPoint type }
      Result := Result + [liIntInRange];
      P.X := Round(XInt);
      P.Y := Round(YInt);
      { Check if the intersection point actually falls on the    }
      { start or end of the division line, and reassign fraction }
      { and numerator to reflect this.                           }
      if (P.X = D1.X) and (P.Y = D1.Y) then begin
        if Fraction > 0.0 then Numerator := -Numerator;
        Fraction := 0.0;
      end
      else if (P.X = D2.X) and (P.Y = D2.Y) then Fraction := 1.0;
    end;

    { Find which side of the partition line the division line is on. }
    case CoordinateMode of
      cmWindows:
        if Numerator > 0.0 then Result := Result + [liDivToRight]
        else if Numerator < 0.0 then Result := Result + [liDivToLeft]
        else begin
          { If numerator=0 we use denominator to find the side. }
          if Denominator < 0.0 then Result := Result + [liDivToRight]
          else if Denominator > 0.0 then Result := Result + [liDivToLeft];
        end;
      cmCartesian:
        if Numerator < 0.0 then Result := Result + [liDivToRight]
        else if Numerator > 0.0 then Result := Result + [liDivToLeft]
        else begin
          { If numerator=0 we use denominator to find the side. }
          if Denominator > 0.0 then Result := Result + [liDivToRight]
          else if Denominator < 0.0 then Result := Result + [liDivToLeft];
        end;
    end;

    if Fraction < 0.0 then Result := Result + [liOffDivStart]
      { Intersection is beyond the start of the division line. }
    else if Fraction = 0.0 then Result := Result + [liAtDivStart]
      { Intersection at the start of the division line. }
    else if Fraction > 1.0 then Result := Result + [liOffDivEnd]
      { Intersection is beyond the end of the division line. }
    else if Fraction = 1.0 then Result := Result + [liAtDivEnd]
      { Intersection at the end of the division line. }
    else Result := Result + [liDivToRight,liDivToLeft];
      { Intersection between the start and the end of the division line. }
  end;
end;