delaunay.pas

等值线插值Pascal程序

        Completes^[ntri] := False;
      End;
      
  end; //the first for

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

  Triangulate := ntri;

  //Free memory
  FreeMem(Completes, sizeof(Completes^));
  FreeMem(Edges, sizeof(Edges^));
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;
begin
  eps:= 0.000001;
  InCircle := False;

  //Check if xc,yc and r have already been calculated
  if  FTriangles^[j].PreCalc=1 then
  begin
    xc := FTriangles^[j].xc;
    yc := FTriangles^[j].yc;
    r  := FTriangles^[j].r;
    rsqr := r*r;
    dx := xp - xc;
    dy := yp - yc;
    drsqr := dx * dx + dy * 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;//else

    dx := x2 - xc;
    dy := y2 - yc;
    rsqr := dx * dx + dy * dy;
    r := Sqrt(rsqr);
    dx := xp - xc;
    dy := yp - yc;
    drsqr := dx * dx + dy * dy;

    //store the xc,yc and r for later use
    FTriangles^[j].PreCalc:=1;
    FTriangles^[j].xc:=xc;
    FTriangles^[j].yc:=yc;
    FTriangles^[j].r:=r;
  end;  //the big else

  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 equation > 0 Then
     WhichSide := -1
  Else If equation = 0 Then
     WhichSide := 0
  Else
     WhichSide := 1;
End;


procedure TDelaunay.Draw;
var
  i: Integer;
begin
  // Clear the form canvas
  ClearBackPage;

  TempBuffer.Canvas.Brush.Color := clwhite;
  //Draw the created triangles
  if (FTriangleCount > 0) then
    For i:= 1 To FTriangleCount do
    begin
      TempBuffer.Canvas.Polygon([Point(Trunc(FVertexs^[FTriangles^[i].vv0].x), Trunc(FVertexs^[FTriangles^[i].vv0].y)),
                               Point(Trunc(FVertexs^[FTriangles^[i].vv1].x), Trunc(FVertexs^[FTriangles^[i].vv1].y)),
                               Point(Trunc(FVertexs^[FTriangles^[i].vv2].x), Trunc(FVertexs^[FTriangles^[i].vv2].y))]);
    end;
  FlipBackPage;
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;



function TDelaunay.GetPointCount: integer;
begin
  Result:= FPointCount-1;
end;


procedure TDelaunay.ScatterContour(ZCount: Integer; Z: Array of Single);
var
  i,j,m: Integer;
  Deside: Integer;
  CastTab : TCastArray;

  sH      : TVectorL3I;
  H,xH,yH : TVectorL3D;
  
  TempD1,TempD2,dMin,dMax: Double ;
  x1,x2,y1,y2: Double; //等值点坐标

  ARecord: PPointPair; //记录点对

  //插值计算
  Function xSec(p1,p2:Integer): Double;
  Begin
    result:= (H[p2]*xH[p1]-H[p1]*xH[p2])/(H[p2]-H[p1]);
  End;

  Function ySec(p1,p2:Integer): Double;
  Begin
    result:= (H[p2]*yH[p1]-H[p1]*yH[p2])/(H[p2]-H[p1]);
  End;

