Punkte in ein Polygon überführen

[Bjoerk]

Ich müßte für eine Triangulation Punkte in ein Polygon überführen. Der Algorithmus trianguliert jeweils die Punkte in einer Grid.Cell (Siehe Anlage). Der Algo benötigt die Punkte jeweils als rechtsdrehendes xor linksdrehendes Polygon. Konvex oder Konkav ist (hier) egal. Deshalb wollte ich Sortieren. Es sind in der Regel nur so 3 bis 6 Punkte je Cell. Das klappt aber nicht (immer). Überlegt hatte ich mir das. Klappt aber irgendwie nicht (immer). Jemand eine Idee? Wäre klasse.

Delphi-Quellcode:

function TFloatPoints.Counterclockwise(const Index: integer): boolean;
var
  A, B, C: TFloatPoint;
  x21, x31, y21, y31: double;
begin
  A := PrevItems[Index];
  B := Items[Index];
  C := NextItems[Index];
  x21 := B.X - A.X;
  x31 := C.X - A.X;
  y21 := B.Y - A.Y;
  y31 := C.Y - A.Y;
  Result := x21 * y31 - x31 * y21 > 0;
end;

procedure TFloatPoints.SortCounterclockwise;
var
  I, J: integer;
begin
  for I := 0 to Count - 2 do
    for J := I + 1 to Count - 1 do
      if not Counterclockwise(J) then
        Exchange(J, Next(J));
  for I := 0 to Count - 1 do
    if not Counterclockwise(I) then
      ShowMessage('not Counterclockwise');
end;

[Jens01]

Bist Du sicher, ob in Counterclockwise die Berechnung richtig ist?
PaulBourke hat das auf dieser Seite gemacht. Ich weiss nicht, ob Du die Formel optimiert hast, jedenfalls multipliziert er dort erst und subt dann.

Clockwise = Summe(Xi*Yi+1 - Xi+1*Yi) < 0

Edit:
Die Formel drüber könnte Deine sein. Die ist aber auch etwas anders, oder?

[Bjoerk]

Hi Bud, Summe Det kann man hier nicht anwenden. Hatten wir das nicht schonmal?

[Jens01]

Zitat:

cross product = ((xi - xi-1),(yi - yi-1)) x ((xi+1 - xi),(yi+1 - yi)
= (xi - xi-1) * (yi+1 - yi) - (yi - yi-1) * (xi+1 - xi)
A positive cross product means we have a counterclockwise polygon.

Ist das jetzt nicht Dein Ansatz?

[Bjoerk]

Ja, das wird wohl das gleiche sein. Jedenfalls immer drei aufeinander folgende Punkte und nicht die Summe aller Punkte. Man muß auf >= 0 prüfen (CompareValue(.., 0) >= 0) und statt J mit Next(J) natürlich I und J tauschen. War wohl gestern Abend schon ziemlich fertig. Dieser Mesh Kram macht einen sowieso völlig fertig. Ich denke mal, ich habs jetzt aber. Die Triangulation einer Grid.Cell mache ich übrigens damit. Den Code muß man allerdings komplett umschreiben. Soll ich's mal anhängen? Wollte mich bei dir übrigens auch noch bedanken, der Interessent vom Statikerforum letztens hat vor einer Woche gekauft (ist ein Bauamt).

[Jens01]

Hier ist mein Testprogramm

Moment, geht nicht..

[Jens01]

Ich lade das Programm gleich noch hoch, es gibt da kleine Probleme bei mir.

Die verlinkte Triangulation ist aber keine Delaunay-Triangulation, oder? Da könnte das Problem auch liegen.

[Jens01]

Hier jetzt mein Test für eine Delaunay-Triangulation.

Das Ding ist nur ein Test, den ich zur Entwicklung genutzt habe.
Man benötigt die Aussenkontur des Polygons und wirft Punkte rein. Die Punkte werden dann zu Dreiecken trianguliert. Wenn Löcher im Polygon sind, werden die entsprechenden Dreiecke einfach rausgelöscht.

[Bjoerk]

Fein. Sieht gut aus. Ich hab hier jedoch i.d.R. mehrere Polygone übereinder und/oder die Polygone können sich überlappen. Außerdem kann ein Polygon rechts- oder linksdrehend sein und aus nur einem oder 2 Punkten bestehen (Pfosten auf Decke, Wand auf Decke, Punktlager, Unterzug ect..). Diese Punkte, Linien und Ränder müssen sich im FE Raster wiederfinden. Und der User gibt eine max. Elementgröße vor. Deshalb der Aufwand mit dem Raster.

Beispiel:

Delphi-Quellcode:

var
  Areas: TAreas;
begin
  Areas := TAreas.Create;
  try
    Areas.Count := 2;

    Areas[0].Add(FloatPoint(0, 0));
    Areas[0].Add(FloatPoint(3, 0));
    Areas[0].Add(FloatPoint(5, 2));
    Areas[0].Add(FloatPoint(5, 0));
    Areas[0].Add(FloatPoint(8.5, 0));
    Areas[0].Add(FloatPoint(8.5, 4));
    Areas[0].Add(FloatPoint(10, 4));
    Areas[0].Add(FloatPoint(10, 1));
    Areas[0].Add(FloatPoint(12.5, 3.5));
    Areas[0].Add(FloatPoint(11, 6));
    Areas[0].Add(FloatPoint(8.5, 6));
    Areas[0].Add(FloatPoint(5, 8));
    Areas[0].Add(FloatPoint(3, 6));
    Areas[0].Add(FloatPoint(5, 6));
    Areas[0].Add(FloatPoint(4, 4));
    Areas[0].Add(FloatPoint(1.5, 3.5));
    Areas[0].Add(FloatPoint(0, 8));
    Areas[0].Add(FloatPoint(0, 3));
    Areas[0].Add(FloatPoint(1.5, 1.5));

    Areas[1].Add(FloatPoint(6, 3));
    Areas[1].Add(FloatPoint(7, 4.5));
    Areas[1].Add(FloatPoint(8, 3));

    FPolygonMesh.TriangleMesh(Areas, 1, 1);

    PaintBox.Invalidate;
  finally
    Areas.Free;
  end;

[Jens01]

Mein Programm ist eine reine Triangulation. Dazu benötigt man mehr als 3 Punkte, sonst ist es ja schon ein Dreieck. Alles andere muss vor der Nutzung berücksichtigt werden.

Ich war mir bei Deiner verlinkten Software nicht so sicher, ob das die beste Lösung für eine Triangulation ist.

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

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.

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

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