Ứng dụng phần mềm mã nguồn mở hỗ trợ giám sát dịch bệnh nguy hiểm

Nguyễn Đức Tuấn, Quách Đồng Thắng

Trung tâm Ứng dụng GIS TP.HCM

TÓM TẮT: Các phần mềm thống kê, phần mềm GIS và mô hình hóa ngày càng trở thành những công cụ hỗ trợ đắc lực, không thể thiếu cho các nhà nghiên cứu dịch tễ học và y tế công cộng. Bên cạnh các phần mềm thương mại, phần mềm mã nguồn mở là một lựa chọn thích hợp nhằm tiết kiệm chi phí đầu tư và nâng cao tính chủ động về công nghệ. Bài viết giới thiệu giải pháp sử dụng phần mềm mã nguồn mở hỗ trợ công tác giám sát dịch bệnh nguy hiểm, cũng như trong nghiên cứu dịch tễ học nói chung.

Từ khóa: Y tế công cộng, dịch tễ học, dịch bệnh, thống kê, hệ thống thông tin địa lý, mô hình tác nhân, R, GAMA, gvSIG.

1          GIỚI THIỆU

Vào năm 1854, trận dịch tả lớn bùng phát ở London đã cướp đi sinh mạng 1/8 dân số của thành phố khi bác sĩ John Snow (1813 – 1858), bằng cách biểu diễn mỗi ca bệnh là một điểm trên bản đồ vẽ tay, đã phát hiện ra nguồn gốc của dịch bệnh là từ các trụ bơm cung cấp nước sinh hoạt bị nhiễm khuẩn trên đường Broad Street. Ông đã đưa ra giải pháp cho ngưng hoạt động các trụ bơm bị nhiễm khuẩn, và đã ngăn được dịch bệnh bùng phát.

Với việc thu thập, vẽ bản đồ, phân tích số liệu và chứng minh rằng nguồn nước bị nhiễm khuẩn chính là nguồn gốc của sự bùng phát dịch tả ở London - chứ không phải do nguyên nhân “khí độc” như quan niệm của các quan chức y tế thời đó, John Snow được xem như cha đẻ của ngành dịch tễ học thực địa và ngành phân tích không gian.

Ngày nay, dịch bệnh ngày càng diễn biến phức tạp và nguy hiểm. Việc nghiên cứu, giám sát và kiểm soát dịch bệnh là rất quan trọng và đòi hỏi kiến thức chuyên môn của các nhà dịch tễ học. Ngoài ra, công tác nghiên cứu dịch tễ học nói chung và giám sát dịch bệnh nguy hiểm nói riêng đòi hỏi có sự hỗ trợ của phần mềm máy tính như phần mềm thống kê, phần mềm GIS và phần mềm mô hình hóa. Sử dụng các công cụ này giúp các nhà nghiên cứu dễ dàng trong việc thu thập, phân tích, thống kê số liệu, thực hiện các phép phân tích không gian, mô hình hóa và mô phỏng dịch bệnh, hỗ trợ công tác giám sát và dự báo tình hình dịch bệnh nhằm đưa ra những biện pháp ứng phó kịp thời và hiệu quả.

2          CÁC PHẦN MỀM MÃ NGUỒN MỞ HỖ TRỢ NGHIÊN CỨU DỊCH TỄ HỌC

2.1        Phần mềm thống kê

Phần mềm R/ RStudio:

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ẽ.  Ngoài các chức năng thống kê nói chung, R hỗ trợ các gói thư viện (package) được xây dựng sẵn để làm việc với dữ liệu dịch tễ nhưsurveilance, R0, RladyBug, EpiEstim, mem, epinet, flubase.

RStudio: Là một môi trường cung cấp giao diện đồ họa giúp người dùng làm việc dễ dàng hơn với R.

Ngoài ra còn có các phần mềm chuyên dụng cho dịch tễ học như:

-         Epi Info.

-         OpenEpi.

2.2        Phần mềm mô hình hóa

Mô hình hóa giúp đánh giá tình hình dịch bệnh theo các kịch bản khác nhau, từ đó đưa ra những dự báo và hỗ trợ công tác kiểm soát, ứng phó với tình hình dịch bệnh. Có nhiều phương pháp mô hình hóa trong nghiên cứu dịch tễ học (nói cách khác là mô hình hóa sự lây lan dịch bệnh trên một dân số nhất định). Bài viết trình bày 02 cách tiếp cận mô hình hóa dịch bệnh: mô hình toán (mathematical model) và mô hình tác nhân (agent-based model).

2.2.1       Mô hình toán SIR

Một dạng mô hình toán thông dụng là mô hình chia ngăn tất định (Deterministic compartmental models). Mô hình này chia dân số ra thành các lớp (ngăn) biểu diễn các giai đoạn/ trạng thái khác  nhau của dịch bệnh.

