{
This unit contains a fpc implementation of the C-code published on
https://web.archive.org/web/20130126163405/http://geomalgorithms.com/a03-_inclusion.html
All original rights are included.
Translated from C to FPC by Ekkehard Domning edo@domis.de
2023-06-08 Changed parameters and bugfix

// Copyright 2000 softSurfer, 2012 Dan Sunday
// This code may be freely used and modified for any purpose
// providing that this copyright notice is included with it.
// SoftSurfer makes no warranty for this code, and cannot be held
// liable for any real or imagined damage resulting from its use.
// Users of this code must verify correctness for their application.
}

unit uWnPnPoly;

{$mode ObjFPC}{$H+}

interface

{
 CrossingNumberPointInPolygon the crossing point test of a polygon.
       Input:
         AX, AY : A test point,
                  V[] = vertex points of a polygon V[n+1] with V[n]=V[0]
         APolygon : vertex points of a polygon V[n+1] with V[n]=V[0]
                  The Polygon is defined by pairs of Double values.
                  The method of passing an untyped reference allows the usage of
                  several datatype like packed records etc.
                  A dynamic Array is passed by referencing the first item
                  The last point of the Polygon must be identical to the first point,
                  to close the polygon.
                  Self-Overlapping of the polygon is *not* allowed
         APolygonPointCount: The Number of points in the polygon =
                  twice the number of doubles in APolygon.
         AFirstValueIsY: If True, the first double value in APolygon is the Y-Value of the coordinate.
                  If False, the first value is the X-Value.
         AYAxisUpPlus : If True the values of Y increases when going in direction top.
                  If false they decreases.
       Return :   True if inside, false if outside
}
function CrossingNumberPointInPolygon(const AX, AY: Double;
                               const APolygon;
                               const APolygonPointCount: Integer;
                               const AFirstValueIsY: Boolean = False;
                               const AYAxisUpPlus : Boolean = False) : Boolean;

{
 WindingPointInPolygon: winding number test for a point in a polygon
      Input:
        AX, AY : A test point,
        APolygon : vertex points of a polygon V[n+1] with V[n]=V[0]
                The Polygon is defined by pairs of Double values.
                The method of passing an untyped reference allows the usage of
                several datatype like packed records etc.
                A dynamic Array is passed by referencing the first item
                The last point of the Polygon must be identical to the first point,
                to close the polygon.
                Self-Overlapping of the polygon is allowed
        APolygonPointCount: The Number of points in the polygon =
                twice the number of doubles in APolygon.
        AFirstValueIsY: If True, the first double value in APolygon is the Y-Value of the coordinate.
                If False, the first value is the X-Value.
        AYAxisUpPlus : If True the values of Y increases when going in direction top.
                 If false they decreases.
      Result:  wn = the winding number (=0 only when P is outside)
}
function WindingPointInPolygon(const AX, AY: Double;
                               const APolygon;
                               const APolygonPointCount: Integer;
                               const AFirstValueIsY: Boolean = False;
                               const AYAxisUpPlus : Boolean = False) : Integer;
{  Testing Example
procedure TForm1.Button1Click(Sender: TObject);
var
  poly : array of Double;
  rc : Integer;
  b : Boolean;
  ax,ay : Double;
  s : String;
begin
  SetLength(poly,5*2);
  //           X               Y
  poly[0] :=  0.0; poly[1] :=   0.0;
  poly[2] := 10.0; poly[3] :=   0.0;
  poly[4] := 10.0; poly[5] :=  -10.0;
  poly[6] :=  0.0; poly[7] :=  -10.0;
  poly[8] :=  0.0; poly[9] :=   0.0;
  ax := 5.0;
  ay := -5.0;
  rc := WindingPointInPolygon(ax,ay,
                              poly[0],
                              Length(poly) div 2,
                              False, // First value is not Y (its X)
                              True); // The Y value increases going up
  b := CrossingNumberPointInPolygon(ax,ay,
                              poly[0],
                              Length(poly) div 2,
                              False, // First value is not Y (its X)
                              True); // The Y value increases going up
  if b then
    s := 'Crossing result: Inside'
  else
    s := 'Crossing result: Outside';
  Label1.Caption := Format('%s , WindingResult: %d',[s,rc]);
end;
}



