Die konkave Hülle

Erstellen einer Clustergrenze mit einem K-nearest Neighbors-Ansatz

Vor ein paar Monaten schrieb ich hier auf Medium einen Artikel über die Kartierung der Verkehrsunfälle in Großbritannien. Ich war hauptsächlich daran interessiert, die Verwendung des DBSCAN-Clustering-Algorithmus für geografische Daten zu veranschaulichen. In dem Artikel habe ich geografische Informationen verwendet, die von der britischen Regierung zu gemeldeten Verkehrsunfällen veröffentlicht wurden. Mein Ziel war es, einen dichtebasierten Clustering-Prozess durchzuführen, um die Bereiche zu finden, in denen Verkehrsunfälle am häufigsten gemeldet werden. Das Endergebnis war die Schaffung einer Reihe von Geo-Zäunen, die diese Unfallherde darstellen.

Wenn Sie alle Punkte in einem bestimmten Cluster sammeln, können Sie eine Vorstellung davon bekommen, wie der Cluster auf einer Karte aussehen würde, aber Ihnen fehlt eine wichtige Information: die äußere Form des Clusters. In diesem Fall handelt es sich um ein geschlossenes Polygon, das auf einer Karte als Geozaun dargestellt werden kann. Es kann davon ausgegangen werden, dass jeder Punkt innerhalb des Geozauns zum Cluster gehört, was diese Form zu einer interessanten Information macht: Sie können sie als Diskriminatorfunktion verwenden. Es kann davon ausgegangen werden, dass alle neu abgetasteten Punkte, die innerhalb des Polygons liegen, zum entsprechenden Cluster gehören. Wie ich in dem Artikel angedeutet habe, können Sie solche Polygone verwenden, um Ihr Fahrrisiko zu behaupten, indem Sie sie verwenden, um Ihren eigenen abgetasteten GPS-Standort zu klassifizieren.

Die Frage ist nun, wie man aus einer Wolke von Punkten, die einen bestimmten Cluster bilden, ein aussagekräftiges Polygon erstellt. Mein Ansatz im ersten Artikel war etwas naïv und spiegelte eine Lösung wider, die ich bereits im Produktionscode verwendet hatte. Diese Lösung beinhaltete das Platzieren eines Kreises mittig an jedem Punkt des Clusters und das Zusammenführen aller Kreise zu einem wolkenförmigen Polygon. Das Ergebnis ist weder sehr schön noch realistisch. Durch die Verwendung von Kreisen als Basisform zum Erstellen des endgültigen Polygons haben diese mehr Punkte als eine stromlinienförmigere Form, wodurch die Speicherkosten erhöht und die Ausführung der Einschlusserkennungsalgorithmen verlangsamt wird.

Andererseits hat dieser Ansatz den Vorteil, dass er (zumindest aus Sicht des Entwicklers) rechnerisch einfach ist, da er die cascaded_union -Funktion von Shapely verwendet, um alle Kreise zusammenzuführen. Ein weiterer Vorteil ist, dass die Form des Polygons implizit unter Verwendung aller Punkte im Cluster definiert wird.

Für anspruchsvollere Ansätze müssen wir irgendwie die Grenzpunkte des Clusters identifizieren, die die Punktwolkenform zu definieren scheinen. Interessanterweise können Sie mit einigen DBSCAN-Implementierungen die Grenzpunkte tatsächlich als Nebenprodukt des Clusterprozesses wiederherstellen. Leider sind diese Informationen (anscheinend) in der Implementierung von SciKit Learn nicht verfügbar, daher müssen wir auskommen.

Der erste Ansatz, der mir in den Sinn kam, war die Berechnung der konvexen Hülle der Punktmenge. Dies ist ein gut verstandener Algorithmus, leidet jedoch unter dem Problem, konkave Formen wie diese nicht zu behandeln:

Diese Form erfasst die Essenz der zugrunde liegenden Punkte nicht korrekt. Würde es als Diskriminator verwendet, würden einige Punkte fälschlicherweise als innerhalb des Clusters klassifiziert, wenn dies nicht der Fall ist. Wir brauchen einen anderen Ansatz.

