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.

[jfheins]

Ist es eine Option, einen vorhandenen Netzgenerator herzunehmen?
Da habe ich zum Beispiel Triangle gefunden. Dem kann man eine passende Datei geben, und er gibt die das Netz zurück.

Ansonsten: Advancing-Frint geht wie folgt vor: Du unterteilst deine Außenkontur zunächst in Knoten und Kanten. Da kannst du ganz stumpf machen, dass du eine Wunschlänge definierst, und dann an jeder Polygonkante schaust, wie viele Unterknoten da denn eingezogen werden müssen.

Advancing-Front geht jetzt schrittweise die Knoten durch, die auf der aktuellen Front liegen und verkleinert die Front (bzw. genauer: das eingeschlossene Gebiet). Wenn die Front leer ist, hast du das Gebiet vollständig vernetzt. Es gibt dabei drei Möglichkeiten:

1. Es wird ein neuer Punkt mit zwei Kanten erzeugt. Bevorzuge ein gleichseitiges Dreieck.
2. Es wird ein Punkt erzeugt, der drei Kanten erhält
3. Es wird eine Kante erzeugt

Welche Fall eintritt, hängt von dem Winkel der aktuellen Kante mit der nächsten Kante zusammen. Guckst du hier: http://www.iue.tuwien.ac.at/phd/fleischmann/node39.html der Winkel alpha.

Das Netz, dass daraus hervoirgeht, kannst du dann natürlich auch noch auf Delaunay überprüfen und Kanten ggf. flippen.

Allgemein würde ich dir aber empfehlen, einen fertigen Netzgenerator zu verwenden.

Edit: Habe gerade noch eine gute Arbeit gefunden: http://elib.uni-stuttgart.de/opus/vo.../geomaindt.pdf
Ab Seite 29.

[Bjoerk]

Ich fürchte du überschätzt (mal wieder) meine mathematischen Fähigkeiten. Ich hab' Anfang der 80iger studiert und spätestens Anfang der Neunziger so ziemlich wieder alles vergessen (mit Ausnahme der TM natürlich). Ich krieg das schon hin, aber anders. Ein PolygnonHatch hab' ich schon, da kann man sogar Winkel und Abstand einstellen (ist mir vorhin eingefallen). Das könnte ich als Punkegenerator verwenden. Ich kann es nur nicht einbauen weil solange meiner Delaunaytriangulation die Umkreisprüfung fehlt werdern die Dreiecke (teilweise) falsch. Deshalb die Frage nach einem eleganten Flip-Algorithmus?

[Bjoerk]

So. Jetzt wirds langsam. Flippen geht aber Hatch ist nicht so prickelnd.

Jens, weißt du noch nach welchem Verfahren du deine Polygone "vorgedreieckt" hast?

[Jens01]

Zitat:

"vorgedreieckt"

Äh, was meinst Du genau?

[Bjoerk]

Zitat von Jens01:

Zitat:

"vorgedreieckt"

Äh, was meinst Du genau?

Damit könnte man Polygone wunderbar Delaunay like vorbereiten und dann easy weiter verfeinern.
(Siehe TDelaunayTriangulation.Triangulate in #18)