Konstruowanie diagramu Voronoi w PostGIS

Próbuję zbudować diagram voronoi z siatki punktów przy użyciu zmodyfikowanego kodu stąd . To jest zapytanie SQL po moich modyfikacjach:

DROP TABLE IF EXISTS example.voronoi;
WITH 
    -- Sample set of points to work with
    Sample AS (SELECT ST_SetSRID(ST_Union(geom), 0) geom FROM example."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, ST_ExteriorRing(ST_ConvexHull(v)))))))).geom, 2180)
INTO example.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)) * 2),ST_Y(x.ct) + ((ST_Y(ST_Centroid(x.edge)) - ST_Y(x.ct)) * 2))
            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

Poniżej - wynik mojego zapytania.

Jak widać, otrzymuję „prawie” prawidłowy schemat voronoi, z wyjątkiem wyróżnionych punktów, które nie mają unikalnego wielokąta. Poniżej jest wytwarzany algorytm QGIS i to, co chciałbym uzyskać z zapytania. Jakieś sugestie, gdzie jest problem z moim kodem?

Chociaż rozwiązuje to bezpośredni problem z zapytaniem o dane, o których mowa, nie jestem z niego zadowolony jako rozwiązanie do ogólnego użytku i wrócę do tej i poprzedniej odpowiedzi, kiedy będę mógł.

Problem polegał na tym, że pierwotne zapytanie wykorzystywało wypukły kadłub na krawędziach voronoi, aby określić zewnętrzną krawędź wielokąta voronoi. Oznaczało to, że niektóre krawędzie voronoi nie docierały na zewnątrz, kiedy powinny. Spojrzałem na użycie funkcji wklęsłego kadłuba, ale tak naprawdę to nie działało tak, jak chciałem.
Jako szybką poprawkę zmieniłem wypukły kadłub, który ma być zbudowany wokół zamkniętego zestawu krawędzi voronoi oraz zderzaka wokół oryginalnych krawędzi. Krawędzie voronoi, które się nie zamykają, są rozciągane na dużą odległość, aby upewnić się, że przecinają się na zewnątrz i służą do budowy wielokątów. Możesz pobawić się parametrami bufora.

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;

Zgodnie z sugestią @ LR1234567 wypróbowania nowej funkcjonalności ST_Voronoi , opracowanej przez @dbaston, oryginalną niesamowitą odpowiedź @MickyT (jak podano w pytaniu OP) i użycie oryginalnych danych można teraz uprościć:

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

co daje wynik, identyczny jak w pytaniu PO.

Jednak cierpi na to ten sam problem, teraz naprawiony w odpowiedzi MickyT, że punkty na wklęsłym kadłubie nie otrzymują otaczającego wielokąta (jak można oczekiwać od algorytmu). Rozwiązałem ten problem za pomocą zapytania z następującą serią kroków.

1. Oblicz wklęsły kadłub punktów wejściowych - punkty na wklęsłym kadłubie mają te, które mają niezwiązane wielokąty na wyjściowym diagramie Voronoi.
2. Znajdź oryginalne punkty na wklęsłym kadłubie (żółte punkty na schemacie 2 poniżej).
3. Buforować wklęsły kadłub (odległość między buforami jest dowolna i być może można go znaleźć bardziej optymalnie na podstawie danych wejściowych?).
4. Znajdź najbliższe punkty w buforze wklęsłego kadłuba najbliższe punktów w kroku 2. Są one pokazane na zielono na poniższym schemacie.
5. Dodaj te punkty do oryginalnego zestawu danych
6. Oblicz diagram Voronoi tego połączonego zestawu danych. Jak widać na trzecim schemacie, punkty na kadłubie mają teraz otaczające wielokąty.

Schemat 2 pokazuje punkty na wklęsłym kadłubie (żółty) i najbliższe punkty buforowania na kadłubie (zielony).

Zapytanie można oczywiście uprościć / skompresować, ale zostawiłem tę formę jako serię CTE, ponieważ łatwiej jest kolejno wykonywać kolejne kroki. To zapytanie działa na oryginalnym zestawie danych w milisekundach (średnio 11 ms na serwerze deweloperskim), podczas gdy odpowiedź MickyT przy użyciu ST_Delauney działa na 4800ms na tym samym serwerze. DBaston twierdzi, że można uzyskać inne zwiększenie prędkości rzędu rzędu od budowania w stosunku do obecnego pnia GEOS, 3,6dev, ze względu na ulepszenia procedur triangulacji.

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;

Schemat 3 pokazujący wszystkie punkty teraz zamknięte w wielokącie

Uwaga: obecnie ST_Voronoi obejmuje budowanie Postgis ze źródła (wersja 2.3 lub trunk) i łączenie z GEOS 3.5 lub nowszym.

Edycja: Właśnie obejrzałem Postgis 2.3, który jest zainstalowany w Amazon Web Services, i wygląda na to, że nazwa funkcji to teraz ST_VoronoiPolygons.

Bez wątpienia to zapytanie / algorytm można ulepszyć. Sugestie mile widziane.

Jeśli masz dostęp do PostGIS 2.3, wypróbuj nową funkcję ST_Voronoi, ostatnio zatwierdzoną:

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

Istnieją wstępnie skompilowane kompilacje dla systemu Windows - http://postgis.net/windows_downloads/