Solución:
Siguiendo una sugerencia de @ LR1234567 para probar la nueva funcionalidad ST_Voronoi que ha sido desarrollada por @dbaston, la sorprendente respuesta original de @MickyT (como se indica en la pregunta de OP) y el uso de los datos originales ahora se puede simplificar a:
WITH voronoi (vor) AS
(SELECT ST_Dump(ST_Voronoi(ST_Collect(geom))) FROM meshpoints)
SELECT (vor).path, (vor).geom FROM voronoi;
lo que da como resultado esta salida, idéntica a la pregunta del OP.
Sin embargo, esto sufre el mismo problema, ahora solucionado en la respuesta de MickyT, que los puntos en el casco cóncavo no obtienen un polígono circundante (como se esperaría del algoritmo). Solucioné este problema con una consulta con la siguiente serie de pasos.
- Calcule el casco cóncavo de los puntos de entrada; los puntos del casco cóncavo son aquellos que tienen polígonos ilimitados en el diagrama de Voronoi de salida.
- Encuentre los puntos originales en el casco cóncavo (puntos amarillos en el diagrama 2 a continuación).
- Proteja el casco cóncavo (¿la distancia del búfer es arbitraria y quizás podría encontrarse de manera más óptima a partir de los datos de entrada?).
- Encuentre los puntos más cercanos en la zona de influencia del casco cóncavo más cercanos a los puntos en el paso 2. Estos se muestran en verde en el diagrama siguiente.
- Agregue estos puntos al conjunto de datos original
- Calcule el diagrama de Voronoi de este conjunto de datos combinados. Como se puede ver en el tercer diagrama, los puntos del casco ahora tienen polígonos circundantes.
Diagrama 2 que muestra los puntos en el casco cóncavo (amarillo) y los puntos más cercanos al amortiguador en el casco (verde).
Obviamente, la consulta podría simplificarse / comprimirse, pero dejé este formulario como una serie de CTE, ya que es más fácil seguir los pasos secuencialmente de esa manera. Esta consulta se ejecuta en el conjunto de datos original en milisegundos (promedio de 11 ms en un servidor de desarrollo), mientras que la respuesta de MickyT usando ST_Delauney se ejecuta en 4800 ms en el mismo servidor. DBaston afirma que se puede obtener otro orden de magnitud de mejora de la velocidad construyendo contra el tronco actual de GEOS, 3.6dev, debido a las mejoras en las rutinas de triangulación.
WITH
conv_hull(geom) AS
(SELECT ST_Concavehull(ST_Union(geom), 1) FROM meshpoints),
edge_points(points) AS
(SELECT mp.geom FROM meshpoints mp, conv_hull ch
WHERE ST_Touches(ch.geom, mp.geom)),
buffered_points(geom) AS
(SELECT ST_Buffer(geom, 100) as geom FROM conv_hull),
closest_points(points) AS
(SELECT
ST_Closestpoint(
ST_Exteriorring(bp.geom), ep.points) as points,
ep.points as epoints
FROM buffered_points bp, edge_points ep),
combined_points(points) AS
(SELECT points FROM closest_points
UNION SELECT geom FROM meshpoints),
voronoi (vor) AS
(SELECT
ST_Dump(
ST_Voronoi(
ST_Collect(points))) as geom
FROM combined_points)
SELECT
(vor).path[1] as id,
(vor).geom
FROM voronoi;
Diagrama 3 que muestra todos los puntos ahora encerrados en un polígono
Nota: Actualmente, ST_Voronoi implica compilar Postgis desde la fuente (versión 2.3 o troncal) y vincularlo con GEOS 3.5 o superior.
Editar: Acabo de ver Postgis 2.3 ya que está instalado en Amazon Web Services, y parece que el nombre de la función ahora es ST_VoronoiPolygons.
Sin duda, esta consulta / algoritmo podría mejorarse. Sugerencias bienvenidas.
Si bien esto soluciona el problema inmediato con la consulta de los datos en cuestión, no estoy satisfecho con ella como una solución para el uso general y volveré a revisar esta y la respuesta anterior cuando pueda.
El problema era que la consulta original utilizaba un casco convexo en los bordes de voronoi para determinar el borde exterior del polígono de voronoi. Esto significaba que algunos de los bordes de los voronoi no llegaban al exterior cuando deberían haberlo hecho. Miré el uso de la funcionalidad del casco cóncavo, pero realmente no funcionó como quería.
Como solución rápida, he cambiado el casco convexo para que se construya alrededor del conjunto cerrado de bordes voronoi más un amortiguador alrededor de los bordes originales. Los bordes de voronoi que no se cierran se extienden una gran distancia para intentar asegurarse de que se crucen con el exterior y se utilicen para construir polígonos. Es posible que desee jugar con los parámetros del búfer.
WITH
-- Sample set of points to work with
Sample AS (SELECT ST_SetSRID(ST_Union(geom), 0) geom FROM 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, (SELECT ST_ExteriorRing(ST_ConvexHull(ST_Union(ST_Union(ST_Buffer(edge,20),ct)))) FROM Edges))))))).geom, 2180) geom
INTO 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)) * 200),ST_Y(x.ct) + ((ST_Y(ST_Centroid(x.edge)) - ST_Y(x.ct)) * 200))
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;
Si tiene acceso a PostGIS 2.3, pruebe la nueva función ST_Voronoi, recientemente comprometida:
http://postgis.net/docs/manual-dev/ST_Voronoi.html
Hay compilaciones precompiladas para Windows: http://postgis.net/windows_downloads/