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

Schnelle Triangulation und konvexe Hülle in 2D

Beitrag von NicTheQuick »

Hallo,

durch Alexi angespornt habe ich mal den Triangulations-Algorithmus von s-hull.org auf eigene Art in Purebasic implementiert.
Heraus gekommen ist folgender Code, der außerdem eine simple GUI beinhaltet.

Zunächst einmal kann man mit #WIDTH und #HEIGHT die Größe der GUI einstellen. Mit #PRECISION kann man wählen zwischen Floats oder Doubles, die dann zum Berechnen verwendet werden sollen. Wenn man #MAKE_DELAUNAY auf #False setzt, dann wird das triangulierte Gitter nicht nach Delaunay optimiert, das heißt es wird kein Edge Flipping ausgeführt.
Am Ende nach der Procedure 'CreateMainWindow()' findet man die Initialisierung der PointCloud2D-Structure.

Code: Alles auswählen

Define.Triangulation::Triangulation pc

CreateRandomizedPointCloud(@pc, 10, 0, 0, #WIDTH - 1, #HEIGHT - 1)
Anhand von 'CreateRandomizedPointCloud()' lässt sich leicht feststellen wie eine Initialisierung abläuft. Im Grunde muss vor einer Neudimensionierung des Arrays 'PointCloud2D\points()' immer zuerst 'Triangulation::clearTriangulation(@pc)' aufgerufen werden, damit die Referenzen zu den Punkten nicht plötzlich auf ungültige Speicheradressen zeigen. Nach einem Aufruf von 'Triangulation::triangulate(@pc, seed)' mit einem beliebigen Punkt als Start, wird die Punktwolke nach dem Verfahren, was auf s-hull.org beschrieben ist, trianguliert.
Nach der Triangulation befinden sich alle Dreiecke in der LinkedList 'pc\triangles()', alle Kanten in 'pc\edges()' und die geordnete Liste von Punkten auf der konvexen Hülle in 'pc\convexHull()'. Das Element 'used.i' in 'Triangulation::ConvexHullPoint2D' wird nur während der Triangulation benötigt und hat danach keinen besonderen Nutzen mehr.
Die Punkte eines Dreiecks sind immer in der selben Richtung angeordnet, in diesem Fall im Uhrzeigersinn (Die Y-Achse muss dafür nach oben zeigen und die X-Achse nach rechts).

Nun zur GUI. Bei 'CreateRandomizedPointCloud()' kann mit dem zweiten Parameter angegeben werden wie viele zufällige Punkte erzeugt werden sollen. Danach startet das Hauptfenster. Mit der linken Maustaste kann man einen Punkt verschieben. Mit einem Rechtsklick ruft man 'Triangulation::triangulate()' auf, wobei der 'seed' hier der Punkt ist, der am nächsten am Mauszeiger war. Hält man die rechte Maustaste gedrückt und bewegt die Maus, wird nach jeder Bewegung erneut trianguliert. Man kann auch beide Tasten gleichzeitig benutzen und einen Punkt bewegen während gleichzeitig trianguliert wird.
In der GUI sieht man in rot die konvexe Hülle und in schwarz das triangulierte Gitter. Im Hintergrund sieht man die Umkreise eines jeden Dreiecks. Wenn #MAKE_DELAUNAY = #True ist, sollten diese Kreise alle grün sein. Ist ein Kreis rot, so bedeutet das, dass dieses Dreieck nicht Delaunay ist. Hat man eine Triangulation erzeugt, die insgesamt Delaunay ist, so kann man durch verschieben eines Punktes mit der linken Maustaste Dreiecke auch Nicht-Delaunay machen. Klickt man anschließend wieder mit der rechten Maustaste, wird wieder alles korrigiert.

Achja, wegen einem Bug in PureBasic, ist an einer Stelle noch ein Workaround nötig.

Versionen:
  • 12.01.2014 Mittags
    Erste Version.
  • 12.01.2014 21:00 Uhr
    Alles als Modul umgebaut.
Aber jetzt erst mal viel Spaß damit!

Code: Alles auswählen

EnableExplicit

#WIDTH = 1600
#HEIGHT = 1000

; If the bug mentioned at http://www.purebasic.fr/english/viewtopic.php?f=23&t=57961
; was solved, set this do #False
#BUG = #True

;-- START OF TRIANGULATION STRUCTURES AND FUNCTIONS

DeclareModule Triangulation
	#PRECISION = #PB_Float ;oder #PB_Double
	#MAKE_DELAUNAY = #True
	
	#DEBUG = #False
	
	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
	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
	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
	
	; Ein Dreieck merkt sich seine drei Eckpunkte und seine drei Kanten
	Structure Triangle2D
		*p.Point2DPC[3]
		*e.Edge2D[3]
	EndStructure
	

	
	Structure Point2DDiff
		*p.Point2D
		diff.prec
	EndStructure
	
	Structure ConvexHullPoint2D
		*p.Point2DPC
		used.i
	EndStructure
	
	Structure CircumCircle
		center.Point2D
		radius.prec
	EndStructure
	
	Structure Triangulation
		n.i
		time.i
		Array points.Point2DPC(1)
		List edges.Edge2D()
		List triangles.Triangle2D()
		List convexHull.ConvexHullPoint2D()
	EndStructure
	
	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)
