Delaunay Triangulation in OpenCascade
eryar@163.com
摘要：本文簡(jiǎn)要介紹了Delaunay三角剖分的基礎(chǔ)理論，并使用OpenCascade的三角剖分算法將邊界BRep表示的幾何體進(jìn)行三角離散化后在OpenSceneGraph中顯示。
關(guān)鍵字：Delaunay Triangulation、OpenCascade、OpenSceneGraph
一、 概述
三角剖分是平面剖分中的一個(gè)重要課題，在數(shù)字圖像處理、計(jì)算機(jī)三維曲面造型、有限元計(jì)算、逆向工程等領(lǐng)域有著廣泛應(yīng)用。由于三角形是平面域中的單純形，與其他平面圖形相比，其有描述方便、處理簡(jiǎn)單等特性，很適合于對(duì)復(fù)雜區(qū)域進(jìn)行簡(jiǎn)化處理。因此，無論在計(jì)算幾何、計(jì)算機(jī)圖形處理、模式識(shí)別、曲面逼近，還有有限元網(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)生、描述、遍歷、簡(jiǎ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來對(duì)圖形進(jìn)行顯示。
三維數(shù)據(jù)交換STL格式文件中保存的都是三角面片的數(shù)據(jù)，STL文件格式是由美國3D System公司開發(fā)，已被工業(yè)界認(rèn)為是目前快速自動(dòng)成型領(lǐng)域的準(zhǔn)標(biāo)準(zhǔn)零件描述文件格式。它對(duì)三維實(shí)體描述的解釋具有惟一性。幾乎所有的幾何造型系統(tǒng)都提供STL文件數(shù)據(jù)交換接口。OpenCascade中的數(shù)據(jù)交換模塊也提供對(duì)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)植物的競(jìng)爭(zhē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ī)劃（在存在障礙物的情況下）；模式識(shí)別——作為尋找物體骨架點(diǎn)的工具；生理學(xué)——用來分析毛細(xì)作用的領(lǐng)域；氣象學(xué)——用來估計(jì)區(qū)域平均降雨量；市場(chǎng)學(xué)——用來建立城市的市場(chǎng)輻射范圍；以及在遙感圖像處理、化學(xué)、地理學(xué)、地質(zhì)學(xué)、冶金學(xué)、數(shù)學(xué)等學(xué)科的應(yīng)用等。
本文只對(duì)OpenCascade中的三角剖分進(jìn)行簡(jiǎn)要介紹，希望對(duì)三角剖分在三維幾何造型方面有興趣的朋友可以對(duì)其深入研究。水平很有限，文中不當(dāng)之處歡迎批評(pí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)
對(duì)于點(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)對(duì)的垂直平分線上的一條線段或半直線，從而為該點(diǎn)對(duì)所在的兩個(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圖的直線對(duì)偶圖是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)。這一特性又稱為空?qǐng)A特性。
3. Delaunay三角剖分
如果點(diǎn)集V的一個(gè)三角剖分T中只包含Delaunay邊，那么該三角剖分稱為Delaunay剖分。
最近點(diǎn)意義下的Voronoi圖的對(duì)偶圖實(shí)際上是點(diǎn)集的一種三角剖分，該三角剖分就是Delaunay剖分（表示為DT(S)），其中每個(gè)三角形的外接圓不包含點(diǎn)集中的其他任何點(diǎn)。因此，在構(gòu)造點(diǎn)集的Voronoi圖之后，再作其對(duì)偶圖，即對(duì)每條Voronoi邊作通過點(diǎn)集中某兩點(diǎn)的垂直平分線，即得到Delaunay三角剖分。
Figure 3.1 Delaunay Triangulation (Generated by MATLAB)
再看幾個(gè)圖片，加深對(duì)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)檫@一特性，對(duì)于給定點(diǎn)集的Delaunay三角剖分總是盡可能避免“瘦長(zhǎng)”三角形，自動(dòng)向等邊三角形逼近；
l 局部?jī)?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)前三角形是局部?jī)?yōu)化的。以英國Bath大學(xué)數(shù)學(xué)分校Bowyer，Green，Sibson為代表的計(jì)算Dirichlet圖的方法屬于加點(diǎn)法，是較早成名的算法之一；以澳大利亞悉尼大學(xué)地學(xué)系Watson為代表的空外接球法也屬于加點(diǎn)法。加點(diǎn)法算法簡(jiǎn)明，是目前應(yīng)用最多的算法，該方法利用了Delaunay空洞性質(zhì)。Bowyer/Watson算法的優(yōu)點(diǎn)是與空間的維數(shù)無關(guān)，并且算法在實(shí)現(xiàn)上比局部變換算法簡(jiǎ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)的邊被加入。最后，在通過換邊方法對(duì)該新點(diǎn)的局部區(qū)域內(nèi)的邊進(jìn)行檢測(cè)和變換，重新維護(hù)網(wǎng)格的Delaunay性質(zhì)。局部變換法的另一個(gè)優(yōu)點(diǎn)是其可以對(duì)已存在的三角網(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ù)主要是用來對(duì)拓?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)過封裝的，便于對(duì)拓?fù)湫螤钸M(jìn)行三角剖分。以下通過一個(gè)簡(jiǎn)單的例子來說明直接使用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三角剖分理論在三維幾何造型中還是比較重要的，通過對(duì)形狀的三角剖分，不僅可以對(duì)其進(jìn)行可視化，還便于對(duì)形狀做進(jìn)一步的處理，如消隱、光照處理等。通過對(duì)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