mobile.365-838.comDelaunay Triangulation in OpenCascade

Delaunay Triangulation in
OpenCascade

eryar@163.com

摘要:本文简要介绍了Delaunay三角剖分的基础理论,并使用OpenCascade的三角形剖分算法将边界BRep表示的几何体实行三角离散化后在OpenSceneGraph中显得。 

关键字:Delaunay Triangulation、OpenCascade、OpenSceneGraph 

一、 概述

三角形剖分是平面剖分中的3个根本课题,在数字图像处理、总结机三个维度曲面造型、有限元总结、逆向工程等领域有所广泛应用。由于三角形是平面域中的单纯形,与别的平面图形相比,其有描述方便、处理大约等特色,很吻合于对复杂区域开展简化处理。因而,无论在盘算几何、计算机图形处理、情势识别、曲面逼近,还有个别许元网格生成上边有广泛的使用。 

虽说曲线、曲面等有标准的方程来表示,不过在在计算机中,只能用离散的方法来逼近。如曲线可用直线段来逼近,而曲面可用多边形或三角形来代表。用多边形网格表示曲面是安顿性中平常利用的样式,能够根据使用须要选取网格的密度。利用三角形面片表示的曲面在计算机图形学中也叫做三角形网格。用三角网格表示曲面须求缓解多少个难点:三角形的发生、描述、遍历、简化和压缩等,那么些标题都是测算几何切磋的规模,相关题材都足以从中找到答案。下图所示的圆柱和立方体是由OpenCascade生成,使用OpenCascade的算法离散成三角网格后在OpenSceneGraph中显得的机能。 

mobile.365-838.com 1

Figure 1.1 Shaded Cylinder and Box 

mobile.365-838.com 2

Figure 1.2 Mesh generated by OpenCascade 

从图中能够看出,平面包车型客车三角形网格功用还不易,曲面包车型大巴三角网格表示只好是相仿表示,能够通过进步网格的密度来充实真实性,但对应渲染的数据量就大了。有人说OpenCascade的来得模块做得不是很好,上述方式则足以只行使OpenCascade的形状模块,再结合OpenSceneGraph来对图片进行体现。 

三维数据沟通STL格式文件中保留的都以三角面片的数目,STL文件格式是由美利哥3D
System集团开支,已被工产业界认为是日前迅猛机动成型领域的准标准零件描述文件格式。它对三维实体描述的分解具有惟一性。大致全部的几何样子系统都提供STL文件数据交换接口。OpenCascade中的数据调换模块也提供对STL格式的协助,由此可知三角网格在几何样子系统中的主要性。 

Voronoi图和Delaunay三角剖分的应用领域11分科学普及:几何建立模型——用来查找三个维度曲面“好的”三角剖分;有限元分析——用来扭转“好的”有限元网格;地理新闻种类——用来进行空中领域分析;结晶学——用来鲜明合金的协会;人类学和考古学——用来规定氏族部落、带头人权威、居住中央或堡垒等的震慑范围;天经济学——用来规定恒星和星系的遍布;生物学生态学和林学——用来分明动植物的竞争;动物学——分析动物的领地;总计学和数码解析——用来分析统计聚合;机器人学——用来展开移动轨迹规划(在存在障碍物的情事下);格局识别——作为寻找物体骨架点的工具;生教育学——用来分析毛细效能的圈子;气象学——用来推测区域平均降水量;市镇学——用来树立城市的市集辐射范围;以及在遥感图像处理、化学、地医学、地质学、冶金学、数学等学科的行使等。 

正文只对OpenCascade中的三角剖分进行简短介绍,希望对三角剖分在三个维度几何样子方面有趣味的对象能够对其深远钻研。水平很单薄,文中不当之处欢迎批评指正、引导,联系邮箱:eryar@163.com。 

二、 Voronoi图

Dirichlet于1850年探讨了平面点的邻域难点,Voronoi于一玖〇玖年将其结果扩张到高维空间。半空间定义Voronoi图:给定平面上n个点集S,S={p一,
p2, …, pn}。定义: 

mobile.365-838.com 3

PiPj连线的垂直平分面将空间分为两半,V(Pi)表示比其余点更类似Pi的点的轨迹是n-三个半平面包车型地铁交,它是3个不多于n-一条边的凸多边形域,称为关联于Pi的Voronoi多边形或关系于Pi的Voronoi域。如下图所示为关联于P一的Voronoi多边形,它是五个4边形,而n=陆. 

