B-Spline Curve Library in Open Cascade
Open Cascade中的B樣條曲線庫
eryar@163.com
摘要Abstract：簡要介紹Open Cascade中的B樣條曲線庫BSplCLib的使用方法，并且結合源程序來對Open Cascade中的B樣條曲線的組成部分如節點矢量、重復度等概念進行介紹，以及通過對計算B樣條基函數的算法進行分析，加深對B樣條曲線概念的理解。
關鍵字Key Word：B Spline Curve、Open Cascade、Knot Vector、Multiplicity
一、 概述 Overview
1946年由Schoenberg提出了B樣條理論，給出了B樣條的差分表達式；1972年de Boor和Cox分別獨立給出了關于B樣條的標準算法。Gordon和Riesenfeld又把B樣條理論用于形狀描述，最終提出了B樣條方法。用B樣條基替代了Bernstein基，構造出B樣條曲線，這種方法繼承了Bezier方法的一切優點，克服了Bezier方法存在的缺點，較成功地解決了局部控制問題，又輕而易舉地在參數連續性基礎上解決了連接問題，從而使自由曲線曲面形狀的描述問題得到較好解決。
p次B樣條曲線的定義為：
其中：
l Pi是控制頂點（control point）；
l Ni,p(u)是定義在非周期節點矢量上的p次B樣條基函數；
有很多方法可以用來定義B樣條基函數以及證明它的一些重要性質。例如，可以采用截尾冪函數的差商定義，開花定義，以及由de Boor和Cox等人提出的遞推公式等來定義。我們這里采用的是遞推定義方法，因為這種方法在計算機實現中是最有效的。
令U={u0,u1,…,um}是一個單調不減的實數序列，即ui<=ui+1，i=0,1,…,m-1。其中，ui稱為節點，U稱為節點矢量，用Ni,p(u)表示第i個p次B樣條基函數，其定義為：
B樣條基有如下性質：
a) 遞推性；
b) 局部支承性；
c) 規范性；
d) 可微性；
根據B樣條曲線定義可知，給定控制頂點Pi（control points），曲線次數p（degree）及節點矢量U（knot vectors），B樣曲線也就確定。對于有理B樣條曲線，還需要參數權重（weights）。
二、 OCC中的B樣條曲線庫 BSplCLib in OCC
在Open Cascade中的工具箱（Toolkit）TKMath中的包（package）BSplCLib是B樣條曲線庫，為B樣條曲線曲面的計算提供了支持。它提供了三方面的功能：
l 對節點矢量（knot vectors）及重復度（multiplicities）的管理；
l 對多維樣條的支持，即B樣條方法中控制頂點的維數可以是任意維數（dimension）；
l 二維和三維樣條曲線的方法；
Open Cascade中的B樣條曲線由下列數據項定義：
|
定義 |
變量類型 |
變量名稱 |
|
控制頂點control points |
TColgp_Array1OfPnt |
Poles |
|
權重weights |
TColStd_Array1OfReal |
Weights |
|
節點knots |
TColStd_Array1OfReal |
Knots |
|
重數multiplicities |
TColStd_Array1OfInteger |
Mults |
|
次數degree |
Standard_Integer |
Degree |
|
周期性periodicity |
Standard_Boolean |
Periodic |
B樣條曲線庫BSplCLib提供了一些基本幾何算法：
l B樣條基函數及其導數的計算BSplCLib::EvalBsplineBasis()；
l 節點插入BSplCLib::InsertKnot()；
l 節點去除BSplCLib::RemoveKnot()；
l 升階BSplCLib::IncreaseDegree()；
l 降階；
結合《The NURBS Book》和Open Cascade中的BSplCLib的源程序，可以高效的學習NURBS。《The NURBS Book》中有詳細的理論推導及算法描述，而Open Cascade中有可以用來實際使用的程序。理論聯系實際，有助于快速理解NURBS的有關概念及其應用。
三、 OCC中B樣條曲線庫的節點和重數Knots and Multiplicity in BSplCLib
由B樣條曲線的可微性可知，節點的重數與B樣條曲線的連續性相關。在節點區間內部，Ni,p(u)是無限次可微的，因在每個節點區間內部，它是一個多項式。在節點處Ni,p(u)是p-k次連續的，其中k是節點的重復度（multiplicity，有時也稱為重數）。因此，增加次數p將提高曲線的連續性，而增加節點的重復度則使連續性降低。
重復度（multiplicity，有時也稱為重數）有兩種不同的理解方式：
l 節點在節點矢量中的重復度；
l 節點相對于一個特定的基函數的重復度；
在Open Cascade中對重復度的理解是前者，即節點在節點矢量中的重復度。下面結合源程序來進行說明。
函數BSplCLib::Knots()用來將給定的節點矢量（節點序列knot sequence）轉換為節點的重復度不大于1的Knots數組和每個節點對應的重復度Mults數組，且數據Knots和Mults的長度必由函數BSplCLib::KnotsLength()得到。Knots()函數的源程序如下所示：
//=======================================================================
//function : Knots
//purpose : Computes the sequence of knots Knots without repetition
// of the knots of multiplicity greater than 1.
// Length of <Knots> and <Mults> must be KnotsLength(KnotSequence,Periodic).
//=======================================================================
void BSplCLib::Knots(const TColStd_Array1OfReal& SeqKnots,
TColStd_Array1OfReal &knots,
TColStd_Array1OfInteger &mult,
// const Standard_Boolean Periodic)
const Standard_Boolean )
{
Standard_Real val = SeqKnots(1);
Standard_Integer kk=1;
knots(kk) = val;
mult(kk) = 1;
for (Standard_Integer jj=2;jj<=SeqKnots.Length();jj++)
{
// test on strict equality on nodes
if (SeqKnots(jj)!=val)
{
val = SeqKnots(jj);
kk++;
knots(kk) = val;
mult(kk) = 1;
}
else
{
mult(kk)++;
}
}
}
從上述代碼可知，直接使用了不等于來判斷兩個節點的值是否相同，而沒有采用誤差處理，即嚴格的相等比較。程序將節點重復度不大于1的節點及其相應的重復度分別保存到knots和mult中。
四、 B樣條曲線的分類 B Spline Curve Type
B樣條曲線一般按定義基函數的節點序列是否等距（均勻）分為均勻B樣條曲線（Uniform B-Spline Curve）和非均勻B樣條曲線（Non Uniform B-Spline Curve）。
B樣條曲線按節點序列中節點分布情況不同，又分為四種類型：均勻B樣條曲線、準均勻B樣條曲線、分段Bezier曲線、一般非均勻B樣條曲線。設給定特征多邊形頂點Vi，i=0,1,…,n，曲線次數k，則有：
l 均勻B樣條曲線（uniform B-Spline curve）：節點序列中節點沿參數軸均勻或等距分布，即所有節點區間長度為大于零的常數（constant）：
l 準均勻B樣曲線（quasi-uniform B-Spline curve）：其節點序列中兩端節點具有重復度k+1，而所有內節點均勻分布，具有重復度1。
l 分段Bezier曲線（piecewise Bezier curve）：其節點序列中兩端節點重復度與準均勻B樣條曲線的相同，所不同的是所有內節點重復度為k。
l 非均勻B樣條曲線（general non-uniform B-Spline curve）：這是對任意分布的節點序列，只要在數學上成立，即節點序列非遞減，都可取。
在基礎類模塊（Module FoundationClasses）的工具箱（Toolkit TKMath）中的包（GeomAbs）中有對B樣樣條曲線類型的定義，源程序如下所示：
//! This enumeration is used to note specific curve form. <br>
enum GeomAbs_BSplKnotDistribution {
GeomAbs_NonUniform,
GeomAbs_Uniform,
GeomAbs_QuasiUniform,
GeomAbs_PiecewiseBezier
};
而類BSplCLib主要是用來管理節點和重復度的，所有將節點和重復度也進行了分類。根據節點矢量是否均勻分布，將節點分配方式（Knot Distribution）分為：均勻（BSplCLib_Uniform）和非均勻（BSplCLib_NonUniform）。源程序如下所示：
enum BSplCLib_KnotDistribution {
BSplCLib_NonUniform,
BSplCLib_Uniform
};
根據重復度數組將重復度的分配方式分為如下三種類型：
n BSplCLib_Constant：重復度都相同；
n BSplCLib_QuasiConstant：首、尾節點的重復度與內部節點的重復度不同；
n BSplClib_NonConstant：其它情況；
源程序如下所示：
enum BSplCLib_MultDistribution {
BSplCLib_NonConstant,
BSplCLib_Constant,
BSplCLib_QuasiConstant
};
判斷節點矢量和重復度矢量類型分別由下列函數實現：
l BSplCLib::KnotForm()；
l BSplCLib::MultForm()；
具體的判斷方法可以查看源程序。
將節點分布方式與重復度的分布方式進行組合，可以得出B樣條曲線的那幾種類型。
五、 B樣條基函數的計算 Evaluate the B-Spline Basis
B樣條基函數的計算主要使用了B樣條基了函數的遞推公式（Cox-deBoor公式）的局部支撐性質，如下所示：
直接由定義可知：
l Ni,0(u)是一個階梯函數，它在半開區間u∈[ui, ui+1)外都為零；
l 當次數p>0時，Ni,p(u)是兩個p-1次基函數的線性組合；
l 計算一組基函數需要事先指定節點矢量U和次數p；
l 半開區間[ui,ui+1)稱為第i個節點區間（knot span），它的長度可以為零，因為相鄰節點可以是相同的；
l 計算p次基函數的過程可以生成一個如下形式的三角形陳列：
B樣條有局部支撐性，即若u不在區間[ui, ui+p+1)，則Ni,p(u)=0。可從下面的三角形中看出N1,3是N1,0、N2,0、N3,0和N4,0的線性組合，而N1,0在區間[u1, u2)上非零，N2,0在區間[u2,u3)上非零，N3,0在區間[u3,u4)上非零，N4,0在區間[u4,u5)上非零，所以N1,3僅在區間[u1,u5)上非零。
在任意給定的節點區間[uj,uj+1)內，最多有p+1個是非零的，它們是Nj-p,p、Nj-p+1、…、Nj,p。例如，在[u3,u4)上，零次基函數中只有N3,0是非零的，一次基函數只有N2,1和N3,1是非零的，非零的三次基函數只有N0,3、N1,3、N2,3、N3,3。這個性質如下圖所示：
上面兩幅圖中右邊的圖中所示的推算過程表明，給定節點序列U及B樣條曲線的次數p，給出任意一個u值，找出其所在的節點區間[ui,ui+1)上，最多有Ni-p,p，Ni-p+1,p，…，Ni,p個非零的基函數。
例如我們根據遞推公式寫出二次基函數的一般形式，如下所示：
當給定的u值在區間[u3,u4)上即（i=3）時，根據上面的三角形，得出下列重要結論：
即這兩項不需要計算。另外一個重要結論就是圖中用相同顏色框中的部分是相同的，也就是下面程序中的變量temp表示的內容。
我們引入下面符號：
由二次基函數推出的三個公式可寫為：
上述推導過程為《The NURBS Book》中的算法，算法代碼如下所示：
理解了變量temp的意義之后，整個程序就很好理解了。
將Open Cascade中計算基函數的算法是不同的，將其源程序摘抄如下所示：
//=======================================================================
//function : Build BSpline Matrix
//purpose : Builds the Bspline Matrix
//=======================================================================
Standard_Integer
BSplCLib::EvalBsplineBasis
//(const Standard_Integer Side, // = 1 rigth side, -1 left side
(const Standard_Integer , // = 1 rigth side, -1 left side
const Standard_Integer DerivativeRequest,
const Standard_Integer Order,
const TColStd_Array1OfReal& FlatKnots,
const Standard_Real Parameter,
Standard_Integer& FirstNonZeroBsplineIndex,
math_Matrix& BsplineBasis)
{
// the matrix must have at least DerivativeRequest + 1
// row and Order columns
// the result are stored in the following way in
// the Bspline matrix
// Let i be the FirstNonZeroBsplineIndex and
// t be the parameter value, k the order of the
// knot vector, r the DerivativeRequest :
//
// B (t) B (t) B (t)
// i i+1 i+k-1
//
// (1) (1) (1)
// B (t) B (t) B (t)
// i i+1 i+k-1
//
//
//
//
// (r) (r) (r)
// B (t) B (t) B (t)
// i i+1 i+k-1
//
Standard_Integer
ReturnCode,
ii,
pp,
qq,
ss,
NumPoles,
LocalRequest ;
// ,Index ;
Standard_Real NewParameter,
Inverse,
Factor,
LocalInverse,
Saved ;
// , *FlatKnotsArray ;
ReturnCode = 0 ;
FirstNonZeroBsplineIndex = 0 ;
LocalRequest = DerivativeRequest ;
if (DerivativeRequest >= Order) {
LocalRequest = Order - 1 ;
}
if (BsplineBasis.LowerCol() != 1 ||
BsplineBasis.UpperCol() < Order ||
BsplineBasis.LowerRow() != 1 ||
BsplineBasis.UpperRow() <= LocalRequest) {
ReturnCode = 1;
goto FINISH ;
}
NumPoles = FlatKnots.Upper() - FlatKnots.Lower() + 1 - Order ;
BSplCLib::LocateParameter(Order - 1,
FlatKnots,
Parameter,
Standard_False,
Order,
NumPoles+1,
ii,
NewParameter) ;
FirstNonZeroBsplineIndex = ii - Order + 1 ;
BsplineBasis(1,1) = 1.0e0 ;
LocalRequest = DerivativeRequest ;
if (DerivativeRequest >= Order) {
LocalRequest = Order - 1 ;
}
for (qq = 2 ; qq <= Order - LocalRequest ; qq++) {
BsplineBasis(1,qq) = 0.0e0 ;
for (pp = 1 ; pp <= qq - 1 ; pp++) {
//
// this should be always invertible if ii is correctly computed
//
Factor = (Parameter - FlatKnots(ii - qq + pp + 1))
/ (FlatKnots(ii + pp) - FlatKnots(ii - qq + pp + 1)) ;
Saved = Factor * BsplineBasis(1,pp) ;
BsplineBasis(1,pp) *= (1.0e0 - Factor) ;
BsplineBasis(1,pp) += BsplineBasis(1,qq) ;
BsplineBasis(1,qq) = Saved ;
}
}
for (qq = Order - LocalRequest + 1 ; qq <= Order ; qq++) {
for (pp = 1 ; pp <= qq - 1 ; pp++) {
BsplineBasis(Order - qq + 2,pp) = BsplineBasis(1,pp) ;
}
BsplineBasis(1,qq) = 0.0e0 ;
for (ss = Order - LocalRequest + 1 ; ss <= qq ; ss++) {
BsplineBasis(Order - ss + 2,qq) = 0.0e0 ;
}
for (pp = 1 ; pp <= qq - 1 ; pp++) {
Inverse = 1.0e0 / (FlatKnots(ii + pp) - FlatKnots(ii - qq + pp + 1)) ;
Factor = (Parameter - FlatKnots(ii - qq + pp + 1)) * Inverse ;
Saved = Factor * BsplineBasis(1,pp) ;
BsplineBasis(1,pp) *= (1.0e0 - Factor) ;
BsplineBasis(1,pp) += BsplineBasis(1,qq) ;
BsplineBasis(1,qq) = Saved ;
LocalInverse = (Standard_Real) (qq - 1) * Inverse ;
for (ss = Order - LocalRequest + 1 ; ss <= qq ; ss++) {
Saved = LocalInverse * BsplineBasis(Order - ss + 2, pp) ;
BsplineBasis(Order - ss + 2, pp) *= - LocalInverse ;
BsplineBasis(Order - ss + 2, pp) += BsplineBasis(Order - ss + 2,qq) ;
BsplineBasis(Order - ss + 2,qq) = Saved ;
}
}
}
FINISH :
return (ReturnCode) ;
}
函數的作用是用來計算所有的基函數及其導數，并將結果以矩陣（數組）的形式保存。結合二次基函數的推導方法，將述代碼寫成公式的形式。函數的參數及其描述如下表所示：
|
變量 |
描述 |
|
DerivativeRequest |
導數的次數 |
|
Order |
B樣條基函數的階數(次數+1) |
|
FlatKnots |
節點矢量 |
|
Parameter |
參數 |
|
FirstNonZeroBspline |
第一個非零基函數的索引值 |
|
BsplineBasis |
基函數值矩陣 |
當導數次數DerivativeRequest大于B樣條基的階數Order時，將計算導數的次數設置為B樣條基的次數（Order-1）。程序代碼如下所示：
LocalRequest = DerivativeRequest ;
if (DerivativeRequest >= Order) {
LocalRequest = Order - 1 ;
}
對B樣條基數計算結果矩陣BsplineBasis存儲空間進行檢查。若存儲空間不足，則會退出，程序代碼如下所示：
if (BsplineBasis.LowerCol() != 1 ||
BsplineBasis.UpperCol() < Order ||
BsplineBasis.LowerRow() != 1 ||
BsplineBasis.UpperRow() <= LocalRequest) {
ReturnCode = 1;
goto FINISH ;
}
確定參數Parameter所在的節點區間的下標（索引值），程序代碼如下所示：
NumPoles = FlatKnots.Upper() - FlatKnots.Lower() + 1 - Order ;
BSplCLib::LocateParameter(Order - 1,
FlatKnots,
Parameter,
Standard_False,
Order,
NumPoles+1,
ii,
NewParameter) ;
確定參數Parameter所在區間的算法是用二分法搜索得到。程序代碼如下所示：
//=======================================================================
//function : Hunt
//purpose :
//=======================================================================
void BSplCLib::Hunt (const Array1OfReal& XX,
const Standard_Real X,
Standard_Integer& Ilc)
{
// replaced by simple dichotomy (RLE)
Ilc = XX.Lower();
const Standard_Real *px = &XX(Ilc);
px -= Ilc;
if (X < px[Ilc]) {
Ilc--;
return;
}
Standard_Integer Ihi = XX.Upper();
if (X > px[Ihi]) {
Ilc = Ihi + 1;
return;
}
Standard_Integer Im;
while (Ihi - Ilc != 1) {
Im = (Ihi + Ilc) >> 1;
if (X > px[Im]) Ilc = Im;
else Ihi = Im;
}
}
確定參數所在區間[ui,ui+1)后，可得到第一個非零基函數的索引值為i-p；
FirstNonZeroBsplineIndex = ii - Order + 1 ;
基函數計算的主要算法代碼如下所示：
BsplineBasis(1,1) = 1.0e0 ;
for (qq = 2 ; qq <= Order - LocalRequest ; qq++) {
BsplineBasis(1,qq) = 0.0e0 ;
for (pp = 1 ; pp <= qq - 1 ; pp++) {
//
// this should be always invertible if ii is correctly computed
//
Factor = (Parameter - FlatKnots(ii - qq + pp + 1))
/ (FlatKnots(ii + pp) - FlatKnots(ii - qq + pp + 1)) ;
Saved = Factor * BsplineBasis(1,pp) ;
BsplineBasis(1,pp) *= (1.0e0 - Factor) ;
BsplineBasis(1,pp) += BsplineBasis(1,qq) ;
BsplineBasis(1,qq) = Saved ;
}
}
其中：
為節點區間[ui,ui+1)上的基函數左邊部分的系數；
為節點區間[ui-1,ui)上的基函數右邊部分的系數；
六、 程序示例 Sample Codes
將上述內容以一個簡單示例程序來驗證，程序代碼如下所示：
/*
* Copyright (c) 2013 eryar All Rights Reserved.
*
* File : Main.cpp
* Author : eryar@163.com
* Date : 2013-03-09
* Version :
*
* Description : Learn the B-Spline Curve library in the Open Cascade.
*
*/
#include <BSplCLib.hxx>
#include <math_Matrix.hxx>
#include <TColStd_Array1OfReal.hxx>
#include <TColStd_Array1OfInteger.hxx>
#include <Geom2d_BSplineCurve.hxx>
#pragma comment(lib, "TKernel.lib")
#pragma comment(lib, "TKMath.lib")
#pragma comment(lib, "TKG2d.lib")
int main(int argc, char* argv[])
{
// Knot vector: [0,0,0,1,2,3,4,4,5,5,5]
TColStd_Array1OfReal knotSeq(1, 11);
knotSeq.Init(0);
knotSeq.SetValue(1, 0);
knotSeq.SetValue(2, 0);
knotSeq.SetValue(3, 0);
knotSeq.SetValue(4, 1);
knotSeq.SetValue(5, 2);
knotSeq.SetValue(6, 3);
knotSeq.SetValue(7, 4);
knotSeq.SetValue(8, 4);
knotSeq.SetValue(9, 5);
knotSeq.SetValue(10, 5);
knotSeq.SetValue(11, 5);
cout<<"Knot Sequence: [ ";
for (Standard_Integer i = 1; i <= knotSeq.Length(); i++)
{
cout<<knotSeq.Value(i)<<" ";
}
cout<<"]"<<endl;
Standard_Integer knotsLen = BSplCLib::KnotsLength(knotSeq);
TColStd_Array1OfReal knots(1, knotsLen);
TColStd_Array1OfInteger mults(1, knotsLen);
// Test Knots, Mults and Knot sequence of BSplCLib.
BSplCLib::Knots(knotSeq, knots, mults);
cout<<"Knots: [ ";
for (Standard_Integer i = 1; i <= knots.Length(); i++)
{
cout<<knots.Value(i)<<" ";
}
cout<<"]"<<endl;
cout<<"Multiplicity: [ ";
for (Standard_Integer i = 1; i <= mults.Length(); i++)
{
cout<<mults.Value(i)<<" ";
}
cout<<"]"<<endl;
if (BSplCLib::KnotForm(knots, 1, knotsLen) == BSplCLib_Uniform)
{
cout<<"Knots is uniform."<<endl;
}
else
{
cout<<"Knots is non-uniform."<<endl;
}
Standard_Real rValue = 2.5;
Standard_Integer iOrder = 2+1;
Standard_Integer iFirstNonZeroIndex = 0;
math_Matrix bSplineBasis(1, 1, 1, iOrder, 0);
BSplCLib::EvalBsplineBasis(1, 0, iOrder, knotSeq, rValue, iFirstNonZeroIndex, bSplineBasis);
cout<<"First Non-Zero Basis index: "<<iFirstNonZeroIndex<<endl;
cout<<bSplineBasis<<endl;
return 0;
}
上述代碼對節點矢量、重復度的概念的驗證，并以一個實例計算所有非零基函數的值。程序輸出為：
Knot Sequence: [ 0 0 0 1 2 3 4 4 5 5 5 ]
Knots: [ 0 1 2 3 4 5 ]
Multiplicity: [ 3 1 1 1 2 3 ]
Knots is uniform.
First Non-Zero Basis index: 3
math_Matrix of RowNumber = 1 and ColNumber = 3
math_Matrix ( 1, 1 ) = 0.125
math_Matrix ( 1, 2 ) = 0.75
math_Matrix ( 1, 3 ) = 0.125
Press any key to continue . . .
七、 結論 Conclusion
通過學習《The NURBS Book》并給合Open Cascade的源程序，理論聯系實際，使對NURBS的學習更輕松。
根據B樣條基的遞推公式，B樣條曲線的局部性是通過節點來具體實現的。與Bezier曲線不同的就是增加了節點這個參數。根據Cox-deBoor遞推公式親自推導出一次、二次、三次B樣條基函數，可以加深對B樣條曲線的理解。
計算給定節點矢量、次數及參數，計算參數所在區間上所有非零基函數算法的步驟為：
l 通過二分法查找出參數所在的節點區間；
l 根據B樣條基的局部支撐性，計算出所在節點區間上所有非零基函數；
八、 致謝 Acknowledgments
感謝曉天的支持與鼓勵。
九、 參考文獻 Bibliography
1． 趙罡，穆國旺，王拉柱譯Les Piegl，Wayne Tiller The NURBS Book(Second Edition) 2010 清華大學出版社
2． 莫容，常智勇 計算機輔助幾何造型技術 2009 科學出版社
3． 王仁宏，李崇君，朱春鋼 計算幾何教程 2008 科學出版社