unit uDelaunay;
interface
uses
SysUtils, Dialogs, Graphics, Types, Math, uUtils, uClasses;
type
TDelaunayTriangle =
record
public
Indices:
array [0..2]
of integer;
// Indices in der Knotenliste FPoints;
procedure Clear;
end;
TDelaunayTriangles =
class
private
procedure SetCount(
const Value: integer);
procedure Flip(
const Index1, Index2, Edge1, Edge2: integer);
public
Item:
array of TDelaunayTriangle;
function Center(
const Index: integer;
var R: double;
Nodes: TFloatPoints): TFLoatPoint;
function IndexOf(
const A, B, C: integer): integer;
function IndexOfPtInCircumcircle(
const NodeIndex: integer;
Nodes: TFloatPoints): integer;
function IndexOfPtInCircumcircleEx(
const NodeIndex, StarIndex: integer;
Nodes: TFloatPoints): integer;
function IndexOfEdge(
const NodeIndex,
Index: integer;
const Nodes: TFloatPoints): integer;
procedure Add(
const A, B, C: integer);
procedure AddItem(
const Value: TDelaunayTriangle);
procedure Delete(
const Index: integer);
procedure Assign(Value: TDelaunayTriangles);
procedure Clear;
function Count: integer;
destructor Destroy;
override;
end;
TDelaunayTriangulation =
class
private
FBorder, FPoints: TFloatPoints;
FTriangles: TDelaunayTriangles;
procedure CheckCircumcircles;
procedure Triangulate;
procedure TriangulatePoint(
const Index: integer);
public
procedure Clear;
procedure Init;
procedure AddPoint(
const Value: TFloatPoint);
procedure AddBorderPoint(
const Value: TFloatPoint);
procedure AddPoints(Value: TFloatPoints);
procedure AddBorderPoints(Value: TFloatPoints);
procedure Draw(Canvas: TCanvas;
const ppMM: double);
property Points: TFloatPoints
read FPoints;
// mm;
property Border: TFloatPoints
read FBorder;
// mm;
property Triangles: TDelaunayTriangles
read FTriangles;
constructor Create;
destructor Destroy;
override;
end;
implementation
function CircumcircleCenter(
const A, B, C: TFloatPoint;
var R: double): TFloatPoint;
var
M1, M2, T1, T2: TFloatPoint;
begin
M1 := FloatPoint((A.X + B.X) / 2, (A.Y + B.Y) / 2);
// Mitte AB;
M2 := FloatPoint((A.X + C.X) / 2, (A.Y + C.Y) / 2);
// Mitte AC;
T1 := Util_RotateFloatPoint(A, M1, 0, 1);
// A um M1 90° drehen; 0 = Cos(90); 1 = Sin(90);
T2 := Util_RotateFloatPoint(A, M2, 0, 1);
// A um M2 90° drehen; ..
if Util_IntersectLinesEx(M1, T1, M2, T2, Result)
then // Schnittpunkt M1T1, M2T2;
R := Util_FloatPointDistance(A, Result)
else
R := 0;
end;
function PtInCircle(
const P, Center: TFloatPoint;
const R: double): boolean;
begin
Result := Util_PtInEllipse(P, Center, R, R);
end;
{ TDelaunayTriangle }
procedure TDelaunayTriangle.Clear;
begin
Indices[0] := -1;
Indices[1] := -1;
Indices[2] := -1;
end;
{ TDelaunayTriangles }
destructor TDelaunayTriangles.Destroy;
begin
Clear;
inherited Destroy;
end;
procedure TDelaunayTriangles.SetCount(
const Value: integer);
begin
SetLength(Item, Value);
end;
procedure TDelaunayTriangles.Assign(Value: TDelaunayTriangles);
var
I: integer;
begin
Clear;
for I := 0
to Value.Count - 1
do
AddItem(Value.Item[I]);
end;
procedure TDelaunayTriangles.Clear;
begin
SetCount(0);
end;
function TDelaunayTriangles.Count: integer;
begin
Result := Length(Item);
end;
procedure TDelaunayTriangles.Add(
const A, B, C: integer);
begin
SetCount(Count + 1);
Item[Count - 1].Indices[0] := A;
Item[Count - 1].Indices[1] := B;
Item[Count - 1].Indices[2] := C;
end;
procedure TDelaunayTriangles.AddItem(
const Value: TDelaunayTriangle);
begin
SetCount(Count + 1);
Item[Count - 1] := Value;
end;
procedure TDelaunayTriangles.Delete(
const Index: integer);
var
I: integer;
begin
for I :=
Index to Count - 2
do
Item[I] := Item[I + 1];
SetCount(Count - 1);
end;
function TDelaunayTriangles.IndexOf(
const A, B, C: integer): integer;
var
I: integer;
B1, B2, B3, B4, B5, B6: boolean;
begin
Result := -1;
for I := 0
to Count - 1
do
begin
B1 := (Item[I].Indices[0] = A)
and (Item[I].Indices[1] = B)
and (Item[I].Indices[2] = C);
B2 := (Item[I].Indices[0] = A)
and (Item[I].Indices[1] = C)
and (Item[I].Indices[2] = B);
B3 := (Item[I].Indices[0] = B)
and (Item[I].Indices[1] = A)
and (Item[I].Indices[2] = C);
B4 := (Item[I].Indices[0] = B)
and (Item[I].Indices[1] = C)
and (Item[I].Indices[2] = A);
B5 := (Item[I].Indices[0] = C)
and (Item[I].Indices[1] = A)
and (Item[I].Indices[2] = B);
B6 := (Item[I].Indices[0] = C)
and (Item[I].Indices[1] = B)
and (Item[I].Indices[2] = A);
if B1
or B2
or B3
or B4
or B5
or B6
then
begin
Result := I;
Break;
end;
end;
end;
function TDelaunayTriangles.IndexOfEdge(
const NodeIndex,
Index: integer;
const Nodes: TFloatPoints): integer;
var
IndexA, IndexB, IndexC: integer;
begin
Result := -1;
try
IndexA := Item[
Index].Indices[0];
IndexB := Item[
Index].Indices[1];
IndexC := Item[
Index].Indices[2];
if Util_SameFloatPoint(Nodes[NodeIndex], Nodes[IndexA])
then
Result := 0
else
if Util_SameFloatPoint(Nodes[NodeIndex], Nodes[IndexB])
then
Result := 1
else
if Util_SameFloatPoint(Nodes[NodeIndex], Nodes[IndexC])
then
Result := 2;
except
on E:
Exception do
ShowMessage('
Fehler vom Typ: ' + E.ClassName + '
, Meldung: ' + E.
Message);
end;
end;
procedure TDelaunayTriangles.Flip(
const Index1, Index2, Edge1, Edge2: integer);
begin
try
except
on E:
Exception do
ShowMessage('
Fehler vom Typ: ' + E.ClassName + '
, Meldung: ' + E.
Message);
end;
end;
function TDelaunayTriangles.IndexOfPtInCircumcircle(
const NodeIndex: integer;
Nodes: TFloatPoints): integer;
begin
Result := IndexOfPtInCircumcircleEx(NodeIndex, 0, Nodes);
end;
function TDelaunayTriangles.Center(
const Index: integer;
var R: double;
Nodes: TFloatPoints): TFLoatPoint;
var
IndexA, IndexB, IndexC: integer;
A, B, C: TFloatPoint;
begin
R := 0;
Result.Clear;
try
IndexA := Item[
Index].Indices[0];
IndexB := Item[
Index].Indices[1];
IndexC := Item[
Index].Indices[2];
A := Nodes[IndexA];
B := Nodes[IndexB];
C := Nodes[IndexC];
Result := CircumcircleCenter(A, B, C, R);
except
on E:
Exception do
ShowMessage('
Fehler vom Typ: ' + E.ClassName + '
, Meldung: ' + E.
Message);
end;
end;
function TDelaunayTriangles.IndexOfPtInCircumcircleEx(
const NodeIndex, StarIndex: integer;
Nodes: TFloatPoints): integer;
var
I: integer;
C: TFloatPoint;
R: double;
begin
Result := -1;
try
for I := StarIndex
to Count - 1
do
begin
C := Center(I, R, Nodes);
if PtInCircle(Nodes[NodeIndex], C, R)
then
Result := I;
// No Break;
end;
except
on E:
Exception do
ShowMessage('
Fehler vom Typ: ' + E.ClassName + '
, Meldung: ' + E.
Message);
end;
end;
{ TDelaunayTriangulation }
constructor TDelaunayTriangulation.Create;
begin
inherited;
FPoints := TFloatPoints.Create;
FBorder := TFloatPoints.Create;
FTriangles := TDelaunayTriangles.Create;
end;
destructor TDelaunayTriangulation.Destroy;
begin
FPoints.Free;
FBorder.Free;
FTriangles.Free;
inherited;
end;
procedure TDelaunayTriangulation.Init;
// need Border;
var
xMin, xMax, yMin, yMax: double;
begin
// 0....1....2;
// . . . .
// . . . .
// . . . .
// 3.........4;
xMin := FBorder.Left;
yMin := FBorder.Top;
xMax := FBorder.Right;
yMax := FBorder.Bottom;
FPoints.Clear;
FPoints.AddXY(xMin, yMin);
// 0;
FPoints.AddXY(xMax / 2, yMin);
// 1;
FPoints.AddXY(xMax, yMin);
// 2;
FPoints.AddXY(xMin, yMax);
// 3;
FPoints.AddXY(xMax, yMax);
// 4;
FTriangles.Clear;
FTriangles.Add(0, 1, 3);
FTriangles.Add(1, 4, 3);
FTriangles.Add(1, 2, 4);
end;
procedure TDelaunayTriangulation.Clear;
begin
FBorder.Clear;
FPoints.Clear;
FTriangles.Clear;
end;
procedure TDelaunayTriangulation.AddPoint(
const Value: TFloatPoint);
// Need Border;
begin
if (FPoints.IndexOf(Value) < 0)
and (FBorder.PtInPolygon(Value))
then
begin
FPoints.Add(Value);
TriangulatePoint(FPoints.Count - 1);
end;
end;
procedure TDelaunayTriangulation.AddBorderPoint(
const Value: TFloatPoint);
begin
if FBorder.IndexOf(Value) < 0
then
FBorder.Add(Value);
end;
procedure TDelaunayTriangulation.AddPoints(Value: TFloatPoints);
var
I: integer;
begin
for I := 0
to Value.Count - 1
do
AddPoint(Value[I]);
end;
procedure TDelaunayTriangulation.AddBorderPoints(Value: TFloatPoints);
var
I: integer;
begin
for I := 0
to Value.Count - 1
do
if FBorder.IndexOf(Value[I]) < 0
then
FBorder.Add(Value[I]);
end;
procedure TDelaunayTriangulation.CheckCircumcircles;
// ???
var
I, J, Index1, Index2, Edge1, Edge2: integer;
begin
EXIT;
I := 0;
J := 0;
while I < FPoints.Count
do
begin
while J < FTriangles.Count - 1
do
begin
Index1 := FTriangles.IndexOfPtInCircumcircleEx(I, J, FPoints);
if Index1 > -1
then
begin
Index2 := FTriangles.IndexOfPtInCircumcircleEx(I, J + 1, FPoints);
if Index2 > -1
then
begin
Edge1 := FTriangles.IndexOfEdge(I, Index1, FPoints);
Edge2 := FTriangles.IndexOfEdge(I, Index2, FPoints);
FTriangles.Flip(Index1, Index2, Edge1, Edge2);
I := 0;
J := -1;
end;
end;
Inc(J);
end;
Inc(I);
end;
end;
procedure TDelaunayTriangulation.Triangulate;
var
TriangleIndex, I, IndexA, IndexB, IndexC: integer;
begin
for I := 0
to FPoints.Count - 1
do
begin
TriangleIndex := FTriangles.IndexOfPtInCircumcircle(I, FPoints);
if TriangleIndex > -1
then
begin
IndexA := FTriangles.Item[TriangleIndex].Indices[0];
IndexB := FTriangles.Item[TriangleIndex].Indices[1];
IndexC := FTriangles.Item[TriangleIndex].Indices[2];
FTriangles.Add(IndexA, IndexB, I);
FTriangles.Add(IndexB, IndexC, I);
FTriangles.Add(IndexC, IndexA, I);
end;
end;
CheckCircumcircles;
end;
procedure TDelaunayTriangulation.TriangulatePoint(
const Index: integer);
var
TriangleIndex, IndexA, IndexB, IndexC: integer;
begin
TriangleIndex := FTriangles.IndexOfPtInCircumcircle(
Index, FPoints);
if TriangleIndex > -1
then
begin
IndexA := FTriangles.Item[TriangleIndex].Indices[0];
IndexB := FTriangles.Item[TriangleIndex].Indices[1];
IndexC := FTriangles.Item[TriangleIndex].Indices[2];
FTriangles.Add(IndexA, IndexB,
Index);
FTriangles.Add(IndexB, IndexC,
Index);
FTriangles.Add(IndexC, IndexA,
Index);
end;
CheckCircumcircles;
end;
procedure TDelaunayTriangulation.Draw(Canvas: TCanvas;
const ppMM: double);
var
R: double;
I, IndexA, IndexB, IndexC: integer;
A, B, C, D: TPoint;
Center: TFloatPoint;
begin
Canvas.Brush.Color := clWhite;
Canvas.Brush.Style := bsSolid;
Canvas.FillRect(Canvas.ClipRect);
For I := 0
To FTriangles.Count - 1
do
begin
IndexA := FTriangles.Item[I].Indices[0];
IndexB := FTriangles.Item[I].Indices[1];
IndexC := FTriangles.Item[I].Indices[2];
A := Util_CanvasPoint(FPoints[IndexA], ppMM);
B := Util_CanvasPoint(FPoints[IndexB], ppMM);
C := Util_CanvasPoint(FPoints[IndexC], ppMM);
Center := FTriangles.Center(I, R, FPoints);
D := Util_CanvasPoint(Center, ppMM);
// ShowMessage(Format('%d %d %d', [IndexA, IndexB, IndexC]));
Canvas.Brush.Style := bsClear;
Canvas.Polygon([A, B, C]);
Canvas.Brush.Color := clWhite;
Canvas.Brush.Style := bsSolid;
Canvas.TextOut(A.X, A.Y, IntToStr(IndexA));
Canvas.TextOut(B.X, B.Y, IntToStr(IndexB));
Canvas.TextOut(C.X, C.Y, IntToStr(IndexC));
// Canvas.TextOut(D.X, D.Y, IntToStr(I));
A := Util_CanvasPoint(FloatPoint(Center.X - R, Center.Y - R), ppMM);
B := Util_CanvasPoint(FloatPoint(Center.X + R, Center.Y + R), ppMM);
Canvas.Brush.Style := bsClear;
// Canvas.Ellipse(A.X, A.Y, B.X, B.Y);
end;
end;
end.