Delaunay Triangulation in OpenCascade

eryar@163.com

摘要：本文簡要介紹了Delaunay三角剖分的基礎理論，并使用OpenCascade的三角剖分算法將邊界BRep表示的幾何體進行三角離散化后在OpenSceneGraph中顯示。

關鍵字：Delaunay Triangulation、OpenCascade、OpenSceneGraph

一、 概述

三角剖分是平面剖分中的一個重要課題，在數字圖像處理、計算機三維曲面造型、有限元計算、逆向工程等領域有著廣泛應用。由于三角形是平面域中的單純形，與其他平面圖形相比，其有描述方便、處理簡單等特性，很適合于對復雜區域進行簡化處理。因此，無論在計算幾何、計算機圖形處理、模式識別、曲面逼近，還有有限元網格生成方面有廣泛的應用。

雖然曲線、曲面等有精確的方程來表示，但是在在計算機中，只能用離散的方式來逼近。如曲線可用直線段來逼近，而曲面可用多邊形或三角形來表示。用多邊形網格表示曲面是設計中經常使用的形式，可以根據應用要求選擇網格的密度。利用三角形面片表示的曲面在計算機圖形學中也稱為三角形網格。用三角形網格表示曲面需要解決幾個問題：三角形的產生、描述、遍歷、簡化和壓縮等，這些問題都是計算幾何研究的范疇，相關問題都可以從中找到答案。下圖所示的圓柱和立方體是由OpenCascade生成，使用OpenCascade的算法離散成三角網格后在OpenSceneGraph中顯示的效果。

Figure 1.1 Shaded Cylinder and Box

Figure 1.2 Mesh generated by OpenCascade

從圖中可以看出，平面的三角形網格效果還不錯，曲面的三角形網格表示只能是近似表示，可以通過提高網格的密度來增加真實性，但相應渲染的數據量就大了。有人說OpenCascade的顯示模塊做得不是很好，上述方法則可以只使用OpenCascade的造型模塊，再結合OpenSceneGraph來對圖形進行顯示。

三維數據交換STL格式文件中保存的都是三角面片的數據，STL文件格式是由美國3D System公司開發，已被工業界認為是目前快速自動成型領域的準標準零件描述文件格式。它對三維實體描述的解釋具有惟一性。幾乎所有的幾何造型系統都提供STL文件數據交換接口。OpenCascade中的數據交換模塊也提供對STL格式的支持，由此可見三角網格在幾何造型系統中的重要性。

Voronoi圖和Delaunay三角剖分的應用領域十分廣泛：幾何建模——用來尋找三維曲面“好的”三角剖分；有限元分析——用來生成“好的”有限元網格；地理信息系統——用來進行空間領域分析；結晶學——用來確定合金的結構；人類學和考古學——用來確定氏族部落、首領權威、居住中心或堡壘等的影響范圍；天文學——用來確定恒星和星系的分布；生物學生態學和林學——用來確定動植物的競爭；動物學——分析動物的領地；統計學和數據分析——用來分析統計聚合；機器人學——用來進行運動軌跡規劃（在存在障礙物的情況下）；模式識別——作為尋找物體骨架點的工具；生理學——用來分析毛細作用的領域；氣象學——用來估計區域平均降雨量；市場學——用來建立城市的市場輻射范圍；以及在遙感圖像處理、化學、地理學、地質學、冶金學、數學等學科的應用等。

本文只對OpenCascade中的三角剖分進行簡要介紹，希望對三角剖分在三維幾何造型方面有興趣的朋友可以對其深入研究。水平很有限，文中不當之處歡迎批評指正、指導，聯系郵箱：eryar@163.com。

二、 Voronoi圖

Dirichlet于1850年研究了平面點的鄰域問題，Voronoi于1908年將其結果擴展到高維空間。半空間定義Voronoi圖：給定平面上n個點集S，S={p1, p2, …, pn}。定義：

