Quadrate erkennen

[Christian+]

Ich habe einige Punkte mit x und y Koordinate und möchte nun herausfinden ob durch beliebige 4 davon ein Quadrat gebildet wird auch gedrehte um dieses dann einzeichnen zu können. Meine Versuche dazu sind bisher alle gescheitert weswegen ich nun hoffe, dass mir vielleicht jemand einen Ansatz geben kann wie ich das lösen könnte.

[gnasen]

Wenn du mindestens 2 Punkte gegeben hast, hast du ja schonmal die Seitenlänge. und den Winkel um den es gedreht wurde.
Damit kannst du die fehlenden 2 Punkte berechnen und abgleichen, ob diese 4 Punkte nun ein Quadrat gleichen.

[Nino]

Interessante Aufgabe.

Wenn du mindestens 2 Punkte gegeben hast, hast du ja schonmal die Seitenlänge. und den Winkel um den es gedreht wurde.

Ich hab's jetzt mal ganz ohne Winkel gemacht, nur mit Abständen zwischen den Punkten.

Mein Code sucht zunächst nach Rauten ( = Rhomben bzw. Salmis ), also nach Vierecken deren Seiten alle gleich lang sind. Das können dann Quadrate sein, müssen aber nicht. Das weiß jeder, der mal ein Viereck aus gleich langen Stäben gebaut hat, die gelenkig miteinander verbunden sind: Die Stäbe lassen sich gegeneinander bewegen, man nennt das daher auch "Gelenkviereck". (Bei einem Dreieck geht das übrigens nicht, das ist immer starr.) Für ein Quadrat braucht es zusätzlich rechte Winkel. Die sind allerdings zwangsläufig vorhanden, wenn beide Diagonalen gleich lang sind.

Das Hauptproblem des folgenden Codes ist -- falls er keine Bugs enthält --, dass er bei vielen Punkten lange dauert! Ob es schneller geht, weiß ich momentan leider nicht.

Code: Alles auswählen

; PB 4.51
; Brute Force-Methode zum Finden aller Quadrate in einer gegebenen Punktmenge.

EnableExplicit

Structure PointD
   x.d
   y.d
EndStructure

Procedure.i AreSame (a.d, b.d)
   ; in : 2 Fließkommazahlen mit doppelter Genauigkeit
   ; out: #True, wenn beide mit der Toleranz #Epsilon gleich sind
   #Epsilon = 0.000001

   If Abs(a-b) < #Epsilon
      ProcedureReturn #True
   Else
      ProcedureReturn #False
   EndIf
EndProcedure

Procedure.d Distance (*p1.PointD, *p2.PointD)
   ; in : 2 Punkte
   ; out: Abstand zwischen den beiden Punkten
   Protected dx.d, dy.d

   dx = *p2\x - *p1\x
   dy = *p2\y - *p1\y
   ProcedureReturn Sqr(dx*dx + dy*dy)
EndProcedure

Macro ShowSquare
   Debug "Quadrat:"
   Debug "A(" + StrD(p(a)\x,1) + "|" + StrD(p(a)\y,1) + ")"
   Debug "B(" + StrD(p(b)\x,1) + "|" + StrD(p(b)\y,1) + ")"
   Debug "C(" + StrD(p(c)\x,1) + "|" + StrD(p(c)\y,1) + ")"
   Debug "D(" + StrD(p(d)\x,1) + "|" + StrD(p(d)\y,1) + ")"
   Debug ""
EndMacro

