OpenCASCADE 參數(shù)曲面面積
eryar@163.com
Abstract. 本文介紹了參數(shù)曲面的第一基本公式，并應(yīng)用曲面的第一基本公式，結(jié)合OpenCASCADE中計(jì)算多重積分的類，對任意參數(shù)曲面的面積進(jìn)行計(jì)算。
Key Words. Parametric Curve, Parametric Surface, Gauss Integration, Global Properties
1.Introduction
我們知道一元函數(shù)y=f(x)的圖像是一條曲線，二元函數(shù)z=f(x,y)的圖像是一張曲面。但是，把曲線曲面表示成參數(shù)方程則更加便利于研究，這種表示方法首先是由歐洲瑞士數(shù)學(xué)家Euler引進(jìn)的。例如，在空間中的一條曲線可以表示為三個(gè)一元函數(shù)：
X=x(t), Y=y(t), Z=z(t)
在向量的概念出現(xiàn)后，空間中的一條曲線可以自然地表示為一個(gè)一元向量函數(shù)：
r=r(t)=(x(t), y(t), z(t))
用向量函數(shù)來表示曲線和曲面后，使曲線曲面一些量的計(jì)算方式比較統(tǒng)一。如曲線可以表示為一元向量函數(shù)，曲面可以表示為二元向量函數(shù)。
本文結(jié)合OpenCASCADE來介紹參數(shù)曲線曲面積分分別計(jì)算曲線弧長和曲面的面積。結(jié)合《微分幾何》來更好地理解曲線曲面相關(guān)知識。
2.Curve Natural Parametric Equations
設(shè)曲線C的參數(shù)方程是r=r(t)，命：
則s是該曲線的一個(gè)不變量，即它與空間中的坐標(biāo)系的選擇無關(guān)，也與該曲線的參數(shù)變換無關(guān)。前者是因?yàn)樵诘芽栔苯亲鴺?biāo)變換下，切向量的長度|r’(t)|是不變的，故s不變。關(guān)于后者可以通過積分的變量的替換來證明，設(shè)參數(shù)變換是：
并且
因此
根據(jù)積分的變量替換公式有：
不變量s的幾何意義就是曲線段的弧長。這說明曲線參數(shù)t可以是任意的，但選擇不同的參數(shù)得到的參數(shù)方程會有不同，但是曲線段的弧長是不變的。以曲線弧長作為曲線方程的參數(shù)，這樣的方程稱為曲線的自然參數(shù)方程N(yùn)atural Parametric Equations。
由曲線的參數(shù)方程可知，曲線弧長的計(jì)算公式為：
幾何意義就是在每個(gè)微元處的切向量的長度求和。
3.Gauss Integration for Arc Length
曲線弧長的計(jì)算就是一元函數(shù)的積分。OpenCASCADE中是如何計(jì)算任意曲線弧長的呢？直接找到相關(guān)的源碼列舉如下：(在類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());
}
上述代碼的意思是直接對曲線的一階導(dǎo)數(shù)的長度求積分，即是弧長。OpenCASCADE的代碼寫得有點(diǎn)難懂，根據(jù)意思把對三維曲線求弧長的代碼改寫下，更便于理解：
//! 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
曲面參數(shù)方程是個(gè)二元向量函數(shù)。根據(jù)《微分幾何》中曲面的第一基本公式（First Fundamental Form of a Surface）可知，曲面上曲線的表達(dá)式為：
r=r(u(t), v(t)) = (x(t), y(t), z(t))
若以s表示曲面上曲線的弧長，則由復(fù)合函數(shù)求導(dǎo)公式可得弧長微分公式：
在古典微分幾何中，上式稱為曲面的第一基本公式，E、F、G稱為第一基本量。在曲面上，每一點(diǎn)的第一基本量與參數(shù)化無關(guān)。
利用曲面第一基本公式可以用于計(jì)算曲面的面積。參數(shù)曲面上與u,v參數(shù)空間的元素dudv對應(yīng)的面積元為：
由參數(shù)曲面法向的計(jì)算可知，曲面的面積元素即為u,v方向上的偏導(dǎo)數(shù)的乘積的模。
其幾何意義可以理解為參數(shù)曲面的面積微元是由u,v方向的偏導(dǎo)數(shù)的向量圍成的一個(gè)四邊形的面積，則整個(gè)曲面的面積即是對面積元素求積分。由于參數(shù)曲面有兩個(gè)參數(shù)，所以若要計(jì)算曲面的面積，只需要對面積元素計(jì)算二重積分即可。
5.Gauss Integration for Area
OpenCASCADE的math包中提供了多重積分的計(jì)算類math_GaussMultipleIntegration，由類名可知積分算法采用了Gauss積分算法。下面通過具體代碼來說明OpenCASCADE中計(jì)算曲面積分的過程。要計(jì)算積分，先要定義被積函數(shù)。因?yàn)閰?shù)曲面與參數(shù)曲線不同，參數(shù)曲線只有一個(gè)參數(shù)，而參數(shù)曲面有兩個(gè)參數(shù)，所以是一個(gè)多元函數(shù)。
//! 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;
};
由于參數(shù)曲面是多元函數(shù)，所以從類math_MultipleVarFunction派生，并在虛函數(shù)NbVariables()中說明有兩個(gè)變量。在虛函數(shù)Value()中計(jì)算面積元素的值，即根據(jù)曲面第一基本公式中面積元素的定義，對參數(shù)曲面求一階導(dǎo)數(shù)，計(jì)算兩個(gè)偏導(dǎo)數(shù)向量的叉乘的模。
有了被積函數(shù)，只需要在定義域?qū)ζ溆?jì)算二重積分，相應(yīng)代碼如下所示：
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算法計(jì)算二重積分，所以需要指定階數(shù)，且階數(shù)越高積分結(jié)果精度越高，這里使用了OpenCASCADE中最高的階數(shù)。
下面通過對基本曲面的面積計(jì)算來驗(yàn)證結(jié)果的正確性，并將計(jì)算結(jié)果和OpenCASCADE中計(jì)算面積的類BRepGProp::SurfaceProperties()結(jié)果進(jìn)行對比。
6.Elementary Surface Area Test
下面通過對OpenCASCADE中幾個(gè)初等曲面的面積進(jìn)行計(jì)算，代碼如下所示：
/*
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;
}
計(jì)算結(jié)果如下圖所示：
上述代碼計(jì)算了曲面的面積，再將曲面轉(zhuǎn)換成B樣條曲面，再使用算法計(jì)算面積。再將曲面和轉(zhuǎn)換的B樣條曲面生成拓樸面，利用OpenCASCADE中計(jì)算曲面面積功能進(jìn)行對比。使用自定義函數(shù)math_AreaFunction利用多重積分類計(jì)算的結(jié)果與OpenCASCADE中計(jì)算曲面面積的值是一致的。當(dāng)把曲面轉(zhuǎn)換成B樣條曲面后，OpenCASCADE計(jì)算的曲面面積偏大。
7.Conclusion
在學(xué)習(xí)《高等數(shù)學(xué)》的積分時(shí)，其主要的一個(gè)應(yīng)用就是計(jì)算弧長、面積和體積等。學(xué)習(xí)高數(shù)抽象概念時(shí)，總會問學(xué)了高數(shù)有什么用？就從計(jì)算機(jī)圖形方面來看，可以利用數(shù)學(xué)工具對任意曲線求弧長，對任意曲面計(jì)算面積等，更具一般性。
通過自定義被積函數(shù)再利用積分算法來計(jì)算任意曲面的面積，將理論與實(shí)踐結(jié)合起來了。即將曲面的第一基本公式與具體的代碼甚至可以利用OpenCASCADE生成對應(yīng)的圖形，這樣抽象的理論就直觀了，更便于理解相應(yīng)的概念。
8.References
1.朱心雄. 自由曲線曲面造型技術(shù). 科學(xué)出版社. 2000
2.陳維桓. 微分幾何. 北京大學(xué)出版社. 2006
3.同濟(jì)大學(xué)數(shù)學(xué)教研室. 高等數(shù)學(xué). 高等教育出版社. 1996