Visulalization Boost Voronoi in OpenSceneGraph
eryar@163.com
Abstract. One of the important features of the boost polygon library is the implementation of the generic sweepline algorithm to construct Voronoi diagrams of points and linear segments in 2D(developed as part of the Google Summer of Code 2010 program). Voronoi diagram data structure has applications in image segmentation, optical character recognition, nearest neighbor queries execution. It is closely related with the other computational geometry conectps: Delaunay triangulation, medial axis, straight skeleton, the largest empty circle. The paper focus on the usage of Boost.Polygon Voronoi Diagram and visualize it in OpenSceneGraph.
Key words. Voronoi, Boost.Polygon, C++, OpenSceneGraph, Visualization
1. Introduction
計算幾何（Computational Geometry）作為一門學科，起源于20世紀70年代，經過近四十多年的發展，其研究內容不斷擴大，涉及Voronoi圖、三角剖分、凸包、直線與多邊形求交、可見性、路徑規劃、多邊形剖分等內容。據相關統計，在數以千計的相關文章中，約有15%是關于Voronoi圖及其對偶（dual）圖Delaunay三角剖分（Delaunay Triangulation）的研究。由于Voronoi圖具有最近性、鄰接性等眾多性質和比較系統的理論體系，如今已經在計算機圖形學、機械工程、地理信息系統、機器人、圖像處理、大數據分析與處理、生物計算及無線傳感網絡等領域得到了廣泛應用，同時也是解決碰撞檢測、路徑規劃、可見性計算、骨架計算以及凸包計算等計算幾何所涉及的其他問題的有效工具。
Voronoi圖的起源最早可以追溯到17世紀。1644年，Descartes用類似Voronoi圖的結構顯示太陽系中物質的分布。數學家G.L. Dirichlet和M.G.Voronoi分別于1850年和1908年在他們的論文中討論了Voronoi圖的概念，所以Voronoi圖又叫Dirichlet tessellation。在其他領域，這個概念也曾獨立地出現，如生物學和生理學中稱之為中軸變換（Medial Axis Transform）或骨架（Skeleton）?；瘜W與物理學中稱之為Wigner-Seitz Zones，氣象學與地理學中稱之為Thiessen多邊形。Voronoi圖最早由Thiessen應用于氣象觀測站中隨機分布的研究。由于M.G. Voronoi從更通用的n維情況對其進行研究和定義，所以Voronoi圖這個名稱為大多數人所使用。
在路徑規劃、機械加工、模式識別、虛擬現實、生物計算等領域，將站點從離散點擴展到線段圓弧等生成Voronoi圖的方式也是非常常見的。
目前可用于生成Voronoi圖的庫有一些，很多是開源庫。像CGAL庫、boost中也提供了生成Voronoi圖的算法。本文根據Boost.Polygon中的Voronoi庫，并用OpenSceneGraph顯示出剖分結果。
2. Boost.Polygon Voronoi Diagram
Boost.Polygon庫提供了構造Voronoi圖的接口，可根據點集、線段集來生成Voronoi圖，如下圖所示：
Figure 2.1 Voronoi Diagram generated by Boost.Polygon Voronoi Algorithms
Boost.Polygon中的基于掃描線算法（sweep-line algorithm）Voronoi庫可以實現如下功能：
v 輸入數據可以是點集和線段；
v 算法的穩定性高及輸出完整的拓樸信息；
v 可以控制輸出的幾何信息的精度；
計算幾何方面以穩定性著稱的CGAL庫中的Voronoi算法只滿足前兩個功能。S-Hull庫以上功能都沒有很好的滿足。下面是一些Boost.Polygon，CGAL，S-Hull庫的對比數據：
Figure 2.2 Construction time for 10 random points
Figure 2.3 Construction time for 100 random points
Figure 2.4 Construction time for 1000 random points
Figure 2.5 Construction time for 10000 random points
Figure 2.6 Memory usage for 100000 random points
Figure 2.7 Logarithmic Execution Time
結論：
v 在輸入上沒有限制這點上CGAL要優于Boost.Polygon；
v Boost.Polygon Voronoi的穩定性要高于S-Hull；
v Boost.Polygon Voronoi和S-Hull的時間復雜度為N*log(N)，而CGAL的不是；
v Boost.Polygon Voronoi的輸出頂點的精度高于CGAL庫；
v Boost.Polygon Voronoi的速度快；
v Boost.Polygon Voronoi根據10000個點或1000個線段來構造Voronoi的時間為0.02秒以內，所以可用來處理有實時性要求的場景；
3. Implementation
Boost.Polygon的Voronoi算法使用簡單，只需要輸入點集或線段集合，就可以直接構造出Voronoi圖了。最簡單的程序示例代碼如下：
/*
* Copyright (c) 2014 eryar All Rights Reserved.
*
* File : Main.cpp
* Author : eryar@163.com
* Date : 2014-05-06 18:28
* Version : V1.0
*
* Description : The Simplest example for boost voronoi library.
* Key words : boost voronoi, C++
*
*/
#include "boost/polygon/voronoi.hpp"
using namespace boost::polygon;
typedef int coordinate_type;
typedef point_data<coordinate_type> Point;
typedef voronoi_diagram<double> VD;
int main(int argc, char* argv[])
{
std::vector<Point> points;
points.push_back(Point(0, 0));
points.push_back(Point(1, 6));
points.push_back(Point(-4, 5));
points.push_back(Point(5, -1));
points.push_back(Point(3, -11));
points.push_back(Point(13, -1));
VD vd;
construct_voronoi(points.begin(), points.end(), &vd);
return 0;
}
且Boost.Polygon的Voronoi算法遍歷Voronoi邊Edges，Voronoi單元cell，Voronoi頂點Vertex也很直接。如下代碼所示為遍歷所有邊，來將剖分結果可視化：
/*
* Copyright (c) 2014 eryar All Rights Reserved.
*
* File : Main.cpp
* Author : eryar@163.com
* Date : 2014-05-06 18:28
* Version : V1.0
*
* Description : VoronoiViewer for boost voronoi library visulization.
* Key words : boost voronoi, C++, OpenSceneGraph
*
*/
#include <osgViewer/Viewer>
#include <osgGA/StateSetManipulator>
#include <osgViewer/ViewerEventHandlers>
#pragma comment(lib, "osgd.lib")
#pragma comment(lib, "osgDBd.lib")
#pragma comment(lib, "osgGAd.lib")
#pragma comment(lib, "osgViewerd.lib")
#include "boost/polygon/voronoi.hpp"
using namespace boost::polygon;
typedef double coordinate_type;
typedef point_data<coordinate_type> Point;
typedef voronoi_diagram<coordinate_type> VD;
osg::Node* BuildVoronoiDiagram(void)
{
srand(static_cast<unsigned int> (time(NULL)));
osg::ref_ptr<osg::Geode> theGeode = new osg::Geode();
osg::ref_ptr<osg::Geometry> theLines = new osg::Geometry();
osg::ref_ptr<osg::Vec3Array> theVertices = new osg::Vec3Array();
VD vd;
std::vector<Point> thePoints;
// Add points for the Voronoi Diagram.
for (int i = 0; i < 100; ++i)
{
int x = rand() % 100;
int y = rand() % 100;
thePoints.push_back(Point(x, y));
// Display the site of the Voronoi Diagram.
theVertices->push_back(osg::Vec3(x - 1, 0.0, y));
theVertices->push_back(osg::Vec3(x + 1, 0.0, y));
theVertices->push_back(osg::Vec3(x, 0.0, y - 1));
theVertices->push_back(osg::Vec3(x, 0.0, y + 1));
}
construct_voronoi(thePoints.begin(), thePoints.end(), &vd);
// Visualize the edge of the Voronoi Diagram.
// Traversing Voronoi edges using edge iterator.
for (VD::const_edge_iterator it = vd.edges().begin(); it != vd.edges().end(); ++it)
{
if (it->is_primary())
{
if (it->is_finite())
{
theVertices->push_back(osg::Vec3(it->vertex0()->x(), 0.0, it->vertex0()->y()));
theVertices->push_back(osg::Vec3(it->vertex1()->x(), 0.0, it->vertex1()->y()));
}
else
{
Point p1 = thePoints[it->cell()->source_index()];
Point p2 = thePoints[it->twin()->cell()->source_index()];
Point origin;
Point direction;
coordinate_type koef = 1.0;
origin.x((p1.x() + p2.x()) * 0.5);
origin.y((p1.y() + p2.y()) * 0.5);
direction.x(p1.y() - p2.y());
direction.y(p2.x() - p1.x());
if (it->vertex0() == NULL)
{
theVertices->push_back(osg::Vec3(
origin.x() - direction.x() * koef,
0.0,
origin.y() - direction.y() * koef));
}
else
{
theVertices->push_back(osg::Vec3(it->vertex0()->x(), 0.0, it->vertex0()->y()));
}
if (it->vertex1() == NULL)
{
theVertices->push_back(osg::Vec3(
origin.x() + direction.x() * koef,
0.0,
origin.y() + direction.y() * koef));
}
else
{
theVertices->push_back(osg::Vec3(it->vertex1()->x(), 0.0, it->vertex1()->y()));
}
}
}
}
theLines->setVertexArray(theVertices);
// Set the colors.
osg::ref_ptr<osg::Vec4Array> theColors = new osg::Vec4Array();
theColors->push_back(osg::Vec4(1.0f, 1.0f, 0.0f, 1.0f));
theLines->setColorArray(theColors);
theLines->setColorBinding(osg::Geometry::BIND_OVERALL);
// Set the normal.
osg::ref_ptr<osg::Vec3Array> theNormals = new osg::Vec3Array();
theNormals->push_back(osg::Vec3(0.0f, -1.0f, 0.0f));
theLines->setNormalArray(theNormals);
theLines->setNormalBinding(osg::Geometry::BIND_OVERALL);
theLines->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::LINES, 0, theVertices->size()));
theGeode->addDrawable(theLines);
return theGeode.release();
}
int main(int argc, char *argv[])
{
osgViewer::Viewer theViewer;
theViewer.setSceneData(BuildVoronoiDiagram());
theViewer.addEventHandler(new osgGA::StateSetManipulator(theViewer.getCamera()->getOrCreateStateSet()));
theViewer.addEventHandler(new osgViewer::StatsHandler);
theViewer.addEventHandler(new osgViewer::WindowSizeHandler);
return theViewer.run();
}
繪制Voronoi的邊時，當邊是有限的finite時，直接可以畫出直線；當邊是infinite時，根據定義計算出了無界邊的方向。顯示結果如下圖所示：
Figure 3.1 Construct Voronoi Diagram by 10 random points by Boost.Polygon
Figure 3.2 Construct Voronoi Diagram by 100 random points by Boost.Polygon
4. Conclusion
Boost.Polygon中的Voronoi圖算法穩定性及性能較高，且可以根據站點查找相關的拓樸信息，如根據站點查找Voronoi單元等；惟一不足的就是默認只處理整數點集。
當輸入有線段時，生成的Voronoi圖中有曲線，曲線的繪制可參考相關例子實現。
Boost.Polygon中的Voronoi算法都以模板實現，編譯時只需要包含相關的頭文件即可，不依賴其他庫，使用還是很方便的。
5. References
1. http://www.boost.org/
2. Boost.Polygon，http://www.boost.org/doc/libs/1_55_0/libs/polygon/doc/index.htm
3. Voronoi Basic Tutorial，\boost_1_54_0\libs\polygon\doc\voronoi_basic_tutorial.htm
4. 汪嘉業, 王文平, 屠長河, 楊承磊. 計算幾何及應用. 科學出版社. 2011
5. 楊承磊, 呂琳, 楊義軍, 孟祥旭. Voronoi圖及其應用. 清華大學出版社. 2013