Punkte in ein Polygon überführen

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