Delaunay-Triangulation

[Jens01]

Ich hab sowas mal vor kurzem gemacht.
Wenn ich mich richtig entsinne, hatte der Code von Bourke das Problem, dass als Ergebnis keine ganzen Dreiecke kamen. Es wurden nur nacheinander Linien gezeichnet. Das brachte mich auf jeden Fall nicht richtig weiter. Also habe ich es mal selbst probiert. Hatte zum Schluß zwar kein schönen Code, aber dann irgendwann ein gutes Ergebnis.

Bei OpenGL gibt es in der GLU-Lib einen 3d-Triangulator. MrMath hat neben seiner guten Matrix-Lib noch diesen Triangulator. Ich habe den aber nie getestet. Der Autor sagte mir mal, dass der Code funktioniert und er den Triangulator für eine Gesichtserkennung nutzt.

[Bjoerk]

Hi Bud,

außerdem ist der Code ja mal so richtig kacke..

Delphi-Quellcode:

  dVertex = record
    x: Double;
    y: Double;
  end;

..

  public
    HowMany: Integer;
    tPoints: Integer
    TargetForm: TForm;
    destructor Destroy;

Gruß
Thomas

[Jens01]

Also noch mal zur Triangulation.
So einfach wie Dejan Vu das meint, war es nicht. Man muß auch irgendwie den Rand/Umriss festlegen. Ebenso die Umrisse der Löcher. Das kann ganz schön verwirrend sein. War es jedenfalls bei mir.
Ich konnte zum Schluß aber Flächenmomente von beliebigen Querschnitten ausrechnen.

Dies Konzept "Fortscheitende Front" von jfheins würde mich auch interessieren. Hatte ich noch nicht gehört. -zumal es sich wie eine rus Militärtaktik anhört-

Wofür brauchst Du das jetzt genau? Finite Elemente?

[Bjoerk]

Ja, FE. Hab zur Zeit nur eine Lösung für Rechtecke. Wie hats du denn die gemeinsamen Lagerpunkte hingekriegt?

[Jens01]

Zitat:

Wie hats du denn die gemeinsamen Lagerpunkte hingekriegt?

Nein, mache kein FE! Brauche das ganze um Löcher in Flächen von Körpern zu schneiden. Z.B. für ein Zapfenloch in einem Pfosten bei der grafischen OpenGL-Darstellung.
Das mit den Flächenmomenten habe ich einfach nur so aus Spaß dazuprogrmmiert, als ich mit dem Programmteil fertig war. War auch kaum noch Arbeit. Ich glaube aber, dass sich Flächenmomente anders einfacher rechnen lassen als mit Triangulationen.

[Bjoerk]

Jens, kannst du mir sagen was der Autor hier macht?

Delphi-Quellcode:

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
  Complete: PComplete;
  Edges: PEdges;
  Nedge: integer;

  //For Super Triangle
  xmin: double;
  xmax: double;
  ymin: double;
  ymax: double;
  xmid: double;
  ymid: double;
  dx: double;
  dy: double;
  dmax: double;

  //General Variables
  i: integer;
  j: integer;
  k: integer;
  ntri: integer;
  xc: double;
  yc: double;
  r: double;
  inc: Boolean;
