delaunay.pas

等值线插值Pascal程序

//Credit to Paul Bourke ([e-mail]) for the original Fortran 77 Program :))
//Conversion to Visual Basic by EluZioN ([e-mail])
//Conversion from VB to Delphi6 by Dr Steve Evans ([e-mail])
///////////////////////////////////////////////////////////////////////////////
//June 2002 Update by Dr Steve Evans ([e-mail]): 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.
///////////////////////////////////////////////////////////////////////////////
//You can use this code however you like providing the above credits remain in tact

//注意：点和三角形都数组都是从1开始计数的

unit Delaunay;

interface

uses Dialogs, Graphics, Forms, Types,classes,math;

//Set these as applicable
Const
  MaxVertices = 500000;
  MaxTriangles = 1000000;
  ExPtTolerance = 0.000001;  //小于这个被认为是同一点


Type
  TCastArray = Array [0..2,0..2,0..2] of Integer;
  TVectorL3D = Array [0..2] of Double;
  TVectorL3I = Array [0..2] of Integer;

  PPointPair = ^TPointPair;
  TPointPair = record
    x1,y1,
    x2,y2: Double
  end;
  //单条等值线
  TLever = record
    FZ: Double;
    Points: TList;
  end;

  //Points (Vertices)
  dVertex = record
    X ,
    Y ,
    Z: Double;
  end;

  //Created Triangles, vv# are the vertex pointers(点的索引)
  dTriangle = record
    vv0: LongInt;
    vv1: LongInt;
    vv2: LongInt;
    PreCalc: Integer;
    xc,yc,r: Double;   //三角形外接圆圆心坐标和半径
  end;

  TDVertexs = array[0..MaxVertices] of dVertex;
  PVertexs = ^TDVertexs;

  TDTriangles = array[0..MaxTriangles] of dTriangle;
  PTriangles = ^TDTriangles;

  TDCompletes = array [0..MaxTriangles] of Boolean;
  PCompletes = ^TDCompletes;

  TDEdges = array[0..2,0..MaxTriangles * 3] of LongInt;
  PEdges = ^TDEdges;

  TDelaunay = class
  private
    { Private declarations }
    FzLow,
    FzHigh: Double;

    FVertexs: PVertexs;
    FTriangles: PTriangles;
    FTriangleCount: Integer;
    FPointCount: Integer;  //Variable for total number of points (vertices)
    procedure QuickSort(var aVertexs: PVertexs; Low,High: Integer);
    function GetPointCount: integer;
    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;

  public
    { Public declarations }
    FLevers: Array of TLever;
    TempBuffer: TBitmap;
    TargetForm: TForm;
    constructor Create;
    destructor Destroy; override;
    procedure Mesh;
    procedure Draw;
    procedure ScatterContour(ZCount: Integer; Z: Array of Single);
    procedure AddPoint(x,y,z: Single);
    procedure ClearBackPage;
    procedure FlipBackPage;
    property zLow: Double read FzLow write FzLow;
    property zHigh: Double read FzHigh write FzHigh;
    property Vertexs: PVertexs read FVertexs;
    property Triangles: PTriangles read FTriangles;
    property TriangleCount: Integer read FTriangleCount;
    property PointCount: Integer read GetPointCount;
  end;

implementation

constructor TDelaunay.Create;
begin
  //Initiate total points to 1, using base 0 causes problems in the functions
  FPointCount := 1;
  FTriangleCount:=0;
  FzLow:= 0;
  FzHigh:= 0;
  TempBuffer:=TBitmap.Create;
  //Allocate memory for arrays
  GetMem(FVertexs, sizeof(FVertexs^));
  GetMem(FTriangles, sizeof(FTriangles^));
end;

destructor TDelaunay.Destroy;
begin
  //Free memory for arrays
  FreeMem(FVertexs, sizeof(FVertexs^));
  FreeMem(FTriangles, sizeof(FTriangles^));
end;

//加入点到FVertexs数组里
procedure TDelaunay.AddPoint(x,y,z: Single);
var
  i: Integer;
  SamePoint: Boolean;
begin
  //Check for duplicate points 检查是否有完全相同的点，
  //如果有则,该点不被加入
  SamePoint := false;
  i := 1;
  while i < FPointCount do
  begin
    If (Abs(x-FVertexs^[i].X) < ExPtTolerance) and
       (Abs(y-FVertexs^[i].Y) < ExPtTolerance) Then
      SamePoint:= true;
    Inc(i);
  end;

  if FzLow > z then
    FzLow:= z
  else if FzHigh < z then
    FzHigh:= z;

  if not SamePoint  then
  begin
    //Set Vertex coordinates
    FVertexs^[FPointCount].X := x;
    FVertexs^[FPointCount].Y := y;
    FVertexs^[FPointCount].Z := z;
    //Increment the total number of points
    //最后得到的点的数目会比实际数目多一个
    FPointCount := FPointCount + 1;
  end;
end;


//构建三角网
procedure TDelaunay.Mesh;
begin
  QuickSort(FVertexs,1,FPointCount-1);
  If FPointCount > 3 Then
  FTriangleCount := Triangulate(FPointCount-1); //'Returns number of triangles created.
end;

