Membangun diagram Voronoi di PostGIS

12

Saya mencoba membuat diagram voronoi dari kisi poin menggunakan kode yang dimodifikasi dari sini . Ini adalah permintaan SQL setelah modifikasi saya:

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

Di bawah - hasil permintaan saya. masukkan deskripsi gambar di sini

Seperti yang Anda lihat, saya mendapatkan "hampir" diagram voronoi yang benar, kecuali titik-titik yang disorot yang tidak memiliki poligon unik. Di bawah ini adalah apa yang dihasilkan algoritma QGIS dan apa yang ingin saya peroleh dari kueri. Adakah saran di mana masalah dengan kode saya?

masukkan deskripsi gambar di sini

DamnBack
sumber
Mungkin Anda juga dapat membandingkan hasil fungsi SpatiaLite "VoronojDiagram" gaia-gis.it/gaia-sins/spatialite-sql-latest.html dan melihat kode sumber di gaia-gis.it/fossil/libspatialite/ indeks .
user30184
Pertanyaan yang bagus, saya telah melihat pertanyaan yang sama dengan referensi Anda, dengan maksud untuk mempercepatnya, tetapi tetap kehabisan waktu. Saya tidak menyadari masalah ini dengan poin luar.
John Powell
5
Untuk apa layaknya kita memiliki ST_Voronoi datang di PostGIS 2.3, Dan Baston akan segera memeriksa kode untuk itu - trac.osgeo.org/postgis/ticket/2259 terlihat cukup banyak dilakukan hanya perlu ditarik. Akan bagus untuk memiliki pengujian orang
LR1234567
Bisakah Anda memposting serangkaian poin yang Anda gunakan? Saya keberatan melakukan sedikit pengujian sendiri
MickyT
@MickyT Ini adalah tautan ke data saya. Data SRID adalah 2180.
DamnBack

Jawaban:

6

Meskipun ini memperbaiki masalah langsung dengan kueri untuk data yang dipermasalahkan, saya tidak senang dengan itu sebagai solusi untuk penggunaan umum dan saya akan meninjau kembali ini dan jawaban sebelumnya ketika saya bisa.

Masalahnya adalah bahwa permintaan asli menggunakan lambung cembung pada tepi voronoi untuk menentukan tepi eksterior untuk poligon voronoi. Ini berarti bahwa beberapa tepi voronoi tidak mencapai bagian luar ketika seharusnya. Saya melihat menggunakan fungsi lambung cekung, tetapi tidak benar-benar berfungsi seperti yang saya inginkan.
Sebagai perbaikan cepat saya telah mengubah lambung cembung yang akan dibangun di sekitar set tepi voronoi yang ditutup ditambah penyangga di sekitar tepi aslinya. Tepi voronoi yang tidak menutup diperpanjang keluar jarak yang besar untuk mencoba dan memastikan bahwa mereka memotong eksterior dan digunakan untuk membangun poligon. Anda mungkin ingin bermain-main dengan parameter buffer.

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;
MickyT
sumber
Terima kasih atas penjelasan dan perbaikan cepat untuk masalahnya! Ini berfungsi dengan data saya (sedikit lebih lambat karena ST_Union(ST_Buffer(geom))), tetapi saya akan melanjutkan pengujian dengan set poin lainnya. Sementara itu saya akan menunggu seperti yang Anda katakan - solusi yang lebih umum. :)
DamnBack
apakah Anda memiliki gambar yang dapat Anda poskan pada hasil akhir Anda?
Jeryl Cook
10

Mengikuti saran dari @ LR1234567 untuk mencoba fungsi ST_Voronoi baru yang telah dikembangkan oleh @dbaston, jawaban luar biasa asli @MickyT (sebagaimana dinyatakan dalam pertanyaan OP) dan menggunakan data asli sekarang dapat disederhanakan menjadi:

WITH voronoi (vor) AS 
     (SELECT ST_Dump(ST_Voronoi(ST_Collect(geom))) FROM meshpoints)
SELECT (vor).path, (vor).geom FROM voronoi;

yang menghasilkan output ini, identik dengan pertanyaan OP.

masukkan deskripsi gambar di sini

Namun, ini menderita dari masalah yang sama, sekarang diperbaiki dalam jawaban MickyT, bahwa poin pada lambung cekung tidak mendapatkan poligon terlampir (seperti yang diharapkan dari algoritma). Saya memperbaiki masalah ini dengan kueri dengan serangkaian langkah berikut.

  1. Hitung lambung cekung dari titik input - titik pada lambung cekung adalah yang memiliki poligon tidak terikat dalam diagram Voronoi keluaran.
  2. Temukan titik asli pada lambung cekung (titik kuning pada diagram 2 di bawah).
  3. Buffer lambung cekung (jarak buffer sewenang-wenang dan mungkin dapat ditemukan lebih optimal dari input data?).
  4. Temukan titik terdekat pada buffer lekukan cekung terdekat dengan titik pada langkah 2. Ini ditampilkan sebagai hijau pada diagram di bawah ini.
  5. Tambahkan titik-titik ini ke kumpulan data asli
  6. Hitung diagram Voronoi dari kumpulan data gabungan ini. Seperti dapat dilihat pada diagram ke-3, titik-titik pada lambung sekarang memiliki poligon terlampir.

