接上
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;
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;
代码2接PART1
/构建三角网
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;
//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;
//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.
//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;
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;
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;