Quách Đồng Thắng
Trung tâm Ứng dụng GIS TP.HCM
TÓM TẮT: Voronoi diagram/ Delaunay triangulation là một trong những vấn đề cơ bản của hình học tính toán, được ứng dụng trong rất nhiều lĩnh vực, đặc biệt ứng dụng trong phân tích không gian. Bài báo trình bày khả năng xây dựng Voronoi diagram/ Delaunay triangulation trong PostgreSQL/PostGIS và trình bày các phép phân tích không gian dựa trên nền tảng Voronoi diagram/ Delaunay triangulation cho các nhà phát triển ứng dụng GIS.
Từ khóa: Phân tích không gian, PostGIS, Voronoi diagram, Delaunay triangulation, PL/Java, PL/R.
1 GIỚI THIỆU
1.1 Giới thiệu chung
Đối với các ứng dụng GIS, lớp bài toán liên quan đến Voronoi diagram/ Delaunay triangulation có rất nhiều ứng dụng thực tế cũng như hỗ trợ các phép phân tích không gian phức tạp khác. Đa số các phần mềm GIS hiện nay đều hỗ trợ công cụ cho người dùng cuối tạo Voronoi diagram/ Delaunay triangulation. Tuy nhiên, đối với các nhà phát triển ứng dụng, việc xây dựng các hàm này ở mức cơ sở dữ liệu là rất cần thiết để có thể dễ dàng kiểm soát và tái sử dụng cho các loại ứng dụng khác nhau.
PostGIS là phần mở rộng của hệ quản trị CSDL mã nguồn mở PostgreSQL cho phép quản lý dữ liệu không gian. Hiện nay PostGIS ver 2.0 chưa hỗ trợ các hàm được xây dựng sẵn để tạo Voronoi diagram/ Delaunay triangulation. Bài báo trình bày cách sử dụng các hàm procedural language như PL/Java và PL/R được xây dựng sẵn bởi bên thứ ba để xây dựng Voronoi diagram/ Delaunay triangulation cho một tập điểm trong PostGIS. Trên cơ sở đó trình bày cách khai thác Voronoi diagram/ Delaunay triangulation để giải quyết các bài toán phân tích không gian liên quan; cũng như đánh giá khả năng phát triển, mở rộng các hàm phân tích không gian trong PostGIS cho các nhà phát triển ứng dụng GIS.
1.2 Voronoi diagram
Về cơ bản, sơ đồ Voronoi được định nghĩa như sau:
- Cho P = {p1,p2,p3,…,pn} là tập n điểm trong không gian Euclidean.
- d(pi,pj): Khoảng cách Euclidean giữa pi và pj.
- Ô Voronoi của pi – kí hiệu V(pi) được định nghĩa:
V(pi) = { q : d(pi, q) < d(pj, q) ∀ j ≠ i }
- Sơ đồ Voronoi của tập điểm P, kí hiệu V(P) là hợp các ô Voronoi của tất cả các điểm thuộc P.
Nói cách khác, sơ đồ Voronoi là một phân hoạch của P thành n vùng. Mỗi vùng ứng với một và chỉ một điểm thuộc P sao cho nếu điểm q thuộc vùng ứng với pi thì khoảng cách từ q đến pi là nhỏ nhất so với các điểm khác thuộc P.
Các thuật toán xây dựng Voronoi diagram:
- Giải thuật vét cạn: Ứng với một điểm pi, lần lượt tính toán giao của tất cả các siêu phẳng của pi và
các điểm pj ≠ pi để xây dựng từng ô Voronoi cho điểm từng điểm pi. Độ phức tạp O(n2 log n).
- Giải thuật chia để trị: độ phức tạp O(n log n)
- Giải thuật dòng quét Fortune: độ phức tạp O(n log n).
1.3 Delaunay triangulation
Delaunay triangulation của một tập điểm P là lưới tam giác với đỉnh là các điểm thuộc P sao cho đường tròn ngoại tiếp của mỗi tam giác không chức bất kì điểm nào thuộc P.
Lưới tam giác Delaunay là đồ thị đối ngẫu của sơ đồ Voronoi. Giải thuật xây dựng Delaunay triagulation cho một tập điểm có độ phức tạp tương đương với Voronoi diagram là O(n logn)
Delaunay triangulation có vai trò quan trọng trong việc xây dựng mô hình TIN và trong biểu diễn các mô hình 3D dạng vector, phục vụ GIS 3D và đồ họa 3D.
1.4 Giới thiệu PostgreSQL/PostGIS
PostGIS là phần mở rộng hỗ trợ quản lý dữ liệu không gian cho hệ quản trị CSDL PostgreSQL. Hiện phiên bản PostGIS ver 2.0 có khoảng 890 hàm tương tác dữ liệu không gian tuân thủ theo chuẩn OGC Simple Feature Specificatins for SQL. Tuy nhiên, đối với lớp các bài toán phân tích không gian phức tạp như: tạo Voronoi diagram, tạo Delaunay triangulation, tìm cặp điểm gần nhất, tìm đường tròn rỗng lớn nhất,.…thì PostGIS chưa hỗ trợ.
Do đó, việc nghiên cứu, cài đặt các phép phân tích không gian ở mức CSDL PostgreSQL/PostGIS là rất cần thiết, vì các lí do sau:
- Nhà phát triển ứng dụng có thể hiện thực các phép phân tích không gian ở mức CSDL và sử dụng cho nhiều loại ứng dụng khác nhau.
- Sử dụng đơn giản thông qua truy vấn SQL.
- Kết quả các phép phân tích không gian có thể dễ dàng chuyển sang các định dạng khác nhau: Well Konwn Text – WKT, Well Known Binary – WKB, shapefile,….
2 THỬ NGHIỆM
Người dùng có thể dễ dàng bổ sung các hàm vào PostgreSQL. PostgreSQL cung cấp 4 loại function:
- Query language functions (hàm viết bằng SQL).
- Procedural language functions (được viết bằng PL/PgSQL, PL/Java, PL/R, PL/Tcl, PL/Perl, và PL/Python).
- Internal functions.
- C-language functions.
Phần sau trình bày cách tạo Voronoi Diagram và Delaunay Triangulation sử dụng procedural language functions, cụ thể như sau:
- Voronoi diagram: Sử dụng gói thư viện deldir (Delaunay Triangulation and Dirichlet (Voronoi) Tessellation) của R bằng ngôn ngữ PL/R.
- Delaunay Triangulation: Sử dụng hàm ST_DelaunayTriangles() được xây dựng sẵn trong Jaspa bằng ngôn ngữ PL/Java.
2.1 Xây dựng Voronoi diagram sử dụng PL/R và deldir Package của R
R là phần mềm mã nguồn mở hỗ trợ phân tích thống kê và hiển thị đồ họa mạnh mẽ.
PL/R là một phần mở rộng ngôn ngữ PostgreSQL cho phép viết các hàm PostgreSQL bằng ngôn ngữ R.
Deldir package (Delaunay Triangulation and Dirichlet (Voronoi) Tessellation) của tác giả Rolf Turner cho phép xây dựng Voronoi diagram/ Delaunay Triangulation từ một tập điểm.
Giải pháp đưa ra là sử dụng function được viết bằng PL/R trong PostgreSQL/PostGIS sử dụng gói deldir trong R để xây dựng Voronoi diagram cho một tập điểm trong PostgreSQL/PostGIS.
Các bước chuẩn bị (Sử dụng PostgreSQL ver 9.1, PostGIS ver 2.0, R ver 2.15.1):
- Cài đặt PostgreSQL/ PostGIS, tạo CSDL với template postgis và bảng dữ liệu “nodes” chứa đối tượng điểm.
- Cài đặt R và deldir package.
- Cài đặt PL/R (For Windows user)
- Download PL/R, copy plr.dll vào thư mục lib của PostgreSQL.
- Tạo biến môi trường R_HOME (C:\Program Files\R\R-2.15.1).
- Bổ sung vào biến Path đường dẫn của R (C:\Program Files\R\R-2.15.1\bin\i386\).
- Khởi động lại dịch vụ PostgreSQL.
Tạo Voronoi diagram sử dụng PL/R:
1. Thực thi script file plr.sql (trong thư mục PL/R).
2. Chạy file voronoi.sql (tác giả Mike Leahy) để tạo hàm voronoi trong PostgreSQL/ PostGIS trên cơ sở sử dụng voronoi function từ package deldir.
2. Chạy file voronoi.sql (tác giả Mike Leahy) để tạo hàm voronoi trong PostgreSQL/ PostGIS trên cơ sở sử dụng voronoi function từ package deldir.
create type voronoi as (id integer, polygon geometry);
create or replace function voronoi(text,text,text) returns setof voronoi as '
library(deldir)
# select the point x/y coordinates into a data frame...
points <- pg.spi.exec(sprintf("select x(%2$s) as x, y(%2$s) as y from %1$s;",arg1,arg2))
# calculate an approprate buffer distance (~10%):
buffer = ((abs(max(points$x)-min(points$x))+abs(max(points$y)-min(points$y)))/2)*(0.10)
# get EWKB for the overall buffer of the convex hull for all points:
buffer <- pg.spi.exec(sprintf("select buffer(convexhull(st_union(%2$s)),%3$.6f) as ewkb from %1$s;",arg1,arg2,buffer))
# the following use of deldir uses high precision and digits to prevent slivers between the output polygons, and uses
# a relatively large bounding box with four dummy points included to ensure that points in the peripheral areas of the
# dataset are appropriately enveloped by their corresponding polygons:
voro = deldir(points$x, points$y, digits=22, frac=0.00000000000000000000000001,list(ndx=2,ndy=2), rw=c(min(points$x)-abs(min(points$x)-max(points$x)), max(points$x)+abs(min(points$x)-max(points$x)), min(points$y)-abs(min(points$y)-max(points$y)), max(points$y)+abs(min(points$y)-max(points$y))))
tiles = tile.list(voro)
poly = array()
id = array()
p = 1
for (i in 1:length(tiles)) {
tile = tiles[[i]]
curpoly = "POLYGON(("
for (j in 1:length(tile$x)) {
curpoly = sprintf("%s %.6f %.6f,",curpoly,tile$x[[j]],tile$y[[j]])
}
curpoly = sprintf("%s %.6f %.6f))",curpoly,tile$x[[1]],tile$y[[1]])
# this bit will find the original point that corresponds to the current polygon, along with its id and the SRID used for the
# point geometry (presumably this is the same for all points)...this will also filter out the extra polygons created for the
# four dummy points, as they will not return a result from this query:
ipoint <- pg.spi.exec(sprintf("select %3$s as id, intersection(''SRID=''||srid(%2$s)||'';%4$s'',''%5$s'') as polygon from %1$s where intersects(%2$s,''SRID=''||srid(%2$s)||'';%4$s'');",arg1,arg2,arg3,curpoly,buffer$ewkb[1]))
if (length(ipoint) > 0)
{
poly[[p]] <- ipoint$polygon[1]
id[[p]] <- ipoint$id[1]
p = (p + 1)
}
}
return(data.frame(id,poly))
' language 'plr';
3. Sử dụng hàm voronoi để tạo voronoi diagram từ bảng nodes.
create table voronoidiagram (gid serial primary key);
SELECT AddGeometryColumn ('public','voronoidiagram','the_geom',4326,'POLYGON',2);
insert into voronoidiagram (gid, the_geom)
SELECT v.id, v.polygon
FROM voronoi('(SELECT gid, the_geom FROM nodes ) as vor', 'vor.the_geom', 'vor.gid') As v
2.2 Xây dựng Delaunay triangulation sử dụng PL/Java và Jaspa
Jaspa (Java Spatial) là phần mở rộng cho phép quản lý dữ liệu không gian trên PostgresQL và H2. Giống như Postgis, Jaspa tạo một template định nghĩa các hàm hỗ trợ quản lý dữ liệu không gian tuân thủ đặc tả OGC Simple Feature Specificatins for SQL. Jaspa quản lý dữ liệu trong schema jaspa, trong đó có hàmST_DelaunayTriangles cho phép xây dựng lưới tam giác Delaunay từ một lớp đối tượng dạng điểm.
Giả sử ta có một database PostGIS với tên gọi voronoi chứa table nodes cần xây dựng Delaunay triagulation. Chúng ta sẽ cài đặt jaspa và sử dụng hàm ST_DelaunayTriangles để xây dựng lưới tam giác cho table nodes trong database PostGIS.
Cài đặt Jaspa:
- Download jaspa và giải nén vào thư mục C:\j aspa4pg . Cấu trúc thư mục như sau:
- C:\jaspa4pg\lib
- C:\jaspa4pg\bin
- C:\jaspa4pg\sql
- C:\jaspa4pg\doc
- C:\jaspa4pg\jre
- Bổ sung vào biến Path trong System Variable:
C:\jaspa4pg\jre\bin;C:\jaspa4pg\jre\bin\client;c:\jaspa4pg\bin;C:\jaspa4pg\jre\bin\server; ; C:\Program Files\PostgreSQL\9.1\bin
- Cấu hình PL/Java: Copy nội dung file C:\jaspa4pg\doc\html\ pgconf-wint.txt và paste vào cuối nội dung file C:\Program Files\PostgreSQL\9.1\data\ postgresql.conf
- Khởi động lại dịch vụ PostgreSQL.
Cài đặt jaspa vào database PostGIS vornoi:
- Chạy script file C:\jaspa4pg\sql\pljavainstall.sql để cài đặt PL/Java.
- Chạy script file C:\jaspa4pg\sql\jaspa.sql để định nghĩa các hàm tương tác dữ liệu không gian trong jaspatemplate.
- Kết quả là schema với tên jaspa được tạo ra với các hàm hỗ trợ quản lý dữ liệu không gian, trong đó có hàm ST_DelaunayTriangles:
Sử dụng hàm ST_DelaunayTriangles của jaspa để xây dựng lưới tam giác Delaunay cho table “nodes” trong PosGIS:
-- Create table delaunay to hold the Delaunay Trianagulation
create table delaunay (gid serial primary key);
SELECT AddGeometryColumn ('public','delaunay','the_geom',4326,'POLYGON',2);
-- Create delaunay triangulation using jaspa.ST_DelaunayTriangles function.
insert into delaunay (the_geom)
SELECT ST_GeomFromEWKB(ST_AsEWKB(jaspa.ST_Dump(jaspa.ST_DelaunayTriangles(jaspa.ST_GeomFromEWKB(ST_AsEWKB(ST_Collect(the_geom))))))) from nodes
3 PHÂN TÍCH KHÔNG GIAN TRONG POSTGIS DỰA TRÊN VORONOI DIAGRAM/ DELAUNAY TRIANGULATION
Xét bài toán tìm cặp điểm gần nhất của một tập điểm cho trước. Nếu sử dụng phương pháp “vét cạn”: lần lượt duyệt qua các cặp điểm, cho đến khi tìm được cặp điểm có khoảng cách gần nhất. Sử dụng truy vấn SQL như sau:
SELECT n1.gid, (select gid from nodes n2 where n1.gid <> n2.gid order by ST_Distance(n1.the_geom, n2.the_geom) limit 1 ),
(select ST_Distance(n1.the_geom, n2.the_geom) from nodes n2 where n1.gid <> n2.gid order by ST_Distance(n1.the_geom, n2.the_geom) limit 1 ) as dist
FROM nodes n1
order by dist
Thử nghiệm trong PostgreSQL, với số điểm là n = 4469, thời gian thực thi là khoảng 60s. Dễ thấy, với phương án “vét cạn”, khi n tăng, thời gian thực thi sẽ tăng lên rất nhiều (cỡ O(n2)). Trong khi đó, lớp bài toán này có thể cài đặt các thuật toán thông minh hơn để giải quyết, hoặc sử dụng kết quả của việc xây dựng sơ đồ Voronoi để truy vấn (ứng với mỗi điểm, chỉ tính khoảng cách đến các điểm khác thuộc ô Voronoi liền kề - hay nói cách khác chính là tìm ra cạnh Delaunay nhỏ nhất), lúc này thời gian và tài nguyên sử dụng trong câu truy vấn được tiết kiệm rất đáng kể.
Ví dụ: tìm cặp điểm gần nhất bằng cách sử dụng cạnh Delaunay (đồ thị đối ngẫu của Voronoi diagram)
SELECT n.gid from nodes n
where ST_touches(n.the_geom, (select the_geom from delaunay_line d order by ST_length(the_geom) limit 1))
Các tính chất của Voronoi diagram được sử dụng để giải quyết nhiều bài toán khác nhau:
- Tìm điểm gần nhất trong một tập điểm P cho trước so với một điểm q: điểm pi thuộc ô Voronoi V(pi) chứa q chính là điểm gần nhất so với p.
- Tìm cặp điểm gần nhất trong một tập điểm P cho trước: ứng với mỗi điểm, chỉ cần tính toán khoảng cách so với các điểm thuộc ô Voronoi liền kề - hay là cạnh của tam giác Delaunay (thay vì tính toán khoảng cách so với n-1 điểm còn lại) để tìm ra cặp điểm gần nhất.
- Tìm k láng giềng gần nhất cho tất cả các điểm trong một tập điểm cho trước.
- Tìm medial axis/ skeleton cho polygon: ứng dụng trong việc tạo tim đường giao thông/ sông hồ dựa trên lớp vùng giao thông/ sông hồ.
- Tìm cây bao trùm nhỏ nhất: Sử dụng đồ thị đối ngẫu của Voronoi diagram là Delaunay triangulation, kết hợp thuật toán Kruskal để tìm cây bao trùm nhỏ nhất.
- Tìm đường tròn rỗng lớn nhất: chỉ cần so sánh các đường tròn có tâm là đỉnh các ô Voronoi.
- Phân cụm (clustering).
- Phân lớp (classification).
Một số ứng dụng thực tế của Voronoi diagram/ Delaunay triangulation:
- Kinh doanh: Một hệ thống siêu thị có nhiều chi nhánh tọa lạc tại các địa điểm khác nhau. Hệ thống siêu thị này cần gửi các poster quảng cáo đến từng địa chỉ nhà của các khách hàng. Vấn đề đặt ra là cần phân vùng cho các chi nhánh sao cho khoảng cách từ chi nhánh đến địa chỉ nhà thuộc phân vùng của nó là nhỏ nhất so với đi từ các chi nhánh khác. Cách phân vùng này chính là xây dựng sơ đồ Voronoi cho các chi nhánh.
- Thủy văn: Sơ đồ Voronoi được ứng dụng để tính lượng mưa bình quân lưu vực (còn được gọi là đa giác Thiessen).
- Dịch tể học: xác định vùng ảnh hưởng của một ca bệnh lây nhiễm.
- Môi trường: Cần xây dựng một khu xử lý rác trong phạm vi một khu dân cư. Yêu cầu đặt ra là chọn vị trí xây dựng sao cho cách xa các hộ dân càng tốt. Bài toán này có thể giải quyết bằng cách tìm đường tròn rỗng lớn nhất trên cơ sở Voronoi diagram.
Một số vấn đề mở rộng của Voronoi diagram:
- Voronoi Diagram cho các điểm xa nhau nhất, được ứng dụng trong bài toán tìm đường tròn bao nhỏ nhất của một tập điểm.
- Voronoi diagram mức 2: ô Voronoi V(pi,pj) chứa tập các điểm gần 2 điểm pi và pj.
- Voronoi diagram mức k: ô Voronoi V(pi,pj, … pk) chứa tập các điểm gần k điểm pi, pj,… pk nhất.
- Voronoi diagram cho tập điểm có trọng số.
- Voronoi diagram cho tập các đoạn thẳng.
- 3D Voronoi diagram.
4 KẾT LUẬN
Bài viết đã giới thiệu một số “thủ thuật” để xây dựng Voronoi diagram/ Delaunay triangulation trong PostGIS và giải pháp khai thác mở rộng một số chức năng phân tích không gian phức tạp trong PostGIS. Do đó, các vấn đề như thuật toán xây dựng và các vấn đề nâng cao của Voronoi diagram/ Delaunay triangulation không được đề cập sâu trong bài viết.
Phần thử nghiệm trong nghiên cứu này cho thấy khả năng sử dụng các hàm được xây dựng sẵn bở bên thứ ba để mở rộng khả năng phân tích không gian trong PostGIS (cũng như của PostgreSQL nói chung). Cụ thể là:
- Giữa PostGIS và R thông qua ngôn ngữ PL/R: sữ dụng gói deldir package (tác giả Rolf Turner) và hàm voronoi sử dụng ngôn ngữ PL/R (tác giả Mike Leahy)
- Giữa PostGIS và jaspa thông qua ngôn ngữ SQL: sử dụng gói jaspa với hàm ST_DelaunayTriangles() được xây dựng sẵn trong jaspa template.
Cách tiếp cận trong bài viết có thể khá phức tạp đối với người dùng cuối nhưng là một cách tiếp cận rất hữu ích và thú vị cho các nhà phát triển ứng dụng GIS trong việc mở rộng các hàm phân tích không gian ở mức cơ sở dữ liệu để sử dụng trong các loại ứng dụng khác nhau.
TÀI LIỆU THAM KHẢO
1) Boissonnat, J.-D. (January 2010). Convex Hulls, Voronoi Diagrams and Delaunay Triangulations. ENS-Lyon.
2) Edwards, R. (December 2, 2010). Determining the Skeleton of a Simple Polygon in (Almost) Linear Time. Oak Ridge, Tennessee.
3) Francis Chin, Jack Snoeyinky, Cao An Wangz. Finding the Medial Axis of a Simple Polygon.
4) Franz Aurenhamme, Rolf Klein. Voronoi Diagrams.
5) Goswami, P. P. Introduction to Computational Geometry.
6) Inkulu, R. Voronoi diagrams: Higher order.
7) Kreveld, M. v. Computational Geometry: its objectives and relation to GIS.
8) Leo Hsu and Regina Obe. (n.d.). bostongis. Retrieved from PL/R and PostGIS: http://www.bostongis.com/?content_name=postgresql_plr_tut02#98
9) Maksym Zavershynskyi, Prof. Evanthia Papadopoulou. Higher-Order Voronoi Diagram of Line Segments.
10) Mark de Berge, Otfried Cheong, Marc vn Kreveld, Mark Overmars. (2008). Computational Geometry - Algorithms and applications. Springer.
11) Martinez, P. J. (n.d.). Retrieved from http://jaspa.upv.es/blog/
12) Nandy, S. C. Voronoi Diagram.
13) Nehab, D. Medial Axis.
14) Obe, L. H. Cross Comparison of Spatially enabled databases: PostGIS, SQL server and Jaspa.
15) Schuster, M. The Largest Empty Circle Problem.
16) Sharifzadeh, M. (2007). Spatial Query Processing using Voronoi diagram.
17) Smid, M. (1997). Closest-Point Problems in Computational Geometry. Magdeburg, Germany.
18) Turner, R. (n.d.). Retrieved from http://cran.r-project.org/web/packages/deldir/index.html
19) Vennemann, K. Beyond PostGIS: new developments in opens source spatial databases.
Great post! Rất tốt cho anh em nghiên cứu và ứng dụng.
ReplyDelete