Comment utiliser ST_DelaunayTriangles pour construire un diagramme de Voronoi?

(modifier 2019) ST_VoronoiPolygons disponible depuis PostGIS v2.3 !

Avec PostGIS 2.1+, nous pouvons utiliser ST_DelaunayTriangles () pour générer une triangulation Delaunay , c'est-à-dire un double graphique de son diagramme de Voronoi , et, en théorie, ils ont une conversion exacte et réversible.

Existe- t -il un script standard SQL sécurisé avec un algorithme optimisé pour cette conversion PostGIS2 Delaunay-à-Voronoi ?

Autres références: 1 , 2

La requête suivante semble faire un ensemble raisonnable de polygones voronoi à partir des triangles de Delaunay.

Je ne suis pas un grand utilisateur de Postgres, donc cela peut probablement être amélioré un peu.

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;

Cela produit l'ensemble de polygones suivant pour les points d'échantillonnage inclus dans la requête

Explication de la requête

Étape 1

Créer les triangles de Delaunay à partir des géométries d'entrée

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

Étape 2

Décomposer les nœuds triangulaires et créer des arêtes. Je pense qu'il devrait y avoir une meilleure façon d'obtenir les bords, mais je n'en ai pas trouvé.

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

Étape 3

Construisez les cercles circonscrits pour chaque triangle et trouvez le centroïde

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

Le EdgesCTE sort chaque front et l'id (chemin) du triangle auquel il appartient.

Étape 4

'Outer Join' la table 'Edges' à elle-même où il y a des bords égaux pour différents triangles (bords intérieurs).

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)

Là où il y a un bord commun, tracez une ligne entre les centroïdes respectifs

Lorsque le bord n'est pas joint (extérieur), tracez une ligne à partir du centre de gravité en passant par le centre du bord. Ne faites cela que si le centre de gravité du cercle se trouve à l'intérieur de l'ensemble de triangles.

Étape 5

Obtenez la coque convexe pour les lignes tracées sous forme de ligne. Rassemblez et fusionnez toutes les lignes. Noeud l'ensemble de lignes de sorte que nous avons un ensemble topologique qui peut être polygonisé.

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