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