unit BeRoTriangulation;
(*****************************************************************************
**                                                                           *
** The TriangulateDelaunayClipping function needs and uses the Clipper       *
** library from phttp://www.angusj.com/delphi/clipper.php for the polygon    *
** clipping work                                                             *
**                                                                           *
******************************************************************************
**        Copyright 2014 Benjamin Rosseaux. All rights reserved.             *
******************************************************************************
**
** Boost Software License - Version 1.0 - August 17th, 2003
**
** Permission is hereby granted, free of charge, to any person or organization
** obtaining a copy of the software and accompanying documentation covered by
** this license (the "Software") to use, reproduce, display, distribute,
** execute, and transmit the Software, and to prepare derivative works of the
** Software, and to permit third-parties to whom the Software is furnished to
** do so, all subject to the following:
**
** The copyright notices in the Software and this entire statement, including
** the above license grant, this restriction and the following disclaimer,
** must be included in all copies of the Software, in whole or in part, and
** all derivative works of the Software, unless such copies or derivative
** works are solely in the form of machine-executable object code generated by
** a source language processor.
**
** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
** IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
** FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
** SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
** FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
** ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
** DEALINGS IN THE SOFTWARE.
**
******************************************************************************
**                                                                           *
** The views and conclusions contained in the software and documentation are *
** those of the authors and should not be interpreted as representing        *
** official policies, either expressed or implied, of Benjamin Rosseaux.     *
**                                                                           *
******************************************************************************)
{$ifdef fpc}
 {$mode delphi}
 {$warnings off}
 {$hints off}
 {$ifdef cpui386}
  {$define cpu386}
  {$asmmode intel}
 {$endif}
 {$ifdef fpc_little_endian}
  {$define little_endian}
 {$else}
  {$ifdef fpc_big_endian}
   {$define big_endian}
  {$endif}
 {$endif}
 {$ifdef fpc_has_internal_sar}
  {$define HasSAR}
 {$endif}
{$else}
 {$define little_endian}
 {$ifndef cpu64}
  {$define cpu64}
 {$endif}
 {$optimization on}
 {$undef HasSAR}
 {$define UseDIV}
{$endif}
{$overflowchecks off}
{$rangechecks off}

{$define FASTER} // For faster results, sacrificing triangulation quality for a bit better speed performance

{$undef FASTONLY} // For only fast seidel-algorithm-based variant

{$undef STANDALONE} // For only to have the TriangulateBeRo variant without any unit dependencies except the RTL itself, of course

interface

type PBeRoTriangulationPoint=^TBeRoTriangulationPoint;
     TBeRoTriangulationPoint=record
      x,y:double;
     end;

     TBeRoTriangulationPolygon=array of TBeRoTriangulationPoint;

     TBeRoTriangulationPolygons=array of TBeRoTriangulationPolygon;

     PBeRoTriangulationTriangle=^TBeRoTriangulationTriangle;
     TBeRoTriangulationTriangle=array[0..2] of TBeRoTriangulationPoint;

     TBeRoTriangulationTriangles=array of TBeRoTriangulationTriangle;

     TBeRoTriangulationPolygonFillRule=(btpfrEVENODD,btpfrNONZERO,btpfrPOSITIVE,btpfrNEGATIVE,btpfrABSGEQTWO);

// The flexible semi-robust variant by my own scanline-rasterization-like sweeping algoritm. Algorithm steps:
// 1. Add edges from input polygons
// 2. Split intersecting edges
// 3. Split edges at Y coordinates of all polygon points including these of intersection edge intersection points
// 4. Construct y-monotone quads with a trapezoid-style shape with scanline-rasterization-like sweeping from
//    top to bottom and from left to right with the wished polygon fill rule
// 5. Optimize and merge the output from step 4. as far as possible
// This variant is still experimental, so do use this function variant on your own risk! :-)
// It works already pretty good, but it is still in progress.
function TriangulateBeRo(const InputPolygons:TBeRoTriangulationPolygons;var OutputTriangles:TBeRoTriangulationTriangles;InputPolygonFillRule:TBeRoTriangulationPolygonFillRule=btpfrEVENODD;Quality:longint=$7fffffff):boolean;

{$ifndef STANDALONE}
{$ifndef FASTONLY}
// Slow bruteforce-style very robust variant with delaunay algorithm, clipping and earclipping for non-realtime usages
// A robust 2D polygon triangulator by combining non-constrained delaunay triangulation as first pass with polygon
// clipping as middle pass and ear clipping as final pass, so it supports even self-intersecting polygon with holes
// as input. It needs and uses the Clipper library from http://www.angusj.com/delphi/clipper.php for the polygon
// clipping work.
// It doesn't support the btpfrABSGEQTWO polygon fill rule, and the output depends on the correctness/accuracy of the
// 3rd-party clipper library
function TriangulateDelaunayClipping(const InputPolygons:TBeRoTriangulationPolygons;var OutputTriangles:TBeRoTriangulationTriangles;InputPolygonFillRule:TBeRoTriangulationPolygonFillRule=btpfrEVENODD):boolean;
{$endif}

// Faster non-robust variant with seidel algorithm for realtime usages
// It supports only the even-odd polygon fill rule in this current implementation, and it can handle only simple polygons
// (with and without holes)
function TriangulateSeidel(const InputPolygons:TBeRoTriangulationPolygons;var OutputTriangles:TBeRoTriangulationTriangles):boolean;
{$endif}

// Faster non-robust variant with ear clipping algorithm for realtime usages
// It can handle only non-self-intersecting simple polygons without holes
function TriangulateEarClipping(const InputPolygon:TBeRoTriangulationPolygon;var OutputTriangles:TBeRoTriangulationTriangles):boolean;

implementation

{$ifndef STANDALONE}
uses Math{$ifndef FASTONLY},Clipper{$endif};
{$endif}

function TriangulateBeRo(const InputPolygons:TBeRoTriangulationPolygons;var OutputTriangles:TBeRoTriangulationTriangles;InputPolygonFillRule:TBeRoTriangulationPolygonFillRule=btpfrEVENODD;Quality:longint=$7fffffff):boolean;
const DoubleResolution=1e-15*1000.0;
      HashBits=8;
      HashSize=1 shl HashBits;
      HashMask=HashSize-1;
      AABBMULTIPLIER=2.0;
      AABBEPSILON=1e-4;
      daabbtNULLNODE=-1;
      TOPLEFT=0;
      TOPRIGHT=1;
      BOTTOMRIGHT=2;
      BOTTOMLEFT=3;
type PGarbageCollectedHeader=^TGarbageCollectedHeader;
     PGarbageCollector=^TGarbageCollector;
     TGarbageCollector=record
      First:PGarbageCollectedHeader;
      Last:PGarbageCollectedHeader;
     end;
     TGarbageCollectedHeader=record
      Previous:PGarbageCollectedHeader;
      Next:PGarbageCollectedHeader;
     end;
     PPPoint=^PPoint;
     PPoint=^TPoint;
     PEdge=^TEdge;
     PEdges=^TEdges;
     PVerticalAATreeNode=^TVerticalAATreeNode;
     PQuad=^TQuad;
     PQuads=^TQuads;
     PTriangle=^TTriangle;
     PTriangles=^TTriangles;
     PPointHashTable=^TPointHashTable;
     PPolygonVertex=^TPolygonVertex;
     PPolygon=^TPolygon;
     PPolygons=^TPolygons;
     PVector=^TVector;
     TVector=record
      x:double;
      y:double;
     end;
     PAABB=^TAABB;
     TAABB=record
      Min:TVector;
      Max:TVector;
     end;
     PDynamicAABBTreeNode=^TDynamicAABBTreeNode;
     TDynamicAABBTreeNode=record
      AABB:TAABB;
      UserData:pointer;
      Children:array[0..1] of longint;
      Height:longint;
      case boolean of
       false:(
        Parent:longint;
       );
       true:(
        Next:longint;
       );
     end;
     PDynamicAABBTreeNodes=^TDynamicAABBTreeNodes;
     TDynamicAABBTreeNodes=array[0..0] of TDynamicAABBTreeNode;
     PDynamicAABBTreeLongintArray=^TDynamicAABBTreeLongintArray;
     TDynamicAABBTreeLongintArray=array[0..65535] of longint;
     PDynamicAABBTree=^TDynamicAABBTree;
     TDynamicAABBTree=record
      Root:longint;
      Nodes:PDynamicAABBTreeNodes;
      NodeCount:longint;
      NodeCapacity:longint;
      FreeList:longint;
      Path:longword;
      InsertionCount:longint;
      Stack:PDynamicAABBTreeLongintArray;
      StackCapacity:longint;
     end;
     TPoint=record
      GarbageCollectedHeader:TGarbageCollectedHeader;
      HashNext:PPoint;
      Hash:longword;
      Previous:PPoint;
      Next:PPoint;
      x:double;
      y:double;
     end;
     TEdge=record
      GarbageCollectedHeader:TGarbageCollectedHeader;
      Previous:PEdge;
      Next:PEdge;
      p:PPoint;
      q:PPoint;
      MinX:double;
      MaxX:double;
      MinY:double;
      MaxY:double;
     end;
     TEdges=record
      GarbageCollectedHeader:TGarbageCollectedHeader;
      First:PEdge;
      Last:PEdge;
     end;
     TVerticalAATreeNode=record
      GarbageCollectedHeader:TGarbageCollectedHeader;
      Parent:PVerticalAATreeNode;
      Left:PVerticalAATreeNode;
      Right:PVerticalAATreeNode;
      Level:longint;
      YCoordinate:double;
      First:PEdge;
      Last:PEdge;
     end;
     TQuad=record
      GarbageCollectedHeader:TGarbageCollectedHeader;
      Previous:PQuad;
      Next:PQuad;
      Points:array[0..3] of PPoint;
      Proxy:longint;
      AABB:TAABB;
      Polygon:PPolygon;
     end;
     TQuads=record
      GarbageCollectedHeader:TGarbageCollectedHeader;
      First:PQuad;
      Last:PQuad;
     end;
     TTriangle=record
      GarbageCollectedHeader:TGarbageCollectedHeader;
      Previous:PTriangle;
      Next:PTriangle;
      Points:array[0..2] of PPoint;
      Proxy:longint;
      AABB:TAABB;
      Polygon:PPolygon;
     end;
     TTriangles=record
      GarbageCollectedHeader:TGarbageCollectedHeader;
      First:PTriangle;
      Last:PTriangle;
     end;
     TPointHashTable=array[0..HashSize-1] of PPoint;
     TPolygonVertex=record
      GarbageCollectedHeader:TGarbageCollectedHeader;
      x:double;
      y:double;
      Previous:PPolygonVertex;
      Next:PPolygonVertex;
      Corresponding:PPolygonVertex;
      Distance:double;
      IsEntry:longbool;
      IsIntersection:longbool;
      Visited:longbool;
     end;
     TPolygon=record
      GarbageCollectedHeader:TGarbageCollectedHeader;
      Previous:PPolygon;
      Next:PPolygon;
      First:PPolygonVertex;
      Count:longint;
     end;
     TPolygons=record
      GarbageCollectedHeader:TGarbageCollectedHeader;
      First:PPolygon;
      Last:PPolygon;
      Count:longint;
     end;
var GarbageCollector:TGarbageCollector;
    Edges:TEdges;
    Quads:TQuads;
    Triangles:TTriangles;
    VerticalAATreeNodeRoot:TVerticalAATreeNode;
    CountOutputTriangles:longint;
    PointHashTable:PPointHashTable;
    QuadDynamicAABBTree:PDynamicAABBTree;
    TriangleDynamicAABBTree:PDynamicAABBTree;
    DoOptimizeTriangles:longbool;
 function HashXY(const x,y:double):longword;
 begin
  result:=(round(x)*73856093) xor (round(y)*19349663);
 end;
 procedure OutputTriangle(const ax,ay,bx,by,cx,cy:double);
 var Triangle:PBeRoTriangulationTriangle;
 begin
  if length(OutputTriangles)<=CountOutputTriangles then begin
   SetLength(OutputTriangles,(CountOutputTriangles+1)*2);
  end;
  Triangle:=@OutputTriangles[CountOutputTriangles];
  inc(CountOutputTriangles);
  if (((bx-ax)*(by+ay))+((cx-bx)*(cy+by))+((ax-cx)*(ay+cy)))<0.0 then begin
   Triangle^[0].x:=ax;
   Triangle^[0].y:=ay;
   Triangle^[1].x:=bx;
   Triangle^[1].y:=by;
   Triangle^[2].x:=cx;
   Triangle^[2].y:=cy;
  end else begin
   Triangle^[0].x:=cx;
   Triangle^[0].y:=cy;
   Triangle^[1].x:=bx;
   Triangle^[1].y:=by;
   Triangle^[2].x:=ax;
   Triangle^[2].y:=ay;
  end;
 end;
 function GarbageCollectorAllocate(Size:longint):pointer;
 begin
  GetMem(result,Size);
  FillChar(result^,Size,0);
  if assigned(GarbageCollector.Last) then begin
   GarbageCollector.Last^.Next:=result;
   PGarbageCollectedHeader(result)^.Previous:=GarbageCollector.Last;
  end else begin
   GarbageCollector.First:=result;
   PGarbageCollectedHeader(result)^.Previous:=nil;
  end;
  GarbageCollector.Last:=result;
  PGarbageCollectedHeader(result)^.Next:=nil;
 end;
 procedure GarbageCollectorFree(Data:pointer);
 begin
  if assigned(PGarbageCollectedHeader(Data)^.Previous) then begin
   PGarbageCollectedHeader(Data)^.Previous^.Next:=PGarbageCollectedHeader(Data)^.Next;
  end else if GarbageCollector.First=Data then begin
   GarbageCollector.First:=PGarbageCollectedHeader(Data)^.Next;
  end;
  if assigned(PGarbageCollectedHeader(Data)^.Next) then begin
   PGarbageCollectedHeader(Data)^.Next^.Previous:=PGarbageCollectedHeader(Data)^.Previous;
  end else if GarbageCollector.Last=Data then begin
   GarbageCollector.Last:=PGarbageCollectedHeader(Data)^.Previous;
  end;
  PGarbageCollectedHeader(Data)^.Previous:=nil;
  PGarbageCollectedHeader(Data)^.Next:=nil;
  FreeMem(Data);
 end;
 procedure GarbageCollectorCleanUp;
 var Current,Next:pointer;
 begin
  Current:=GarbageCollector.First;
  GarbageCollector.First:=nil;
  GarbageCollector.Last:=nil;
  while assigned(Current) do begin
   Next:=PGarbageCollectedHeader(Current)^.Next;
   FreeMem(Current);
   Current:=Next;
  end;
 end;
 function Min(const a,b:double):double; overload;
 begin
  if a<b then begin
   result:=a;
  end else begin
   result:=b;
  end;
 end;
 function Min(const a,b:longint):longint; overload;
 begin
  if a<b then begin
   result:=a;
  end else begin
   result:=b;
  end;
 end;
 function Max(const a,b:double):double; overload;
 begin
  if a<b then begin
   result:=b;
  end else begin
   result:=a;
  end;
 end;
 function Max(const a,b:longint):longint; overload;
 begin
  if a<b then begin
   result:=b;
  end else begin
   result:=a;
  end;
 end;
 function Sign(const a:double):longint; overload;
 begin
  if a<0.0 then begin
   result:=-1;
  end else if a>0.0 then begin
   result:=1;
  end else begin
   result:=0;
  end;
 end;
 function Sign(const a:longint):longint; overload;
 begin
  if a<0.0 then begin
   result:=-1;
  end else if a>0.0 then begin
   result:=1;
  end else begin
   result:=0;
  end;
 end;
 function IsZero(const a:double):boolean;
 begin
  Result:=abs(a)<=DoubleResolution;
 end;
 function SameValue(const a,b:double):boolean;
 var Epsilon:double;
 begin
  Epsilon:=Max(Min(abs(a),abs(b))*DoubleResolution,DoubleResolution);
  if a>b then begin
   result:=(a-b)<=Epsilon;
  end else begin
   result:=(b-a)<=Epsilon;
  end;
 end;
 function ExactSameValue(const a,b:double):boolean;
 begin
  result:=(a=b) or (int64(pointer(@a)^)=int64(pointer(@b)^));
 end;
 function LerpValue(const a,b,t:double):double;
 begin
  if t<=0.0 then begin
   result:=a;
  end else if t>=1.0 then begin
   result:=b;
  end else begin
   result:=(a*(1.0-t))+(b*t);
  end;
 end;
 function AABBCombine(const AABB,WithAABB:TAABB):TAABB;
 begin
  result.Min.x:=Min(AABB.Min.x,WithAABB.Min.x);
  result.Min.y:=Min(AABB.Min.y,WithAABB.Min.y);
  result.Max.x:=Max(AABB.Max.x,WithAABB.Max.x);
  result.Max.y:=Max(AABB.Max.y,WithAABB.Max.y);
 end;
 function AABBCost(const AABB:TAABB):double;
 begin
//result:=(AABB.Max.x-AABB.Min.x)+(AABB.Max.y-AABB.Min.y); // Manhattan distance
  result:=(AABB.Max.x-AABB.Min.x)*(AABB.Max.y-AABB.Min.y); // Volume
 end;
 function AABBContains(const InAABB,AABB:TAABB):boolean;
 begin
  result:=(InAABB.Min.x<=AABB.Min.x) and (InAABB.Min.y<=AABB.Min.y) and
          (InAABB.Max.x>=AABB.Min.x) and (InAABB.Max.y>=AABB.Min.y) and
          (InAABB.Min.x<=AABB.Max.x) and (InAABB.Min.y<=AABB.Max.y) and
          (InAABB.Max.x>=AABB.Max.x) and (InAABB.Max.y>=AABB.Max.y);
 end;
 function AABBIntersect(const AABB,WithAABB:TAABB):boolean;
 begin
  result:=((AABB.Max.x>=WithAABB.Min.x) and (AABB.Min.x<=WithAABB.Max.x)) and
          ((AABB.Max.y>=WithAABB.Min.y) and (AABB.Min.y<=WithAABB.Max.y));