begin
  //Allocate memory
  GetMem(Complete, sizeof(Complete^));
  GetMem(Edges, sizeof(Edges^));

  //Find the maximum and minimum vertex bounds.
  //This is to allow calculation of the bounding triangle
  xmin := Vertex^[1].X;
  ymin := Vertex^[1].Y;
  xmax := xmin;
  ymax := ymin;
  For i := 2 To nvert do
  begin
    if Vertex^[i].X < xmin then
      xmin := Vertex^[i].X;
    if Vertex^[i].X > xmax then
      xmax := Vertex^[i].X;
    if Vertex^[i].Y < ymin then
      ymin := Vertex^[i].Y;
    if Vertex^[i].Y > ymax then
      ymax := Vertex^[i].Y;
  end;

  dx := xmax - xmin;
  dy := ymax - ymin;
  if dx > dy then
    dmax := dx
  Else
    dmax := dy;

  xmid := Trunc((xmax + xmin) / 2);
  ymid := Trunc((ymax + ymin) / 2);

  //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.

  Vertex^[nvert + 1].X := (xmid - 2 * dmax);
  Vertex^[nvert + 1].Y := (ymid - dmax);
  Vertex^[nvert + 2].X := xmid;
  Vertex^[nvert + 2].Y := (ymid + 2 * dmax);
  Vertex^[nvert + 3].X := (xmid + 2 * dmax);
  Vertex^[nvert + 3].Y := (ymid - dmax);
  Triangle^[1].vv0 := nvert + 1;
  Triangle^[1].vv1 := nvert + 2;
  Triangle^[1].vv2 := nvert + 3;
  Triangle^[1].Precalc := 0;

  Complete[1] := False;
  ntri := 1;

  //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 Complete^[j] <> True then
      begin
        inc := InCircle(Vertex^[i].X, Vertex^[i].Y, Vertex^[Triangle^[j].vv0].X,
          Vertex^[Triangle^[j].vv0].Y, Vertex^[Triangle^[j].vv1].X,
          Vertex^[Triangle^[j].vv1].Y, Vertex^[Triangle^[j].vv2].X,
          Vertex^[Triangle^[j].vv2].Y, xc, yc, r, j);
        //Include this if points are sorted by X
        if (xc + r) < Vertex[i].X then
          complete[j] := True
        Else
          if inc then
          begin
            Edges^[1, Nedge + 1] := Triangle^[j].vv0;
            Edges^[2, Nedge + 1] := Triangle^[j].vv1;
            Edges^[1, Nedge + 2] := Triangle^[j].vv1;
            Edges^[2, Nedge + 2] := Triangle^[j].vv2;
            Edges^[1, Nedge + 3] := Triangle^[j].vv2;
            Edges^[2, Nedge + 3] := Triangle^[j].vv0;
            Nedge := Nedge + 3;
            Triangle^[j].vv0 := Triangle^[ntri].vv0;
            Triangle^[j].vv1 := Triangle^[ntri].vv1;
            Triangle^[j].vv2 := Triangle^[ntri].vv2;
            Triangle^[j].PreCalc := Triangle^[ntri].PreCalc;
            Triangle^[j].xc := Triangle^[ntri].xc;
            Triangle^[j].yc := Triangle^[ntri].yc;
            Triangle^[j].r := Triangle^[ntri].r;
            Triangle^[ntri].PreCalc := 0;
            Complete^[j] := Complete^[ntri];
            j := j - 1;
            ntri := ntri - 1;
          end;
      end;
    until j >= ntri;

    // 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
    begin
      if Not (Edges^[1, j] = 0) And Not (Edges^[2, j] = 0) then
      begin
        For k := j + 1 To Nedge do
        begin
          if Not (Edges^[1, k] = 0) And Not (Edges^[2, k] = 0) then
          begin
            if Edges^[1, j] = Edges^[2, k] then
            begin
              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;
            end;
          end;
        end;
      end;
    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
    begin
      if Not (Edges^[1, j] = 0) And Not (Edges^[2, j] = 0) then
      begin
        ntri := ntri + 1;
        Triangle^[ntri].vv0 := Edges^[1, j];
        Triangle^[ntri].vv1 := Edges^[2, j];
        Triangle^[ntri].vv2 := i;
        Triangle^[ntri].PreCalc := 0;
        Complete^[ntri] := False;
      end;
    end;
  end;

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

  Triangulate := ntri;

  //Free memory
  FreeMem(Complete, sizeof(Complete^));
  FreeMem(Edges, sizeof(Edges^));
end;

[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. Ansonsten wird er da irgendwie den Algo durchlaufen. Das sehe ich auch nicht so auf den ersten Blick.

[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?