Polynomial Library in OpenCascade
摘要Abstract：分析冪基曲線即多項式曲線在OpenCascade中的計算方法，以及利用OpenSceneGraph來顯示OpenCascade的計算結果，加深對多項式曲線的理解。?
關鍵字Key Words：OpenCascade、PLib、OpenSceneGraph、Polynomial Library?
??
一、 概述 Overview
CAGD（Computer Aided Geometry Design）中，基表示的參數(shù)矢函數(shù)形式已經成為形狀數(shù)學描述的標準形式。人們首先注意到在各類函數(shù)中，多項式函數(shù)（Polynomial Function）能較好地滿足要求。它表示形式簡單，又無窮次可微，且容易計算函數(shù)值及各階導數(shù)值。采用多項式函數(shù)作為基函數(shù)即多項式基，相應可以得到參數(shù)多項式曲線曲面。?
人們很早就注意到這樣一個問題：設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]上的任何連續(xù)函數(shù)都可以用多項式來一致逼近。該定理實際上正好解決了利用多項式組成的函數(shù)來表示連續(xù)函數(shù)的問題。?
冪（又稱單項式monomial）基uj(j=0,1,,,n)是最簡單的多項式基。相應多項式的全體構成n次多項式空間。n次多項式空間中任一組n+1個線性無關的多項式都可以作為一組基，因此就有無窮多組基。不同組基之間僅僅相差一個線性變換。?
如果允許坐標函數(shù)x(u)，y(u)，z(u)是任意的，則可以得到范圍很廣的曲線。但是在實際開發(fā)一個幾何造型系統(tǒng)時，我們需要一些折中，理想的情況是將坐標函數(shù)限制在滿足下述條件的一類函數(shù)中：?
l 能夠精確地表示用戶需要的所有曲線；?
l 在計算機中能夠被方便、高效、精確地處理，特別是可以高效地計算曲線上的點及各階導矢；函數(shù)的數(shù)值計算對浮點數(shù)舍入誤差不敏感；函數(shù)所需要的存儲量較小；?
l 比較簡單，在數(shù)學上易于理解。?
一類被廣泛使用的函數(shù)就是多項式。盡管它們滿足上標準中的后兩項，但有很多類型的重要曲線（及曲面）不能用多項式精確表示，在系統(tǒng)中，這些曲線只能用多項式逼近。同一參數(shù)多項式曲線可以采用不同的基表示，如冪基表示和Bezier表示，由些決定了它們具有不同的性質，因而就有不同的優(yōu)點。?
?
二、 冪基曲線的計算 Calculate Polynomial Function
一條n次曲線的冪基表示形式是：?
，其中
是矢量。?
給定u0，計算冪基曲線上的點C（u0）的最有效算法是英國數(shù)學家W.G.Horner提出的Horner方法。Horner算法是遞歸概念的一個典型實例，它采用最少的乘法來進行多項式求值，使計算由X^n問題轉化為O（n）的問題。?
l 當次數(shù)＝1時：
l 當次數(shù)＝2時：
l ……?
l 當次數(shù)＝n時：
用數(shù)組來直接計算：
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]區(qū)間上相應曲線上的點，將這些點連成線就構成了冪基曲線的近似表示。在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）的數(shù)學計算工具箱（TKMath Toolkit）中有個PLib包，用以對多項式進行基本的計算。PLib庫中的函數(shù)都是靜態(tài)函數(shù)，所以都是類函數(shù)，可以用類名加函數(shù)名直接調用。?
PLib可對多項式進行如下計算：?
l 計算多項式的值：EvalPolynomial；?
l 計算Lagrange插值：EvalLagrange；?
l 計算Hermite插值：EvalCubicHermite；?
其中計算多項式值的方法也是用的Horner方法。?
包PLib中提供了計算幾何的數(shù)學基礎中多項式插值中大部分插值計算。結合書籍《計算幾何教程》科學出版社中第一章的理論內容及OpenCascade的源程序，可以方便計算幾何的數(shù)學基礎知識的學習。?
??
五、 使用PLib Apply PLib Class
因為包PLib中的類PLib都是靜態(tài)函數(shù)，所以函數(shù)傳入的參數(shù)比較多，若要使用這些計算函數(shù)，需要對其函數(shù)參數(shù)進行了解。為了對不同維數(shù)多項式進行計算，類PLib中把空間點轉換成了實數(shù)數(shù)組，并提供了相互轉換的函數(shù)。以計算多項式值為例，來說明使用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;
}
函數(shù)PLib::SetPoles可以將空間坐標點轉換為實數(shù)數(shù)組。在調用無微分計算多項式的函數(shù)時，將坐標點的實數(shù)數(shù)組的首地址作為參數(shù)之一傳入。?
為了與前面的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為了統(tǒng)一多項式的計算，將空間點都轉換成數(shù)組后再進行計算，在這其中大量使用了指針，代碼可讀性也不是很好，需要仔細、耐心。?
??
七、 參考資料?
1. 王仁宏、李崇君、朱春鋼 計算幾何教程 科學出版社 2008.6?
2. 趙罡、穆國旺、王拉柱 非均勻有理B樣條《The NURBS Book》 清華大學出版社 2010.12?
3.? OpenCascade source code?
?
PDF Version and Source Code: Polynomial Library in OpenCascade
更多文章、技術交流、商務合作、聯(lián)系博主
微信掃碼或搜索：z360901061
微信掃一掃加我為好友
QQ號聯(lián)系： 360901061
您的支持是博主寫作最大的動力，如果您喜歡我的文章，感覺我的文章對您有幫助，請用微信掃描下面二維碼支持博主2元、5元、10元、20元等您想捐的金額吧，狠狠點擊下面給點支持吧，站長非常感激您！手機微信長按不能支付解決辦法：請將微信支付二維碼保存到相冊，切換到微信，然后點擊微信右上角掃一掃功能，選擇支付二維碼完成支付。
【本文對您有幫助就好】元