PiPj連線的垂直平分面將空間分為兩半，V(Pi)表示比其他點更接近Pi的點的軌跡是n-1個半平面的交，它是一個不多于n-1條邊的凸多邊形域，稱為關聯于Pi的Voronoi多邊形或關聯于Pi的Voronoi域。如下圖所示為關聯于P1的Voronoi多邊形，它是一個四邊形，而n=6.

Figure 2.1 n=6時的一種V(p1)

對于點集S中的每個點都可以做一個Voronoi多邊形，這樣的n個Voronoi多邊形組成的圖稱為Voronoi圖，記為Vor(S)。如下圖所示：

Figure 2.2 Voronoi diagram for 10 randomly points (Generated by MATLAB)

圖中的頂點和邊分別稱為Voronoi頂點和Voronoi邊。顯然，|S|=n時，Vor(S)劃分平面成n個多邊形域，每個多邊形域V(Pi)包含S中的一個點而且只包含S中的一個點，Vor(S)的邊是S中某點對的垂直平分線上的一條線段或半直線，從而為該點對所在的兩個多邊形域所共有。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圖有如下性質：

l n個點的點集S的Voronoi圖至多有2n-5個頂點和3n-6條邊；

l 每個Voronoi點恰好是三條Voronoi邊的交點；

l 設v是Vor(S)的頂點，則圓C(v)內不含S的其他點；

l 點集S中點Pi的每一個最近鄰近點確定V(Pi)的一條邊；

l Voronoi圖的直線對偶圖是S的一個三角剖分；

l 如果Pi,Pj屬于S，并且通過Pi,Pj有一個不包含S中其他點的圓，那么線段PiPj是點集S三角剖分的一條邊，反之亦成立。

三、 Delaunay三角剖分 

1. 二維實數域上的三角剖分

假設V是二維實數域上的有限點集，邊e是由點集中的點作端點構成的封閉線段，E為e的集合，那么該點集V的一個三角剖分T＝（V，E）是一個平面圖：

l 除了端點，平面圖中的邊不包含點集中的任何點；

l 沒有相交邊；

l 平面圖中所有的面都是三角面，且所有三角面的合集是點集V的凸包。

2. Delaunay邊

假設E中的一條邊（兩端點a,b），e滿足下列條件，則稱為Delaunay邊：存在一個圓經過a，b兩點，圓內不包含點集V中的任何的點。這一特性又稱為空圓特性。

3. Delaunay三角剖分

如果點集V的一個三角剖分T中只包含Delaunay邊，那么該三角剖分稱為Delaunay剖分。

最近點意義下的Voronoi圖的對偶圖實際上是點集的一種三角剖分，該三角剖分就是Delaunay剖分（表示為DT(S)），其中每個三角形的外接圓不包含點集中的其他任何點。因此，在構造點集的Voronoi圖之后，再作其對偶圖，即對每條Voronoi邊作通過點集中某兩點的垂直平分線，即得到Delaunay三角剖分。

Figure 3.1 Delaunay Triangulation (Generated by MATLAB)

再看幾個圖片，加深對Delaunay三角剖分的理解：

Figure 3.2 Delaunay Edge 

Figure 3.3 Illustrate Delaunay Edge

Figure 3.4 Delaunay Edge

4. Delaunay三角剖分的特性

l 1978年Sibson證明了在二維的情況下，在點集的所有三角剖分中，Delaunay三角剖分使得生成的三角形的最小角達到最大（max-min angle）。因為這一特性，對于給定點集的Delaunay三角剖分總是盡可能避免“瘦長”三角形，自動向等邊三角形逼近；

l 局部優化與整體優化（locally optimal and globally optimal）；

l Delaunay空洞（cavity）與局部重連（local reconnection）；

5. 經典的Delaunay三角剖分算法

目前常用的算法分為幾種，有掃描線法（Sweepline）、隨機增量法（Incremental）、分治法（Divide and Conquer）等。

經典的Delaunay三角剖分算法主要有兩類：Bowyer/Watson算法和局部變換法。

