I'm trying to construct voronoi diagram from grid of points using modified code from here. This is SQL query after my modifications:
DROP TABLE IF EXISTS example.voronoi;
WITH
-- Sample set of points to work with
Sample AS (SELECT ST_SetSRID(ST_Union(geom), 0) geom FROM example."MeshPoints2d"),
-- 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_SetSRID((ST_Dump(ST_Polygonize(ST_Node(ST_LineMerge(ST_Union(v, ST_ExteriorRing(ST_ConvexHull(v)))))))).geom, 2180)
INTO example.voronoi
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
As you can see I get "almost" correct voronoi diagram, except highlighted points which don't have unique polygon. Below is what QGIS algorithm produce and what I would like to obtain from query. Any sugestions where is the problem with my code?
Best Answer
Following a suggestion by @LR1234567 to try out the new ST_Voronoi functionality that has been developed by @dbaston, @MickyT 's original amazing answer (as stated in OP's question) and using the original data can now be simplified to:
which results in this output, identical to the OP's question.
However, this suffers from the same issue, now fixed in MickyT's answer, that points on the concave hull do not get an enclosing polygon (as would be expected from the algorithm). I fixed this issue with a query with the following series of steps.
Diagram 2 showing points on concave hull (yellow) and closest points to buffer on hull (green).
The query could obviously be simplified/compressed, but I left it is this form as a series of CTEs, as it is easier to follow the steps sequentially that way. This query runs on the original data set in milliseconds (11ms average on a dev server) whereas MickyT's answer using ST_Delauney runs in 4800ms on the same server. DBaston claims that another order of magnitude speed enhancement can be gained from building against the current trunk of GEOS, 3.6dev, due to enhancements in triangulation routines.
Diagram 3 showing all points now enclosed in a polygon
Note: Currently ST_Voronoi involves building Postgis from source (version 2.3, or trunk) and linking against GEOS 3.5 or above.
Edit: I've just looked at Postgis 2.3 as it is installed on Amazon Web Services, and it seems the function name is now ST_VoronoiPolygons.
No doubt this query/algorithm could be improved. Suggestions welcome.