Delaunay-Triangulation

[Bjoerk]

Zitat von Medium:

Ich habe es zwar schon gemacht, aber so komplett ohne Ansatz und Ausgangslage ist sinnvolle Hilfe nicht möglich. Fang an, und bei konkreten Problemen sind wir sicher gerne dabei - ich zumindest, weil ich finde solche Themen spannend

Zitat von Jens01:

Ich glaube, er macht erst so ein Supertriangle. Das soll wahrscheinlich den äußeren Rand beschreiben, was aber ein Problem ist, wenn Du eine bestimmte Außenkontur haben willst und Punkte für ein Loch hast.[..]

Genau so sieht‘s aus. Wie baust du denn die Ränder ein hast und hast du eine Idee für einen Flip-Algorithmus? Das Einarbeiten der Ränder und überhaupt erst mal ein Raster an Punkten zu bekommen das dann trianguliert wird scheint mir hier die Hauptarbeit zu werden? Hab den Algo übrigens völlig neu gemacht (vom Link oben keine einzige Zeile übernommen). Hier mal soweit wie ich bis jetzt bin.

Delphi-Quellcode:

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.