begin
  //分配记录等值线的数组
  for i:= 0 to Length(FLevers)-1 do
    if Assigned(FLevers[i].Points) then
      FLevers[i].Points.Free;
  SetLength(FLevers,ZCount);
  for i:= 0 to ZCount-1 do
  begin
    FLevers[i].FZ:= Z[i];
    FLevers[i].Points:= TList.Create;
  end;

  //每个三角行内出现等值点的情况映射，有27种情况
  //这27种情况是根据三角形的三个顶点高程与等值点
  //的大小比较得来得，每个点有三种情况：大、小、等
  //0..19 为 对各种情况的处理方法，有20种
  CastTab[0,0,0]:= 0; CastTab[0,0,1]:= 0; CastTab[0,0,2]:= 1;
  CastTab[0,1,0]:= 0; CastTab[0,1,1]:= 2; CastTab[0,1,2]:= 3;
  CastTab[0,2,0]:= 4; CastTab[0,2,1]:= 5; CastTab[0,2,2]:= 6;

  CastTab[1,0,0]:= 0; CastTab[1,0,1]:= 7; CastTab[1,0,2]:= 8;
  CastTab[1,1,0]:= 9; CastTab[1,1,1]:= 10; CastTab[1,1,2]:= 9;
  CastTab[1,2,0]:= 8; CastTab[1,2,1]:= 7; CastTab[1,2,2]:= 0;

  CastTab[2,0,0]:= 6; CastTab[2,0,1]:= 5; CastTab[2,0,2]:= 4;
  CastTab[2,1,0]:= 3; CastTab[2,1,1]:= 2; CastTab[2,1,2]:= 0;
  CastTab[2,2,0]:= 1; CastTab[2,2,1]:= 0; CastTab[2,2,2]:= 0;

  for i:= 1 to TriangleCount do
  begin

    //获得三角形三个顶点中的最小值和最大值
    TempD1:= min(FVertexs^[FTriangles^[i].vv0].Z,FVertexs^[FTriangles^[i].vv1].Z);
    TempD2:= min(FVertexs^[FTriangles^[i].vv1].Z,FVertexs^[FTriangles^[i].vv2].Z);
    dMin:= min(TempD1,TempD2);
    TempD1:= max(FVertexs^[FTriangles^[i].vv0].Z,FVertexs^[FTriangles^[i].vv1].Z);
    TempD2:= max(FVertexs^[FTriangles^[i].vv1].Z,FVertexs^[FTriangles^[i].vv2].Z);
    dMax:= max(TempD1,TempD2);

    for j:= 0 to ZCount-1 do
      if (Z[j] >= dMin) And (Z[j] <= dMax) Then
      begin

        H[0] := FVertexs^[FTriangles^[i].vv0].Z-Z[j];
        xH[0]:= FVertexs^[FTriangles^[i].vv0].X;
        yH[0]:= FVertexs^[FTriangles^[i].vv0].Y;
        H[1] := FVertexs^[FTriangles^[i].vv1].Z-Z[j];
        xH[1]:= FVertexs^[FTriangles^[i].vv1].X;
        yH[1]:= FVertexs^[FTriangles^[i].vv1].Y;
        H[2] := FVertexs^[FTriangles^[i].vv2].Z-Z[j];
        xH[2]:= FVertexs^[FTriangles^[i].vv2].X;
        yH[2]:= FVertexs^[FTriangles^[i].vv2].Y;

        for m:= 0 to 2 do
          If H[m] > 0 Then
            sH[m]:= 1
          Else If H[m]<0 Then
            sH[m]:= -1
          Else
            sH[m]:= 0;

        Deside := CastTab[sH[0]+1 ,sH[1]+1, sH[2]+1];

        If NOT(deside = 0) Then // 0的情况不处理
        begin
          Case deside Of
                  1:  begin
                        x1:= xSec(0,2);
                        y1:= ySec(0,2);
                        x2:= xSec(1,2);
                        y2:= ySec(1,2);
                      end;
                  2:  begin
                        x1:= xH[1];
                        y1:= yH[1];
                        x2:= xH[2];
                        y2:= yH[2];
                      end;
                  3:  begin
                        x1:= xH[1];
                        y1:= yH[1];
                        x2:= xSec(0,2);
                        y2:= ySec(0,2);
                      end;
                  4:  begin
                        x1:= xSec(0,1);
                        y1:= ySec(0,1);
                        x2:= xSec(1,2);
                        y2:= ySec(1,2);
                      end;
                  5:  Begin
                        x1:= xH[2];
                        y1:= yH[2];
                        x2:= xSec(0,1);
                        y2:= ySec(0,1);
                      End;
                  6:  Begin
                        x1:= xSec(0,1);
                        y1:= ySec(0,1);
                        x2:= xSec(0,2);
                        y2:= ySec(0,2);
                      End;
                  7:  Begin
                        x1:= xH[0];
                        y1:= yH[0];
                        x2:= xH[2];
                        y2:= yH[2];
                      End;
                  8:  Begin
                        x1:= xH[0];
                        y1:= yH[0];
                        x2:= xSec(1,2);
                        y2:= ySec(1,2);
                      End;
                  9:  Begin
                        x1:= xH[0];
                        y1:= yH[0];
                        x2:= xH[1];
                        y2:= yH[1];
                      End;
                  10: begin         //there is some argument here
                        x1:= xH[0];
                        y1:= yH[0];
                        x2:= xH[2];
                        y2:= yH[2];
                      end;
                end;//----case

          //此处获得该三角形内的等值点
          New(ARecord);
          ARecord^.x1:= x1;
          ARecord^.y1:= y1;
          ARecord^.x2:= x2;
          ARecord^.y2:= y2;
          FLevers[j].Points.Add(ARecord);
        end; //if not(deside)
      end;// if Z[]
  end;
end;

end.