Mô hình chia ngăn tất định tiêu biểu là mô hình SIR (Suscceptible – Infectious – Recover/ immune/ expired) được Kermack và McKendrick đề xuất vào năm 1927.

Mô hình SIR giả sử rằng tất cả các cá thể trong cộng đồng đều có khả năng bị nhiễm bệnh như nhau và hoàn toàn miễn dịch sau giai đoạn bị lây nhiễm. Dân số được chia thành 3 lớp riêng biệt: S – người có khả năng bị nhiễm bệnh, I – người đã bị lây nhiễm và có thể lây cho người khác, R – người đã từng bị nhiễm và hiện tại miễn nhiễm (hồi phục hoặc tử vong).

Cho S(t), I(t) và R(t) lần lượt là số lượng người có khả năng bị nhiễm, người bị nhiễm và người miễn nhiễm tại thời điểm t và N là cỡ dân số:

N = S(t) + I(t) + R(t)

Giả sử tốc độ người khỏe mạnh mắc bệnh khi tiếp xúc với người bị nhiễm là α, và người bị nhiễm sẽ được miễn nhiễm với tốc độ là β. Mô hình SIR được biểu diễn dưới dạng các phương trình vi phân sau:

Ngoài ra còn có các biến thể từ mô hình SIR như mô hình SI, SIS, SIRS, SEIR (Susceptible – Exposed –  Infectious –Recovered), MSEIR (Immunized – Susceptible – Exposed – Infectious – Recovered). Tuy nhiên, các mô hình toán học có nhược điểm là không tính đến yếu tố thay đổi theo không gian và thời gian như sự thay đổi cấu trúc dân số và tính tương tác giữa các cá thể.

2.2.2       Mô hình tác nhân

Phần mềm mô hình hóa, đặc biệt là mô hình tác nhân Agent-based Modeling (ABM) được ứng dụng trong rất nhiều lĩnh vực, trong đó có xây dựng mô hình và mô phỏng dịch bệnh. Trong khi mô hình chia ngăn không thể cài đặt các hành vi phức tạp như tính di động và tính tương tác giữa các cá thể, mô hình tác nhân có thể khắc phục nhược điểm này bằng cách cài đặt các hành vi đặc trưng của từng cá thể - đóng vai trò là một tác nhân trong mô hình. Ví dụ như người nhiễm bệnh có các mối quan hệ xã hội, giao tiếp với người ở khu vực địa lý nào, từ đó có thể dự báo tình hình lây lan dịch bệnh thông qua các hành vi của từng ca bệnh.

Một số phần mềm mô hình tác nhân mã nguồn mở thông dụng:

-         Repast / Agent Analysist

-         GAMA

-         NetLogo

2.3        Phần mềm GIS

Các ca bệnh luôn gắn với một vị trí không gian trong mối tương quan với các đối tượng khác, và là một yếu tố quan trọng trong các nghiên cứu dịch tễ. Với sự hỗ trợ của GIS, đặc biệt là các phép phân tích không gian, các nhà dịch tễ học sẽ có công cụ hữu hiệu trong điều tra, thu thập, hiển thị dữ liệu và kiểm soát dịch bệnh.

Một số khái niệm:

-     Geographical Information System (GIS): Hệ thống thu thập, lưu trữ, xử lý và kết xuất dữ liệu không gian.

-     Dữ liệu không gian: Thông tin về vị trí, hình dạng và các thuộc tính mô tả của đối tượng.

-     Luật Tobler 1 về địa lý: “Mọi vật đều có quan hệ với những vật khác, tuy nhiên những vật gần nhau có quan hệ nhiều hơn những vật xa nhau.” (“Everything is related to everything else, but near things are more related than distant things.’’ – Tobler, 1970). Luật Tobler 1 rất có ý nghĩa trong nghiên cứu dịch tễ học.

Nguồn dữ liệu GIS hỗ trợ các phân tích dịch tễ học:

-     Dữ liệu dân số trong mối tương quan với các vấn đề kinh tế - xã hội.

-     Dữ liệu môi trường, sinh thái: tình trạng ô nhiễm, thực phủ (từ ảnh viễn thám).

-     Dữ liệu địa hình, thủy văn, khí hậu.

-     Dữ liệu sử dụng đất và cơ sở hạ tầng kỹ thuật: giao thông, trường học, hệ thống cấp thoát nước.

-     Dữ liệu cơ sở y tế và dữ liệu dịch tễ học: các ca nhiễm bệnh, các ca tử vong,…

-     Phân bố dịch bệnh và mạng lưới các cơ sở y tế.