//点按X坐标从小到大排序
procedure TDelaunay.QuickSort(var aVertexs: PVertexs; Low,High: Integer);
//Sort all points by x
  procedure DoQuickSort(var aVertexs: PVertexs; iLo, iHi: Integer);
  var
    Lo, Hi: Integer;
    Mid: Double;
    T: dVertex;
  begin
    Lo := iLo;
    Hi := iHi;
    Mid := aVertexs^[(Lo + Hi) div 2].X;
    repeat
      while aVertexs^[Lo].x < Mid do Inc(Lo);
      while aVertexs^[Hi].x > Mid do Dec(Hi);
      if Lo <= Hi then
      begin
        T := aVertexs^[Lo];
        aVertexs^[Lo] := aVertexs^[Hi];
        aVertexs^[Hi] := T;
        Inc(Lo);
        Dec(Hi);
      end;
    until Lo > Hi;
    if Hi > iLo then DoQuickSort(aVertexs, iLo, Hi);
    if Lo < iHi then DoQuickSort(aVertexs, Lo, iHi);
  end;
begin
  DoQuickSort(aVertexs, Low, High);
end;


//真正构建三角网(nVert：点的个数)
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
  Completes: PCompletes;
  Edges: PEdges;
  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;  //是否在外接圆中
begin
  //Allocate memory
  GetMem(Completes, sizeof(Completes^));
  GetMem(Edges, sizeof(Edges^));

  //Find the maximum and minimum vertex bounds.
  //This is to allow calculation of the bounding triangle
  xmin := FVertexs^[1].x;
  ymin := FVertexs^[1].y;
  xmax := xmin;
  ymax := ymin;
  For i := 2 To nvert do
  begin
  If FVertexs^[i].x < xmin Then xmin := FVertexs^[i].x;
  If FVertexs^[i].x > xmax Then xmax := FVertexs^[i].x;
  If FVertexs^[i].y < ymin Then ymin := FVertexs^[i].y;
  If FVertexs^[i].y > ymax Then ymax := FVertexs^[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.

  FVertexs^[nvert + 1].x := (xmid - 2 * dmax);
  FVertexs^[nvert + 1].y := (ymid - dmax);
  FVertexs^[nvert + 2].x := xmid;
  FVertexs^[nvert + 2].y := (ymid + 2 * dmax);
  FVertexs^[nvert + 3].x := (xmid + 2 * dmax);
  FVertexs^[nvert + 3].y := (ymid - dmax);
  FTriangles^[1].vv0 := nvert + 1;
  FTriangles^[1].vv1 := nvert + 2;
  FTriangles^[1].vv2 := nvert + 3;
  FTriangles^[1].Precalc := 0;

  Completes[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 Completes^[j] <> True Then
      begin
        inc := InCircle(FVertexs^[i].x, FVertexs^[i].y, FVertexs^[FTriangles^[j].vv0].x,
                        FVertexs^[FTriangles^[j].vv0].y, FVertexs^[FTriangles^[j].vv1].x,
                        FVertexs^[FTriangles^[j].vv1].y, FVertexs^[FTriangles^[j].vv2].x,
                        FVertexs^[FTriangles^[j].vv2].y, xc, yc, r,j);
            //Include this if points are sorted by X
        If (xc + r) < FVertexs[i].x Then  //
          completes[j] := True          //
        Else If inc Then
        begin
          Edges^[1, Nedge + 1] := FTriangles^[j].vv0;
          Edges^[2, Nedge + 1] := FTriangles^[j].vv1;
          Edges^[1, Nedge + 2] := FTriangles^[j].vv1;
          Edges^[2, Nedge + 2] := FTriangles^[j].vv2;
          Edges^[1, Nedge + 3] := FTriangles^[j].vv2;
          Edges^[2, Nedge + 3] := FTriangles^[j].vv0;
          Nedge := Nedge + 3;
          FTriangles^[j].vv0 := FTriangles^[ntri].vv0;
          FTriangles^[j].vv1 := FTriangles^[ntri].vv1;
          FTriangles^[j].vv2 := FTriangles^[ntri].vv2;
          FTriangles^[j].PreCalc:=FTriangles^[ntri].PreCalc;
          FTriangles^[j].xc:=FTriangles^[ntri].xc;
          FTriangles^[j].yc:=FTriangles^[ntri].yc;
          FTriangles^[j].r:=FTriangles^[ntri].r;
          FTriangles^[ntri].PreCalc:=0;
          Completes^[j] := Completes^[ntri];
          j := j - 1;
          ntri := ntri - 1;
        End;//else
      End;  //if
    until j>=ntri; //repeat

    // 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
      If Not (Edges^[1, j] = 0) And Not (Edges^[2, j] = 0) Then
        For k := j + 1 To Nedge do
          If Not (Edges^[1, k] = 0) And Not (Edges^[2, k] = 0) Then
            If Edges^[1, j] = Edges^[2, k] Then
              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;

    //  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
      If Not (Edges^[1, j] = 0) And Not (Edges^[2, j] = 0) Then
      begin
        ntri := ntri + 1;
        FTriangles^[ntri].vv0 := Edges^[1, j];
        FTriangles^[ntri].vv1 := Edges^[2, j];
        FTriangles^[ntri].vv2 := i;
        FTriangles^[ntri].PreCalc:=0;