Jak użyć ST_DelaunayTriangles do skonstruowania diagramu Voronoi?

(edycja 2019) ST_VoronoiPolygons dostępne od PostGIS v2.3 !

W PostGIS 2.1+ możemy użyć ST_DelaunayTriangles () do wygenerowania triangulacji Delaunaya , która jest podwójnym wykresem jej diagramu Voronoi , i teoretycznie mają one dokładną i odwracalną konwersję.

Czy istnieje jakiś bezpieczny skrypt w standardzie SQL ze zoptymalizowanym algorytmem dla tej konwersji PostGIS2 Delaunay na Voronoi ?

Inne referencje: 1 , 2

Poniższe zapytanie wydaje się wykonywać rozsądny zestaw wielokątów voronoi, zaczynając od trójkątów Delaunay.

Nie jestem dużym użytkownikiem Postgres, więc prawdopodobnie można go nieco poprawić.

WITH 
    -- Sample set of points to work with
    Sample AS (SELECT ST_GeomFromText('MULTIPOINT (12 5, 5 7, 2 5, 19 6, 19 13, 15 18, 10 20, 4 18, 0 13, 0 6, 4 1, 10 0, 15 1, 19 6)') geom),
    -- 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_Polygonize(ST_Node(ST_LineMerge(ST_Union(v, ST_ExteriorRing(ST_ConvexHull(v))))))
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;

W ten sposób powstaje następujący zestaw wielokątów dla przykładowych punktów zawartych w zapytaniu

Wyjaśnienie zapytania

Krok 1

Utwórz trójkąty Delaunay z geometrii wejściowych

SELECT (gd).Path id, ST_ExteriorRing((gd).Geom) g -- ID and make triangle a linestring
FROM (SELECT (ST_Dump(ST_DelaunayTriangles(geom))) gd FROM Sample) a -- Get Delaunay Triangles

Krok 2

Rozłóż trójkątne węzły i wykonaj krawędzie. Myślę, że powinien istnieć lepszy sposób na uzyskanie krawędzi, ale go nie znalazłem.

SELECT ...
        ST_MakeLine(p1,p2) ,
        ST_MakeLine(p2,p3) ,
        ST_MakeLine(p3,p1)
        ...
FROM    (
    -- Decompose to points
    SELECT id,
        ST_PointN(g,1) p1,
        ST_PointN(g,2) p2,
        ST_PointN(g,3) p3
    FROM    (
        ... Step 1...
        )b
    ) c

Krok 3

Zbuduj opisane okręgi dla każdego trójkąta i znajdź środek ciężkości

SELECT ... Step 2 ...
    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    (
        ... Step 1...
        )b
    ) c

EdgesCTE wyprowadza każdą krawędź i ID (ścieżki) trójkąta, do której należy.

Krok 4

„Outer Join” stół „Edges” do siebie, w którym są równe krawędzie dla różnych trójkątów (krawędzie wewnętrzne).

SELECT  
    ...
    ST_MakeLine(
    x.ct, -- Circumscribed Circle centroid
    CASE 
    WHEN y.id IS NULL THEN
        CASE WHEN ST_Within( -- Don't draw lines back towards the original set
            x.ct,
            (SELECT ST_ConvexHull(geom) FROM sample)) THEN
            -- 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),
                T_Y(x.ct) + ((ST_Y(ST_Centroid(x.edge)) - ST_Y(x.ct)) * 2)
            )
        END
    ELSE 
        y.ct -- Centroid of triangle with common edge
    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)

Tam, gdzie jest wspólna krawędź, narysuj linię między odpowiednimi centroidami

Tam, gdzie krawędź nie jest połączona (zewnętrzna), narysuj linię od środka ciężkości przez środek krawędzi. Zrób to tylko, jeśli środek koła znajduje się w zestawie trójkątów.

Krok 5

Uzyskaj wypukły kadłub dla narysowanych linii jako linię. Połącz i połącz wszystkie linie. Węzeł zestawu linii, abyśmy mieli zestaw topologiczny, który można poligonizować.

SELECT ST_Polygonize(ST_Node(ST_LineMerge(ST_Union(v, ST_ExteriorRing(ST_ConvexHull(v))))))