Delaunay-Triangulation

[Bjoerk]

Ich hab leider gar keinen Plan von diesem Algo. Hat jemand mal was in der Richtung gemacht oder eine Idee wie man sowas angehen könnte?

http://de.wikipedia.org/wiki/Delaunay-Triangulation

[Dejan Vu]

Wo ist denn dein Problem? Ich habe es auch noch nicht gemacht, aber es sieht doch nicht so schwer aus.

Delphi-Quellcode:

Delaunay := TDelaunay.Create;
Delaunay.Initialize();
For P in Points do Delaunay.Add(P);

Initialize ist einfach und was das 'Add' macht, steht da doch auch.

Edit: Ok, ein wenig hochnäsig. Mach mal 'stepwise refinement' bzw. top down programming(*)
also

Delphi-Quellcode:

Procedure TDelaunay.Add (P : TPoint);
var
  i : Integer;
  tmp : TTriangle;

begin
  for i:=0 to Triangles.Length-1 do 
    if Triangles[i].Contains(P) then begin
      tmp := Triangles[i];
      Triangles.Remove(i);
      Triangles.Add(TTriangle.create(tmp.a,p,tmp.b));
      Triangles.Add(TTriangle.create(tmp.b,p,tmp.c));
      Triangles.Add(TTriangle.create(tmp.c,p,tmp.a));
      CheckBoundariesWithFlipAlgorithm;
      exit;
     end;

  raise exception.Create('Point not in triangle');
end;

Ok. Bleibt noch 'CheckBoundariesWithFlipAlgorithm'. Der geht einfach alle Dreiecke durch und prüft die Bedingungen mit diesem Umkreis. Oder nur für die drei Dreiecke, die dazugekommen sind, wäre vielleicht sinniger.

(*) Alles auskodieren, was einem klar ist und alles, was noch nicht klar ist, als Methode aufrufen. Danach die neuen Methoden nach dem gleichen Schema ausformulieren. Klappt meistens, ist gut lesbar und fast schon ein wenig 'clean code'.

[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 :)

[pelzig]

Zitat:

Zitat von Dejan Vu (Beitrag 1271252)

Wo ist denn dein Problem? Ich habe es auch noch nicht gemacht, aber es sieht doch nicht so schwer aus.

Delphi-Quellcode:

Delaunay := TDelaunay.Create;
Delaunay.Initialize();
For P in Points do Delaunay.Add(P);

Initialize ist einfach und was das 'Add' macht, steht da doch auch.

Ich habe zwar keinen Führerschein, aber:

- Auto kaufen
- Losfahren

sollte doch wohl kein Problem in einer Millionenstadt darstellen!

Meint wohl Dejan Vu?

MfG

[Dejan Vu]

Genau so. Wenn Du das nicht verstehst, dein Problem. Wie man fährt, habe ich oben ja schon gezeigt.
Der Algorithmus wird erklärt und man muss ihn nur in Delphi übersetzen. Wenn man die Sprache beherscht, geht das doch recht flott. Wie man sieht.

[pelzig]

Ganz genau so in #2:

Geändert von Dejan Vu (Heute um 17:37 Uhr)

Hinterher ist man immer schlauer und altklüger ... :-D

MfG

[Bjoerk]

Ok. Ich versuch mal bissl näher zu erläutern. Und gleich vorab, so einfach wie sich das Dejan Vu in seinem judenglichen Leichtsinn vorstellt ist es natürlich nicht. Stellen wir uns eine beliebige Wohnung vor. Diese Wohnung hat Zimmer. Wir stehen in einem dieser Zimmer und schauen nach unten auf den Boden. Was wir dann sehen nennt sich der Grundriss dieses Zimmers. Wir bekommen vom Architekten einen Grundriss der gesamten Wohnung. Die Zimmer nennen wir Polygone und verlegen die gesamte Wohnung mit dreiecks- bzw. vierecksförmigen Parkett. Das Parkett nennen wir finite Elemente. Wir sollen das Parkett wie folgt verlegen: Es dürfen beliebige drei- und viereckige Teile verwendet werden, die Teile sollen jedoch in etwa gleich groß sein. Wir sollen nun in jedes Polygon ein Netz reinlegen ("ParkettFloodFill"), so daß jedes Polygonnetz die Randknoten der benachbarten Polygonnetze teilt. Dann gibt es ggf. Treppenlöcher. Die nennen wir Aussparungen. In diesen Aussparungen gibt es kein Netz. In etwa klar? :gruebel:

[jfheins]

Wenn es schon um FEM geht: Das Verfahren "Fortscheitende Front" führt meistens zu schöneren Netzen.
Morgen kann ich dir vielleicht etwas ausführlicher helfen.

[Bjoerk]

Gerade eben entdeckt. Das hat wohl ein netter Mensch nach Delphi übersetzt. Schau ich mir näher an.

[Bjoerk]

Wow jfheins, das sieht mega aus. So hätt ich's gern..

[Jens01]

Ich hab sowas mal vor kurzem gemacht.8-)
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?

:-D 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:

Zitat von Medium (Beitrag 1271255)

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:

Zitat von Jens01 (Beitrag 1271329)

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. :oops: 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? :gruebel:

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

[Dejan Vu]

Zitat:

Zitat von Bjoerk (Beitrag 1271269)

Ok. Ich versuch mal bissl näher zu erläutern...

