Как использовать ST_DelaunayTriangles для построения диаграммы Вороного?

13

(редактировать 2019) ST_VoronoiPolygons доступны с PostGIS v2.3 !


С PostGIS 2.1+ мы можем использовать ST_DelaunayTriangles () для генерации триангуляции Делоне , которая является двойственным графом ее диаграммы Вороного , и, теоретически, они имеют точное и обратимое преобразование.

Существует ли какой-либо безопасный сценарий SQL-стандарта с оптимизированным алгоритмом для этого преобразования PostGIS2 Делоне-в Вороной ?


Другие ссылки: 1 , 2

Питер Краусс
источник
Является gist.github.com/djq/4714788 вид вещи вы после?
MickyT
Я думаю, что он хочет чисто SQL-реализацию с использованием ST_DelaunayTriangles ()
Raphael
Смотрите этот ответ для установки ST_DelaunayTrianglesв Linux Debian Stable .
Питер Краусс
! ST_VoronoiPolygons доступны с PostGIS 2.3
Питер Краусс

Ответы:

23

Следующий запрос, по-видимому, делает разумный набор полигонов вороной, начиная с треугольников Делоне.

Я не большой пользователь Postgres, так что, возможно, его можно немного улучшить.

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;

Это создает следующий набор полигонов для точек выборки, включенных в запрос введите описание изображения здесь

Объяснение запроса

Шаг 1

Создайте треугольники Делоне из входных геометрий

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

Шаг 2

Разложить узлы треугольника и сделать ребра можно. Я думаю, что должен быть лучший способ получить края, но я не нашел один.

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

введите описание изображения здесь

Шаг 3

Постройте описанные круги для каждого треугольника и найдите центроид

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

введите описание изображения здесь

EdgesКТР выводит каждое ребро и идентификатор (путь) треугольника он принадлежит.

Шаг 4

«Внешнее» Соедините таблицу «Края» с собой, где имеются равные ребра для разных треугольников (внутренних ребер).

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)

Там, где есть общий край, проведите линию между соответствующими центроидами.

введите описание изображения здесь

Там, где край не соединен (снаружи), проведите линию от центра тяжести через центр края. Делайте это только в том случае, если центр тяжести круга находится внутри набора треугольников.

введите описание изображения здесь

Шаг 5

Получите выпуклый корпус для нарисованных линий в виде линии. Объедините и объедините все линии. Узел линии установлен так, чтобы у нас был топологический набор, который можно полигонизировать.

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

введите описание изображения здесь

MickyT
источник
Хорошая подсказка, возможно решение (!). Мне нужно проверить, но не могу сейчас ... Анализ: вы используете ST_ConvexHullи ST_Centroidвместо этого "перпендикулярные биссектрисы", как в прямом алгоритме, предложенном моим ref1 / Kenneth Sloa ... Почему не прямое решение?
Питер Краусс
Я в основном делаю перпендикулярные биссектрисы для внешних краев, просто без всякой математики :) Я добавлю объяснение шагов, которые я предпринял к ответу
MickyT
Хорошие иллюстрации и объяснения, очень дидактические!   Вы отправили все, что мне нужно (!), Но в эти дни у меня нет Postgis2.1 для тестирования ... Могу ли я проверить здесь (в качестве комментария) некоторые вопросы, на которые любой может ответить тестированием?   1) ST_Polygonize «создает коллекцию GeometryCollection, содержащую возможные полигоны», все они являются ячейками Вороного, верно?   2) Что касается производительности, вы думаете, что ваше решение на базе центроида имеет схожее время процессора, чем «все вычисления перпендикулярных биссектрис»?
Питер Краусс
@PeterKrauss 1) ST_polygonize создает клетки вороной из линейной работы. Хитрость заключается в том, чтобы убедиться, что вся работа линии разбита на узлы. 2) Я не думаю, что будет большая разница между вычислением деления пополам и использованием ST_Centroid на линии. Но это должно быть проверено.
MickyT
Смотрите этот ответ для установки ST_DelaunayTrianglesв Linux Debian Stable .
Питер Краусс