Como usar ST_DelaunayTriangles para construir um diagrama de Voronoi?

12

(editar 2019) ST_VoronoiPolygons disponíveis desde o PostGIS v2.3 !


Com o PostGIS 2.1+, podemos usar ST_DelaunayTriangles () para gerar uma triangulação de Delaunay , que é um gráfico duplo do diagrama de Voronoi e, em teoria, eles têm uma conversão exata e reversível.

Existe algum script padrão do SQL seguro com um algoritmo otimizado para esta conversão PostGIS2 Delaunay em Voronoi ?


Outras refs: 1 , 2

Peter Krauss
fonte
Gist.github.com/djq/4714788 é o tipo de coisa que você procura ?
MickyT
Acho que ele quer uma implementação puramente SQL usando ST_DelaunayTriangles ()
raphael
Veja esta resposta para instalar ST_DelaunayTrianglesno Linux Debian Stable .
Peter Krauss
! ST_VoronoiPolygons disponíveis desde o PostGIS 2.3
Peter Krauss

Respostas:

22

A consulta a seguir parece fazer um conjunto razoável de polígonos de voronoi a partir dos triângulos de Delaunay.

Como não sou um grande usuário do Postgres, é provável que possa ser melhorado bastante.

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;

Isso produz o seguinte conjunto de polígonos para os pontos de amostra incluídos na consulta insira a descrição da imagem aqui

Explicação da consulta

Passo 1

Crie os triângulos de Delaunay a partir das geometrias de entrada

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

Passo 2

Decomponha os nós do triângulo e faça arestas. Acho que deveria haver uma maneira melhor de obter as bordas, mas não encontrei uma.

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

insira a descrição da imagem aqui

etapa 3

Construa os círculos circunscritos para cada triângulo e encontre o centróide

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

insira a descrição da imagem aqui

O EdgesCTE gera cada aresta e o id (caminho) do triângulo ao qual pertence.

Passo 4

'Junção externa' a tabela 'Bordas', onde existem arestas iguais para triângulos diferentes (arestas internas).

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)

Onde existe uma aresta comum, desenhe uma linha entre os respectivos centróides

insira a descrição da imagem aqui

Onde a aresta não estiver unida (exterior), desenhe uma linha do centróide através do centro da aresta. Faça isso apenas se o centróide do círculo estiver dentro do conjunto de triângulos.

insira a descrição da imagem aqui

Etapa 5

Obtenha o casco convexo para as linhas desenhadas como uma linha. Unir e mesclar todas as linhas. Nó o conjunto de linhas para que tenhamos um conjunto topológico que possa ser poligonalizado.

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

insira a descrição da imagem aqui

MickyT
fonte
Boa pista, talvez uma solução (!). Preciso testar, mas não posso agora ... Analisar: você usa ST_ConvexHulle, em ST_Centroidvez disso, "bissetores perpendiculares", como no algoritmo direto sugerido pelo meu ref1 / Kenneth Sloa ... Por que não a solução direta?
Peter Krauss
Eu praticamente estão fazendo mediatrizes para as bordas exteriores, só que sem toda a matemática :) Vou acrescentar uma explicação sobre os passos que tomou para a resposta
MickyT
Boas ilustrações e explicações, muito didáticas!   Você postou tudo o que eu preciso (!), Mas hoje em dia não tenho o Postgis2.1 para testar ... Posso verificar aqui (como comentário) algumas perguntas que qualquer um pode responder testando?   1) o ST_Polygonize "cria uma GeometryCollection contendo possíveis polígonos", são todas as células de Voronoi, correto?   2) sobre desempenho, você acha que sua solução baseada em centróide tem um tempo de CPU semelhante a "toda a matemática do cálculo de bissetores perpendiculares"?
Peter Krauss
@PeterKrauss 1) ST_polygonize cria as células voronoi a partir da linha de trabalho. O truque é garantir que todo o trabalho da linha seja dividido nos nós. 2) Acho que não haveria muita diferença entre calcular a bissecção e usar o ST_Centroid na linha. Mas isso precisaria ser testado.
MickyT
Veja esta resposta para instalar ST_DelaunayTrianglesno Linux Debian Stable .
Peter Krauss