Die konkave Rumpfalternative

Glücklicherweise gibt es Alternativen zu diesem Zustand: Wir können einen konkaven Rumpf berechnen. So sieht die konkave Hülle aus, wenn sie auf denselben Punktesatz wie im vorherigen Bild angewendet wird:

Oder vielleicht dieser:

Wie Sie sehen können, gibt es im Gegensatz zum konvexen Rumpf keine einheitliche Definition dessen, was der konkave Rumpf einer Menge von Punkten ist. Mit dem Algorithmus, den ich hier vorstelle, wird die Wahl, wie konkav Ihre Rümpfe sein sollen, über einen einzigen Parameter getroffen: k — die Anzahl der nächsten Nachbarn, die bei der Rumpfberechnung berücksichtigt werden. Mal sehen, wie das funktioniert.

Der Algorithmus

Der Algorithmus, den ich hier vorstelle, wurde vor über einem Jahrzehnt von Adriano Moreira und Maribel Yasmina Santos von der Universität Minho, Portugal, beschrieben . Aus dem Abstract:

Dieses Papier beschreibt einen Algorithmus zur Berechnung der Hüllkurve eines Satzes von Punkten in einer Ebene, der konvexe auf nicht konvexe Rümpfe erzeugt, die den von den gegebenen Punkten belegten Bereich darstellen. Der vorgeschlagene Algorithmus basiert auf einem k-Nearest Neighbors-Ansatz, bei dem der Wert von k, dem einzigen Algorithmusparameter, verwendet wird, um die “Glätte” der endgültigen Lösung zu steuern.

Da ich diesen Algorithmus auf geographische Informationen anwenden werde, mussten einige Änderungen vorgenommen werden, nämlich bei der Berechnung von Winkeln und Entfernungen . Diese ändern jedoch in keiner Weise den Kern des Algorithmus, der durch die folgenden Schritte allgemein beschrieben werden kann:

1. Suchen Sie den Punkt mit der niedrigsten y-Koordinate (Breitengrad) und machen Sie ihn zum aktuellen Punkt.
2. Finden Sie die k-nächstgelegenen Punkte zum aktuellen Punkt.
3. Wählen Sie aus den k-nächsten Punkten den Punkt aus, der der größten Rechtskurve aus dem vorherigen Winkel entspricht. Hier verwenden wir das Konzept der Peilung und beginnen mit einem Winkel von 270 Grad (genau nach Westen).
4. Überprüfen Sie, ob sich der neue Punkt nicht schneidet, indem Sie den neuen Punkt zur wachsenden Linienzeichenfolge hinzufügen. Wenn dies der Fall ist, wählen Sie einen anderen Punkt aus der k-Liste aus oder starten Sie ihn mit einem größeren Wert von k neu.
5. Machen Sie den neuen Punkt zum aktuellen Punkt und entfernen Sie ihn aus der Liste.
6. Fügen Sie nach k Iterationen den ersten Punkt wieder zur Liste hinzu.
7. Schleife zu Nummer 2.

Der Algorithmus scheint recht einfach zu sein, aber es gibt eine Reihe von Details, die beachtet werden müssen, insbesondere weil es sich um geografische Koordinaten handelt. Entfernungen und Winkel werden auf andere Weise gemessen.

Der Code

Hier veröffentliche ich eine angepasste Version des Codes des vorherigen Artikels. Sie finden immer noch denselben Clustercode und denselben wolkenförmigen Clustergenerator. Die aktualisierte Version enthält jetzt ein Paket mit dem Namen geomath.hulls, in dem Sie die Klasse ConcaveHull finden. Um Ihre konkaven Rümpfe zu erstellen, gehen Sie wie folgt vor:

Im obigen Code ist points ein Array von Dimensionen (N, 2), wobei die Zeilen die beobachteten Punkte und die Spalten die geografischen Koordinaten (Längen- und Breitengrad) enthalten. Das resultierende Array hat genau die gleiche Struktur, enthält jedoch nur die Punkte, die zur polygonalen Form des Clusters gehören. Eine Art Filter.

