Schnelle Triangulation und konvexe Hülle in 2D

Hier könnt Ihr gute, von Euch geschriebene Codes posten. Sie müssen auf jeden Fall funktionieren und sollten möglichst effizient, elegant und beispielhaft oder einfach nur cool sein.
Lambda
Beiträge: 526
Registriert: 16.06.2011 14:38

Re: Schnelle Triangulation und konvexe Hülle in 2D

Beitrag von Lambda »

Ja danke. :wink: Sollte ich so oder so optional von außerhalb anwenden, teils ist der Effekt auch gewünscht.
Benutzeravatar
NicTheQuick
Ein Admin
Beiträge: 8675
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

Beitrag von NicTheQuick »

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.
Bild
Benutzeravatar
NicTheQuick
Ein Admin
Beiträge: 8675
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

Beitrag von NicTheQuick »

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.

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
Bild
Antworten