Saltar al contenido

Construyendo el diagrama de Voronoi en PostGIS

Solución:

Siguiendo una sugerencia de @ LR1234567 para probar la nueva funcionalidad ST_Voronoi que ha sido desarrollada por @dbaston, la sorprendente respuesta original de @MickyT (como se indica en la pregunta de OP) y el uso de los datos originales ahora se puede simplificar a:

WITH voronoi (vor) AS 
     (SELECT ST_Dump(ST_Voronoi(ST_Collect(geom))) FROM meshpoints)
SELECT (vor).path, (vor).geom FROM voronoi;

lo que da como resultado esta salida, idéntica a la pregunta del OP.

ingrese la descripción de la imagen aquí

Sin embargo, esto sufre el mismo problema, ahora solucionado en la respuesta de MickyT, que los puntos en el casco cóncavo no obtienen un polígono circundante (como se esperaría del algoritmo). Solucioné este problema con una consulta con la siguiente serie de pasos.

  1. Calcule el casco cóncavo de los puntos de entrada; los puntos del casco cóncavo son aquellos que tienen polígonos ilimitados en el diagrama de Voronoi de salida.
  2. Encuentre los puntos originales en el casco cóncavo (puntos amarillos en el diagrama 2 a continuación).
  3. Proteja el casco cóncavo (¿la distancia del búfer es arbitraria y quizás podría encontrarse de manera más óptima a partir de los datos de entrada?).
  4. Encuentre los puntos más cercanos en la zona de influencia del casco cóncavo más cercanos a los puntos en el paso 2. Estos se muestran en verde en el diagrama siguiente.
  5. Agregue estos puntos al conjunto de datos original
  6. Calcule el diagrama de Voronoi de este conjunto de datos combinados. Como se puede ver en el tercer diagrama, los puntos del casco ahora tienen polígonos circundantes.

Diagrama 2 que muestra los puntos en el casco cóncavo (amarillo) y los puntos más cercanos al amortiguador en el casco (verde).
Diagrama 2.

Obviamente, la consulta podría simplificarse / comprimirse, pero dejé este formulario como una serie de CTE, ya que es más fácil seguir los pasos secuencialmente de esa manera. Esta consulta se ejecuta en el conjunto de datos original en milisegundos (promedio de 11 ms en un servidor de desarrollo), mientras que la respuesta de MickyT usando ST_Delauney se ejecuta en 4800 ms en el mismo servidor. DBaston afirma que se puede obtener otro orden de magnitud de mejora de la velocidad construyendo contra el tronco actual de GEOS, 3.6dev, debido a las mejoras en las rutinas de triangulación.

WITH 
  conv_hull(geom) AS 
        (SELECT ST_Concavehull(ST_Union(geom), 1) FROM meshpoints), 
  edge_points(points) AS 
        (SELECT mp.geom FROM meshpoints mp, conv_hull ch 
        WHERE ST_Touches(ch.geom, mp.geom)), 
  buffered_points(geom) AS
        (SELECT ST_Buffer(geom, 100) as geom FROM conv_hull),
  closest_points(points) AS
        (SELECT 
              ST_Closestpoint(
                   ST_Exteriorring(bp.geom), ep.points) as points,
             ep.points as epoints 
         FROM buffered_points bp, edge_points ep),
  combined_points(points) AS
        (SELECT points FROM closest_points 
        UNION SELECT geom FROM meshpoints),
  voronoi (vor) AS 
       (SELECT 
            ST_Dump(
                  ST_Voronoi(
                    ST_Collect(points))) as geom 
        FROM combined_points)
 SELECT 
     (vor).path[1] as id, 
     (vor).geom 
 FROM voronoi;

Diagrama 3 que muestra todos los puntos ahora encerrados en un polígono
diagrama 3

Nota: Actualmente, ST_Voronoi implica compilar Postgis desde la fuente (versión 2.3 o troncal) y vincularlo con GEOS 3.5 o superior.

Editar: Acabo de ver Postgis 2.3 ya que está instalado en Amazon Web Services, y parece que el nombre de la función ahora es ST_VoronoiPolygons.

Sin duda, esta consulta / algoritmo podría mejorarse. Sugerencias bienvenidas.