Da wir mit Arrays umgehen, ist es nur natürlich, NumPy in den Kampf zu bringen. Alle Berechnungen wurden nach Möglichkeit ordnungsgemäß vektorisiert, und es wurden Anstrengungen unternommen, um die Leistung beim Hinzufügen und Entfernen von Elementen aus Arrays zu verbessern (Spoiler: Sie werden überhaupt nicht verschoben). Eine der fehlenden Verbesserungen ist die Code-Parallelisierung. Aber das kann warten.

Ich habe den Code wie in der Arbeit beschrieben um den Algorithmus herum organisiert, obwohl während der Übersetzung einige Optimierungen vorgenommen wurden. Der Algorithmus basiert auf einer Reihe von Unterprogrammen, die durch das Papier eindeutig identifiziert werden, also lassen Sie uns diese jetzt aus dem Weg räumen. Für Ihren Lesekomfort, Ich werde die gleichen Namen wie in der Zeitung verwenden.

CleanList -Die Reinigung der Liste der Punkte wird im Klassenkonstruktor durchgeführt:

Wie Sie sehen können, wird die Liste der Punkte aus Leistungsgründen als NumPy-Array implementiert. Die Bereinigung der Liste erfolgt in Zeile 10, in der nur die eindeutigen Punkte gespeichert werden. Das Datensatzarray ist mit Beobachtungen in Zeilen und geografischen Koordinaten in den beiden Spalten organisiert. Beachten Sie, dass ich in Zeile 13 auch ein boolesches Array erstelle, das zum Indizieren in das Hauptdatensatzarray verwendet wird, um das Entfernen von Elementen und das gelegentliche Hinzufügen von Elementen zu erleichtern. Ich habe diese Technik in der NumPy-Dokumentation als “Maske” bezeichnet und sie ist sehr leistungsfähig. Was die Primzahlen betrifft, werde ich sie später besprechen.

FindMinYPoint – Dies erfordert eine kleine Funktion:

Diese Funktion wird mit dem Dataset-Array als Argument aufgerufen und gibt den Index des Punktes mit dem niedrigsten Breitengrad zurück. Beachten Sie, dass Zeilen mit dem Längengrad in der ersten Spalte und dem Breitengrad in der zweiten Spalte codiert sind.

RemovePoint
AddPoint – Dies sind aufgrund der Verwendung des Arrays indices keine Probleme. Dieses Array wird verwendet, um die aktiven Indizes im Hauptdatensatzarray zu speichern, sodass das Entfernen von Elementen aus dem Datensatz ein Kinderspiel ist.

Obwohl der in der Arbeit beschriebene Algorithmus das Hinzufügen eines Punkts zu dem Array erfordert, aus dem der Rumpf besteht, wird dies tatsächlich als implementiert:

Später wird die Variable test_hull wieder hull zugewiesen, wenn die Linienzeichenfolge als nicht überschneidend angesehen wird. Aber ich bin immer vor dem Spiel hier. Das Entfernen eines Punktes aus dem Dataset-Array ist so einfach wie:

self.indices = False

Wenn Sie es wieder hinzufügen, müssen Sie nur den Array-Wert im selben Index auf true zurücksetzen. Aber all diese Bequemlichkeit kommt mit dem Preis, dass wir die Indizes im Auge behalten müssen. Dazu später mehr.

NearestPoints – Hier wird es interessant, weil wir es nicht mit planaren Koordinaten zu tun haben, also raus mit Pythagoras und rein mit Haversine:

Beachten Sie, dass der zweite und dritte Parameter Arrays im Datensatzformat sind: Längengrad in der ersten Spalte und Breitengrad in der zweiten Spalte. Wie Sie sehen, gibt diese Funktion ein Array von Entfernungen in Metern zwischen dem Punkt im zweiten Argument und den Punkten im dritten Argument zurück. Sobald wir diese haben, können wir die k-nächsten Nachbarn auf einfache Weise erhalten. Dafür gibt es jedoch eine spezielle Funktion, die einige Erklärungen verdient:

