Wie benutze ich ST_DelaunayTriangles, um ein Voronoi-Diagramm zu erstellen?

(edit 2019) ST_VoronoiPolygons verfügbar seit PostGIS v2.3 !

Mit PostGIS 2.1+ können wir ST_DelaunayTriangles () verwenden , um eine Delaunay-Triangulation zu generieren , die ein duales Diagramm des Voronoi-Diagramms ist und theoretisch eine genaue und reversible Konvertierung aufweist.

Gibt es ein sicheres SQL-Standard-Skript mit einem optimierten Algorithmus für diese PostGIS2-Konvertierung von Delaunay nach Voronoi ?

Andere Referenzen: 1 , 2

Die folgende Abfrage scheint eine sinnvolle Menge von Voronoi-Polygonen zu erstellen, die von den Delaunay-Dreiecken ausgehen.

Ich bin kein großer Postgres-Benutzer, daher kann es wahrscheinlich einiges verbessert werden.

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;

Dadurch wird die folgende Gruppe von Polygonen für die in der Abfrage enthaltenen Beispielpunkte erstellt

Abfrage Erläuterung

Schritt 1

Erstellen Sie die Delaunay-Dreiecke aus den Eingabegeometrien

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

Schritt 2

Zerlegen Sie die Dreiecksknoten, und erstellen Sie Kanten. Ich denke, es sollte einen besseren Weg geben, um die Kanten zu bekommen, aber ich habe keinen gefunden.

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

Schritt 3

Bilde für jedes Dreieck die umschriebenen Kreise und finde den Schwerpunkt

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

Der EdgesCTE gibt jede Kante und die ID (Pfad) des Dreiecks aus, zu dem er gehört.

Schritt 4

"Außen" Verbinden Sie die Tabelle "Kanten" mit sich selbst, wobei es gleiche Kanten für verschiedene Dreiecke (Innenkanten) gibt.

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)

Wo es eine gemeinsame Kante gibt, ziehen Sie eine Linie zwischen den jeweiligen Schwerpunkten

Wenn die Kante nicht verbunden ist (außen), ziehen Sie eine Linie vom Schwerpunkt durch die Mitte der Kante. Tun Sie dies nur, wenn der Schwerpunkt des Kreises innerhalb der Dreiecke liegt.

Schritt 5

Holen Sie sich die konvexe Hülle für die gezeichneten Linien als Linie. Vereinige alle Linien und füge sie zusammen. Knoten Sie den Liniensatz, damit wir einen topologischen Satz haben, der polygonisiert werden kann.

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