//Credit to Paul Bourke (pbourke@swin.edu.au) for the original Fortran 77 Program :))
//Conversion to Visual Basic by EluZioN (EluZioN@casesladder.com)
//Conversion from VB to Delphi6 by Dr Steve Evans (steve@lociuk.com)
///////////////////////////////////////////////////////////////////////////////
//June 2002 Update by Dr Steve Evans (steve@lociuk.com): Heap memory allocation
//added to prevent stack overflow when MaxVertices and MaxTriangles are very large.
//Additional Updates in June 2002:
//Bug in InCircle function fixed. Radius r := Sqrt(rsqr).
//Check for duplicate points added when inserting new point.
//For speed, all points pre-sorted in x direction using quicksort algorithm and
//triangles flagged when no longer needed. The circumcircle centre and radius of
//the triangles are now stored to improve calculation time.
//
//June 2005 (More) Dynamic memory allocation and a few optimizations added by
//Gunnar Blumert (gunnar@blumert.de)
///////////////////////////////////////////////////////////////////////////////
//You can use this code however you like providing the above credits remain in tact

unit Delaunay;

interface

uses Windows, SysUtils, Dialogs, Graphics, Forms, Types, Math;

// Set these values as appropriate for your application
// If these values are not sufficient, more memory will
// be allocated as needed
const InitialVertexCount = 10000;
      Vertex_Increment   = 1000;
      ExPtTolerance = 0.000001;

//Points (Vertices)
Type dVertex = record
    x,y : double;
end;

//Created Triangles, vv# are the vertex pointers
Type dTriangle = record
    vv0: LongInt;
    vv1: LongInt;
    vv2: LongInt;
    PreCalc: Integer;
    xc,yc,r: Double;
end;

type TDVertex = array of dVertex;

TDVertexClass = class
private
  fData: TDVertex;
  procedure SetItem(Index: integer; Value: dVertex);
  function GetItem(Index: integer): dVertex;
public
  constructor Create(Capacity: integer);
  destructor Destroy; override;
  property Items[index: integer]: dVertex read getitem write setitem; default;
end;

type TDTriangle = array of dTriangle;

TDTriangleClass = class
private
  fData: TDTriangle;
  procedure SetItem(Index: integer; value: dTriangle);
  function GetItem(Index: integer): dTriangle;
public
  constructor Create(Capacity: integer);
  destructor destroy; override;
  property Items[index: integer]: dTriangle read getitem write setitem; default;
end;

TByteArray = array of byte;
TBoolArrayClass = class
private
  fData: TByteArray;
  function GetItem(Index: integer): boolean;
  procedure SetItem(Index: integer; Value: boolean);
  function ByteIndex(Index: integer): integer;
  function BitMask(Index: integer): byte;
public
  constructor create(Capacity: integer);
  destructor Destroy; override;
  property Items[index: integer]: boolean read GetItem write Setitem; default;
end;

TIntArray = array of integer;
TIntArrayClass = class
private
  fData: TIntArray;
  function GetItem(Index: integer): integer;
  procedure SetItem(Index: integer; Value: integer);
public
  constructor create(Capacity: integer);
  destructor destroy; override;
  property Items[index: integer]: integer read getitem write setitem; default;
end;

TDEdges = array[0..2] of TIntArrayClass;

type
  TDelaunay = class
  private
    { Private declarations }
    Vertex: TDVertexClass;
    Triangle: TDTriangleClass;
    function InCircle(xp, yp, x1, y1, x2, y2, x3, y3: Double;
             var xc: Double; var yc: Double; var r: Double; j: Integer): Boolean;
    Function WhichSide(xp, yp, x1, y1, x2, y2: Double): Integer;
    Function Triangulate(nvert: Integer): Integer;
    function FindPoint(_x, _y: Double; out index: integer): boolean;
  public
    { Public declarations }
    TempBuffer: TBitmap;
    HowMany: Integer;
    tPoints: Integer; //Variable for total number of points (vertices)
    TargetForm: TForm;
    constructor Create;
    destructor Destroy; override;
    procedure Mesh;
    procedure Draw;
    procedure AddPoint(x,y: Double);
    procedure ClearBackPage;
    procedure FlipBackPage;
    procedure QuickSort(Low,High: Integer);
  end;

implementation

constructor TDVertexClass.Create(Capacity: integer);
begin
  inherited create;
  SetLength(fData,Capacity);
end;