Die Funktion beginnt mit der Erstellung eines Arrays mit den Basisindizes. Dies sind die Indizes der Punkte, die nicht aus dem Datensatzarray entfernt wurden. Wenn wir zum Beispiel in einem Zehn-Punkte-Cluster mit dem Entfernen des ersten Punkts beginnen würden, wäre das Array der Basisindizes . Als nächstes berechnen wir die Abstände und sortieren die resultierenden Array-Indizes. Die ersten k werden extrahiert und dann als Maske verwendet, um die Basisindizes abzurufen. Es ist irgendwie verzerrt, funktioniert aber. Wie Sie sehen, gibt die Funktion kein Array von Koordinaten zurück, sondern ein Array von Indizes in das Datensatzarray.

SortByAngle – Hier gibt es mehr Probleme, da wir keine einfachen Winkel berechnen, sondern Lager. Diese werden als Null Grad genau nach Norden gemessen, wobei die Winkel im Uhrzeigersinn zunehmen. Hier ist der Kern des Codes, der die Lager berechnet:

Die Funktion gibt ein Array von Peilungen zurück, gemessen von dem Punkt, dessen Index im ersten Argument steht, bis zu den Punkten, deren Indizes im dritten Argument stehen. Sortieren ist einfach:

Zu diesem Zeitpunkt enthält das Kandidatenarray die Indizes der k-nächsten Punkte, sortiert nach absteigender Reihenfolge der Peilung.

IntersectQ – Anstatt meine eigenen Linienkreuzungsfunktionen zu rollen, wandte ich mich an Shapely, um Hilfe zu erhalten. Tatsächlich behandeln wir beim Erstellen des Polygons im Wesentlichen eine Linienzeichenfolge und hängen Segmente an, die sich nicht mit den vorherigen schneiden. Das Testen ist einfach: Wir nehmen das im Bau befindliche Rumpf-Array auf, konvertieren es in ein formschönes Linien-String-Objekt und testen, ob es einfach ist (nicht selbstschneidend) oder nicht.

Kurz gesagt, eine formschöne Linienzeichenfolge wird komplex, wenn sie sich selbst kreuzt, sodass das Prädikat is_simple falsch wird. Einfach.

PointInPolygon — Dies erwies sich als das am schwierigsten zu implementierende. Erlauben Sie mir, dies anhand des Codes zu erklären, der die endgültige Validierung des Hull-Polygons durchführt (überprüft, ob alle Punkte des Clusters im Polygon enthalten sind):

Shapelys Funktionen zum Testen auf Schnittpunkt und Einbeziehung hätten ausreichen sollen, um zu überprüfen, ob das endgültige Hüllenpolygon alle Punkte des Clusters überlappt, dies war jedoch nicht der Fall. Warum? Shapely ist koordinatenunabhängig und behandelt geografische Koordinaten, die in Breiten- und Längengraden ausgedrückt werden, genauso wie Koordinaten auf einer kartesischen Ebene. Aber die Welt verhält sich anders, wenn Sie auf einer Kugel leben, und Winkel (oder Lager) sind entlang einer geodätischen nicht konstant. Das Beispiel einer geodätischen Linie, die Bagdad mit Osaka verbindet, veranschaulicht dies perfekt. Es kommt also vor, dass der Algorithmus unter bestimmten Umständen einen Punkt basierend auf dem Peilungskriterium enthalten kann, aber später, unter Verwendung der planaren Algorithmen von Shapely, als etwas außerhalb des Polygons angesehen wird. Das ist es, was die kleine Entfernungskorrektur dort macht.

Ich habe eine Weile gebraucht, um das herauszufinden. Meine Debugging-Hilfe war QGIS, ein großartiges Stück freier Software. Bei jedem Schritt der verdächtigen Berechnungen würde ich die Daten im WKT-Format in eine CSV-Datei ausgeben, um sie als Layer einzulesen. Ein echter Lebensretter!

Wenn das Polygon nicht alle Punkte des Clusters abdeckt, besteht die einzige Möglichkeit darin, k zu erhöhen und es erneut zu versuchen. Hier habe ich ein bisschen meiner eigenen Intuition hinzugefügt.

Prime k

