Polynomial Library in OpenCascade
摘要Abstract：分析冪基曲線即多項式曲線在OpenCascade中的計算方法，以及利用OpenSceneGraph來顯示OpenCascade的計算結果，加深對多項式曲線的理解。?
關鍵字Key Words：OpenCascade、PLib、OpenSceneGraph、Polynomial Library?
??
一、 概述 Overview
CAGD（Computer Aided Geometry Design）中，基表示的參數矢函數形式已經成為形狀數學描述的標準形式。人們首先注意到在各類函數中，多項式函數（Polynomial Function）能較好地滿足要求。它表示形式簡單，又無窮次可微，且容易計算函數值及各階導數值。采用多項式函數作為基函數即多項式基，相應可以得到參數多項式曲線曲面。?
人們很早就注意到這樣一個問題：設f(x)∈C[a,b]，對于任意給的ε>0，是否存在這樣的多項式p(x)，使得?
對一切a≤x≤b一致成立？?
1885年，Weierstrass證明了p(x)的存在性；?
1912年，Bernstein將滿足條件的多項式構造出來。?
Weierstrass定理：設f(x)∈C[a,b]，那么對于任意給的ε>0，都存在這樣的多項式p(x)，使得?
Weierstrass定理說明，[a,b]上的任何連續函數都可以用多項式來一致逼近。該定理實際上正好解決了利用多項式組成的函數來表示連續函數的問題。?
冪（又稱單項式monomial）基uj(j=0,1,,,n)是最簡單的多項式基。相應多項式的全體構成n次多項式空間。n次多項式空間中任一組n+1個線性無關的多項式都可以作為一組基，因此就有無窮多組基。不同組基之間僅僅相差一個線性變換。?
如果允許坐標函數x(u)，y(u)，z(u)是任意的，則可以得到范圍很廣的曲線。但是在實際開發一個幾何造型系統時，我們需要一些折中，理想的情況是將坐標函數限制在滿足下述條件的一類函數中：?
l 能夠精確地表示用戶需要的所有曲線；?
l 在計算機中能夠被方便、高效、精確地處理，特別是可以高效地計算曲線上的點及各階導矢；函數的數值計算對浮點數舍入誤差不敏感；函數所需要的存儲量較小；?
l 比較簡單，在數學上易于理解。?
一類被廣泛使用的函數就是多項式。盡管它們滿足上標準中的后兩項，但有很多類型的重要曲線（及曲面）不能用多項式精確表示，在系統中，這些曲線只能用多項式逼近。同一參數多項式曲線可以采用不同的基表示，如冪基表示和Bezier表示，由些決定了它們具有不同的性質，因而就有不同的優點。?
?
二、 冪基曲線的計算 Calculate Polynomial Function
一條n次曲線的冪基表示形式是：?
，其中
是矢量。?
給定u0，計算冪基曲線上的點C（u0）的最有效算法是英國數學家W.G.Horner提出的Horner方法。Horner算法是遞歸概念的一個典型實例，它采用最少的乘法來進行多項式求值，使計算由X^n問題轉化為O（n）的問題。?
l 當次數＝1時：
l 當次數＝2時：
l ……?
l 當次數＝n時：
用數組來直接計算：
void
Horner1(a, n, u0, C)
{
C
=
a[n];
for
(
int
i = n-
1
; i >=
0
; i--
)
{
C
= C * u0 +
a[i];
}
}
?在OpenSceneGraph中來計算：??
void
Horner(osg::Vec3Array* a,
double
u, osg::Vec3&
p)
{
int
n = a->size() -
1
;
p
= a->
at(n);
for
(
int
i = n-
1
; i >=
0
; i--
)
{
p
= p * u + a->
at(i);
}
}
?
三、 冪基曲線的顯示 Render Power Basis Curve
利用Horner方法可以計算[0,1]區間上相應曲線上的點，將這些點連成線就構成了冪基曲線的近似表示。在OpenSceneGraph中顯示冪基曲線程序如下所示：?
#include <osg/Vec3>
#include
<osg/Array>
#include
<osg/Geode>
#include
<osg/Group>
#include
<osgGA/StateSetManipulator>
#include
<osgViewer/Viewer>
#include
<osgViewer/ViewerEventHandlers>
#pragma
comment(lib, "osgd.lib")
#pragma
comment(lib, "osgGAd.lib")
#pragma
comment(lib, "osgViewerd.lib")
/*
* @breif Compute point on power basis curve.
* @param [in] a:
* @param [in] x:
* @return: Point on power basis curve.
*/
void
Horner(osg::Vec3Array* a,
double
u, osg::Vec3&
p)
{
int
n = a->size() -
1
;
if
(-
1
==
n)
{
return
;
}
p
= a->
at(n);
for
(
int
i = n-
1
; i >=
0
; i--
)
{
p
= p * u + a->
at(i);
}
}
osg::Node
*
RenderPowerBasisCurve()
{
const
int
nStep =
100
;
osg::Geode
* curveNode =
new
osg::Geode();
osg::ref_ptr
<osg::Geometry> curveGeom =
new
osg::Geometry();
osg::ref_ptr
<osg::Vec3Array> curvePnts =
new
osg::Vec3Array();
//
Test to compuate point on power basis curve.
osg::ref_ptr<osg::Vec3Array> ctrlPnts =
new
osg::Vec3Array;
ctrlPnts
->push_back(osg::Vec3(
0
,
0
,
6
));
ctrlPnts
->push_back(osg::Vec3(
3
,
0
,
6
));
ctrlPnts
->push_back(osg::Vec3(
6
,
0
,
3
));
ctrlPnts
->push_back(osg::Vec3(
6
,
0
,
0
));
for
(
int
i =
0
; i < nStep; i++
)
{
osg::Vec3 pnt;
Horner(ctrlPnts, i
*
1.0
/
nStep, pnt);
curvePnts
->
push_back(pnt);
}
curveGeom
->
setVertexArray(curvePnts);
curveGeom
->addPrimitiveSet(
new
osg::DrawArrays(osg::PrimitiveSet::LINE_STRIP,
0
, curvePnts->
size()));
curveNode
->
addDrawable(curveGeom);
return
curveNode;
}
int
main(
int
argc,
char
*
argv[])
{
osgViewer::Viewer myViewer;
osg::ref_ptr
<osg::Group> root =
new
osg::Group();
root
->
addChild(RenderPowerBasisCurve());
myViewer.setSceneData(root);
myViewer.addEventHandler(
new
osgGA::StateSetManipulator(myViewer.getCamera()->
getOrCreateStateSet()));
myViewer.addEventHandler(
new
osgViewer::StatsHandler);
myViewer.addEventHandler(
new
osgViewer::WindowSizeHandler);
return
myViewer.run();
}
設置不同的控制點ctrlPnts，就得到不同的曲線。?
當n=1時，有兩個控制點a0, a1，表示由a0到a0+a1的直線段，如圖3.1所示：?
Figure 3.1 連接兩點（0,0,6）到（3,0,6）的直線?
當n=2時，曲線是一段由a0到a0+a1+a2的拋物線弧，如圖3.2所示：?
Figure 3.2 拋物線弧（0,0,6）（3,0,6）（6,0,3）?
??
四、 occ中的多項式計算庫The PLib in OCC
在OpenCascade中的基礎模塊（FoundationClasses）的數學計算工具箱（TKMath Toolkit）中有個PLib包，用以對多項式進行基本的計算。PLib庫中的函數都是靜態函數，所以都是類函數，可以用類名加函數名直接調用。?
PLib可對多項式進行如下計算：?
l 計算多項式的值：EvalPolynomial；?
l 計算Lagrange插值：EvalLagrange；?
l 計算Hermite插值：EvalCubicHermite；?
其中計算多項式值的方法也是用的Horner方法。?
包PLib中提供了計算幾何的數學基礎中多項式插值中大部分插值計算。結合書籍《計算幾何教程》科學出版社中第一章的理論內容及OpenCascade的源程序，可以方便計算幾何的數學基礎知識的學習。?
??
五、 使用PLib Apply PLib Class
因為包PLib中的類PLib都是靜態函數，所以函數傳入的參數比較多，若要使用這些計算函數，需要對其函數參數進行了解。為了對不同維數多項式進行計算，類PLib中把空間點轉換成了實數數組，并提供了相互轉換的函數。以計算多項式值為例，來說明使用PLib的方法。程序代碼如下所示：?
osg::Node* TestPLib(
void
)
{
const
int
nStep =
100
;
osg::Geode
* curveNode =
new
osg::Geode();
osg::ref_ptr
<osg::Geometry> curveGeom =
new
osg::Geometry();
osg::ref_ptr
<osg::Vec3Array> curvePnts =
new
osg::Vec3Array();
TColgp_Array1OfPnt poles(
1
,
3
);
TColStd_Array1OfReal fp(
1
, poles.Length() *
3
);
TColStd_Array1OfReal points(
0
,
1
,
3
);
Standard_Real
* polynomialCoeff = (Standard_Real*) &
(fp(fp.Lower()));
Standard_Real
* result = (Standard_Real*)&
(points(points.Lower()));
poles.SetValue(
1
, gp_Pnt(
0
,
0
,
6
));
poles.SetValue(
2
, gp_Pnt(
3
,
0
,
6
));
poles.SetValue(
3
, gp_Pnt(
6
,
0
,
3
));
//
Change poles to flat array.
PLib::SetPoles(poles, fp);
//
Three control point, so degree is 3-1=2.
Standard_Integer nDegree =
3
-
1
;
//
Because point are 3 Dimension.
Standard_Integer nDimension =
3
;
for
(
int
i =
0
; i < nStep; i++
)
{
PLib::NoDerivativeEvalPolynomial(
i
*
1.0
/
nStep,
nDegree,
nDimension,
nDegree
*
nDimension,
polynomialCoeff[
0
],
result[
0
]);
//
curvePnts->push_back(osg::Vec3(result[
0
], result[
1
], result[
2
]));
}
curveGeom
->
setVertexArray(curvePnts);
curveGeom
->addPrimitiveSet(
new
osg::DrawArrays(osg::PrimitiveSet::LINE_STRIP,
0
, curvePnts->
size()));
curveNode
->
addDrawable(curveGeom);
return
curveNode;
}
函數PLib::SetPoles可以將空間坐標點轉換為實數數組。在調用無微分計算多項式的函數時，將坐標點的實數數組的首地址作為參數之一傳入。?
為了與前面的Horner方法計算多項式的結果進行對比，將OpenCascade對相同點計算的結果也顯示出來。如下圖5.1所示：?
Figure 5.1 PLib compute result VS. Previous Horner method?
由上圖可知，PLib的計算結果與前面的Horner方法計算結果相同。查看OpenCascade的源程序，得其多項式計算方法也是采用的Horner方法。?
void
PLib::NoDerivativeEvalPolynomial(
const
Standard_Real Par,
const
Standard_Integer Degree,
const
Standard_Integer Dimension,
const
Standard_Integer DegreeDimension,
Standard_Real
&
PolynomialCoeff,
Standard_Real
&
Results)
{
Standard_Integer jj;
Standard_Real
*RA = &
Results ;
Standard_Real
*PA = &
PolynomialCoeff ;
Standard_Real
*tmpRA =
RA;
Standard_Real
*tmpPA = PA +
DegreeDimension;
switch
(Dimension) {
case
1
: {
*tmpRA = *
tmpPA;
for
(jj = Degree ; jj >
0
; jj--
) {
tmpPA
--
;
*tmpRA = Par * (*tmpRA) + (*
tmpPA);
}
break
;
}
}
從上述計算一維多項式的代碼可以看出，計算方法與前面的Horner方法相同。?
??
六、 結論?
學習使用Horner方法來計算多項式的值，并將計算結果在OpenSceneGraph中顯示。通過使用OpenCascade的多項式庫PLib來計算多項式的值，并結合其源程序來理解如何使用庫PLib。庫PLib為了統一多項式的計算，將空間點都轉換成數組后再進行計算，在這其中大量使用了指針，代碼可讀性也不是很好，需要仔細、耐心。?
??
七、 參考資料?
1. 王仁宏、李崇君、朱春鋼 計算幾何教程 科學出版社 2008.6?
2. 趙罡、穆國旺、王拉柱 非均勻有理B樣條《The NURBS Book》 清華大學出版社 2010.12?
3.? OpenCascade source code?
?
PDF Version and Source Code: Polynomial Library in OpenCascade
更多文章、技術交流、商務合作、聯系博主
微信掃碼或搜索：z360901061
微信掃一掃加我為好友
QQ號聯系： 360901061
您的支持是博主寫作最大的動力，如果您喜歡我的文章，感覺我的文章對您有幫助，請用微信掃描下面二維碼支持博主2元、5元、10元、20元等您想捐的金額吧，狠狠點擊下面給點支持吧，站長非常感激您！手機微信長按不能支付解決辦法：請將微信支付二維碼保存到相冊，切換到微信，然后點擊微信右上角掃一掃功能，選擇支付二維碼完成支付。
【本文對您有幫助就好】元
