Delaunay Triangulation in OpenCascade

Size: px

Start display at page:

Download "Delaunay Triangulation in OpenCascade"

含鲁
5 years ago
Views:

Delaunay Triangulation in OpenCascade eryar@163.

1 Delaunay Triangulation in OpenCascade 摘要 : 本文简要介绍了 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

2 从图中可以看出, 平面的三角形网格效果还不错, 曲面的三角形网格表示只能是近似表示, 可以通过提高网格的密度来增加真实性, 但相应渲染的数据量就大了有人说 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} 定义 : V ( p i ) = Ι i j H( p, p PiPj 连线的垂直平分面将空间分为两半,V(Pi) 表示比其他点更接近 Pi 的点的轨迹是

3 二 Voronoi 图 Dirichlet 于 1850 年研究了平面点的邻域问题,Voronoi 于 1908 年将其结果扩展到高维空间半空间定义 Voronoi 图 : 给定平面上 n 个点集 S,S={p1, p2,, pn} 定义 : V ( p i ) = Ι i j H( p, p PiPj 连线的垂直平分面将空间分为两半,V(Pi) 表示比其他点更接近 Pi 的点的轨迹是 n-1 个半平面的交, 它是一个不多于 n-1 条边的凸多边形域, 称为关联于 Pi 的 Voronoi 多边形或关联于 Pi 的 Voronoi 域如下图所示为关联于 P1 的 Voronoi 多边形, 它是一个四边形, 而 n=6. i j ) 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.

gov/voro++/) Voronoi 图有如下性质 : n 个点的点集 S 的 Voronoi 图至多有 2n-5 个顶点和 3n-6 条边 ; 每个 Voronoi 点恰好是三条 Voronoi 边的交点 ; 设 v 是 Vor(S)

4 共有 Vor(S) 中有的多边形域是无界的 Figure 2.3 Ten shops in a flat city and their Voronoi cells ( Figure 2.4 Voronoi tessellation in a cylinder (Voro++ library: Voronoi 图有如下性质 : n 个点的点集 S 的 Voronoi 图至多有 2n-5 个顶点和 3n-6 条边 ; 每个 Voronoi 点恰好是三条 Voronoi 边的交点 ; 设 v 是 Vor(S) 的顶点, 则圆 C(v) 内不含 S 的其他点 ; 点集 S 中点 Pi 的每一个最近邻近点确定 V(Pi) 的一条边 ; Voronoi 图的直线对偶图是 S 的一个三角剖分 ; 如果 Pi,Pj 属于 S, 并且通过 Pi,Pj 有一个不包含 S 中其他点的圆, 那么线段 PiPj 是点集 S 三角剖分的一条边, 反之亦成立

三 Delaunay 三角剖分 1. 二维实数域上的三角剖分假设 V 是二维实数域上的有限点集, 边 e 是由点集中的点作端点构成的封闭线段,E 为 e 的集合, 那么该点集 V 的一个三角剖分 T=(V,E) 是一个平面图 : 除了端点, 平面图中的边不包含点集中的任何点 ; 没有相交边 ; 平面图中所有的面都是三角面, 且所有三角面的合集是点集 V 的凸包 2.

5 三 Delaunay 三角剖分 1. 二维实数域上的三角剖分假设 V 是二维实数域上的有限点集, 边 e 是由点集中的点作端点构成的封闭线段,E 为 e 的集合, 那么该点集 V 的一个三角剖分 T=(V,E) 是一个平面图 : 除了端点, 平面图中的边不包含点集中的任何点 ; 没有相交边 ; 平面图中所有的面都是三角面, 且所有三角面的合集是点集 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 三角剖分的理解 :

6 Figure 3.2 Delaunay Edge Figure 3.3 Illustrate Delaunay Edge

7 Figure 3.4 Delaunay Edge 4. Delaunay 三角剖分的特性 1978 年 Sibson 证明了在二维的情况下, 在点集的所有三角剖分中,Delaunay 三角剖分使得生成的三角形的最小角达到最大 (max-min angle) 因为这一特性, 对于给定点集的 Delaunay 三角剖分总是尽可能避免瘦长三角形, 自动向等边三角形逼近 ; 局部优化与整体优化 (locally optimal and globally optimal); Delaunay 空洞 (cavity) 与局部重连 (local reconnection); 5. 经典的 Delaunay 三角剖分算法目前常用的算法分为几种, 有扫描线法 (Sweepline) 随机增量法 (Incremental) 分治法 (Divide and Conquer) 等经典的 Delaunay 三角剖分算法主要有两类 :Bowyer/Watson 算法和局部变换法 Bowyer/Watson 算法又称为 Delaunay 空洞算法或加点法, 以 Bowyer 和 Watson 算法为代表从一个三角形开始, 每次加一个点, 保证每一步得到的当前三角形是局部优化的以英国 Bath 大学数学分校 Bowyer,Green,Sibson 为代表的计算 Dirichlet 图的方法属于加点法, 是较早成名的算法之一 ; 以澳大利亚悉尼大学地学系 Watson 为代表的空外接球法也属于加点法加点法算法简明, 是目前应用最多的算法, 该方法利用了 Delaunay 空洞性质 Bowyer/Watson 算法的优点是与空间的维数无关, 并且算法在实现上比局部变换算法简单该算法在新点加入到 Delaunay 网格时, 部分外接球包含新点的三角形单元不再符合 Delaunay 属性, 则这些三角形单元被删除, 形成 Delaunay 空洞, 然后算法将新点与组成空洞的每一个顶点相连生成一个新边, 根据空球属性可以证明这些新边都是局部 Delaunay 的, 因此新生成的三角网格仍是 Delaunay 的

8 Figure 3.5 Illustration of 2D Bowyer/Watson algorithm for Delaunay Triangulation 局部变换法又称为换边换面法当利用局部变换法实现增量式点集的 Delaunay 三角剖分时, 首先定位新加入点所在的三角形, 然后在网格中加入三个新的连接该三角形顶点与新顶点的边, 若该新点位于某条边上, 则该边被删除, 四条连接该新点的边被加入最后, 在通过换边方法对该新点的局部区域内的边进行检测和变换, 重新维护网格的 Delaunay 性质局部变换法的另一个优点是其可以对已存在的三角网格进行优化, 使其变换成为 Delaunay 三角网格, 该方法的缺点则是当算法扩展到高维空间时变得较为复杂

9 四 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 : * 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")

10 #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;

11 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())); } }

