{
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

// 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.
       Return:  True if inside, false if outside
}
function CrossingNumberPointInPolygon(const AX, AY: Double;
                               const APolygon;
                               const APolygonPointCount: Integer;
                               const AFirstValueIsY: Boolean) : 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.
      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) : Integer;

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;
}
//===================================================================
*)

// 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 AX, AY, BX, BY, CX, CY : Double) : Double;
begin
  Result := ((BX-AX) * (CY-AY)) - ((CX-AX) * (BY-AY));
end;

// 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) : Boolean;
var
  pts : TStaticDoubleArray absolute APolygon;
  cn : Integer;
  i : Integer;
  yofs, xofs : Integer;
  x0, x1, y0, y1 : Double;
  vt : Double;
begin
  cn := 0;    // the  crossing number counter
  Result := False;
  if APolygonPointCount < 4 then Exit; // At least a closed triangle is needed to have any point inside
  if AFirstValueIsY then
  begin
    yofs := 0;
    xofs := 1;
  end
  else
  begin
    yofs := 1;
    xofs := 0;
  end;
  x0 := pts[0+xofs];
  y0 := pts[0+yofs];
  // 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];
    if (((y0 <= AY) and (y1 > AY)) or     // an upward crossing
        ((y0 > AY) and (y1 <=  AY))) then // a downward crossing
    begin
      // compute  the actual edge-ray intersect x-coordinate
      vt := (AY  - 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) : Integer;
var
  wn : Integer;
  pts : TStaticDoubleArray absolute APolygon;
  i : Integer;
  yofs, xofs : Integer;
  x0, x1, y0, y1 : Double;
begin
  Result := 0;
  if APolygonPointCount < 4 then Exit; // At least a closed triangle is needed to have any point inside
  wn := 0;
  if AFirstValueIsY then
  begin
    yofs := 0;
    xofs := 1;
  end
  else
  begin
    yofs := 1;
    xofs := 0;
  end;
  x0 := pts[0+xofs];
  y0 := pts[0+yofs];
  // 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];
    if (y0 <= AY) then      // start y <= P.y
    begin
      if (y1 > AY) then     // an upward crossing
      begin
        if (y0 <= AY) and   // start y <= P.y
           (IsLeft(x0,y0,x1,y1,AX,AY) > 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 <= AY) and   // a downward crossing
           (IsLeft(x0,y0,x1,y1,AX,AY) < 0) then  // P right of edge
          Dec(wn);          // have a valid down intersect
      end;
    end;
    x0 := x1;
    y0 := y1;
  end;
  Result := wn;
end;

end.

