# How to use ST_DelaunayTriangles to construct a Voronoi diagram?

(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

• Is gist.github.com/djq/4714788 the sort of thing you are after? Commented Sep 30, 2014 at 19:04
• I think he wants a purely SQL implementation using ST_DelaunayTriangles() Commented Sep 30, 2014 at 20:01
• See this answer to install `ST_DelaunayTriangles` in Linux Debian Stable. Commented Oct 9, 2014 at 0:57
• ! ST_VoronoiPolygons available since PostGIS 2.3 Commented Mar 9, 2018 at 1:12

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

## 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
``````

### 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
``````

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

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.

### 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))))))
``````

• 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? Commented Sep 30, 2014 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 Commented Sep 30, 2014 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"? Commented Oct 1, 2014 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. Commented Oct 1, 2014 at 18:03
• See this answer to install `ST_DelaunayTriangles` in Linux Debian Stable. Commented Oct 9, 2014 at 0:56