12 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 new osgga:: ::StateSetManipulator(myViewer.getCamera()-> ()->getorcreatestateset() )); myviewer.addeventhandler(new new osgviewer:: ::StatsHandler); myviewer.addeventhandler(new new osgviewer:: ::WindowSizeHandler); return myviewer.run(); } 结果如下图所示 : Figure 4.1 Cylinder mesh generated by BRepMesh::Mesh

13 BRepMesh::Mesh 是经过封装的, 便于对拓扑形状进行三角剖分以下通过一个简单的例子来说明直接使用 BRepMesh_Delaun 的方法 : /** * Copyright (c) 2013 eryar All Rights Reserved. * * File : Main.cpp * Author : eryar@163.com * Date : * 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()) {

14 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;

15 } 上述程序是以一个正方形为例, 使用 BRepMesh_Delaun 三角剖分的结果为两个三角形, 如下所示 : Iterate Mesh Triangles: , 1 0, 0 1, 0 ======== , 1 0, 1 0, 0 ======== 以上结果都是二维空间上的, 三维空间中的使用方法可以参考类 : BRepMesh_FastDiscretFace 这个类说明了如何将一个面进行网格划分

16 五结论 Delaunay 三角剖分理论在三维几何造型中还是比较重要的, 通过对形状的三角剖分, 不仅可以对其进行可视化, 还便于对形状做进一步的处理, 如消隐光照处理等通过对 OpenCascade 中三角剖分算法的使用, 以进一步了解三角剖分理论应用及其算法实现

17 六参考资料 1. 周培德. 计算几何算法设计与分析. 清华大学出版社, 李海生. Delaunay 三角剖分理论及可视化应用研究. 哈尔滨工业大学出版社, 何援军. 计算机图形学. 机械工业出版社, 周元峰, 孙峰, 王文平, 汪嘉业, 张彩明. 基于局部修复的移动数据点 Delaunay 三角化快速更新方法. 计算机辅助设计与图形学学报, 2011, 12:

绘制OpenCascade中的曲线

绘制OpenCascade中的曲线在 OpenSceneGraph 中绘制 OpenCascade 的曲线 Draw OpenCascade Geometry Curves in OpenSceneGraph eryar@163.com 摘要 Abstract: 本文简要说明 OpenCascade 中几何曲线的数据, 并将这些几何曲线在 OpenSceneGraph 中绘制出来关键字 KeyWords:OpenCascade Geometry