EndDeclareModule

Module Triangulation

	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
	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\radius = Sqr(Pow(*cc\center\x - *a\x, 2) + Pow(*cc\center\y - *a\y, 2))
	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 < 0.0)
	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)
				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) ;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 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
			EndWith
			ProcedureReturn @*pc\triangles()
		EndIf
		
		ProcedureReturn #False
	EndProcedure
	
	Procedure.i isEdgeDelaunay(*edge.Edge2D) ;correct
		Protected cc.CircumCircle, i.i, j.i, *p.Point2D
		
		With *edge
			If (\t[1])
				For i = 0 To 1	;Iterate over adjacent triangles
					getCircumCircle(\t[i]\p[0], \t[i]\p[1], \t[i]\p[2], @cc)
					
					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 - cc\center\x, 2) + Pow(*p\y - cc\center\y, 2) - cc\radius * cc\radius < #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) ;correct
		Protected NewList *ndEdges.Edge2D()
		Protected cc.CircumCircle, *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
		
		ForEach *pc\edges()
			If (*pc\edges()\t[1] And (Not isEdgeDelaunay(*pc\edges())))
				If (AddElement(*ndEdges()))
					*ndEdges() = @*pc\edges()
				EndIf
				*pc\edges()\mark = #True
			Else
				*pc\edges()\mark = #False
			EndIf
		Next
		
		While FirstElement(*ndEdges())
			*actualEdge = *ndEdges()
			DeleteElement(*ndEdges())
			With *actualEdge
				\mark = #False
				doFlip = #False
				
				For i = 0 To 1	;Iterate over adjacent triangles
					getCircumCircle(\t[i]\p[0], \t[i]\p[1], \t[i]\p[2], @cc)
					
					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 - cc\center\x, 2) + Pow(*p\y - cc\center\y, 2) - cc\radius * cc\radius < #EPSILON)
								p_i(1 - i) = j
								doFlip = #True
							EndIf
						EndIf
					Next
				Next
				
				If (doFlip)
					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 " + i + " is not in the right order! [" + \t[i] + "]"
							error = #True
						EndIf
					Next
					If (error)
						Debug "Error on " + *actualEdge
						ProcedureReturn #False
					EndIf
					Debug "flip on " + *actualEdge + " with [" + \t[0] + "] and [" + \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
					
					;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])
								If (AddElement(*ndEdges()))
									*ndEdges() = *newEdge
									*newEdge\mark = #True
								EndIf
							EndIf
						Next
					Next
				EndIf
			EndWith
		Wend
		
	EndProcedure
	
	Procedure triangulate(*pc.Triangulation, seed.i = 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)
			
			;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\radius = Infinity()
			For i = 2 To \n - 1
				getCircumCircle(*x0, *xj, sortedPoints(i)\p, @C)
				If (C\radius < bestC\radius)
					bestIndex = i
					bestC = C
				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
			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)
			
			;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 prppagating 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
			
			Protected *p.Point2D, *last.ConvexHullPoint2D, alreadyAdded.i, convexHullSize.i, j.i
			For i = 3 To \n - 1
				*p = sortedPoints(i)\p
				
				alreadyAdded = #False
				convexHullSize = ListSize(\convexHull())
				FirstElement(\convexHull())
				*last = @\convexHull()
				NextElement(\convexHull())
				For j = 0 To convexHullSize - 1
					If (isRightHand(*last\p, *p, \convexHull()\p))
						addTriangle(*pc, *last\p, *p, \convexHull()\p)
						\convexHull()\used + 1
						*last\used + 1
						
						If (Not alreadyAdded)
							InsertElement(\convexHull())
							\convexHull()\p = *p
							NextElement(\convexHull())
							convexHullSize + 1
							alreadyAdded = #True
						EndIf
					EndIf
					*last = @\convexHull()
					If (Not NextElement(\convexHull()))
						FirstElement(\convexHull())
					EndIf
				Next
				
				ForEach \convexHull()
					If (\convexHull()\used = 2)
						DeleteElement(\convexHull(), 1)
					Else
						\convexHull()\used = 0
					EndIf
				Next
				
			Next
			
			;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)
			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
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
	
	With *pc
		Triangulation::clearTriangulation(*pc)
		\n = n
		ReDim \points(\n - 1)
		
		For i = 0 To \n - 1
			\points(i)\x = (Random(#MAX_INTEGER) * (maxX - minX)) / #MAX_INTEGER + minX
			\points(i)\y = (Random(#MAX_INTEGER) * (maxY - minY)) / #MAX_INTEGER + minY
		Next
	EndWith
EndProcedure

Procedure DrawPoints(*main.Window)
	Protected i.i
	
	With *main\pc
		If StartDrawing(CanvasOutput(*main\canvasId))
			DrawingMode(#PB_2DDrawing_Default)
			
			Box(0, 0, GadgetWidth(*main\canvasId), GadgetHeight(*main\canvasId), $ffffff)
			
			DrawingMode(#PB_2DDrawing_Outlined)
			Protected cc.Triangulation::CircumCircle, color.i
			color = $cfcfff
			ForEach \triangles()
				Triangulation::getCircumCircle(\triangles()\p[0], \triangles()\p[1], \triangles()\p[2], @cc)
				If (Triangulation::isTriangleDelaunay(@\triangles()))
					color = $cfffcf
				Else
					color = $cfcfff
				EndIf
				Circle(cc\center\x, cc\center\y, cc\radius, color)
			Next
			
			ForEach \edges()
				CompilerIf #BUG
					Protected.Triangulation::Point2D *p0, *p1
					*p0 = \edges()\p[0]
					*p1 = \edges()\p[1]
					LineXY(*p0\x, *p0\y, *p1\x, *p1\y, $7f7f7f)
				CompilerElse
					LineXY(\edges()\p[0]\x, \edges()\p[0]\y, \edges()\p[1]\x, \edges()\p[1]\y, $7f7f7f)
				CompilerEndIf
			Next
			
			DrawingMode(#PB_2DDrawing_Transparent)
			For i = 0 To \n - 1
				Circle(\points(i)\x, \points(i)\y, 1, $000000)
				;Plot(\points(i)\x, \points(i)\y, $000000)
			Next
			
			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
			
			DrawingMode(#PB_2DDrawing_Default)
			DrawText(0, 0, " Points: " + Str(\n) + "  " +
			               "Edges: " + Str(ListSize(\edges())) + "  " +
			               "Triangles: " + Str(ListSize(\triangles())) + "  " +
			               "Time: " + Str(\time) + " ms ", $0000ff, $ffffff)
			StopDrawing()
		EndIf
	EndWith
	
	ProcedureReturn #True
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_LeftButtonDown
				\leftDown = #True
				\clicked = Triangulation::getPoint(\pc, x, y)
			
			Case #PB_EventType_RightButtonDown
				\rightDown = #True
				Triangulation::triangulate(\pc, Triangulation::getPoint(\pc, x, y))
				DrawPoints(*main)
			
			Case #PB_EventType_MouseMove
				If (\leftDown)
					Triangulation::setPoint(\pc, \clicked, x, y)
				EndIf
				If (\rightDown)
					If (\clicked >= 0)
						seed = \clicked
					Else
						seed = Triangulation::getPoint(\pc, x, y)
					EndIf
					Triangulation::triangulate(\pc, seed)
				EndIf
				DrawPoints(*main)
				
			Case #PB_EventType_LeftButtonUp
				\leftDown = #False
				\clicked = -1
			
			Case #PB_EventType_RightButtonUp
				\rightDown = #False
				\clicked = -1
		EndSelect
	EndWith
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
		
		\canvasId = CanvasGadget(#PB_Any, 0, 0, \width, \height)
		\clicked = -1
		SetGadgetData(\canvasId, *main)
		If (Not \canvasId)
			CloseWindow(\id)
			ProcedureReturn #False
		EndIf
		
		BindGadgetEvent(\canvasId, @CanvasEvent())
	
		ProcedureReturn \id
	EndWith
EndProcedure

Define.Triangulation::Triangulation pc

CreateRandomizedPointCloud(@pc, 10, 0, 0, #WIDTH - 1, #HEIGHT - 1)

Define.Window main
main\width = #WIDTH
main\height = #HEIGHT
main\title = "Triangulation Test"
main\pc = @pc

CreateMainWindow(@main)

DrawPoints(@main)

Repeat : Until WaitWindowEvent() = #PB_Event_CloseWindow
Bild
Benutzeravatar
mk-soft
Beiträge: 3695
Registriert: 24.11.2004 13:12
Wohnort: Germany

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

Beitrag von mk-soft »

5 von 5 Sterne :allright:
Alles ist möglich, fragt sich nur wie...
Projekte ThreadToGUI / EventDesigner V3 / OOP-BaseClass-Modul
Downloads auf MyWebspace / OneDrive
Lambda
Beiträge: 526
Registriert: 16.06.2011 14:38

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

Beitrag von Lambda »

Spitze Nic, und nochmals Danke für deine Mühe. :allright:

Einen drauf setzen könntest du nur noch mit der dritten Dimension, volumetrisches Gelände oder Metaball ähnliche Wolken generieren. ^^
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 »

Prinzipiell weiß ich wie man aus einer 3D-Punktwolke ein dreidimensionales Gitter macht. Das wird ja zum Beispiel auch bei der Kinect gemacht, die als Rohdaten ebenfalls nur Punktwolken liefert. Allerdings ist die Umsetzung dann doch etwas schwieriger. :wink:
Bild
Lambda
Beiträge: 526
Registriert: 16.06.2011 14:38

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

Beitrag von Lambda »

An sich, wie bei Delaunay wäre es die zusätzliche Achse beim "Kreisen". Allerdings kommt noch einiges hinzu wenn bspw. innere Punkte ausgeschlossen werden würden, da hast du recht. Wäre das gegeben, wäre es tatsächlich eine Metaball Implementierung :? , da freut sich die CPU. :D

Musste eben ein Problem mit makeDelaunay feststellen. Ohne ist das Ergebnis wie gehabt, mit Delaunay ist bei:

Code: Alles auswählen

getCircumCircle(*a.SV2, *b.SV2, *c.SV2, *cc.CircumCircle)
*b Null. :| (im Projekt angewandt)
Edit: nur ein eigener Fehler bei der Einspeisung
Zuletzt geändert von Lambda am 12.01.2014 17:57, insgesamt 1-mal geändert.
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 »

Das ist merkwürdig. Aber die Zeile hilft mir so nicht weiter. Da bräuchte ich schon Daten, mit denen ich das testen kann. Wieso hast du die Strukturnamen eigentlich geändert?
Und hast du vorher immer PC_ClearTriangulation() benutzt?
Bild
Lambda
Beiträge: 526
Registriert: 16.06.2011 14:38

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

Beitrag von Lambda »

Nur um Namenskollisionen zu umgehen und weniger Schreibarbeit. Gehe dem erst noch genauer auf den Grund, vermute bei der Einspeisung der Punkte. Ein Beispiel wäre schlecht möglich, ist aber doch recht sicher kein Fehler der Triangulation.

Edit: Pardo, war nur ein Fehler bei der Einspeisung. :wink:
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 »

Wichtig ist, dass man wirklich die Structure 'Point2DPC' nimmt, die pro Punkt eine Liste verknüpfter Kanten zur Verfügung stellt. Da muss natürlich dann auch jedes Element mit 'InitializeStructure()' initialisiert werden, wenn man nicht wie ich gleich auf Standard-Arrays zurückgreift, die das automatisch initialisieren.
Bild
Lambda
Beiträge: 526
Registriert: 16.06.2011 14:38

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

Beitrag von Lambda »

An welcher Stelle wird denn ein Punkt zur Triangulation gewählt? Ich würde meine Konturen noch mit einwerfen damit diese nicht überschritten werden. Ist eine Zelle größer hat der Algo. zu großen Spielraum, was in engen Kurven etwas aus den Konturen wirkt.

Kein Fehler, es ist natürlich eine freie Punktwolke. Nur beim UV mapping wird z.B den Konturen nach gerechnet, was die Textur dann beim überschreiten kachelt. Ich würde daher vor Verbindung jede Linie auf Schnitt mit der Kontur prüfen.
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 »

Alexi hat geschrieben:An welcher Stelle wird denn ein Punkt zur Triangulation gewählt?
Wie meinst du das genau? Es wird nur ein Punkt gewählt, nämlich der, wo die Triangulation beginnen soll.

In 'PC_triangulate()' gibt es das Array 'sortedPoints()'. Zuerst werden in den Schritten 2 bis 4 die ersten drei Punkte für das initiale Dreieck berechnet. Dann werden in Schritt 5 die restlichen Punkte im selben Array, also ab Position 3, neu sortiert anhand der Entfernung zum Umkreismittelpunkt des initialen Dreiecks. In Schritt 7 werden schließlich die sortierten Punkte nacheinander hinzugefügt.

Hilft dir das weiter?
Bild
Antworten