unit TeeTriangulate;
{$I TeeDefs.inc}

// Adapted from Andre Bester's algorithm (anb@iafrica.com)

interface

uses
  {$IFDEF FMX}
  FMXTee.Canvas
  {$ELSE}
  TeCanvas
  {$ENDIF}
  ;

type
  TTriValue=Double; // Single

  TGetValueProc=function(Index:Integer):TTriValue of object;
  TSetValueProc=procedure(Dest,Source:Integer) of object;

  TTriangulator=class
  private
    function IDxchg(const I1,I2,I3,I4:Integer):Integer;
  public
    IPL : packed Array of Integer;
    IPT : packed Array of Integer;

    GetX : TGetValueProc;
    GetZ : TGetValueProc;
    SetXYZ : TSetValueProc;

    NumLines     : Integer;
    NumTriangles : Integer;
    ImprovedTriangles : Boolean;

    Destructor Destroy; override;

    function CreateTriangles(NumPoints:Integer):Boolean;
  end;

implementation

uses
  Math;

Destructor TTriangulator.Destroy;
begin
  IPL:=nil;
  IPT:=nil;

  inherited;
end;

function TTriangulator.IDxchg(const I1,I2,I3,I4:Integer):Integer;
var x1,x2,x3,x4 : TTriValue;
    y1,y2,y3,y4 : TTriValue;
    u1,u2,u3,u4 : TTriValue;
    A1,A2,B1,B2 : TTriValue;
    C1,C2       : TTriValue;
    S1,S2,S3,S4 : TTriValue;
begin
  result:=0;

  x1:=GetX(I1);
  y1:=GetZ(I1);
  x2:=GetX(I2);
  y2:=GetZ(I2);
  x3:=GetX(I3);
  y3:=GetZ(I3);
  x4:=GetX(I4);
  y4:=GetZ(I4);

  u3:=(y2-y3)*(x1-x3)-(x2-x3)*(y1-y3);
  u4:=(y1-y4)*(x2-x4)-(x1-x4)*(y2-y4);

  if (u3*u4)>0 then
  begin
    u1:=(y3-y1)*(x4-x1)-(x3-x1)*(y4-y1);
    u2:=(y4-y2)*(x3-x2)-(x4-x2)*(y3-y2);

    A1:=sqr(x1-x3)+sqr(y1-y3);
    B1:=sqr(x4-x1)+sqr(y4-y1);
    C1:=sqr(x3-x4)+sqr(y3-y4);
    A2:=sqr(x2-x4)+sqr(y2-y4);
    B2:=sqr(x3-x2)+sqr(y3-y2);
    C2:=sqr(x2-x1)+sqr(y2-y1);

    S1:=Sqr(u1)/(C1*Math.Max(A1,B1));
    S2:=Sqr(u2)/(C1*Math.Max(A2,B2));
    S3:=Sqr(u3)/(C2*Math.Max(B2,A1));
    S4:=Sqr(u4)/(C2*Math.Max(B1,A2));

    if Math.Min(S1,S2) < Math.Min(S3,S4) then
       result:=1;
  end;
end;

