Geekly articles weekly

Fortune-Algorithmus, Implementierungsdetails

In den letzten Wochen habe ich an der Implementierung des Fortune-Algorithmus in C ++ gearbeitet. Dieser Algorithmus verwendet viele 2D-Punkte und erstellt daraus ein Voronoi-Diagramm . Wenn Sie nicht wissen, was ein Voronoi-Diagramm ist, sehen Sie sich die Abbildung an:


Für jeden Einstiegspunkt, der als „Site“ bezeichnet wird, müssen wir viele Punkte finden, die näher an diesem Ort liegen als an allen anderen. Solche Punktsätze bilden die im obigen Bild gezeigten Zellen.

Es ist bemerkenswert in Fortunes Algorithmus, dass er solche Diagramme rechtzeitig erstellt O(n logn) (was für einen Vergleichsalgorithmus optimal ist), wobei n Ist die Anzahl der Plätze.

Ich schreibe diesen Artikel, weil ich die Implementierung dieses Algorithmus für eine sehr schwierige Aufgabe halte. Im Moment ist dies der komplizierteste der Algorithmen, die ich implementieren musste. Daher möchte ich die Probleme, auf die ich gestoßen bin, und deren Lösung teilen.

Wie üblich wird der Code auf github veröffentlicht , und alle von mir verwendeten Referenzmaterialien sind am Ende des Artikels aufgeführt.

Beschreibung des Fortune-Algorithmus


Ich werde nicht erklären, wie der Algorithmus funktioniert, weil andere Leute es bereits gut gemacht haben. Ich kann empfehlen, diese beiden Artikel zu studieren: hier und hier . Das zweite ist sehr interessant - der Autor hat eine interaktive Demo in Javascript geschrieben, die zum Verständnis der Funktionsweise des Algorithmus hilfreich ist. Wenn Sie einen formelleren Ansatz mit allen Beweisen benötigen, empfehle ich, Kapitel 7 der Computational Geometry, 3. Ausgabe, zu lesen.

Außerdem beschäftige ich mich lieber mit Implementierungsdetails, die nicht gut dokumentiert sind. Und sie machen die Implementierung des Algorithmus so komplex. Insbesondere werde ich mich auf die verwendeten Datenstrukturen konzentrieren.

Ich habe gerade einen Pseudocode des Algorithmus geschrieben, damit Sie eine Vorstellung von seiner globalen Struktur bekommen:

  Fügen Sie der Ereigniswarteschlange für jeden Ort ein Ortsereignis hinzu
 bis die Ereigniswarteschlange leer ist
     Rufen Sie das Top-Ereignis ab
     wenn das Ereignis ein Ortsereignis ist
         Fügen Sie einen neuen Bogen in die Küste ein
         Suchen Sie nach neuen Kreisereignissen
     sonst
         Erstellen Sie einen Scheitelpunkt im Diagramm
         Wir entfernen einen festgezogenen Bogen von der Küste
         ungültige Ereignisse löschen
         Suchen Sie nach neuen Kreisereignissen

Diagrammdatenstruktur


Das erste Problem, auf das ich stieß, war die Wahl der Art und Weise, wie das Voronoi-Diagramm gespeichert werden soll.

Ich entschied mich für eine Datenstruktur, die in der Computergeometrie weit verbreitet ist und als doppelt verbundene Kantenliste (DCEL) bezeichnet wird.