mobile.365-838.com 4

Figure 2.一 n=陆时的一种V(p壹) 

对此点集S中的各个点都能够做3个Voronoi多边形,那样的n个Voronoi多边形组成的图称为Voronoi图,记为Vor(S)。如下图所示: 

mobile.365-838.com 5

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

图中的顶点和边分外号字为Voronoi顶点和Voronoi边。鲜明,|S|=n时,Vor(S)划分平面成n个多边形域,每个多边形域V(Pi)包含S中的一个点同时只包括S中的3个点,Vor(S)的边是S中某点对的垂直平分线上的一条线条或半直线,从而为该点对所在的多少个多边形域所共有。Vor(S)中部分多边形域是无界的。 

mobile.365-838.com 6

Figure 2.3 Ten shops in a flat city and their Voronoi cells 

(http://en.wikipedia.org/wiki/Voronoi\_diagram

mobile.365-838.com 7

Figure 2.4 Voronoi tessellation in a cylinder (Voro++ library:
http://math.lbl.gov/voro++/

Voronoi图有如下性质: 

l n个点的点集S的Voronoi图至多有二n-四个极端和3n-陆条边; 

l 每一个Voronoi点恰好是3条Voronoi边的交点; 

l 设v是Vor(S)的终极,则圆C(v)内不含S的别的点; 

l 点集S中式点心Pi的每1个以来走近点明确V(Pi)的一条边; 

l Voronoi图的直线对偶图是S的三个三角剖分; 

l
假如Pi,Pj属于S,并且经过Pi,Pj有贰个不分包S中其余点的圆,那么线段PiPj是点集S三角剖分的一条边,反之亦建立。 

三、 Delaunay三角剖分 

一. 二维实数域上的三角剖分

倘诺V是贰维实数域上的有限点集,边e是由点集中的点作端点构成的封闭线段,E为e的聚众,那么该点集V的1个三角形剖分T=(V,E)是三个平面图: 

l 除了端点,平面图中的边不含有点集中的任何点; 

l 未有相交边; 

l 平面图中具有的面都是三角面,且富有三角面包车型大巴合集是点集V的凸包。 

2. Delaunay边

如果E中的一条边(两端点a,b),e满足下列原则,则称为Delaunay边:存在三个圆经过a,b两点,圆内不含有点集V中的任何的点。这一风味又称之为空圆性格。 

三. Delaunay三角剖分

如果点集V的三个三角剖分T中只包括Delaunay边,那么该三角剖分称为Delaunay剖分。 

近些年点意思下的Voronoi图的双双图实际上是点集的壹种三角剖分,该三角剖分正是Delaunay剖分(表示为DT(S)),当中各类三角形的外接圆不包罗点集中的别的任何点。由此,在构造点集的Voronoi图之后,再作其对偶图,即对每条Voronoi边作通过点集中某两点的垂直平分线,即获取Delaunay三角剖分。 

mobile.365-838.com 8

Figure 3.1 Delaunay Triangulation (Generated by MATLAB) 

再看多少个图片,加深对Delaunay三角剖分的知情: 

mobile.365-838.com 9

Figure 3.2 Delaunay Edge  

mobile.365-838.com 10

Figure 3.3 Illustrate Delaunay Edge 

mobile.365-838.com 11

Figure 3.4 Delaunay Edge 

肆. Delaunay三角剖分的表征

l
一9七陆年Sibson注解了在二维的场合下,在点集的保有三角剖分中,Delaunay三角剖分使得生成的三角形的最小角达到最大(max-min
angle)。因为那1特色,对于给确定地点集的Delaunay三角剖分总是竭尽制止“瘦长”三角形,自动向等边三角形逼近; 

l 局地优化与完整优化(locally optimal and globally optimal); 

l Delaunay空洞(cavity)与一些重连(local reconnection); 

  1. 经典的Delaunay三角剖分算法 

当下常用的算法分为二种,有扫描线法(Sweepline)、随机增量法(Incremental)、分治法(Divide
and Conquer)等。 

经文的Delaunay三角剖分算法主要有两类:Bowyer/沃特son算法和有些变换法。 

l
Bowyer/沃特son算法又称之为Delaunay空洞算法或加点法,以Bowyer和Watson算法为表示。从二个三角形最先,每一趟加3个点,保证每一步得到的当下三角形是①对优化的。以大不列颠及英格兰联合王国Bath大学数学分校Bowyer,格林,Sibson为代表的盘算Dirichlet图的艺术属于加点法,是较早成名的算法之1;以澳大哈利法克斯联邦(Commonwealth of Australia)莫斯科赫鲁高校学地球科学系沃特son为代表的空外接球法也属于加点法。加点法算法简明,是当下采取最多的算法,该措施应用了Delaunay空洞性质。Bowyer/沃特son算法的帮助和益处是与上空的维数非亲非故,并且算法在达成上比部分变换算法不难。该算法在新点进入到Delaunay网格时,部极度接球包括新点的三角形单元不再符合Delaunay属性,则那几个三角形单元被删去,形成Delaunay空洞,然后算法将新点与组合空洞的每三个巅峰相连生成多少个新边,依照空球属性能够评释这么些新边都以1对Delaunay的,因而新生成的三角形网格仍是Delaunay的。 

mobile.365-838.com 12

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

l
局地变换法又称作换边、换面法。当使用1些变换法完成增量式点集的Delaunay三角剖分时,首先定位新加入点所在的三角形,然后在网格中参预三个新的总是该三角形顶点与新顶点的边,若该新点位于某条边上,则该边被删去,4条连接该新点的边被到场。最后,在通过换边方法对该新点的一部分区域内的边进行检验和更换,重新维护网格的Delaunay性质。局地变换法的另多少个独到之处是其得以对已存在的三角网格进行优化,使其更换来为Delaunay三角网格,该办法的症结则是当算法扩充到高维空间时变得比较复杂。 

肆、 Delaunay三角剖分在OpenCascade的行使

OpenCascade中网格剖分的包首要有BRepMesh、MeshAlgo、MeshVS,在那之中,类MeshAlgo_Delaunay使用算法沃特son来进展Delaunay三角剖分。从类StlTransfer中的注释The
triangulation is computed with the Delaunay algorithm implemented in
package
BRepMesh.能够见见包BRepMesh就是Delaunay三角剖分的切切实实贯彻。使用格局如下: 

BRepMesh::Mesh (aShape, Deflection); 

本条函数主若是用来对拓扑形状举办三角剖分。以下通过将2个圆柱三角剖分为例表达什么将八个拓扑形状举办三角剖分并将结果开始展览可视化。 

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

结果如下图所示: 

mobile.365-838.com 13

Figure 4.1 Cylinder mesh generated by BRepMesh::Mesh 

BRepMesh::Mesh是由此包装的,便于对拓扑形状举办三角剖分。以下通过一个回顾的事例来注解直接利用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;
}

上述顺序是以三个星型为例,使用BRepMesh_mobile.365-838.com,Delaun三角剖分的结果为五个三角形,如下所示: 

Iterate Mesh Triangles: 
——– 
1 , 1 
0 , 0 
1 , 0 
======== 
——– 
1 , 1 
0 , 1 
0 , 0 
======== 

 以上结果都是二维空间上的,三个维度空间中的使用办法可以参考类:BRepMesh_法斯特DiscretFace。那一个类表达了哪些将多少个面拓展网格划分。 

五、 结论

Delaunay三角剖分理论在三个维度几何样子中照旧比较首要的,通过对造型的三角剖分,不仅能够对其举行可视化,还有利于对造型做尤其的处理,如消隐、光照处理等。通过对OpenCascade中三角剖分算法的利用,以进一步精晓三角剖分理论运用及其算法完结。 

六、 参考资料

  1. 周培德. 总结几何—算法设计与分析. 清华东军事和政院学出版社, 201一 

  2. 李海生. Delaunay三角剖分理论及可视化应用研究. 奇瓦瓦药科高校出版社,
    2010 

  3. 何援军. 计算机图形学. 机械工业出版社, 20十 

  4. 周元峰, 孙峰, 王文平, 汪嘉业, 张彩明.
    基于局地修复的移位数据点Delaunay三角化急忙翻新方法.
    计算机帮忙设计与图片学学报, 201一, 1二: 200陆-101二 

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

 

PDF Version: Delaunay Triangulation in
OpenCascade

admin

网站地图xml地图