function TTriangulator.CreateTriangles(NumPoints:Integer):Boolean;
var NLT3 : Integer;
    ITF  : packed Array[1..2] of Integer;
    ipl1 : Integer;
    ipl2 : Integer;
    IPTI1: Integer;
    IPTI2: Integer;

  Procedure CalcBorder;
  var i    : Integer;
      JLT3 : Integer;
      IPLJ1: Integer;
      IPLJ2: Integer;
  begin
    for i:=1 to NLT3 div 3 do
    begin
      JLT3:=i*3;
      IPLJ1:=IPL[JLT3-2];
      IPLJ2:=IPL[JLT3-1];

      if ((IPLJ1=ipl1) and (IPLJ2=IPTI2)) or
         ((IPLJ2=ipl1) and (IPLJ1=IPTI2)) then
         IPL[JLT3]:=ITF[1];

      if ((IPLJ1=ipl2) and (IPLJ2=IPTI1)) or
         ((IPLJ2=ipl2) and (IPLJ1=IPTI1)) then
         IPL[JLT3]:=ITF[2];
    end;
  end;

  Var DSQMN : TTriValue;
      IPMN1 : Integer;
      IPMN2 : Integer;
      ip1   : Integer;
      ip2   : Integer;
      NDPM1 : Integer;
      xd1   : TTriValue;
      xd2   : TTriValue;
      yd1   : TTriValue;
      yd2   : TTriValue;
      NDP0  : Integer;
      DSQI  : TTriValue;
   tmpCount : Integer;

  { find closest pair and their midpoint }
  Function FindClosestPair:Boolean;
  begin
    result:=False;
    DSQMN:=Sqr(GetX(1)-GetX(0))+Sqr(GetZ(1)-GetZ(0));

    IPMN1:=1;
    IPMN2:=2;
    ip1:=1;
    while (not result) and (ip1<=NDPM1) do
    begin
      xd1:=GetX(IP1);
      yd1:=GetZ(IP1);

      ip2:=ip1+1;
      while (not result) and (ip2<=NDP0) do
      begin
        xd2:=GetX(IP2);
        yd2:=GetZ(IP2);
        DSQI:=Sqr(xd2-xd1)+Sqr(yd2-yd1);

        if DSQI=0.0 then
        begin
          SetXYZ(ip2,NDP0);

          Dec(tmpCount);
          Dec(NDP0);
          Dec(NDPM1);
          Dec(ip2);
        end
        else
        if DSQI<DSQMN then
        begin
          DSQMN:=DSQI;
          IPMN1:=ip1;
          IPMN2:=ip2;
        end;

        Inc(ip2);
      end;

      Inc(ip1);
    end;
  end;

var
  JPMN : Integer;
  IWL  : Array of Integer;
  IWP  : Array of Integer;
  WK   : Array of TTriValue;

  Procedure SortRest;
  Var XDMP : TTriValue;
      YDMP : TTriValue;
      jp1  : Integer;
      jp2  : Integer;
      tmpip1 : Integer;
      tmp  : Integer;
  begin
    XDMP:=(GetX(IPMN1)+GetX(IPMN2))*0.5;
    YDMP:=(GetZ(IPMN1)+GetZ(IPMN2))*0.5;

    // sort other (NDP-2) datapoints in ascending order of
    // distance from midpoint and stores datapoint numbers
    // in IWP array

    jp1:=2;
    for tmpip1:=1 to NDP0 do
    if (tmpip1<>IPMN1) and (tmpip1<>IPMN2) then
    begin
      Inc(jp1);
      IWP[jp1]:=tmpip1;
      WK[jp1]:=Sqr(GetX(tmpIP1)-XDMP)+Sqr(GetZ(tmpIP1)-YDMP);
    end;

    for jp1:=3 to NDPM1 do
    begin
      DSQMN:=WK[jp1];
      JPMN:=jp1;
      for jp2:=jp1 to NDP0 do
      begin
        // optimized...
        if WK[jp2]<DSQMN then
        begin
          DSQMN:=WK[jp2];
          JPMN:=jp2;
        end;
      end;

      tmp:=IWP[jp1];
      IWP[jp1]:=IWP[JPMN];
      IWP[JPMN]:=tmp;
      WK[JPMN]:=WK[jp1];
    end;
  end;

Var DSQ12 : TTriValue;
Const Ratio = 1.0E-6;
    NRep  = 100;

  Procedure CheckColinear;
  Var AR       : TTriValue;
      dx21     : TTriValue;
      dy21     : TTriValue;
      CoLinear : TTriValue;
      jp       : Integer;
      ip       : Integer;
      jpmx     : Integer;
      i        : Integer;
  begin
    // if necessary modifies ordering so that first
    // three datapoints are not colinear
    AR:=DSQ12*ratio;

    xd1:=GetX(IPMN1);
    yd1:=GetZ(IPMN1);

    dx21:=GetX(IPMN2)-xd1;
    dy21:=GetZ(IPMN2)-yd1;

    ip:=0;
    jp:=3;
    CoLinear:=0.0;

    while (jp<=NDP0) and (colinear<=AR) do
    begin
      ip:=IWP[jp];
      CoLinear:=Abs((GetZ(IP)-yd1)*dx21-(GetX(IP)-xd1)*dy21);
      Inc(jp);
    end;

    Dec(jp);

