Punkte in ein Polygon überführen

[Bjoerk]

Das nette an dem Link oben ist daß dieser Algo konvexe und Konkave Polygone trianguliert..

Deinen Ansatz über Delauny mit dem nachträglichen rauslöschen finde ich aber auch nicht schlecht. Solange es nur Flächen sind. Genial.

Delphi-Quellcode:

procedure TDelaunyTriangulationForm.AssignButtonClick(Sender: TObject);
var
  I, A, B, C: integer;
  P1, P2, P3: TFloatPoint;
  Points: TFloatPoints;
begin
  Points := TFloatPoints.Create;
  try
    Points.Add(FloatPoint(0, 0));
    Points.Add(FloatPoint(3, 0));
    Points.Add(FloatPoint(5, 2));
    Points.Add(FloatPoint(5, 0));
    Points.Add(FloatPoint(8.5, 0));
    Points.Add(FloatPoint(8.5, 4));
    Points.Add(FloatPoint(10, 4));
    Points.Add(FloatPoint(10, 1));
    Points.Add(FloatPoint(12.5, 3.5));
    Points.Add(FloatPoint(11, 6));
    Points.Add(FloatPoint(8.5, 6));
    Points.Add(FloatPoint(5, 8));
    Points.Add(FloatPoint(3, 6));
    Points.Add(FloatPoint(5, 6));
    Points.Add(FloatPoint(4, 4));
    Points.Add(FloatPoint(1.5, 3.5));
    Points.Add(FloatPoint(0, 8));
    Points.Add(FloatPoint(0, 3));
    Points.Add(FloatPoint(1.5, 1.5));
    FMesh.Assign(Points);
    FMesh.Mesh;
    for I := FMesh.TrianglesCount downto 1 do
    begin
      A := FMesh.A[I];
      B := FMesh.B[I];
      C := FMesh.C[I];
      P1 := FloatPoint(FMesh.X[A], FMesh.Y[A]);
      P2 := FloatPoint(FMesh.X[B], FMesh.Y[B]);
      P3 := FloatPoint(FMesh.X[C], FMesh.Y[C]);
      if not Points.PtIn(TriangleIncircleCenter(P1, P2, P3), true) then
        FMesh.DeleteTriangle(I);
    end;
    PaintBox.Invalidate;
  finally
    Points.Free;
  end;
end;

[Bjoerk]

Jens, du hast völlig recht. Das geht mit dem Delauny auch für meine Zwecke.

Delphi-Quellcode:

procedure TPBDelaunyTriangulation.PolygonsMesh(Value: TAreas; const dX, dY: double); // PB = Paul Bourke;
var
  Line, Horz, Vert: TFloatLine;
  I, J, Row, Col: integer;
  P: TFloatPoint;
  MaxX, MinX, MaxY, MinY: double;
  AxisHorz, AxisVert: TFloats;
  Nodes: TFloatPoints;
