OpenCASCADE 参数曲面面积
eryar@163.com
Abstract. 本文介绍了参数曲面的第一基本公式，并应用曲面的第一基本公式，结合OpenCASCADE中计算多重积分的类，对任意参数曲面的面积进行计算。
Key Words. Parametric Curve, Parametric Surface, Gauss Integration, Global Properties
1.Introduction
我们知道一元函数y=f(x)的图像是一条曲线，二元函数z=f(x,y)的图像是一张曲面。但是，把曲线曲面表示成参数方程则更加便利于研究，这种表示方法首先是由欧洲瑞士数学家Euler引进的。例如，在空间中的一条曲线可以表示为三个一元函数：
X=x(t), Y=y(t), Z=z(t)
在向量的概念出现后，空间中的一条曲线可以自然地表示为一个一元向量函数：
r=r(t)=(x(t), y(t), z(t))
用向量函数来表示曲线和曲面后，使曲线曲面一些量的计算方式比较统一。如曲线可以表示为一元向量函数，曲面可以表示为二元向量函数。
本文结合OpenCASCADE来介绍参数曲线曲面积分分别计算曲线弧长和曲面的面积。结合《微分几何》来更好地理解曲线曲面相关知识。
2.Curve Natural Parametric Equations
设曲线C的参数方程是r=r(t)，命：
则s是该曲线的一个不变量，即它与空间中的坐标系的选择无关，也与该曲线的参数变换无关。前者是因为在笛卡尔直角坐标变换下，切向量的长度|r’(t)|是不变的，故s不变。关于后者可以通过积分的变量的替换来证明，设参数变换是：
并且
因此
根据积分的变量替换公式有：
不变量s的几何意义就是曲线段的弧长。这说明曲线参数t可以是任意的，但选择不同的参数得到的参数方程会有不同，但是曲线段的弧长是不变的。以曲线弧长作为曲线方程的参数，这样的方程称为曲线的自然参数方程Natural Parametric Equations。
由曲线的参数方程可知，曲线弧长的计算公式为：
几何意义就是在每个微元处的切向量的长度求和。
3.Gauss Integration for Arc Length
曲线弧长的计算就是一元函数的积分。OpenCASCADE中是如何计算任意曲线弧长的呢？直接找到相关的源码列举如下：(在类CPnts_AbscissaPoint中)
// auxiliary functions to compute the length of the derivative
static Standard_Real f3d(const Standard_Real X, const Standard_Address C)
{
gp_Pnt P;
gp_Vec V;
((Adaptor3d_Curve*)C)->D1(X,P,V);
return V.Magnitude();
}
static Standard_Real f2d(const Standard_Real X, const Standard_Address C)
{
gp_Pnt2d P;
gp_Vec2d V;
((Adaptor2d_Curve2d*)C)->D1(X,P,V);
return V.Magnitude();
}
//==================================================================
//function : Length
//purpose : 3d with parameters
//==================================================================
Standard_Real CPnts_AbscissaPoint::Length(const Adaptor3d_Curve& C,
const Standard_Real U1,
const Standard_Real U2)
{
CPnts_MyGaussFunction FG;
//POP pout WNT
CPnts_RealFunction rf = f3d;
FG.Init(rf,(Standard_Address)&C);
// FG.Init(f3d,(Standard_Address)&C);
math_GaussSingleIntegration TheLength(FG, U1, U2, order(C));
if (!TheLength.IsDone()) {
throw Standard_ConstructionError();
}
return Abs(TheLength.Value());
}
//==================================================================
//function : Length
//purpose : 2d with parameters
//==================================================================
Standard_Real CPnts_AbscissaPoint::Length(const Adaptor2d_Curve2d& C,
const Standard_Real U1,
const Standard_Real U2)
{
CPnts_MyGaussFunction FG;
//POP pout WNT
CPnts_RealFunction rf = f2d;
FG.Init(rf,(Standard_Address)&C);
// FG.Init(f2d,(Standard_Address)&C);
math_GaussSingleIntegration TheLength(FG, U1, U2, order(C));
if (!TheLength.IsDone()) {
throw Standard_ConstructionError();
}
return Abs(TheLength.Value());
}
上述代码的意思是直接对曲线的一阶导数的长度求积分，即是弧长。OpenCASCADE的代码写得有点难懂，根据意思把对三维曲线求弧长的代码改写下，更便于理解：
//! Function for curve length evaluation.
class math_LengthFunction : public math_Function
{
public:
math_LengthFunction(const Handle(Geom_Curve)& theCurve)
: myCurve(theCurve)
{
}
virtual Standard_Boolean Value(const Standard_Real X, Standard_Real& F)
{
gp_Pnt aP;
gp_Vec aV;
myCurve->D1(X, aP, aV);
F = aV.Magnitude();
return Standard_True;
}
private:
Handle(Geom_Curve) myCurve;
};
4.First Fundamental Form of a Surface
曲面参数方程是个二元向量函数。根据《微分几何》中曲面的第一基本公式（First Fundamental Form of a Surface）可知，曲面上曲线的表达式为：
r=r(u(t), v(t)) = (x(t), y(t), z(t))
若以s表示曲面上曲线的弧长，则由复合函数求导公式可得弧长微分公式：
在古典微分几何中，上式称为曲面的第一基本公式，E、F、G称为第一基本量。在曲面上，每一点的第一基本量与参数化无关。
利用曲面第一基本公式可以用于计算曲面的面积。参数曲面上与u,v参数空间的元素dudv对应的面积元为：
由参数曲面法向的计算可知，曲面的面积元素即为u,v方向上的偏导数的乘积的模。
其几何意义可以理解为参数曲面的面积微元是由u,v方向的偏导数的向量围成的一个四边形的面积，则整个曲面的面积即是对面积元素求积分。由于参数曲面有两个参数，所以若要计算曲面的面积，只需要对面积元素计算二重积分即可。
5.Gauss Integration for Area
OpenCASCADE的math包中提供了多重积分的计算类math_GaussMultipleIntegration，由类名可知积分算法采用了Gauss积分算法。下面通过具体代码来说明OpenCASCADE中计算曲面积分的过程。要计算积分，先要定义被积函数。因为参数曲面与参数曲线不同，参数曲线只有一个参数，而参数曲面有两个参数，所以是一个多元函数。
//! 2D variable function for surface area evaluation.
class math_AreaFunction : public math_MultipleVarFunction
{
public:
math_AreaFunction(const Handle(Geom_Surface)& theSurface)
: mySurface(theSurface)
{
}
virtual Standard_Integer NbVariables() const
{
return 2;
}
virtual Standard_Boolean Value(const math_Vector& X, Standard_Real& Y)
{
gp_Pnt aP;
gp_Vec aDu;
gp_Vec aDv;
mySurface->D1(X(1), X(2), aP, aDu, aDv);
Y = aDu.Crossed(aDv).Magnitude();
return Standard_True;
}
private:
Handle(Geom_Surface) mySurface;
};
由于参数曲面是多元函数，所以从类math_MultipleVarFunction派生，并在虚函数NbVariables()中说明有两个变量。在虚函数Value()中计算面积元素的值，即根据曲面第一基本公式中面积元素的定义，对参数曲面求一阶导数，计算两个偏导数向量的叉乘的模。
有了被积函数，只需要在定义域对其计算二重积分，相应代码如下所示：
void evalArea(const Handle(Geom_Surface)& theSurface, const math_Vector& theLower, const math_Vector& theUpper)
{
math_IntegerVector aOrder(1, 2, math::GaussPointsMax());
math_AreaFunction aFunction(theSurface);
math_GaussMultipleIntegration anIntegral(aFunction, theLower, theUpper, aOrder);
if (anIntegral.IsDone())
{
anIntegral.Dump(std::cout);
}
}
通过theLower和theUpper指定定义域，由于采用了Gauss-Legendre算法计算二重积分，所以需要指定阶数，且阶数越高积分结果精度越高，这里使用了OpenCASCADE中最高的阶数。
下面通过对基本曲面的面积计算来验证结果的正确性，并将计算结果和OpenCASCADE中计算面积的类BRepGProp::SurfaceProperties()结果进行对比。
6.Elementary Surface Area Test
下面通过对OpenCASCADE中几个初等曲面的面积进行计算，代码如下所示：
/*
Copyright(C) 2017 Shing Liu(eryar@163.com)
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files(the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and / or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions :
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/
#include <gp_Pnt.hxx>
#include <gp_Vec.hxx>
#include <math.hxx>
#include <math_Function.hxx>
#include <math_MultipleVarFunction.hxx>
#include <math_GaussMultipleIntegration.hxx>
#include <Geom_Plane.hxx>
#include <Geom_ConicalSurface.hxx>
#include <Geom_CylindricalSurface.hxx>
#include <Geom_SphericalSurface.hxx>
#include <Geom_ToroidalSurface.hxx>
#include <Geom_BSplineSurface.hxx>
#include <Geom_RectangularTrimmedSurface.hxx>
#include <GeomConvert.hxx>
#include <GProp_GProps.hxx>
#include <TopoDS_Face.hxx>
#include <BRepGProp.hxx>
#include <BRepBuilderAPI_MakeFace.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")
#pragma comment(lib, "TKBRep.lib")
#pragma comment(lib, "TKTopAlgo.lib")
//! 2D variable function for surface area evaluation.
class math_AreaFunction : public math_MultipleVarFunction
{
public:
math_AreaFunction(const Handle(Geom_Surface)& theSurface)
: mySurface(theSurface)
{
}
virtual Standard_Integer NbVariables() const
{
return 2;
}
virtual Standard_Boolean Value(const math_Vector& X, Standard_Real& Y)
{
gp_Pnt aP;
gp_Vec aDu;
gp_Vec aDv;
Standard_Real E = 0.0;
Standard_Real F = 0.0;
Standard_Real G = 0.0;
mySurface->D1(X(1), X(2), aP, aDu, aDv);
E = aDu.Dot(aDu);
F = aDu.Dot(aDv);
G = aDv.Dot(aDv);
Y = Sqrt(E * G - F * F);
//Y = aDu.Crossed(aDv).Magnitude();
return Standard_True;
}
private:
Handle(Geom_Surface) mySurface;
};
void evalArea(const Handle(Geom_Surface)& theSurface, const math_Vector& theLower, const math_Vector& theUpper)
{
math_IntegerVector aOrder(1, 2, math::GaussPointsMax());
math_AreaFunction aFunction(theSurface);
math_GaussMultipleIntegration anIntegral(aFunction, theLower, theUpper, aOrder);
if (anIntegral.IsDone())
{
anIntegral.Dump(std::cout);
}
}
void evalArea(const Handle(Geom_BoundedSurface)& theSurface)
{
math_IntegerVector aOrder(1, 2, math::GaussPointsMax());
math_Vector aLower(1, 2, 0.0);
math_Vector aUpper(1, 2, 0.0);
theSurface->Bounds(aLower(1), aUpper(1), aLower(2), aUpper(2));
math_AreaFunction aFunction(theSurface);
math_GaussMultipleIntegration anIntegral(aFunction, aLower, aUpper, aOrder);
if (anIntegral.IsDone())
{
anIntegral.Dump(std::cout);
}
}
void testFace(const TopoDS_Shape& theFace)
{
GProp_GProps aSurfaceProps;
BRepGProp::SurfaceProperties(theFace, aSurfaceProps);
std::cout << "Face area: " << aSurfaceProps.Mass() << std::endl;
}
void testPlane()
{
std::cout << "====== Test Plane Area =====" << std::endl;
Handle(Geom_Plane) aPlaneSurface = new Geom_Plane(gp::XOY());
math_Vector aLower(1, 2);
math_Vector aUpper(1, 2);
// Parameter U range.
aLower(1) = 0.0;
aUpper(1) = 2.0;
// Parameter V range.
aLower(2) = 0.0;
aUpper(2) = 3.0;
evalArea(aPlaneSurface, aLower, aUpper);
// Convert to BSpline Surface.
Handle(Geom_RectangularTrimmedSurface) aTrimmedSurface =
new Geom_RectangularTrimmedSurface(aPlaneSurface, aLower(1), aUpper(1), aLower(2), aUpper(2));
Handle(Geom_BSplineSurface) aBSplineSurface = GeomConvert::SurfaceToBSplineSurface(aTrimmedSurface);
evalArea(aBSplineSurface);
// Test Face.
TopoDS_Face aFace = BRepBuilderAPI_MakeFace(aTrimmedSurface, Precision::Confusion()).Face();
testFace(aFace);
aFace = BRepBuilderAPI_MakeFace(aBSplineSurface, Precision::Confusion()).Face();
testFace(aFace);
}
void testCylinder()
{
std::cout << "====== Test Cylinder Area =====" << std::endl;
Handle(Geom_CylindricalSurface) aCylindrialSurface = new Geom_CylindricalSurface(gp::XOY(), 1.0);
math_Vector aLower(1, 2);
math_Vector aUpper(1, 2);
aLower(1) = 0.0;
aUpper(1) = M_PI * 2.0;
aLower(2) = 0.0;
aUpper(2) = 3.0;
evalArea(aCylindrialSurface, aLower, aUpper);
// Convert to BSpline Surface.
Handle(Geom_RectangularTrimmedSurface) aTrimmedSurface =
new Geom_RectangularTrimmedSurface(aCylindrialSurface, aLower(1), aUpper(1), aLower(2), aUpper(2));
Handle(Geom_BSplineSurface) aBSplineSurface = GeomConvert::SurfaceToBSplineSurface(aTrimmedSurface);
evalArea(aBSplineSurface);
// Test Face.
TopoDS_Face aFace = BRepBuilderAPI_MakeFace(aTrimmedSurface, Precision::Confusion()).Face();
testFace(aFace);
aFace = BRepBuilderAPI_MakeFace(aBSplineSurface, Precision::Confusion()).Face();
testFace(aFace);
}
void testSphere()
{
std::cout << "====== Test Sphere Area =====" << std::endl;
Handle(Geom_SphericalSurface) aSphericalSurface = new Geom_SphericalSurface(gp::XOY(), 1.0);
math_Vector aLower(1, 2);
math_Vector aUpper(1, 2);
aSphericalSurface->Bounds(aLower(1), aUpper(1), aLower(2), aUpper(2));
evalArea(aSphericalSurface, aLower, aUpper);
// Convert to BSpline Surface.
Handle(Geom_BSplineSurface) aBSplineSurface = GeomConvert::SurfaceToBSplineSurface(aSphericalSurface);
evalArea(aBSplineSurface);
// Test Face.
TopoDS_Face aFace = BRepBuilderAPI_MakeFace(aSphericalSurface, Precision::Confusion()).Face();
testFace(aFace);
aFace = BRepBuilderAPI_MakeFace(aBSplineSurface, Precision::Confusion()).Face();
testFace(aFace);
}
void test()
{
testPlane();
testSphere();
testCylinder();
}
int main(int argc, char* argv[])
{
test();
return 0;
}
计算结果如下图所示：
上述代码计算了曲面的面积，再将曲面转换成B样条曲面，再使用算法计算面积。再将曲面和转换的B样条曲面生成拓朴面，利用OpenCASCADE中计算曲面面积功能进行对比。使用自定义函数math_AreaFunction利用多重积分类计算的结果与OpenCASCADE中计算曲面面积的值是一致的。当把曲面转换成B样条曲面后，OpenCASCADE计算的曲面面积偏大。
7.Conclusion
在学习《高等数学》的积分时，其主要的一个应用就是计算弧长、面积和体积等。学习高数抽象概念时，总会问学了高数有什么用？就从计算机图形方面来看，可以利用数学工具对任意曲线求弧长，对任意曲面计算面积等，更具一般性。
通过自定义被积函数再利用积分算法来计算任意曲面的面积，将理论与实践结合起来了。即将曲面的第一基本公式与具体的代码甚至可以利用OpenCASCADE生成对应的图形，这样抽象的理论就直观了，更便于理解相应的概念。
8.References
1.朱心雄. 自由曲线曲面造型技术. 科学出版社. 2000
2.陈维桓. 微分几何. 北京大学出版社. 2006
3.同济大学数学教研室. 高等数学. 高等教育出版社. 1996