OpenCascade中的Delaunay三角剖分

<a href="mailto:[email protected]">[email protected]</a>

摘要：本文簡要介紹了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三角剖分的應用領域十分廣泛：幾何模組化——用來尋找三維曲面“好的”三角剖分；有限元分析——用來生成“好的”有限元網格；地理資訊系統——用來進行空間領域分析；結晶學——用來确定合金的結構；人類學和考古學——用來确定氏族部落、首領權威、居住中心或堡壘等的影響範圍；天文學——用來确定恒星和星系的分布；生物學生态學和林學——用來确定動植物的競争；動物學——分析動物的領地；統計學和資料分析——用來分析統計聚合；機器人學——用來進行運動軌迹規劃（在存在障礙物的情況下）；模式識别——作為尋找物體骨架點的工具；生理學——用來分析毛細作用的領域；氣象學——用來估計區域平均降雨量；市場學——用來建立城市的市場輻射範圍；以及在遙感圖像處理、化學、地理學、地質學、冶金學、數學等學科的應用等。

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三角剖分的一條邊，反之亦成立。

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

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

l 沒有相交邊；

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

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

如果點集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

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三角網格，該方法的缺點則是當算法擴充到高維空間時變得較為複雜。

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  : [email protected] 

*        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&gt; 

#include gp_Pln.hxx&gt; 

#include BRep_Tool.hxx&gt; 

#include TopoDS.hxx&gt; 

#include TopoDS_Edge.hxx&gt; 

#include TopoDS_Wire.hxx&gt; 

#include TopoDS_Face.hxx&gt; 

#include BRepBuilderAPI_MakeEdge.hxx&gt; 

#include BRepBuilderAPI_MakeWire.hxx&gt; 

#include BRepBuilderAPI_MakeFace.hxx&gt; 

#include BRepPrimAPI_MakeBox.hxx&gt; 

#include BRepPrimAPI_MakeCone.hxx&gt; 

#include BRepPrimAPI_MakeCylinder.hxx&gt; 

#include BRepPrimApI_MakeSphere.hxx&gt; 

#include BRepMesh.hxx&gt; 

#include TopExp_Explorer.hxx&gt; 

#include Poly_Triangulation.hxx&gt; 

#include TShort_Array1OfShortReal.hxx&gt; 

#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&gt; 

#include osgViewer/Viewer&gt; 

#include osgViewer/ViewerEventHandlers&gt; 

#include osgGA/StateSetManipulator&gt; 

#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&amp; aShape) 

{ 

    osg::ref_ptrosg::Group&gt; root = new osg::Group(); 

    osg::ref_ptrosg::Geode&gt; geode = new osg::Geode(); 

    osg::ref_ptrosg::Geometry&gt; triGeom = new osg::Geometry(); 

    osg::ref_ptrosg::Vec3Array&gt; vertices = new osg::Vec3Array(); 

    osg::ref_ptrosg::Vec3Array&gt; 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-&gt;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-&gt;NbNodes()); 

        Poly_Array1OfTriangle triangles(1, triFace-&gt;NbTriangles()); 

        nodes = triFace-&gt;Nodes(); 

        triangles = triFace-&gt;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 &gt; gp::Resolution()) 

                normal.Normalize(); 

} 

else 

                normal.SetCoord(0., 0., 0.); 

            vertices-&gt;push_back(osg::Vec3(vertex1.X(), vertex1.Y(), vertex1.Z())); 

            vertices-&gt;push_back(osg::Vec3(vertex2.X(), vertex2.Y(), vertex2.Z())); 

            vertices-&gt;push_back(osg::Vec3(vertex3.X(), vertex3.Y(), vertex3.Z())); 

            normals-&gt;push_back(osg::Vec3(normal.X(), normal.Y(), normal.Z())); 

    triGeom-&gt;setVertexArray(vertices.get()); 

    triGeom-&gt;addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::TRIANGLES, 0, vertices-&gt;size())); 

    triGeom-&gt;setNormalArray(normals); 

    triGeom-&gt;setNormalBinding(osg::Geometry::BIND_PER_PRIMITIVE); 

    geode-&gt;addDrawable(triGeom); 

    root-&gt;addChild(geode); 

return root.release(); 

int main(int argc, char* argv[]) 

    osgViewer::Viewer myViewer; 

    root-&gt;addChild(BuildShapeMesh(BRepPrimAPI_MakeCylinder(.6, 1))); 

    myViewer.setSceneData(root); 

    myViewer.addEventHandler(new osgGA::StateSetManipulator(myViewer.getCamera()-&gt;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的方法：

/** 

#include BRepMesh_Edge.hxx&gt; 

#include BRepMesh_Delaun.hxx&gt; 

#include BRepMesh_Array1OfVertexOfDelaun.hxx&gt; 

#include TColStd_MapIteratorOfMapOfInteger.hxx&gt; 

    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-&gt;ElemOfDomain()); triDom.More(); triDom.Next()) 

        Standard_Integer triId = triDom.Key(); 

const BRepMesh_Triangle&amp; curTri = meshData-&gt;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&amp; edge1 = meshData-&gt;GetLink(edgeIndex1); 

const BRepMesh_Edge&amp; edge2 = meshData-&gt;GetLink(edgeIndex2); 

const BRepMesh_Edge&amp; edge3 = meshData-&gt;GetLink(edgeIndex3); 

        vertexIndex1 = (o1? edge1.FirstNode(): edge1.LastNode()); 

        vertexIndex2 = (o1? edge1.LastNode() : edge1.FirstNode()); 

        vertexIndex3 = (o2? edge2.LastNode() : edge2.FirstNode()); 

const BRepMesh_Vertex&amp; vertex1 = meshData-&gt;GetNode(vertexIndex1); 

const BRepMesh_Vertex&amp; vertex2 = meshData-&gt;GetNode(vertexIndex2); 

const BRepMesh_Vertex&amp; vertex3 = meshData-&gt;GetNode(vertexIndex3); 

const gp_XY&amp; p1 = vertex1.Coord(); 

const gp_XY&amp; p2 = vertex2.Coord(); 

const gp_XY&amp; p3 = vertex3.Coord(); 

        std::cout"--------"std::endl; 

        std::coutp1.X()" , "p1.Y()std::endl; 

        std::coutp2.X()" , "p2.Y()std::endl; 

        std::coutp3.X()" , "p3.Y()std::endl; 

        std::cout"========"std::endl; 

return 0; 

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

Iterate Mesh Triangles: 

-------- 

1 , 1 

0 , 0 

1 , 0 

======== 

0 , 1 

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

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

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

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

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

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