l Bowyer/Watson算法又稱為Delaunay空洞算法或加點法，以Bowyer和Watson算法為代表。從一個三角形開始，每次加一個點，保證每一步得到的當前三角形是局部優化的。以英國Bath大學數學分校Bowyer，Green，Sibson為代表的計算Dirichlet圖的方法屬于加點法，是較早成名的算法之一；以澳大利亞悉尼大學地學系Watson為代表的空外接球法也屬于加點法。加點法算法簡明，是目前應用最多的算法，該方法利用了Delaunay空洞性質。Bowyer/Watson算法的優點是與空間的維數無關，并且算法在實現上比局部變換算法簡單。該算法在新點加入到Delaunay網格時，部分外接球包含新點的三角形單元不再符合Delaunay屬性，則這些三角形單元被刪除，形成Delaunay空洞，然后算法將新點與組成空洞的每一個頂點相連生成一個新邊，根據空球屬性可以證明這些新邊都是局部Delaunay的，因此新生成的三角網格仍是Delaunay的。

Figure 3.5 Illustration of 2D Bowyer/Watson algorithm for Delaunay Triangulation

l 局部變換法又稱為換邊、換面法。當利用局部變換法實現增量式點集的Delaunay三角剖分時，首先定位新加入點所在的三角形，然后在網格中加入三個新的連接該三角形頂點與新頂點的邊，若該新點位于某條邊上，則該邊被刪除，四條連接該新點的邊被加入。最后，在通過換邊方法對該新點的局部區域內的邊進行檢測和變換，重新維護網格的Delaunay性質。局部變換法的另一個優點是其可以對已存在的三角網格進行優化，使其變換成為Delaunay三角網格，該方法的缺點則是當算法擴展到高維空間時變得較為復雜。

四、 Delaunay三角剖分在OpenCascade的應用

OpenCascade中網格剖分的包主要有BRepMesh、MeshAlgo、MeshVS，其中，類MeshAlgo_Delaunay使用算法Watson來進行Delaunay三角剖分。從類StlTransfer中的注釋The triangulation is computed with the Delaunay algorithm implemented in package BRepMesh.可以看出包BRepMesh就是Delaunay三角剖分的具體實現。使用方法如下：

BRepMesh::Mesh (aShape, Deflection);

這個函數主要是用來對拓撲形狀進行三角剖分。以下通過將一個圓柱三角剖分為例說明如何將一個拓撲形狀進行三角剖分并將結果進行可視化。

/**//** 

*    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(); 

} 

結果如下圖所示：

Figure 4.1 Cylinder mesh generated by BRepMesh::Mesh

BRepMesh::Mesh是經過封裝的，便于對拓撲形狀進行三角剖分。以下通過一個簡單的例子來說明直接使用BRepMesh_Delaun的方法：

上述程序是以一個正方形為例，使用BRepMesh_Delaun三角剖分的結果為兩個三角形，如下所示：

以上結果都是二維空間上的，三維空間中的使用方法可以參考類：BRepMesh_FastDiscretFace。這個類說明了如何將一個面進行網格劃分。

五、 結論

Delaunay三角剖分理論在三維幾何造型中還是比較重要的，通過對形狀的三角剖分，不僅可以對其進行可視化，還便于對形狀做進一步的處理，如消隱、光照處理等。通過對OpenCascade中三角剖分算法的使用，以進一步了解三角剖分理論應用及其算法實現。

六、 參考資料

1. 周培德. 計算幾何—算法設計與分析. 清華大學出版社, 2011

2. 李海生. Delaunay三角剖分理論及可視化應用研究. 哈爾濱工業大學出版社, 2010

3. 何援軍. 計算機圖形學. 機械工業出版社, 2010

4. 周元峰, 孫峰, 王文平, 汪嘉業, 張彩明. 基于局部修復的移動數據點Delaunay三角化快速更新方法. 計算機輔助設計與圖形學學報, 2011, 12: 2006-1012

5. http://en.wikipedia.org/wiki/Voronoi_diagram