//      if jp=NDP0 then
//         raise ETriSurfaceException.Create(TeeMsg_TriSurfaceAllColinear);

    if jp<>3 then
    if Length(IWP)>3 then
    begin
      jpmx:=jp;
      jp:=jpmx+1;

      for i:=4 to jpmx do
      begin
        Dec(jp);
        IWP[jp]:=IWP[jp-1];
      end;

      IWP[3]:=ip;
    end;
  end;

Var NTT3 : Integer;

  // forms first triangle-vertices in IPT array and border
  // line segments and triangle number in IPL array
  Procedure AddFirst;

    function Side(Const u1,v1,u2,v2,u3,v3:TTriValue):TTriValue;
    begin
      result:=(v3-v1)*(u2-u1)-(u3-u1)*(v2-v1);
    end;

  Var ip3 : Integer;
  begin
    ip1:=IPMN1;
    ip2:=IPMN2;

    if Length(IWP)>3 then
       ip3:=IWP[3]
    else
       ip3:=0;

    if Side(GetX(IP1),GetZ(IP1),GetX(IP2),GetZ(IP2),GetX(IP3),GetZ(IP3))<0 then
    begin
      ip1:=IPMN2;
      ip2:=IPMN1;
    end;

    NumTriangles:=1;
    NTT3:=3;

    { first triangle }
    IPT[1]:=ip1;
    IPT[2]:=ip2;
    IPT[3]:=ip3;

    NumLines:=3;
    NLT3:=9;
    IPL[1]:=ip1;
    IPL[2]:=ip2;
    IPL[3]:=1;
    IPL[4]:=ip2;
    IPL[5]:=ip3;
    IPL[6]:=1;
    IPL[7]:=ip3;
    IPL[8]:=ip1;
    IPL[9]:=1;
  end;

  Procedure CalcTriangle(jp1:Integer);
  Var DXMN : TTriValue;
      DYMN : TTriValue;
      ARMN,
      dxmx,dymx,
      dsqmx,armx : TTriValue;
      NSH,JWL : Integer;
      NLN,NLNT3,
      ITT3,
      NLF : Integer;
      tmp,
      jpmx : Integer;

    Procedure Part1;
    var jp2 : Integer;
        AR  : TTriValue;
        DX  : TTriValue;
        DY  : TTriValue;
    begin
      for jp2:=2 to NumLines do
      begin
        ip2:=IPL[3*jp2-2];

        xd2:=GetX(IP2);
        yd2:=GetZ(IP2);

        DX:=xd2-xd1;
        DY:=yd2-yd1;

        AR:=DY*DXMN-DX*DYMN;
        if AR<=ARMN then
        begin
          DSQI:=Sqr(DX)+Sqr(DY);

          if (AR<-ARMN) or (DSQI<DSQMN) then
          begin
            JPMN:=jp2;
            DXMN:=DX;
            DYMN:=DY;
            DSQMN:=DSQI;
            ARMN:=DSQMN*ratio;
          end;
        end;

        AR:=DY*DXMX-DX*DYMX;
        if AR>=-ARMX then
        begin
          DSQI:=Sqr(DX)+Sqr(DY);

          if (AR>ARMX) or (DSQI<DSQMX) then
          begin
            JPMX:=jp2;
            DXMX:=DX;
            DYMX:=DY;
            DSQMX:=DSQI;
            ARMX:=DSQMX*ratio;
          end;
        end;
      end;
    end;

    Procedure ShiftIPLArray;
    var i : Integer;
        tmpSource : Integer;
    begin
      // shifts the IPL array to have invisible border
      // line segments contained in 1st part of array
      for i:=1 to NSH do
      begin
        tmp:=i*3;
        tmpSource:=tmp+NLT3;
        IPL[tmpSource-2]:=IPL[tmp-2];
        IPL[tmpSource-1]:=IPL[tmp-1];
        IPL[tmpSource]  :=IPL[tmp];
      end;

      for i:=1 to NLT3 div 3 do
      begin
        tmp:=i*3;
        tmpSource:=tmp+(NSH*3);
        IPL[tmp-2]:=IPL[tmpSource-2];
        IPL[tmp-1]:=IPL[tmpSource-1];
        IPL[tmp]  :=IPL[tmpSource];
      end;

      Dec(JPMX,NSH);
    end;

    Procedure AddTriangles;
    var jp2   : Integer;
        IPTI  : Integer;
        IT    : Integer;
        jp2t3 : Integer;
    begin
      // adds triangles to IPT array, updates border line
      // segments in IPL array and sets flags for the border
      // line segments to be reexamined in the iwl array
      JWL:=0;
      NLNT3:=0;
      for jp2:=JPMX to NumLines do
      begin
        jp2t3:=jp2*3;
        ipl1:=IPL[jp2t3-2];
        ipl2:=IPL[jp2t3-1];
        IT:=IPL[jp2t3];

        // add triangle to IPT array
        Inc(NumTriangles);
        Inc(NTT3,3);
        IPT[NTT3-2]:=ipl2;
        IPT[NTT3-1]:=ipl1;
        IPT[NTT3]:=ip1;

        // updates borderline segments in ipl array
        if jp2=JPMX then
        begin
          IPL[jp2t3-1]:=ip1;
          IPL[jp2t3]:=NumTriangles;
        end;

        if jp2=NumLines then
        begin
          NLN:=JPMX+1;
          NLNT3:=NLN*3;
          IPL[NLNT3-2]:=ip1;
          IPL[NLNT3-1]:=IPL[1];
          IPL[NLNT3]:=NumTriangles;
        end;

        // determine vertex that is not on borderline segments
        ITT3:=IT*3;
        IPTI:=IPT[ITT3-2];

        if (IPTI=ipl1) or (IPTI=ipl2) then
        begin
          IPTI:=IPT[ITT3-1];

          if (IPTI=ipl1) or (IPTI=ipl2) then
             IPTI:=IPT[ITT3];
        end;

        // checks if exchange is necessary
        if IDxchg(ip1,IPTI,ipl1,ipl2)<>0 then
        begin
          // modifies ipt array if necessary
          IPT[ITT3-2]:=IPTI;
          IPT[ITT3-1]:=ipl1;
          IPT[ITT3]:=ip1;
          IPT[NTT3-1]:=IPTI;

          if jp2=JPMX then IPL[jp2t3]:=IT;

          if (jp2=NumLines) and (IPL[3]=IT) then
             IPL[3]:=NumTriangles;

          // set flags in IWL array
          JWL:=JWL+4;
          IWL[JWL-3]:=ipl1;
          IWL[JWL-2]:=IPTI;
          IWL[JWL-1]:=IPTI;
          IWL[JWL]:=ipl2;
        end;
      end;
    end;

    Procedure ImproveTriangles;
    Var ILF    : Integer;
        tmpNLF : Integer;
        IPT1   : Integer;
        IPT2   : Integer;
        IPT3   : Integer;
        IREP   : Integer;
        IT1T3  : Integer;
        IT2T3  : Integer;
        LoopFlag : Boolean;
        NTF    : Integer;
        NTT3P3 : Integer;
    begin
      // improve triangulation
      NTT3P3:=NTT3+3;
      IREP:=1;
      while IREP<=NREP do
      begin
        for ILF:=1 to NLF do
        begin
          ipl1:=IWL[ILF*2-1];
          ipl2:=IWL[ILF*2];

          // locates in ipt array two triangles on
          // both sides of flagged line segment
          NTF:=0;
          LoopFlag:=True;
          tmp:=3;
            
          while LoopFlag and (tmp<=NTT3) do
          begin
            ITT3:=NTT3P3-tmp;

            IPT1:=IPT[ITT3-2];
            IPT2:=IPT[ITT3-1];
            IPT3:=IPT[ITT3];

            if (ipl1=IPT1) or (ipl1=IPT2) or (ipl1=IPT3) then // todo: optimize?
            begin
              if (ipl2=IPT1) or (ipl2=IPT2) or (ipl2=IPT3) then
              begin
                Inc(NTF);
                ITF[NTF]:=ITT3 div 3;

                if NTF=2 then
                   LoopFlag:=False;
              end;
            end;

            Inc(tmp,3);
          end;

          if NTF>=2 then
          begin
            IT1T3:=ITF[1]*3;
            IPTI1:=IPT[IT1T3-2];

            if (IPTI1=ipl1) or (IPTI1=ipl2) then
            begin
              IPTI1:=IPT[IT1T3-1];
              if (IPTI1=ipl1) or (IPTI1=ipl2) then IPTI1:=IPT[IT1T3];
            end;

            IT2T3:=ITF[2]*3;
            IPTI2:=IPT[IT2T3-2];

            if (IPTI2=ipl1) or (IPTI2=ipl2) then
            begin
              IPTI2:=IPT[IT2T3-1];
              if (IPTI2=ipl1) or (IPTI2=ipl2) then IPTI2:=IPT[IT2T3];
            end;

             // checks if exchange necessary
            if IDxchg(IPTI1,IPTI2,ipl1,ipl2)<>0 then
            begin
               IPT[IT1T3-2]:=IPTI1;
               IPT[IT1T3-1]:=IPTI2;
               IPT[IT1T3]:=ipl1;

               IPT[IT2T3-2]:=IPTI2;
               IPT[IT2T3-1]:=IPTI1;
               IPT[IT2T3]:=ipl2;

               JWL:=JWL+8;

               IWL[JWL-7]:=ipl1;
               IWL[JWL-6]:=IPTI1;
               IWL[JWL-5]:=IPTI1;
               IWL[JWL-4]:=ipl2;
               IWL[JWL-3]:=ipl2;
               IWL[JWL-2]:=IPTI2;
               IWL[JWL-1]:=IPTI2;
               IWL[JWL]  :=ipl1;

               CalcBorder;
            end;
          end;
        end;

        tmp:=NLF;
        NLF:=JWL div 2;
        if NLF=tmp then break
        else
        begin // reset IWL array for next round
          JWL:=0;
          tmp:=(tmp+1)*2;
          tmpNLF:=2*NLF;
          while tmp<=tmpNLF do
          begin
            Inc(JWL,2);
            IWL[JWL-1]:=IWL[tmp-1];
            IWL[JWL]  :=IWL[tmp];
            Inc(tmp,2);
          end;
          NLF:=JWL div 2;
        end;

        Inc(IREP);
      end;
    end;

  begin
    ip1:=IWP[jp1];

    xd1:=GetX(IP1);
    yd1:=GetZ(IP1);

    // determine visible borderline segments
    ip2:=IPL[1];
    JPMN:=1;

    xd2:=GetX(IP2);
    yd2:=GetZ(IP2);

    DXMN:=xd2-xd1;
    DYMN:=yd2-yd1;
    DSQMN:=Sqr(DXMN)+Sqr(DYMN);

    ARMN:=DSQMN*Ratio;
    jpmx:=1;
    dxmx:=DXMN;
    dymx:=DYMN;
    dsqmx:=DSQMN;
    armx:=ARMN;

    Part1;

    if jpmx<jpmn then
       Inc(jpmx,NumLines);

    NSH:=JPMN-1;
    if NSH>0 then
       ShiftIPLArray;

    AddTriangles;

    NumLines:=NLN;
    NLT3:=NLNT3;
    NLF:=JWL div 2;

    if (NLF<>0) and ImprovedTriangles then
       ImproveTriangles;
  end;

var t : Integer;
begin
  result:=False;

  tmpCount:=NumPoints;
  NDP0:=tmpCount-1;
  NDPM1:=NDP0-1;

  SetLength(IPT,1+(6*tmpCount-15));
  SetLength(IPL,6*tmpCount);
  SetLength(IWL,18*tmpCount);
  SetLength(IWP,tmpCount);
  SetLength(WK,tmpCount);
  try
    if not FindClosestPair then
    begin
      DSQ12:=DSQMN;
      SortRest;
      CheckColinear;
      AddFirst;

      // Add the remaining NDP-3 data points one by one:
      for t:=4 to NDP0 do
          CalcTriangle(t);

      result:=True;
    end;
  finally
    IWL:=nil;
    IWP:=nil;
    WK:=nil;
  end;
end;

end.
