Delaunay Triangulation in OpenCascade
eryar@163.com
摘要：本文簡要介紹了Delaunay三角剖分的基礎(chǔ)理論，并使用OpenCascade的三角剖分算法將邊界BRep表示的幾何體進(jìn)行三角離散化后在OpenSceneGraph中顯示。
關(guān)鍵字：Delaunay Triangulation、OpenCascade、OpenSceneGraph
一、 概述
三角剖分是平面剖分中的一個(gè)重要課題，在數(shù)字圖像處理、計(jì)算機(jī)三維曲面造型、有限元計(jì)算、逆向工程等領(lǐng)域有著廣泛應(yīng)用。由于三角形是平面域中的單純形，與其他平面圖形相比，其有描述方便、處理簡單等特性，很適合于對復(fù)雜區(qū)域進(jìn)行簡化處理。因此，無論在計(jì)算幾何、計(jì)算機(jī)圖形處理、模式識別、曲面逼近，還有有限元網(wǎng)格生成方面有廣泛的應(yīng)用。
雖然曲線、曲面等有精確的方程來表示，但是在在計(jì)算機(jī)中，只能用離散的方式來逼近。如曲線可用直線段來逼近，而曲面可用多邊形或三角形來表示。用多邊形網(wǎng)格表示曲面是設(shè)計(jì)中經(jīng)常使用的形式，可以根據(jù)應(yīng)用要求選擇網(wǎng)格的密度。利用三角形面片表示的曲面在計(jì)算機(jī)圖形學(xué)中也稱為三角形網(wǎng)格。用三角形網(wǎng)格表示曲面需要解決幾個(gè)問題：三角形的產(chǎn)生、描述、遍歷、簡化和壓縮等，這些問題都是計(jì)算幾何研究的范疇，相關(guān)問題都可以從中找到答案。下圖所示的圓柱和立方體是由OpenCascade生成，使用OpenCascade的算法離散成三角網(wǎng)格后在OpenSceneGraph中顯示的效果。
Figure 1.1 Shaded Cylinder and Box
Figure 1.2 Mesh generated by OpenCascade
從圖中可以看出，平面的三角形網(wǎng)格效果還不錯(cuò)，曲面的三角形網(wǎng)格表示只能是近似表示，可以通過提高網(wǎng)格的密度來增加真實(shí)性，但相應(yīng)渲染的數(shù)據(jù)量就大了。有人說OpenCascade的顯示模塊做得不是很好，上述方法則可以只使用OpenCascade的造型模塊，再結(jié)合OpenSceneGraph來對圖形進(jìn)行顯示。
三維數(shù)據(jù)交換STL格式文件中保存的都是三角面片的數(shù)據(jù)，STL文件格式是由美國3D System公司開發(fā)，已被工業(yè)界認(rèn)為是目前快速自動(dòng)成型領(lǐng)域的準(zhǔn)標(biāo)準(zhǔn)零件描述文件格式。它對三維實(shí)體描述的解釋具有惟一性。幾乎所有的幾何造型系統(tǒng)都提供STL文件數(shù)據(jù)交換接口。OpenCascade中的數(shù)據(jù)交換模塊也提供對STL格式的支持，由此可見三角網(wǎng)格在幾何造型系統(tǒng)中的重要性。
Voronoi圖和Delaunay三角剖分的應(yīng)用領(lǐng)域十分廣泛：幾何建模——用來尋找三維曲面“好的”三角剖分；有限元分析——用來生成“好的”有限元網(wǎng)格；地理信息系統(tǒng)——用來進(jìn)行空間領(lǐng)域分析；結(jié)晶學(xué)——用來確定合金的結(jié)構(gòu)；人類學(xué)和考古學(xué)——用來確定氏族部落、首領(lǐng)權(quán)威、居住中心或堡壘等的影響范圍；天文學(xué)——用來確定恒星和星系的分布；生物學(xué)生態(tài)學(xué)和林學(xué)——用來確定動(dòng)植物的競爭；動(dòng)物學(xué)——分析動(dòng)物的領(lǐng)地；統(tǒng)計(jì)學(xué)和數(shù)據(jù)分析——用來分析統(tǒng)計(jì)聚合；機(jī)器人學(xué)——用來進(jìn)行運(yùn)動(dòng)軌跡規(guī)劃（在存在障礙物的情況下）；模式識別——作為尋找物體骨架點(diǎn)的工具；生理學(xué)——用來分析毛細(xì)作用的領(lǐng)域；氣象學(xué)——用來估計(jì)區(qū)域平均降雨量；市場學(xué)——用來建立城市的市場輻射范圍；以及在遙感圖像處理、化學(xué)、地理學(xué)、地質(zhì)學(xué)、冶金學(xué)、數(shù)學(xué)等學(xué)科的應(yīng)用等。
本文只對OpenCascade中的三角剖分進(jìn)行簡要介紹，希望對三角剖分在三維幾何造型方面有興趣的朋友可以對其深入研究。水平很有限，文中不當(dāng)之處歡迎批評指正、指導(dǎo)，聯(lián)系郵箱：eryar@163.com。
二、 Voronoi圖
Dirichlet于1850年研究了平面點(diǎn)的鄰域問題，Voronoi于1908年將其結(jié)果擴(kuò)展到高維空間。半空間定義Voronoi圖：給定平面上n個(gè)點(diǎn)集S，S={p1, p2, …, pn}。定義：
PiPj連線的垂直平分面將空間分為兩半，V(Pi)表示比其他點(diǎn)更接近Pi的點(diǎn)的軌跡是n-1個(gè)半平面的交，它是一個(gè)不多于n-1條邊的凸多邊形域，稱為關(guān)聯(lián)于Pi的Voronoi多邊形或關(guān)聯(lián)于Pi的Voronoi域。如下圖所示為關(guān)聯(lián)于P1的Voronoi多邊形，它是一個(gè)四邊形，而n=6.
Figure 2.1 n=6時(shí)的一種V(p1)
對于點(diǎn)集S中的每個(gè)點(diǎn)都可以做一個(gè)Voronoi多邊形，這樣的n個(gè)Voronoi多邊形組成的圖稱為Voronoi圖，記為Vor(S)。如下圖所示：
Figure 2.2 Voronoi diagram for 10 randomly points (Generated by MATLAB)
圖中的頂點(diǎn)和邊分別稱為Voronoi頂點(diǎn)和Voronoi邊。顯然，|S|=n時(shí)，Vor(S)劃分平面成n個(gè)多邊形域，每個(gè)多邊形域V(Pi)包含S中的一個(gè)點(diǎn)而且只包含S中的一個(gè)點(diǎn)，Vor(S)的邊是S中某點(diǎn)對的垂直平分線上的一條線段或半直線，從而為該點(diǎn)對所在的兩個(gè)多邊形域所共有。Vor(S)中有的多邊形域是無界的。
Figure 2.3 Ten shops in a flat city and their Voronoi cells
(http://en.wikipedia.org/wiki/Voronoi_diagram)
Figure 2.4 Voronoi tessellation in a cylinder (Voro++ library: http://math.lbl.gov/voro++/)
Voronoi圖有如下性質(zhì)：
l n個(gè)點(diǎn)的點(diǎn)集S的Voronoi圖至多有2n-5個(gè)頂點(diǎn)和3n-6條邊；
l 每個(gè)Voronoi點(diǎn)恰好是三條Voronoi邊的交點(diǎn)；
l 設(shè)v是Vor(S)的頂點(diǎn)，則圓C(v)內(nèi)不含S的其他點(diǎn)；
l 點(diǎn)集S中點(diǎn)Pi的每一個(gè)最近鄰近點(diǎn)確定V(Pi)的一條邊；
l Voronoi圖的直線對偶圖是S的一個(gè)三角剖分；
l 如果Pi,Pj屬于S，并且通過Pi,Pj有一個(gè)不包含S中其他點(diǎn)的圓，那么線段PiPj是點(diǎn)集S三角剖分的一條邊，反之亦成立。
三、 Delaunay三角剖分
1. 二維實(shí)數(shù)域上的三角剖分
假設(shè)V是二維實(shí)數(shù)域上的有限點(diǎn)集，邊e是由點(diǎn)集中的點(diǎn)作端點(diǎn)構(gòu)成的封閉線段，E為e的集合，那么該點(diǎn)集V的一個(gè)三角剖分T＝（V，E）是一個(gè)平面圖：
l 除了端點(diǎn)，平面圖中的邊不包含點(diǎn)集中的任何點(diǎn)；
l 沒有相交邊；
l 平面圖中所有的面都是三角面，且所有三角面的合集是點(diǎn)集V的凸包。
2. Delaunay邊
假設(shè)E中的一條邊（兩端點(diǎn)a,b），e滿足下列條件，則稱為Delaunay邊：存在一個(gè)圓經(jīng)過a，b兩點(diǎn)，圓內(nèi)不包含點(diǎn)集V中的任何的點(diǎn)。這一特性又稱為空圓特性。
3. Delaunay三角剖分
如果點(diǎn)集V的一個(gè)三角剖分T中只包含Delaunay邊，那么該三角剖分稱為Delaunay剖分。
最近點(diǎn)意義下的Voronoi圖的對偶圖實(shí)際上是點(diǎn)集的一種三角剖分，該三角剖分就是Delaunay剖分（表示為DT(S)），其中每個(gè)三角形的外接圓不包含點(diǎn)集中的其他任何點(diǎn)。因此，在構(gòu)造點(diǎn)集的Voronoi圖之后，再作其對偶圖，即對每條Voronoi邊作通過點(diǎn)集中某兩點(diǎn)的垂直平分線，即得到Delaunay三角剖分。
Figure 3.1 Delaunay Triangulation (Generated by MATLAB)
再看幾個(gè)圖片，加深對Delaunay三角剖分的理解：
Figure 3.2 Delaunay Edge
Figure 3.3 Illustrate Delaunay Edge
Figure 3.4 Delaunay Edge
4. Delaunay三角剖分的特性
l 1978年Sibson證明了在二維的情況下，在點(diǎn)集的所有三角剖分中，Delaunay三角剖分使得生成的三角形的最小角達(dá)到最大（max-min angle）。因?yàn)檫@一特性，對于給定點(diǎn)集的Delaunay三角剖分總是盡可能避免“瘦長”三角形，自動(dòng)向等邊三角形逼近；
l 局部優(yōu)化與整體優(yōu)化（locally optimal and globally optimal）；
l Delaunay空洞（cavity）與局部重連（local reconnection）；
5. 經(jīng)典的Delaunay三角剖分算法
目前常用的算法分為幾種，有掃描線法（Sweepline）、隨機(jī)增量法（Incremental）、分治法（Divide and Conquer）等。
經(jīng)典的Delaunay三角剖分算法主要有兩類：Bowyer/Watson算法和局部變換法。
l Bowyer/Watson算法又稱為Delaunay空洞算法或加點(diǎn)法，以Bowyer和Watson算法為代表。從一個(gè)三角形開始，每次加一個(gè)點(diǎn)，保證每一步得到的當(dāng)前三角形是局部優(yōu)化的。以英國Bath大學(xué)數(shù)學(xué)分校Bowyer，Green，Sibson為代表的計(jì)算Dirichlet圖的方法屬于加點(diǎn)法，是較早成名的算法之一；以澳大利亞悉尼大學(xué)地學(xué)系Watson為代表的空外接球法也屬于加點(diǎn)法。加點(diǎn)法算法簡明，是目前應(yīng)用最多的算法，該方法利用了Delaunay空洞性質(zhì)。Bowyer/Watson算法的優(yōu)點(diǎn)是與空間的維數(shù)無關(guān)，并且算法在實(shí)現(xiàn)上比局部變換算法簡單。該算法在新點(diǎn)加入到Delaunay網(wǎng)格時(shí)，部分外接球包含新點(diǎn)的三角形單元不再符合Delaunay屬性，則這些三角形單元被刪除，形成Delaunay空洞，然后算法將新點(diǎn)與組成空洞的每一個(gè)頂點(diǎn)相連生成一個(gè)新邊，根據(jù)空球?qū)傩钥梢宰C明這些新邊都是局部Delaunay的，因此新生成的三角網(wǎng)格仍是Delaunay的。
Figure 3.5 Illustration of 2D Bowyer/Watson algorithm for Delaunay Triangulation
l 局部變換法又稱為換邊、換面法。當(dāng)利用局部變換法實(shí)現(xiàn)增量式點(diǎn)集的Delaunay三角剖分時(shí)，首先定位新加入點(diǎn)所在的三角形，然后在網(wǎng)格中加入三個(gè)新的連接該三角形頂點(diǎn)與新頂點(diǎn)的邊，若該新點(diǎn)位于某條邊上，則該邊被刪除，四條連接該新點(diǎn)的邊被加入。最后，在通過換邊方法對該新點(diǎn)的局部區(qū)域內(nèi)的邊進(jìn)行檢測和變換，重新維護(hù)網(wǎng)格的Delaunay性質(zhì)。局部變換法的另一個(gè)優(yōu)點(diǎn)是其可以對已存在的三角網(wǎng)格進(jìn)行優(yōu)化，使其變換成為Delaunay三角網(wǎng)格，該方法的缺點(diǎn)則是當(dāng)算法擴(kuò)展到高維空間時(shí)變得較為復(fù)雜。
四、 Delaunay三角剖分在OpenCascade的應(yīng)用
OpenCascade中網(wǎng)格剖分的包主要有BRepMesh、MeshAlgo、MeshVS，其中，類MeshAlgo_Delaunay使用算法Watson來進(jìn)行Delaunay三角剖分。從類StlTransfer中的注釋The triangulation is computed with the Delaunay algorithm implemented in package BRepMesh.可以看出包BRepMesh就是Delaunay三角剖分的具體實(shí)現(xiàn)。使用方法如下：
BRepMesh::Mesh (aShape, Deflection);
這個(gè)函數(shù)主要是用來對拓?fù)湫螤钸M(jìn)行三角剖分。以下通過將一個(gè)圓柱三角剖分為例說明如何將一個(gè)拓?fù)湫螤钸M(jìn)行三角剖分并將結(jié)果進(jìn)行可視化。
/**//**
* Copyright (c) 2013 eryar All Rights Reserved.
*
* File : Main.cpp
* Author : eryar@163.com
* Date : 2013-05-26
* Version : 0.1
*
* Description : Use BRepMesh_Delaun class to learn
* Delaunay's triangulation algorithm.
*
*/
// Open Cascade library.
#include <gp_Pnt.hxx>
#include <gp_Pln.hxx>
#include <BRep_Tool.hxx>
#include <TopoDS.hxx>
#include <TopoDS_Edge.hxx>
#include <TopoDS_Wire.hxx>
#include <TopoDS_Face.hxx>
#include <BRepBuilderAPI_MakeEdge.hxx>
#include <BRepBuilderAPI_MakeWire.hxx>
#include <BRepBuilderAPI_MakeFace.hxx>
#include <BRepPrimAPI_MakeBox.hxx>
#include <BRepPrimAPI_MakeCone.hxx>
#include <BRepPrimAPI_MakeCylinder.hxx>
#include <BRepPrimApI_MakeSphere.hxx>
#include <BRepMesh.hxx>
#include <TopExp_Explorer.hxx>
#include <Poly_Triangulation.hxx>
#include <TShort_Array1OfShortReal.hxx>
#pragma comment(lib, "TKernel.lib")
#pragma comment(lib, "TKMath.lib")
#pragma comment(lib, "TKBRep.lib")
#pragma comment(lib, "TKPrim.lib")
#pragma comment(lib, "TKMesh.lib")
#pragma comment(lib, "TKTopAlgo.lib")
// OpenSceneGraph library.
#include <osgDB/ReadFile>
#include <osgViewer/Viewer>
#include <osgViewer/ViewerEventHandlers>
#include <osgGA/StateSetManipulator>
#pragma comment(lib, "osgd.lib")
#pragma comment(lib, "osgDbd.lib")
#pragma comment(lib, "osgGAd.lib")
#pragma comment(lib, "osgViewerd.lib")
osg::Node* BuildShapeMesh(const TopoDS_Shape& aShape)
{
osg::ref_ptr<osg::Group> root = new osg::Group();
osg::ref_ptr<osg::Geode> geode = new osg::Geode();
osg::ref_ptr<osg::Geometry> triGeom = new osg::Geometry();
osg::ref_ptr<osg::Vec3Array> vertices = new osg::Vec3Array();
osg::ref_ptr<osg::Vec3Array> normals = new osg::Vec3Array();
BRepMesh::Mesh(aShape, 1);
TopExp_Explorer faceExplorer;
for (faceExplorer.Init(aShape, TopAbs_FACE); faceExplorer.More(); faceExplorer.Next())
{
TopLoc_Location loc;
TopoDS_Face aFace = TopoDS::Face(faceExplorer.Current());
Handle_Poly_Triangulation triFace = BRep_Tool::Triangulation(aFace, loc);
Standard_Integer nTriangles = triFace->NbTriangles();
gp_Pnt vertex1;
gp_Pnt vertex2;
gp_Pnt vertex3;
Standard_Integer nVertexIndex1 = 0;
Standard_Integer nVertexIndex2 = 0;
Standard_Integer nVertexIndex3 = 0;
TColgp_Array1OfPnt nodes(1, triFace->NbNodes());
Poly_Array1OfTriangle triangles(1, triFace->NbTriangles());
nodes = triFace->Nodes();
triangles = triFace->Triangles();
for (Standard_Integer i = 1; i <= nTriangles; i++)
{
Poly_Triangle aTriangle = triangles.Value(i);
aTriangle.Get(nVertexIndex1, nVertexIndex2, nVertexIndex3);
vertex1 = nodes.Value(nVertexIndex1);
vertex2 = nodes.Value(nVertexIndex2);
vertex3 = nodes.Value(nVertexIndex3);
gp_XYZ vector12(vertex2.XYZ() - vertex1.XYZ());
gp_XYZ vector13(vertex3.XYZ() - vertex1.XYZ());
gp_XYZ normal = vector12.Crossed(vector13);
Standard_Real rModulus = normal.Modulus();
if (rModulus > gp::Resolution())
{
normal.Normalize();
}
else
{
normal.SetCoord(0., 0., 0.);
}
vertices->push_back(osg::Vec3(vertex1.X(), vertex1.Y(), vertex1.Z()));
vertices->push_back(osg::Vec3(vertex2.X(), vertex2.Y(), vertex2.Z()));
vertices->push_back(osg::Vec3(vertex3.X(), vertex3.Y(), vertex3.Z()));
normals->push_back(osg::Vec3(normal.X(), normal.Y(), normal.Z()));
}
}
triGeom->setVertexArray(vertices.get());
triGeom->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::TRIANGLES, 0, vertices->size()));
triGeom->setNormalArray(normals);
triGeom->setNormalBinding(osg::Geometry::BIND_PER_PRIMITIVE);
geode->addDrawable(triGeom);
root->addChild(geode);
return root.release();
}
int main(int argc, char* argv[])
{
osgViewer::Viewer myViewer;
osg::ref_ptr<osg::Group> root = new osg::Group();
root->addChild(BuildShapeMesh(BRepPrimAPI_MakeCylinder(.6, 1)));
myViewer.setSceneData(root);
myViewer.addEventHandler(new osgGA::StateSetManipulator(myViewer.getCamera()->getOrCreateStateSet()));
myViewer.addEventHandler(new osgViewer::StatsHandler);
myViewer.addEventHandler(new osgViewer::WindowSizeHandler);
return myViewer.run();
}
結(jié)果如下圖所示：
Figure 4.1 Cylinder mesh generated by BRepMesh::Mesh
BRepMesh::Mesh是經(jīng)過封裝的，便于對拓?fù)湫螤钸M(jìn)行三角剖分。以下通過一個(gè)簡單的例子來說明直接使用BRepMesh_Delaun的方法：
/**
* Copyright (c) 2013 eryar All Rights Reserved.
*
* File : Main.cpp
* Author : eryar@163.com
* Date : 2013-05-26
* Version : 0.1
*
* Description : Use BRepMesh_Delaun class to learn
* Delaunay's triangulation algorithm.
*
*/
#include <BRepMesh_Edge.hxx>
#include <BRepMesh_Delaun.hxx>
#include <BRepMesh_Array1OfVertexOfDelaun.hxx>
#include <TColStd_MapIteratorOfMapOfInteger.hxx>
#pragma comment(lib, "TKernel.lib")
#pragma comment(lib, "TKMesh.lib")
int main(int argc, char* argv[])
{
BRepMesh_Array1OfVertexOfDelaun vertices(1, 4);
vertices.SetValue(1, BRepMesh_Vertex(0, 0, MeshDS_Free));
vertices.SetValue(2, BRepMesh_Vertex(1, 0, MeshDS_Free));
vertices.SetValue(3, BRepMesh_Vertex(1, 1, MeshDS_Free));
vertices.SetValue(4, BRepMesh_Vertex(0, 1, MeshDS_Free));
BRepMesh_Delaun triangulation(vertices);
//triangulation.AddVertex(BRepMesh_Vertex(0.5, 0.5, MeshDS_OnSurface));
Handle_BRepMesh_DataStructureOfDelaun meshData = triangulation.Result();
std::cout<<"Iterate Mesh Triangles:"<<std::endl;
MeshDS_MapOfInteger::Iterator triDom;
for (triDom.Initialize(meshData->ElemOfDomain()); triDom.More(); triDom.Next())
{
Standard_Integer triId = triDom.Key();
const BRepMesh_Triangle& curTri = meshData->GetElement(triId);
Standard_Integer vertexIndex1 = 0;
Standard_Integer vertexIndex2 = 0;
Standard_Integer vertexIndex3 = 0;
Standard_Integer edgeIndex1 = 0;
Standard_Integer edgeIndex2 = 0;
Standard_Integer edgeIndex3 = 0;
Standard_Boolean o1 = Standard_False;
Standard_Boolean o2 = Standard_False;
Standard_Boolean o3 = Standard_False;
curTri.Edges(edgeIndex1, edgeIndex2, edgeIndex3, o1, o2, o3);
const BRepMesh_Edge& edge1 = meshData->GetLink(edgeIndex1);
const BRepMesh_Edge& edge2 = meshData->GetLink(edgeIndex2);
const BRepMesh_Edge& edge3 = meshData->GetLink(edgeIndex3);
vertexIndex1 = (o1? edge1.FirstNode(): edge1.LastNode());
vertexIndex2 = (o1? edge1.LastNode() : edge1.FirstNode());
vertexIndex3 = (o2? edge2.LastNode() : edge2.FirstNode());
const BRepMesh_Vertex& vertex1 = meshData->GetNode(vertexIndex1);
const BRepMesh_Vertex& vertex2 = meshData->GetNode(vertexIndex2);
const BRepMesh_Vertex& vertex3 = meshData->GetNode(vertexIndex3);
const gp_XY& p1 = vertex1.Coord();
const gp_XY& p2 = vertex2.Coord();
const gp_XY& p3 = vertex3.Coord();
std::cout<<"--------"<<std::endl;
std::cout<<p1.X()<<" , "<<p1.Y()<<std::endl;
std::cout<<p2.X()<<" , "<<p2.Y()<<std::endl;
std::cout<<p3.X()<<" , "<<p3.Y()<<std::endl;
std::cout<<"========"<<std::endl;
}
return 0;
}
上述程序是以一個(gè)正方形為例，使用BRepMesh_Delaun三角剖分的結(jié)果為兩個(gè)三角形，如下所示：
Iterate Mesh Triangles:
--------
1 , 1
0 , 0
1 , 0
========
--------
1 , 1
0 , 1
0 , 0
========
以上結(jié)果都是二維空間上的，三維空間中的使用方法可以參考類：BRepMesh_FastDiscretFace。這個(gè)類說明了如何將一個(gè)面進(jìn)行網(wǎng)格劃分。
五、 結(jié)論
Delaunay三角剖分理論在三維幾何造型中還是比較重要的，通過對形狀的三角剖分，不僅可以對其進(jìn)行可視化，還便于對形狀做進(jìn)一步的處理，如消隱、光照處理等。通過對OpenCascade中三角剖分算法的使用，以進(jìn)一步了解三角剖分理論應(yīng)用及其算法實(shí)現(xiàn)。
六、 參考資料
1. 周培德. 計(jì)算幾何—算法設(shè)計(jì)與分析. 清華大學(xué)出版社, 2011
2. 李海生. Delaunay三角剖分理論及可視化應(yīng)用研究. 哈爾濱工業(yè)大學(xué)出版社, 2010
3. 何援軍. 計(jì)算機(jī)圖形學(xué). 機(jī)械工業(yè)出版社, 2010
4. 周元峰, 孫峰, 王文平, 汪嘉業(yè), 張彩明. 基于局部修復(fù)的移動(dòng)數(shù)據(jù)點(diǎn)Delaunay三角化快速更新方法. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 2011, 12: 2006-1012
5. http://en.wikipedia.org/wiki/Voronoi_diagram