Procedure FindSquares (Array p.PointD(1))
   Protected.d seitenlaenge
   Protected.i a, b, c, d    ; Indizes der Eckpunkte eines zu prüfenden Vierecks
   Protected.i n = ArraySize(p())
   
   For a = 1 To n - 3
      For b = a + 1 To n - 2
         seitenlaenge = Distance(@p(a), @p(b))
         For c = b + 1 To n - 1
            If AreSame(Distance(@p(c),@p(b)), seitenlaenge)  ; C ist möglicher Quadrat-Nachbar von B.
               For d = c + 1 To n
                  If AreSame(Distance(@p(c),@p(d)), seitenlaenge) And AreSame(Distance(@p(a),@p(d)), seitenlaenge)
                     ; Die Punkte A, B, C und D bilden schonmal eine Raute (Rhombus).
                     If AreSame(Distance(@p(a),@p(c)), Distance(@p(b),@p(d)))
                        ; Beide Diagonalen sind gleich lang => Quadrat.
                        ShowSquare
                     EndIf
                  EndIf
               Next
            ElseIf AreSame(Distance(@p(c),@p(a)), seitenlaenge)  ; C ist möglicher Quadrat-Nachbar von A.
               For d = c + 1 To n
                  If AreSame(Distance(@p(c),@p(d)), seitenlaenge) And AreSame(Distance(@p(b),@p(d)), seitenlaenge)
                     ; Die Punkte A, B, C und D bilden schonmal eine Raute (Rhombus).
                     If AreSame(Distance(@p(a),@p(d)), Distance(@p(b),@p(c)))
                        ; Beide Diagonalen sind gleich lang => Quadrat.
                        ShowSquare
                     EndIf
                  EndIf
               Next
            EndIf
         Next
      Next
   Next
EndProcedure


;-- Demo
Dim p.PointD(10)

; Array mit Punkten füllen
p( 1)\x = 3.5 : p( 1)\y = 6.0
p( 2)\x = 1.5 : p( 2)\y = 1.5
p( 3)\x = 5.5 : p( 3)\y = 1.5
p( 4)\x = 2.5 : p( 4)\y = 2.5
p( 5)\x = 1.5 : p( 5)\y = 3.5
p( 6)\x = 3.5 : p( 6)\y = 1.5
p( 7)\x = 5.5 : p( 7)\y = 6.0
p( 8)\x = 3.5 : p( 8)\y = 3.5
p( 9)\x = 3.5 : p( 9)\y = 5.0
p(10)\x = 5.5 : p(10)\y = 3.5

FindSquares(p())

Debug "Fertig."

Grüße, Nino

[Christian+]

Danke die Ideen von euch sind sehr gut ich habe gleich mal versucht drauf aufzubauen. Dazu hab ich mal versucht noch einen Weg zu finden die Idee von gnasen mit dem berechnen von 2 Punkten irgendwie umzusetzen da ich denke dass das eventuell schneller sein sollte auch wenn das nicht so wichtig ist da allzu große Mengen an Punkten ich der Funktion nicht als Input geben werde. Dabei kam der folgende Quellecode raus ich muss den nur nochmal mit verschieden Testdaten füttern um sicherzustellen das er auch immer richtig funktioniert.

Code: Alles auswählen

EnableExplicit

Macro ShowSquare
  LineXY(p(a)\x, p(a)\y, p(b)\x, p(b)\y, RGB(0,0,255))
  LineXY(p(b)\x, p(b)\y, p(c)\x, p(c)\y, RGB(0,0,255))
  LineXY(p(c)\x, p(c)\y, p(d)\x, p(d)\y, RGB(0,0,255))
  LineXY(p(d)\x, p(d)\y, p(a)\x, p(a)\y, RGB(0,0,255))
  Debug "Quadrat:"
  Debug "A(" + Str(p(a)\x) + "|" + Str(p(a)\y) + ")"
  Debug "B(" + Str(p(b)\x) + "|" + Str(p(b)\y) + ")"
  Debug "C(" + Str(p(c)\x) + "|" + Str(p(c)\y) + ")"
  Debug "D(" + Str(p(d)\x) + "|" + Str(p(d)\y) + ")"
  Debug ""
EndMacro

Procedure FindSquares(Array p.Point(1))
  
  Protected.i a, b, c, d ;Punkte des Quadrat
  Protected.i vektor.Point, z.Point
  Protected.i n = ArraySize(p()) ;Anzahl Punkte
  
  For a = 0 To n ;für alle Punkte
    For b = a+1 To n ;für alle noch nicht geprüften Punkte
      
      ;Vektor zwischen a und b ermitteln
      vektor\x = p(b)\x - p(a)\x
      vektor\y = p(b)\y - p(a)\y
      
      ;schauen ob Punkt c im Array
      For c = a+1 To n 
        
        ;Punkt c berechnen
        z\x = p(b)\x-vektor\y
        z\y = p(b)\y+vektor\x
        
        If p(c)\x = z\x And p(c)\y = z\y
          
          ;Punkt d berechnen
          z\x = p(a)\x-vektor\y
          z\y = p(a)\y+vektor\x
          
          ;schauen ob Punkt d im Array
          For d = a+1 To n 
            If p(d)\x = z\x And p(d)\y = z\y
              
              ShowSquare
              Break 2
              
            EndIf
          Next d
          
        EndIf
      Next c
      
    Next b
  Next a
  
