Punkte in ein Polygon überführen

[Jens01]

Zitat:

Der Bourke produziert aber keine Linien sondern echte Dreiecke.

Hmm, hatte ich anders in Erinnerung. Wieso habe ich das denn noch mal neu gemacht? Wahrscheinlich, weil der Code so undurchsichtig war....

[Bjoerk]

Japp. Wie gesagt, ich hab den Code völlig umgeschrieben. Ich häng ihn mal an, weil: Da du einen eigenen Delauny geschrieben hast, bist du wohl einer der wenigen Menschen auf diesem Planeten, der mir eventuell sagen könnte, was Bourke in seiner Triangulate treibt? Warum der Umweg über die Ränder und wieso die beiden Hilfsvariablen (hab sie InCircleCalculated und TriangleComplete umbenannt). Nur falls du Zeit und Lust hast..

Ich dachte eigentlich, daß der Delauny sich einfach für jeden Punkt das Dreieck mit dem minimalsten Abstand zu dessen Umkreismittelpunkt sucht und dann den Punkt mit den drei Punkten dieses Dreiecks verbindet, also so die 3 neuen Dreiecke entstehen?

[Bjoerk]

Bin mir ziemlich sicher, daß wir die fünf Hilfsvariablen gar nicht brauchen!?

Delphi-Quellcode:

function TPBDelaunyTriangulation.InCircle(const NodeIndex, TriangleIndex: integer): boolean;
var
  A, B, C: integer;
  xC, yC, m1, m2, mx1, mx2, my1, my2, SqrR1, SqrR2: double;