-    Các nguồn dữ liệu hỗ trợ khác: ảnh viễn thám, GPS (phục vụ thu thập dữ liệu về vị trí các ca bệnh tại thực địa).

3          GIỚI THIỆU GIẢI PHÁP CỦA TRUNG TÂM ỨNG DỤNG GIS TP.HCM

3.1        Sử dụng phần mềm mã nguồn mở

-    Phần mềm thống kê: R/ RStudio và các package hỗ trợ nghiên cứu dịch tễ.

-   Phần mềm GIS: gvSIG, gvSIG mobile, Quantum GIS: thu thập, lưu trữ, phân tích, kết xuất dữ liệu không gian như bản đồ phân bố dịch bệnh, hệ thống giao thông, ranh giới hành chính,…

-    Phần mềm mô hình hóa tác nhân: GAMA được dùng để xây dựng mô hình và mô phỏng dịch bệnh, có thể tích hợp dữ liệu GIS.

-    Hệ quản trị CSDL: PostgreSQL/ PostGIS để quản lý cơ sở dữ liệu khi triển khai hệ thống đa người dùng.

-    Dữ liệu GIS: CSDL GIS nền TP.HCM do Trung tâm Ứng dụng GIS xây dựng và cập nhật, gồm các lớp ranh giới hành chính, giao thông, thủy hệ, thửa đất,...

Các hình sau minh họa bản đồ dịch tả của John Snow đã được số hóa trong gvSIG và Quantum GIS: các trụ nước và ca bệnh tả được số hóa dưới dạng điểm trên nền lớp dữ liệu giao thông.

Hình 3: Bản đồ phân bố dịch bệnh tả của John Snow trên gvSIG

Với sự hỗ trợ của các phép phân tích không gian trong phần mềm GIS, việc phân tích dữ liệu dịch tễ được thực hiện dễ dàng, trực quan và chính xác. Ví dụ sau cho thấy biểu đồ Voronoi được xây dựng để phân hoạch “vùng ảnh hưởng” của từng trụ bơm nước. Khi lớp “vùng ảnh hưởng” này được chồng lên vị trí các ca bệnh sẽ cho phép xác định nhanh chóng các ca bệnh nào gần trụ bơm nước nào nhất (do đó có khả năng bị “ảnh hưởng” từ trụ nước đó nhất).

3.2        Thử nghiệm xây dựng mô hình SIR bằng phần mềm GAMA

GAMA là phần mềm mã nguồn mở được phát triển trên nền Java, được ứng dụng trong rất nhiều lĩnh vực, trong đó có lĩnh vực nghiên cứu dịch tễ học và mô phỏng dịch bệnh. Người dùng có thể tự xây dựng, định nghĩa hành vi của các tác nhân trong mô hình và mô phỏng (chạy mô hình) thông qua ngôn ngữ mô hình hóa GAML trong môi trường Eclipse IDE. GAMA hỗ trợ làm việc với dữ liệu GIS dạng shapefile, giúp người dùng có thể đưa các dữ liệu thực tế vào mô hình như hệ thống giao thông, sông hồ, ranh giới hành chính, bản đồ phân bố dịch bệnh,…

GAMA giúp các nhà nghiên cứu dễ dàng trong việc:

-    Mô hình hóa các hệ thống phức tạp một cách trực quan và linh hoạt.

-    Mô phỏng (chạy mô hình).

-    Phân tích mô hình và các kết quả mô phỏng.

-   Thao tác dữ liệu GIS một cách đơn giản, dễ dàng chuyển dữ liệu GIS thành các tác nhân trong mô hình

-   Có thể kế thừa các mô hình khác đã được xây dựng sẵn.

Phần sau trình bày cách mô hình hóa mô hình toán học SIR trong phần mềm GAMA như là một cách giới thiệu cách sử dụng GAMA trong việc xây dựng các mô hình, mô phỏng, dự báo tình hình dịch bệnh.

Trong GAMA, mô hình được xây dựng bằng ngôn ngữ GAML (file *.gaml). Cấu trúc một GAMA project gồm các thư mục như doc, images, includes và models.

Cấu trúc một mô hình (file*.gaml) được xây dựng trong GAMA gồm:

-   Global: Định nghĩa các biến dùng chung, các hàm khởi tạo khi chạy mô hình.

-   Environment: Định nghĩa phạm vi khung nhìn của mô hình – hữu dụng khi tích hợp dữ liệu GIS vào mô hình.

-   Entities: Định nghĩa các tác nhân trong mô hình qua từ khóa species.

-   Experiment: Định nghĩa các hiển thị khi chạy mô hình.