destructor TDVertexClass.Destroy;
begin
  finalize(fData);
  inherited;
end;

function TDVertexClass.GetItem(Index: integer): dVertex;
begin
  if High(fData) < Index then SetLength(fData,Index+Vertex_Increment);
  Result := fData[Index]
end;

procedure TDVertexClass.SetItem(Index: integer; Value: dVertex);
begin
  if High(fData) < Index then SetLength(fData,Index+Vertex_Increment);
  fData[Index] := Value
end;

constructor TDTriangleClass.Create(Capacity: integer);
begin
  inherited create;
  SetLength(fData,Capacity);
end;

destructor TDTriangleClass.destroy;
begin
  Finalize(fData);
  inherited
end;

function TDTriangleClass.GetItem(Index: integer): dTriangle;
begin
  if High(fData) < Index then SetLength(fData,Index+Vertex_Increment*2);
  Result := fData[index]
end;

procedure TDTriangleClass.SetItem(Index: integer; value: dTriangle);
begin
  if High(fData) < Index then SetLength(fData,Index+Vertex_Increment*2);
  fData[index] := value
end;

constructor TBoolArrayClass.Create(Capacity: integer);
begin
  inherited create;
  SetLength(fData,Capacity shr 3);
end;

destructor TBoolArrayClass.Destroy;
begin
  Finalize(fData);
  inherited;
end;

function TBoolArrayClass.ByteIndex(Index: integer): integer;
begin
  Result := index shr 3;
end;

function TBoolArrayClass.BitMask(Index: integer): byte;
begin
  Result := 1 shl (Index and $7)
end;

function TBoolArrayClass.GetItem(Index: integer): boolean;
begin
  Result := fData[ByteIndex(Index)] and BitMask(Index) <> 0;
end;

procedure TBoolArrayClass.SetItem(Index: integer; Value: boolean);
var bi: integer;
begin
  bi := ByteIndex(Index);
  if bi > High(fData) then SetLength(fData,bi+(Vertex_Increment shr 3));
  if Value then fData[bi] := fData[bi] or BitMask(Index) else fData[bi] := fData[bi] and (not BitMask(Index));
end;

constructor TIntArrayClass.create(Capacity: integer);
begin
  inherited create;
  SetLength(fData,Capacity);
end;

destructor TIntArrayClass.Destroy;
begin
  Finalize(fData);
  inherited
end;

function TIntArrayClass.GetItem(Index: integer): integer;
begin
  if High(fData) < Index then SetLength(fDAta,Index+Vertex_Increment);
  Result := fData[index]
end;

procedure TIntArrayClass.SetItem(Index: integer; Value: integer);
begin
  if High(fData) < Index then SetLength(fData,Index+Vertex_Increment);
  fData[Index] := Value;
end;

constructor TDelaunay.Create;
begin
  inherited;
  //Initiate total points to 1, using base 0 causes problems in the functions
  tPoints := 1;
  HowMany:=0;
  TempBuffer:=TBitmap.Create;
  Vertex := TDVertexClass.Create(InitialVertexCount);
  Triangle := TDTriangleClass.Create(InitialVertexCount*2);
end;

destructor TDelaunay.Destroy;
begin
  FreeAndNIL(Vertex);
  FreeAndNIL(Triangle);
  inherited;
end;

function TDelaunay.FindPoint(_x,_y: double; out index: integer): boolean;
var i: integer;
    dx, dy: double;
begin
  Result := false;
  for i := 1 to pred(tPoints) do with Vertex[i] do begin
    dx := _x-x;
    dy := _y-y;
    if IsZero(dx,ExPtTolerance) and IsZero(dy,ExPtTolerance) then begin
      Result := true;
      index := i;
      exit
    end else if x > _x then exit;
  end;
end;

function TDelaunay.InCircle(xp, yp, x1, y1, x2, y2, x3, y3: Double;
    var xc: Double; var yc: Double; var r: Double; j: Integer): Boolean;
//Return TRUE if the point (xp,yp) lies inside the circumcircle
//made up by points (x1,y1) (x2,y2) (x3,y3)
//The circumcircle centre is returned in (xc,yc) and the radius r
//NOTE: A point on the edge is inside the circumcircle
var
  eps: Double;
  m1: Double;
  m2: Double;
  mx1: Double;
  mx2: Double;
  my1: Double;
  my2: Double;
  dx: Double;
  dy: Double;
  rsqr: Double;
  drsqr: Double;
  T: dTriangle;
