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]

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

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