Si bien esto soluciona el problema inmediato con la consulta de los datos en cuestión, no estoy satisfecho con ella como una solución para el uso general y volveré a revisar esta y la respuesta anterior cuando pueda.

El problema era que la consulta original utilizaba un casco convexo en los bordes de voronoi para determinar el borde exterior del polígono de voronoi. Esto significaba que algunos de los bordes de los voronoi no llegaban al exterior cuando deberían haberlo hecho. Miré el uso de la funcionalidad del casco cóncavo, pero realmente no funcionó como quería.
Como solución rápida, he cambiado el casco convexo para que se construya alrededor del conjunto cerrado de bordes voronoi más un amortiguador alrededor de los bordes originales. Los bordes de voronoi que no se cierran se extienden una gran distancia para intentar asegurarse de que se crucen con el exterior y se utilicen para construir polígonos. Es posible que desee jugar con los parámetros del búfer.

WITH 
    -- Sample set of points to work with
    Sample AS (SELECT ST_SetSRID(ST_Union(geom), 0) geom FROM MeshPoints2d),
    -- Build edges and circumscribe points to generate a centroid
    Edges AS (
    SELECT id,
        UNNEST(ARRAY['e1','e2','e3']) EdgeName,
        UNNEST(ARRAY[
            ST_MakeLine(p1,p2) ,
            ST_MakeLine(p2,p3) ,
            ST_MakeLine(p3,p1)]) Edge,
        ST_Centroid(ST_ConvexHull(ST_Union(-- Done this way due to issues I had with LineToCurve
            ST_CurveToLine(REPLACE(ST_AsText(ST_LineMerge(ST_Union(ST_MakeLine(p1,p2),ST_MakeLine(p2,p3)))),'LINE','CIRCULAR'),15),
            ST_CurveToLine(REPLACE(ST_AsText(ST_LineMerge(ST_Union(ST_MakeLine(p2,p3),ST_MakeLine(p3,p1)))),'LINE','CIRCULAR'),15)
        ))) ct      
    FROM    (
        -- Decompose to points
        SELECT id,
            ST_PointN(g,1) p1,
            ST_PointN(g,2) p2,
            ST_PointN(g,3) p3
        FROM    (
            SELECT (gd).Path id, ST_ExteriorRing((gd).geom) g -- ID andmake triangle a linestring
            FROM (SELECT (ST_Dump(ST_DelaunayTriangles(geom))) gd FROM Sample) a -- Get Delaunay Triangles
            )b
        ) c
    )
SELECT ST_SetSRID((ST_Dump(ST_Polygonize(ST_Node(ST_LineMerge(ST_Union(v, (SELECT ST_ExteriorRing(ST_ConvexHull(ST_Union(ST_Union(ST_Buffer(edge,20),ct)))) FROM Edges))))))).geom, 2180) geom
INTO voronoi
FROM (
    SELECT  -- Create voronoi edges and reduce to a multilinestring
        ST_LineMerge(ST_Union(ST_MakeLine(
        x.ct,
        CASE 
        WHEN y.id IS NULL THEN
            CASE WHEN ST_Within(
                x.ct,
                (SELECT ST_ConvexHull(geom) FROM sample)) THEN -- Don't draw lines back towards the original set
                -- Project line out twice the distance from convex hull
                ST_MakePoint(ST_X(x.ct) + ((ST_X(ST_Centroid(x.edge)) - ST_X(x.ct)) * 200),ST_Y(x.ct) + ((ST_Y(ST_Centroid(x.edge)) - ST_Y(x.ct)) * 200))
            END
        ELSE 
            y.ct
        END
        ))) v
    FROM    Edges x 
        LEFT OUTER JOIN -- Self Join based on edges
        Edges y ON x.id <> y.id AND ST_Equals(x.edge,y.edge)
    ) z;

Si tiene acceso a PostGIS 2.3, pruebe la nueva función ST_Voronoi, recientemente comprometida:

http://postgis.net/docs/manual-dev/ST_Voronoi.html

Hay compilaciones precompiladas para Windows: http://postgis.net/windows_downloads/

¡Haz clic para puntuar esta entrada!
(Votos: 0 Promedio: 0)



Utiliza Nuestro Buscador

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *