OpenCASCADE曲面求交之网格离散法1

eryar@163.com

1 Introduction

由朱心雄等著《自由曲线曲面造型技术》书中对曲面求交之网格离散法描述如下：该法的基本思想是先将曲面离散为由小平面片组成的网格，当网格足够密时，可以认为已经非常接近真实曲面，对分别表示不同曲面的两张网格，利用平面片求交法求得的交线，并以此交线近似代表曲面间的交线。这种方法原理简明，便于实现，适用范围广，任意参数曲面均可利用该法求交。但为获取精确地交线，则必须生成非常细密的网格，这将导致占用内存多，计算花费大。因此，实际工作中很少单一使用离散网格法，而常将其与其他方法结合使用。

OpenCASCADE中对于曲面求交也提供离散网格法，其中曲面的离散网格由类IntPatch_Polyhedron表示，两个网格面求交使用类IntPatch_InterferencePolyhedron。本文主要介绍曲面的网格表示类IntPatch_Polyhedron。

2 Polyhedron

OpenCASCADE用于曲面求交的网格离散算法相对BRepMesh中的算法要简单很多，主要思路是根据参数U，V方向上的采样点数量来计算曲面上的点，再根据固定公式将采样点连成三角形。其中生成采样点代码如下所示：

成员变量CMyPnts是采样点数组，CMyU和CMyV是采样点在曲面上的参数。将采样点连成三角形函数如下图所示：

根据上述生成采样点及三角形函数，对于平面生成的三角形如下图所示：

其中Triangle()函数中变量line表示参数u方向上第几条线，代入具体的索引Index来看规律：

当参数索引 Index为1时，line为1，得到的三角形为1-4-5；

当参数索引Index为2时，line为1，得到的三角形为1-5-2；

当参数索引Index为3时，line为1，得到的三角形为2-5-6；

当参数索引Index为4时，line为1，得到的三角形为2-6-3；

当参数索引Index为5时，line为2，得到的三角形为4-7-8；

当参数索引Index为6时，line为2，得到的三角形为4-8-5；

从上可以得到生成三角形的规律，即根据索引Index计算正在处理的三角形是参数u方向上第几条线line，生成这条线上在参数v方向上的所有的三角形。生成的三角形都是逆时针的。

下面我们看看对于一些基本曲面，这种离散网格算法生成的网格效果：

球面的离散网格

圆柱面的离散网格

圆环面的离散网格

B样条曲面

3 Conclusion

综上所述，类IntPatch_Polyhedron中生成网格的算法主要依赖曲面在参数U，V上的采样点数量。默认采样点数量是根据函数NbSamplesV()和NbSamplesU()生成。

也可以指定采样点数量，当采样点数量越多，则生成的三角形越多，网格越密。当然这种方式也可用来生成曲面的显示数据，生成速度很快，唯一的缺陷是生成显示用网格的精度只能通过采样点数量来控制，对于曲率变化大的曲面，若指定多的采样点，则会生成大量三角形，占用大量内存空间。

附上测试代码：

