12

(edit 2019) ST_VoronoiPolygons available since PostGIS v2.3!


With PostGIS 2.1+ we can use ST_DelaunayTriangles() to generate a Delaunay triangulation, that is a dual graph of its Voronoi diagram, and, in theory, they have an exact and reversible conversion.

Does any safe SQL-standard script with an optimized algorithm exist for this PostGIS2 Delaunay-to-Voronoi conversion?


Other refs: 1, 2

21
+50

The following query appears to do a reasonable set of voronoi polygons starting from the Delaunay Triangles.

I'm not a big Postgres user, so it can probably be improved quite a bit.

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;

This produces the following set of polygons for the sample points included in the query enter image description here

Query Explanation

Step 1

Create the Delaunay Triangles from the input geometries

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

Step 2

Decompose the triangle nodes and make edges can be made. I think there should be a better way to get the edges, but I didn't find one.

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

enter image description here

Step 3

Build the circumscribed circles for each triangle and find the centroid

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

enter image description here

The Edges CTE outputs each edge and the id(path) of the triangle it belongs to.

Step 4

'Outer Join' the 'Edges' table to itself where there are equal edges for different triangles (interior edges).

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)

Where there is a common edge draw a line between the respective centroids

enter image description here

Where the edge is not joined (exterior) draw a line from the centroid through the centre of the edge. Only do this if the centroid of the circle is inside the set of triangles.

enter image description here

Step 5

Get the convex hull for the drawn lines as a line. Union up and merge all the lines. Node the line set so that we have a topological set that can be polygonized.

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

enter image description here

  • Good clue, perhaps a solution (!). I need to test, but can't now... Analysing: you use ST_ConvexHull and ST_Centroid instead "perpendicular bisectors" as in the direct algorithm suggested by my ref1/Kenneth Sloa... Why not the direct solution? – Peter Krauss Sep 30 '14 at 20:06
  • I pretty much are doing perpendicular bisectors for the exterior edges, just without all the math:) I'll add an explanation of the steps I took to the answer – MickyT Sep 30 '14 at 20:53
  • Good illustrations and explanation, very didactic!   You posted all that I need (!), but this days I not have Postgis2.1 to test... Can I check here (as comment) some questions that any one can answer by testing?   1) the ST_Polygonize "creates a GeometryCollection containing possible polygons", they are all the Voronoi cells, correct?   2) about performance, you think your centroid-based solution have similar CPU-time than "all the math of perpendicular bisectors calculation"? – Peter Krauss Oct 1 '14 at 11:21
  • @PeterKrauss 1) ST_polygonize does create the voronoi cells from the line work. The trick with it is to make sure all the line work is split at the nodes. 2) I don't think there would be a lot of difference between calculating the bisection and using ST_Centroid on the line. But it would need to be tested. – MickyT Oct 1 '14 at 18:03
  • See this answer to install ST_DelaunayTriangles in Linux Debian Stable. – Peter Krauss Oct 9 '14 at 0:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.