EndProcedure

Dim p.Point(29)

p(0)\x = 100 : p(0)\y =  80
p(1)\x =  80 : p(1)\y = 180
p(2)\x =  80 : p(2)\y = 220
p(3)\x = 260 : p(3)\y = 160
p(4)\x = 240 : p(4)\y = 180
p(5)\x = 160 : p(5)\y = 240
p(6)\x = 100 : p(6)\y =  40
p(7)\x =  80 : p(7)\y = 200
p(8)\x = 180 : p(8)\y = 140
p(9)\x =  80 : p(9)\y =  60

p(10)\x = 120 : p(10)\y = 140
p(11)\x = 260 : p(11)\y =  80
p(12)\x =  40 : p(12)\y = 260
p(13)\x = 220 : p(13)\y =  60
p(14)\x = 200 : p(14)\y = 260
p(15)\x =  40 : p(15)\y = 160
p(16)\x =  60 : p(16)\y = 200
p(17)\x = 260 : p(17)\y = 200
p(18)\x = 160 : p(18)\y = 160
p(19)\x = 100 : p(19)\y = 200

p(20)\x =  80 : p(20)\y = 260
p(21)\x = 140 : p(21)\y = 160
p(22)\x =  40 : p(22)\y = 220
p(23)\x = 200 : p(23)\y = 220
p(24)\x = 260 : p(24)\y = 180
p(25)\x = 160 : p(25)\y = 220
p(26)\x = 160 : p(26)\y = 200
p(27)\x = 220 : p(27)\y =  80
p(28)\x =  60 : p(28)\y = 160
p(29)\x =  40 : p(29)\y = 100

Define.i image, i