#include <TColgp_Array2OfPnt.hxx>
#include <Geom_Plane.hxx>
#include <Geom_CylindricalSurface.hxx>
#include <Geom_ConicalSurface.hxx>
#include <Geom_SphericalSurface.hxx>
#include <Geom_ToroidalSurface.hxx>
#include <Geom_BSplineSurface.hxx>
#include <GeomAdaptor_Surface.hxx>
#include <GeomAPI_PointsToBSplineSurface.hxx>
#include <IntPatch_Polyhedron.hxx>
#include <IntPatch_InterferencePolyhedron.hxx>
#pragma comment(lib, "TKernel.lib")
#pragma comment(lib, "TKMath.lib")
#pragma comment(lib, "TKG2d.lib")
#pragma comment(lib, "TKG3d.lib")
#pragma comment(lib, "TKGeomBase.lib")
#pragma comment(lib, "TKGeomAlgo.lib")
void makeSurface(Handle(Geom_BSplineSurface)& theSurface)
{
    TColgp_Array2OfPnt aPoints(1, 5, 1, 5);
    aPoints.SetValue(1, 1, gp_Pnt(-4, -4, 5));
    aPoints.SetValue(1, 2, gp_Pnt(-4, -2, 5));
    aPoints.SetValue(1, 3, gp_Pnt(-4, 0, 4));
    aPoints.SetValue(1, 4, gp_Pnt(-4, 2, 5));
    aPoints.SetValue(1, 5, gp_Pnt(-4, 4, 5));
    aPoints.SetValue(2, 1, gp_Pnt(-2, -4, 4));
    aPoints.SetValue(2, 2, gp_Pnt(-2, -2, 4));
    aPoints.SetValue(2, 3, gp_Pnt(-2, 0, 4));
    aPoints.SetValue(2, 4, gp_Pnt(-2, 2, 4));
    aPoints.SetValue(2, 5, gp_Pnt(-2, 5, 4));
    aPoints.SetValue(3, 1, gp_Pnt(0, -4, 3.5));
    aPoints.SetValue(3, 2, gp_Pnt(0, -2, 3.5));
    aPoints.SetValue(3, 3, gp_Pnt(0, 0, 3.5));
    aPoints.SetValue(3, 4, gp_Pnt(0, 2, 3.5));
    aPoints.SetValue(3, 5, gp_Pnt(0, 5, 3.5));
    aPoints.SetValue(4, 1, gp_Pnt(2, -4, 4));
    aPoints.SetValue(4, 2, gp_Pnt(2, -2, 4));
    aPoints.SetValue(4, 3, gp_Pnt(2, 0, 3.5));
    aPoints.SetValue(4, 4, gp_Pnt(2, 2, 5));
    aPoints.SetValue(4, 5, gp_Pnt(2, 5, 4));
    aPoints.SetValue(5, 1, gp_Pnt(4, -4, 5));
    aPoints.SetValue(5, 2, gp_Pnt(4, -2, 5));
    aPoints.SetValue(5, 3, gp_Pnt(4, 0, 5));
    aPoints.SetValue(5, 4, gp_Pnt(4, 2, 6));
    aPoints.SetValue(5, 5, gp_Pnt(4, 5, 5));
    theSurface = GeomAPI_PointsToBSplineSurface(aPoints).Surface();
}
void writeStl(const IntPatch_Polyhedron& thePolyhedron, const std::string& theFileName)
{
    // Dump surface polyhedron to STL file.
    std::ofstream aStlFile(theFileName);
    aStlFile << "solid polyhedron" << std::endl;
    // Dump triangles.
    for (Standard_Integer t = 1; t <= thePolyhedron.NbTriangles(); ++t)
    {
        Standard_Integer aPi1 = 0;
        Standard_Integer aPi2 = 0;
        Standard_Integer aPi3 = 0;
        thePolyhedron.Triangle(t, aPi1, aPi2, aPi3);
        const gp_Pnt& aP1 = thePolyhedron.Point(aPi1);
        const gp_Pnt& aP2 = thePolyhedron.Point(aPi2);
        const gp_Pnt& aP3 = thePolyhedron.Point(aPi3);
        aStlFile << "facet" << std::endl;
        aStlFile << "outer loop" << std::endl;
        aStlFile << "vertex " << aP1.X() << " " << aP1.Y() << " " << aP1.Z() << std::endl;
        aStlFile << "vertex " << aP2.X() << " " << aP2.Y() << " " << aP2.Z() << std::endl;
        aStlFile << "vertex " << aP3.X() << " " << aP3.Y() << " " << aP3.Z() << std::endl;
        aStlFile << "endloop" << std::endl;
        aStlFile << "endfacet" << std::endl;
    }
    aStlFile << "endsolid polyhedron" << std::endl;
    aStlFile.close();
}
void testPolyhedron()
{
    // Plane surface polyhedron.
    Handle(Geom_Plane) aPlane = new Geom_Plane(gp::XOY());
    Handle(GeomAdaptor_Surface) aSurfaceAdaptor = new GeomAdaptor_Surface(aPlane, 0.0, 10.0, 0.0, 20.0);
    IntPatch_Polyhedron aPlanePolyhedron(aSurfaceAdaptor);
    writeStl(aPlanePolyhedron, "d:/plane.stl");
    // Spherical surface polyhedron.
    Handle(Geom_SphericalSurface) aSphericalSurface = new Geom_SphericalSurface(gp::XOY(), 3.0);
    aSurfaceAdaptor = new GeomAdaptor_Surface(aSphericalSurface);
    IntPatch_Polyhedron aSphericalPolyhedron(aSurfaceAdaptor);
    writeStl(aSphericalPolyhedron, "d:/spherical.stl");
    // Cylindrical surface polyhedron.
    Handle(Geom_CylindricalSurface) aCylindricalSurface = new Geom_CylindricalSurface(gp::XOY(), 5.0);
    aSurfaceAdaptor = new GeomAdaptor_Surface(aCylindricalSurface, 0.0, M_PI, 0.0, 8.0);
    IntPatch_Polyhedron aCylindricalPolyhedron(aSurfaceAdaptor);
    writeStl(aCylindricalPolyhedron, "d:/cylindrical.stl");
    // Toroidal Surface polyhedron.
    Handle(Geom_ToroidalSurface) aToroidalSurface = new Geom_ToroidalSurface(gp::XOY(), 10.0, 3.0);
    aSurfaceAdaptor = new GeomAdaptor_Surface(aToroidalSurface);
    IntPatch_Polyhedron aToroidalPolyhedron(aSurfaceAdaptor);
    writeStl(aToroidalPolyhedron, "d:/toroidal.stl");
    // BSpline surface polyhedron.
    Handle(Geom_BSplineSurface) aBSplineSurface;
    makeSurface(aBSplineSurface);
    aSurfaceAdaptor = new GeomAdaptor_Surface(aBSplineSurface);
    IntPatch_Polyhedron aPolyhedron(aSurfaceAdaptor);
    writeStl(aPolyhedron, "d:/bspline.stl");
}
int main(int argc, char* argv[])
{
    testPolyhedron();
    return 0;
}