begin
  // TAreas = Liste von Polygonen; Polygon[Index].Count = 1 .. N;
  // - Polygon[Index].Count = 1; -> Punkt;
  // - Polygon[Index].Count = 2; -> Linie;
  // - Polygon[Index].Count > 2 and Counterclockwise; -> pos. Fläche;
  // - Polygon[Index].Count > 2 and not Counterclockwise; -> neg. Fläche;
  AxisHorz := TFloats.Create;
  AxisVert := TFloats.Create;
  Nodes := TFloatPoints.Create;
  try
    // I von III: Nodes ermitteln;
    for I := 0 to Value.Count - 1 do // Alle Flächen, Punkte und Linien;
      for J := 0 to Value[I].Count - 1 do
      begin
        AxisHorz.Add(Value[I].Items[J].Y);
        AxisVert.Add(Value[I].Items[J].X);
      end;
    AxisHorz.RemoveSameValues; // Sortieren und Doppel rauslöschen;
    AxisVert.RemoveSameValues;
    AxisHorz.Refine(dX); // ggf. Zwischenwerte einfügen, so daß der Abstand..
    AxisVert.Refine(dY); // ..zwischen 2 Achsen <= dX bzw. dY wird;
    MinX := Value.MinX; // Unten/Links;
    MinY := Value.MinY;
    MaxX := Value.MaxX; // Oben/Rechts;
    MaxY := Value.MaxY;
    // Schnittpunkte ermitten;
    // *** Kriterium für TAreas.PtIn:
    // - Wenn in einer pos. Fläche und nicht in einer neg. Fläche;
    // - Polygonlinien können sich berühren:
    // - Value.PtIn führt für pos. Flächen Inflate(+1mm) und für neg. Flächen Inflate(-1mm) durch;
    // - Value.PtIn gibt den Index des Polygons zurück, in dem sich der Punkt befindet;
    // *** Kriterium für Line.Intersect;
    // - Gibt nur Schnittpunkte zurück, die sich innerhalb der Strecken A1A2 und B1B2 befinden;
    for Row := 0 to AxisVert.Count - 1 do
    begin
      Vert.P1 := FloatPoint(AxisVert[Row], MinY);
      Vert.P2 := FloatPoint(AxisVert[Row], MaxY);
      for I := 0 to Value.Count - 1 do
        for J := 0 to Value[I].Count - 1 do
        begin
          Line.P1 := Value[I].Items[J];
          Line.P2 := Value[I].NextItems[J];
          if Line.Intersect(Vert, false, P) then // 1 von 3: Schnittpunkte VertLines / PolygonLines;
            if Value.PtIn(P) > -1 then Nodes.Add(P);
      end;
      for Col := 0 to AxisHorz.Count - 1 do
      begin
        Horz.P1 := FloatPoint(MinX, AxisHorz[Col]);
        Horz.P2 := FloatPoint(MaxX, AxisHorz[Col]);
        for I := 0 to Value.Count - 1 do
          for J := 0 to Value[I].Count - 1 do
          begin
            Line.P1 := Value[I].Items[J];
            Line.P2 := Value[I].NextItems[J];
            if Line.Intersect(Horz, false, P) then // 2 von 3: Schnittpunkte HorzLines / PolygonLines;
              if Value.PtIn(P) > -1 then Nodes.Add(P);
          end;
        if Vert.Intersect(Horz, false, P) then // 3 von 3: Schnittpunkte HorzLines / VertLines;
          if Value.PtIn(P) > -1  then Nodes.Add(P);
      end;
    end;
    Nodes.RemoveSameValues; // Doppel rauslöschen;
    // II von III: Delauny;
    Clear;
    for I := 0 to Nodes.Count - 1 do
      Add(Nodes.X[I], Nodes.Y[I]);
    Mesh;
    // III von III: Nicht vorhandene Knoten und Elemente rauslöschen;
    for I := FTrianglesCount downto 1 do // Delauny ist 1-basiert;
      if Value.PtIn(IncircleCenter[I]) < 0 then
        DeleteTriangle(I);
    for I := FNodesCount downto 1 do // ..
      if not NeedNode[I] then
        DeleteNode(I);
    RefreshCapacity;
    Assert(CheckMesh, 'PolygonsMesh.CheckMesh');
    // Ggf. die Dreiecke in Tetragons überführen -> QuadMesh;
  finally
    AxisHorz.Free;
    AxisVert.Free;
    Nodes.Free;
  end;
end;

[Jens01]

Zitat:

Das nette an dem Link oben ist daß dieser Algo konvexe und Konkave Polygone trianguliert..

Das ist bei meinem Algo auch egal. Er benötigt aber die Außenkontur, damit nicht Punkte verbunden werden, die nicht zusammen gehören.
Bei dem Delaunay von Bourke und vielen anderen, hatte ich immer das Problem, dass da immer nur Linien produziert werden, was zum Zeichnen ausreichend ist, aber für die Weiterverarbeitung in einem Mesh nicht. Dazu braucht man dann richtige Dreiecke.
Ich habe mir den Algo deshalb noch mal neu entwickelt. In meinem Testprogramm sind eigentlich echte Dreiecke zu sehen.

[Jens01]

Sooo, eins noch.

Zitat:

Solange es nur Flächen sind. Genial.

Darüber habe ich noch mal nachgedacht. Mit dem Delaunay lassen sich aber auch Volumen zerlegen. Man sollte beliebige Volumen in Pyramiden aufteilen können.
Im 2d prüft man, ob ein Punkt im Kreis liegt. Im 3d müßte man prüfen, ob ein Punkt in einer aufgespannten Kugel liegt.

[Bjoerk]

Mit dem Bourke hab ich mich auch länger beschäftigt. Der Bourke produziert aber keine Linien sondern echte Dreiecke. Der veröffentlichte Delphi Code ist allerdings die totale Katastrophe. Man muß ihn komplett umschreiben. Bourke macht ja den Trick mit dem Supertriangle und löscht die Dreiecke zu dem am Schluß wieder raus, (genau wie wir mit den Dreiecken, die nicht im Polygon liegen). Anyway..

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