OpenWindow(#PB_Any, 100, 100, 300, 300, "Quadrate erkennen")

image = CreateImage(#PB_Any, 300, 300)

StartDrawing(ImageOutput(image))

  Box(0, 0, 300, 300)
  
  FindSquares(p())
  
  For i = 0 To ArraySize(p())
    Circle(p(i)\x, p(i)\y, 1, RGB(0,0,0))
  Next
 
StopDrawing()

ImageGadget(0, 0, 0, 300, 300, ImageID(image))
  
Repeat
  
Until WaitWindowEvent() = #PB_Event_CloseWindow

End

[DarkDragon]

Dabei kam der folgende Quellecode raus ich muss den nur nochmal mit verschieden Testdaten füttern um sicherzustellen das er auch immer richtig funktioniert.

Du weißt aber schon, dass das kein wissenschaftlich verwertbarer Beweis ist? Dann funktioniert das immernoch nur für deine Testfälle und nicht immer.

[Sirius-2337]

Mir ist aufgefallen, dass 4 auf der selben stelle liegenden Punkte auch aus Quadrat erkannt werden, aber sonst siehts ganz gut aus.

[Christian+]

Vielen Dank für das Feedback.

@DarkDragon wenn das Programm fehlerfrei ist sollte ja alles passen beruht ja auf meinen mathematischen Kenntnissen die hoffentlich richtig angewandt sind und um sicherzustellen das ich nicht doch irgendeinen blöden Fehler gemacht habe fällt mir nichts besseres ein wie es einfach gründlich zu testen lass mich aber gerne einer besseren Methode belehren.

@Sirius-2337 das mit den 4 Punkten auf derselben Stelle ist mir auch schon aufgefallen da das für mich aber egal ist habe ich da keine Abfrage eingebaut werde ich aber glaube mal noch machen muss ich ja nur überprüfen ob die Länge des Vektors ungleich 0 ist.

[Nino]

Dabei kam der folgende Quellecode raus

Gefällt mir besser als mein eigener. Ich habe noch ein paar kleine Änderungsvorschläge für die Prozedur FindSquares().

Eine kosmetische Änderung:

Code: Alles auswählen

Protected.i vektor.Point, z.Point

Hier würde ich nicht den Standardtyp auf .i setzen, weil anschließend die folgenden Variablen ohnehin einzeln als .Point deklariert werden. Meines Erachtens besser lesbar:

Code: Alles auswählen

Protected vektor.Point, z.Point

oder

Code: Alles auswählen

Protected.Point vektor, z

Hier ist eine kleine Zeitersparnis möglich:

Code: Alles auswählen

      ;schauen ob Punkt c im Array
      For c = a+1 To n 
        
        ;Punkt c berechnen
        z\x = p(b)\x-vektor\y
        z\y = p(b)\y+vektor\x

Die Berechnung muss nicht bei jedem Schleifendurchlauf neu erfolgen, sie braucht nur 1x vor Beginn der For/Next-Schleife ausgeführt zu werden. (Weil z aber auch noch später für den Vergleich mit p(d) verwendet wird, braucht man dafür eine zusätzliche Variable.)

Insgesamt sieht das bei mir dann so aus:

Code: Alles auswählen

EnableExplicit

Macro ShowSquare
  LineXY(p(a)\x, p(a)\y, p(b)\x, p(b)\y, RGB(0,0,255))
  LineXY(p(b)\x, p(b)\y, p(c)\x, p(c)\y, RGB(0,0,255))
  LineXY(p(c)\x, p(c)\y, p(d)\x, p(d)\y, RGB(0,0,255))
  LineXY(p(d)\x, p(d)\y, p(a)\x, p(a)\y, RGB(0,0,255))
  Debug "Quadrat:"
  Debug "A(" + Str(p(a)\x) + "|" + Str(p(a)\y) + ")"
  Debug "B(" + Str(p(b)\x) + "|" + Str(p(b)\y) + ")"
  Debug "C(" + Str(p(c)\x) + "|" + Str(p(c)\y) + ")"
  Debug "D(" + Str(p(d)\x) + "|" + Str(p(d)\y) + ")"
  Debug ""
EndMacro

Procedure FindSquares (Array p.Point(1))
   Protected.i a, b, c, d          ; Punkte eines Quadrats
   Protected.Point vektor, c_, d_
   Protected.i n = ArraySize(p())
   
   For a = 0 To n-1
      For b = a+1 To n  ; für alle noch nicht geprüften Punkte
         ; Vektor von a nach b ermitteln
         vektor\x = p(b)\x - p(a)\x
         vektor\y = p(b)\y - p(a)\y
         
         ; Punkt c_ berechnen
         c_\x = p(b)\x - vektor\y
         c_\y = p(b)\y + vektor\x
            
         ; schauen ob Punkt c_ im Array ist
         For c = a+1 To n
            If p(c)\x = c_\x And p(c)\y = c_\y
               ; Punkt d_ berechnen
               d_\x = p(a)\x - vektor\y
               d_\y = p(a)\y + vektor\x
               
               ; schauen ob Punkt d_ im Array ist
               For d = a+1 To n
                  If p(d)\x = d_\x And p(d)\y = d_\y
                     ShowSquare
                     Break 2     ; nächsten Punkt B prüfen
                  EndIf
               Next d
            EndIf
         Next c
         
      Next b
   Next a
EndProcedure


Define.i image, i
Dim p.Point(29)

p(0)\x = 100 : p(0)\y =  80
p(1)\x =  80 : p(1)\y = 180
p(2)\x =  80 : p(2)\y = 220
p(3)\x = 260 : p(3)\y = 160
p(4)\x = 240 : p(4)\y = 180
p(5)\x = 160 : p(5)\y = 240
p(6)\x = 100 : p(6)\y =  40
p(7)\x =  80 : p(7)\y = 200
p(8)\x = 180 : p(8)\y = 140
p(9)\x =  80 : p(9)\y =  60

p(10)\x = 120 : p(10)\y = 140
p(11)\x = 260 : p(11)\y =  80
p(12)\x =  40 : p(12)\y = 260
p(13)\x = 220 : p(13)\y =  60
p(14)\x = 200 : p(14)\y = 260
p(15)\x =  40 : p(15)\y = 160
p(16)\x =  60 : p(16)\y = 200
p(17)\x = 260 : p(17)\y = 200
p(18)\x = 160 : p(18)\y = 160
p(19)\x = 100 : p(19)\y = 200

p(20)\x =  80 : p(20)\y = 260
p(21)\x = 140 : p(21)\y = 160
p(22)\x =  40 : p(22)\y = 220
p(23)\x = 200 : p(23)\y = 220
p(24)\x = 260 : p(24)\y = 180
p(25)\x = 160 : p(25)\y = 220
p(26)\x = 160 : p(26)\y = 200
p(27)\x = 220 : p(27)\y =  80
p(28)\x =  60 : p(28)\y = 160
p(29)\x =  40 : p(29)\y = 100

OpenWindow(#PB_Any, 100, 100, 300, 300, "Quadrate finden")

image = CreateImage(#PB_Any, 300, 300)

StartDrawing(ImageOutput(image))
Box(0, 0, 300, 300)

FindSquares(p())

For i = 0 To ArraySize(p())
   Circle(p(i)\x, p(i)\y, 1, RGB(0,0,0))
Next
StopDrawing()

ImageGadget(0, 0, 0, 300, 300, ImageID(image))

Repeat
Until WaitWindowEvent() = #PB_Event_CloseWindow

da ich denke dass das eventuell schneller sein sollte auch wenn das nicht so wichtig ist da allzu große Mengen an Punkten ich der Funktion nicht als Input geben werde

Ja OK, aber vielleicht hat später mal jemand anderes Bedarf für solchen Code, und derjenige muss evtl. doch größere Mengen an Punkten verarbeiten. Immerhin steigt momentan die Laufzeit ungefähr mit der 4. Potenz der Punktanzahl, das ist ziemlich ungünstig. Ich habe noch eine Idee zur Beschleunigung des Codes, aber die will ich erstmal prüfen.

Grüße, Nino

[NicTheQuick]

Mit einem Kd-Tree könnt ihr in logarithmischer Zeit Punkte finden. Damit sollte noch einmal ein erheblicher Geschwindigkeitszuwachs möglich sein. der Kd-Tree ist allerdings nur für statische Dinge gut. Also dann, wenn die Punkte von Anfang an fest stehen und sich während der Prüfung auf Quadrate nicht verändern.

[Nino]

Sorry wegen der späten Antwort ... Leider hatte ich zwischendurch wenig Zeit für die PB-Programmierung.

Mit einem Kd-Tree könnt ihr in logarithmischer Zeit Punkte finden. Damit sollte noch einmal ein erheblicher Geschwindigkeitszuwachs möglich sein.

Genau an sowas hatte ich gedacht, ohne dass ich allerdings den Namen dafür gekannt hätte ... vielen Dank für den Hinweis. Ich habe in Sedgewick (1991): Algorithmen, S. 434 ff. etwas über mehrdimensionale Suchbäume gefunden.

Zusätzlich habe ich noch eine Version mit einer Map geschrieben. Der folgende Code vergleicht die Geschwindigketen der Quadrat-Erkennung, wenn lauter unterschiedliche zufällige Punkte in einem Array, in einer Map oder in einem KD-Baum gespeichert sind. Vorher habe ich natürlich geprüft, dass alle 3 Algorithmen immer die selben Ergebnisse liefern.

//edit 22.09.2013:

• Kleinen Bug in KDtree_FindPoint() behoben.
• Der Code war für PB 4.51. Er ist jetzt so angepasst, dass er auch unter den neueren PB-Versionen läuft.
• Einige andere kleine Verbesserungen.

Code: Alles auswählen

; <http://www.purebasic.fr/german/viewtopic.php?f=4&t=24473>, 22.09.2013

CompilerIf #PB_Compiler_Debugger
   Debug "Für Zeitmessungen den Debugger ausschalten!"
   End
CompilerEndIf

EnableExplicit

CompilerIf #PB_Compiler_Version < 510
   Macro Bool (_expression_)
      ((_expression_) Or #False)
   EndMacro
CompilerEndIf


Procedure.i CreateUniqueRandomPoints (numPoints.i, maxX.i, maxY.i, Array a.Point(1))
   ; -- Array aus *unterschiedlichen* zufälligen Punkten erzeugen
   ; in : numPoints: Anzahl der Punkte
   ;      maxX     : maximaler x-Wert (Minimum ist 0)
   ;      maxY     : maximaler y-Wert (Minimum ist 0)
   ;      a()      : Array für die Punkte
   ; out: a()          : mit unterschiedlichen zufälligen Punkten gefülltes Array
   ;      Funktionswert: 0 bei Fehler, sonst 1
   Protected NewMap test()
   Protected point$
   Protected.i x, y, i=0
   
   If numPoints > (maxX+1)*(maxY+1)  ; numPoints ist zu groß für
      ProcedureReturn 0              ; den erlaubten Bereich.
   EndIf
   
   Dim a.Point(numPoints-1)
   
   While i < numPoints
      x = Random(maxX)
      y = Random(maxY)
      
      point$ = Str(x) + "|" + Str(y)
      If FindMapElement(test(), point$) = 0
         AddMapElement(test(), point$, #PB_Map_NoElementCheck)
         a(i)\x = x
         a(i)\y = y
         i + 1
      EndIf
   Wend
   
   ProcedureReturn 1
EndProcedure


#MaxX = 100
#MaxY = 100
#NumPoints = 750

Define.i i, s1, s2a, s2b, s3a, s3b, t, t1, t2a, t2b, t3a, t3b
Define.i NumSquares1, NumSquares2, NumSquares3
Dim RandomPoint.Point(0)

;-- lauter *unterschiedliche* zufällige Punkte erzeugen (Zeit wird nicht mitgezählt)
If CreateUniqueRandomPoints(#NumPoints, #MaxX, #MaxY, RandomPoint()) = 0
   MessageRequester("Fehler", "#NumPoints ist zu groß.")
   End
EndIf

;------------------------------------------------------------------------------

Procedure FindSquares_Array (Array p.Point(1))
   Protected.i a, b, c, d          ; Punkte eines Quadrats
   Protected.Point vektor, c_, d_
   Protected.i count=0, last=ArraySize(p())
   
   For a = 0 To last-1
      For b = a+1 To last  ; für alle noch nicht geprüften Punkte
         ; Vektor von a nach b ermitteln
         vektor\x = p(b)\x - p(a)\x
         vektor\y = p(b)\y - p(a)\y
         
         ; Punkt c_ berechnen
         c_\x = p(b)\x - vektor\y
         c_\y = p(b)\y + vektor\x
         
         ; schauen ob Punkt c_ im Array ist
         For c = a+1 To last
            If p(c)\x = c_\x And p(c)\y = c_\y
               ; Punkt d_ berechnen
               d_\x = p(a)\x - vektor\y
               d_\y = p(a)\y + vektor\x
               
               ; schauen ob Punkt d_ im Array ist
               For d = a+1 To last
                  If p(d)\x = d_\x And p(d)\y = d_\y
                     count + 1
                     Break 2     ; nächsten Punkt B prüfen
                  EndIf
               Next d
            EndIf
         Next c
         
      Next b
   Next a
   
   ProcedureReturn count
EndProcedure


s1 = ElapsedMilliseconds()
NumSquares1 = FindSquares_Array(RandomPoint())
t1 = ElapsedMilliseconds() - s1

;------------------------------------------------------------------------------

Structure Node
   p.Point
   f.i
EndStructure


Procedure FindSquares_Map (Map p.Node())
   Protected.Point a, b, c, d    ; Punkte eines Quadrats
   Protected.Point vektor
   Protected a_cur$, b_cur$
   Protected.i count=0
   
   ResetMap(p())
   While NextMapElement(p())     ; für alle Punkte
      p()\f = #True
      a_cur$ = MapKey(p())       ; Schlüssel des aktuellen Elements merken
      a = p()\p
      While NextMapElement(p())  ; für alle noch nicht geprüften Punkte
         b_cur$ = MapKey(p())    ; Schlüssel des aktuellen Elements merken
         b = p()\p
         
         ; Vektor von a nach b ermitteln
         vektor\x = b\x - a\x
         vektor\y = b\y - a\y
         
         ; Punkt c berechnen
         c\x = b\x - vektor\y
         c\y = b\y + vektor\x
         
         ; schauen ob Punkt c in der Map ist und noch nicht "dran" war
         If FindMapElement(p(), Str(c\x)+"|"+Str(c\y)) And p()\f = #False
            ; Punkt d berechnen
            d\x = a\x - vektor\y
            d\y = a\y + vektor\x
            
            ; schauen ob Punkt d in der Map ist und noch nicht "dran" war
            If FindMapElement(p(), Str(d\x)+"|"+Str(d\y)) And p()\f = #False
               count + 1
            EndIf
            
         EndIf
         FindMapElement(p(), b_cur$)  ; zurück zum aktuellen Map-Element der inneren Schleife
      Wend
      FindMapElement(p(), a_cur$)     ; zurück zum aktuellen Map-Element der äußeren Schleife
   Wend
   
   ProcedureReturn count
EndProcedure


NewMap p2.Node()
; NewMap p2.Node(2048)  ; bringt nicht viel, der KD-Baum ist trotzdem deutlich schneller

s2a = ElapsedMilliseconds()
For i = 0 To #NumPoints-1
   With RandomPoint(i)
      p2(Str(\x)+"|"+Str(\y))\p\x = \x
      p2()\p\y = \y
   EndWith
Next

s2b = ElapsedMilliseconds()
NumSquares2 = FindSquares_Map(p2())
t = ElapsedMilliseconds()

t2a = t - s2a
t2b = t - s2b

;------------------------------------------------------------------------------

Structure KDtree_Node
   p.Point
   *l.KDtree_Node      ; Zeiger auf linken  Nachfolger
   *r.KDtree_Node      ; Zeiger auf rechten Nachfolger
   found.i             ; #True/#False
EndStructure


Procedure KDtree_Insert (List tree.KDtree_Node(), *t.KDtree_Node, *p.Point)
   ; -- Hilfsroutine für KDtree_Create()
   ; in : tree(): KD-Baum
   ;      *t: Zeiger auf Anfangsknoten des KD-Baumes
   ;      *p: Zeiger auf einzufügenden Punkt
   ; out: tree(): KD-Baum mit neuem Knoten
   Protected.KDtree_Node *pred
   Protected.i left, d=#True
   
   ; Position für neuen Punkt suchen
   While *t
      *pred = *t
      
      If d
         left = Bool(*t\p\x > *p\x)
      Else
         left = Bool(*t\p\y > *p\y)
      EndIf
      If left
         *t = *t\l
      Else
         *t = *t\r
      EndIf
      d = Bool(Not d)
   Wend
   
   ; Knoten mit neuem Punkt in Liste einfügen
   If left
      *pred\l = AddElement(tree())
   Else
      *pred\r = AddElement(tree())
   EndIf
   CopyStructure(*p, @tree()\p, Point)
   tree()\l = 0
   tree()\r = 0
EndProcedure


Procedure.i KDtree_Create (List tree.KDtree_Node(), Array a.Point(1), x_min.i=0, y_min.i=0)
   ; in : a()  : Array mit 2-dimensionalen Punkten
   ;      x_min: kleinster x-Wert in a()
   ;      y_min: kleinster y-Wert in a()
   ; out: tree()       : KD-Baum aus diesen Punkten
   ;      Funktionswert: Zeiger auf Kopfknoten des KD-Baumes
   Protected.KDtree_Node *head
   Protected.i i, last=ArraySize(a())
   
   ClearList(tree())
   *head = AddElement(tree())
   tree()\p\x = x_min
   tree()\p\y = y_min
   tree()\l = 0
   tree()\r = 0
   
   For i = 0 To last
      KDtree_Insert(tree(), *head, @a(i))
   Next
   
   ProcedureReturn *head
EndProcedure


Procedure.i KDtree_FindPoint (*t.KDtree_Node, *p.Point)
   ; in : *t: Zeiger auf Anfangsknoten des KD-Baumes
   ;      *p: Zeiger auf zu suchenden Punkt
   ; out: Funktionswert: Zeiger auf Knoten mit gefundenem Punkt,
   ;                     oder 0 wenn nicht gefunden
   ;                     (Head wird bei der Suche übersprungen)
   Protected.i left, d=#True
   
   If *t
      Repeat
         If d
            left = Bool(*t\p\x > *p\x)
         Else
            left = Bool(*t\p\y > *p\y)
         EndIf
         If left
            *t = *t\l
         Else
            *t = *t\r
         EndIf
         d = Bool(Not d)
      Until *t = 0 Or (*t\p\x = *p\x And *t\p\y = *p\y)
   EndIf
   
   ProcedureReturn *t
EndProcedure


Procedure.i FindSquares_Tree (List tree.KDtree_Node(), *head.KDtree_Node)
   Protected.KDtree_Node *a, *t
   Protected.Point vektor, b, c, d
   Protected.i count=0
   
   FirstElement(tree())          ; Head überspringen
   While NextElement(tree())     ; für alle Punkte
      tree()\found = #True
      *a = @tree()               ; aktuelles Element der äußeren Schleife
      While NextElement(tree())  ; für alle noch nicht geprüften Punkte
         b = tree()\p
         
         ; Vektor von a nach b ermitteln
         vektor\x = b\x - *a\p\x
         vektor\y = b\y - *a\p\y
         
         ; Punkt c berechnen
         c\x = b\x - vektor\y
         c\y = b\y + vektor\x
         
         ; schauen ob Punkt c in der Liste ist
         *t = KDtree_FindPoint(*head, c)
         If *t And (*t\found = #False)
            ; Punkt d berechnen
            d\x = *a\p\x - vektor\y
            d\y = *a\p\y + vektor\x
            
            ; schauen ob Punkt d in der Liste ist
            *t = KDtree_FindPoint(*head, d)
            If *t And (*t\found = #False)
               count + 1
            EndIf
         EndIf
      Wend
      ChangeCurrentElement(tree(), *a)        ; zurück zum aktuellen Element der äußeren Schleife
   Wend
   
   ProcedureReturn count
EndProcedure


Define *head.KDtree_Node
NewList tree.KDtree_Node()

s3a = ElapsedMilliseconds()
*head = KDtree_Create(tree(), RandomPoint())

s3b = ElapsedMilliseconds()
NumSquares3 = FindSquares_Tree(tree(), *head)
t = ElapsedMilliseconds()

t3a = t - s3a
t3b = t - s3b

;------------------------------------------------------------------------------

Define msg$

;-- Zeiten anzeigen:
msg$ = "FindSquares_Array(): "  + Str(NumSquares1) + " Quadrate in " + Str(t1) + " ms" + #LF$
msg$ + "FindSquares_Map()  : "  + Str(NumSquares2) + " Quadrate in " + Str(t2a) + "/" + Str(t2b) + " ms" + #LF$
msg$ + "FindSquares_Tree() : "  + Str(NumSquares3) + " Quadrate in " + Str(t3a) + "/" + Str(t3b) + " ms"

MessageRequester("Ergebnis", msg$)

Die Initiaisierung der Map bzw. des KD-Baums fällt übrigens zeitlich kaum ins Gewicht, bei meinen Versuchen waren die "Brutto"- und die "Netto"-Zeiten immer gleich. Ich habe nacheinander 500, 750, 1000, 1250 und 1500 Punkte zufällig aus 100x100 Punkten ausgewählt (und das Programm mit jeder Punktzahl 6x laufen lassen). Daraus resultieren folgende Diagramme:





Der Algorithmus mit der Map kann (bei vielen Punkten) etwas beschleunigt werden, wenn man die Anzahl der Slots erhöht. Das bewirkt aber nur einen graduellen Gewinn, keinen wesentlichen Unterschied in diesem Vergleich, bei dem der KD-Baum immer deutlich am schnellsten ist.

Grüße, Nino