implementation
type
  TStaticDoubleArray = array[0..(MaxInt div SizeOf(Double))-1] of Double;

(*
// Copyright 2000 softSurfer, 2012 Dan Sunday
// This code may be freely used and modified for any purpose
// providing that this copyright notice is included with it.
// SoftSurfer makes no warranty for this code, and cannot be held
// liable for any real or imagined damage resulting from its use.
// Users of this code must verify correctness for their application.


// a Point is defined by its coordinates {int x, y;}
//===================================================================


// isLeft(): tests if a point is Left|On|Right of an infinite line.
//    Input:  three points P0, P1, and P2
//    Return: >0 for P2 left of the line through P0 and P1
//            =0 for P2  on the line
//            <0 for P2  right of the line
//    See: Algorithm 1 "Area of Triangles and Polygons"
inline int
isLeft( Point P0, Point P1, Point P2 )
{
    return ( (P1.x - P0.x) * (P2.y - P0.y)
            - (P2.x -  P0.x) * (P1.y - P0.y) );
}
//===================================================================


// cn_PnPoly(): crossing number test for a point in a polygon
//      Input:   P = a point,
//               V[] = vertex points of a polygon V[n+1] with V[n]=V[0]
//      Return:  0 = outside, 1 = inside
// This code is patterned after [Franklin, 2000]
int
cn_PnPoly( Point P, Point* V, int n )
{
    int    cn = 0;    // the  crossing number counter

    // loop through all edges of the polygon
    for (int i=0; i<n; i++) {    // edge from V[i]  to V[i+1]
       if (((V[i].y <= P.y) && (V[i+1].y > P.y))     // an upward crossing
        || ((V[i].y > P.y) && (V[i+1].y <=  P.y))) { // a downward crossing
            // compute  the actual edge-ray intersect x-coordinate
            float vt = (float)(P.y  - V[i].y) / (V[i+1].y - V[i].y);
            if (P.x <  V[i].x + vt * (V[i+1].x - V[i].x)) // P.x < intersect
                 ++cn;   // a valid crossing of y=P.y right of P.x
        }
    }
    return (cn&1);    // 0 if even (out), and 1 if  odd (in)

}
//===================================================================


// wn_PnPoly(): winding number test for a point in a polygon
//      Input:   P = a point,
//               V[] = vertex points of a polygon V[n+1] with V[n]=V[0]
//      Return:  wn = the winding number (=0 only when P is outside)
int
wn_PnPoly( Point P, Point* V, int n )
{
    int    wn = 0;    // the  winding number counter

    // loop through all edges of the polygon
    for (int i=0; i<n; i++) {   // edge from V[i] to  V[i+1]
        if (V[i].y <= P.y) {          // start y <= P.y
            if (V[i+1].y  > P.y)      // an upward crossing
                 if (isLeft( V[i], V[i+1], P) > 0)  // P left of  edge
                     ++wn;            // have  a valid up intersect
        }
        else {                        // start y > P.y (no test needed)
            if (V[i+1].y  <= P.y)     // a downward crossing
                 if (isLeft( V[i], V[i+1], P) < 0)  // P right of  edge
                     --wn;            // have  a valid down intersect
        }
    }
    return wn;
}
//===================================================================
*)


// cn_PnPoly(): crossing number test for a point in a polygon
//      Input:   P = a point,
//               V[] = vertex points of a polygon V[n+1] with V[n]=V[0]
//      Return:  0 = outside, 1 = inside
// This code is patterned after [Franklin, 2000]
function CrossingNumberPointInPolygon(const AX, AY: Double; const APolygon;
  const APolygonPointCount: Integer; const AFirstValueIsY: Boolean;
  const AYAxisUpPlus: Boolean): Boolean;
var
  pts : TStaticDoubleArray absolute APolygon;
  cn : Integer;
  i : Integer;
  yofs, xofs : Integer;
  x0, x1, y0, y1 : Double;
  vt : Double;
  ysign : Double;
  lY : Double;