begin
  A := FA[TriangleIndex];
  B := FB[TriangleIndex];
  C := FC[TriangleIndex];
  if SameValue(FY[B], FY[A]) then
  begin
    m2 := -(FX[C] - FX[B]) / (FY[C] - FY[B]);
    mx2 := (FX[B] + FX[C]) / 2;
    my2 := (FY[B] + FY[C]) / 2;
    xC := (FX[B] + FX[A]) / 2;
    yC := m2 * (xC - mx2) + my2;
  end
  else
    if SameValue(FY[C], FY[B]) then
    begin
      m1 := -(FX[B] - FX[A]) / (FY[B] - FY[A]);
      mx1 := (FX[A] + FX[B]) / 2;
      my1 := (FY[A] + FY[B]) / 2;
      xC := (FX[C] + FX[B]) / 2;
      yC := m1 * (xC - mx1) + my1;
    end
    else
    begin
      m1 := -(FX[B] - FX[A]) / (FY[B] - FY[A]);
      m2 := -(FX[C] - FX[B]) / (FY[C] - FY[B]);
      mx1 := (FX[A] + FX[B]) / 2;
      mx2 := (FX[B] + FX[C]) / 2;
      my1 := (FY[A] + FY[B]) / 2;
      my2 := (FY[B] + FY[C]) / 2;
      if not SameValue(m1 - m2, 0) then
      begin
        xC := (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
        yC := m1 * (xC - mx1) + my1;
      end
      else
      begin
        xC := (FX[A] + FX[B] + FX[C]) / 3;
        yC := (FY[A] + FY[B] + FY[C]) / 3;
      end;
    end;
  SqrR1 := Sqr(FX[NodeIndex] - xC) + Sqr(FY[NodeIndex] - yC);
  SqrR2 := Sqr(FX[B] - xC) + Sqr(FY[B] - yC);
  Result := CompareValue(SqrR1, SqrR2) <= 0; // = PtInCirlc(X, Y, xC, yC, R2);
end;

procedure TPBDelaunyTriangulation.RemoveInvalidEdges;
var
  I, J: integer;
begin
  for I := 1 to FEdgesCount - 1 do
    if (FLeft[I] <> 0) and (FRight[I] <> 0) then
      for J := I + 1 to FEdgesCount do
        if (FLeft[J] <> 0) and (FRight[J] <> 0) then
          if (FLeft[I] = FRight[J]) and (FRight[I] = FLeft[J]) then
          begin
            FLeft[I] := 0;
            FRight[I] := 0;
            FLeft[J] := 0;
            FRight[J] := 0;
          end;
end;

function TPBDelaunyTriangulation.Triangulate: integer;
var
  Triangle, Node, I: integer;
begin
  SetSuperTriangle;
  TriangulateClear;
  Result := 1;
  try
    for Node := 1 to FNodesCount do
    begin
      FEdgesCount := 0;
      Triangle := 0;
      while Triangle < Result do
      begin
        Inc(Triangle);
        if InCircle(Node, Triangle) then
        begin
          FLeft[FEdgesCount + 1] := FA[Triangle];
          FRight[FEdgesCount + 1] := FB[Triangle];
          FLeft[FEdgesCount + 2] := FB[Triangle];
          FRight[FEdgesCount + 2] := FC[Triangle];
          FLeft[FEdgesCount + 3] := FC[Triangle];
          FRight[FEdgesCount + 3] := FA[Triangle];
          Inc(FEdgesCount, 3);
          FA[Triangle] := FA[Result];
          FB[Triangle] := FB[Result];
          FC[Triangle] := FC[Result];
          Dec(Triangle);
          Dec(Result);
        end;
      end;
      RemoveInvalidEdges;
      for I := 1 to FEdgesCount do
        if (FLeft[I] <> 0) and (FRight[I] <> 0) then
        begin
          Inc(Result);
          FA[Result] := FLeft[I];
          FB[Result] := FRight[I];
          FC[Result] := Node;
        end;
    end;
  finally
    FTrianglesCount := DeleteSuperTriangle(Result);
  end;
end;

[Bjoerk]

So. Jetzt. Ja. Hab bei Bourke auf der Hompage mal gesucht. Und siehe da, dort steht völlig klar und eindeutig, wie er vorgeht.

Delphi-Quellcode:

procedure TPBDelaunyTriangulation.Triangulate;
var
  Triangle, Node, I: integer;
  Left, Right: TPBIntegers; // Edges;
begin
  // input : vertex list;
  // output : triangle list;
  //    initialize the triangle list;
  //    determine the supertriangle;
  //    add supertriangle vertices to the end of the vertex list;
  //    add the supertriangle to the triangle list;
  //    for each sample point in the vertex list;
  //       initialize the edge buffer;
  //       for each triangle currently in the triangle list;
  //          calculate the triangle circumcircle center and radius;
  //          if the point lies in the triangle circumcircle then;
  //             add the three triangle edges to the edge buffer;
  //             remove the triangle from the triangle list;
  //          endif;
  //       endfor;
  //       delete all doubly specified edges from the edge buffer;
  //          this leaves the edges of the enclosing polygon only;
  //       add to the triangle list all triangles formed between the point;
  //          and the edges of the enclosing polygon;
  //    endfor;
  //    remove any triangles from the triangle list that use the supertriangle vertices;
  //    remove the supertriangle vertices from the vertex list;
  // end;
  Left := TPBIntegers.Create;
  Right := TPBIntegers.Create;
  try
    AddSuperTriangle;
    for Node := 1 to NodesCount - 3 do
    begin
      Left.Clear;
      Right.Clear;
      for Triangle := TrianglesCount downto 1 do
      begin
        if InCircle(Node, Triangle) then
        begin
          Left.Add(A[Triangle]);
          Right.Add(B[Triangle]);
          Left.Add(B[Triangle]);
          Right.Add(C[Triangle]);
          Left.Add(C[Triangle]);
          Right.Add(A[Triangle]);
          FA.Delete(Triangle);
          FB.Delete(Triangle);
          FC.Delete(Triangle);
        end;
      end;
      DeleteInvalidEdges(Left, Right);
      for I := 1 to Left.Count do
      begin
        FA.Add(Left[I]);
        FB.Add(Right[I]);
        FC.Add(Node);
      end;
    end;
    DeleteSuperTriangle;
  finally
    Left.Free;
    Right.Free;
  end;
end;

Was mich an dem Code doch sehr störte, war diese katastrophale Triangulate, so daß man kein Items.Add ect. verwenden konnte, obwohl der Algo sozusagen ein Paradebeispiel für Listen ist. Das 1 basierte hab ich aber gelassen, stört mich nicht weiter (hab ja extra Listen dazu geschrieben). Werd die unit auch noch der Bourke Gemeinde rüberschicken. Hab noch eine kleine grafische Klasse ergänzt. Bitte den oben geposteten Code nicht mehr verwenden, sondern ggf. diesen. Schöne Pfingsten. LG Thomas

[Jens01]

Zitat:

Da du einen eigenen Delauny geschrieben hast, bist du wohl einer der wenigen Menschen auf diesem Planeten, der mir eventuell sagen könnte, was Bourke in seiner Triangulate treibt?

Eigentlich ist das ganz einfach.

Wenn ich richtig erinnere:
Du nimmst Deine Punktwolke verbindest 3 Punkte, schlägst einen Kreis um die Punkte, so dass alle Punkte auf dem Kreis liegen und überprüfst, dass sonst kein anderer Punkt innerhalb des Kreises liegt.

Von der Taktik her, musst Du halt so vorgehen, dass Du mit den 2 Punkten beginnst, die am dichtesten zusammen liegen. Von den 2 Punkten aus nimmst Du Dir systematisch (per Schleife) alle anderen Punkte vor, schlägst einen Kreis und prüfst, ob noch ein Punkt drin liegt. Wenn kein weiterer Punkt drin liegt, machst Du die 3 Punkte zum Dreieck.
Von den jeweiligen Seiten dieses Dreiecks arbeitest Du weiter.

[Bjoerk]

Jens, vielen Dank für deine Mühe. Ich habs jetzt nach Bourke (Siehe Code). Dann ist es vergleichsweise einfach. Kann man so 1 zu 1 programmieren. In Delphi, Siehe #19 oder angehängte unit dort.

Code:

subroutine triangulate
input : vertex list
output : triangle list
   initialize the triangle list
   determine the supertriangle
   add supertriangle vertices to the end of the vertex list
   add the supertriangle to the triangle list
   for each sample point in the vertex list
      initialize the edge buffer
      for each triangle currently in the triangle list
         calculate the triangle circumcircle center and radius
         if the point lies in the triangle circumcircle then
            add the three triangle edges to the edge buffer
            remove the triangle from the triangle list
         endif
      endfor
      delete all doubly specified edges from the edge buffer
         this leaves the edges of the enclosing polygon only
      add to the triangle list all triangles formed between the point 
         and the edges of the enclosing polygon
   endfor
   remove any triangles from the triangle list that use the supertriangle vertices
   remove the supertriangle vertices from the vertex list
end