Der Artikel schlägt vor, den Wert von k um eins zu erhöhen und den Algorithmus erneut von Grund auf neu auszuführen. Meine frühen Tests mit dieser Option waren nicht sehr zufriedenstellend: Die Laufzeiten auf problematischen Clustern wären langsam. Dies war auf den langsamen Anstieg von k zurückzuführen, daher beschloss ich, einen anderen Erhöhungsplan zu verwenden: eine Tabelle mit Primzahlen. Der Algorithmus beginnt bereits mit k = 3 , daher war es eine einfache Erweiterung, ihn auf einer Liste von Primzahlen zu entwickeln. Dies ist, was Sie im rekursiven Aufruf sehen:

Ich habe eine Sache für Primzahlen, wissen Sie …

Blow Up

Die von diesem Algorithmus erzeugten konkaven Hüllenpolygone müssen noch weiter verarbeitet werden, da sie nur Punkte innerhalb des Rumpfes unterscheiden, aber nicht in der Nähe davon. Die Lösung besteht darin, diesen dünnen Clustern etwas Polsterung hinzuzufügen. Hier verwende ich genau die gleiche Technik wie zuvor, und so sieht es aus:

Hier habe ich Shapelys buffer -Funktion verwendet, um den Trick zu machen.

Die Funktion akzeptiert ein formschönes Polygon und gibt eine aufgeblasene Version von sich selbst zurück. Der zweite Parameter ist der Radius in Metern der hinzugefügten Auffüllung.

Ausführen des Codes

Ziehen Sie zunächst den Code aus dem GitHub-Repository auf Ihren lokalen Computer. Die Datei, die Sie ausführen möchten, befindet sich ShowHotSpots.py im Hauptverzeichnis. Bei der ersten Ausführung liest der Code die britischen Verkehrsunfalldaten von 2013 bis 2016 ein und clustert sie. Die Ergebnisse werden dann als CSV-Datei für nachfolgende Läufe zwischengespeichert.

Sie erhalten dann zwei Karten: Die erste wird mit den wolkenförmigen Clustern generiert, während die zweite den hier diskutierten konkaven Clustering-Algorithmus verwendet. Während der Polygongenerierungscode ausgeführt wird, werden möglicherweise einige Fehler gemeldet. Um zu verstehen, warum der Algorithmus keine konkave Hülle erstellt, schreibt der Code die Cluster in CSV-Dateien in das Verzeichnis data/out/failed/. Wie üblich können Sie diese Dateien mit QGIS als Layer importieren.

Im Wesentlichen schlägt dieser Algorithmus fehl, wenn er nicht genügend Punkte findet, um die Form zu “umgehen”, ohne sich selbst zu schneiden. Dies bedeutet, dass Sie bereit sein müssen, diese Cluster entweder zu verwerfen oder eine andere Behandlung auf sie anzuwenden (konvexe Hülle oder koaleszierte Blasen).

Konkavität

Es ist ein Wrap. In diesem Artikel habe ich eine Methode zur Nachbearbeitung von DBSCAN-generierten geografischen Clustern in konkave Formen vorgestellt. Diese Methode kann im Vergleich zu anderen Alternativen ein besser passendes Außenpolygon für die Cluster bereitstellen.

Vielen Dank fürs Lesen und viel Spaß beim Basteln mit dem Code!

Kryszkiewicz M., Lasek P. (2010) TI-DBSCAN: Clustering mit DBSCAN mittels der Dreiecksungleichung. In: Szczuka M., Kryszkiewicz M., Ramanna S., Jensen R., Hu Q. (Hrsg.) Rough Sets and Current Trends in Computing. RSCTC 2010. Vorlesungsnotizen in der Informatik, vol 6086. Springer, Berlin, Heidelberg

Scikit-learn: Maschinelles Lernen in Python, Pedregosa et al., JMLR 12, S. 2825-2830, 2011

Moreira, A. und Santos, M.Y., 2007, Concave Hull: A K-nearest neighbors approach for the computation of the region occupied by a set of points

Berechnen Sie Abstand, Peilung und mehr zwischen Breiten- /Längengradpunkten

GitHub repository