Ich habe es allerdings anders gemacht als die. Das merkt man deutlich, wenn man die Codes vergleicht.
Das einzige, was noch fehlt, ist das Edge-Flipping von Delauney. Das mache ich dann vielleicht morgen oder so.
Kurze Erklärung:
Es wird eine Punktwolke aus 100 Punkten erstellt (Zeile 434). Dann kann man mit der linken Maustaste die Punkte noch zurecht rücken, wenn mal will. Mit einem Rechtsklick auf einen Punkt beginnt dann die Triangulation von dem gewählten Punkt aus. Man kann den Punkt auch dabei bewegen und das ganze in Echtzeit betrachten. In rot sieht man dann noch die konvexe Hülle, die zufällig mit berechnet wird. Das ganze geht bisher recht fix. Bei 10000 Punkten dauert es bei mir so 200 ms.
Viel Spaß.
Code: Alles auswählen
EnableExplicit
#WIDTH = 800
#HEIGHT = 800
#MAX_INTEGER = 1 << (SizeOf(Integer) * 8 - 1) - 1
#MIN_INTEGER = ~#MAX_INTEGER
#PRECISION = #PB_Float
CompilerIf #PRECISION = #PB_Float
Macro prec
f
EndMacro
CompilerElse
Macro prec
d
EndMacro
CompilerEndIf
Structure Window
id.i
width.i
height.i
title.s
canvasId.i
*pc.PointCloud2D
clicked.i
calcOnMove.i
EndStructure
Structure Line2D
*p.Point2D[2]
EndStructure
Structure Point2D
x.prec
y.prec
s.s
List *lines.Line2D()
EndStructure
Structure Vector2D Extends Point2D
EndStructure
Structure Point2DDiff
*p.Point2D
diff.prec
EndStructure
Structure ConvexHullPoint2D
*p.Point2D
used.i
EndStructure
Structure Triangle2D
*l.Line2D[3]
EndStructure
Structure PointCloud2D
min.Point2D
max.Point2D
n.i
time.i
Array points.Point2D(1)
List lines.Line2D()
List triangles.Triangle2D()
List convexHull.ConvexHullPoint2D()
EndStructure
Procedure PC_clearLines(*pc.PointCloud2D)
Protected i.i
ClearList(*pc\triangles())
For i = 0 To *pc\n - 1
ClearList(*pc\points(i)\lines())
Next
ClearList(*pc\lines())
ClearList(*pc\convexHull())
EndProcedure
;Erstellt zufällige Punkte
Procedure PC_RandomizePointCloud(*pc.PointCloud2D)
Protected i.i, *p.Point2D
With *pc
ReDim \points(\n - 1)
ClearList(\lines())
For i = 0 To \n - 1
*p = @\points(i)
*p\x = (Random(#MAX_INTEGER) * (\max\x - \min\x)) / #MAX_INTEGER + \min\x
*p\y = (Random(#MAX_INTEGER) * (\max\y - \min\y)) / #MAX_INTEGER + \min\y
Next
EndWith
EndProcedure
Procedure PC_DrawPoints(*main.Window, *pc.PointCloud2D = 0)
Protected i.i
If (*pc = 0)
*pc = *main\pc
If (*pc = 0)
ProcedureReturn #False
EndIf
EndIf
With *pc
If StartDrawing(CanvasOutput(*main\canvasId))
DrawingMode(#PB_2DDrawing_Default)
Box(0, 0, GadgetWidth(*main\canvasId), GadgetHeight(*main\canvasId), $ffffff)
DrawingMode(#PB_2DDrawing_Transparent)
For i = 0 To \n - 1
Circle(\points(i)\x, \points(i)\y, 2, $000000)
;DrawText(\points(i)\x, \points(i)\y, \points(i)\s, $ff0000)
Next
ForEach \lines()
LineXY(\lines()\p[0]\x, \lines()\p[0]\y, \lines()\p[1]\x, \lines()\p[1]\y, $dfdfdf)
Next
Protected ConvexHullSize.i = ListSize(\convexHull())
If (ConvexHullSize > 2)
Protected *last.ConvexHullPoint2D = 0
LastElement(\convexHull())
*last = @\convexHull()
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)
*last = @\convexHull()
Next
EndIf
DrawText(0, 0, "Time: " + Str(\time) + " ms", 0)
StopDrawing()
EndIf
EndWith
ProcedureReturn #True
EndProcedure
Procedure PC_getPoint(*pc.PointCloud2D, x.prec, y.prec)
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 PC_setPoint(*pc.PointCloud2D, i.i, x.prec, y.prec)
If (i >= 0 And i < *pc\n)
*pc\points(i)\x = x
*pc\points(i)\y = y
EndIf
EndProcedure
Structure CircumCircle
center.Point2D
radius.prec
EndStructure
Procedure PC_getCircumCircle(*a.Point2D, *b.Point2D, *c.Point2D, *cc.CircumCircle)
Protected a.prec = Sqr(Pow(*c\x - *b\x, 2) + Pow(*c\y - *b\y, 2))
Protected b.prec = Sqr(Pow(*c\x - *a\x, 2) + Pow(*c\y - *a\y, 2))
Protected c.prec = Sqr(Pow(*b\x - *a\x, 2) + Pow(*b\y - *a\y, 2))
;http://en.wikipedia.org/wiki/Circumcenter
Protected s.prec = 0.5 * (a + b + c)
*cc\radius = a * b * c / Sqr(s * (s - a) * (s - b) * (s - c))
Protected d.prec = 2 * (*a\x * (*b\y - *c\y) +
*b\x * (*c\y - *a\y) +
*c\x * (*a\y - *b\y))
*cc\center\x = (a * a * (*b\y - *c\y) +
b * b * (*c\y - *a\y) +
c * c * (*a\y - *b\y)) / d
*cc\center\y = (a * a * (*c\x - *b\x) +
b * b * (*a\x - *c\x) +
c * c * (*b\x - *a\x)) / d
*cc\center\x = (*a\x + *b\x + *c\x) / 3.0
*cc\center\y = (*a\y + *b\y + *c\y) / 3.0
EndProcedure
Procedure.i PC_isRightHand(*a.Point2D, *b.Point2D, *c.Point2D)
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 PC_orderRightHand(*a.Point2D, *p_b.Integer, *p_c.Integer)
If (Not PC_isRightHand(*a, *p_b\i, *p_c\i))
Swap *p_b\i, *p_c\i
EndIf
EndProcedure
Procedure PC_addLine(*pc.PointCloud2D, *a.Point2D, *b.Point2D)
Protected lineExistsAlready.i = #False
If (*a = *b)
ProcedureReturn #False
EndIf
ForEach *a\lines()
If (*a\lines()\p[0] = *b Or *a\lines()\p[1] = *b)
ProcedureReturn *a\lines()
EndIf
Next
If AddElement(*pc\lines())
*pc\lines()\p[0] = *a
*pc\lines()\p[1] = *b
AddElement(*a\lines())
*a\lines() = @*pc\lines()
AddElement(*b\lines())
*b\lines() = @*pc\lines()
ProcedureReturn @*pc\lines()
EndIf
EndProcedure
Procedure PC_addTriangle(*pc.PointCloud2D, *a.Point2D, *b.Point2D, *c.Point2D)
;Assume AddElements() never fails!
;PC_orderRightHand(*a, @*b, @*c)
AddElement(*pc\triangles())
*pc\triangles()\l[0] = PC_addLine(*pc, *a, *b)
*pc\triangles()\l[1] = PC_addLine(*pc, *b, *c)
*pc\triangles()\l[2] = PC_addLine(*pc, *c, *a)
ProcedureReturn @*pc\triangles()
EndProcedure
Procedure PC_triangulate(*pc.PointCloud2D, seed.i)
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.Point2D = @\points(seed)
;2. Sort according to |x_i - x_0|^2.
For i = 0 To \n - 1
sortedPoints(i)\p = @\points(i)
\points(i)\s = ""
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.Point2D = 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 minRadius.prec = Infinity()
Protected bestIndex.i
Protected cc.CircumCircle
Protected C.Point2D
For i = 2 To \n - 1
PC_getCircumCircle(*x0, *xj, sortedPoints(i)\p, @cc)
If (cc\radius < minRadius)
;Debug Str(cc\center\x) + ", " + Str(cc\center\y) + " (" + StrF(cc\radius) + ")"
minRadius = cc\radius
bestIndex = i
C = cc\center
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 = 0 To \n - 1
sortedPoints(i)\diff = Pow(C\x - sortedPoints(i)\p\x, 2) + Pow(C\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.Point2D = sortedPoints(2)\p
PC_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.
PC_clearLines(*pc)
PC_addTriangle(*pc, *x0, *xk, *xj)
ClearList(\convexHull())
AddElement(\convexHull()) : \convexHull()\p = *x0 : *x0\s = "x_0"
AddElement(\convexHull()) : \convexHull()\p = *xk : *xk\s = "x_k"
AddElement(\convexHull()) : \convexHull()\p = *xj : *xj\s = "x_j"
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 (PC_isRightHand(*last\p, *p, \convexHull()\p))
PC_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).
;TODO
\time = ElapsedMilliseconds() - \time
EndWith
EndProcedure
Procedure CanvasEvent()
Protected x.i, y.i, gadgetId.i, *main.Window
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
\clicked = PC_getPoint(\pc, x, y)
Case #PB_EventType_RightButtonDown
\clicked = PC_getPoint(\pc, x, y)
\calcOnMove = #True
Case #PB_EventType_MouseMove
PC_setPoint(\pc, \clicked, x, y)
If (\calcOnMove)
PC_triangulate(\pc, \clicked)
EndIf
PC_DrawPoints(*main)
Case #PB_EventType_LeftButtonUp
\clicked = -1
Case #PB_EventType_RightButtonUp
\calcOnMove = #False
PC_triangulate(\pc, \clicked)
PC_DrawPoints(*main)
\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.PointCloud2D pc
pc\min\x = 0
pc\min\y = 0
pc\max\x = #WIDTH - 1
pc\max\y = #HEIGHT - 1
pc\n = 1000
PC_RandomizePointCloud(@pc)
Define.Window main
main\width = #WIDTH
main\height = #HEIGHT
main\title = "Triangulation Test"
main\pc = @pc
CreateMainWindow(@main)
PC_DrawPoints(@main)
Repeat : Until WaitWindowEvent() = #PB_Event_CloseWindow