Meine VoronoiDiagram Klasse verwendet vier Container als Felder:

 class VoronoiDiagram { public: // ... private: std::vector<Site> mSites; std::vector<Face> mFaces; std::list<Vertex> mVertices; std::list<HalfEdge> mHalfEdges; }

Ich werde über jeden von ihnen ausführlich sprechen.

Die Site Klasse beschreibt den Einstiegspunkt. Jeder Ort hat einen Index, der nützlich ist, um ihn in die Warteschlange zu stellen, Koordinaten und einen Zeiger auf die Zelle ( face ):

 struct Site { std::size_t index; Vector2 point; Face* face; };

Die Eckpunkte der Zelle werden durch die Vertex Klasse dargestellt, sie haben nur ein Koordinatenfeld:

 struct Vertex { Vector2 point; };

Hier ist die Implementierung der Halbkanten:

 struct HalfEdge { Vertex* origin; Vertex* destination; HalfEdge* twin; Face* incidentFace; HalfEdge* prev; HalfEdge* next; };

Sie fragen sich vielleicht, was ist eine halbe Rippe? Eine Kante im Voronoi-Diagramm ist zwei benachbarten Zellen gemeinsam. In der DCEL-Datenstruktur teilen wir diese Kanten in zwei Halbkanten, eine für jede Zelle, und sie sind durch einen Doppelzeiger verbunden. Darüber hinaus hat die Halbkante einen Start- und einen Endpunkt. Das incidentFace Feld gibt die Fläche an, zu der die Halbkante gehört. Zellen in DCEL werden als zyklische doppelt verknüpfte Liste von Halbkanten implementiert, in der benachbarte Halbkanten miteinander verbunden sind. Daher geben die Felder prev und next die vorherige und die nächste Halbkante in der Zelle an.

Die folgende Abbildung zeigt alle diese Felder für die rote Halbkante:


Schließlich definiert die Face Klasse die Zelle. Es enthält einfach einen Zeiger auf seinen Platz und einen anderen auf eine seiner Halbrippen. Es spielt keine Rolle, welche der Halbkanten ausgewählt ist, da die Zelle ein geschlossenes Polygon ist. Auf diese Weise erhalten wir Zugriff auf alle Halbkanten, während wir eine zyklisch verknüpfte Liste durchlaufen.

 struct Face { Site* site; HalfEdge* outerComponent; };

Ereigniswarteschlange


Die Standardmethode zum Implementieren einer Ereigniswarteschlange ist eine Prioritätswarteschlange. Bei der Verarbeitung von Orts- und Kreisereignissen müssen wir möglicherweise Kreisereignisse aus der Warteschlange entfernen, da sie nicht mehr gültig sind. Bei den meisten Implementierungen von Warteschlangen mit Standardpriorität können Sie jedoch kein Element löschen, das sich nicht oben befindet. Dies gilt insbesondere für std::priority_queue .

Es gibt zwei Möglichkeiten, um dieses Problem zu lösen. Die erste, einfachere Möglichkeit besteht darin, Ereignissen ein valid Flag hinzuzufügen. valid wird zunächst auf true . Anstatt das Kreisereignis aus der Warteschlange zu entfernen, können wir sein Flag einfach auf false . Wenn wir schließlich alle Ereignisse in der Hauptschleife verarbeiten, prüfen wir, ob der valid Flag-Wert des Ereignisses false ist. Wenn ja, überspringen Sie ihn einfach und verarbeiten Sie den nächsten.

Die zweite Methode, die ich angewendet habe, war, nicht std::priority_queue . Stattdessen habe ich meine eigene Prioritätswarteschlange implementiert, die das Entfernen aller darin enthaltenen Elemente unterstützt. Die Implementierung einer solchen Warteschlange ist recht einfach. Ich habe diese Methode gewählt, weil sie den Algorithmuscode sauberer macht.

Küste


Die Küstenlinien-Datenstruktur ist ein komplexer Teil des Algorithmus. Bei falscher Implementierung gibt es keine Garantie dafür, dass der Algorithmus in ausgeführt wird O(n logn) . Der Schlüssel, um diese zeitliche Komplexität zu erreichen, ist die Verwendung eines selbstausgleichenden Baums. Aber es ist leichter gesagt als getan!

Den meisten Ressourcen, die ich untersucht habe (die beiden oben genannten Artikel und das Buch Computational Geometry ), wird empfohlen, die Küstenlinie als Baum zu implementieren, in dem interne Knoten Bruchpunkte und Blätter Bögen anzeigen. Aber sie sagen nichts darüber, wie man einen Baum balanciert. Ich denke, dass ein solches Modell nicht das beste ist, und hier ist der Grund:

  • Es enthält redundante Informationen: Wir wissen, dass es einen Unterbrechungspunkt zwischen zwei benachbarten Bögen gibt, daher ist es nicht erforderlich, diese Punkte als Knoten darzustellen
  • es ist nicht ausreichend für den Selbstausgleich: Nur der durch Bruchstellen gebildete Teilbaum kann ausgeglichen werden. Wir können wirklich nicht den gesamten Baum ausbalancieren, da die Bögen sonst zu internen Knoten und Blättern von Haltepunkten werden können. Das Schreiben eines Algorithmus, um nur den von internen Knoten gebildeten Teilbaum auszugleichen, scheint mir ein Albtraum zu sein.

Deshalb habe ich beschlossen, die Küste anders darzustellen. In meiner Implementierung ist die Küste auch ein Baum, aber alle Knoten sind Bögen. Ein solches Modell weist keinen der aufgeführten Nachteile auf.

Hier ist die Definition von Arc Arc in meiner Implementierung:

 struct Arc { enum class Color{RED, BLACK}; // Hierarchy Arc* parent; Arc* left; Arc* right; // Diagram VoronoiDiagram::Site* site; VoronoiDiagram::HalfEdge* leftHalfEdge; VoronoiDiagram::HalfEdge* rightHalfEdge; Event* event; // Optimizations Arc* prev; Arc* next; // Only for balancing Color color; };

Die ersten drei Felder dienen zur Strukturierung des Baums. Das Feld leftHalfEdge gibt die Halbkante an, die vom leftHalfEdge linken Punkt des Bogens gezeichnet wird. Und rightHalfEdge befindet sich an der halben Kante, die vom äußersten rechten Punkt gezeichnet wird. Zwei Zeiger, prev und next werden verwendet, um direkten Zugriff auf den vorherigen und nächsten Bogen der Küste zu erhalten. Darüber hinaus können Sie die Küste als doppelt verknüpfte Liste umgehen. Schließlich hat jeder Bogen eine Farbe, mit der die Küstenlinie ausgeglichen wird.

Um die Küste auszugleichen, entschied ich mich für ein rot-schwarzes Schema . Beim Schreiben von Code hat mich das Buch Einführung in Algorithmen inspiriert. Kapitel 13 beschreibt zwei interessante Algorithmen, insertFixup und deleteFixup , die den Baum nach dem Einfügen oder Löschen ausgleichen.

Ich kann die im Buch gezeigte insert jedoch nicht verwenden, da die Schlüssel verwendet werden, um die Einfügemarke des Knotens zu finden. Der Fortune-Algorithmus enthält keine Schlüssel. Wir wissen nur, dass wir einen Bogen vor oder nach dem anderen in die Küste einfügen müssen. Um dies zu implementieren, habe ich die insertAfter insertBefore und insertAfter :

 void Beachline::insertBefore(Arc* x, Arc* y) { // Find the right place if (isNil(x->left)) { x->left = y; y->parent = x; } else { x->prev->right = y; y->parent = x->prev; } // Set the pointers y->prev = x->prev; if (!isNil(y->prev)) y->prev->next = y; y->next = x; x->prev = y; // Balance the tree insertFixup(y); }

Das Einfügen von y vor x erfolgt in drei Schritten:

  1. Suchen Sie einen Platz zum Einfügen eines neuen Knotens. Zu diesem Zweck habe ich die folgende Beobachtung verwendet: Das linke Kind x oder das rechte Kind x->prev ist Nil , und dasjenige, das Nil ist, steht vor x und nach x->prev .
  2. Innerhalb der Küste behalten wir die Struktur einer doppelt verknüpften Liste bei, daher müssen wir die prev und next Zeiger der Elemente x->prev , y und x entsprechend aktualisieren.
  3. Schließlich rufen wir einfach die im Buch beschriebene Methode insertFixup auf, um den Baum auszugleichen.

insertAfter wird ähnlich implementiert.

Die Entfernungsmethode aus dem Buch kann ohne Änderungen implementiert werden.

Diagrammlimit


Hier ist die Ausgabe des oben beschriebenen Fortune-Algorithmus:


Es gibt ein kleines Problem mit einigen Kanten von Zellen am Rand des Bildes: Sie werden nicht gezeichnet, weil sie unendlich sind.

Schlimmer noch, eine Zelle ist möglicherweise kein einzelnes Fragment. Wenn wir beispielsweise drei Punkte auf derselben Linie nehmen, hat der Mittelpunkt zwei unendliche Halbkanten, die nicht miteinander verbunden sind. Dies passt nicht sehr gut zu uns, da wir nicht auf eine der halben Kanten zugreifen können, da die Zelle eine verknüpfte Liste von Kanten ist.

Um diese Probleme zu lösen, werden wir das Diagramm einschränken. Damit meine ich, dass wir jede Zelle des Diagramms so einschränken, dass sie keine unendlichen Kanten mehr hat und jede Zelle ein geschlossenes Polygon ist.

Glücklicherweise können wir mit dem Fortune-Algorithmus schnell endlose Kanten finden: Sie entsprechen Halbkanten, die sich am Ende des Algorithmus noch an der Küste befinden.

Mein Restriktionsalgorithmus erhält eine Box als Eingabe und besteht aus drei Schritten:

  1. Es ermöglicht die Platzierung jedes Scheitelpunkts des Diagramms innerhalb des Rechtecks.
  2. Schneiden Sie jede unendliche Kante ab.
  3. Schließt Zellen.

Stufe 1 ist trivial - wir müssen das Rechteck nur erweitern, wenn es keinen Scheitelpunkt enthält.

Stufe 2 ist ebenfalls recht einfach: Sie besteht aus der Berechnung der Schnittpunkte zwischen den Strahlen und dem Rechteck.

Stufe 3 ist auch nicht sehr kompliziert, erfordert nur Aufmerksamkeit. Ich führe es in zwei Stufen durch. Zuerst füge ich die Eckpunkte des Rechtecks ​​zu den Zellen hinzu, an deren Eckpunkten sie sein sollten. Dann stelle ich sicher, dass alle Eckpunkte der Zelle durch Halbkanten verbunden sind.

Ich empfehle Ihnen, den Code zu studieren und Fragen zu stellen, wenn Sie Details zu diesem Teil benötigen.

Hier ist das Ausgabediagramm des Begrenzungsalgorithmus:


Jetzt sehen wir, dass alle Kanten gezeichnet sind. Und wenn Sie verkleinern, können Sie sicherstellen, dass alle Zellen geschlossen sind:


Schnittpunkt mit Rechteck


Großartig! Aber das erste Bild vom Anfang des Artikels ist besser, oder?

In vielen Anwendungen ist es nützlich, den Schnittpunkt zwischen dem Voronoi-Diagramm und dem Rechteck zu haben, wie im ersten Bild gezeigt.

Das Gute ist, dass es nach dem Einschränken des Diagramms viel einfacher ist, dies zu tun. Die schlechte Nachricht ist, dass wir vorsichtig sein müssen, obwohl der Algorithmus nicht sehr kompliziert ist.

Die Idee ist folgende: Wir gehen um die halbe Kante jeder Zelle herum und überprüfen den Schnittpunkt zwischen der halben Kante und dem Rechteck. Fünf Fälle sind möglich:

  1. Die Halbrippe befindet sich vollständig innerhalb des Rechtecks: Wir speichern eine solche Halbrippe
  2. Die Halbrippe befindet sich vollständig außerhalb des Rechtecks: Wir werfen eine solche Halbrippe weg
  3. Die Halbrippe geht außerhalb des Rechtecks: Wir schneiden die Halbrippe ab und speichern sie als die letzte Halbrippe, die ausgeht .
  4. Die Halbrippe geht in das Rechteck: Wir schneiden die Halbrippe ab, um sie mit der letzten Halbrippe zu verbinden , die ausgegangen ist (wir speichern sie in Fall 3 oder 5).
  5. Die Halbrippe kreuzt das Rechteck zweimal: Wir schneiden die Halbrippe ab und fügen eine Halbrippe hinzu, um sie mit der letzten Halbrippe zu verknüpfen , die ausgeht , und speichern sie dann als die neue letzte Halbrippe, die ausgeht .

Ja, es gab viele Fälle. Ich habe ein Bild erstellt, um sie alle zu zeigen:


Das orangefarbene Polygon ist die ursprüngliche Zelle und das rote ist die abgeschnittene Zelle. Abgeschnittene Halbrippen sind rot markiert. Grüne Rippen wurden hinzugefügt, um die in das Rechteck eintretenden Halbrippen mit den herauskommenden Halbrippen zu verbinden.

Wenn wir diesen Algorithmus auf ein begrenztes Diagramm anwenden, erhalten wir das erwartete Ergebnis:


Fazit


Der Artikel erwies sich als ziemlich lang. Und ich bin sicher, dass Ihnen viele Aspekte immer noch nicht klar sind. Trotzdem hoffe ich, dass es Ihnen nützlich sein wird. Untersuchen Sie den Code auf Details.

Um zusammenzufassen und sicherzustellen, dass wir keine Zeit umsonst verschwendeten, habe ich auf meinem (billigen) Laptop die Zeit gemessen, um das Voronoi-Diagramm für eine andere Anzahl von Orten zu berechnen:

  • n=1000 : 2 ms
  • n=$10.00 : 33 ms
  • n=$100.00 : 450 ms
  • n=$1.000.00 : 6600 ms

Ich habe nichts mit diesen Indikatoren zu vergleichen, aber es scheint, dass es unglaublich schnell ist!

Referenzen


Source: https://habr.com/ru/post/de430628/

More articles:


All Articles