Diagram 2 menunjukkan titik pada lambung cekung (kuning) dan titik terdekat ke buffer pada lambung (hijau). Diagram 2.

Query jelas dapat disederhanakan / dikompresi, tapi saya meninggalkannya sebagai bentuk CTE, karena lebih mudah untuk mengikuti langkah-langkah berurutan seperti itu. Kueri ini berjalan pada data asli yang diatur dalam milidetik (rata-rata 11ms pada server dev) sedangkan jawaban MickyT menggunakan ST_Delauney berjalan pada 4800 ms pada server yang sama. DBaston mengklaim bahwa urutan peningkatan kecepatan magnitudo lainnya dapat diperoleh dari pembangunan terhadap batang GEOS saat ini, 3.6dev, karena peningkatan dalam rutinitas triangulasi.

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;

Diagram 3 menunjukkan semua titik yang sekarang terlampir dalam poligon diagram 3

Catatan: Saat ini ST_Voronoi melibatkan pembuatan Postgis dari sumber (versi 2.3, atau trunk) dan menghubungkannya dengan GEOS 3.5 atau lebih tinggi.

Sunting: Saya baru saja melihat Postgis 2.3 saat diinstal di Amazon Web Services, dan sepertinya nama fungsinya sekarang adalah ST_VoronoiPolygons.

Tidak diragukan lagi kueri / algoritme ini dapat ditingkatkan. Saran diterima.

John Powell
sumber
@daston. Ingin tahu apakah Anda memiliki komentar tentang pendekatan ini?
John Powell
1
baik, semua poin yang mendapatkan poligon melampirkan, hanya saja itu tidak proporsional besar. Jika dan bagaimana ini harus dibuat lebih kecil cukup subyektif, dan tanpa tahu persis apa yang diinginkan untuk poligon luar sulit untuk mengetahui apa cara "terbaik". Bagiku metode yang bagus bagiku. Saya telah menggunakan metode yang kurang canggih yang serupa dengan semangat Anda, menjatuhkan poin tambahan di sepanjang batas cembung buffered pada jarak tetap yang ditentukan oleh kerapatan titik rata-rata.
dbaston
@daston. Terima kasih, hanya memastikan saya tidak melewatkan sesuatu yang jelas. Algoritma untuk mengecilkan poligon luar ke sesuatu yang lebih sesuai dengan ukuran yang dalam (dalam kasus kode pos saya) adalah sesuatu yang saya harus memikirkan lebih lanjut.
John Powell
@ John Barça Terima kasih atas solusi hebat lainnya. Kecepatan perhitungan lebih dari memuaskan dengan pendekatan ini. Sayangnya saya ingin menggunakan algoritma ini di dalam plugin QGIS saya, dan itu harus bekerja dengan PostGIS 2.1+ di luar kotak. Tapi yang pasti saya akan menggunakan solusi ini setelah rilis resmi PostGIS 2.3. Bagaimanapun terima kasih atas jawaban yang komprehensif. :)
DamnBack
@ DamnBack. Terima kasih banyak Saya membutuhkan ini untuk bekerja dan pertanyaan Anda sebenarnya banyak membantu saya, karena saya tidak tahu tentang ST_Voronoi keluar, dan solusi yang lebih lama jauh lebih lambat (seperti yang Anda perhatikan). Itu menyenangkan mencari tahu juga :-)
John Powell
3

Jika Anda memiliki akses ke PostGIS 2.3, cobalah fungsi ST_Voronoi baru, yang baru-baru ini dilakukan:

http://postgis.net/docs/manual-dev/ST_Voronoi.html

Ada build yang telah dikompilasi untuk windows - http://postgis.net/windows_downloads/

LR1234567
sumber
Terima kasih atas info bahwa ada fungsi ST_Voronoi bawaan yang baru - saya akan memeriksanya. Sayangnya saya memerlukan solusi yang berfungsi pada versi PostGIS 2.1+, jadi permintaan @MickyT adalah yang paling dekat dengan kebutuhan saya saat ini.
DamnBack
@ LR1234567. Apakah ini memerlukan versi GEOS tertentu. Saya punya waktu besok untuk menguji 2.3 dan ST_Voronoi.
John Powell
Membutuhkan GEOS 3.5
LR1234567