begin
  Result := False;
  if APolygonPointCount < 4 then Exit; // At least a closed triangle is needed to have any point inside
  if AYAxisUpPlus then
    ysign := -1.0
  else
    ysign := 1.0;
  if AFirstValueIsY then
  begin
    yofs := 0;
    xofs := 1;
  end
  else
  begin
    yofs := 1;
    xofs := 0;
  end;
  lY := AY * ysign; // The local copy of AY, sign adjusted
  x0 := pts[0+xofs];
  y0 := pts[0+yofs] * ysign;
  cn := 0;    // the  crossing number counter
  // loop through all edges of the polygon
  for i := 1 to APolygonPointCount-1 do
  begin
    // edge from V[i-1]  to V[i]
    x1 := pts[(i shl 1)+xofs];
    y1 := pts[(i shl 1)+yofs] * ysign;
    if (((y0 <= lY) and (y1 > lY)) or     // an upward crossing
        ((y0 > lY) and (y1 <=  lY))) then // a downward crossing
    begin
      // compute  the actual edge-ray intersect x-coordinate
      vt := (lY  - y0) / (y1 - y0);
      if AX < (x0 + vt * (x1 - x0)) then // P.x < intersect
        Inc(cn);  // a valid crossing of y=P.y right of P.x
    end;
    x0 := x1;
    y0 := y1;
  end;
  Result := Odd(cn); // False if even (out), and True if odd (in)
end;

// wn_PnPoly(): winding number test for a point in a polygon
//      Input:   P = a point,
//               V[] = vertex points of a polygon V[n+1] with V[n]=V[0]
//      Return:  wn = the winding number (=0 only when P is outside)
function WindingPointInPolygon(const AX, AY: Double; const APolygon;
  const APolygonPointCount: Integer; const AFirstValueIsY: Boolean;
  const AYAxisUpPlus: Boolean): Integer;

  // isLeft(): tests if a point is Left|On|Right of an infinite line.
  //    Input:  three points P0, P1, and P2
  //    Return: >0 for P2 left of the line through P0 and P1
  //            =0 for P2  on the line
  //            <0 for P2  right of the line
  //    See: Algorithm 1 "Area of Triangles and Polygons"

  function IsLeft(const X0, Y0, X1, Y1, X2, Y2 : Double) : Double;inline;
  begin
    Result := ((X1-X0) * (Y2-Y0)) - ((X2-X0) * (Y1-Y0));
  end;

var
  wn : Integer;
  pts : TStaticDoubleArray absolute APolygon;
  i : Integer;
  yofs, xofs : Integer;
  x0, x1, y0, y1 : Double;
  ysign : Double;
  lY : Double;
begin
  Result := 0;
  if APolygonPointCount < 4 then Exit; // At least a closed triangle is needed to have any point inside
  if AYAxisUpPlus then
    ysign := -1.0
  else
    ysign := 1.0;
  if AFirstValueIsY then
  begin
    yofs := 0;
    xofs := 1;
  end
  else
  begin
    yofs := 1;
    xofs := 0;
  end;
  lY := AY * ysign; //Local copy of AY, sign adjusted
  x0 := pts[0+xofs];
  y0 := pts[0+yofs] * ysign;
  wn := 0; // The winding number

  // loop through all edges of the polygon
  for i := 1 to APolygonPointCount-1 do
  begin
    // edge from V[i-1]  to V[i]
    x1 := pts[(i shl 1)+xofs];
    y1 := pts[(i shl 1)+yofs] * ysign;
    if (y0 <= lY) then      // start y <= P.y
    begin
      if (y1 > lY) and    // an upward crossing
         (IsLeft(x0,y0,x1,y1,AX,lY) > 0) then  // P left of edge
         Inc(wn);          // have  a valid up intersect
    end
    else                  // start y > P.y (no test needed)
    begin
      if (y1 <= lY) and   // a downward crossing
         (IsLeft(x0,y0,x1,y1,AX,lY) < 0) then  // P right of edge
        Dec(wn);          // have a valid down intersect
    end;
    x0 := x1;
    y0 := y1;
  end;
  Result := wn;
end;

end.