Espesor mínimo de pared de un polígono no convexo con agujeros

Question

¿Cuál es la forma más eficiente de encontrar el espesor mínimo de pared (valor y ubicación) de un área poligonal compleja no convexa que incluye agujeros? Véase el ejemplo de un polígono en azul, con el espesor de pared mínimo en rojo, aunque en este caso la ubicación es ambigua, si las dos líneas adyacentes son paralelas.

Hasta ahora, lo hemos intentado:

• Subdivisión de las líneas del polígono y búsqueda de las líneas mínimas de puntos dentro del polígono (fuerza bruta, no es eficiente para polígonos complejos con >10'000 puntos)
• Triangulación de Delaunay y búsqueda de aristas mínimas dentro del polígono. No es lo suficientemente preciso, sólo es factible si se combina con la subdivisión de las líneas del polígono primero. Aquí está el ejemplo (nº 3), donde la triangulación de Delaunay no encontraría las aristas simplex en rojo, sino que se perdería el grosor mínimo de la pared en la caja verde:
  
• Incremento iterativo del buffer de erosión para encontrar la inserción mínima, donde el polígono de erosión se rompe en múltiples partes = la mitad del espesor de pared mínimo. El problema es encontrar después la ubicación del espesor de pared mínimo con este enfoque. Además, la erosión no siempre se rompe en múltiples partes y pasa por alto los "callejones sin salida". Aquí hay un ejemplo (nº 2) que erosiona hasta una línea y da el grosor mínimo de pared equivocado:
  
• Primero se encuentra el eje medial y luego se busca el círculo mínimo en el eje medial que cubra el área del polígono pero que no se superponga. Edición: El problema es que hay muchos "candidatos erróneos" en el eje medio: Por ejemplo, el círculo A (nº 1) sería erróneo, el círculo B denotaría el grosor de pared mínimo correcto:
  

Answer 1

Dos ideas más para echar al bote:

1. Rasteriza tu polígono y utiliza una transformación de distancia (devuelve la imagen de la distancia más corta de cada píxel no nulo a un píxel nulo). Esqueletar su imagen rasterizada, luego tomar los valores de la imagen transformada de distancia a lo largo del esqueleto. De ese conjunto tendrá unos mínimos que deberían corresponder a su anchura más estrecha. Este conjunto puede ser utilizado como puntos de búsqueda iniciales para luego implementar su enfoque de fuerza bruta. Hay que tener en cuenta que el esqueleto bisecará las esquinas de los objetos, y en esos lugares la transformada de distancia se acercará a cero (a medida que se acerque al límite del objeto). Esto puede ser problemático, pero representa una cuestión con su problema - ¿por qué no puede ser el ancho más pequeño en una esquina (y ser esencialmente cero)? Usted podría abordar este problema mediante el establecimiento de un umbral en la distancia más corta alrededor del perímetro entre los dos puntos (si se encuentran en el mismo objeto). Puede utilizar una transformada de distancia geodésica en el conjunto de píxeles del perímetro para encontrar rápidamente ese valor.

   Este método requiere que se tome una decisión sobre la resolución del polígono rasterizado, lo que introduce cierta dependencia de la escala. Y si se elige una resolución demasiado alta, la transformación de distancia puede llevar mucho tiempo. Pero en general son bastante rápidas. Este método podría no darle la precisión que desea, pero al menos podría darle un conjunto mucho más pequeño de ubicaciones que necesitan ser comprobadas.

2. Tu método de fuerza bruta no es un mal punto de partida. Tuve un problema similar en el que tenía que encontrar todas las intersecciones de una línea (larga) con ella misma, y pude acelerar enormemente el tiempo de búsqueda utilizando un algoritmo de búsqueda de árbol kd (en ese momento usaba rangesearch en Matlab) para encontrar primero los puntos dentro de una vecindad. De este modo, sólo se fuerza bruscamente un pequeño subconjunto del número total de puntos.

Answer 2

Uno de los métodos más eficientes para encontrar el espesor de pared mínimo (valor y ubicación) de un área de polígono complejo, no convexo, incluyendo agujeros, podría ser mediante el uso de una capa regularmente espaciada (o aleatoria) de puntos para determinar, primero, el segmento más cercano con el contexto para cada punto y, a continuación, el punto de intersección entre el segmento incremental y el polígono de lado opuesto; basado en directores cosenos.

Se pueden utilizar distancias incrementales hasta que el primer segmento alcance e intersecte algún polígono lateral (el espesor mínimo de la pared).

Para probar mi enfoque cloné su polígono con agujeros y creé una capa de puntos aleatorios dentro del polígono con 100 puntos; como se puede observar en la siguiente imagen:

El código utilizado por PyQGIS tiene el siguiente aspecto:

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

y produce una capa de memoria de distancias incrementales (sólo para fines de visualización) e imprime el espesor mínimo de la pared en formato WKT.

Después de ejecutar el código en la consola de Python de QGIS obtuve el resultado de la siguiente imagen:

Se puede observar que sólo una distancia incremental llegó primero al lado opuesto en el área esperada.

El formato WKT impreso (para el espesor mínimo de la pared) se utiliza con el plugin QuickWKT de QGIS para visualizar ese segmento en la siguiente imagen:

La ligera inclinación se produjo porque el "segmento más cercano con contexto" se unió a un vértice; en lugar de polígono lateral. Sin embargo, se puede evitar con una excepción de código o más puntos.

Answer 3

El enfoque del eje medial es correcto, sólo se necesita un criterio para ignorar los círculos malos: Cada círculo en el eje medial toca la superficie en dos (o más) puntos. Imagina vectores desde el centro del círculo a estos puntos de la superficie y observa que el ángulo entre ellos es de 180° para el círculo bueno B y sólo de 90° para el círculo malo A.

Answer 4

Un algoritmo genérico de dibujo de polígonos funciona ordenando los segmentos de línea de arriba a abajo, y trabajando a través de ellos para dibujar líneas ortogonales de píxeles. Como esta forma de descomponer los polígonos (sin curvaturas) es muy rápida, puede utilizarse como base. En lugar de ir sólo de arriba a abajo, puedes ir a 0, 30, 60 y 90 grados, y encontrar tu sección ortogonal más corta (¡espesor mínimo de la pared!), sólo tienes que calcular eso una vez por punto y no para cualquier tipo de "resolución de píxeles".

Ver computer-graphics-scan-line-polygon-fill-algorithm