Suite

Comment utiliser ST_DelaunayTriangles pour construire un diagramme de Voronoi ?

Comment utiliser ST_DelaunayTriangles pour construire un diagramme de Voronoi ?


(modifier 2019) ST_VoronoiPolygones disponibles depuis PostGIS v2.3!


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

Est-ce que tout coffre-fort Script standard SQL avec un algorithme optimisé existent pour cette conversion PostGIS2 Delaunay à Voronoi?


Autres références : 1, 2


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

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

WITH -- Exemple d'ensemble de points à utiliser avec 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), -- Construire des arêtes et circonscrire des points pour générer un centre de gravité 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(-- Fait de cette façon en raison de problèmes avec 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 ( -- Décomposer en 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 Tr iangles )b ) c ) SELECT ST_Polygonize(ST_Node(ST_LineMerge(ST_Union(v, ST_ExteriorRing(ST_ConvexHull(v)))))) FROM ( SELECT -- Créer des arêtes voronoi et réduire en une chaîne multiligne 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 -- Ne pas tracer de lignes vers l'ensemble d'origine -- Projeter la ligne deux fois plus loin que convex coque 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 -- Auto-jointure basée sur les arêtes Edges y ON x.id <> y.id AND ST_Equals( x.bord,y.bord) ) 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 et faire du triangle une chaîne de lignes FROM (SELECT (ST_Dump(ST_DelaunayTriangles(geom))) gd FROM Sample) a -- Get Delaunay Triangles

Étape 2

Décomposez les nœuds du triangle et créez 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 ( -- Décompose en points SELECT id, ST_PointN(g,1) p1, ST_PointN(g,2) p2, ST_PointN (g,3) p3 DE (… Étape 1… )b ) c

Étape 3

Construisez les cercles circonscrits pour chaque triangle et trouvez le centre de gravité

SELECT… Étape 2… ST_Centroid(ST_ConvexHull(ST_Union(-- Fait de cette façon en raison de problèmes que j'ai eus avec 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 ( -- Décomposer en points SELECT id, ST_PointN(g,1) p1, ST_PointN(g,2) p2, ST_PointN(g,3) p3 FROM (… Étape 1… )b ) c

leBordsCTE génère chaque bord et l'identifiant (chemin) du triangle auquel il appartient.

Étape 4

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

SELECT… ST_MakeLine( x.ct, -- Centroïde du cercle circonscrit CASE WHEN y.id IS NULL THEN CASE WHEN ST_Within( -- Ne dessine pas de lignes vers l'ensemble d'origine x.ct, (SELECT ST_ConvexHull(geom) FROM sample) ) ALORS -- Ligne de projet sur deux fois la distance de l'enveloppe convexe 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 -- Centroïde du triangle avec arête commune END ))) v FROM Edges x LEFT OUTER JOIN - - Auto-jointure basée sur les bords Bords y ON x.id <> y.id ET 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 du centre de gravité au centre du bord. Ne le faites que si le centre de gravité du cercle est à l'intérieur de l'ensemble des triangles.

Étape 5

Obtenez l'enveloppe convexe des lignes dessinées sous forme de ligne. Union vers le haut et fusionner toutes les lignes. Noeud l'ensemble de lignes de manière à avoir un ensemble topologique pouvant être polygonisé.

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


Voir la vidéo: Delaunay Triangulation and Voronoi Diagram Demo