end;
 function DynamicAABBTreeCreate:PDynamicAABBTree;
 var i:longint;
 begin
  GetMem(result,SizeOf(TDynamicAABBTree));
  FillChar(result^,SizeOf(TDynamicAABBTree),#0);
  result^.Root:=daabbtNULLNODE;
  result^.NodeCount:=0;
  result^.NodeCapacity:=16;
  GetMem(result^.Nodes,result^.NodeCapacity*SizeOf(TDynamicAABBTreeNode));
  FillChar(result^.Nodes^,result^.NodeCapacity*SizeOf(TDynamicAABBTreeNode),0);
  for i:=0 to result^.NodeCapacity-2 do begin
   result^.Nodes^[i].Next:=i+1;
   result^.Nodes^[i].Height:=-1;
  end;
  result^.Nodes^[result^.NodeCapacity-1].Next:=daabbtNULLNODE;
  result^.Nodes^[result^.NodeCapacity-1].Height:=-1;
  result^.FreeList:=0;
  result^.Path:=0;
  result^.InsertionCount:=0;
  result^.StackCapacity:=16;
  GetMem(result^.Stack,result^.StackCapacity*SizeOf(longint));
 end;
 procedure DynamicAABBTreeDestroy(const DynamicAABBTree:PDynamicAABBTree);
 begin
  FreeMem(DynamicAABBTree^.Nodes);
  FreeMem(DynamicAABBTree^.Stack);
  FreeMem(DynamicAABBTree);
 end;
 function DynamicAABBTreeAllocateNode(const DynamicAABBTree:PDynamicAABBTree):longint;
 var Node:PDynamicAABBTreeNode;
     i:longint;
 begin
  if DynamicAABBTree^.FreeList=daabbtNULLNODE then begin
   inc(DynamicAABBTree^.NodeCapacity,DynamicAABBTree^.NodeCapacity);
   ReallocMem(DynamicAABBTree^.Nodes,DynamicAABBTree^.NodeCapacity*SizeOf(TDynamicAABBTreeNode));
   FillChar(DynamicAABBTree^.Nodes^[DynamicAABBTree^.NodeCount],(DynamicAABBTree^.NodeCapacity-DynamicAABBTree^.NodeCount)*SizeOf(TDynamicAABBTreeNode),#0);
   for i:=DynamicAABBTree^.NodeCount to DynamicAABBTree^.NodeCapacity-2 do begin
    DynamicAABBTree^.Nodes^[i].Next:=i+1;
    DynamicAABBTree^.Nodes^[i].Height:=-1;
   end;
   DynamicAABBTree^.Nodes^[DynamicAABBTree^.NodeCapacity-1].Next:=daabbtNULLNODE;
   DynamicAABBTree^.Nodes^[DynamicAABBTree^.NodeCapacity-1].Height:=-1;
   DynamicAABBTree^.FreeList:=DynamicAABBTree^.NodeCount;
  end;
  result:=DynamicAABBTree^.FreeList;
  DynamicAABBTree^.FreeList:=DynamicAABBTree^.Nodes^[result].Next;
  Node:=@DynamicAABBTree^.Nodes^[result];
  Node^.Parent:=daabbtNULLNODE;
  Node^.Children[0]:=daabbtNULLNODE;
  Node^.Children[1]:=daabbtNULLNODE;
  Node^.Height:=0;
  Node^.UserData:=nil;
  inc(DynamicAABBTree^.NodeCount);
 end;
 procedure DynamicAABBTreeFreeNode(const DynamicAABBTree:PDynamicAABBTree;NodeID:longint);
 var Node:PDynamicAABBTreeNode;
 begin
  Node:=@DynamicAABBTree^.Nodes^[NodeID];
  Node^.Next:=DynamicAABBTree^.FreeList;
  Node^.Height:=-1;
  DynamicAABBTree^.FreeList:=NodeID;
  dec(DynamicAABBTree^.NodeCount);
 end;
 function DynamicAABBTreeBalance(const DynamicAABBTree:PDynamicAABBTree;NodeAID:longint):longint;
 var NodeA,NodeB,NodeC,NodeD,NodeE,NodeF,NodeG:PDynamicAABBTreeNode;
     NodeBID,NodeCID,NodeDID,NodeEID,NodeFID,NodeGID,NodeBalance:longint;
 begin
  NodeA:=@DynamicAABBTree^.Nodes^[NodeAID];
  if (NodeA.Children[0]<0) or (NodeA^.Height<2) then begin
   result:=NodeAID;
  end else begin
   NodeBID:=NodeA^.Children[0];
   NodeCID:=NodeA^.Children[1];
   NodeB:=@DynamicAABBTree^.Nodes^[NodeBID];
   NodeC:=@DynamicAABBTree^.Nodes^[NodeCID];
   NodeBalance:=NodeC^.Height-NodeB^.Height;
   if NodeBalance>1 then begin
    NodeFID:=NodeC^.Children[0];
    NodeGID:=NodeC^.Children[1];
    NodeF:=@DynamicAABBTree^.Nodes^[NodeFID];
    NodeG:=@DynamicAABBTree^.Nodes^[NodeGID];
    NodeC^.Children[0]:=NodeAID;
    NodeC^.Parent:=NodeA^.Parent;
    NodeA^.Parent:=NodeCID;
    if NodeC^.Parent>=0 then begin
     if DynamicAABBTree^.Nodes^[NodeC^.Parent].Children[0]=NodeAID then begin
      DynamicAABBTree^.Nodes^[NodeC^.Parent].Children[0]:=NodeCID;
     end else begin
      DynamicAABBTree^.Nodes^[NodeC^.Parent].Children[1]:=NodeCID;
     end;
    end else begin
     DynamicAABBTree^.Root:=NodeCID;
    end;
    if NodeF^.Height>NodeG^.Height then begin
     NodeC^.Children[1]:=NodeFID;
     NodeA^.Children[1]:=NodeGID;
     NodeG^.Parent:=NodeAID;
     NodeA^.AABB:=AABBCombine(NodeB^.AABB,NodeG^.AABB);
     NodeC^.AABB:=AABBCombine(NodeA^.AABB,NodeF^.AABB);
     NodeA^.Height:=1+Max(NodeB^.Height,NodeG^.Height);
     NodeC^.Height:=1+Max(NodeA^.Height,NodeF^.Height);
    end else begin
     NodeC^.Children[1]:=NodeGID;
     NodeA^.Children[1]:=NodeFID;
     NodeF^.Parent:=NodeAID;
     NodeA^.AABB:=AABBCombine(NodeB^.AABB,NodeF^.AABB);
     NodeC^.AABB:=AABBCombine(NodeA^.AABB,NodeG^.AABB);
     NodeA^.Height:=1+Max(NodeB^.Height,NodeF^.Height);
     NodeC^.Height:=1+Max(NodeA^.Height,NodeG^.Height);
    end;
    result:=NodeCID;
   end else if NodeBalance<-1 then begin
    NodeDID:=NodeB^.Children[0];
    NodeEID:=NodeB^.Children[1];
    NodeD:=@DynamicAABBTree^.Nodes^[NodeDID];
    NodeE:=@DynamicAABBTree^.Nodes^[NodeEID];
    NodeB^.Children[0]:=NodeAID;
    NodeB^.Parent:=NodeA^.Parent;
    NodeA^.Parent:=NodeBID;
    if NodeB^.Parent>=0 then begin
     if DynamicAABBTree^.Nodes^[NodeB^.Parent].Children[0]=NodeAID then begin
      DynamicAABBTree^.Nodes^[NodeB^.Parent].Children[0]:=NodeBID;
     end else begin
      DynamicAABBTree^.Nodes^[NodeB^.Parent].Children[1]:=NodeBID;
     end;
    end else begin
     DynamicAABBTree^.Root:=NodeBID;
    end;
    if NodeD^.Height>NodeE^.Height then begin
     NodeB^.Children[1]:=NodeDID;
     NodeA^.Children[0]:=NodeEID;
     NodeE^.Parent:=NodeAID;
     NodeA^.AABB:=AABBCombine(NodeC^.AABB,NodeE^.AABB);
     NodeB^.AABB:=AABBCombine(NodeA^.AABB,NodeD^.AABB);
     NodeA^.Height:=1+Max(NodeC^.Height,NodeE^.Height);
     NodeB^.Height:=1+Max(NodeA^.Height,NodeD^.Height);
    end else begin
     NodeB^.Children[1]:=NodeEID;
     NodeA^.Children[0]:=NodeDID;
     NodeD^.Parent:=NodeAID;
     NodeA^.AABB:=AABBCombine(NodeC^.AABB,NodeD^.AABB);
     NodeB^.AABB:=AABBCombine(NodeA^.AABB,NodeE^.AABB);
     NodeA^.Height:=1+Max(NodeC^.Height,NodeD^.Height);
     NodeB^.Height:=1+Max(NodeA^.Height,NodeE^.Height);
    end;
    result:=NodeBID;
   end else begin
    result:=NodeAID;
   end;
  end;
 end;
 procedure DynamicAABBTreeInsertLeaf(const DynamicAABBTree:PDynamicAABBTree;Leaf:longint);
 var Node:PDynamicAABBTreeNode;
     LeafAABB,CombinedAABB,AABB:TAABB;
     Index,Sibling,OldParent,NewParent:longint;
     Children:array[0..1] of longint;
     CombinedCost,Cost,InheritanceCost:double;
     Costs:array[0..1] of double;
 begin
  inc(DynamicAABBTree^.InsertionCount);
  if DynamicAABBTree^.Root<0 then begin
   DynamicAABBTree^.Root:=Leaf;
   DynamicAABBTree^.Nodes^[Leaf].Parent:=daabbtNULLNODE;
  end else begin
   LeafAABB:=DynamicAABBTree^.Nodes^[Leaf].AABB;
   Index:=DynamicAABBTree^.Root;
   while DynamicAABBTree^.Nodes^[Index].Children[0]>=0 do begin
    Children[0]:=DynamicAABBTree^.Nodes^[Index].Children[0];
    Children[1]:=DynamicAABBTree^.Nodes^[Index].Children[1];
    CombinedAABB:=AABBCombine(DynamicAABBTree^.Nodes^[Index].AABB,LeafAABB);
    CombinedCost:=AABBCost(CombinedAABB);
    Cost:=CombinedCost*2.0;
    InheritanceCost:=2.0*(CombinedCost-AABBCost(DynamicAABBTree^.Nodes^[Index].AABB));
    AABB:=AABBCombine(LeafAABB,DynamicAABBTree^.Nodes^[Children[0]].AABB);
    if DynamicAABBTree^.Nodes^[Children[0]].Children[0]<0 then begin
     Costs[0]:=AABBCost(AABB)+InheritanceCost;
    end else begin
     Costs[0]:=(AABBCost(AABB)-AABBCost(DynamicAABBTree^.Nodes^[Children[0]].AABB))+InheritanceCost;
    end;
    AABB:=AABBCombine(LeafAABB,DynamicAABBTree^.Nodes^[Children[1]].AABB);
    if DynamicAABBTree^.Nodes^[Children[1]].Children[1]<0 then begin
     Costs[1]:=AABBCost(AABB)+InheritanceCost;
    end else begin
     Costs[1]:=(AABBCost(AABB)-AABBCost(DynamicAABBTree^.Nodes^[Children[1]].AABB))+InheritanceCost;
    end;
    if (Cost<Costs[0]) and (Cost<Costs[1]) then begin
     break;
    end else begin
     if Costs[0]<Costs[1] then begin
      Index:=Children[0];
     end else begin
      Index:=Children[1];
     end;
    end;
   end;
   Sibling:=Index;
   OldParent:=DynamicAABBTree^.Nodes^[Sibling].Parent;
   NewParent:=DynamicAABBTreeAllocateNode(DynamicAABBTree);
   DynamicAABBTree^.Nodes^[NewParent].Parent:=OldParent;
   DynamicAABBTree^.Nodes^[NewParent].UserData:=nil;
   DynamicAABBTree^.Nodes^[NewParent].AABB:=AABBCombine(LeafAABB,DynamicAABBTree^.Nodes^[Sibling].AABB);
   DynamicAABBTree^.Nodes^[NewParent].Height:=DynamicAABBTree^.Nodes^[Sibling].Height+1;
   if OldParent>=0 then begin
    if DynamicAABBTree^.Nodes^[OldParent].Children[0]=Sibling then begin
     DynamicAABBTree^.Nodes^[OldParent].Children[0]:=NewParent;
    end else begin
     DynamicAABBTree^.Nodes^[OldParent].Children[1]:=NewParent;
    end;
    DynamicAABBTree^.Nodes^[NewParent].Children[0]:=Sibling;
    DynamicAABBTree^.Nodes^[NewParent].Children[1]:=Leaf;
    DynamicAABBTree^.Nodes^[Sibling].Parent:=NewParent;
    DynamicAABBTree^.Nodes^[Leaf].Parent:=NewParent;
   end else begin
    DynamicAABBTree^.Nodes^[NewParent].Children[0]:=Sibling;
    DynamicAABBTree^.Nodes^[NewParent].Children[1]:=Leaf;
    DynamicAABBTree^.Nodes^[Sibling].Parent:=NewParent;
    DynamicAABBTree^.Nodes^[Leaf].Parent:=NewParent;
    DynamicAABBTree^.Root:=NewParent;
   end;
   Index:=DynamicAABBTree^.Nodes^[Leaf].Parent;
   while Index>=0 do begin
    Index:=DynamicAABBTreeBalance(DynamicAABBTree,Index);
    Node:=@DynamicAABBTree^.Nodes^[Index];
    Node^.AABB:=AABBCombine(DynamicAABBTree^.Nodes^[Node^.Children[0]].AABB,DynamicAABBTree^.Nodes^[Node^.Children[1]].AABB);
    Node^.Height:=1+Max(DynamicAABBTree^.Nodes^[Node^.Children[0]].Height,DynamicAABBTree^.Nodes^[Node^.Children[1]].Height);
    Index:=Node^.Parent;
   end;  
  end;
 end;
 procedure DynamicAABBTreeRemoveLeaf(const DynamicAABBTree:PDynamicAABBTree;Leaf:longint);
 var Node:PDynamicAABBTreeNode;
     Parent,GrandParent,Sibling,Index:longint;
 begin
  if DynamicAABBTree^.Root=Leaf then begin
   DynamicAABBTree^.Root:=daabbtNULLNODE;
  end else begin
   Parent:=DynamicAABBTree^.Nodes^[Leaf].Parent;
   GrandParent:=DynamicAABBTree^.Nodes^[Parent].Parent;
   if DynamicAABBTree^.Nodes^[Parent].Children[0]=Leaf then begin
    Sibling:=DynamicAABBTree^.Nodes^[Parent].Children[1];
   end else begin
    Sibling:=DynamicAABBTree^.Nodes^[Parent].Children[0];
   end;
   if GrandParent>=0 then begin
    if DynamicAABBTree^.Nodes^[GrandParent].Children[0]=Parent then begin
     DynamicAABBTree^.Nodes^[GrandParent].Children[0]:=Sibling;
    end else begin
     DynamicAABBTree^.Nodes^[GrandParent].Children[1]:=Sibling;
    end;
    DynamicAABBTree^.Nodes^[Sibling].Parent:=GrandParent;
    DynamicAABBTreeFreeNode(DynamicAABBTree,Parent);
    Index:=GrandParent;
    while Index>=0 do begin
     Index:=DynamicAABBTreeBalance(DynamicAABBTree,Index);
     Node:=@DynamicAABBTree^.Nodes^[Index];
     Node^.AABB:=AABBCombine(DynamicAABBTree^.Nodes^[Node^.Children[0]].AABB,DynamicAABBTree^.Nodes^[Node^.Children[1]].AABB);
     Node^.Height:=1+Max(DynamicAABBTree^.Nodes^[Node^.Children[0]].Height,DynamicAABBTree^.Nodes^[Node^.Children[1]].Height);
     Index:=Node^.Parent;
    end;
   end else begin
    DynamicAABBTree^.Root:=Sibling;
    DynamicAABBTree^.Nodes^[Sibling].Parent:=daabbtNULLNODE;
    DynamicAABBTreeFreeNode(DynamicAABBTree,Parent);
   end;
  end;
 end;
 function DynamicAABBTreeCreateProxy(const DynamicAABBTree:PDynamicAABBTree;const AABB:TAABB;UserData:pointer):longint;
 var Node:PDynamicAABBTreeNode;
 begin
  result:=DynamicAABBTreeAllocateNode(DynamicAABBTree);
  Node:=@DynamicAABBTree^.Nodes^[result];
  Node^.AABB.Min.x:=AABB.Min.x-AABBEPSILON;
  Node^.AABB.Min.y:=AABB.Min.y-AABBEPSILON;
  Node^.AABB.Max.x:=AABB.Max.x+AABBEPSILON;
  Node^.AABB.Max.y:=AABB.Max.y+AABBEPSILON;
  Node^.UserData:=UserData;
  Node^.Height:=0;
  DynamicAABBTreeInsertLeaf(DynamicAABBTree,result);
 end;
 procedure DynamicAABBTreeDestroyProxy(const DynamicAABBTree:PDynamicAABBTree;NodeID:longint);
 begin
  DynamicAABBTreeRemoveLeaf(DynamicAABBTree,NodeID);
  DynamicAABBTreeFreeNode(DynamicAABBTree,NodeID);
 end;
 procedure VerticalAATreeClearNode(Root:PVerticalAATreeNode;var Node:PVerticalAATreeNode);
 begin
  if assigned(Node) then begin
   VerticalAATreeClearNode(Root,Node^.Left);
   VerticalAATreeClearNode(Root,Node^.Right);
   if Node<>Root then begin
    FreeMem(Node);
    Node:=nil;
   end;
  end;
 end;
 procedure VerticalAATreeInit(Root:PVerticalAATreeNode);
 begin
  FillChar(Root^,SizeOf(TVerticalAATreeNode),0);
  Root^.Level:=$7fffffff;
 end;
 procedure VerticalAATreeDone(Root:PVerticalAATreeNode);
 begin
  VerticalAATreeClearNode(Root,Root);
  Root^.Level:=$7fffffff;
 end;
 function VerticalAATreeFirst(Root:PVerticalAATreeNode):PVerticalAATreeNode;
 begin
  if assigned(Root^.Left) then begin
   result:=Root;
   while assigned(result^.Left) do begin
    result:=result^.Left;
   end;
  end else begin
   result:=nil;
  end;
 end;
 function VerticalAATreePrevious(Root,n:PVerticalAATreeNode):PVerticalAATreeNode;
 begin
  if assigned(n^.Left) then begin
   result:=n^.Left;
   while assigned(result^.Right) do begin
    result:=result^.Right;
   end;
  end else begin
   while assigned(n^.Parent) and (n^.Parent^.Left=n) do begin
    n:=n^.Parent;
   end;
   n:=n^.Parent;
   if assigned(n) and (n<>Root) then begin
    result:=n;
   end else begin
    result:=nil;
   end;
  end;
 end;
 function VerticalAATreeNext(Root,n:PVerticalAATreeNode):PVerticalAATreeNode;
 begin
  if assigned(n^.Right) then begin
   result:=n^.Right;
   while assigned(result^.Left) do begin
    result:=result^.Left;
   end;
  end else begin
   while assigned(n^.Parent) and (n^.Parent^.Right=n) do begin
    n:=n^.Parent;
   end;
   n:=n^.Parent;
   if assigned(n) and (n<>Root) then begin
    result:=n;
   end else begin
    result:=nil;
   end;
  end;
 end;
 procedure VerticalAATreeSkew(Root,OldParent:PVerticalAATreeNode);
 var NewParent:PVerticalAATreeNode;
 begin
  Assert(assigned(OldParent));
  NewParent:=OldParent^.Left;
  Assert(assigned(NewParent));
  if OldParent^.Parent^.Left=OldParent then begin
   OldParent^.Parent^.Left:=NewParent;
  end else begin
   OldParent^.Parent^.Right:=NewParent;
  end;
  NewParent^.Parent:=OldParent^.Parent;
  OldParent^.Parent:=NewParent;
  OldParent^.Left:=NewParent^.Right;
  if assigned(OldParent^.Left) then begin
   OldParent^.Left^.Parent:=OldParent;
  end;
  NewParent^.Right:=OldParent;
  if assigned(OldParent^.Left) then begin
   OldParent^.Level:=OldParent^.Left^.Level+1;
  end else begin
   OldParent^.Level:=1;
  end;
 end;
 function VerticalAATreeSplit(Root,OldParent:PVerticalAATreeNode):boolean;
 var NewParent:PVerticalAATreeNode;
 begin
  Assert(assigned(OldParent));
  NewParent:=OldParent^.Right;
  if assigned(NewParent) and assigned(NewParent^.Right) and (NewParent^.Right^.Level=OldParent^.Level) then begin
   if OldParent^.Parent^.Left=OldParent then begin
    OldParent^.Parent^.Left:=NewParent;
   end else begin
    OldParent^.Parent^.Right:=NewParent;
   end;
   NewParent^.Parent:=OldParent^.Parent;
   OldParent^.Parent:=NewParent;
   OldParent^.Right:=NewParent^.Left;
   if assigned(OldParent^.Right) then begin
    OldParent^.Right^.Parent:=OldParent;
   end;
   NewParent^.Left:=OldParent;
   NewParent^.Level:=OldParent^.Level+1;
   result:=true;
  end else begin
   result:=false;
  end;
 end;
 procedure VerticalAATreeRebalanceAfterLeafAdd(Root,n:PVerticalAATreeNode);
 begin
  n^.Level:=1;
  n^.Left:=nil;
  n^.Right:=nil;
  n:=n^.Parent;
  while n<>Root do begin
   if (assigned(n^.Left) and (n^.Level<>(n^.Left^.Level+1))) or ((not assigned(n^.Left)) and (n^.Level<>1)) then begin
    VerticalAATreeSkew(Root,n);
    if (not assigned(n^.Right)) or (n^.Level<>n^.Right^.Level) then begin
     n:=n^.Parent;
    end;
   end;
   if not VerticalAATreeSplit(Root,n^.Parent) then begin
    break;
   end;
   n:=n^.Parent;
  end;
 end;
 function VerticalAATreeFindNode(Root:PVerticalAATreeNode;YCoordinate:double):PVerticalAATreeNode;
 var n:PVerticalAATreeNode;
 begin
  result:=nil;
  n:=Root^.Left;
  while assigned(n) do begin
   if YCoordinate=n^.YCoordinate then begin
    result:=n;
    break;
   end else if YCoordinate<n^.YCoordinate then begin
    n:=n^.Left;
   end else begin
    n:=n^.Right;
   end;
  end;
 end;
 function VerticalAATreeInsertNode(Root:PVerticalAATreeNode;YCoordinate:double):PVerticalAATreeNode;
 var s:PVerticalAATreeNode;
     LessThan:boolean;
 begin
  result:=VerticalAATreeFindNode(Root,YCoordinate);
  if not assigned(result) then begin
   GetMem(result,SizeOf(TVerticalAATreeNode));
   FillChar(result^,SizeOf(TVerticalAATreeNode),0);
   result^.YCoordinate:=YCoordinate;
   result^.First:=nil;
   result^.Last:=nil;
   s:=Root;
   LessThan:=true;
   while (LessThan and assigned(s^.Left)) or ((not LessThan) and assigned(s^.Right)) do begin
    if LessThan then begin
     s:=s^.Left;
    end else begin
     s:=s^.Right;
    end;
    LessThan:=YCoordinate<s^.YCoordinate;
   end;
   if LessThan then begin
    s^.Left:=result;
   end else begin
    s^.Right:=result;
   end;
   result^.Parent:=s;
   VerticalAATreeRebalanceAfterLeafAdd(Root,result);
  end;
 end;
 function PointCreate(const x,y:double):PPoint;
 var Hash,HashIndex:longword;
 begin
  Hash:=HashXY(x,y);
  HashIndex:=Hash and HashMask;
  result:=PointHashTable[HashIndex];
  while assigned(result) and ((result^.Hash<>Hash) or not (ExactSameValue(result^.x,x) and ExactSameValue(result^.y,y))) do begin
   result:=result^.HashNext;
  end;
  if not assigned(result) then begin
   result:=GarbageCollectorAllocate(SizeOf(TPoint));
   result^.HashNext:=PointHashTable[HashIndex];
   PointHashTable[HashIndex]:=result;
   result^.Hash:=Hash;
   result^.Previous:=nil;
   result^.Next:=nil;
   result^.x:=x;
   result^.y:=y;
  end;
 end;
 function PointEquals(const a,b:PPoint):boolean;
 begin
  result:=(a=b) or (SameValue(a^.x,b^.x) and SameValue(a^.y,b^.y));
 end;
 function PointStrictlyEquals(const a,b:PPoint):boolean;
 begin
  result:=(a=b) or (ExactSameValue(a^.x,b^.x) and ExactSameValue(a^.y,b^.y));
 end;
 procedure EdgeUpdateMinMax(const Edge:PEdge);
 begin
  Edge^.MinX:=Min(Edge^.p^.x,Edge^.q^.x);
  Edge^.MaxX:=Max(Edge^.p^.x,Edge^.q^.x);
  Edge^.MinY:=Min(Edge^.p^.y,Edge^.q^.y);
  Edge^.MaxY:=Max(Edge^.p^.y,Edge^.q^.y);
 end;
 function EdgeCreate(const p,q:PPoint;DoAdd:boolean=true):PEdge;
 begin
  result:=GarbageCollectorAllocate(SizeOf(TEdge));
  if DoAdd then begin
   if assigned(Edges.Last) then begin
    Edges.Last^.Next:=result;
    result^.Previous:=Edges.Last;
   end else begin
    Edges.First:=result;
    result^.Previous:=nil;
   end;
   Edges.Last:=result;
   result^.Next:=nil;
  end else begin
   result^.Previous:=nil;
   result^.Next:=nil;
  end;
  result^.p:=p;
  result^.q:=q;
  EdgeUpdateMinMax(result);
 end;
 procedure EdgesAddEdge(const Edge:PEdge);
 begin
  if assigned(Edges.Last) then begin
   Edges.Last^.Next:=Edge;
   Edge^.Previous:=Edges.Last;
  end else begin
   Edges.First:=Edge;
   Edge^.Previous:=nil;
  end;
  Edges.Last:=Edge;
  Edge^.Next:=nil;
 end;               
 function EdgeCompare(const a,b:PEdge):longint;
 begin
  if SameValue(a^.MinY,b^.MinY) then begin
   if SameValue(a^.MinX,b^.MinX) then begin
    if SameValue(a^.MaxX,b^.MaxX) then begin
     result:=Sign(a^.MaxY-b^.MaxY);
    end else begin
     result:=Sign(a^.MaxX-b^.MaxX);
    end;
   end else begin
    result:=Sign(a^.MinX-b^.MinX);
   end;
  end else begin
   result:=Sign(a^.MinY-b^.MinY);
  end;
 end;
 procedure EdgesInsertAtEdge(AtEdge,Edge:PEdge);
 begin
  if EdgeCompare(AtEdge,Edge)<=0 then begin
   Edge^.Next:=AtEdge^.Next;
   Edge^.Previous:=AtEdge;
   if assigned(AtEdge^.Next) then begin
    AtEdge^.Next^.Previous:=Edge;
   end else begin
    Edges.Last:=Edge;
   end;
   AtEdge^.Next:=Edge;
  end else begin
   Edge^.Previous:=AtEdge^.Previous;
   Edge^.Next:=AtEdge;
   if assigned(AtEdge^.Previous) then begin
    AtEdge^.Previous^.Next:=Edge;
   end else begin
    Edges.First:=Edge;
   end;
   AtEdge^.Previous:=Edge;
  end;
 end;
 function QuadCreate:PQuad;
 begin
  result:=GarbageCollectorAllocate(SizeOf(TQuad));
  if assigned(Quads.Last) then begin
   Quads.Last^.Next:=result;
   result^.Previous:=Quads.Last;
  end else begin
   Quads.First:=result;
   result^.Previous:=nil;
  end;
  Quads.Last:=result;
  result^.Next:=nil;
 end;
 procedure QuadDestroy(const Quad:PQuad);
 begin
  if assigned(Quad^.Previous) then begin
   Quad^.Previous^.Next:=Quad^.Next;
  end else if Quads.First=Quad then begin
   Quads.First:=Quad^.Next;
  end;
  if assigned(Quad^.Next) then begin
   Quad^.Next^.Previous:=Quad^.Previous;
  end else if Quads.Last=Quad then begin
   Quads.Last:=Quad^.Previous;
  end;
  Quad^.Previous:=nil;
  Quad^.Next:=nil;
  GarbageCollectorFree(Quad);
 end;
 procedure QuadUpdateMinMax(const Quad:PQuad);
 begin
  Quad^.AABB.Min.x:=Min(Min(Min(Quad^.Points[0]^.x,Quad^.Points[1]^.x),Quad^.Points[2]^.x),Quad^.Points[3]^.x);
  Quad^.AABB.Max.x:=Max(Max(Max(Quad^.Points[0]^.x,Quad^.Points[1]^.x),Quad^.Points[2]^.x),Quad^.Points[3]^.x);
  Quad^.AABB.Min.y:=Min(Min(Min(Quad^.Points[0]^.y,Quad^.Points[1]^.y),Quad^.Points[2]^.y),Quad^.Points[3]^.y);
  Quad^.AABB.Max.y:=Max(Max(Max(Quad^.Points[0]^.y,Quad^.Points[1]^.y),Quad^.Points[2]^.y),Quad^.Points[3]^.y);
 end;
 function TriangleCreate:PTriangle;
 begin
  result:=GarbageCollectorAllocate(SizeOf(TTriangle));
  if assigned(Triangles.Last) then begin
   Triangles.Last^.Next:=result;
   result^.Previous:=Triangles.Last;
  end else begin
   Triangles.First:=result;
   result^.Previous:=nil;
  end;
  Triangles.Last:=result;
  result^.Next:=nil;
 end;
 procedure TriangleDestroy(const Triangle:PTriangle);
 begin
  if assigned(Triangle^.Previous) then begin
   Triangle^.Previous^.Next:=Triangle^.Next;
  end else if Triangles.First=Triangle then begin
   Triangles.First:=Triangle^.Next;
  end;
  if assigned(Triangle^.Next) then begin
   Triangle^.Next^.Previous:=Triangle^.Previous;
  end else if Triangles.Last=Triangle then begin
   Triangles.Last:=Triangle^.Previous;
  end;
  Triangle^.Previous:=nil;
  Triangle^.Next:=nil;
  GarbageCollectorFree(Triangle);
 end;
 procedure TriangleUpdateMinMax(const Triangle:PTriangle);
 begin
  Triangle^.AABB.Min.x:=Min(Min(Triangle^.Points[0]^.x,Triangle^.Points[1]^.x),Triangle^.Points[2]^.x);
  Triangle^.AABB.Max.x:=Max(Max(Triangle^.Points[0]^.x,Triangle^.Points[1]^.x),Triangle^.Points[2]^.x);
  Triangle^.AABB.Min.y:=Min(Min(Triangle^.Points[0]^.y,Triangle^.Points[1]^.y),Triangle^.Points[2]^.y);
  Triangle^.AABB.Max.y:=Max(Max(Triangle^.Points[0]^.y,Triangle^.Points[1]^.y),Triangle^.Points[2]^.y);
 end;
 procedure AddEdges;
 var i,j,Count:longint;
     LastPoint,Point:PPoint;
     InputPolygon:TBeRoTriangulationPolygon;
     LastInputPolygonPoint,InputPolygonPoint:PBeRoTriangulationPoint;
     Points:array of TBeRoTriangulationPoint;
 begin
  Points:=nil;
  try
   for i:=0 to length(InputPolygons)-1 do begin
    InputPolygon:=InputPolygons[i];
    if length(InputPolygon)>0 then begin
     Count:=0;
     LastInputPolygonPoint:=@InputPolygon[length(InputPolygon)-1];
     if length(Points)<length(InputPolygon) then begin
      SetLength(Points,length(InputPolygon)*2);
     end;
     for j:=0 to length(InputPolygon)-1 do begin
      InputPolygonPoint:=@InputPolygon[j];
      if not (SameValue(LastInputPolygonPoint^.x,InputPolygonPoint^.x) and SameValue(LastInputPolygonPoint^.y,InputPolygonPoint^.y)) then begin
       Points[Count]:=InputPolygonPoint^;
       inc(Count);
       LastInputPolygonPoint:=InputPolygonPoint;
      end;
     end;
     if Count>2 then begin
      LastInputPolygonPoint:=@Points[Count-1];
      LastPoint:=PointCreate(LastInputPolygonPoint^.x,LastInputPolygonPoint^.y);
      VerticalAATreeInsertNode(@VerticalAATreeNodeRoot,LastInputPolygonPoint^.y);
      for j:=0 to Count-1 do begin
       InputPolygonPoint:=@Points[j];
       VerticalAATreeInsertNode(@VerticalAATreeNodeRoot,InputPolygonPoint^.y);
       Point:=PointCreate(InputPolygonPoint^.x,InputPolygonPoint^.y);
       EdgeCreate(LastPoint,Point);
       LastPoint:=Point;
       LastInputPolygonPoint:=InputPolygonPoint;
      end;
     end;
    end;
   end;
  finally
   SetLength(Points,0);
  end;
 end;
 procedure MergeSortEdges;
 var PartA,PartB,CurrentEdge:PEdge;
     InSize,PartASize,PartBSize,Merges:longint;
 begin
  if assigned(Edges.First) then begin
   InSize:=1;
   while true do begin
    PartA:=Edges.First;
    Edges.First:=nil;
    Edges.Last:=nil;
    Merges:=0;
    while assigned(PartA) do begin
     inc(Merges);
     PartB:=PartA;
     PartASize:=0;
     while PartASize<InSize do begin
      inc(PartASize);
      PartB:=PartB^.Next;
      if not assigned(PartB) then begin
       break;
      end;
     end;
     PartBSize:=InSize;
     while (PartASize>0) or ((PartBSize>0) and assigned(PartB)) do begin
      if PartASize=0 then begin
       CurrentEdge:=PartB;
       PartB:=PartB^.Next;
       dec(PartBSize);
      end else if (PartBSize=0) or not assigned(PartB) then begin
       CurrentEdge:=PartA;
       PartA:=PartA^.Next;
       dec(PartASize);
      end else if EdgeCompare(PartA,PartB)<=0 then begin
       CurrentEdge:=PartA;
       PartA:=PartA^.Next;
       dec(PartASize);
      end else begin
       CurrentEdge:=PartB;
       PartB:=PartB^.Next;
       dec(PartBSize);
      end;
      if assigned(Edges.Last) then begin
       Edges.Last^.Next:=CurrentEdge;
      end else begin
       Edges.First:=CurrentEdge;
      end;
      CurrentEdge^.Previous:=Edges.Last;
      Edges.Last:=CurrentEdge;
     end;
     PartA:=PartB;
    end;
    Edges.Last^.Next:=nil;
    if Merges<=1 then begin
     break;
    end;
    inc(InSize,InSize);
   end;
  end;
 end;
 procedure InsertionSortEdges;
 var CurrentEdge,WorkEdge,TemporaryEdge,PreviousEdge,NextEdge:PEdge;
 begin
  CurrentEdge:=Edges.First;
  if assigned(CurrentEdge) then begin
   NextEdge:=CurrentEdge^.Next;
   if assigned(NextEdge) then begin
    CurrentEdge:=NextEdge;
    while assigned(CurrentEdge) do begin
     WorkEdge:=CurrentEdge;
     TemporaryEdge:=WorkEdge^.Previous;
     CurrentEdge:=CurrentEdge^.Next;
     if assigned(TemporaryEdge) and (EdgeCompare(TemporaryEdge,WorkEdge)>0) then begin
      repeat
       TemporaryEdge:=TemporaryEdge^.Previous;
      until not (assigned(TemporaryEdge) and (EdgeCompare(TemporaryEdge,WorkEdge)>0));
      PreviousEdge:=WorkEdge^.Previous;
      NextEdge:=WorkEdge^.Next;
      PreviousEdge^.Next:=NextEdge;
      if assigned(NextEdge) then begin
       NextEdge^.Previous:=PreviousEdge;
      end else if Edges.Last=WorkEdge then begin
       Edges.Last:=PreviousEdge;
      end;
      if assigned(TemporaryEdge) then begin
       PreviousEdge:=TemporaryEdge;
       TemporaryEdge:=TemporaryEdge^.Next;
       PreviousEdge^.Next:=WorkEdge;
       WorkEdge^.Previous:=PreviousEdge;
       TemporaryEdge^.Previous:=WorkEdge;
       WorkEdge^.Next:=TemporaryEdge;
      end else begin
       TemporaryEdge:=Edges.First;
       WorkEdge^.Previous:=nil;
       WorkEdge^.Next:=TemporaryEdge;
       TemporaryEdge^.Previous:=WorkEdge;
       Edges.First:=WorkEdge;
      end;
     end;
    end;
   end;
  end;
 end;
 procedure SplitEdgesAtIntersections;
 var EdgeA,EdgeB,SplittedEdge,NewEdges:PEdge;
     a0,a1,b0,b1,IntersectionPoint:PPoint;
     a10x,a10y,b10x,b10y,ab0x,ab0y,d,ai,bi,xi,yi:double;
     TryAgain:boolean;
 begin
  repeat
   NewEdges:=nil;
   EdgeA:=Edges.First;
   while assigned(EdgeA) do begin
    EdgeB:=EdgeA^.Next;
    while assigned(EdgeB) and (EdgeB^.MinY<=EdgeA^.MaxY) do begin
     if (EdgeA<>EdgeB) and ((EdgeA^.MinX<=EdgeB^.MaxX) and (EdgeB^.MinX<=EdgeA^.MaxX)) then begin
      a0:=EdgeA^.p;
      a1:=EdgeA^.q;
      b0:=EdgeB^.p;
      b1:=EdgeB^.q;
      a10x:=a1^.x-a0^.x;
      a10y:=a1^.y-a0^.y;
      b10x:=b1^.x-b0^.x;
      b10y:=b1^.y-b0^.y;
      d:=(a10x*b10y)-(b10x*a10y);
      if not IsZero(d) then begin
       ab0x:=a0^.x-b0^.x;
       ab0y:=a0^.y-b0^.y;
       ai:=((b10x*ab0y)-(b10y*ab0x))/d;
       if (ai>=0.0) and (ai<=1.0) then begin
        bi:=((a10x*ab0y)-(a10y*ab0x))/d;
        if (bi>=0.0) and (bi<=1.0) then begin
         d:=1.0-ai;
         xi:=(a0^.x*d)+(a1^.x*ai);
         yi:=(a0^.y*d)+(a1^.y*ai);
         if not ((SameValue(a0^.x,xi) and SameValue(a0^.y,yi)) or
                 (SameValue(a1^.x,xi) and SameValue(a1^.y,yi)) or
                 (SameValue(b0^.x,xi) and SameValue(b0^.y,yi)) or
                 (SameValue(b1^.x,xi) and SameValue(b1^.y,yi))) then begin
          VerticalAATreeInsertNode(@VerticalAATreeNodeRoot,yi);
          IntersectionPoint:=PointCreate(xi,yi);
          begin
           SplittedEdge:=EdgeCreate(IntersectionPoint,EdgeA^.q,false);
           EdgeA^.q:=SplittedEdge^.p;
           EdgeUpdateMinMax(EdgeA);
           SplittedEdge^.Next:=NewEdges;
           NewEdges:=SplittedEdge;
          end;
          begin
           SplittedEdge:=EdgeCreate(IntersectionPoint,EdgeB^.q,false);
           EdgeB^.q:=SplittedEdge^.p;
           EdgeUpdateMinMax(EdgeB);
           SplittedEdge^.Next:=NewEdges;
           NewEdges:=SplittedEdge;
          end;
         end;
        end;
       end;
      end;
     end;
     EdgeB:=EdgeB^.Next;
    end;
    EdgeA:=EdgeA^.Next;
   end;
   TryAgain:=assigned(NewEdges);
   while assigned(NewEdges) do begin
    SplittedEdge:=NewEdges;
    NewEdges:=SplittedEdge^.Next;
    EdgesAddEdge(SplittedEdge);
   end;
   InsertionSortEdges;
  until not TryAgain;
 end;
 procedure SplitEdgesAtYCoordinates;
 var Current,Start:PVerticalAATreeNode;
     Edge,SplittedEdge:PEdge;
     p0,p1:PPoint;
     TryAgain,EdgeTryAgain:boolean;
 begin
  repeat
   TryAgain:=false;
   Edge:=Edges.First;
   while assigned(Edge) do begin
    if Edge^.p^.y<Edge^.q^.y then begin
     p0:=Edge^.p;
     p1:=Edge^.q;
    end else begin
     p0:=Edge^.q;
     p1:=Edge^.p;
    end;
    Start:=VerticalAATreeNodeRoot.Left;
    while assigned(Start) and (Start^.YCoordinate>p0^.y) do begin
     Start:=VerticalAATreePrevious(@VerticalAATreeNodeRoot,Start);
    end;
    while assigned(Start) and (p0^.y<Start^.YCoordinate) do begin
     Start:=VerticalAATreeNext(@VerticalAATreeNodeRoot,Start);
    end;
    if not assigned(Start) then begin
     Start:=VerticalAATreeFirst(@VerticalAATreeNodeRoot);
    end;
    repeat
     Current:=Start;
     EdgeTryAgain:=false;
     while assigned(Current) and (Current^.YCoordinate<p1^.y) do begin
      if (p0^.y<p1^.y) and ((p0^.y<Current^.YCoordinate) and (Current^.YCoordinate<p1^.y)) and not (ExactSameValue(p0^.y,Current^.YCoordinate) or ExactSameValue(p1^.y,Current^.YCoordinate)) then begin
       SplittedEdge:=EdgeCreate(PointCreate(LerpValue(p0^.x,p1^.x,(Current^.YCoordinate-p0^.y)/(p1^.y-p0^.y)),Current^.YCoordinate),Edge^.q,false);
       Edge^.q:=SplittedEdge^.p;
       EdgeUpdateMinMax(Edge);
       EdgesInsertAtEdge(Edge,SplittedEdge);
       TryAgain:=true;
       EdgeTryAgain:=true;
       if Edge^.p^.y<Edge^.q^.y then begin
        p0:=Edge^.p;
        p1:=Edge^.q;
       end else begin
        p0:=Edge^.q;
        p1:=Edge^.p;
       end;
       break;
      end;
      Current:=VerticalAATreeNext(@VerticalAATreeNodeRoot,Current);
     end;
    until not EdgeTryAgain;
    Edge:=Edge^.Next;
   end;
   InsertionSortEdges;
  until not TryAgain;
 end;
 procedure SweepVertical;
  procedure SweepHorizontal(const MinY,MaxY:double;Edge:PEDge);
  var CurrentEdge,LastEdge:PEdge;
      Winding:longint;
      Flag,Draw:boolean;
      Quad:PQuad;
      AABB:TAABB;
  begin
   Winding:=0;
   Flag:=false;
   CurrentEdge:=Edge;
   LastEdge:=nil;
   while assigned(CurrentEdge) and (CurrentEdge^.MinY<MaxY) do begin
    if (CurrentEdge^.MinY<CurrentEdge^.MaxY) and ((CurrentEdge^.MinY<MaxY) and (MinY<CurrentEdge^.MaxY)) then begin
     if assigned(LastEdge) then begin
      case InputPolygonFillRule of
       btpfrEVENODD:begin
        Flag:=not Flag;
        Draw:=Flag;
       end;
       else begin
        if LastEdge^.p^.y<=LastEdge^.q^.y then begin
         inc(Winding);
        end else begin
         dec(Winding);
        end;
        case InputPolygonFillRule of
         btpfrNONZERO:begin
          Draw:=Winding<>0;
         end;
         btpfrPOSITIVE:begin
          Draw:=Winding>0;
         end;
         btpfrNEGATIVE:begin
          Draw:=Winding<0;
         end;
         else {btpfrABSGEQTWO:}begin
          Draw:=abs(Winding)>=2;
         end
        end;
       end;
      end;
      if Draw then begin    
       Quad:=QuadCreate;
       if LastEdge^.p^.y<LastEdge^.q^.y then begin
        Quad^.Points[TOPLEFT]:=LastEdge^.p;
        Quad^.Points[BOTTOMLEFT]:=LastEdge^.q;
       end else begin
        Quad^.Points[TOPLEFT]:=LastEdge^.q;
        Quad^.Points[BOTTOMLEFT]:=LastEdge^.p;
       end;
       if CurrentEdge^.p^.y<CurrentEdge^.q^.y then begin
        Quad^.Points[TOPRIGHT]:=CurrentEdge^.p;
        Quad^.Points[BOTTOMRIGHT]:=CurrentEdge^.q;
       end else begin
        Quad^.Points[TOPRIGHT]:=CurrentEdge^.q;
        Quad^.Points[BOTTOMRIGHT]:=CurrentEdge^.p;
       end;
       QuadUpdateMinMax(Quad);
       AABB.Min.x:=Quad^.AABB.Min.x-1.0;
       AABB.Min.y:=Quad^.AABB.Min.y-1.0;
       AABB.Max.x:=Quad^.AABB.Max.x+1.0;
       AABB.Max.y:=Quad^.AABB.Max.y+1.0;
       Quad^.Proxy:=DynamicAABBTreeCreateProxy(QuadDynamicAABBTree,AABB,Quad);
      end;
     end;
     LastEdge:=CurrentEdge;
    end;
    CurrentEdge:=CurrentEdge^.Next;
   end;
  end;
 var Current,Last:PVerticalAATreeNode;
     Edge:PEdge;
 begin
  Edge:=Edges.First;
  Last:=nil;
  Current:=VerticalAATreeFirst(@VerticalAATreeNodeRoot);
  while assigned(Current) do begin
   if assigned(Last) then begin
    while assigned(Edge) and (Edge^.MinY<Last^.YCoordinate) do begin
     Edge:=Edge^.Next;
    end;
    SweepHorizontal(Last^.YCoordinate,Current^.YCoordinate,Edge);
   end;
   Last:=Current;
   Current:=VerticalAATreeNext(@VerticalAATreeNodeRoot,Current);
  end;
 end;
 procedure AddTriangle(const ax,ay,bx,by,cx,cy:double); overload;
 var Triangle:PTriangle;
     AABB:TAABB;
 begin
  if DoOptimizeTriangles then begin
   Triangle:=TriangleCreate;
   Triangle^.Points[0]:=PointCreate(ax,ay);
   Triangle^.Points[1]:=PointCreate(bx,by);
   Triangle^.Points[2]:=PointCreate(cx,cy);
   TriangleUpdateMinMax(Triangle);
   AABB.Min.x:=Triangle^.AABB.Min.x-1.0;
   AABB.Min.y:=Triangle^.AABB.Min.y-1.0;
   AABB.Max.x:=Triangle^.AABB.Max.x+1.0;
   AABB.Max.y:=Triangle^.AABB.Max.y+1.0;
   Triangle^.Proxy:=DynamicAABBTreeCreateProxy(TriangleDynamicAABBTree,AABB,Triangle);
  end else begin
   OutputTriangle(ax,ay,bx,by,cx,cy);
  end;
 end;
 procedure AddTriangle(const a,b,c:PPoint); overload;
 var Triangle:PTriangle;
     AABB:TAABB;
 begin
  if DoOptimizeTriangles then begin
   Triangle:=TriangleCreate;
   Triangle^.Points[0]:=a;
   Triangle^.Points[1]:=b;
   Triangle^.Points[2]:=c;
   TriangleUpdateMinMax(Triangle);
   AABB.Min.x:=Triangle^.AABB.Min.x-1.0;
   AABB.Min.y:=Triangle^.AABB.Min.y-1.0;
   AABB.Max.x:=Triangle^.AABB.Max.x+1.0;
   AABB.Max.y:=Triangle^.AABB.Max.y+1.0;
   Triangle^.Proxy:=DynamicAABBTreeCreateProxy(TriangleDynamicAABBTree,AABB,Triangle);
  end else begin
   OutputTriangle(a^.x,a^.y,b^.x,b^.y,c^.x,c^.y);
  end;
 end;
 function SameDirection(const a0,a1,b0,b1:PPoint):boolean; overload;
 var ax,ay,al,bx,by,bl:double;
 begin
  result:=false;
  ax:=a1^.x-a0^.x;
  ay:=a1^.y-a0^.y;
  al:=sqr(ax)+sqr(ay);
  if not IsZero(al) then begin
   bx:=b1^.x-b0^.x;
   by:=b1^.y-b0^.y;
   bl:=sqr(bx)+sqr(by);
   if not IsZero(bl) then begin
    al:=sqrt(al);
    ax:=ax/al;
    ay:=ay/al;
    bl:=sqrt(bl);
    bx:=bx/bl;
    by:=by/bl;
    result:=SameValue(ax,bx) and SameValue(ay,by);
   end;
  end;
 end;
 function SameDirection(const a0,a1,b0,b1:PPolygonVertex):boolean; overload;
 var ax,ay,al,bx,by,bl:double;
 begin
  result:=false;
  ax:=a1^.x-a0^.x;
  ay:=a1^.y-a0^.y;
  al:=sqr(ax)+sqr(ay);
  if not IsZero(al) then begin
   bx:=b1^.x-b0^.x;
   by:=b1^.y-b0^.y;
   bl:=sqr(bx)+sqr(by);
   if not IsZero(bl) then begin
    al:=sqrt(al);
    ax:=ax/al;
    ay:=ay/al;
    bl:=sqrt(bl);
    bx:=bx/bl;
    by:=by/bl;
    result:=SameValue(ax,bx) and SameValue(ay,by);
   end;
  end;
 end;
 procedure OptimizeQuads;
  function MergeQuads(Quad,OtherQuad:PQuad):boolean;
  begin
   result:=false;
   if (PointEquals(Quad.Points[BOTTOMLEFT],OtherQuad.Points[TOPLEFT]) and PointEquals(Quad.Points[BOTTOMRIGHT],OtherQuad.Points[TOPRIGHT])) and
      (SameDirection(Quad.Points[TOPLEFT],Quad.Points[BOTTOMLEFT],OtherQuad.Points[TOPLEFT],OtherQuad.Points[BOTTOMLEFT]) and
       SameDirection(Quad.Points[TOPRIGHT],Quad.Points[BOTTOMRIGHT],OtherQuad.Points[TOPRIGHT],OtherQuad.Points[BOTTOMRIGHT])) then begin
    Quad.Points[BOTTOMLEFT]:=OtherQuad.Points[BOTTOMLEFT];
    Quad.Points[BOTTOMRIGHT]:=OtherQuad.Points[BOTTOMRIGHT];
    result:=true;
   end else if (PointEquals(Quad.Points[TOPLEFT],OtherQuad.Points[BOTTOMLEFT]) and PointEquals(Quad.Points[TOPRIGHT],OtherQuad.Points[BOTTOMRIGHT])) and
               (SameDirection(Quad.Points[BOTTOMLEFT],Quad.Points[TOPLEFT],OtherQuad.Points[BOTTOMLEFT],OtherQuad.Points[TOPLEFT]) and
                SameDirection(Quad.Points[BOTTOMRIGHT],Quad.Points[TOPRIGHT],OtherQuad.Points[BOTTOMRIGHT],OtherQuad.Points[TOPRIGHT])) then begin
    Quad.Points[TOPLEFT]:=OtherQuad.Points[TOPLEFT];
    Quad.Points[TOPRIGHT]:=OtherQuad.Points[TOPRIGHT];
    result:=true;
   end else if (PointEquals(Quad.Points[TOPRIGHT],OtherQuad.Points[TOPLEFT]) and PointEquals(Quad.Points[BOTTOMRIGHT],OtherQuad.Points[BOTTOMLEFT])) and
               (SameDirection(Quad.Points[TOPLEFT],Quad.Points[TOPRIGHT],OtherQuad.Points[TOPLEFT],OtherQuad.Points[TOPRIGHT]) and
                SameDirection(Quad.Points[BOTTOMLEFT],Quad.Points[BOTTOMRIGHT],OtherQuad.Points[BOTTOMLEFT],OtherQuad.Points[BOTTOMRIGHT])) then begin
    Quad.Points[TOPRIGHT]:=OtherQuad.Points[TOPRIGHT];
    Quad.Points[BOTTOMRIGHT]:=OtherQuad.Points[BOTTOMRIGHT];
    result:=true;
   end else if (PointEquals(Quad.Points[TOPLEFT],OtherQuad.Points[TOPRIGHT]) and PointEquals(Quad.Points[BOTTOMLEFT],OtherQuad.Points[BOTTOMRIGHT])) and
               (SameDirection(Quad.Points[TOPRIGHT],Quad.Points[TOPLEFT],OtherQuad.Points[TOPRIGHT],OtherQuad.Points[TOPLEFT]) and
                SameDirection(Quad.Points[BOTTOMRIGHT],Quad.Points[BOTTOMLEFT],OtherQuad.Points[BOTTOMRIGHT],OtherQuad.Points[BOTTOMLEFT])) then begin
    Quad.Points[TOPLEFT]:=OtherQuad.Points[TOPLEFT];
    Quad.Points[BOTTOMLEFT]:=OtherQuad.Points[BOTTOMLEFT];
    result:=true;
   end;
  end;
 var Quad,OtherQuad:PQuad;
     LocalStack:PDynamicAABBTreeLongintArray;
     LocalStackPointer,NodeID:longint;
     Node:PDynamicAABBTreeNode;
     TryAgain:boolean;
 begin
  Quad:=Quads.First;
  while assigned(Quad) do begin
   repeat
    TryAgain:=false;
    LocalStack:=QuadDynamicAABBTree^.Stack;
    LocalStack^[0]:=QuadDynamicAABBTree^.Root;
    LocalStackPointer:=1;
    while LocalStackPointer>0 do begin
     dec(LocalStackPointer);
     NodeID:=LocalStack^[LocalStackPointer];
     if NodeID>=0 then begin
      Node:=@QuadDynamicAABBTree^.Nodes[NodeID];
      if AABBIntersect(Node^.AABB,Quad^.AABB) then begin
       if Node^.Children[0]<0 then begin
        OtherQuad:=Node^.UserData;
        if Quad<>OtherQuad then begin
         if MergeQuads(Quad,OtherQuad) then begin
          DynamicAABBTreeDestroyProxy(QuadDynamicAABBTree,OtherQuad^.Proxy);
          QuadDestroy(OtherQuad);
          QuadUpdateMinMax(Quad);
          DynamicAABBTreeRemoveLeaf(QuadDynamicAABBTree,Quad^.Proxy);
          Node^.AABB.Min.x:=Quad^.AABB.Min.x-1.0;
          Node^.AABB.Min.y:=Quad^.AABB.Min.y-1.0;
          Node^.AABB.Max.x:=Quad^.AABB.Max.x+1.0;
          Node^.AABB.Max.y:=Quad^.AABB.Max.y+1.0;
          DynamicAABBTreeInsertLeaf(QuadDynamicAABBTree,Quad^.Proxy);
          TryAgain:=true;
          break;
         end;
        end;
       end else begin
        if QuadDynamicAABBTree^.StackCapacity<=(LocalStackPointer+2) then begin
         QuadDynamicAABBTree^.StackCapacity:=(LocalStackPointer+2)*2;
         ReallocMem(QuadDynamicAABBTree^.Stack,QuadDynamicAABBTree^.StackCapacity*SizeOf(longint));
         LocalStack:=QuadDynamicAABBTree^.Stack;
        end;
        LocalStack^[LocalStackPointer+0]:=Node^.Children[0];
        LocalStack^[LocalStackPointer+1]:=Node^.Children[1];
        inc(LocalStackPointer,2);
       end;
      end;
     end;
    end;
   until not TryAgain;
   Quad:=Quad^.Next;
  end;
 end;
 procedure TriangulateQuads;
  procedure SortQuadPoints(const Quad:PQuad);
   function CrossProduct(const a,b:PPoint):double;
   begin
    result:=(a^.x*b^.y)-(a^.y*b^.x);
   end;
   function Orientation(const a,b,c:PPoint):double;
   begin
    result:=CrossProduct(a,b)+CrossProduct(b,c)+CrossProduct(c,a);
   end;
   procedure Swap(a,b:PPPoint);
   var t:PPoint;
   begin
    t:=a^;
    a^:=b^;
    b^:=t;
   end;
  var a,b,c,d:PPPoint;
  begin        
   a:=@Quad.Points[0];
   b:=@Quad.Points[1];
   c:=@Quad.Points[2];
   d:=@Quad.Points[3];
   if Orientation(a^,b^,c^)<0.0 then begin
    if Orientation(a^,c^,d^)<0.0 then begin
     exit;
    end else if Orientation(a^,b^,d^)<0.0 then begin
     Swap(d,c);
    end else begin
     Swap(a,d);
    end;
   end else if Orientation(a^,c^,d^)<0.0 then begin
    if Orientation(a^,b^,d^)<0.0 then begin
     Swap(b,c);
    end else begin
     Swap(a,b);
    end;
   end else begin
    Swap(a,c);
   end;
  end;
 var Quad,Next:PQuad;
 begin
  Quad:=Quads.First;
  while assigned(Quad) do begin
   Next:=Quad^.Next;
   SortQuadPoints(Quad);
   AddTriangle(Quad^.Points[0],Quad^.Points[1],Quad^.Points[2]);
   AddTriangle(Quad^.Points[0],Quad^.Points[2],Quad^.Points[3]);
   QuadDestroy(Quad);
   Quad:=Next;
  end;
 end;
 procedure OptimizeTriangles;
  function MergeTriangles(Triangle,OtherTriangle:PTriangle):boolean;
  const ModuloNextThree:array[0..2] of longint=(1,2,0);
        RemainingVertex:array[0..2,0..2] of longint=((-1,2,1),
                                                     (2,-1,0),
                                                     (1,0,-1));
  var i,j,k,l,cb0,cb1,pa,pb,pc,SharedVertices:longint;
      Coincident,Parallel:boolean;
  begin
   result:=false;
   if Triangle<>OtherTriangle then begin
    Coincident:=false;
    Parallel:=false;
    cb0:=-1;
    cb1:=-1;
    pa:=-1;
    pb:=-1;
    SharedVertices:=0;
    for i:=0 to 2 do begin
     j:=ModuloNextThree[i];
     for k:=0 to 2 do begin
      l:=ModuloNextThree[k];
      if PointEquals(Triangle^.Points[i],OtherTriangle^.Points[k]) then begin
       inc(SharedVertices);
       if SharedVertices>2 then begin
        exit;
       end;
      end;
      if ((PointEquals(Triangle^.Points[i],OtherTriangle^.Points[k]) and PointEquals(Triangle^.Points[j],OtherTriangle^.Points[l]))) or
         ((PointEquals(Triangle^.Points[i],OtherTriangle^.Points[l]) and PointEquals(Triangle^.Points[j],OtherTriangle^.Points[k]))) then begin
       cb0:=k;
       cb1:=l;
       Coincident:=true;
      end else if PointEquals(Triangle^.Points[j],OtherTriangle^.Points[k]) and
                  SameDirection(Triangle^.Points[i],Triangle^.Points[j],OtherTriangle^.Points[k],OtherTriangle^.Points[l]) then begin
       pa:=j;
       pb:=k;
       Parallel:=true;
      end else if PointEquals(Triangle^.Points[j],OtherTriangle^.Points[l]) and
                  SameDirection(Triangle^.Points[i],Triangle^.Points[j],OtherTriangle^.Points[l],OtherTriangle^.Points[k]) then begin
       pa:=j;
       pb:=l;
       Parallel:=true;
      end else if PointEquals(Triangle^.Points[i],OtherTriangle^.Points[k]) and
                  SameDirection(Triangle^.Points[j],Triangle^.Points[i],OtherTriangle^.Points[k],OtherTriangle^.Points[l]) then begin
       pa:=i;
       pb:=k;
       Parallel:=true;
      end else if PointEquals(Triangle^.Points[i],OtherTriangle^.Points[l]) and
                  SameDirection(Triangle^.Points[j],Triangle^.Points[i],OtherTriangle^.Points[l],OtherTriangle^.Points[k]) then begin
       pa:=i;
       pb:=l;
       Parallel:=true;
      end;
     end;
    end;
    if (SharedVertices=2) and Coincident and Parallel then begin
     if pb=cb0 then begin
      pc:=RemainingVertex[pb,cb1];
     end else if pb=cb1 then begin
      pc:=RemainingVertex[pb,cb0];
     end else begin
      pc:=-1;
     end;
     if pc>=0 then begin
      Triangle^.Points[pa]:=OtherTriangle^.Points[pc];
      result:=true;
      exit;
     end;
    end;
   end;
  end;
 var Triangle,OtherTriangle:PTriangle;
     LocalStack:PDynamicAABBTreeLongintArray;
     LocalStackPointer,NodeID:longint;
     Node:PDynamicAABBTreeNode;
     TryAgain:boolean;
 begin
  Triangle:=Triangles.First;
  while assigned(Triangle) do begin
   repeat
    TryAgain:=false;
    LocalStack:=TriangleDynamicAABBTree^.Stack;
    LocalStack^[0]:=TriangleDynamicAABBTree^.Root;
    LocalStackPointer:=1;
    while LocalStackPointer>0 do begin
     dec(LocalStackPointer);
     NodeID:=LocalStack^[LocalStackPointer];
     if NodeID>=0 then begin
      Node:=@TriangleDynamicAABBTree^.Nodes[NodeID];
      if AABBIntersect(Node^.AABB,Triangle^.AABB) then begin
       if Node^.Children[0]<0 then begin
        OtherTriangle:=Node^.UserData;
        if Triangle<>OtherTriangle then begin
         if MergeTriangles(Triangle,OtherTriangle) then begin
          DynamicAABBTreeDestroyProxy(TriangleDynamicAABBTree,OtherTriangle^.Proxy);
          TriangleDestroy(OtherTriangle);
          TriangleUpdateMinMax(Triangle);
          DynamicAABBTreeRemoveLeaf(TriangleDynamicAABBTree,Triangle^.Proxy);
          Node^.AABB.Min.x:=Triangle^.AABB.Min.x-1.0;
          Node^.AABB.Min.y:=Triangle^.AABB.Min.y-1.0;
          Node^.AABB.Max.x:=Triangle^.AABB.Max.x+1.0;
          Node^.AABB.Max.y:=Triangle^.AABB.Max.y+1.0;
          DynamicAABBTreeInsertLeaf(TriangleDynamicAABBTree,Triangle^.Proxy);
          TryAgain:=true;
          //exit;
          break;
         end;
        end;
       end else begin
        if TriangleDynamicAABBTree^.StackCapacity<=(LocalStackPointer+2) then begin
         TriangleDynamicAABBTree^.StackCapacity:=(LocalStackPointer+2)*2;
         ReallocMem(TriangleDynamicAABBTree^.Stack,TriangleDynamicAABBTree^.StackCapacity*SizeOf(longint));
         LocalStack:=TriangleDynamicAABBTree^.Stack;
        end;
        LocalStack^[LocalStackPointer+0]:=Node^.Children[0];
        LocalStack^[LocalStackPointer+1]:=Node^.Children[1];
        inc(LocalStackPointer,2);
       end;
      end;
     end;
    end;
   until not TryAgain;
   Triangle:=Triangle^.Next;
  end;
 end;
 procedure PushOutputTriangles;
 var Triangle,Next:PTriangle;
 begin
  Triangle:=Triangles.First;
  while assigned(Triangle) do begin
   Next:=Triangle^.Next;
   OutputTriangle(Triangle^.Points[0]^.x,Triangle^.Points[0]^.y,Triangle^.Points[1]^.x,Triangle^.Points[1]^.y,Triangle^.Points[2]^.x,Triangle^.Points[2]^.y);
   Triangle:=Next;
  end;
 end;
 function PolygonVertexCreate(const x,y:double):PPolygonVertex;
 begin
  result:=GarbageCollectorAllocate(SizeOf(TPolygonVertex));
  result^.x:=x;
  result^.y:=y;
  result^.Previous:=nil;
  result^.Next:=nil;
  result^.Corresponding:=nil;
  result^.Distance:=0.0;
  result^.IsEntry:=true;
  result^.IsIntersection:=false;
  result^.Visited:=false;
 end;
 function PolygonVertexCreateIntersection(const x,y,Distance:double):PPolygonVertex;
 begin
  result:=PolygonVertexCreate(x,y);
  result^.Distance:=Distance;
  result^.IsEntry:=false;
  result^.IsIntersection:=true;
 end;
 procedure PolygonVertexMarkIntersection(var PolygonVertex:PPolygonVertex);
 begin
  PolygonVertex^.Distance:=0.0;
  PolygonVertex^.IsEntry:=false;
  PolygonVertex^.IsIntersection:=true;
 end;
 procedure PolygonVertexVisit(const PolygonVertex:PPolygonVertex);
 begin
  PolygonVertex^.Visited:=true;
  if assigned(PolygonVertex^.Corresponding) and not PolygonVertex^.Corresponding.Visited then begin
   PolygonVertexVisit(PolygonVertex^.Corresponding);
  end;
 end;
 function PolygonVertexEquals(const a,b:PPolygonVertex):boolean;
 begin
  result:=(a=b) or (SameValue(a^.x,b^.x) and SameValue(a^.y,b^.y));
 end;
 function PolygonVertexIsInside(const PolygonVertex:PPolygonVertex;const Polygon:PPolygon):boolean;
 var Vertex,Next:PPolygonVertex;
     x,y:double;
 begin
  result:=false;
  Vertex:=Polygon^.First;
  Next:=Vertex^.Next;
  x:=PolygonVertex^.x;
  y:=PolygonVertex^.y;
  repeat
   if (((Vertex^.y<y) and (Next^.y>=y)) or ((Next^.y<y) and (Vertex^.y>=y))) and ((Vertex^.x<=x) or (Next^.x<=x)) then begin
    result:=result xor ((Vertex^.x+((Next^.x-Vertex^.x)*((y-Vertex^.y)/(Next^.y-Vertex^.y))))<x);
   end;
   Vertex:=Vertex^.Next;
   Next:=Vertex^.Next;
   if not assigned(Next) then begin
    Next:=Polygon^.First;
   end;
  until Vertex=Polygon^.First;
 end;
 function PolygonCreate:PPolygon;
 begin
  result:=GarbageCollectorAllocate(SizeOf(TPolygon));
  result^.First:=nil;
  result^.Count:=0;
 end;
 function PolygonsCreate:PPolygons;
 begin
  result:=GarbageCollectorAllocate(SizeOf(TPolygons));
  result^.First:=nil;
  result^.Last:=nil;
  result^.Count:=0;
 end;
 procedure PolygonVertexDestroy(PolygonVertex:PPolygonVertex);
 begin
  GarbageCollectorFree(PolygonVertex);
 end;
 procedure PolygonDestroy(Polygon:PPolygon);
 var Current,Next:PPolygonVertex;
 begin
  Current:=Polygon^.First;
  while assigned(Current) do begin
   Next:=Current^.Next;
   PolygonVertexDestroy(Current);
   if Next=Polygon^.First then begin
    break;
   end else begin
    Current:=Next;
   end;
  end;
  Polygon^.First:=nil;
  GarbageCollectorFree(Polygon);
 end;
 procedure PolygonsDestroy(Polygons:PPolygons);
 var Current,Next:PPolygon;
 begin
  Current:=Polygons^.First;
  while assigned(Current) do begin
   Next:=Current^.Next;
   PolygonDestroy(Current);
   Current:=Next;
  end;
  Polygons^.First:=nil;
  Polygons^.Last:=nil;
  GarbageCollectorFree(Polygons);
 end;
 procedure PolygonsAdd(const Polygons:PPolygons;const Polygon:PPolygon);
 begin
  if assigned(Polygons^.Last) then begin
   Polygons^.Last^.Next:=Polygon;
   Polygon^.Previous:=Polygons^.Last;
  end else begin
   Polygons^.First:=Polygon;
   Polygon^.Previous:=nil;
  end;
  Polygons^.Last:=Polygon;
  Polygon^.Next:=nil;
  inc(Polygons^.Count);
 end;
 procedure PolygonAddVertex(const Polygon:PPolygon;const Vertex:PPolygonVertex);
 var Previous,Next:PPolygonVertex;
 begin
  if assigned(Polygon^.First) then begin
   Next:=Polygon^.First;
   Previous:=Next^.Previous;
   Next^.Previous:=Vertex;
   Vertex^.Next:=Next;
   Vertex^.Previous:=Previous;
   Previous^.Next:=Vertex;
  end else begin
   Polygon^.First:=Vertex;
   Polygon^.First^.Next:=Vertex;
   Polygon^.First^.Previous:=Vertex;
  end;
  inc(Polygon^.Count);
 end;
 function PolygonClone(const Polygon:PPolygon):PPolygon;
 var Current:PPolygonVertex;
 begin
  result:=PolygonCreate;
  Current:=Polygon^.First;
  if assigned(Current) then begin
   repeat
    PolygonAddVertex(result,PolygonVertexCreate(Current^.x,Current^.y));
    Current:=Current^.Next;
   until Current=Polygon^.First;
  end;
 end;
 procedure PolygonInsertVertex(const Polygon:PPolygon;const Vertex,Start,Stop:PPolygonVertex);
 var Previous,Current:PPolygonVertex;
 begin
  Current:=Start;
  while (assigned(Start) and (Current<>Stop)) and (Current^.Distance<Vertex^.Distance) do begin
   Current:=Current^.Next;
  end;
  Vertex^.Next:=Current;
  Previous:=Current^.Previous;
  Vertex^.Previous:=Previous;
  Previous^.Next:=Vertex;
  Current^.Previous:=Vertex;
  inc(Polygon^.Count);
 end;
 function PolygonGetNext(const Polygon:PPolygon;const Vertex:PPolygonVertex):PPolygonVertex;
 begin
  result:=Vertex;
  while assigned(result) and result^.IsIntersection do begin
   result:=result^.Next;
  end;
 end;
 function PolygonGetFirstIntersect(const Polygon:PPolygon):PPolygonVertex;
 begin
  result:=Polygon^.First;
  repeat
   if result^.IsIntersection and not result^.Visited then begin
    exit;
   end;
   result:=result^.Next;
  until result=Polygon^.First;
 end;
 function PolygonHasUnprocessed(const Polygon:PPolygon):boolean;
 var Vertex:PPolygonVertex;
 begin
  result:=false;
  Vertex:=Polygon^.First;
  repeat
   if Vertex^.IsIntersection and not Vertex^.Visited then begin
    result:=true;
    exit;
   end;
   Vertex:=Vertex^.Next;
  until Vertex=Polygon^.First;
 end;
 procedure PolygonClean(const Polygon:PPolygon);
 var Vertex,Previous,Next:PPolygonVertex;
 begin
  Vertex:=Polygon^.First;
  if assigned(Vertex) then begin
   Vertex:=Vertex^.Next;
   while Vertex<>Polygon^.First do begin
    Previous:=Vertex^.Previous;
    Next:=Vertex^.Next;
    if (Vertex<>Next) and (SameValue(Vertex^.x,Next^.x) and SameValue(Vertex^.y,Next^.y)) then begin
     dec(Polygon^.Count);
     if Polygon^.First=Vertex then begin
      Polygon^.First:=Next;
     end;
     Next^.Previous:=Previous;
     Previous^.Next:=Next;
     PolygonVertexDestroy(Vertex);
    end;
    Vertex:=Next;
   end;
  end;
 end;
 procedure PolygonOptimize(const Polygon:PPolygon);
 var Vertex,Previous,Next:PPolygonVertex;
 begin
  Vertex:=Polygon^.First;
  if assigned(Vertex) then begin
   Vertex:=Vertex^.Next;
   while Vertex<>Polygon^.First do begin
    Previous:=Vertex^.Previous;
    Next:=Vertex^.Next;
    if (Vertex<>Next) and
       (((SameValue(Vertex^.x,Next^.x) and SameValue(Vertex^.y,Next^.y))) or
       (((Vertex<>Previous) and (Previous<>Next)) and SameDirection(Previous,Vertex,Vertex,Next))) then begin
     dec(Polygon^.Count);
     if Polygon^.First=Vertex then begin
      Polygon^.First:=Next;
     end;
     Next^.Previous:=Previous;
     Previous^.Next:=Next;
     PolygonVertexDestroy(Vertex);
    end;
    Vertex:=Next;
   end;
  end;
 end;
 function PolygonTouching(const Source,Clip:PPolygon):boolean;
  function Distance(const px,py,ax,ay,bx,by:double):double;
  var l,t,x,y:double;
  begin
   l:=sqr(ax-bx)+sqr(ay-by);
   if IsZero(l) then begin
    result:=sqr(px-ax)+sqr(py-ay);
   end else begin
    t:=(((px-ax)*(bx*ax))+((py-ay)*(by*ay)))/l;
    if t<0.0 then begin
     result:=sqr(px-ax)+sqr(py-ay);
    end else if t>1.0 then begin
     result:=sqr(px-bx)+sqr(py-by);
    end else begin
     l:=1.0-t;
     x:=(ax*l)+(bx*t);
     y:=(ay*l)+(by*t);
     result:=sqr(px-x)+sqr(py-y);
    end;
   end;
   if not IsZero(result) then begin
    result:=sqrt(result);
   end;
  end;
 var SourceVertex,NextSourceVertex,ClipVertex,NextClipVertex:PPolygonVertex;
 begin
  result:=false;
  SourceVertex:=Source^.First;
  ClipVertex:=Clip^.First;
  repeat
   NextSourceVertex:=SourceVertex^.Next;
   repeat
    NextClipVertex:=ClipVertex^.Next;
    if ((IsZero(Distance(SourceVertex^.x,SourceVertex^.y,ClipVertex^.x,ClipVertex^.y,NextClipVertex^.x,NextClipVertex^.y)) and not
        ((SameValue(SourceVertex^.x,ClipVertex^.x) and SameValue(SourceVertex^.y,ClipVertex^.y)) or
         (SameValue(SourceVertex^.x,NextClipVertex^.x) and SameValue(SourceVertex^.y,NextClipVertex^.y))))) or
       ((IsZero(Distance(ClipVertex^.x,ClipVertex^.y,SourceVertex^.x,SourceVertex^.y,NextSourceVertex^.x,NextSourceVertex^.y)) and not
        ((SameValue(ClipVertex^.x,SourceVertex^.x) and SameValue(ClipVertex^.y,SourceVertex^.y)) or
         (SameValue(ClipVertex^.x,NextSourceVertex^.x) and SameValue(ClipVertex^.y,NextSourceVertex^.y))))) then begin
     result:=true;
     exit;
    end else if ((SameValue(SourceVertex^.x,ClipVertex^.x) and SameValue(SourceVertex^.y,ClipVertex^.y)) and
                 (SameValue(NextSourceVertex^.x,NextClipVertex^.x) and SameValue(NextSourceVertex^.y,NextClipVertex^.y))) or
                ((SameValue(SourceVertex^.x,NextClipVertex^.x) and SameValue(SourceVertex^.y,NextClipVertex^.y)) and
                 (SameValue(NextSourceVertex^.x,ClipVertex^.x) and SameValue(NextSourceVertex^.y,ClipVertex^.y))) then begin
     result:=true;
     exit;
    end;
    ClipVertex:=NextClipVertex;
   until ClipVertex=Clip^.First;
   SourceVertex:=NextSourceVertex;
  until SourceVertex=Source^.First;
 end;
 function PolygonIntersecting(const Source,Clip:PPolygon):boolean;
 var SourceVertex,ClipVertex,s0,s1,c0,c1:PPolygonVertex;
     s10x,s10y,c10x,c10y,sc0x,sc0y,d,ToSource,ToClip,IntersectionX,IntersectionY:double;
 begin
  result:=false;
  SourceVertex:=Source^.First;
  ClipVertex:=Clip^.First;
  repeat
   if not SourceVertex^.IsIntersection then begin
    repeat
     if not ClipVertex^.IsIntersection then begin
      s0:=SourceVertex;
      s1:=PolygonGetNext(Source,SourceVertex^.Next);
      c0:=ClipVertex;
      c1:=PolygonGetNext(Clip,ClipVertex^.Next);
      s10x:=s1^.x-s0^.x;
      s10y:=s1^.y-s0^.y;
      c10x:=c1^.x-c0^.x;
      c10y:=c1^.y-c0^.y;
      d:=(s10x*c10y)-(c10x*s10y);
      if not IsZero(d) then begin
       sc0x:=s0^.x-c0^.x;
       sc0y:=s0^.y-c0^.y;
       ToSource:=((c10x*sc0y)-(sc0x*c10y))/d;
       ToClip:=((s10x*sc0y)-(sc0x*s10y))/d;
       if ((ToSource>0.0) and (ToSource<1.0)) and ((ToClip>0.0) and (ToClip<1.0)) then begin
        d:=1.0-ToSource;
        IntersectionX:=(s0.x*d)+(s1.x*ToSource);
        IntersectionY:=(s0.y*d)+(s1.y*ToSource);
        if not ((SameValue(IntersectionX,s0^.x) and SameValue(IntersectionY,s0^.y)) or
                (SameValue(IntersectionX,s1^.x) and SameValue(IntersectionY,s1^.y)) or
                (SameValue(IntersectionX,c0^.x) and SameValue(IntersectionY,c0^.y)) or
                (SameValue(IntersectionX,c1^.x) and SameValue(IntersectionY,c1^.y))) then begin
         result:=true;
         exit;
        end;
       end;
      end;
     end;
     ClipVertex:=ClipVertex^.Next;
    until ClipVertex=Clip^.First;
   end;
   SourceVertex:=SourceVertex^.Next;
  until SourceVertex=Source^.First;
 end;
 function PolygonSelfIntersecting(const Source:PPolygon):boolean;
 begin
  result:=PolygonIntersecting(Source,Source);
 end;
 function PolygonClip(const Source,Clip:PPolygon;SourceForwards,ClipForwards:boolean):PPolygons;
 var SourceVertex,ClipVertex,SourceIntersection,ClipIntersection,Current,s0,s1,c0,c1:PPolygonVertex;
     Clipped:PPolygon;
     s10x,s10y,c10x,c10y,sc0x,sc0y,d,ToSource,ToClip,IntersectionX,IntersectionY:double;
 begin
  result:=nil;
  SourceVertex:=Source^.First;
  repeat
   if not SourceVertex^.IsIntersection then begin
    ClipVertex:=Clip^.First;
    repeat
     if not ClipVertex^.IsIntersection then begin
      s0:=SourceVertex;
      s1:=PolygonGetNext(Source,SourceVertex^.Next);
      c0:=ClipVertex;
      c1:=PolygonGetNext(Clip,ClipVertex^.Next);
      if (SameValue(s0^.x,c0^.x) and SameValue(s0^.y,c0^.y)) and (SameValue(s1^.x,c1^.x) and SameValue(s1^.y,c1^.y)) then begin
       PolygonVertexMarkIntersection(s0);
       PolygonVertexMarkIntersection(c0);
       PolygonVertexMarkIntersection(s1);
       PolygonVertexMarkIntersection(c1);
       s0^.Corresponding:=c0;
       c0^.Corresponding:=s0;
       s1^.Corresponding:=c1;
       c1^.Corresponding:=s1;
      end else if (SameValue(s0^.x,c1^.x) and SameValue(s0^.y,c1^.y)) and (SameValue(s1^.x,c0^.x) and SameValue(s1^.y,c0^.y)) then begin
       PolygonVertexMarkIntersection(s0);
       PolygonVertexMarkIntersection(c0);
       PolygonVertexMarkIntersection(s1);
       PolygonVertexMarkIntersection(c1);
       s0^.Corresponding:=c1;
       c0^.Corresponding:=s1;
       s1^.Corresponding:=c0;
       c1^.Corresponding:=s0;
      end else begin
       s10x:=s1^.x-s0^.x;
       s10y:=s1^.y-s0^.y;
       c10x:=c1^.x-c0^.x;
       c10y:=c1^.y-c0^.y;
       d:=(s10x*c10y)-(c10x*s10y);
       if not IsZero(d) then begin
        sc0x:=s0^.x-c0^.x;
        sc0y:=s0^.y-c0^.y;
        ToSource:=((c10x*sc0y)-(sc0x*c10y))/d;
        ToClip:=((s10x*sc0y)-(sc0x*s10y))/d;
        if ((ToSource>0.0) and (ToSource<1.0)) and ((ToClip>0.0) and (ToClip<1.0)) then begin
         d:=1.0-ToSource;
         IntersectionX:=(s0.x*d)+(s1.x*ToSource);
         IntersectionY:=(s0.y*d)+(s1.y*ToSource);
         if not ((SameValue(IntersectionX,s0^.x) and SameValue(IntersectionY,s0^.y)) or
                 (SameValue(IntersectionX,s1^.x) and SameValue(IntersectionY,s1^.y)) or
                 (SameValue(IntersectionX,c0^.x) and SameValue(IntersectionY,c0^.y)) or
                 (SameValue(IntersectionX,c1^.x) and SameValue(IntersectionY,c1^.y))) then begin
          SourceIntersection:=PolygonVertexCreateIntersection(IntersectionX,IntersectionY,ToSource);
          ClipIntersection:=PolygonVertexCreateIntersection(IntersectionX,IntersectionY,ToClip);
          SourceIntersection^.Corresponding:=ClipIntersection;
          ClipIntersection^.Corresponding:=SourceIntersection;
          PolygonInsertVertex(Source,SourceIntersection,SourceVertex,PolygonGetNext(Source,SourceVertex^.Next));
          PolygonInsertVertex(Clip,ClipIntersection,ClipVertex,PolygonGetNext(Clip,ClipVertex^.Next));
         end;
        end else if ((IsZero(ToSource) or SameValue(ToSource,1.0)) and ((ToClip>0.0) and (ToClip<1.0))) or
                    (((ToSource>0.0) and (ToSource<1.0)) and (IsZero(ToClip) or SameValue(ToClip,1.0))) then begin
         exit;
        end;
       end;
      end;
     end;
     ClipVertex:=ClipVertex^.Next;
    until ClipVertex=Clip^.First;
   end;
   SourceVertex:=SourceVertex^.Next;
  until SourceVertex=Source^.First;
  SourceVertex:=Source^.First;
  ClipVertex:=Clip^.First;
  SourceForwards:=SourceForwards xor PolygonVertexIsInside(SourceVertex,Clip);
  ClipForwards:=ClipForwards xor PolygonVertexIsInside(ClipVertex,Source);
  repeat
   if SourceVertex^.IsIntersection then begin
    SourceVertex^.IsEntry:=SourceForwards;
    SourceForwards:=not SourceForwards;
   end;
   SourceVertex:=SourceVertex^.Next;
  until SourceVertex=Source^.First;
  repeat
   if ClipVertex^.IsIntersection then begin
    ClipVertex^.IsEntry:=ClipForwards;
    ClipForwards:=not ClipForwards;
   end;
   ClipVertex:=ClipVertex^.Next;
  until ClipVertex=Clip^.First;
  while PolygonHasUnprocessed(Source) do begin
   Current:=PolygonGetFirstIntersect(Source);
   Clipped:=PolygonCreate;
   PolygonAddVertex(Clipped,PolygonVertexCreate(Current^.x,Current^.y));
   repeat                                    
    PolygonVertexVisit(Current);
    if Current^.IsEntry then begin
     repeat
      Current:=Current^.Next;
      PolygonAddVertex(Clipped,PolygonVertexCreate(Current^.x,Current^.y));
     until Current^.IsIntersection;
    end else begin
     repeat
      Current:=Current^.Previous;
      PolygonAddVertex(Clipped,PolygonVertexCreate(Current^.x,Current^.y));
     until Current^.IsIntersection;
    end;
    Current:=Current^.Corresponding;
   until Current^.Visited;
   if not assigned(result) then begin
    result:=PolygonsCreate;
   end;
   PolygonsAdd(result,Clipped);
  end;
 end;
 function PolygonUnion(const Source,Clip:PPolygon):PPolygons;
 begin
  result:=PolygonClip(Source,Clip,false,false);
 end;
 function PolygonIntersection(const Source,Clip:PPolygon):PPolygons;
 begin
  result:=PolygonClip(Source,Clip,true,true);
 end;
 function PolygonDifference(const Source,Clip:PPolygon):PPolygons;
 begin
  result:=PolygonClip(Source,Clip,false,true);
 end;
 function PolygonTryMerge(var a:PPolygon;const b:PPolygon):boolean;
 var Polygons:PPolygons;
     ca,cb:PPolygon;
 begin
  result:=false;
  ca:=PolygonClone(a);
  cb:=PolygonClone(b);
  if assigned(ca) and assigned(cb) then begin
   PolygonClean(ca);
   PolygonClean(cb);
   if PolygonTouching(ca,cb) and not PolygonIntersecting(ca,cb) then begin
    Polygons:=PolygonUnion(ca,cb);
    if assigned(Polygons) then begin
     if assigned(Polygons^.First) and (Polygons^.First=Polygons^.Last) and not PolygonSelfIntersecting(Polygons^.First) then begin
      PolygonDestroy(a);
      a:=Polygons^.First;
      Polygons^.First:=nil;
      Polygons^.Last:=nil;
      result:=true;
     end;
     PolygonsDestroy(Polygons);
    end;
   end;
  end;
  if assigned(ca) then begin
   PolygonDestroy(ca);
  end;
  if assigned(cb) then begin
   PolygonDestroy(cb);
  end;
 end;
 function MergeQuadsToPolygons:PPolygons;
 var Quad,OtherQuad:PQuad;
     LocalStack:PDynamicAABBTreeLongintArray;
     LocalStackPointer,NodeID:longint;
     Node:PDynamicAABBTreeNode;
     TryAgain:boolean;
     CurrentPolygon,OtherPolygon:PPolygon;
     AABB:TAABB;
 begin
  result:=PolygonsCreate;
  Quad:=Quads.First;
  while assigned(Quad) do begin
   Quad^.Polygon:=PolygonCreate;
   PolygonAddVertex(Quad^.Polygon,PolygonVertexCreate(Quad^.Points[0]^.x,Quad^.Points[0]^.y));
   PolygonAddVertex(Quad^.Polygon,PolygonVertexCreate(Quad^.Points[1]^.x,Quad^.Points[1]^.y));
   PolygonAddVertex(Quad^.Polygon,PolygonVertexCreate(Quad^.Points[2]^.x,Quad^.Points[2]^.y));
   PolygonAddVertex(Quad^.Polygon,PolygonVertexCreate(Quad^.Points[3]^.x,Quad^.Points[3]^.y));
   Quad:=Quad^.Next;
  end;
  Quad:=Quads.First;
  while assigned(Quad) do begin
   AABB:=Quad^.AABB;
   CurrentPolygon:=Quad^.Polygon;
   Quad^.Polygon:=nil;
   DynamicAABBTreeDestroyProxy(QuadDynamicAABBTree,Quad^.Proxy);
   QuadDestroy(Quad);
   repeat
    TryAgain:=false;
    LocalStack:=QuadDynamicAABBTree^.Stack;
    LocalStack^[0]:=QuadDynamicAABBTree^.Root;
    LocalStackPointer:=1;
    while LocalStackPointer>0 do begin
     dec(LocalStackPointer);
     NodeID:=LocalStack^[LocalStackPointer];
     if NodeID>=0 then begin
      Node:=@QuadDynamicAABBTree^.Nodes[NodeID];
      if AABBIntersect(Node^.AABB,AABB) then begin
       if Node^.Children[0]<0 then begin
        OtherQuad:=Node^.UserData;
        if assigned(OtherQuad) then begin
         OtherPolygon:=OtherQuad^.Polygon;
         if PolygonTryMerge(CurrentPolygon,OtherPolygon) then begin
          PolygonDestroy(OtherPolygon);
          AABB:=AABBCombine(AABB,Node^.AABB);
          OtherQuad^.Polygon:=nil;
          DynamicAABBTreeDestroyProxy(QuadDynamicAABBTree,OtherQuad^.Proxy);
          QuadDestroy(OtherQuad);
          TryAgain:=true;
          break;
         end;
        end;
       end else begin
        if QuadDynamicAABBTree^.StackCapacity<=(LocalStackPointer+2) then begin
         QuadDynamicAABBTree^.StackCapacity:=(LocalStackPointer+2)*2;
         ReallocMem(QuadDynamicAABBTree^.Stack,QuadDynamicAABBTree^.StackCapacity*SizeOf(longint));
         LocalStack:=QuadDynamicAABBTree^.Stack;
        end;
        LocalStack^[LocalStackPointer+0]:=Node^.Children[0];
        LocalStack^[LocalStackPointer+1]:=Node^.Children[1];
        inc(LocalStackPointer,2);
       end;
      end;
     end;
    end;
   until not TryAgain;
   PolygonsAdd(result,CurrentPolygon);
   Quad:=Quads.First;
  end;
 end;
 procedure PolygonTriangulate(const Polygon:PPolygon);
 const EPSILON=1e-14;
 type PVertex=^TVertex;
      TVertex=record
       ArrayPrevious:PVertex;
       ArrayNext:PVertex;
       Previous:PVertex;
       Next:PVertex;
       x,y:double;
      end;
      TVertices=array of TVertex;
 var Vertices:TVertices;
     i:longint;
     PolygonVertex,NextPolygonVertex:PPolygonVertex;
     Vertex,Current,Start,Stop,Next:PVertex;
     a,dx0,dy0,dx1,dy1,px,py,vx,vy,nx,ny,vnx,vny,npx,npy,pvx,pvy,cx,cy,pcx,pcy,vcx,vcy,ncx,ncy:double;
     IsEar,Convex,FirstSign:boolean;
 begin
  PolygonOptimize(Polygon);
  if Polygon^.Count=3 then begin
   AddTriangle(Polygon^.First^.x,Polygon^.First^.y,Polygon^.First^.Next^.x,Polygon^.First^.Next^.y,Polygon^.First^.Next^.Next^.x,Polygon^.First^.Next^.Next^.y);
  end else if Polygon^.Count>3 then begin
   Convex:=true;
   FirstSign:=false;
   PolygonVertex:=Polygon^.First;
   repeat
    NextPolygonVertex:=PolygonVertex^.Next;
    dx0:=NextPolygonVertex^.Next^.x-NextPolygonVertex^.x;
    dy0:=NextPolygonVertex^.Next^.y-NextPolygonVertex^.y;
    dx1:=PolygonVertex^.x-NextPolygonVertex^.x;
    dy1:=PolygonVertex^.y-NextPolygonVertex^.y;
    a:=(dx0*dy1)-(dx1*dy0);
    if PolygonVertex=Polygon^.First then begin
     FirstSign:=a<0.0;
    end else begin
     if FirstSign<>(a<0.0) then begin
      Convex:=false;
      break;
     end;
    end;
    PolygonVertex:=NextPolygonVertex;
   until PolygonVertex=Polygon^.First;
   if Convex then begin
    PolygonVertex:=Polygon^.First^.Next;
    repeat
     NextPolygonVertex:=PolygonVertex^.Next;
     AddTriangle(Polygon^.First^.x,Polygon^.First^.y,PolygonVertex^.x,PolygonVertex^.y,NextPolygonVertex^.x,NextPolygonVertex^.y);
     PolygonVertex:=NextPolygonVertex;
    until PolygonVertex^.Next=Polygon^.First;
   end else begin
    Vertices:=nil;
    try
     SetLength(Vertices,Polygon^.Count);
     PolygonVertex:=Polygon^.First;
     for i:=0 to length(Vertices)-1 do begin
      Vertices[i].x:=PolygonVertex^.x;
      Vertices[i].y:=PolygonVertex^.y;
      PolygonVertex:=PolygonVertex^.Next;
     end;
     for i:=1 to length(Vertices)-1 do begin
      Vertices[i].ArrayPrevious:=@Vertices[i-1];
     end;
     for i:=0 to length(Vertices)-2 do begin
      Vertices[i].ArrayNext:=@Vertices[i+1];
     end;
     Vertices[length(Vertices)-1].ArrayNext:=@Vertices[0];
     Vertices[0].ArrayPrevious:=@Vertices[length(Vertices)-1];
     Start:=@Vertices[0];
     Vertex:=Start;
     Start:=Vertex;
     Current:=Vertex;
     a:=0.0;
     repeat
      Next:=Current^.ArrayNext;
      a:=a+((Current^.x-Next^.x)*(Current^.y+Next^.y));
      Current:=Next;
     until Current=Start;
     if a<0.0 then begin
      repeat
       Vertex^.Previous:=Vertex^.ArrayNext;
       Vertex^.Next:=Vertex^.ArrayPrevious;
       Vertex:=Vertex^.ArrayPrevious;
      until Vertex=Start;
     end else begin
      repeat
       Vertex^.Previous:=Vertex^.ArrayPrevious;
       Vertex^.Next:=Vertex^.ArrayNext;
       Vertex:=Vertex^.ArrayNext;
      until Vertex=Start;
     end;
     Vertex:=Start;
     while Vertex^.Previous<>Vertex^.Next do begin
      px:=Vertex^.Previous^.x;
      py:=Vertex^.Previous^.y;
      vx:=Vertex^.x;
      vy:=Vertex^.y;
      nx:=Vertex^.Next^.x;
      ny:=Vertex^.Next^.y;
      vnx:=nx-vx;
      vny:=ny-vy;
      npx:=px-nx;
      npy:=py-ny;
      pvx:=vx-px;
      pvy:=vy-py;
      IsEar:=((pvx*vny)-(vnx*pvy))>EPSILON;
      if IsEar then begin
       Current:=Vertex^.Next;
       Stop:=Vertex^.Previous^.Previous;
       while Current<>Stop do begin
        Current:=Current^.Next;
        cx:=Current^.x;
        cy:=Current^.y;
        pcx:=cx-px;
        pcy:=cy-py;
        if (pcx<>0.0) or (pcy<>0.0) then begin
         vcx:=cx-vx;
         vcy:=cy-vy;
         if (vcx<>0.0) or (vcy<>0.0) then begin
          ncx:=cx-nx;
          ncy:=cy-ny;
          if (ncx<>0.0) or (ncy<>0.0) then begin
           if (((vnx*vcy)-(vny*vcx))>=(-EPSILON)) and (((pvx*pcy)-(pvy*pcx))>=(-EPSILON)) and (((npx*ncy)-(npy*ncx))>=(-EPSILON)) then begin
            IsEar:=false;
            break;
           end;
          end;
         end;
        end;
       end;
      end;
      if IsEar then begin
       AddTriangle(Vertex^.Previous^.x,Vertex^.Previous^.y,Vertex^.x,Vertex^.y,Vertex^.Next^.x,Vertex^.Next^.y);
       Vertex^.Previous^.Next:=Vertex^.Next;
       Vertex^.Next^.Previous:=Vertex^.Previous;
       Vertex:=Vertex^.Next;
       Start:=Vertex;
      end else begin
       Vertex:=Vertex^.Next;
       if Vertex=Start then begin
        SetLength(Vertices,0);
        exit;
       end;
      end;
     end;
    finally
     SetLength(Vertices,0);
    end;
   end;
  end;
 end;
 procedure PolygonsTriangulate(const Polygons:PPolygons);
 var Polygon:PPolygon;
 begin
  if assigned(Polygons) then begin
   Polygon:=Polygons^.First;
   while assigned(Polygon) do begin
    PolygonTriangulate(Polygon);
    Polygon:=Polygon^.Next;
   end;
  end;
 end;
 function IsSimple(const Polygon:TBeRoTriangulationPolygon):boolean;
 var a1i,b1i,k:longint;
     a0p,a1p,b0p,b1p:PBeRoTriangulationPoint;
     a10x,a10y,b10x,b10y,ab0x,ab0y,d,ai,bi,xi,yi:double;
 begin
  k:=length(Polygon);
  result:=k>2;
  if result then begin
   a0p:=@Polygon[k-1];
   for a1i:=0 to k-1 do begin
    a1p:=@Polygon[a1i];
    if SameValue(a0p^.x,a1p^.x) and SameValue(a0p^.y,a1p^.y) then begin
     result:=false;
     exit;
    end;
    a0p:=a1p;
   end;
   a0p:=@Polygon[k-1];
   for a1i:=0 to k-1 do begin
    a1p:=@Polygon[a1i];
    b0p:=a1p;
    for b1i:=a1i+1 to k-1 do begin
     b1p:=@Polygon[b1i];
     if a1i<>b1i then begin
      a10x:=a1p^.x-a0p^.x;
      a10y:=a1p^.y-a0p^.y;
      b10x:=b1p^.x-b0p^.x;
      b10y:=b1p^.y-b0p^.y;
      d:=(a10x*b10y)-(b10x*a10y);
      if not IsZero(d) then begin
       ab0x:=a0p^.x-b0p^.x;
       ab0y:=a0p^.y-b0p^.y;
       ai:=((b10x*ab0y)-(b10y*ab0x))/d;
       if (ai>=0.0) and (ai<=1.0) then begin
        bi:=((a10x*ab0y)-(a10y*ab0x))/d;
        if (bi>=0.0) and (bi<=1.0) then begin
         d:=1.0-ai;
         xi:=(a0p^.x*d)+(a1p^.x*ai);
         yi:=(a0p^.y*d)+(a1p^.y*ai);
         if not ((SameValue(a0p^.x,xi) and SameValue(a0p^.y,yi)) or
                 (SameValue(a1p^.x,xi) and SameValue(a1p^.y,yi)) or
                 (SameValue(b0p^.x,xi) and SameValue(b0p^.y,yi)) or
                 (SameValue(b1p^.x,xi) and SameValue(b1p^.y,yi))) then begin
          result:=false;
          exit;
         end;
        end;
       end;
      end;
     end;
     b0p:=b1p;
    end;
    a0p:=a1p;
   end;
  end;
 end;
var i,j,k:longint;
    Polygons:PPolygons;
    Polygon:PPolygon;
begin
 result:=false;
 GarbageCollector.First:=nil;
 GarbageCollector.Last:=nil;
 try
  k:=0;
  for i:=0 to length(InputPolygons)-1 do begin
   j:=length(InputPolygons[i]);
   if j<3 then begin
    exit;
   end;
   k:=k+j;
  end;
  if k>2 then begin
   CountOutputTriangles:=0;
   if (length(InputPolygons)=1) and (length(InputPolygons[0])=3) then begin
    OutputTriangle(InputPolygons[0,0].x,InputPolygons[0,0].y,InputPolygons[0,1].x,InputPolygons[0,1].y,InputPolygons[0,2].x,InputPolygons[0,2].y);
   end else begin
    if Quality<0 then begin
     Quality:=0;
    end else if Quality>4 then begin
     Quality:=4;
    end;
    DoOptimizeTriangles:=Quality in [2,4];
    GetMem(PointHashTable,SizeOf(TPointHashTable));
    try
     FillChar(PointHashTable^,SizeOf(TPointHashTable),0);
     if DoOptimizeTriangles then begin
      FillChar(Triangles,SizeOf(TTriangles),0);
      TriangleDynamicAABBTree:=DynamicAABBTreeCreate;
     end;
     try
      if (Quality>3) and ((length(InputPolygons)=1) and IsSimple(InputPolygons[0])) then begin
       Polygon:=PolygonCreate;
       for i:=0 to length(InputPolygons[0])-1 do begin
        PolygonAddVertex(Polygon,PolygonVertexCreate(InputPolygons[0,i].x,InputPolygons[0,i].y));
       end;
       PolygonTriangulate(Polygon);
       PolygonDestroy(Polygon);
      end else begin
       FillChar(Edges,SizeOf(TEdges),0);
       FillChar(Quads,SizeOf(TQuads),0);
       FillChar(Triangles,SizeOf(TTriangles),0);
       VerticalAATreeInit(@VerticalAATreeNodeRoot);
       try
        QuadDynamicAABBTree:=DynamicAABBTreeCreate;
        try
         AddEdges;
         MergeSortEdges;
         SplitEdgesAtIntersections;
         SplitEdgesAtYCoordinates;
         SweepVertical;
         case Quality of
          1,2:begin
           OptimizeQuads;
           TriangulateQuads;
          end;
          3,4:begin
           OptimizeQuads;
           Polygons:=MergeQuadsToPolygons;
           if assigned(Polygons) then begin
            PolygonsTriangulate(Polygons);
           end else begin
            TriangulateQuads;
           end;
          end;
          else {0:}begin
           TriangulateQuads;
          end;
         end;
        finally
         DynamicAABBTreeDestroy(QuadDynamicAABBTree);
        end;
       finally
        VerticalAATreeDone(@VerticalAATreeNodeRoot);
       end;
      end;
      if DoOptimizeTriangles then begin
       OptimizeTriangles;
       PushOutputTriangles;
      end;
     finally
      if DoOptimizeTriangles then begin
       DynamicAABBTreeDestroy(TriangleDynamicAABBTree);
      end;
     end;
    finally
     FreeMem(PointHashTable);
    end;
   end;
   SetLength(OutputTriangles,CountOutputTriangles);
   result:=CountOutputTriangles>0;
  end;
 finally
  GarbageCollectorCleanUp;
 end;
end;

{$ifndef STANDALONE}
{$ifndef FASTONLY}
function TriangulateDelaunayClipping(const InputPolygons:TBeRoTriangulationPolygons;var OutputTriangles:TBeRoTriangulationTriangles;InputPolygonFillRule:TBeRoTriangulationPolygonFillRule=btpfrEVENODD):boolean;
const ScaleFactor=16777216.0;
type PDelaunay2DVector=^TDelaunay2DVector;
     TDelaunay2DVector=record
      x,y:double;
     end;
     TDelaunay2DVertices=array of TDelaunay2DVector;
     TDelaunay2DTriangleIndices=array[0..3] of longint;
     PDelaunay2DTriangle=^TDelaunay2DTriangle;
     TDelaunay2DTriangle=TDelaunay2DTriangleIndices;
     TDelaunay2DTriangles=array of TDelaunay2DTriangle;
     TPDelaunay2DTriangles=array of PDelaunay2DTriangle;
var CountOutputTriangles:longint;
    Paths:TPaths;
 procedure PushTriangle(const ax,ay,bx,by,cx,cy:double);
 var Triangle:PBeRoTriangulationTriangle;
 begin
  if length(OutputTriangles)<=CountOutputTriangles then begin
   SetLength(OutputTriangles,(CountOutputTriangles+1)*2);
  end;
  Triangle:=@OutputTriangles[CountOutputTriangles];
  inc(CountOutputTriangles);
  if (((bx-ax)*(by+ay))+((cx-bx)*(cy+by))+((ax-cx)*(ay+cy)))<0.0 then begin
   Triangle^[0].x:=ax;
   Triangle^[0].y:=ay;
   Triangle^[1].x:=bx;
   Triangle^[1].y:=by;
   Triangle^[2].x:=cx;
   Triangle^[2].y:=cy;
  end else begin
   Triangle^[0].x:=cx;
   Triangle^[0].y:=cy;
   Triangle^[1].x:=bx;
   Triangle^[1].y:=by;
   Triangle^[2].x:=ax;
   Triangle^[2].y:=ay;
  end;
 end;
 procedure PushPathByEarClipping(const Path:TPath);
 const EPSILON=1e-14;
 type PVertex=^TVertex;
      TVertex=record
       ArrayPrevious:PVertex;
       ArrayNext:PVertex;
       Previous:PVertex;
       Next:PVertex;
       x,y:double;
      end;
      TVertices=array of TVertex;
 var Vertices:TVertices;
     i:longint;
     Vertex,Current,Start,Stop,Next:PVertex;
     a,px,py,vx,vy,nx,ny,vnx,vny,npx,npy,pvx,pvy,cx,cy,pcx,pcy,vcx,vcy,ncx,ncy:double;
     IsEar:boolean;
 begin
  if length(Path)=3 then begin
   PushTriangle(Path[0].x/ScaleFactor,Path[0].y/ScaleFactor,
                Path[1].x/ScaleFactor,Path[1].y/ScaleFactor,
                Path[2].x/ScaleFactor,Path[2].y/ScaleFactor);
  end else if length(Path)>3 then begin
   Vertices:=nil;
   try
    SetLength(Vertices,length(Path));
    for i:=0 to length(Vertices)-1 do begin
     Vertices[i].x:=Path[i].x/ScaleFactor;
     Vertices[i].y:=Path[i].y/ScaleFactor;
    end;
    for i:=1 to length(Vertices)-1 do begin
     Vertices[i].ArrayPrevious:=@Vertices[i-1];
    end;
    for i:=0 to length(Vertices)-2 do begin
     Vertices[i].ArrayNext:=@Vertices[i+1];
    end;
    Vertices[length(Vertices)-1].ArrayNext:=@Vertices[0];
    Vertices[0].ArrayPrevious:=@Vertices[length(Vertices)-1];
    Start:=@Vertices[0];
    Vertex:=Start;
    Start:=Vertex;
    Current:=Vertex;
    a:=0.0;
    repeat
     Next:=Current^.ArrayNext;
     a:=a+((Current^.x-Next^.x)*(Current^.y+Next^.y));
     Current:=Next;
    until Current=Start;
    if a<0.0 then begin
     repeat
      Vertex^.Previous:=Vertex^.ArrayNext;
      Vertex^.Next:=Vertex^.ArrayPrevious;
      Vertex:=Vertex^.ArrayPrevious;
     until Vertex=Start;
    end else begin
     repeat
      Vertex^.Previous:=Vertex^.ArrayPrevious;
      Vertex^.Next:=Vertex^.ArrayNext;
      Vertex:=Vertex^.ArrayNext;
     until Vertex=Start;
    end;
    Vertex:=Start;
    while Vertex^.Previous<>Vertex^.Next do begin
     px:=Vertex^.Previous^.x;
     py:=Vertex^.Previous^.y;
     vx:=Vertex^.x;
     vy:=Vertex^.y;
     nx:=Vertex^.Next^.x;
     ny:=Vertex^.Next^.y;
     vnx:=nx-vx;
     vny:=ny-vy;
     npx:=px-nx;
     npy:=py-ny;
     pvx:=vx-px;
     pvy:=vy-py;
     IsEar:=((pvx*vny)-(vnx*pvy))>EPSILON;
     if IsEar then begin
      Current:=Vertex^.Next;
      Stop:=Vertex^.Previous^.Previous;
      while Current<>Stop do begin
       Current:=Current^.Next;
       cx:=Current^.x;
       cy:=Current^.y;
       pcx:=cx-px;
       pcy:=cy-py;
       if (pcx<>0.0) or (pcy<>0.0) then begin
        vcx:=cx-vx;
        vcy:=cy-vy;
        if (vcx<>0.0) or (vcy<>0.0) then begin
         ncx:=cx-nx;
         ncy:=cy-ny;
         if (ncx<>0.0) or (ncy<>0.0) then begin
          if (((vnx*vcy)-(vny*vcx))>=(-EPSILON)) and (((pvx*pcy)-(pvy*pcx))>=(-EPSILON)) and (((npx*ncy)-(npy*ncx))>=(-EPSILON)) then begin
           IsEar:=false;
           break;
          end;
         end;
        end;
       end;
      end;
     end;
     if IsEar then begin
      PushTriangle(Vertex^.Previous^.x,Vertex^.Previous^.y,Vertex^.x,Vertex^.y,Vertex^.Next^.x,Vertex^.Next^.y);
      Vertex^.Previous^.Next:=Vertex^.Next;
      Vertex^.Next^.Previous:=Vertex^.Previous;
      Vertex:=Vertex^.Next;
      Start:=Vertex;
     end else begin
      Vertex:=Vertex^.Next;
      if Vertex=Start then begin
       SetLength(Vertices,0);
       exit;
      end;
     end;
    end;
   finally
    SetLength(Vertices,0);
   end;
  end;
 end;
 function DelaunayTriangulationWatson(const Points:TDelaunay2DVertices;var Triangles:TDelaunay2DTriangles):boolean;
 // The Delaunay triangulation algorithm in this sub function is based on nnsort.c by David F. Watson
 label Corner1,Corner2,Corner3;
 const Dim=2;
       Dim1=Dim+1;
       BIGNUM=1e+37;
       EPSILON=0.00001;
       TSIZE=75;
       RANGE=10.0;
       dm=dim-1;
 var Vertices:array of array of longint;
     nrs,i0,i1,i2,i3,i4,i5,i6,i7,i8,i9,i11,nts,tsz,i,ntris,a,b,c,d,itri:longint;
     xx,yy,bgs:double;
     ii:array[0..2] of longint;
     mxy:array[0..2,0..Dim-1] of double;
     wrk:array[0..Dim-1,0..Dim1-1] of double;
     pts:array of array[0..Dim-1] of double;
     tmp:array of array[0..Dim-1] of longint;
     id:array of longint;
     a3s:array of array[0..Dim1-1] of longint;
     ccr:array of array[0..Dim1-1] of double;
     nverts:array of longint;
 begin
  Vertices:=nil;
  pts:=nil;
  tmp:=nil;
  id:=nil;
  a3s:=nil;
  ccr:=nil;
  nverts:=nil;
  try
   nrs:=length(Points);
   for i0:=0 to Dim-1 do begin
    mxy[0,i0]:=-BIGNUM;
    mxy[1,i0]:=BIGNUM;
   end;
   for i0:=0 to Dim-1 do begin
    for i1:=0 to Dim1-1 do begin
     wrk[i0,i1]:=-RANGE;
    end;
   end;
   for i0:=0 to Dim-1 do begin
    wrk[i0,i0]:=RANGE*((3*Dim)-1);
   end;
   SetLength(pts,nrs+Dim1);
   for i0:=0 to nrs-1 do begin
    pts[i0,0]:=Points[i0].x;
    pts[i0,1]:=Points[i0].y;
    for i1:=0 to Dim-1 do begin
     if mxy[0,i1]<pts[i0,i1] then begin
      mxy[0,i1]:=pts[i0,i1];
     end;
     if mxy[1,i1]>pts[i0,i1] then begin
      mxy[1,i1]:=pts[i0,i1];
     end;
    end;
   end;
   bgs:=0;
   for i0:=0 to Dim-1 do begin
    mxy[0,i0]:=mxy[0,i0]-mxy[1,i0];
    if bgs<mxy[0,i0] then begin
     bgs:=mxy[0,i0];
    end;
   end;
   bgs:=bgs*EPSILON;
   RandSeed:=357;
   for i0:=0 to nrs-1 do begin
    for i1:=0 to Dim-1 do begin
     pts[i0,i1]:=pts[i0,i1]+(bgs*(0.5-random));
    end;
   end;
   for i0:=0 to Dim1-1 do begin
    for i1:=0 to Dim-1 do begin
     pts[nrs+i0,i1]:=mxy[1,i1]+(wrk[i1,i0]*mxy[0,i1]);
    end;
   end;
   i1:=1;
   for i0:=2 to dim1-1 do begin
    i1:=i1*i0;
   end;
   tsz:=TSIZE*i1;
   SetLength(tmp,tsz+1);
   i1:=i1*(nrs+(50*i1));
   SetLength(id,i1);
   for i0:=0 to i1-1 do begin
    id[i0]:=i0;
   end;
   SetLength(a3s,i1);
   SetLength(ccr,i1);
   a3s[0,0]:=nrs;
   for i0:=1 to dim1-1 do begin
    a3s[0,i0]:=a3s[0,i0-1]+1;
   end;
   ccr[0,Dim]:=BIGNUM;
   for i0:=0 to Dim-1 do begin
    ccr[0,i0]:=0;
   end;
   nts:=1;
   i4:=1;
   for i0:=0 to nrs-1 do begin
    i1:=-1;
    i7:=-1;
    i9:=0;
    for i11:=0 to nts-1 do begin
     inc(i1);
     while a3s[i1,0]<0 do begin
      inc(i1);
     end;
     xx:=ccr[i1,Dim];
     for i2:=0 to Dim-1 do begin
      xx:=xx-sqr(pts[i0,i2]-ccr[i1,i2]);
      if xx<0.0 then begin
       goto Corner3;
      end;
     end;
     dec(i9);
     dec(i4);
     id[i4]:=i1;
     for i2:=0 to Dim1-1 do begin
      ii[0]:=0;
      if ii[0]=i2 then begin
       inc(ii[0]);
      end;
      for i3:=1 to Dim-1 do begin
       ii[i3]:=ii[i3-1]+1;
       if ii[i3]=i2 then begin
        inc(ii[i3]);
       end;
      end;
      if i7>dm then begin
       i8:=i7;
       for i3:=0 to i8 do begin
        for i5:=0 to dim-1 do begin
         if a3s[i1,ii[i5]]<>tmp[i3,i5] then begin
          goto Corner1;
         end;
        end;
        for i6:=0 to Dim-1 do begin
         tmp[i3,i6]:=tmp[i8,i6];
        end;
        dec(i7);
        goto Corner2;
        Corner1:
       end;
      end;
      inc(i7);
      if (i7+1)>tsz then begin
       tsz:=(i7+1)*2;
       SetLength(tmp,tsz);
      end;
      for i3:=0 to Dim-1 do begin
       tmp[i7,i3]:=a3s[i1,ii[i3]];
      end;
      Corner2:
     end;
     a3s[i1,0]:=-1;
     Corner3:
    end;
    for i1:=0 to i7 do begin
     for i2:=0 to Dim-1 do begin
      wrk[i2,Dim]:=0;
      for i3:=0 to Dim-1 do begin
       wrk[i2,i3]:=pts[tmp[i1,i2],i3]-pts[i0,i3];
       wrk[i2,dim]:=wrk[i2,dim]+(wrk[i2,i3]*(pts[tmp[i1,i2],i3]+pts[i0,i3])*0.5);
      end;
     end;
     xx:=(wrk[0,0]*wrk[1,1])-(wrk[1,0]*wrk[0,1]);
     ccr[id[i4],0]:=((wrk[0,2]*wrk[1,1])-(wrk[1,2]*wrk[0,1]))/xx;
     ccr[id[i4],1]:=((wrk[0,0]*wrk[1,2])-(wrk[1,0]*wrk[0,2]))/xx;
     ccr[id[i4],Dim]:=0;
     for i2:=0 to Dim-1 do begin
      ccr[id[i4],dim]:=ccr[id[i4],dim]+sqr(pts[i0,i2]-ccr[id[i4],i2]);
      a3s[id[i4],i2]:=tmp[i1,i2];
     end;
     a3s[id[i4],dim]:=i0;
     inc(i4);
     inc(i9);
    end;
    inc(nts,i9);
   end;
   SetLength(nverts,nrs);
   for i:=0 to nrs-1 do begin
    nverts[i]:=0;
   end;
   ntris:=0;
   i0:=-1;
   for i11:=0 to nts-1 do begin
    inc(i0);
    while a3s[i0,0]<0 do begin
     inc(i0);
    end;
    if a3s[i0,0]<nrs then begin
     inc(ntris);
     inc(nverts[a3s[i0,0]]);
     inc(nverts[a3s[i0,1]]);
     inc(nverts[a3s[i0,2]]);
    end;
   end;
   SetLength(Vertices,nrs);
   for i:=0 to nrs-1 do begin
    SetLength(Vertices[i],nverts[i]);
    nverts[i]:=0;
   end;
   SetLength(Triangles,ntris);
   itri:=0;
   i0:=-1;
   for i11:=0 to nts-1 do begin
    inc(i0);
    while a3s[i0,0]<0 do begin
     inc(i0);
    end;
    if a3s[i0,0]<nrs then begin
     a:=a3s[i0,0];
     b:=a3s[i0,1];
     c:=a3s[i0,2];
     Vertices[a,nverts[a]]:=itri;
     Vertices[b,nverts[b]]:=itri;
     Vertices[c,nverts[c]]:=itri;
     inc(nverts[a]);
     inc(nverts[b]);
     inc(nverts[c]);
     Triangles[itri,0]:=a;
     Triangles[itri,1]:=b;
     Triangles[itri,2]:=c;
     inc(itri);
    end;
   end;
   result:=length(Triangles)>0;
  finally
   SetLength(nverts,0);
   SetLength(ccr,0);
   SetLength(a3s,0);
   SetLength(id,0);
   SetLength(tmp,0);
   SetLength(pts,0);
   SetLength(Vertices,0);
  end;
 end;
 procedure PushPathByRecursiveDelaunayAndClipping(const Path:TPath;Depth:longint);
 var Delaunay2DVertices:TDelaunay2DVertices;
     Delaunay2DTriangles:TDelaunay2DTriangles;
     Clipper:TClipper;
     i,j:longint;
     TempPaths:TPaths;
 begin
  if length(Path)=3 then begin
   PushTriangle(Path[0].x/ScaleFactor,Path[0].y/ScaleFactor,
                Path[1].x/ScaleFactor,Path[1].y/ScaleFactor,
                Path[2].x/ScaleFactor,Path[2].y/ScaleFactor);
  end else if length(Path)>3 then begin
   if Depth>0 then begin
    Delaunay2DVertices:=nil;
    Delaunay2DTriangles:=nil;
    TempPaths:=nil;
    try
     SetLength(Delaunay2DVertices,length(Path));
     for i:=0 to length(Path)-1 do begin
      Delaunay2DVertices[i].x:=Path[i].x;
      Delaunay2DVertices[i].y:=Path[i].y;
     end;
     if DelaunayTriangulationWatson(Delaunay2DVertices,Delaunay2DTriangles) then begin
      Clipper:=TClipper.Create;
      try
       for i:=0 to length(Delaunay2DTriangles)-1 do begin
        SetLength(TempPaths,1,3);
        SetLength(TempPaths[0],3);
        TempPaths[0,0].x:=round(Delaunay2DVertices[Delaunay2DTriangles[i,0]].x);
        TempPaths[0,0].y:=round(Delaunay2DVertices[Delaunay2DTriangles[i,0]].y);
        TempPaths[0,1].x:=round(Delaunay2DVertices[Delaunay2DTriangles[i,1]].x);
        TempPaths[0,1].y:=round(Delaunay2DVertices[Delaunay2DTriangles[i,1]].y);
        TempPaths[0,2].x:=round(Delaunay2DVertices[Delaunay2DTriangles[i,2]].x);
        TempPaths[0,2].y:=round(Delaunay2DVertices[Delaunay2DTriangles[i,2]].y);
        Clipper.Clear;
        Clipper.StrictlySimple:=true;
        Clipper.AddPaths(TempPaths,ptSubject,true);
        Clipper.AddPath(Path,ptClip,true);
        Clipper.Execute(ctIntersection,TempPaths,pftEvenOdd,pftEvenOdd);
        for j:=0 to length(TempPaths)-1 do begin
         PushPathByRecursiveDelaunayAndClipping(TempPaths[j],Depth-1);
        end;
       end;
      finally
       Clipper.Free;
      end;
     end else begin
      PushPathByEarClipping(Path);
     end;
    finally
     SetLength(Delaunay2DVertices,0);
     SetLength(Delaunay2DTriangles,0);
     SetLength(TempPaths,0);
    end;
   end else begin
    PushPathByEarClipping(Path);
   end;
  end;
 end;
var i,j,k:longint;
    Delaunay2DVertices:TDelaunay2DVertices;
    Delaunay2DTriangles:TDelaunay2DTriangles;
    TempPaths:TPaths;
    Clipper:TClipper;
    PolyFillType:TPolyFillType;
begin
 result:=false;
 OutputTriangles:=nil;
 Delaunay2DVertices:=nil;
 Delaunay2DTriangles:=nil;
 Paths:=nil;
 TempPaths:=nil;
 try
  k:=0;
  for i:=0 to length(InputPolygons)-1 do begin
   j:=length(InputPolygons[i]);
   if j<3 then begin
    exit;
   end;
   k:=k+j;
  end;
  if k>2 then begin
   CountOutputTriangles:=0;
   if (length(InputPolygons)=1) and (length(InputPolygons[0])=3) then begin
    PushTriangle(InputPolygons[0,0].x,InputPolygons[0,0].y,InputPolygons[0,1].x,InputPolygons[0,1].y,InputPolygons[0,2].x,InputPolygons[0,2].y);
   end else begin
    SetLength(Delaunay2DVertices,k);
    k:=0;
    SetLength(Paths,length(InputPolygons));
    for i:=0 to length(InputPolygons)-1 do begin
     SetLength(Paths[i],length(InputPolygons[i]));
     for j:=0 to length(InputPolygons[i])-1 do begin
      Paths[i,j].x:=round(InputPolygons[i,j].x*ScaleFactor);
      Paths[i,j].y:=round(InputPolygons[i,j].y*ScaleFactor);
      Delaunay2DVertices[k].x:=Paths[i,j].x;
      Delaunay2DVertices[k].y:=Paths[i,j].y;
      inc(k);
     end;
    end;
    if DelaunayTriangulationWatson(Delaunay2DVertices,Delaunay2DTriangles) then begin
     Clipper:=TClipper.Create;
     try
      for i:=0 to length(Delaunay2DTriangles)-1 do begin
       SetLength(TempPaths,1,3);
       TempPaths[0,0].x:=round(Delaunay2DVertices[Delaunay2DTriangles[i,0]].x);
       TempPaths[0,0].y:=round(Delaunay2DVertices[Delaunay2DTriangles[i,0]].y);
       TempPaths[0,1].x:=round(Delaunay2DVertices[Delaunay2DTriangles[i,1]].x);
       TempPaths[0,1].y:=round(Delaunay2DVertices[Delaunay2DTriangles[i,1]].y);
       TempPaths[0,2].x:=round(Delaunay2DVertices[Delaunay2DTriangles[i,2]].x);
       TempPaths[0,2].y:=round(Delaunay2DVertices[Delaunay2DTriangles[i,2]].y);
       Clipper.Clear;
       Clipper.StrictlySimple:=true;
       Clipper.AddPaths(TempPaths,ptSubject,true);
       Clipper.AddPaths(Paths,ptClip,true);
       case InputPolygonFillRule of
        btpfrEVENODD:begin
         PolyFillType:=pftEvenOdd;
        end;
        btpfrNONZERO:begin
         PolyFillType:=pftNonZero;
        end;
        btpfrPOSITIVE:begin
         PolyFillType:=pftPositive;
        end;
        btpfrNEGATIVE:begin
         PolyFillType:=pftNegative;
        end;
        else begin
         PolyFillType:=pftEvenOdd;
        end;
       end;
       Clipper.Execute(ctIntersection,TempPaths,pftEvenOdd,PolyFillType);
       for j:=0 to length(TempPaths)-1 do begin
{$ifdef FASTER}
        // For faster results, sacrificing triangulation quality for raw speed
        PushPathByEarClipping(TempPaths[j]);
{$else}
        // For more optimal results
        PushPathByRecursiveDelaunayAndClipping(TempPaths[j],3);
{$endif}
       end;
      end;
     finally
      Clipper.Free;
     end;
    end;
   end;
   SetLength(OutputTriangles,CountOutputTriangles);
   result:=CountOutputTriangles>0;
  end;
 finally
  SetLength(Delaunay2DVertices,0);
  SetLength(Delaunay2DTriangles,0);
  SetLength(Paths,0);
  SetLength(TempPaths,0);
 end;
end;
{$endif}

// Faster non-robust-variant with seidel algorithm for realtime usages
function TriangulateSeidel(const InputPolygons:TBeRoTriangulationPolygons;var OutputTriangles:TBeRoTriangulationTriangles):boolean;
type PGarbageCollectedHeader=^TGarbageCollectedHeader;
     TGarbageCollectedHeader=record
      Next:PGarbageCollectedHeader;
     end;
     PPoint=^TPoint;
     PEdge=^TEdge;
     PPolygon=^TPolygon;
     PPolygonNode=^TPolygonNode;
     PNode=^TNode;
     PTrapezoid=^TTrapezoid;
     PQueryGraph=^TQueryGraph;
     PTrapezoidalMap=^TTrapezoidalMap;
     TPoint=record
      GarbageCollectedHeader:TGarbageCollectedHeader;
      x:double;
      y:double;
     end;
     TEdge=record
      GarbageCollectedHeader:TGarbageCollectedHeader;
      p:PPoint;
      q:PPoint;
      Slope:double;
      Polygon:PPolygon;
      Below:PTrapezoid;
      Above:PTrapezoid;
     end;
     TPolygon=record
      GarbageCollectedHeader:TGarbageCollectedHeader;
      Length:longint;
      First:PPolygonNode;
      Last:PPolygonNode;
     end;
     TPolygonNode=record
      GarbageCollectedHeader:TGarbageCollectedHeader;
      Point:PPoint;
      Next:PPolygonNode;
      Previous:PPolygonNode;
      Ear:longbool;
     end;
     TNode=record
      GarbageCollectedHeader:TGarbageCollectedHeader;
      Point:PPoint;
      Edge:PEdge;
      Left:PNode;
      Right:PNode;
      Trapezoid:PTrapezoid;
     end;
     TTrapezoid=record
      GarbageCollectedHeader:TGarbageCollectedHeader;
      Previous:PTrapezoid;
      Next:PTrapezoid;
      LeftPoint:PPoint;
      RightPoint:PPoint;
      Top:PEdge;
      Bottom:PEdge;
      Inside:longbool;
      Removed:longbool;
      Sink:PNode;
      UpperLeft:PTrapezoid;
      UpperRight:PTrapezoid;
      LowerLeft:PTrapezoid;
      LowerRight:PTrapezoid;
     end;
     TQueryGraph=record
      GarbageCollectedHeader:TGarbageCollectedHeader;
      Head:PNode;
     end;
     TTrapezoidalMap=record
      GarbageCollectedHeader:TGarbageCollectedHeader;
      Root:PTrapezoid;
      First:PTrapezoid;
      Last:PTrapezoid;
      Count:longint;
      QueryGraph:PQueryGraph;
      bcross:PEdge;
      tcross:PEdge;
     end;
var GarbageCollectorList:pointer;
    CountOutputTriangles:longint;
    FalseTrapezoid:PTrapezoid;
 function GarbageCollectorAllocate(Size:longint):pointer;
 begin
  GetMem(result,Size);
  FillChar(result^,Size,0);
  PGarbageCollectedHeader(result)^.Next:=GarbageCollectorList;
  GarbageCollectorList:=result;
 end;
 procedure GarbageCollectorCleanUp;
 var Current,Next:pointer;
 begin
  Current:=GarbageCollectorList;
  GarbageCollectorList:=nil;
  while assigned(Current) do begin
   Next:=PGarbageCollectedHeader(Current)^.Next;
   FreeMem(Current);
   Current:=Next;
  end;
 end;
 function NewPoint(const x,y:double):PPoint;
 begin
  result:=GarbageCollectorAllocate(SizeOf(TPoint));
  result^.x:=x;
  result^.y:=y;
 end;
 function NewEdge(const p,q:PPoint):PEdge;
 begin
  result:=GarbageCollectorAllocate(SizeOf(TEdge));
  result^.p:=p;
  result^.q:=q;
  if IsInfinite(p^.x) or IsInfinite(p^.y) or IsInfinite(q^.x) or IsInfinite(q^.y) or IsZero(q^.x-p^.x) then begin
   result^.Slope:=Infinity;
  end else begin
   result^.Slope:=(q^.y-p^.y)/(q^.x-p^.x);
  end;
  result^.Polygon:=nil;
  result^.Below:=nil;
  result^.Above:=nil;
 end;
 function NewPolygonNode(const Point:PPoint):PPolygonNode;
 begin
  result:=GarbageCollectorAllocate(SizeOf(TPolygonNode));
  result^.Point:=Point;
  result^.Next:=nil;
  result^.Previous:=nil;
  result^.Ear:=false;
 end;
 function NewPolygon:PPolygon;
 begin
  result:=GarbageCollectorAllocate(SizeOf(TPolygon));
  result^.Length:=0;
  result^.First:=nil;
  result^.Last:=nil;
 end;
 procedure PolygonAdd(const Polygon:PPolygon;const Point:PPoint);
 var Node:PPolygonNode;
 begin
  Node:=NewPolygonNode(Point);
  if Polygon^.Length=0 then begin
   Polygon^.First:=Node;
   Polygon^.Last:=Node;
  end else begin
   Polygon^.Last^.Next:=Node;
   Node^.Previous:=Polygon^.Last;
   Polygon^.Last:=Node;
  end;
  inc(Polygon^.Length);
 end;
 procedure PolygonRemove(const Polygon:PPolygon;const Node:PPolygonNode);
 begin
  if Polygon^.Length>0 then begin
   if Polygon^.First=Node then begin
    Polygon^.First:=Polygon^.First^.Next;
    if assigned(Polygon^.First) then begin
     Polygon^.First^.Previous:=nil;
    end else begin
     Polygon^.Last:=nil;
    end;
   end else if Polygon^.Last=Node then begin
    Polygon^.Last:=Polygon^.Last^.Previous;
    Polygon^.Last^.Next:=nil;
   end else begin
    Node^.Previous^.Next:=Node^.Next;
    Node^.Next^.Previous:=Node^.Previous;
   end;
   Node^.Next:=nil;
   Node^.Previous:=nil;
   dec(Polygon^.Length);
  end;
 end;
 procedure PolygonInsertBefore(const Polygon:PPolygon;const Point:PPoint;const Node:PPolygonNode);
 var NewNode:PPolygonNode;
 begin
  NewNode:=NewPolygonNode(Point);
  NewNode^.Previous:=Node^.Previous;
  NewNode^.Next:=Node;
  if assigned(Node^.Previous) then begin
   Node^.Previous^.Next:=NewNode;
  end else begin
   Polygon^.First:=NewNode;
  end;
  Node^.Previous:=NewNode;
  inc(Polygon^.Length);
 end;
 function neg(const p1,p2:PPoint):boolean;
 begin
  result:=(p1^.x<>p2^.x) or (p1^.y<>p2^.y);
 end;
 function Clone(const p:PPoint):PPoint;
 begin
  result:=NewPoint(p^.x,p^.y);
 end;
 function Cross(const a,b,c:PPoint):double;
 begin
  result:=((a^.x-c^.x)*(b^.y-c^.y))-((a^.y-c^.y)*(b^.x-c^.x));
 end;
 function EdgeOrient(const Edge:PEdge;const Point:PPoint):double;
 begin
  result:=Cross(Edge^.p,Edge^.q,Point);
 end;
 function EdgeAbove(const Edge:PEdge;const Point:PPoint):boolean;
 begin
  result:=Cross(Edge^.p,Edge^.q,Point)<0.0;
 end;
 function EdgeBelow(const Edge:PEdge;const Point:PPoint):boolean;
 begin
  result:=Cross(Edge^.p,Edge^.q,Point)>0.0;
 end;
 procedure PushTriangle(const a,b,c:PPoint);
 var Triangle:PBeRoTriangulationTriangle;
 begin
  if length(OutputTriangles)<=CountOutputTriangles then begin
   SetLength(OutputTriangles,(CountOutputTriangles+1)*2);
  end;
  Triangle:=@OutputTriangles[CountOutputTriangles];
  inc(CountOutputTriangles);
  if (((b^.x-a^.x)*(b^.y+a^.y))+((c^.x-b^.x)*(c^.y+b^.y))+((a^.x-c^.x)*(a^.y+c^.y)))<0.0 then begin
   Triangle^[0].x:=a^.x;
   Triangle^[0].y:=a^.y;
   Triangle^[1].x:=b^.x;
   Triangle^[1].y:=b^.y;
   Triangle^[2].x:=c^.x;
   Triangle^[2].y:=c^.y;
  end else begin
   Triangle^[0].x:=c^.x;
   Triangle^[0].y:=c^.y;
   Triangle^[1].x:=b^.x;
   Triangle^[1].y:=b^.y;
   Triangle^[2].x:=a^.x;
   Triangle^[2].y:=a^.y;
  end;
 end;
 procedure TriangulateMountain(const Edge:PEdge);
 var ConvexPoints:array of PPolygonNode;
     CountConvexPoints,ConvexPointIndex:longint;
     a,b:PPoint;
     Polygon:PPolygon;
     p,Ear,Previous,Next:PPolygonNode;
     Positive:boolean;
  procedure AddEar(Node:PPolygonNode);
  begin
   if (not Node^.Ear) and (Node<>Polygon^.First) and (Node<>Polygon^.Last) and
       ((Cross(Node^.Next^.Point,Node^.Previous^.Point,Node^.Point)>0.0)=Positive) then begin
    Node^.Ear:=true;
    inc(CountConvexPoints);
    if length(ConvexPoints)<CountConvexPoints then begin
     SetLength(ConvexPoints,CountConvexPoints*2);
    end;
    ConvexPoints[CountConvexPoints-1]:=Node;
   end;
  end;
 begin
  ConvexPoints:=nil;
  CountConvexPoints:=0;
  try
   a:=Edge^.p;
   b:=Edge^.q;
   Polygon:=Edge^.Polygon;
   p:=Polygon^.First^.Next;
   if Polygon^.Length>2 then begin
    if Polygon^.Length=3 then begin
     PushTriangle(a,p^.Point,b);
    end else begin
     Positive:=Cross(p^.Point,b,a)>0.0;
     while p<>Polygon^.Last do begin
      AddEar(p);
      p:=p^.Next;
     end;
     ConvexPointIndex:=0;
     while ConvexPointIndex<CountConvexPoints do begin
      Ear:=ConvexPoints[ConvexPointIndex];
      inc(ConvexPointIndex);
      Previous:=Ear^.Previous;
      Next:=Ear^.Next;
      PushTriangle(Previous^.Point,Ear^.Point,Next^.Point);
      PolygonRemove(Polygon,Ear);
      AddEar(Previous);
      AddEar(Next);
     end;
    end;
   end;
  finally
   SetLength(ConvexPoints,0);
  end;
 end;
 function NewTrapezoid(const LeftPoint,RightPoint:PPoint;const Top,Bottom:PEdge):PTrapezoid;
 begin
  result:=GarbageCollectorAllocate(SizeOf(TTrapezoid));
  result^.Previous:=nil;
  result^.Next:=nil;
  result^.LeftPoint:=LeftPoint;
  result^.RightPoint:=RightPoint;
  result^.Top:=Top;
  result^.Bottom:=Bottom;
  result^.Inside:=false;
  result^.Removed:=false;
  result^.Sink:=nil;
  result^.UpperLeft:=nil;
  result^.UpperRight:=nil;
  result^.LowerLeft:=nil;
  result^.LowerRight:=nil;
 end;
 function NewNode:PNode;
 begin
  result:=GarbageCollectorAllocate(SizeOf(TNode));
  result^.Point:=nil;
  result^.Edge:=nil;
  result^.Left:=nil;
  result^.Right:=nil;
  result^.Trapezoid:=nil;
 end;
 function SetYNode(const Node:PNode;const Edge:PEdge;const Left,Right:PNode):PNode;
 begin
  result:=Node;
  result^.Edge:=Edge;
  result^.Left:=Left;
  result^.Right:=Right;
  result^.Trapezoid:=nil;
 end;
 function SetXNode(const Node:PNode;const Point:PPoint;const Left,Right:PNode):PNode;
 begin
  result:=Node;
  result^.Point:=Point;
  result^.Left:=Left;
  result^.Right:=Right;
  result^.Trapezoid:=nil;
 end;
 function SetSink(const Node:PNode;const Trapezoid:PTrapezoid):PNode;
 begin
  result:=Node;
  result^.Trapezoid:=Trapezoid;
  Trapezoid^.Sink:=Node;
 end;
 function GetSink(Trapezoid:PTrapezoid):PNode;
 begin
  if assigned(Trapezoid^.Sink) then begin
   result:=Trapezoid^.Sink;
  end else begin
   result:=SetSink(NewNode,Trapezoid);
  end;
 end;
 function NewQueryGraph(Head:PTrapezoid):PQueryGraph;
 begin
  result:=GarbageCollectorAllocate(SizeOf(TQueryGraph));
  result^.Head:=GetSink(Head);
 end;
 function QueryGraphLocate(const QueryGraph:PQueryGraph;const Point:PPoint;const Slope:double):PTrapezoid;
 var Node:PNode;
     Orient:double;
 begin
  result:=nil;
  Node:=QueryGraph^.Head;
  while assigned(Node) do begin
   if assigned(Node^.Trapezoid) then begin
    result:=Node^.Trapezoid;
    exit;
   end else begin
    if assigned(Node^.Point) then begin
     if Point^.x>=Node^.Point^.x then begin
      Node:=Node^.Right;
     end else begin
      Node:=Node^.Left;
     end;
    end else if assigned(Node^.Edge) then begin
     Orient:=EdgeOrient(Node^.Edge,Point);
     if Orient<0.0 then begin
      Node:=Node^.Right;
     end else if Orient>0.0 then begin
      Node:=Node^.Left;
     end else if Slope<Node^.Edge^.Slope then begin
      Node:=Node^.Right;
     end else begin
      Node:=Node^.Left;
     end;
    end else begin
     Assert(false);
     exit;
    end;
   end;
  end;
 end;
 procedure QueryGraphCase1(const QueryGraph:PQueryGraph;const Sink:PNode;const Edge:PEdge;const t1,t2,t3,t4:PTrapezoid);
 var yNode,qNode:PNode;
 begin
  yNode:=SetYNode(NewNode,Edge,GetSink(t2),GetSink(t3));
  qNode:=SetXNode(NewNode,Edge^.q,yNode,GetSink(t4));
  SetXNode(Sink,Edge^.p,GetSink(t1),qNode);
 end;
 procedure QueryGraphCase2(const QueryGraph:PQueryGraph;const Sink:PNode;const Edge:PEdge;const t1,t2,t3:PTrapezoid);
 var yNode:PNode;
 begin
  yNode:=SetYNode(NewNode,Edge,GetSink(t2),GetSink(t3));
  SetXNode(Sink,Edge^.p,GetSink(t1),yNode);
 end;
 procedure QueryGraphCase3(const QueryGraph:PQueryGraph;const Sink:PNode;const Edge:PEdge;const t1,t2:PTrapezoid);
 begin
  SetYNode(Sink,Edge,GetSink(t1),GetSink(t2));
 end;
 procedure QueryGraphCase4(const QueryGraph:PQueryGraph;const Sink:PNode;const Edge:PEdge;const t1,t2,t3:PTrapezoid);
 var yNode:PNode;
 begin
  yNode:=SetYNode(NewNode,Edge,GetSink(t1),GetSink(t2));
  SetXNode(Sink,Edge^.q,yNode,GetSink(t3));
 end;
 procedure TrapezoidUpdateLeft(const Trapezoid,ul,ll:PTrapezoid);
 begin
  Trapezoid^.UpperLeft:=ul;
  if assigned(ul) then begin
   ul^.UpperRight:=Trapezoid;
  end;
  Trapezoid^.LowerLeft:=ll;
  if assigned(ll) then begin
   ll^.LowerRight:=Trapezoid;
  end;
 end;
 procedure TrapezoidUpdateRight(const Trapezoid,ur,lr:PTrapezoid);
 begin
  Trapezoid^.UpperRight:=ur;
  if assigned(ur) then begin
   ur^.UpperLeft:=Trapezoid;
  end;
  Trapezoid^.LowerRight:=lr;
  if assigned(lr) then begin
   lr^.LowerLeft:=Trapezoid;
  end;
 end;
 procedure TrapezoidUpdateLeftRight(const Trapezoid,ul,ll,ur,lr:PTrapezoid);
 begin
  TrapezoidUpdateLeft(Trapezoid,ul,ll);
  TrapezoidUpdateRight(Trapezoid,ur,lr);
 end;
 procedure TrapezoidMarkInside(const Trapezoid:PTrapezoid);
 var Stack:array of PTrapezoid;
     StackPointer:longint;
     t:PTrapezoid;
 begin
  Stack:=nil;
  try
   SetLength(Stack,16);
   Stack[0]:=Trapezoid;
   StackPointer:=1;
   while StackPointer>0 do begin
    dec(StackPointer);
    t:=Stack[StackPointer];
    if not t^.Inside then begin
     t^.Inside:=true;
     if assigned(t^.UpperLeft) then begin
      if length(Stack)<=StackPointer then begin
       SetLength(Stack,(StackPointer+1)*2);
      end;
      Stack[StackPointer]:=t^.UpperLeft;
      inc(StackPointer);
     end;
     if assigned(t^.LowerLeft) then begin
      if length(Stack)<=StackPointer then begin
       SetLength(Stack,(StackPointer+1)*2);
      end;
      Stack[StackPointer]:=t^.LowerLeft;
      inc(StackPointer);
     end;
     if assigned(t^.UpperRight) then begin
      if length(Stack)<=StackPointer then begin
       SetLength(Stack,(StackPointer+1)*2);
      end;
      Stack[StackPointer]:=t^.UpperRight;
      inc(StackPointer);
     end;
     if assigned(t^.LowerRight) then begin
      if length(Stack)<=StackPointer then begin
       SetLength(Stack,(StackPointer+1)*2);
      end;
      Stack[StackPointer]:=t^.LowerRight;
      inc(StackPointer);
     end;
    end;
   end;
  finally
   SetLength(Stack,0);
  end;
 end;
 function TrapezoidContains(const Trapezoid:PTrapezoid;const Point:PPoint):boolean;
 begin
  result:=(Point^.x>Trapezoid^.LeftPoint^.x) and
          (Point^.x<Trapezoid^.RightPoint^.x) and
          EdgeAbove(Trapezoid^.Top,Point) and
          EdgeBelow(Trapezoid^.Bottom,Point);
 end;
 procedure TrapezoidAddPoint(const Trapezoid:PTrapezoid;const Edge:PEdge;const Point:PPoint);
 var Polygon:PPolygon;
     v:PPolygonNode;
 begin
  Polygon:=Edge^.Polygon;
  if assigned(Polygon) then begin
   v:=Polygon^.First;
   while assigned(v) do begin
    if not neg(Point,v^.Point) then begin
     exit;
    end;
    if Point^.x<v^.Point^.x then begin
     PolygonInsertBefore(Polygon,Point,v);
     exit;
    end;
    v:=v^.Next;
   end;
   PolygonAdd(Polygon,Point);
  end else begin
   if neg(Point,Edge^.p) and neg(Point,Edge^.q) then begin
    Polygon:=NewPolygon;
    Edge^.Polygon:=Polygon;
    PolygonAdd(Polygon,Edge^.p);
    PolygonAdd(Polygon,Point);
    PolygonAdd(Polygon,Edge^.q);
   end;
  end;
 end;
 procedure TrapezoidAddPoints(const Trapezoid:PTrapezoid);
 begin
  if Trapezoid^.LeftPoint<>Trapezoid^.Bottom^.p then begin
   TrapezoidAddPoint(Trapezoid,Trapezoid^.Bottom,Trapezoid^.LeftPoint);
  end;
  if Trapezoid^.RightPoint<>Trapezoid^.Bottom^.q then begin
   TrapezoidAddPoint(Trapezoid,Trapezoid^.Bottom,Trapezoid^.RightPoint);
  end;
  if Trapezoid^.LeftPoint<>Trapezoid^.Top^.p then begin
   TrapezoidAddPoint(Trapezoid,Trapezoid^.Top,Trapezoid^.LeftPoint);
  end;
  if Trapezoid^.RightPoint<>Trapezoid^.Top^.q then begin
   TrapezoidAddPoint(Trapezoid,Trapezoid^.Top,Trapezoid^.RightPoint);
  end;
 end;
 function NewTrapezoidalMap:PTrapezoidalMap;
 var Top,Bottom:PEdge;
 begin
  Top:=NewEdge(NewPoint(-Infinity,Infinity),NewPoint(Infinity,Infinity));
  Bottom:=NewEdge(NewPoint(-Infinity,-Infinity),NewPoint(Infinity,-Infinity));
  result:=GarbageCollectorAllocate(SizeOf(TTrapezoidalMap));
  result^.Root:=NewTrapezoid(Bottom.p,Top.q,Top,Bottom);
  result^.First:=result^.Root;
  result^.Last:=result^.Root;
  result^.Count:=1;
  result^.QueryGraph:=NewQueryGraph(result^.Root);
 end;
 procedure TrapezoidalMapPush(const TrapezoidalMap:PTrapezoidalMap;const Trapezoid:PTrapezoid);
 begin
  if assigned(TrapezoidalMap^.Last) then begin
   TrapezoidalMap^.Last^.Next:=Trapezoid;
   Trapezoid^.Previous:=TrapezoidalMap^.Last;
  end else begin
   TrapezoidalMap^.First:=Trapezoid;
   Trapezoid^.Previous:=nil;
  end;
  TrapezoidalMap^.Last:=Trapezoid;
  Trapezoid^.Next:=nil;
 end;
 function TrapezoidalMapNextTrapezoid(const TrapezoidalMap:PTrapezoidalMap;const t:PTrapezoid;const Edge:PEdge):PTrapezoid;
 begin
  if Edge^.q^.x<=t^.RightPoint^.x then begin
   result:=FalseTrapezoid;
  end else if EdgeAbove(Edge,t^.RightPoint) then begin
   result:=t^.UpperRight;
  end else begin
   result:=t^.LowerRight;
  end;
 end;
 function TrapezoidalMapCase1(const TrapezoidalMap:PTrapezoidalMap;const t:PTrapezoid;const e:PEdge):PTrapezoid;
 var t2,t3,t4:PTrapezoid;
     Sink:PNode;
 begin
  t2:=NewTrapezoid(e^.p,e^.q,t^.Top,e);
  t3:=NewTrapezoid(e^.p,e^.q,e,t^.Bottom);
  t4:=NewTrapezoid(e^.q,t^.RightPoint,t^.Top,t^.Bottom);
  TrapezoidUpdateRight(t4,t^.UpperRight,t^.LowerRight);
  TrapezoidUpdateLeft(t4,t2,t3);
  t^.RightPoint:=e^.p;
  TrapezoidUpdateRight(t,t2,t3);
  Sink:=t^.Sink;
  t^.Sink:=nil;
  QueryGraphCase1(TrapezoidalMap^.QueryGraph,Sink,e,t,t2,t3,t4);
  TrapezoidalMapPush(TrapezoidalMap,t2);
  TrapezoidalMapPush(TrapezoidalMap,t3);
  TrapezoidalMapPush(TrapezoidalMap,t4);
  result:=FalseTrapezoid;
 end;
 function TrapezoidalMapCase2(const TrapezoidalMap:PTrapezoidalMap;const t:PTrapezoid;const e:PEdge):PTrapezoid;
 var Next,t2,t3:PTrapezoid;
     Sink:PNode;
 begin
  Next:=TrapezoidalMapNextTrapezoid(TrapezoidalMap,t,e);
  t2:=NewTrapezoid(e^.p,t^.RightPoint,t^.Top,e);
  t3:=NewTrapezoid(e^.p,t^.RightPoint,e,t^.Bottom);
  t^.RightPoint:=e^.p;
  TrapezoidUpdateLeft(t,t^.UpperLeft,t^.LowerLeft);
  TrapezoidUpdateLeftRight(t2,t,nil,t^.UpperRight,nil);
  TrapezoidUpdateLeftRight(t3,nil,t,nil,t^.LowerRight);
  TrapezoidalMap^.bcross:=t^.Bottom;
  TrapezoidalMap^.tcross:=t^.Top;
  e^.Above:=t2;
  e^.Below:=t3;
  Sink:=t^.Sink;
  t^.Sink:=nil;
  QueryGraphCase2(TrapezoidalMap^.QueryGraph,Sink,e,t,t2,t3);
  TrapezoidalMapPush(TrapezoidalMap,t2);
  TrapezoidalMapPush(TrapezoidalMap,t3);
  result:=Next;
 end;
 function TrapezoidalMapCase3(const TrapezoidalMap:PTrapezoidalMap;const t:PTrapezoid;const e:PEdge):PTrapezoid;
 var Next,t1,t2,LowerRight,LowerLeft:PTrapezoid;
     Bottom,Top:PEdge;
     Sink:PNode;
 begin
  Next:=TrapezoidalMapNextTrapezoid(TrapezoidalMap,t,e);
  Bottom:=t^.Bottom;
  LowerRight:=t^.LowerRight;
  LowerLeft:=t^.LowerLeft;
  Top:=t^.Top;
  if TrapezoidalMap^.tcross=t^.Top then begin
   t1:=t^.UpperLeft;
   TrapezoidUpdateRight(t1,t^.UpperRight,nil);
   t1^.RightPoint:=t^.RightPoint;
  end else begin
   t1:=t;
   t1^.Bottom:=e;
   t1^.LowerLeft:=e^.Above;
   if assigned(e^.Above) then begin
    e^.Above^.LowerRight:=t1;
   end;
   t1^.LowerRight:=nil;
  end;
  if TrapezoidalMap^.bcross=Bottom then begin
   t2:=LowerLeft;
   TrapezoidUpdateRight(t2,nil,LowerRight);
   t2^.RightPoint:=t^.RightPoint;
  end else if t1=t then begin
   t2:=NewTrapezoid(t^.LeftPoint,t^.RightPoint,e,Bottom);
   TrapezoidUpdateLeftRight(t2,e^.Below,LowerLeft,nil,LowerRight);
   TrapezoidalMapPush(TrapezoidalMap,t2);
  end else begin
   t2:=t;
   t2^.Top:=e;
   t2^.UpperLeft:=e^.Below;
   if assigned(e^.Below) then begin
    e^.Below^.UpperRight:=t2;
   end;
   t2^.UpperRight:=nil;
  end;
  if (t<>t1) and (t<>t2) then begin
   t^.Removed:=true;
  end;
  TrapezoidalMap^.bcross:=Bottom;
  TrapezoidalMap^.tcross:=Top;
  e^.Above:=t1;
  e^.Below:=t2;
  Sink:=t^.Sink;
  t^.Sink:=nil;
  QueryGraphCase3(TrapezoidalMap^.QueryGraph,Sink,e,t1,t2);
  result:=Next;
 end;
 function TrapezoidalMapCase4(const TrapezoidalMap:PTrapezoidalMap;const t:PTrapezoid;const e:PEdge):PTrapezoid;
 var Next,t1,t2:PTrapezoid;
     Sink:PNode;
 begin
  Next:=TrapezoidalMapNextTrapezoid(TrapezoidalMap,t,e);
  if TrapezoidalMap^.tcross=t^.Top then begin
   t1:=t^.UpperLeft;
   t1^.RightPoint:=e^.q;
  end else begin
   t1:=NewTrapezoid(t^.LeftPoint,e^.q,t^.Top,e);
   TrapezoidUpdateLeft(t1,t^.UpperLeft,e^.Above);
   TrapezoidalMapPush(TrapezoidalMap,t1);
  end;
  if TrapezoidalMap^.bcross=t^.Bottom then begin
   t2:=t^.LowerLeft;
   t2^.RightPoint:=e^.q;
  end else begin
   t2:=NewTrapezoid(t^.LeftPoint,e^.q,e,t^.Bottom);
   TrapezoidUpdateLeft(t2,e^.Below,t^.LowerLeft);
   TrapezoidalMapPush(TrapezoidalMap,t2);
  end;
  t^.LeftPoint:=e^.q;
  TrapezoidUpdateLeft(t,t1,t2);
  Sink:=t^.Sink;
  t^.Sink:=nil;
  QueryGraphCase4(TrapezoidalMap^.QueryGraph,Sink,e,t1,t2,t);
  result:=Next;
 end;
 function TrapezoidalMapAddEdge(const TrapezoidalMap:PTrapezoidalMap;const Edge:PEdge):boolean;
 var t:PTrapezoid;
     cp,cq:boolean;
 begin
  result:=false;
  t:=QueryGraphLocate(TrapezoidalMap^.QueryGraph,Edge^.p,Edge^.Slope);
  if assigned(t) then begin
   cp:=false;
   cq:=false;
   while assigned(t) and (t<>FalseTrapezoid) do begin
    if cp then begin
     cp:=false;
    end else begin
     cp:=TrapezoidContains(t,Edge^.p);
    end;
    if cq then begin
     cq:=false;
    end else begin
     cq:=TrapezoidContains(t,Edge^.q);
    end;
    if cp and cq then begin
     t:=TrapezoidalMapCase1(TrapezoidalMap,t,Edge);
     break;
    end else if cp and not cq then begin
     t:=TrapezoidalMapCase2(TrapezoidalMap,t,Edge);
    end else if (not cp) and (not cq) then begin
     t:=TrapezoidalMapCase3(TrapezoidalMap,t,Edge);
    end else begin
     t:=TrapezoidalMapCase4(TrapezoidalMap,t,Edge);
    end;
    if not assigned(t) then begin
     exit;
    end;
   end;
   TrapezoidalMap^.bcross:=nil;
   TrapezoidalMap^.tcross:=nil;
   result:=true;
  end;
 end;
 function TrapezoidalMapCollectPoints(const TrapezoidalMap:PTrapezoidalMap):boolean;
 var Current:PTrapezoid;
 begin
  result:=false;
  Current:=TrapezoidalMap.First;
  while assigned(Current) do begin
   if not Current^.Removed then begin
    if (Current^.Top=TrapezoidalMap^.Root^.Top) and (assigned(Current^.Bottom^.Below) and not Current^.Bottom^.Below^.Removed) then begin
     TrapezoidMarkInside(Current^.Bottom^.Below);
     break;
    end;
    if (Current^.Bottom=TrapezoidalMap^.Root^.Bottom) and (assigned(Current^.Top^.Above) and not Current^.Top^.Above^.Removed) then begin
     TrapezoidMarkInside(Current^.Top^.Above);
     break;
    end;
   end;
   Current:=Current^.Next;
  end;
  Current:=TrapezoidalMap.First;
  while assigned(Current) do begin
   if Current^.Inside and not Current^.Removed then begin
    if IsInfinite(Current^.Top^.p^.y) then begin
     exit;
    end;
    TrapezoidAddPoints(Current);
   end;
   Current:=Current^.Next;
  end;
  result:=true;
 end;
 function ShearTransform(const p:TBeRoTriangulationPoint):PPoint;
 const SHEAR=1e-10;
 begin
  result:=NewPoint(p.x+(p.y*SHEAR),p.y);
 end;
var i,j,k,h,l:longint;
    InputPolygon:TBeRoTriangulationPolygon;
    p,q:PPoint;
    Edges:array of PEdge;
    Map:PTrapezoidalMap;
begin
 GarbageCollectorList:=nil;
 try
  result:=false;
  Edges:=nil;
  try
   h:=0;
   for i:=0 to length(InputPolygons)-1 do begin
    l:=length(InputPolygons[i]);
    if l<3 then begin
     exit;
    end;
    h:=h+l;
   end;
   if h>2 then begin
    SetLength(Edges,h);
    h:=0;
    for i:=0 to length(InputPolygons)-1 do begin
     InputPolygon:=InputPolygons[i];
     l:=length(InputPolygon);
     for j:=0 to l-1 do begin
      k:=j+1;
      if k>=l then begin
       k:=0;
      end;
      if j<>0 then begin
       p:=q;
      end else begin
       p:=ShearTransform(InputPolygon[j]);
      end;
      q:=ShearTransform(InputPolygon[k]);
      if p^.x>q^.x then begin
       Edges[h]:=NewEdge(q,p);
      end else begin
       Edges[h]:=NewEdge(p,q);
      end;
      inc(h);
     end;
    end;
    Map:=NewTrapezoidalMap;
    FalseTrapezoid:=NewTrapezoid(nil,nil,nil,nil);
    for i:=0 to h-1 do begin
     if not TrapezoidalMapAddEdge(Map,Edges[i]) then begin
      exit;
     end;
    end;
    if not TrapezoidalMapCollectPoints(Map) then begin
     exit;
    end;
    CountOutputTriangles:=0;
    for i:=0 to h-1 do begin
     if assigned(Edges[i]^.Polygon) and (Edges[i]^.Polygon^.Length>0) then begin
      TriangulateMountain(Edges[i]);
     end;
    end;
    SetLength(OutputTriangles,CountOutputTriangles);
    result:=length(OutputTriangles)>0;
   end;
  finally
   SetLength(Edges,0);
  end;
 finally
  GarbageCollectorCleanUp;
 end;
end;
{$endif}

function TriangulateEarClipping(const InputPolygon:TBeRoTriangulationPolygon;var OutputTriangles:TBeRoTriangulationTriangles):boolean;
const EPSILON=1e-14;
type PVertex=^TVertex;
     TVertex=record
      ArrayPrevious:PVertex;
      ArrayNext:PVertex;
      Previous:PVertex;
      Next:PVertex;
      x,y:double;
     end;
     TVertices=array of TVertex;
var Vertices:TVertices;
    i,j,k,CountOutputTriangles:longint;
    Vertex,Current,Start,Stop,Next:PVertex;
    a,dx0,dy0,dx1,dy1,px,py,vx,vy,nx,ny,vnx,vny,npx,npy,pvx,pvy,cx,cy,pcx,pcy,vcx,vcy,ncx,ncy:double;
    IsEar,Convex,FirstSign:boolean;
 procedure PushTriangle(const ax,ay,bx,by,cx,cy:double);
 var Triangle:PBeRoTriangulationTriangle;
 begin
  if length(OutputTriangles)<=CountOutputTriangles then begin
   SetLength(OutputTriangles,(CountOutputTriangles+1)*2);
  end;
  Triangle:=@OutputTriangles[CountOutputTriangles];
  inc(CountOutputTriangles);
  if (((bx-ax)*(by+ay))+((cx-bx)*(cy+by))+((ax-cx)*(ay+cy)))<0.0 then begin
   Triangle^[0].x:=ax;
   Triangle^[0].y:=ay;
   Triangle^[1].x:=bx;
   Triangle^[1].y:=by;
   Triangle^[2].x:=cx;
   Triangle^[2].y:=cy;
  end else begin
   Triangle^[0].x:=cx;
   Triangle^[0].y:=cy;
   Triangle^[1].x:=bx;
   Triangle^[1].y:=by;
   Triangle^[2].x:=ax;
   Triangle^[2].y:=ay;
  end;
 end;
begin
 CountOutputTriangles:=0;
 if length(InputPolygon)=3 then begin
  PushTriangle(InputPolygon[0].x,InputPolygon[0].y,InputPolygon[1].x,InputPolygon[1].y,InputPolygon[2].x,InputPolygon[2].y);
 end else if length(InputPolygon)>3 then begin
  Convex:=true;
  FirstSign:=false;
  for i:=0 to length(InputPolygon)-1 do begin
   j:=i+1;
   if j>=length(InputPolygon) then begin
    j:=0;
   end;
   k:=j+1;
   if k>=length(InputPolygon) then begin
    k:=0;
   end;
   dx0:=InputPolygon[k].x-InputPolygon[j].x;
   dy0:=InputPolygon[k].y-InputPolygon[j].y;
   dx1:=InputPolygon[i].x-InputPolygon[j].x;
   dy1:=InputPolygon[i].y-InputPolygon[j].y;
   a:=(dx0*dy1)-(dx1*dy0);
   if i=0 then begin
    FirstSign:=a<0.0;
   end else begin
    if FirstSign<>(a<0.0) then begin
     Convex:=false;
     break;
    end;
   end;
  end;
  if Convex then begin
   SetLength(OutputTriangles,length(InputPolygon)-3);
   for i:=1 to length(InputPolygon)-2 do begin
    PushTriangle(InputPolygon[0].x,InputPolygon[0].y,InputPolygon[i].x,InputPolygon[i].y,InputPolygon[i+1].x,InputPolygon[i+1].y);
   end;
  end else begin
   Vertices:=nil;
   try
    SetLength(Vertices,length(InputPolygon));
    for i:=0 to length(Vertices)-1 do begin
     Vertices[i].x:=InputPolygon[i].x;
     Vertices[i].y:=InputPolygon[i].y;
    end;
    for i:=1 to length(Vertices)-1 do begin
     Vertices[i].ArrayPrevious:=@Vertices[i-1];
    end;
    for i:=0 to length(Vertices)-2 do begin
     Vertices[i].ArrayNext:=@Vertices[i+1];
    end;
    Vertices[length(Vertices)-1].ArrayNext:=@Vertices[0];
    Vertices[0].ArrayPrevious:=@Vertices[length(Vertices)-1];
    Start:=@Vertices[0];
    Vertex:=Start;
    Start:=Vertex;
    Current:=Vertex;
    a:=0.0;
    repeat
     Next:=Current^.ArrayNext;
     a:=a+((Current^.x-Next^.x)*(Current^.y+Next^.y));
     Current:=Next;
    until Current=Start;
    if a<0.0 then begin
     repeat
      Vertex^.Previous:=Vertex^.ArrayNext;
      Vertex^.Next:=Vertex^.ArrayPrevious;
      Vertex:=Vertex^.ArrayPrevious;
     until Vertex=Start;
    end else begin
     repeat
      Vertex^.Previous:=Vertex^.ArrayPrevious;
      Vertex^.Next:=Vertex^.ArrayNext;
      Vertex:=Vertex^.ArrayNext;
     until Vertex=Start;
    end;
    Vertex:=Start;
    while Vertex^.Previous<>Vertex^.Next do begin
     px:=Vertex^.Previous^.x;
     py:=Vertex^.Previous^.y;
     vx:=Vertex^.x;
     vy:=Vertex^.y;
     nx:=Vertex^.Next^.x;
     ny:=Vertex^.Next^.y;
     vnx:=nx-vx;
     vny:=ny-vy;
     npx:=px-nx;
     npy:=py-ny;
     pvx:=vx-px;
     pvy:=vy-py;
     IsEar:=((pvx*vny)-(vnx*pvy))>EPSILON;
     if IsEar then begin
      Current:=Vertex^.Next;
      Stop:=Vertex^.Previous^.Previous;
      while Current<>Stop do begin
       Current:=Current^.Next;
       cx:=Current^.x;
       cy:=Current^.y;
       pcx:=cx-px;
       pcy:=cy-py;
       if (pcx<>0.0) or (pcy<>0.0) then begin
        vcx:=cx-vx;
        vcy:=cy-vy;
        if (vcx<>0.0) or (vcy<>0.0) then begin
         ncx:=cx-nx;
         ncy:=cy-ny;
         if (ncx<>0.0) or (ncy<>0.0) then begin
          if (((vnx*vcy)-(vny*vcx))>=(-EPSILON)) and (((pvx*pcy)-(pvy*pcx))>=(-EPSILON)) and (((npx*ncy)-(npy*ncx))>=(-EPSILON)) then begin
           IsEar:=false;
           break;
          end;
         end;
        end;
       end;
      end;
     end;
     if IsEar then begin
      PushTriangle(Vertex^.Previous^.x,Vertex^.Previous^.y,Vertex^.x,Vertex^.y,Vertex^.Next^.x,Vertex^.Next^.y);
      Vertex^.Previous^.Next:=Vertex^.Next;
      Vertex^.Next^.Previous:=Vertex^.Previous;
      Vertex:=Vertex^.Next;
      Start:=Vertex;
     end else begin
      Vertex:=Vertex^.Next;
      if Vertex=Start then begin
       SetLength(Vertices,0);
       exit;
      end;
     end;
    end;
   finally
    SetLength(Vertices,0);
   end;
  end;
 end;
 SetLength(OutputTriangles,CountOutputTriangles);
 result:=CountOutputTriangles>0;
end;

end.