Hat das dann noch etwas mit der Delaunay-Triangulation zu tun?

[Bjoerk]

Zitat:

Zitat von Jens01 (Beitrag 1271864)

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)

[Bjoerk]

Kaum vergeht ein Jahr, schon bekommt man eine Idee wie man's umsetzen könnte. :-D

Delphi-Quellcode:

function TPolygonTriangulation.AdvancingFrontTriangulate(const RefineCount: integer): boolean;
var
  A, B, C, I: integer; // Dreieck ABC (B = Current)
  D: TFloatPoint;
  Beta: double;
  Triangle: TPolygonTriangle;
  AdvancingFront: TFloatPoints;
begin
  Result := true;
  AdvancingFront := TFloatPoints.Create;
  try
    AdvancingFront.Assign(FPoints);
    AdvancingFront.RefinePolygon(RefineCount);
    FTriangles.Clear;
    B := 0;
    while Result and (AdvancingFront.Count > 3) do
    begin
      A := AdvancingFront.Prev(B);
      C := AdvancingFront.Next(B);
      Triangle.A := AdvancingFront[A];
      Triangle.B := AdvancingFront[B];
      Triangle.C := AdvancingFront[C];
      Beta := Triangle.Beta;
      if CompareFloat(Beta, 0.5 * Pi) < 0 then // Konvex;
      begin
        AdvancingFront.Delete(B);
        FTriangles.Add(Triangle);
      end
      else
        if (CompareFloat(Beta, 0.5 * Pi) >= 0) and (CompareFloat(Beta, Pi) <= 0) then // Konvex;
        begin
          D.X := 0.5 * (Triangle.A.X + Triangle.C.X) + 0.2 * (Triangle.A.Y - Triangle.C.Y);
          D.Y := 0.5 * (Triangle.A.Y + Triangle.C.Y) + 0.2 * (Triangle.C.X - Triangle.A.X);
          I := FTriangles.IndexOfPtIn(D);
          if I < 0 then
          begin
            AdvancingFront[B] := D;
            FTriangles.AddTriangle(Triangle.A, Triangle.B, D);
            FTriangles.AddTriangle(Triangle.B, Triangle.C, D);
          end
          else
          begin
           // ???
           // Result := false;
          end;
        end
        else
        begin
          Result := false; // Konkav; to do ..
          FTriangles.Clear;
        end;
      if Assigned(FOnTriangulate) then
        FOnTriangulate(50, 100);
      B := AdvancingFront.Next(B);
    end;
    if Result then
    begin
      FTriangles.AddTriangle(AdvancingFront[0], AdvancingFront[1], AdvancingFront[2]);
      FTriangles.CircumcircleFlip;
    end;
  finally
    AdvancingFront.Free;
  end;
end;

Hab mal bissl weitergemacht. Macht einen echt fertig.. Ich kriegs jetzt bis soweit hin, bis ein neuer Punkt innerhalb eines bereits bestehenden Dreiecks zu liegen kommt. Aber, wie geht es dann weiter? Jemand eine Idee dazu?

Anlage

[Jens01]

Ich habe mal bei H.Herrmann gelesen, dass das das richtige Konzept für FEM-Netze ist.

[Bjoerk]

Danke für den Link. Der Nelson sucht sich den neuen Punkt nach gewissen Kriterien. Damit hat man aber noch nicht die eigentlich Triangulation (wenn ich's richtig gelesen hab). Hab für dieses Jahr auch schon wieder die Schnauze voll von diesem Thema. :shock: :cyclops: :drunken::wall: Ist wohl nur was für Männer ohne Nerven.. Rufe ggf. mal den Autor dieser Diss an (wohnt grad um die Ecke). LG Thomas

[Jens01]

Also die reine Triangulation ist doch einfach:shock:
Ich habe mir vor einiger Zeit einen fertig gemacht. Unterm Strich ist die finale Lösung ganz einfach und ein überschaubarer Code.
u.a. gibt es hier auch einen Triangulator. Den habe ich aber nicht probiert.

[Bjoerk]

Jens, den Delauny kann ich dir auch schicken..

Stell dir aber eine Geschoßdecke vor. Da gibt es verschiedene Polygone (Zimmer) die aneinander stoßen. Die Übergänge müssen geglättet werden. Dann gibt es i.d.R. ein Treppenloch. Dann kann es vorgegebene Achsen geben (Unterzüge). Dann gibt es ggf. Einzellasten aus Dachpfosten auf der Decke. Dann gibt es evtl. Stützen. Das heiß, die Lage all dieser Knoten ist fest und das Raster muß diese in etwa treffen und das Treppenloch muß ausgespart werden. Usw. usw..

Ich hab das bis jetzt ja auch in meiner Software, nur daß der User ein Raster angeben muß. Es sind nur Koordinaten zulässig, die einen Rasterpunkt darstellen. Und Schrägen werden abgetreppt (da Rechteckselemente). Ich wollte halt zusätzlich eine neues FEM Programm entwickeln (Dreiecke).

Apropos Dreiecke, hast du einen Tipp für eine GUTE Steifigkeitsmatrix für Dreieckselemente?

[Jens01]

Okay, dann habe ich Dich falsch verstanden.

[ikdeveloper]

I cannot find the uUtils, uClasses units for the TDelaunayTriangulation class. Could you help.