begin
  eps:= 0.000001;
  InCircle := False;

  //Check if xc,yc and r have already been calculated
  if  Triangle[j].PreCalc=1 then
  begin
    xc := Triangle[j].xc;
    yc := Triangle[j].yc;
    r  := Triangle[j].r;
    rsqr := sqr(r);
    dx := xp - xc;
    dy := yp - yc;
    drsqr := sqr(dx) + sqr(dy);
  end else
  begin



  If (Abs(y1 - y2) < eps) And (Abs(y2 - y3) < eps) Then
  begin
  ShowMessage('INCIRCUM - F - Points are coincident !!');
  Exit;
  end;

  If Abs(y2 - y1) < eps Then
  begin
  m2 := -(x3 - x2) / (y3 - y2);
  mx2 := (x2 + x3) / 2;
  my2 := (y2 + y3) / 2;
  xc := (x2 + x1) / 2;
  yc := m2 * (xc - mx2) + my2;
  end
  Else If Abs(y3 - y2) < eps Then
  begin
  m1 := -(x2 - x1) / (y2 - y1);
  mx1 := (x1 + x2) / 2;
  my1 := (y1 + y2) / 2;
  xc := (x3 + x2) / 2;
  yc := m1 * (xc - mx1) + my1;
  end
  Else
  begin
  m1 := -(x2 - x1) / (y2 - y1);
  m2 := -(x3 - x2) / (y3 - y2);
  mx1 := (x1 + x2) / 2;
  mx2 := (x2 + x3) / 2;
  my1 := (y1 + y2) / 2;
  my2 := (y2 + y3) / 2;
    if (m1-m2)<>0 then  //se
    begin
    xc := (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
    yc := m1 * (xc - mx1) + my1;
    end else
    begin
    xc:= (x1+x2+x3)/3;
    yc:= (y1+y2+y3)/3;
    end;

  end;

  dx := x2 - xc;
  dy := y2 - yc;
  rsqr := dx * dx + dy * dy;
  r := Sqrt(rsqr);
  dx := xp - xc;
  dy := yp - yc;
  drsqr := sqr(dx) + sqr(dy);
  //store the xc,yc and r for later use
  T := Triangle[j];
  T.PreCalc := 1;
  T.xc := xc;
  T.yc := yc;
  T.r := r;
  Triangle[j] := T;

  end;

  If drsqr <= rsqr Then InCircle := True;

end;



Function TDelaunay.WhichSide(xp, yp, x1, y1, x2, y2: Double): Integer;
//Determines which side of a line the point (xp,yp) lies.
//The line goes from (x1,y1) to (x2,y2)
//Returns -1 for a point to the left
//         0 for a point on the line
//        +1 for a point to the right
var
 equation: Double;
begin
  equation := ((yp - y1) * (x2 - x1)) - ((y2 - y1) * (xp - x1));

  if IsZero(Equation) then Result := 0 else
  if Equation > 0 then Result := -1
  else Result := 1;

(*
  If equation > 0 Then
     WhichSide := -1
  Else If equation = 0 Then
     WhichSide := 0
  Else
     WhichSide := 1;
*)
End;



Function TDelaunay.Triangulate(nvert: Integer): Integer;
//Takes as input NVERT vertices in arrays Vertex()
//Returned is a list of NTRI triangular faces in the array
//Triangle(). These triangles are arranged in clockwise order.
var

//  Complete: TDComplete;
  Complete: TBoolArrayClass;
  Edges: TDEdges;
  Nedge: LongInt;

  //For Super Triangle
  xmin: Double;
  xmax: Double;
  ymin: Double;
  ymax: Double;
  xmid: Double;
  ymid: Double;
  dx: Double;
  dy: Double;
  dmax: Double;

  //General Variables
  i : Integer;
  j : Integer;
  k : Integer;
  ntri : Integer;
  xc : Double;
  yc : Double;
  r : Double;
  inc : Boolean;

  dv: DVertex;
  dt: dTriangle;
begin

      //Allocate memory
  Complete := TBoolArrayClass.create(tPoints shr 3);
  try
    for i := 0 to 2 do Edges[i] := TIntArrayClass.create(tPoints);
    try
      //Find the maximum and minimum vertex bounds.
      //This is to allow calculation of the bounding triangle
      xmin := Vertex[1].x;
      ymin := Vertex[1].y;
      xmax := xmin;
      ymax := ymin;
        For i := 2 To nvert do
        begin
        If Vertex[i].x < xmin Then xmin := Vertex[i].x;
        If Vertex[i].x > xmax Then xmax := Vertex[i].x;
        If Vertex[i].y < ymin Then ymin := Vertex[i].y;
        If Vertex[i].y > ymax Then ymax := Vertex[i].y;
        end;

      dx := xmax - xmin;
      dy := ymax - ymin;
      If dx > dy Then
          dmax := dx
      Else
          dmax := dy;

      xmid := Trunc((xmax + xmin) / 2);
      ymid := Trunc((ymax + ymin) / 2);

      //Set up the supertriangle
      //This is a triangle which encompasses all the sample points.
      //The supertriangle coordinates are added to the end of the
      //vertex list. The supertriangle is the first triangle in
      //the triangle list.
      dv.x := round(xmid - 2 * dmax);
      dv.y := round(ymid - dmax);
      Vertex[nvert+1] := dv;
      dv.x := round(xmid);
      dv.y := round(ymid + 2 * dmax);
      Vertex[nvert+2] := dv;
      dv.x := round(xmid + 2 * dmax);
      dv.y := round(ymid - dmax);
      Vertex[nvert+3] := dv;

      dt := Triangle[1];
      with dt do begin
        vv0 := nvert+1;
        vv1 := nvert+2;
        vv2 := nvert+3;
        Precalc := 0;
      end;
      Triangle[1] := dt;

      Complete[1] := False;
      ntri := 1;

      //Include each point one at a time into the existing mesh
      For i := 1 To nvert do
      begin
        Nedge := 0;
          //Set up the edge buffer.
          //If the point (Vertex(i).x,Vertex(i).y) lies inside the circumcircle then the
          //three edges of that triangle are added to the edge buffer.
          j := 0;
            repeat
              j := j + 1;
              If Complete[j] <> True Then
              begin
                  inc := InCircle(Vertex[i].x, Vertex[i].y, Vertex[Triangle[j].vv0].x,
                                  Vertex[Triangle[j].vv0].y, Vertex[Triangle[j].vv1].x,
                                  Vertex[Triangle[j].vv1].y, Vertex[Triangle[j].vv2].x,
                                  Vertex[Triangle[j].vv2].y, xc, yc, r,j);
                  //Include this if points are sorted by X
                  If (xc + r) < Vertex[i].x Then  //
                     complete[j] := True          //
                  Else                            //
                  If inc Then
                  begin
                      Edges[1][Nedge+1] := Triangle[j].vv0;
                      Edges[1][Nedge + 1] := Triangle[j].vv0;
                      Edges[2][Nedge + 1] := Triangle[j].vv1;
                      Edges[1][Nedge + 2] := Triangle[j].vv1;
                      Edges[2][Nedge + 2] := Triangle[j].vv2;
                      Edges[1][Nedge + 3] := Triangle[j].vv2;
                      Edges[2][Nedge + 3] := Triangle[j].vv0;
                      Nedge := Nedge + 3;
                      Triangle[j] := Triangle[ntri];
                      dt := Triangle[ntri];
                      dt.PreCalc := 0;
                      Triangle[ntri] := dt;
                      Complete[j] := Complete[ntri];
                      j := j - 1;
                      ntri := ntri - 1;
                  End;
              End;
          until j>=ntri;

      // Tag multiple edges
      // Note: if all triangles are specified anticlockwise then all
      // interior edges are opposite pointing in direction.
          For j := 1 To Nedge - 1 do
          begin
//              If Not (Edges[1][j] = 0) And Not (Edges[2][j] = 0) Then
              if (Edges[1][j] <> 0) or (Edges[2][j] <> 0) then
              begin
                  For k := j + 1 To Nedge do
                  begin
//                      If Not (Edges[1][k] = 0) And Not (Edges[2][k] = 0) Then
                      if (Edges[1][k] <> 0) or (Edges[2][k] <> 0) then
                      begin
                          If Edges[1][j] = Edges[2][k] Then
                          begin
                              If Edges[2][j] = Edges[1][k] Then
                               begin
                                  Edges[1][j] := 0;
                                  Edges[2][j] := 0;
                                  Edges[1][k] := 0;
                                  Edges[2][k] := 0;
                               End;
                           End;
                      End;
                  end;
              End;
          end;

      //  Form new triangles for the current point
      //  Skipping over any tagged edges.
      //  All edges are arranged in clockwise order.
          For j := 1 To Nedge do
          begin
//                  If Not (Edges[1][j] = 0) And Not (Edges[2][j] = 0) Then
                  if (Edges[1][j] <> 0) or (Edges[2][j] <> 0) then
                  begin
                      ntri := ntri + 1;
                      dt := Triangle[ntri];
                      with dt do begin
                        vv0 := Edges[1][j];
                        vv1 := Edges[2][j];
                        vv2 := i;
                        Precalc := 0;
                      end;
                      Triangle[ntri] := dt;
                      Complete[ntri] := False;
                  End;
          end;
      end;

      //Remove triangles with supertriangle vertices
      //These are triangles which have a vertex number greater than NVERT
      i:= 0;
      repeat
        i := i + 1;
        If (Triangle[i].vv0 > nvert) Or (Triangle[i].vv1 > nvert) Or (Triangle[i].vv2 > nvert) Then
        begin
          dt := Triangle[i];
          dt.vv0 := Triangle[ntri].vv0;
          dt.vv1 := Triangle[ntri].vv1;
          dt.vv2 := Triangle[ntri].vv2;
          Triangle[i] := dt;
           i := i - 1;
           ntri := ntri - 1;
        End;
      until i>=ntri;

      Triangulate := ntri;
    finally
      for i := 0 to 2 do FreeAndNIL(Edges[i]);
    end
  finally
    FreeAndNIL(Complete);
  end
End;


procedure TDelaunay.Mesh;
begin
  QuickSort(1,pred(tPoints));
  If tPoints > 3 Then
  HowMany := Triangulate(tPoints-1); //'Returns number of triangles created.
end;

procedure TDelaunay.Draw;
var
  //variable to hold how many triangles are created by the triangulate function
  i: Integer;
begin
  // Clear the form canvas
  ClearBackPage;

  TempBuffer.Canvas.Brush.Color := clTeal;
  //Draw the created triangles
  if (HowMany > 0) then
  begin
    For i:= 1 To HowMany do
    begin
    TempBuffer.Canvas.Polygon([Point(Trunc(Vertex[Triangle[i].vv0].x), -Trunc(Vertex[Triangle[i].vv0].y)),
                                  Point(Trunc(Vertex[Triangle[i].vv1].x), -Trunc(Vertex[Triangle[i].vv1].y)),
                                  Point(Trunc(Vertex[Triangle[i].vv2].x), -Trunc(Vertex[Triangle[i].vv2].y))]);
    end;
  end;

  FlipBackPage;
end;

procedure TDelaunay.AddPoint(x,y: Double);
var dv: DVertex;
    i: integer;
begin
  if FindPoint(x,y,i) then exit;
  dv.x := x;
  dv.y := -y;
  Vertex[tPoints] := dv;
  inc(tPoints);
end;

procedure TDelaunay.ClearBackPage;
begin
  TempBuffer.Height:=TargetForm.Height;
  TempBuffer.Width:=TargetForm.Width;
  TempBuffer.Canvas.Brush.Color := clSilver;
  TempBuffer.Canvas.FillRect(Rect(0,0,TargetForm.Width,TargetForm.Height));
end;

procedure TDelaunay.FlipBackPage;
var
  ARect : TRect;
begin
  ARect := Rect(0,0,TargetForm.Width,TargetForm.Height);
  TargetForm.Canvas.CopyRect(ARect, TempBuffer.Canvas, ARect);
end;


procedure TDelaunay.QuickSort(Low,High: Integer);
//Sort all points by x
  procedure DoQuickSort(iLo, iHi: Integer);
  var
    Lo, Hi: Integer;
    Mid: Double;
    T: dVertex;
    A: TDVertex;
  begin
    A:=Vertex.fData;
    Lo := iLo;
    Hi := iHi;
    Mid := A[(Lo + Hi) div 2].x;
    repeat
      while A[Lo].x < Mid do Inc(Lo);
      while A[Hi].x > Mid do Dec(Hi);
      if Lo <= Hi then
      begin
        T := A[Lo];
        A[Lo] := A[Hi];
        A[Hi] := T;
        Inc(Lo);
        Dec(Hi);
      end;
    until Lo > Hi;
    if Hi > iLo then DoQuickSort(iLo, Hi);
    if Lo < iHi then DoQuickSort(Lo, iHi);
  end;
begin
  DoQuickSort(Low, High);
end;

end.