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时：
用数组来直接计算：
2 {
3 C = a[n];
4
5 for ( int i = n - 1 ; i >= 0 ; i -- )
6 {
7 C = C * u0 + a[i];
8 }
9 }
10
在OpenSceneGraph中来计算：
2 {
3 int n = a -> size() - 1 ;
4
5 p = a -> at(n);
6
7 for ( int i = n - 1 ; i >= 0 ; i -- )
8 {
9 p = p * u + a -> at(i);
10 }
11 }
三、 幂基曲线的显示 Render Power Basis Curve
利用Horner方法可以计算[0,1]区间上相应曲线上的点，将这些点连成线就构成了幂基曲线的近似表示。在OpenSceneGraph中显示幂基曲线程序如下所示：
2 * Copyright (c) 2013 eryar All Rights Reserved.
3 *
4 * File : Main.cpp
5 * Author : eryar@163.com
6 * Date : 2013-05-06
7 * Version : 0.1
8 *
9 * Description : PLib means Polynomial functions library.
10 * The PLib in occ provides basic computation functions
11 * for polynomial functions.
12 *
13 */
14
15 // OpenSceneGraph
16 #include < osg / Vec3 >
17 #include < osg / Array >
18 #include < osg / Geode >
19 #include < osg / Group >
20 #include < osg / MatrixTransform >
21 #include < osgGA / StateSetManipulator >
22 #include < osgViewer / Viewer >
23 #include < osgViewer / ViewerEventHandlers >
24
25 #pragma comment(lib, " osgd.lib " )
26 #pragma comment(lib, " osgGAd.lib " )
27 #pragma comment(lib, " osgViewerd.lib " )
28
29 // OpenCascade
30 #include < PLib.hxx >
31 #include < TColgp_Array1OfPnt.hxx >
32 #include < TColStd_Array1OfReal.hxx >
33
34 #pragma comment(lib, " TKernel.lib " )
35 #pragma comment(lib, " TKMath.lib " )
36
37 /*
38 * @breif Compute point on power basis curve.
39 * @param [in] a:
40 * @param [in] x:
41 * @return: Point on power basis curve.
42 */
43 void Horner(osg::Vec3Array * a, double u, osg::Vec3 & p)
44 {
45 int n = a -> size() - 1 ;
46
47 if ( - 1 == n)
48 {
49 return ;
50 }
51
52 p = a -> at(n);
53
54 for ( int i = n - 1 ; i >= 0 ; i -- )
55 {
56 p = p * u + a -> at(i);
57 }
58 }
59
60 osg::Node * RenderPowerBasisCurve()
61 {
62 const int nStep = 100 ;
63 osg::Geode * curveNode = new osg::Geode();
64 osg::ref_ptr < osg::Geometry > curveGeom = new osg::Geometry();
65 osg::ref_ptr < osg::Vec3Array > curvePnts = new osg::Vec3Array();
66
67 // Test to compuate point on power basis curve.
68 osg::ref_ptr < osg::Vec3Array > ctrlPnts = new osg::Vec3Array;
69 ctrlPnts -> push_back(osg::Vec3( 0 , 0 , 6 ));
70 ctrlPnts -> push_back(osg::Vec3( 3 , 0 , 6 ));
71 ctrlPnts -> push_back(osg::Vec3( 6 , 0 , 3 ));
72 // ctrlPnts->push_back(osg::Vec3(6, 0, 0));
73
74 for ( int i = 0 ; i < nStep; i ++ )
75 {
76 osg::Vec3 pnt;
77 Horner(ctrlPnts, i * 1.0 / nStep, pnt);
78
79 curvePnts -> push_back(pnt);
80 }
81
82 curveGeom -> setVertexArray(curvePnts);
83 curveGeom -> addPrimitiveSet( new osg::DrawArrays(osg::PrimitiveSet::LINE_STRIP, 0 , curvePnts -> size()));
84
85 curveNode -> addDrawable(curveGeom);
86
87 return curveNode;
88 }
89
90 osg::Node * TestPLib( void )
91 {
92 const int nStep = 100 ;
93 osg::Geode * curveNode = new osg::Geode();
94 osg::ref_ptr < osg::Geometry > curveGeom = new osg::Geometry();
95 osg::ref_ptr < osg::Vec3Array > curvePnts = new osg::Vec3Array();
96
97 TColgp_Array1OfPnt poles( 1 , 3 );
98 TColStd_Array1OfReal fp( 1 , poles.Length() * 3 );
99 TColStd_Array1OfReal points( 0 , 1 , 3 );
100
101 Standard_Real * polynomialCoeff = (Standard_Real * ) & (fp(fp.Lower()));
102 Standard_Real * result = (Standard_Real * ) & (points(points.Lower()));
103
104 poles.SetValue( 1 , gp_Pnt( 0 , 0 , 6 ));
105 poles.SetValue( 2 , gp_Pnt( 3 , 0 , 6 ));
106 poles.SetValue( 3 , gp_Pnt( 6 , 0 , 3 ));
107
108 // Change poles to flat array.
109 PLib::SetPoles(poles, fp);
110
111 // Three control point, so degree is 3-1=2.
112 Standard_Integer nDegree = 3 - 1 ;
113
114 // Because point are 3 Dimension.
115 Standard_Integer nDimension = 3 ;
116
117 for ( int i = 0 ; i < nStep; i ++ )
118 {
119 PLib::NoDerivativeEvalPolynomial(
120 i * 1.0 / nStep,
121 nDegree,
122 nDimension,
123 nDegree * nDimension,
124 polynomialCoeff[ 0 ],
125 result[ 0 ]);
126
127 //
128 curvePnts -> push_back(osg::Vec3(result[ 0 ], result[ 1 ], result[ 2 ]));
129 }
130
131 curveGeom -> setVertexArray(curvePnts);
132 curveGeom -> addPrimitiveSet( new osg::DrawArrays(osg::PrimitiveSet::LINE_STRIP, 0 , curvePnts -> size()));
133
134 curveNode -> addDrawable(curveGeom);
135
136 return curveNode;
137 }
138
139 int main( int argc, char * argv[])
140 {
141 osgViewer::Viewer myViewer;
142 osg::ref_ptr < osg::Group > root = new osg::Group();
143
144 root -> addChild(RenderPowerBasisCurve());
145
146 // Translate along x axis.
147 osg::ref_ptr < osg::MatrixTransform > transform = new osg::MatrixTransform();
148 transform -> setMatrix(osg::Matrix::translate( 1 , 0 , 0 ));
149 transform -> addChild(TestPLib());
150
151 root -> addChild(transform);
152
153 myViewer.setSceneData(root);
154
155 myViewer.addEventHandler( new osgGA::StateSetManipulator(myViewer.getCamera() -> getOrCreateStateSet()));
156 myViewer.addEventHandler( new osgViewer::StatsHandler);
157 myViewer.addEventHandler( new osgViewer::WindowSizeHandler);
158
159 return myViewer.run();
160 }
161
设置不同的控制点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的方法。程序代码如下所示：
2 {
3 const int nStep = 100 ;
4 osg::Geode * curveNode = new osg::Geode();
5 osg::ref_ptr < osg::Geometry > curveGeom = new osg::Geometry();
6 osg::ref_ptr < osg::Vec3Array > curvePnts = new osg::Vec3Array();
7
8 TColgp_Array1OfPnt poles( 1 , 3 );
9 TColStd_Array1OfReal fp( 1 , poles.Length() * 3 );
10 TColStd_Array1OfReal points( 0 , 1 , 3 );
11
12 Standard_Real * polynomialCoeff = (Standard_Real * ) & (fp(fp.Lower()));
13 Standard_Real * result = (Standard_Real * ) & (points(points.Lower()));
14
15 poles.SetValue( 1 , gp_Pnt( 0 , 0 , 6 ));
16 poles.SetValue( 2 , gp_Pnt( 3 , 0 , 6 ));
17 poles.SetValue( 3 , gp_Pnt( 6 , 0 , 3 ));
18
19 // Change poles to flat array.
20 PLib::SetPoles(poles, fp);
21
22 // Three control point, so degree is 3-1=2.
23 Standard_Integer nDegree = 3 - 1 ;
24
25 // Because point are 3 Dimension.
26 Standard_Integer nDimension = 3 ;
27
28 for ( int i = 0 ; i < nStep; i ++ )
29 {
30 PLib::NoDerivativeEvalPolynomial(
31 i * 1.0 / nStep,
32 nDegree,
33 nDimension,
34 nDegree * nDimension,
35 polynomialCoeff[ 0 ],
36 result[ 0 ]);
37
38 //
39 curvePnts -> push_back(osg::Vec3(result[ 0 ], result[ 1 ], result[ 2 ]));
40 }
41
42 curveGeom -> setVertexArray(curvePnts);
43 curveGeom -> addPrimitiveSet( new osg::DrawArrays(osg::PrimitiveSet::LINE_STRIP, 0 , curvePnts -> size()));
44
45 curveNode -> addDrawable(curveGeom);
46
47 return curveNode;
48 }
函数PLib::SetPoles可以将空间坐标点转换为实数数组。在调用无微分计算多项式的函数时，将坐标点的实数数组的首地址作为参数之一传入。
为了与前面的Horner方法计算多项式的结果进行对比，将OpenCascade对相同点计算的结果也显示出来。如下图5.1所示：
Figure 5.1 PLib compute result VS. Previous Horner method
由上图可知，PLib的计算结果与前面的Horner方法计算结果相同。查看OpenCascade的源程序，得其多项式计算方法也是采用的Horner方法。
2
3 const Standard_Integer Degree,
4
5 const Standard_Integer Dimension,
6
7 const Standard_Integer DegreeDimension,
8
9 Standard_Real & PolynomialCoeff,
10
11 Standard_Real & Results)
12
13 {
14
15 Standard_Integer jj;
16
17 Standard_Real * RA = & Results ;
18
19 Standard_Real * PA = & PolynomialCoeff ;
20
21 Standard_Real * tmpRA = RA;
22
23 Standard_Real * tmpPA = PA + DegreeDimension;
24
25 switch (Dimension) {
26
27 case 1 : {
28
29 * tmpRA = * tmpPA;
30
31 for (jj = Degree ; jj > 0 ; jj -- ) {
32
33 tmpPA -- ;
34
35 * tmpRA = Par * ( * tmpRA) + ( * tmpPA);
36
37 }
38
39 break ;
40
41 }
42
43 }
44
45
从上述计算一维多项式的代码可以看出，计算方法与前面的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