Thử nghiệm xây dựng mô hình và mô phỏng SIR bằng cách tạo file SIR.gaml định nghĩa các nhóm S, I, R, các màn hình hiển thị kết quả mô phỏng và các tham số của mô hình. Các tham số có thể hiệu chỉnh khi chạy mô hình là:

-   Number of agents: Cỡ dân số.

-   Number of Infected: Số ca nhiễm.

-   Hệ số α, β.

Kết quả chạy mô hình SIR trong GAMA

Đây là một ví dụ rất đơn giản về cách sử dụng GAMA để xây dựng mô hình và mô phỏng dịch bệnh. Trong thực tế, để xây dựng mô hình một cách hiệu quả và gần với thế giới thực, cần kết hợp sử dụng các nguồn dữ liệu GIS như vị trí các ca bệnh được thu thập bằng GPS và số hóa thành dữ liệu GIS, sử dụng thành thạo công cụ GAMA để có thể định nghĩa các tác nhân tham gia vào mô hình, cũng như thực hiện các kịch bản mô phỏng và phân tích kết quả mô phỏng.

Hình sau cho thấy việc tích hợp dữ liệu GIS vào mô hình SIR trong GAMA để mô phỏng lây lan dịch bệnh trong một khu vực địa lý nhất định.

4          KẾT LUẬN

Ứng dụng các phần mềm, đặc biệt là các phần mềm mã nguồn mở thống kê, GIS và mô hình hóa trong nghiên cứu dịch tễ học là rất cần thiết nhằm tiết kiệm thời gian (thực hiện các phép tính thống kê, chạy mô hình) và tiết kiệm chi phí đầu tư.

Bài viết đã giới thiệu một số phần mềm mã nguồn mở hữu ích, đã được sử dụng rất rộng rãi trong nhiều lĩnh vực, trong đó có lĩnh vực nghiên cứu dịch tễ học và giám sát dịch bệnh.

Việc kết hợp mô hình tác nhân và dữ liệu GIS giúp việc mô hình hóa dịch bệnh mang tính trực quan và tính thực tiễn. Bên cạnh các kiến thức chuyên môn về dịch tễ học, nhà nghiên cứu dịch tễ cần đầu tư các nghiên cứu sâu hơn về các phần mềm mã nguồn mở được đề xuất trong bài viết như tìm hiểu ngôn ngữ GAML trong phần mềm mô hình tác nhân GAMA, sử dụng các phần mềm gvSIG/ gvSIG Mobile/ Quantum GIS, các kỹ thuật cải thiện tốc độ mô hình hóa và mô phỏng với một khối lượng dữ liệu lớn.

5         TÀI LIỆU THAM KHẢO

1)      al., A. D. (n.d.). GAMA platform. Retrieved from http://code.google.com/p/gama-platform/

2)      Boulos, M. N. Geographic Informatics in Health. Bath BA2 7AY, UK.

3)      Dionne C.G. Law, Rachel A. Wilfert. FOCUS on Field epidemiology. North Carolina Center for Public Health Preparedness.

4)      Enrique Frías-Martínez, G. W.-M. An Agent-Based Model of Epidemic Spread using Human Mobility and Social Network Information.

5)      Greenley, B. An Agent-Based Model of Recurring Epidemics in a Population with Quarantine Capabilities .

6)      Iannelli, M. The mathematicalmodeling of epidemics.

7)      JJ Allaire, Joe Cheng, Josh Paulson, Paul DiCristina. (n.d.). RStudio. Retrieved from http://rstudio.org/

8)      Khaled M. Khalil, M. Abdel-Aziz, Taymour T. Nazmy, Abdel-Badeeh M. Salem. An Agent-Based Modeling for Pandemic Influenza in Egypt.

9)      Peter Dawson, Alex Skvortsov, Russell Connell, Ralph Gailis. (2007). Epidemic Spread Modelling: Alignment of Agent-based Simulation with a SIR Mathematical Model.

10)  Russell Connell, Peter Dawson and Alex Skvortsov. Comparison of an Agent-based Model of Disease.

11)  Tassier, T. (2005). SIR Model of Epidemics.

12)  Wu, J. T. (2010). Basic epidemic theory.

Phân tích không gian trong PostGIS dựa trên Voronoi diagram/ Delaunay triangulation

Xem trên Asia GeoSpatial Forum 2012: Full paper (pdf), Presentation (pps)

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 = {p₁,p₂,p₃,…,p_n} là tập n điểm trong không gian Euclidean.

-       d(p_i,p_j): Khoảng cách Euclidean giữa p_i và p_j.

-       Ô Voronoi của p_i – kí hiệu V(p_i) được định nghĩa:

V(p_i) = { q : d(p_i, q) < d(p_j, 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 p_i thì khoảng cách từ q đến p_i 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(n² 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.

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(n²)). 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.