Mindestwandstärke eines nicht konvexen Polygons mit Löchern

Was ist der effizienteste Weg, um die minimale Wandstärke (Wert und Position) eines komplexen, nicht konvexen Polygonbereichs einschließlich Löchern zu ermitteln? Siehe Beispiel eines Polygons in Blau mit der minimalen Wandstärke in Rot, obwohl in diesem Fall die Position nicht eindeutig ist, wenn die beiden benachbarten Linien parallel sind.

Bisher haben wir versucht:

• Unterteilen von Polygonlinien und Finden von minimalen Punkt-Punkt-Linien innerhalb des Polygons (Brute Force, nicht effizient für komplexe Polygone mit> 10'000 Punkten)
• Verzögern Sie die Triangulation und finden Sie minimale Kanten innerhalb des Polygons. Nicht präzise genug, nur machbar, wenn zuerst eine Unterteilung der Polygonlinien kombiniert wird. Hier ist ein Beispiel (Nr. 3), in dem die Delaunay-Triangulation keine Simplexkanten in Rot findet, aber die minimale Wandstärke in der grünen Box verpasst:
  
• Iterativ zunehmender Erosionspuffer, um den minimalen Einschub zu finden, bei dem das Erosionspolygon in mehrere Teile zerfällt = die Hälfte der minimalen Wandstärke. Problem ist es, bei diesem Ansatz anschließend den Ort der minimalen Wandstärke zu finden. Darüber hinaus zerfällt die Erosion nicht immer in mehrere Teile und übersieht "Sackgassen". Hier ist ein Beispiel (Nr. 2), das zu einer Linie erodiert und die falsche minimale Wandstärke ergibt:
  
• Suchen Sie zuerst die Mittelachse und suchen Sie dann nach dem Mindestkreis auf der Mittelachse, der den Polygonbereich bedeckt, aber nicht überlappt. Edit: Problematisch sind die vielen "falschen Kandidaten" auf der Mittelachse: ZB. (Nr. 1) Kreis A wäre falsch, Kreis B würde die korrekte Mindestwandstärke bezeichnen:
  

Eine der effizientesten Methoden, um die minimale Wandstärke (Wert und Position) eines komplexen, nicht konvexen Polygonbereichs einschließlich Löchern zu ermitteln, könnte die Verwendung einer regelmäßig beabstandeten Schicht (oder zufälligen) von Punkten zur Bestimmung des ersten nächstgelegenen Segments sein mit Kontext für jeden Punkt und als nächstes den Schnittpunkt zwischen dem inkrementellen Segment und dem Polygon der gegenüberliegenden Seite; basiert auf Direktoren Cosinus.

Inkrementelle Abstände können verwendet werden, bis das erste Segment ein Seitenpolygon (die minimale Wandstärke) erreicht und schneidet.

Um meinen Ansatz auszuprobieren, habe ich Ihr Polygon mit Löchern geklont und eine zufällige Punktebene innerhalb des Polygons mit 100 Punkten erstellt. wie es auf folgendem Bild zu sehen ist:

Der von PyQGIS verwendete Code sieht wie folgt aus:

import math

def azimuth(point1, point2):
   return point1.azimuth(point2) #in degrees

def cosdir_azim(azim):
   azim = math.radians(azim)
   cosa = math.sin(azim)
   cosb = math.cos(azim)
   return cosa,cosb

registry = QgsMapLayerRegistry.instance()    
polygon = registry.mapLayersByName('polygon_with_holes')
point_layer = registry.mapLayersByName('Random_points')
points = [ feat.geometry().asPoint() for feat in point_layer[0].getFeatures() ]
feat_polygon = polygon[0].getFeatures().next()

#producing rings polygons
rings_polygon = feat_polygon.geometry().asPolygon()
segments = []
epsg = point_layer[0].crs().authid()
uri = "LineString?crs=" + epsg + "&field=id:integer""&index=yes"
mem_layer = QgsVectorLayer(uri,
                           'increasing_segments',
                           'memory')
prov = mem_layer.dataProvider()
length = 10
pt2 = 0 
k = 0
while pt2 == 0:
    for i, point in enumerate(points):
        #determining closest distance to vertex or side polygon
        dist1 = feat_polygon.geometry().closestSegmentWithContext(point)[0]
        #determining point with closest distance to vertex or side polygon
        pt = feat_polygon.geometry().closestSegmentWithContext(point)[1]
        cosa, cosb = cosdir_azim(azimuth(pt, point))
        #extending segment in opposite direction based in director cosine and length 
        op_pt  = QgsPoint(point.x() + (length*cosa), point.y() + (length*cosb))
        segments.append([pt,op_pt])
        geom = QgsGeometry.fromPolyline([point,op_pt])
        for ring in rings_polygon:
            geom_ring = QgsGeometry.fromPolyline(ring)
            if geom.intersects(geom_ring):
                pt3 = geom.intersection(geom_ring)
                pt2 = pt3.distance(QgsGeometry.fromPoint(point))
                ms = [pt3.asPoint(), pt]
    length += 100
    k += 1
new_segments = segments[len(segments) -1 - len(segments)/k: len(segments) - 1]

feats = [ QgsFeature() for i in range(len(new_segments)) ]
for i,feat in enumerate(feats):
    feat.setAttributes([i])
    geom = QgsGeometry.fromPolyline(new_segments[i])
    feat.setGeometry(geom)

prov.addFeatures(feats)
QgsMapLayerRegistry.instance().addMapLayer(mem_layer)
minimum_segment = QgsGeometry().fromPolyline(ms).exportToWkt()
print minimum_segment, k

und es erzeugt eine Speicherschicht mit inkrementellen Abständen (nur zu Visualisierungszwecken) und druckt eine minimale Wandstärke im WKT-Format.

Nachdem ich den Code in der Python Console von QGIS ausgeführt habe, habe ich das folgende Bild erhalten:

Es kann beobachtet werden, dass nur ein inkrementeller Abstand zuerst die gegenüberliegende Seite im erwarteten Bereich erreichte.

Das gedruckte WKT-Format (für minimale Wandstärke) wird mit dem QuickWKT-Plugin von QGIS verwendet, um dieses Segment in der folgenden Abbildung zu visualisieren:

Die leichte Neigung wurde erzeugt, weil "engstes Segment mit Kontext" mit einem Scheitelpunkt verbunden war; stattdessen Seitenpolygon. Dies kann jedoch mit einer Code-Ausnahme oder mehr Punkten vermieden werden.

Zwei weitere Ideen, die Sie in den Topf werfen sollten:

1. Rasteren Sie Ihr Polygon und verwenden Sie eine Entfernungstransformation (gibt das Bild der kürzesten Entfernung von jedem Pixel ungleich Null auf ein Nullpixel zurück). Skelettieren Sie Ihr gerastertes Bild und nehmen Sie dann die Werte des entfernungstransformierten Bildes entlang des Skeletts. Von diesem Satz haben Sie einige Mindestwerte, die Ihrer engsten Breite entsprechen sollten. Dieses Set kann als anfänglicher Suchpunkt verwendet werden, um dann Ihren Brute-Force-Ansatz zu implementieren. Ich sollte beachten, dass das Skelett die Ecken von Objekten halbiert und an diesen Stellen die Entfernungstransformation gegen Null geht (wenn Sie sich der Objektgrenze nähern). Dies mag problematisch sein, stellt jedoch ein Problem mit Ihrem Problem dar - warum kann ' t die kleinste Breite an einer Ecke sein (und im Wesentlichen Null sein)? Sie können dieses Problem beheben, indem Sie einen Schwellenwert für den kürzesten Abstand um den Umfang zwischen den beiden Punkten festlegen (wenn sie auf demselben Objekt liegen). Sie können eine geodätische Entfernungstransformation für den Satz von Umfangspixeln verwenden, um diesen Wert schnell zu finden.

   Bei dieser Methode müssen Sie eine Entscheidung über die Auflösung des gerasterten Polygons treffen, was zu einer gewissen Skalierungsabhängigkeit führt. Wenn Sie eine zu hohe Auflösung wählen, kann die Entfernungstransformation zeitaufwändig sein. Aber im Allgemeinen sind sie ziemlich schnell. Diese Methode bietet Ihnen möglicherweise nicht die gewünschte Präzision, aber zumindest eine viel kleinere Anzahl von Positionen, die überprüft werden müssen.

2. Ihre Brute-Force-Methode ist kein schlechter Ausgangspunkt. Ich hatte ein ähnliches Problem, bei dem ich alle Schnittpunkte einer (langen) Linie mit sich selbst finden musste, und ich konnte die Suchzeit mithilfe eines kd-Baum-Suchalgorithmus (ich habe damals die Bereichssuche in Matlab verwendet) erheblich beschleunigen Punkte innerhalb einer Nachbarschaft zuerst. Auf diese Weise erzwingen Sie nur eine kleine Teilmenge der Gesamtzahl der Punkte.

Der Ansatz der Mittelachse ist korrekt. Sie benötigen lediglich ein Kriterium, um fehlerhafte Kreise zu ignorieren: Jeder Kreis auf der Mittelachse berührt die Oberfläche an zwei (oder mehr) Punkten. Stellen Sie sich Vektoren vom Mittelpunkt des Kreises zu diesen Punkten auf der Oberfläche vor und beobachten Sie, dass der Winkel zwischen 180 ° für den guten Kreis B und nur 90 ° für den schlechten Kreis A beträgt.

Ein generischer Polygonzeichnungsalgorithmus sortiert Liniensegmente von oben nach unten und zeichnet orthogonale Pixellinien durch. Da diese Art des Abbaus von Polygonen (ohne Krümmungen) sehr schnell ist, kann sie als Grundlage verwendet werden. Anstatt nur von oben nach unten zu gehen, können Sie 0, 30, 60 und 90 Grad gehen und Ihren kürzesten orthogonalen Abschnitt (= minimale Wandstärke!) Finden. Sie müssen dies nur einmal pro Punkt berechnen und nicht für jede Art von "Pixelauflösung".

Siehe Computergrafik-Scan-Line-Polygon-Fill-Algorithmus