OpenCascade与三角剖分 笔记

笔记

摘自 http://www.cppblog.com/eryar/archive/2013/05/26/200605.html

概述

三角剖分是平面剖分中的一个重要课题，在数字图像处理、计算机三维曲面造型、有限元计算、逆向工程等领域有着广泛应用。由于三角形是平面域中的单纯形，与其他平面图形相比，其有描述方便、处理简单等特性，很适合于对复杂区域进行简化处理。因此，无论在计算几何、计算机图形处理、模式识别、曲面逼近，还有有限元网格生成方面有广泛的应用。

虽然曲线、曲面等有精确的方程来表示，但是在在计算机中，只能用离散的方式来逼近。如曲线可用直线段来逼近，而曲面可用多边形或三角形来表示。用多边形网格表示曲面是设计中经常使用的形式，可以根据应用要求选择网格的密度。利用三角形面片表示的曲面在计算机图形学中也称为三角形网格。用三角形网格表示曲面需要解决几个问题：三角形的产生、描述、遍历、简化和压缩等，这些问题都是计算几何研究的范畴，相关问题都可以从中找到答案。下图所示的圆柱和立方体是由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格式的支持，由此可见三角网格在几何造型系统中的重要性。

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