Schnelle Triangulation und konvexe Hülle in 2D
Re: Schnelle Triangulation und konvexe Hülle in 2D
Ja danke. Sollte ich so oder so optional von außerhalb anwenden, teils ist der Effekt auch gewünscht.
- NicTheQuick
- Ein Admin
- Beiträge: 8679
- Registriert: 29.08.2004 20:20
- Computerausstattung: Ryzen 7 5800X, 32 GB DDR4-3200
Ubuntu 22.04.3 LTS
GeForce RTX 3080 Ti - Wohnort: Saarbrücken
- Kontaktdaten:
Re: Schnelle Triangulation und konvexe Hülle in 2D
Falls du es noch nicht gesehen hast, gibt es jetzt im ersten Post auch eine Version als Modul, dann gibt es auch keine Probleme mehr mit gleichnamigen Structures.
- NicTheQuick
- Ein Admin
- Beiträge: 8679
- Registriert: 29.08.2004 20:20
- Computerausstattung: Ryzen 7 5800X, 32 GB DDR4-3200
Ubuntu 22.04.3 LTS
GeForce RTX 3080 Ti - Wohnort: Saarbrücken
- Kontaktdaten:
Re: Schnelle Triangulation und konvexe Hülle in 2D
Ich poste mal für alle, die interessiert sind, meine aktuelle Debug-Version.
Die Funktion CreateRandomizedPointCloud() erstellt momentan eine etwas schräge, aber regelmäßige Punktwolke. Damit hatte ich die letzte Zeit die meisten Probleme. Ihr könnt damit herumspielen wie ihr wollt.
Startet man das Programm, ist die Punktwolkte zu sehen. Mit Linksklick und Drag'n'Drop kann man die Punkte verschieben.
Mit Rechtsklick kann man die Triangulation starten. Nach dem ersten Klick sieht man die initiale konvexe Hülle in Form eines Dreiecks, die momentan immer ab Punkt 0 erstellt wird.
Mit weiteren Rechtsklicks oder der Taste f kann man im Algorithmus vorwärts gehen. Dabei poppen auch ein paar Ausgaben im Debugfenster auf. Diese zeigen an, ob das gerade eben markierte Dreieck die konvexe Hülle erweitert (ADD) oder nicht (not), dahinter stehen die IDs der Punkt, die neu hinzugefügt wurden. Manchmal scheint man dieses Dreieck nicht zu sehen, das liegt aber daran, dass dann zufällig gerade alle drei Punkte auf einer Linie liegen. Und das ist auch das Problem, mit dem ich zu kämpfen hatte. Zwischendurch sieht man auch ein paar sehr große Zahlen. Das ist genau dann der Fall, wenn die Dreieckspunkte auf einer Linie liegen. Das versuche ich abzufangen. In ganz ungünstigen Situationen kann die Triangulation aber auch gar nicht abgeschlossen werden. Das habe ich noch nicht abgefangen, aber vielleicht findet ihr ja so einen Fall. Konstruieren kann man ihn auf jeden Fall.
Wer übrigens anfängt Punkte zu verschieben während die Triangulation läuft, kann sie sich damit kaputt machen. Einen Abbruch-Button gibt es nicht. Man muss dann einmal bis zum Ende durchlaufen (f gedrückt halten) oder das Programm neu starten. Wie gesagt, das ist eh alles nur zum Debuggen, vor allem auch das ganze grafische Zeugs. Aber ich dachte das könnte euch gefallen.
Die Funktion CreateRandomizedPointCloud() erstellt momentan eine etwas schräge, aber regelmäßige Punktwolke. Damit hatte ich die letzte Zeit die meisten Probleme. Ihr könnt damit herumspielen wie ihr wollt.
Startet man das Programm, ist die Punktwolkte zu sehen. Mit Linksklick und Drag'n'Drop kann man die Punkte verschieben.
Mit Rechtsklick kann man die Triangulation starten. Nach dem ersten Klick sieht man die initiale konvexe Hülle in Form eines Dreiecks, die momentan immer ab Punkt 0 erstellt wird.
Mit weiteren Rechtsklicks oder der Taste f kann man im Algorithmus vorwärts gehen. Dabei poppen auch ein paar Ausgaben im Debugfenster auf. Diese zeigen an, ob das gerade eben markierte Dreieck die konvexe Hülle erweitert (ADD) oder nicht (not), dahinter stehen die IDs der Punkt, die neu hinzugefügt wurden. Manchmal scheint man dieses Dreieck nicht zu sehen, das liegt aber daran, dass dann zufällig gerade alle drei Punkte auf einer Linie liegen. Und das ist auch das Problem, mit dem ich zu kämpfen hatte. Zwischendurch sieht man auch ein paar sehr große Zahlen. Das ist genau dann der Fall, wenn die Dreieckspunkte auf einer Linie liegen. Das versuche ich abzufangen. In ganz ungünstigen Situationen kann die Triangulation aber auch gar nicht abgeschlossen werden. Das habe ich noch nicht abgefangen, aber vielleicht findet ihr ja so einen Fall. Konstruieren kann man ihn auf jeden Fall.
Wer übrigens anfängt Punkte zu verschieben während die Triangulation läuft, kann sie sich damit kaputt machen. Einen Abbruch-Button gibt es nicht. Man muss dann einmal bis zum Ende durchlaufen (f gedrückt halten) oder das Programm neu starten. Wie gesagt, das ist eh alles nur zum Debuggen, vor allem auch das ganze grafische Zeugs. Aber ich dachte das könnte euch gefallen.
Code: Alles auswählen
CompilerIf Not #PB_Compiler_Thread
CompilerError "Threadsicherer Modus muss aktiviert sein."
CompilerEndIf
#WIDTH = 1400
#HEIGHT = 800
;-- START OF TRIANGULATION STRUCTURES AND FUNCTIONS
DeclareModule Triangulation
EnableExplicit
#PRECISION = #PB_Float ;oder #PB_Double
#MAKE_DELAUNAY = #True
#IGNORE_SECOND_SORTING = #False
#MAX_FLIPPING = 0
#DEBUG = #True
CompilerIf #PRECISION = #PB_Float
Macro prec
f
EndMacro
#EPSILON = 0.00001
CompilerElseIf #PRECISION = #PB_Double
Macro prec
d
EndMacro
#EPSILON = 0.00000001
CompilerElse
CompilerError "Precision not supported."
CompilerEndIf
Structure Point2D
x.prec
y.prec
id.i
EndStructure
Structure Vector2D Extends Point2D
EndStructure
; Eine Kante merkt sich seine Endpunkte und die maximal zwei Dreicke, zu denen sie gehört
Structure Edge2D
*p.Point2DPC[2]
*t.Triangle2D[2]
mark.i ; Nur für makeDelauney. #True, wenn die Kante in der ToFlip-Liste ist.
flipped.i ; Nur für makeDelauney. Zähler wie häufig eine Kante geflipt wurde.
EndStructure
; Ein Punkt merkt sich natürlich seine Koordinaten und die Kanten und Dreiecken, zu denen er gehört.
Structure Point2DPC Extends Point2D
List *edges.Edge2D()
EndStructure
Structure CircumCircle
center.Point2D
;radius.prec
radiusSq.prec
EndStructure
; Ein Dreieck merkt sich seine drei Eckpunkte und seine drei Kanten
Structure Triangle2D
*p.Point2DPC[3]
*e.Edge2D[3]
cc.CircumCircle
EndStructure
Structure Point2DDiff
*p.Point2D
diff.prec
EndStructure
Structure ConvexHullPoint2D
*p.Point2DPC
used.i
EndStructure
Structure Triangulation
n.i
time.i
Array points.Point2DPC(1)
List edges.Edge2D()
List triangles.Triangle2D()
List convexHull.ConvexHullPoint2D()
nextTriangle.Triangle2D[2]
EndStructure
Prototype StepCallback(*pc.Triangulation, *userData)
Declare clearTriangulation(*pc.Triangulation)
Declare getCircumCircle(*a.Point2D, *b.Point2D, *c.Point2D, *cc.CircumCircle)
Declare.i getPoint(*pc.Triangulation, x.prec, y.prec)
Declare setPoint(*pc.Triangulation, i.i, x.prec, y.prec)
Declare.i isRightHand(*a.Point2D, *b.Point2D, *c.Point2D)
Declare.i isEdgeDelaunay(*edge.Edge2D)
Declare.i isTriangleDelaunay(*triangle.Triangle2D)
Declare triangulate(*pc.Triangulation, seed.i = 0, *cb.StepCallback = 0, *userData = 0)
EndDeclareModule
Module Triangulation
Procedure.s strT(*t.Triangle2D)
ProcedureReturn "T" + *t\p[0]\id + "-" + *t\p[1]\id + "-" + *t\p[2]\id
EndProcedure
Procedure clearTriangulation(*pc.Triangulation) ;correct
Protected i.i
ClearList(*pc\triangles())
ClearList(*pc\edges())
ClearList(*pc\convexHull())
For i = 0 To *pc\n - 1
ClearList(*pc\points(i)\edges())
Next
For i = 0 To 5
*pc\nextTriangle[i % 2]\p[i / 2] = 0
Next
EndProcedure
Procedure.i getPoint(*pc.Triangulation, x.prec, y.prec) ;correct
Protected i.i, diffSq.prec = Infinity(), actualDiffSq.prec
Protected nearest.i = -1
With *pc
For i = 0 To \n - 1
actualDiffSq = Pow(\points(i)\x - x, 2) + Pow(\points(i)\y - y, 2)
If (actualDiffSq < diffSq)
nearest = i
diffSq = actualDiffSq
EndIf
Next
EndWith
ProcedureReturn nearest
EndProcedure
Procedure setPoint(*pc.Triangulation, i.i, x.prec, y.prec) ;correct
If (i >= 0 And i < *pc\n)
*pc\points(i)\x = x
*pc\points(i)\y = y
EndIf
EndProcedure
Procedure getCircumCircle(*a.Point2D, *b.Point2D, *c.Point2D, *cc.CircumCircle) ;correct
Protected i.i
For i = 0 To 1
Protected m.Point2D
m\x = 0.5 * (*b\x - *c\x)
m\y = 0.5 * (*b\y - *c\y)
Protected v1.Vector2D
v1\x = *b\x - *a\x
v1\y = *b\y - *a\y
Protected v2.Vector2D
v2\x = *c\x - *a\x
v2\y = *c\y - *a\y
Protected lambda.prec = NaN()
Protected t.prec = v2\y * v1\x - v1\y * v2\x
If (Abs(t) > #EPSILON)
lambda = (v2\x * m\x + v2\y * m\y) / t
Break
EndIf
Swap *a, *b
Next
*cc\center\x = lambda * v1\y + 0.5 * (*a\x + *b\x)
*cc\center\y = - lambda * v1\x + 0.5 * (*a\y + *b\y)
*cc\radiusSq = Pow(*cc\center\x - *a\x, 2) + Pow(*cc\center\y - *a\y, 2)
If *cc\radiusSq >= 1.0 / (#EPSILON * #EPSILON)
Debug *cc\radiusSq
EndIf
ProcedureReturn Bool(*cc\radiusSq < 1.0 / (#EPSILON * #EPSILON))
EndProcedure
Procedure.i isRightHand(*a.Point2D, *b.Point2D, *c.Point2D) ;correct
Protected ab.Vector2D, bc.Vector2D
ab\x = *b\x - *a\x
ab\y = *b\y - *a\y
bc\x = *c\x - *b\x
bc\y = *c\y - *b\y
Protected z.prec = ab\x * bc\y - ab\y * bc\x
ProcedureReturn Bool(z < -#EPSILON)
EndProcedure
Procedure orderRightHand(*a.Point2D, *p_b.Integer, *p_c.Integer) ;correct
If (Not isRightHand(*a, *p_b\i, *p_c\i))
Swap *p_b\i, *p_c\i
EndIf
EndProcedure
Procedure.i addEdge(*pc.Triangulation, *a.Point2DPC, *b.Point2DPC) ;correct
;Zwei identische Punkte ergeben keine Linie
If (*a = *b)
ProcedureReturn #False
EndIf
;Existiert bereits eine Linie mit diesen beiden Punkten?
ForEach *a\edges()
;If (*a\edges()\p[0] = *b Or *a\edges()\p[1] = *b) ;changed
If (*a\edges()\p[0] = *a And *a\edges()\p[1] = *b) XOr (*a\edges()\p[0] = *b And *a\edges()\p[1] = *a)
ProcedureReturn *a\edges()
EndIf
Next
;Füge die Linie hinzu und speichere ihre Referenz in den Punkten
If AddElement(*pc\edges())
*pc\edges()\p[0] = *a
*pc\edges()\p[1] = *b
*pc\edges()\t[0] = 0
*pc\edges()\t[1] = 0
AddElement(*a\edges())
*a\edges() = @*pc\edges()
AddElement(*b\edges())
*b\edges() = @*pc\edges()
ProcedureReturn @*pc\edges()
EndIf
ProcedureReturn #False
EndProcedure
Procedure.i addTriangle(*pc.Triangulation, *a.Point2DPC, *b.Point2DPC, *c.Point2DPC, force.i = #False) ;correct
;Zwei identische Punkte ergeben kein Dreieck
If (*a = *b Or *b = *c Or *a = *c)
ProcedureReturn #False
EndIf
;Make sure the triangle's points are ordered in right hand order
orderRightHand(*a, @*b, @*c)
Protected cc.CircumCircle
If Not getCircumCircle(*a, *b, *c, cc)
If Not force
ProcedureReturn #False
EndIf
EndIf
Protected Dim *l.Edge2D(2)
*l(0) = addEdge(*pc, *a, *b)
*l(1) = addEdge(*pc, *b, *c)
*l(2) = addEdge(*pc, *c, *a)
;Prüfe, ob die drei Linien zufällig schon ein Dreieck bilden.
;Falls ja, dann gib einfach das zurück.
Protected j.i, k.i, *triangle.Triangle2D, used.i
For j = 0 To 1
used = 0
*triangle = *l(0)\t[j]
If (*triangle)
For k = 0 To 1
used + Bool(*l(1)\t[k] = *triangle)
used + Bool(*l(2)\t[k] = *triangle)
Next
If (used = 2)
Debug "Doppeltes Dreieck?"
ProcedureReturn *triangle
EndIf
EndIf
Next
If (AddElement(*pc\triangles()))
With *pc\triangles()
For j = 0 To 2
\e[j] = *l(j)
\e[j]\t[Bool(\e[j]\t[0])] = @*pc\triangles()
Next
\p[0] = *a
\p[1] = *b
\p[2] = *c
\cc = cc
EndWith
ProcedureReturn @*pc\triangles()
EndIf
ProcedureReturn #False
EndProcedure
Procedure.i isEdgeDelaunay(*edge.Edge2D) ;correct
Protected i.i, j.i, *p.Point2D
With *edge
If (\t[1])
For i = 0 To 1
For j = 0 To 2 ;Iterate over points of other triangle
*p = \t[1 - i]\p[j]
If (*p <> \p[0] And *p <> \p[1]) ;Is the actual point not belonging to the actual edge?
; Is this point from the other triangle within the circum circle of the actual triangle?
If (Pow(*p\x - \t[i]\cc\center\x, 2) + Pow(*p\y - \t[i]\cc\center\y, 2) - \t[i]\cc\radiusSq < -#EPSILON)
ProcedureReturn #False
EndIf
EndIf
Next
Next
EndIf
EndWith
ProcedureReturn #True
EndProcedure
Procedure.i isTriangleDelaunay(*triangle.Triangle2D) ;correct
Protected cc.CircumCircle, i.i, j.i, k.i, *edge.Edge2D, *p.Point2D
For k = 0 To 2
If (Not isEdgeDelaunay(*triangle\e[k]))
ProcedureReturn #False
EndIf
Next
ProcedureReturn #True
EndProcedure
Procedure makeDelaunay(*pc.Triangulation, *cb.StepCallback = 0, *userData = 0) ;correct
Protected NewList *ndEdges.Edge2D()
Protected Dim cc.CircumCircle(1), *p.Point2D
Protected Dim p_i.i(1) ;p_i(x) : Index to point of triangle x which not belongs to the actual edge
Protected i.i, j.i, doFlip.i
Protected *actualEdge.Edge2D, *newEdge.Edge2D
; Iteriere über alle Kanten
ForEach *pc\edges()
; Gehört die Kante zu zwei Dreiecken, prüfe, ob sie Delauney ist
If (*pc\edges()\t[1] And (Not isEdgeDelaunay(*pc\edges())))
; Wenn sie es nicht ist, füge sie zu der Prüfliste hinzu
If (AddElement(*ndEdges()))
*ndEdges() = @*pc\edges()
EndIf
; Und markiere sie
*pc\edges()\mark = #True
Else
*pc\edges()\mark = #False
EndIf
*pc\edges()\flipped = 0
Next
If *cb
*cb(*pc, *userData)
EndIf
While FirstElement(*ndEdges())
*actualEdge = *ndEdges()
DeleteElement(*ndEdges())
With *actualEdge
\mark = #False
doFlip = #False
For i = 0 To 1
For j = 0 To 2 ;Iterate over points of other triangle
*p = \t[1 - i]\p[j]
If (*p <> \p[0] And *p <> \p[1]) ;Is the actual point not belonging to the actual edge?
p_i(1 - i) = j
; Is this point from the other triangle within the circum circle of the actual triangle?
If Pow(*p\x - \t[i]\cc\center\x, 2) + Pow(*p\y - \t[i]\cc\center\y, 2) < \t[i]\cc\radiusSq - #EPSILON * #EPSILON
doFlip = #True
EndIf
EndIf
Next
Next
If (doFlip)
CopyStructure(\t[0], *pc\nextTriangle[0], Triangle2D)
CopyStructure(\t[1], *pc\nextTriangle[1], Triangle2D)
If *cb
*cb(*pc, *userData)
EndIf
CompilerIf #DEBUG
Protected error.i = #False
;Be sure p_i is correct
For i = 0 To 1
If (\t[i]\p[p_i(i)] = \p[0] Or \t[i]\p[p_i(i)] = \p[1])
Debug "p_i(" + i + ") ist nicht korrekt! (Points) [" + \t[i] + "]"
error = #True
EndIf
If (\t[i]\e[(p_i(i) + 1) % 3] <> *actualEdge)
Debug "p_i(" + i + ") ist nicht korrekt! (Edges) [" + \t[i] + "]"
error = #True
EndIf
If (Not isRightHand(\t[i]\p[0], \t[i]\p[1], \t[i]\p[2]))
Debug "Triangle " + strT(\t[i]) + " is not in the right order! [" + \t[i] + "]"
error = #True
EndIf
Next
If (error)
Debug "Error on " + *actualEdge\p[0]\id + "-" + *actualEdge\p[1]\id
ProcedureReturn #False
EndIf
Debug "flip on " + *actualEdge\p[0]\id + "-" + *actualEdge\p[1]\id + " With [" + strT(\t[0]) + "] And [" + strT(\t[1]) + "]"
CompilerEndIf
;Delete Edge from Points-Array
For i = 0 To 1
ForEach \p[i]\edges()
If (\p[i]\edges() = *actualEdge)
DeleteElement(\p[i]\edges())
Break
EndIf
Next
Next
;Swap points to create new triangles.
;This loop runs without error if the triangle's points are ordered in right hand order
;and the edge's order correlates to the points.
For i = 0 To 1
;Give actual edge the new coordinates
\p[i] = \t[i]\p[p_i(i)]
;Make the edge known to the point
AddElement(\p[i]\edges())
\p[i]\edges() = *actualEdge
;Change third point of triangle i to first point of triangle 1 - i
\t[i]\p[(p_i(i) + 2) % 3] = \t[1 - i]\p[p_i(1 - i)]
;Change second edge of triangle i to third edge of triangle 1 - i
\t[i]\e[(p_i(i) + 1) % 3] = \t[1 - i]\e[(p_i(1 - i) + 2) % 3]
;Correct neighbours of new second edge of triangle i
Protected *e.Edge2D
*e = \t[i]\e[(p_i(i) + 1) % 3]
For j = 0 To 1
If (*e\t[j] = \t[1 - i])
*e\t[j] = \t[i]
EndIf
Next
Next
For i = 0 To 1
;Change third edge of triangle i to actual edge
\t[i]\e[(p_i(i) + 2) % 3] = *actualEdge
Next
CompilerIf #DEBUG
If (addEdge(*pc, \t[0]\p[p_i(0)], \t[1]\p[p_i(1)]) <> *actualEdge)
Debug "actualEdge problem!"
EndIf
For i = 0 To 1
If (Not isRightHand(\t[i]\p[0], \t[i]\p[1], \t[i]\p[2]))
Debug "Triangle " + i + " is no more in the right order!"
EndIf
Next
CompilerEndIf
\flipped + 1
;Edge is now Delauney. Add adjacent Edges to List.
LastElement(*ndEdges())
For i = 0 To 1
For j = 0 To 1
*newEdge = \t[i]\e[(p_i(i) + j) % 3]
If ((Not *newEdge\mark) And *newEdge\t[1] And (#MAX_FLIPPING = 0 Or *newEdge\flipped < #MAX_FLIPPING))
If (AddElement(*ndEdges()))
*ndEdges() = *newEdge
*newEdge\mark = #True
EndIf
EndIf
Next
Next
; Update CircumCirles
For i = 0 To 1
If Not getCircumCircle(\t[i]\p[0], \t[i]\p[1], \t[i]\p[2], \t[i]\cc)
Debug "oh oh"
EndIf
Next
EndIf
EndWith
Wend
*pc\nextTriangle[0]\p[0] = 0
*pc\nextTriangle[0]\p[1] = 0
*pc\nextTriangle[0]\p[2] = 0
*pc\nextTriangle[1]\p[0] = 0
*pc\nextTriangle[1]\p[1] = 0
*pc\nextTriangle[1]\p[2] = 0
If *cb
*cb(*pc, *userData)
EndIf
EndProcedure
Procedure triangulate(*pc.Triangulation, seed.i = 0, *cb.StepCallback = 0, *userData = 0) ;correct
Protected Dim sortedPoints.Point2DDiff(*pc\n - 1)
Protected i.i
With *pc
\time = ElapsedMilliseconds()
If (\n < 3)
ProcedureReturn #False
EndIf
;1. Select a seed point x_0 from x_i.
If (seed < 0 Or seed >= \n)
ProcedureReturn #False
EndIf
Protected *x0.Point2DPC = @\points(seed)
clearTriangulation(*pc)
;TODO Hier würde eine lineare Suche reichen, da in Schrift 5 sowieso wieder sortiert wird
;2. Sort according to |x_i - x_0|^2.
For i = 0 To \n - 1
sortedPoints(i)\p = @\points(i)
sortedPoints(i)\diff = Pow(*x0\x - \points(i)\x, 2) + Pow(*x0\y - \points(i)\y, 2)
Next
SortStructuredArray(sortedPoints(), #PB_Sort_Ascending, OffsetOf(Point2DDiff\diff), #PRECISION)
;3. Find the point x_j closest to x_0.
Protected *xj.Point2DPC = sortedPoints(1)\p
;4. Find the point x_k that creates the smallest circum-circle
; with x_0 and x_j and record the center of the circum-circle C.
Protected bestIndex.i
Protected C.CircumCircle, bestC.CircumCircle
bestC\radiusSq = Infinity()
For i = 2 To \n - 1
getCircumCircle(*x0, *xj, sortedPoints(i)\p, @C)
If (C\radiusSq < bestC\radiusSq)
bestIndex = i
bestC = C
EndIf
; Wir suchen nur solange bis der (sortedPoints(i)-x_0)^2 < (bestC\radius * 2) ^ 2 ist
If sortedPoints(i)\diff > bestC\radiusSq * 4
Break
EndIf
Next
;5. Resort the remaining points according to |x_i - C|^2 to
; give points s_i.
If (bestIndex <> 2)
Swap sortedPoints(2)\p, sortedPoints(bestIndex)\p
EndIf
CompilerIf Not #IGNORE_SECOND_SORTING
For i = 3 To \n - 1
sortedPoints(i)\diff = Pow(bestC\center\x - sortedPoints(i)\p\x, 2) + Pow(bestC\center\y - sortedPoints(i)\p\y, 2)
Next
SortStructuredArray(sortedPoints(), #PB_Sort_Ascending, OffsetOf(Point2DDiff\diff), #PRECISION, 3, \n - 1)
CompilerEndIf
;6. Order point x_0, x_j, x_k to give a right handed system.
; This is the initial seed convex hull.
Protected *xk.Point2DPC = sortedPoints(2)\p
orderRightHand(*x0, @*xk, @*xj)
;7. Sequentially add the points s_i to the propagating 2D convex
; hull that is seeded with the triangle formed from x_0, x_j, x_k.
; As a new point is added the facets of the 2D-hull that are visible
; to it form new triangles.
addTriangle(*pc, *x0, *xk, *xj)
ClearList(\convexHull())
AddElement(\convexHull()) : \convexHull()\p = *x0
AddElement(\convexHull()) : \convexHull()\p = *xk
AddElement(\convexHull()) : \convexHull()\p = *xj
If *cb
*cb(*pc, *userData)
EndIf
Protected *p.Point2D, *last.ConvexHullPoint2D, alreadyAdded.i, convexHullSize.i, j.i
Protected *cEdges, *cTriangles, *cConvexHull
;Protected *bestLast.ConvexHullPoint2D, *best.ConvexHullPoint2D, bestRadius.prec = Infinity(), cc.CircumCircle
For i = 3 To \n - 1
*p = sortedPoints(i)\p
Debug "Punkt " + *p\id
alreadyAdded = #False
convexHullSize = ListSize(\convexHull())
FirstElement(\convexHull())
*last = @\convexHull()
NextElement(\convexHull())
For j = 0 To convexHullSize - 1
\nextTriangle[0]\p[0] = *last\p
\nextTriangle[0]\p[1] = *p
\nextTriangle[0]\p[2] = \convexHull()\p
If *cb
*cConvexHull = @\convexHull()
*cEdges = @\edges()
*cTriangles = @\triangles()
*cb(*pc, *userData)
ChangeCurrentElement(\convexHull(), *cConvexHull)
ChangeCurrentElement(\edges(), *cEdges)
ChangeCurrentElement(\triangles(), *cTriangles)
EndIf
If (isRightHand(*last\p, *p, \convexHull()\p))
If addTriangle(*pc, *last\p, *p, \convexHull()\p)
Debug " + ADD " + *last\p\id + "-" + *p\id + "-" + \convexHull()\p\id
\convexHull()\used + 1
*last\used + 1
If (Not alreadyAdded)
InsertElement(\convexHull())
\convexHull()\p = *p
NextElement(\convexHull())
convexHullSize + 1
alreadyAdded = #True
EndIf
Else
Debug " - not"
EndIf
Else
Debug " - not " + *last\p\id + "-" + *p\id + "-" + \convexHull()\p\id
EndIf
*last = @\convexHull()
If (Not NextElement(\convexHull()))
FirstElement(\convexHull())
EndIf
Next
;TODO vielleicht muss das schon viel früher ausgeführt werden.
ForEach \convexHull()
If (\convexHull()\used >= 2)
DeleteElement(\convexHull(), 1)
Else
\convexHull()\used = 0
EndIf
Next
Next
\nextTriangle[0]\p[0] = 0
\nextTriangle[0]\p[1] = 0
\nextTriangle[0]\p[2] = 0
;8. A non-overlapping triangulation of the set of points is created.
; (This is an extremely fast method for creating an non-overlapping
; triangualtion of a 2D point set).
;9: Adjacent pairs of triangles of this triangulation must be 'flipped'
; to create a Delaunay triangulation from the initial non-overlapping
; triangulation.
CompilerIf #MAKE_DELAUNAY
makeDelaunay(*pc, *cb, *userData)
CompilerEndIf
\time = ElapsedMilliseconds() - \time
EndWith
EndProcedure
EndModule
;-- END OF TRIANGULATION FUNCTIONS
Structure Window
id.i
width.i
height.i
title.s
canvasId.i
*pc.Triangulation::Triangulation
clicked.i
leftDown.i
rightDown.i
mutex.i
thread.i
newThread.i
EndStructure
#MAX_INTEGER = 1 << (SizeOf(Integer) * 8 - 1) - 1
#MIN_INTEGER = ~#MAX_INTEGER
;Erstellt zufällige Punkte
Procedure CreateRandomizedPointCloud(*pc.Triangulation::Triangulation, n.i, minX.d = 0.0, minY.d = 0.0, maxX.d = 1.0, maxY.d = 1.0)
Protected i.i, w.i
With *pc
Triangulation::clearTriangulation(*pc)
\n = n
ReDim \points(\n - 1)
w = Round(Sqr(\n), #PB_Round_Up)
For i = 0 To \n - 1
\points(i)\id = i
\points(i)\x = (Random(#MAX_INTEGER) * (maxX - minX)) / #MAX_INTEGER + minX
\points(i)\y = (Random(#MAX_INTEGER) * (maxY - minY)) / #MAX_INTEGER + minY
\points(i)\x = Mod(i * (maxX - minX) * (maxY - minY) / \n, maxX - minX)
\points(i)\y = Mod((i * (maxX - minX) * (maxY - minY) / \n) / (maxX - minX), maxY - minY)
;\points(i)\x = (0.5 + Mod(i, w)) * (maxX - minX) / w
;\points(i)\y = (0.5 + Int(i / w)) * (maxY - minY) / w
Next
EndWith
EndProcedure
#VECTORDRAWING = #True
Procedure DrawPoints(*main.Window)
Protected i.i, width.i, height.i
Static font.i
With *main\pc
CompilerIf #VECTORDRAWING
If Not IsFont(font)
font = LoadFont(#PB_Any, "Courier New", 16, #PB_Font_Bold)
EndIf
If StartVectorDrawing(CanvasVectorOutput(*main\canvasId))
VectorFont(FontID(font), 16)
CompilerElse
If StartDrawing(CanvasOutput(*main\canvasId))
DrawingMode(#PB_2DDrawing_Default)
CompilerEndIf
width = GadgetWidth(*main\canvasId)
height = GadgetHeight(*main\canvasId)
CompilerIf #VECTORDRAWING
VectorSourceColor(RGBA(255, 255, 255, 255))
FillVectorOutput()
CompilerElse
Box(0, 0, GadgetWidth(*main\canvasId), GadgetHeight(*main\canvasId), $ffffff)
DrawingMode(#PB_2DDrawing_Outlined)
CompilerEndIf
CompilerIf #VECTORDRAWING
If LastElement(\convexHull())
MovePathCursor(\convexHull()\p\x, \convexHull()\p\y)
VectorSourceColor(RGBA($ff, 0, 0, $3f))
ForEach \convexHull()
AddPathLine(\convexHull()\p\x, \convexHull()\p\y)
Next
DashPath(5, 10, #PB_Path_RoundCorner)
EndIf
CompilerElse
Protected ConvexHullSize.i = ListSize(\convexHull())
If (ConvexHullSize > 2)
Protected *last.Triangulation::ConvexHullPoint2D = 0
LastElement(\convexHull())
*last = @\convexHull()
i = 0
ForEach \convexHull()
LineXY(*last\p\x, *last\p\y, \convexHull()\p\x, \convexHull()\p\y, $0000ff)
DrawText((*last\p\x + \convexHull()\p\x) / 2, (*last\p\y + \convexHull()\p\y) / 2, Str(i), 0)
i + 1
*last = @\convexHull()
Next
EndIf
CompilerEndIf
Protected cc.Triangulation::CircumCircle, color.i
ForEach \triangles()
If (Triangulation::isTriangleDelaunay(@\triangles()))
color = RGBA(0, $ff, 0, $5f)
Else
color = RGBA($ff, 0, 0, $5f)
EndIf
Triangulation::getCircumCircle(\triangles()\p[0], \triangles()\p[1], \triangles()\p[2], @cc)
CompilerIf #VECTORDRAWING
VectorSourceColor(color)
AddPathCircle(cc\center\x, cc\center\y, Sqr(cc\radiusSq))
StrokePath(1)
CompilerElse
; Ohne das hier, kann Kreise malen ganz schön langsam sein
If cc\center\x >= 0 And cc\center\x < width And cc\center\y >= 0 And cc\center\y < height
Circle(cc\center\x, cc\center\y, Sqr(cc\radiusSq), color)
EndIf
CompilerEndIf
Next
ForEach \edges()
If \edges()\mark
color = RGBA(0, 0, $ff, $ff)
Else
color = RGBA($7f, $7f, $7f, $ff)
EndIf
CompilerIf #VECTORDRAWING
MovePathCursor(\edges()\p[0]\x, \edges()\p[0]\y)
VectorSourceColor(color)
AddPathLine(\edges()\p[1]\x, \edges()\p[1]\y)
StrokePath(1 + \edges()\mark * 2)
CompilerElse
LineXY(\edges()\p[0]\x, \edges()\p[0]\y, \edges()\p[1]\x, \edges()\p[1]\y, color)
CompilerEndIf
Next
CompilerIf Not #VECTORDRAWING
DrawingMode(#PB_2DDrawing_Transparent)
CompilerEndIf
For i = 0 To \n - 1
Protected s.s = Str(\points(i)\id)
CompilerIf #VECTORDRAWING
If \points(i) = \nextTriangle\p[1]
VectorSourceColor(RGBA($7f, $7f, $ff, $ff))
AddPathCircle(\points(i)\x, \points(i)\y, 10)
FillPath()
color = RGBA($ff, $3f, $3f, $ff)
Else
color = RGBA(0, 0, 0, $ff)
EndIf
VectorSourceColor(color)
AddPathCircle(\points(i)\x, \points(i)\y, 8)
FillPath()
MovePathCursor(\points(i)\x - VectorTextWidth(s) / 2 + 1, \points(i)\y - VectorTextHeight(s) / 2 + 1)
VectorSourceColor(RGBA($ff, $ff, $ff, $ff))
DrawVectorText(s)
CompilerElse
Circle(\points(i)\x, \points(i)\y, 1, $000000)
DrawText(\points(i)\x, \points(i)\y, s, $ff0000)
CompilerEndIf
Next
CompilerIf #VECTORDRAWING
For i = 0 To 1
If \nextTriangle[i]\p[0] And \nextTriangle[i]\p[1] And \nextTriangle[i]\p[2]
If Triangulation::isRightHand(\nextTriangle[i]\p[0], \nextTriangle[i]\p[1], \nextTriangle[i]\p[2])
color = RGBA($7f, $ff, $7f, $7f)
Else
color = RGBA($ff, $7f, $7f, $7f)
EndIf
MovePathCursor(\nextTriangle[i]\p[0]\x, \nextTriangle[i]\p[0]\y)
AddPathLine(\nextTriangle[i]\p[1]\x, \nextTriangle[i]\p[1]\y)
AddPathLine(\nextTriangle[i]\p[2]\x, \nextTriangle[i]\p[2]\y)
ClosePath()
VectorSourceColor(color)
FillPath()
EndIf
Next
CompilerEndIf
CompilerIf Not #VECTORDRAWING
DrawingMode(#PB_2DDrawing_Default)
DrawText(0, height - 20, " Points: " + Str(\n) + " " +
"Edges: " + Str(ListSize(\edges())) + " " +
"Triangles: " + Str(ListSize(\triangles())) + " " +
"Time: " + Str(\time) + " ms ", $0000ff, $ffffff)
CompilerEndIf
CompilerIf #VECTORDRAWING
StopVectorDrawing()
CompilerElse
StopDrawing()
CompilerEndIf
EndIf
EndWith
ProcedureReturn #True
EndProcedure
Procedure callback(*pc, *main.Window)
PostEvent(#PB_Event_FirstCustomValue, *main\id, 0)
WaitSemaphore(*main\mutex)
EndProcedure
Procedure triangulate(*main.Window)
Triangulation::triangulate(*main\pc, 0, @callback(), *main)
EndProcedure
Procedure CanvasEvent()
Protected x.i, y.i, gadgetId.i, *main.Window, seed.i
gadgetId = EventGadget()
*main = GetGadgetData(gadgetId)
With *main
x = GetGadgetAttribute(\canvasId, #PB_Canvas_MouseX)
y = GetGadgetAttribute(\canvasId, #PB_Canvas_MouseY)
Select (EventType())
Case #PB_EventType_RightClick
If Not \newThread
If IsThread(\thread)
SignalSemaphore(\mutex)
EndIf
EndIf
\newThread = #False
Case #PB_EventType_LeftButtonDown
\leftDown = #True
\clicked = Triangulation::getPoint(\pc, x, y)
Case #PB_EventType_Input
If GetGadgetAttribute(\canvasId, #PB_Canvas_Input) = 'f'
If IsThread(\thread)
SignalSemaphore(\mutex)
EndIf
EndIf
Case #PB_EventType_RightButtonDown
\rightDown = #True
If Not IsThread(\thread)
\thread = CreateThread(@triangulate(), *main)
\newThread = #True
EndIf
;Triangulation::triangulate(\pc, Triangulation::getPoint(\pc, x, y), @callback(), *main)
;DrawPoints(*main)
Case #PB_EventType_MouseMove
If (\leftDown)
Triangulation::setPoint(\pc, \clicked, x, y)
DrawPoints(*main)
EndIf
If Not IsThread(\thread)
If (\rightDown)
If (\clicked >= 0)
seed = \clicked
Else
seed = Triangulation::getPoint(\pc, x, y)
EndIf
Triangulation::triangulate(\pc, seed)
EndIf
;DrawPoints(*main)
EndIf
Case #PB_EventType_LeftButtonUp
\leftDown = #False
\clicked = -1
Case #PB_EventType_RightButtonUp
\rightDown = #False
\clicked = -1
EndSelect
EndWith
EndProcedure
Procedure Redraw()
Protected *main.Window = GetWindowData(EventWindow())
If Event() = #PB_Event_FirstCustomValue
DrawPoints(*main)
EndIf
EndProcedure
Procedure CreateMainWindow(*main.Window)
With *main
\id = OpenWindow(#PB_Any, 0, 0, \width, \height, \title, #PB_Window_MinimizeGadget | #PB_Window_SystemMenu | #PB_Window_ScreenCentered)
If (Not \id) : ProcedureReturn #False : EndIf
SetWindowData(\id, *main)
\canvasId = CanvasGadget(#PB_Any, 0, 0, \width, \height, #PB_Canvas_Keyboard)
\clicked = -1
SetGadgetData(\canvasId, *main)
If (Not \canvasId)
CloseWindow(\id)
ProcedureReturn #False
EndIf
BindGadgetEvent(\canvasId, @CanvasEvent())
BindEvent(#PB_Event_FirstCustomValue, @Redraw(), \id)
ProcedureReturn \id
EndWith
EndProcedure
Define.Triangulation::Triangulation pc
CreateRandomizedPointCloud(@pc, 53, 0, 0, #WIDTH - 1, #HEIGHT - 1)
Define.Window main
main\width = #WIDTH
main\height = #HEIGHT
main\title = "Triangulation Test"
main\pc = @pc
main\mutex = CreateSemaphore(0)
CreateMainWindow(@main)
DrawPoints(@main)
Repeat : Until WaitWindowEvent() = #PB_Event_CloseWindow