Bagaimana cara menggunakan ST_DelaunayTriangles untuk membuat diagram Voronoi?

13

(sunting 2019) ST_VoronoiPolygons tersedia sejak PostGIS v2.3 !


Dengan PostGIS 2.1+ kita dapat menggunakan ST_DelaunayTriangles () untuk menghasilkan triangulasi Delaunay , yaitu grafik ganda dari diagram Voronoi -nya , dan, secara teori, mereka memiliki konversi yang tepat dan dapat dibalik.

Apakah ada skrip SQL- aman yang aman dengan algoritma yang dioptimalkan ada untuk konversi PostGIS2 Delaunay-to-Voronoi ini ?


Referensi lainnya: 1 , 2

Peter Krauss
sumber
Apakah gist.github.com/djq/4714788 adalah hal yang Anda cari?
MickyT
Saya pikir dia ingin implementasi SQL murni menggunakan ST_DelaunayTriangles ()
raphael
Lihat jawaban ini untuk menginstal ST_DelaunayTrianglesdi Linux Debian Stable .
Peter Krauss
! ST_VoronoiPolygons tersedia sejak PostGIS 2.3
Peter Krauss

Jawaban:

23

Kueri berikut tampaknya melakukan serangkaian voronoi poligon yang wajar mulai dari Delaunay Segitiga.

Saya bukan pengguna Postgres besar, jadi mungkin bisa ditingkatkan sedikit.

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;

Ini menghasilkan set poligon berikut untuk titik sampel yang termasuk dalam kueri masukkan deskripsi gambar di sini

Penjelasan Pertanyaan

Langkah 1

Buat Segitiga Delaunay dari geometri input

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

Langkah 2

Membusuk node segitiga dan membuat tepi dapat dibuat. Saya pikir seharusnya ada cara yang lebih baik untuk mendapatkan tepian, tetapi saya tidak menemukannya.

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

masukkan deskripsi gambar di sini

Langkah 3

Bangun lingkaran yang dibatasi untuk setiap segitiga dan temukan 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

masukkan deskripsi gambar di sini

The EdgesCTE output setiap tepi dan id (path) dari segitiga itu milik.

Langkah 4

'Outer Join' the 'Edges' table untuk dirinya sendiri di mana ada tepi yang sama untuk segitiga yang berbeda (tepi interior).

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)

Di mana ada tepi umum, buat garis di antara centroid masing-masing

masukkan deskripsi gambar di sini

Di mana ujung tidak bergabung (eksterior) gambar garis dari pusat massa melalui pusat tepi. Hanya lakukan ini jika pusat massa dari lingkaran berada di dalam himpunan segitiga.

masukkan deskripsi gambar di sini

Langkah 5

Dapatkan cembung cembung untuk garis yang ditarik sebagai garis. Menyatukan dan menggabungkan semua lini. Node set garis sehingga kita memiliki set topologi yang dapat dipoligonisasi.

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

masukkan deskripsi gambar di sini

MickyT
sumber
Petunjuk bagus, mungkin solusi (!). Saya perlu menguji, tetapi tidak bisa sekarang ... Menganalisis: Anda menggunakan ST_ConvexHulldan ST_Centroidsebagai gantinya "garis-garis tegak lurus" seperti dalam algoritma langsung yang disarankan oleh ref1 / Kenneth Sloa saya ... Mengapa bukan solusi langsung?
Peter Krauss
Saya cukup banyak melakukan garis-garis tegak lurus untuk tepi luar, hanya tanpa semua matematika :) Saya akan menambahkan penjelasan tentang langkah-langkah yang saya ambil untuk jawabannya
MickyT
Ilustrasi dan penjelasan yang baik, sangat didaktik!   Anda memposting semua yang saya butuhkan (!), Tetapi hari ini saya tidak memiliki Postgis2.1 untuk diuji ... Dapatkah saya memeriksa di sini (sebagai komentar) beberapa pertanyaan yang dapat dijawab oleh seseorang dengan menguji?   1) ST_Polygonize "membuat GeometryCollection berisi kemungkinan poligon", semuanya adalah sel Voronoi, benar?   2) tentang kinerja, apakah menurut Anda solusi berbasis centroid Anda memiliki waktu CPU yang sama dengan "semua perhitungan matematika garis lurus"?
Peter Krauss
@PeterKrauss 1) ST_polygonize tidak membuat sel voronoi dari pekerjaan garis. Trik dengan itu adalah untuk memastikan semua pekerjaan garis terpecah di node. 2) Saya tidak berpikir akan ada banyak perbedaan antara menghitung pembelahan dua dan menggunakan ST_Centroid di telepon. Tapi itu perlu diuji.
MickyT
Lihat jawaban ini untuk menginstal ST_DelaunayTrianglesdi Linux Debian Stable .
Peter Krauss