OPEN CASCADE BSpline Curve Interpolation
eryar@163.com
Abstract. Global curve interpolation to point data is a way to construct curves. The paper focus on the interpolate algorithm in OPEN CASCADE, and give a simple example to demonstrate the usage of the GeomAPI_Interpolate class.
Key Words. Interpolate, NURBS, BSpline, OPEN CASCADE
1.Introduction
曲線曲面擬合Curve and Surface Fitting的方式可以分為兩類：插值interpolation和逼近approximation。采用插值的方式時，所創(chuàng)建的曲線或曲面必須精確地滿足所給的數(shù)據(jù)條件，例如曲線通過所給的插值點。采用逼近的方式時，創(chuàng)建的曲線或曲面不必精確地滿足所給的數(shù)據(jù)條件，只要在一定的誤差范圍內(nèi)接近即可，如下圖所示：
Figure 1.1 A curve interpolating five points and end derivatives(The NURBS Book)
Figure 1.2 A curve approximating points(The NURBS Book)
本文先簡要介紹B樣條插值的原理，再結合OPEN CASCADE源碼說明如何對給定點插值B樣條曲線及OPEN CASCADE中插值曲線的一些注意事項。
2.Global Interpolation
假設給定一組數(shù)據(jù)點{Qk}，k=0,1,…n，我們想要用一條p次非有理B樣條曲線插值于這些點。如果我們?yōu)槊總€點Qk指定了一個參數(shù)值uk，并且選定了一個合適的節(jié)點矢量U，我們就可以建立一個系數(shù)矩陣為(n+1)x(n+1)的線性方程組：
n+1個控制點Pi是未知量。剩下的問題是如何確定Qk對應的參數(shù)值uk及節(jié)點矢量U，這將影響到曲線的形狀和參數(shù)化。常見的選取uk的方法有均勻參數(shù)化法、弦長參數(shù)化法和向心參數(shù)化法。
2.1 Equally Spaced 均勻參數(shù)法
假設參數(shù)限定在[0，1]范圍內(nèi)，那么
當參數(shù)范圍為[a,b]時，
均勻參數(shù)化法是最簡單的構造參數(shù)的方法，但是不推薦采用這種方法，因為當數(shù)據(jù)點分布步均勻時，會產(chǎn)生很奇怪的形狀，如打圈自交。
Figure 2.1 B-Spline curve interpolation with the uniformly spaced method[1]
2.2 Chord Length 弦長參數(shù)法
令d為總弦長，且把參數(shù)范圍限定在[0，1]之間，則：
這是最常用的方法，并且一般用它就足夠了，考慮到弦長參數(shù)化接近曲線的均勻參數(shù)化，在這種意義下，它給出了曲線的一個“好”的參數(shù)化。
2.3 Centripetal Method 向心參數(shù)法
這是一個更新的方法，當數(shù)據(jù)點急轉彎變化時，這個方法能得到比弦長參數(shù)化更好的結果。
3.BSpline Interpolation in OPEN CASCADE
OPEN CASCADE對曲線的插值是通過GeomAPI包中的GeomAPI_Interpolate實現(xiàn)的。由其代碼注釋可知，這個類的功能是可以對一系列點進行插值得到C2連續(xù)的B樣條曲線，當對插值點處的切矢不作要求時。對點直接插值的構造函數(shù)為：
GeomAPI_Interpolate (const Handle< TColgp_HArray1OfPnt > &Points, const Standard_Boolean PeriodicFlag, const Standard_Real Tolerance)
其中參數(shù)Points為插值點，PeriodicFlag為是否周期曲線，Tolerance是對插值點進行檢查用的容差。Tolerance容易產(chǎn)生誤解，根據(jù)插值曲線的定義，插值曲線是要求通過插值點的，所以不存在插值得到的曲線和插值點之間的容差。
經(jīng)過查看OPEN CASCADE中插值曲線的源碼，可以得出對曲線進行插值的步驟如下：
v 檢查是否有重復的插值點CheckPoints；
v 生成參數(shù)BuildParameters；
v 使用BSplCLib::Interpolate()進行插值；
v 根據(jù)參數(shù)及次數(shù)生成系數(shù)矩陣，再結合插值點，對系數(shù)矩陣和插值點組成的方程組進行求解。
檢查插值點代碼如下：
static Standard_Boolean CheckPoints(const TColgp_Array1OfPnt& PointArray,
const Standard_Real Tolerance)
{
Standard_Integer ii ;
Standard_Real tolerance_squared = Tolerance * Tolerance,
distance_squared ;
Standard_Boolean result = Standard_True ;
for (ii = PointArray.Lower() ; result && ii < PointArray.Upper() ; ii++) {
distance_squared =
PointArray.Value(ii).SquareDistance(PointArray.Value(ii+1)) ;
result = (distance_squared >= tolerance_squared) ;
}
return result ;
}
由上述代碼可知，Tolerance主要是用于檢查插值點是否在容差范圍內(nèi)有重合現(xiàn)象。生成參數(shù)代碼如下：
static void BuildParameters(const Standard_Boolean PeriodicFlag,
const TColgp_Array1OfPnt& PointsArray,
Handle(TColStd_HArray1OfReal)& ParametersPtr)
{
Standard_Integer ii,
index ;
Standard_Real distance ;
Standard_Integer
num_parameters = PointsArray.Length() ;
if (PeriodicFlag) {
num_parameters += 1 ;
}
ParametersPtr =
new TColStd_HArray1OfReal(1,
num_parameters) ;
ParametersPtr->SetValue(1,0.0e0) ;
index = 2 ;
for (ii = PointsArray.Lower() ; ii < PointsArray.Upper() ; ii++) {
distance =
PointsArray.Value(ii).Distance(PointsArray.Value(ii+1)) ;
ParametersPtr->SetValue(index,
ParametersPtr->Value(ii) + distance) ;
index += 1 ;
}
if (PeriodicFlag) {
distance =
PointsArray.Value(PointsArray.Upper()).
Distance(PointsArray.Value(PointsArray.Lower())) ;
ParametersPtr->SetValue(index,
ParametersPtr->Value(ii) + distance) ;
}
}
由上述代碼可知，OPEN CASCADE插值生成參數(shù)的方法如下：
不是上述三種常用方法的之一，和弦長參數(shù)化法類似，但是沒有去除以總弦長。生成節(jié)點矢量之前為了得到曲線的次數(shù)，做了如下處理：
if (num_poles == 2 && !myTangentRequest) {
degree = 1 ;
}
else if (num_poles == 3 && !myTangentRequest) {
degree = 2 ;
num_distinct_knots = 2 ;
}
else {
degree = 3 ;
num_poles += 2 ;
if (myTangentRequest)
for (ii = myTangentFlags->Lower() + 1 ;
ii < myTangentFlags->Upper() ; ii++) {
if (myTangentFlags->Value(ii)) {
num_poles += 1 ;
}
}
}
由上述代碼可知，插值要求至少有兩個插值點。當只有兩個插值點時，插值曲線次數(shù)為1，即為直線；當有三個插值點且沒有切矢的要求時，插值曲線次數(shù)為2次；當插值點數(shù)多于3個時，插值曲線次數(shù)為3。即對于多于三個點進行插值時，最高只能得到3次曲線，也即C2連續(xù)的曲線。進行B樣條插值的代碼如下：
void BSplCLib::Interpolate(const Standard_Integer Degree,
const TColStd_Array1OfReal& FlatKnots,
const TColStd_Array1OfReal& Parameters,
const TColStd_Array1OfInteger& ContactOrderArray,
const Standard_Integer ArrayDimension,
Standard_Real& Poles,
Standard_Integer& InversionProblem)
{
Standard_Integer ErrorCode,
UpperBandWidth,
LowerBandWidth ;
// Standard_Real *PolesArray = &Poles ;
math_Matrix InterpolationMatrix(1, Parameters.Length(),
1, 2 * Degree + 1) ;
ErrorCode =
BSplCLib::BuildBSpMatrix(Parameters,
ContactOrderArray,
FlatKnots,
Degree,
InterpolationMatrix,
UpperBandWidth,
LowerBandWidth) ;
if(ErrorCode)
Standard_OutOfRange::Raise("BSplCLib::Interpolate");
ErrorCode =
BSplCLib::FactorBandedMatrix(InterpolationMatrix,
UpperBandWidth,
LowerBandWidth,
InversionProblem) ;
if(ErrorCode)
Standard_OutOfRange::Raise("BSplCLib::Interpolate");
ErrorCode =
BSplCLib::SolveBandedSystem(InterpolationMatrix,
UpperBandWidth,
LowerBandWidth,
ArrayDimension,
Poles) ;
if(ErrorCode)
Standard_OutOfRange::Raise("BSplCLib::Interpolate");
}
先是根據(jù)參數(shù)及插值曲線次數(shù)生成系數(shù)矩陣，再對系數(shù)矩陣和插值點構成的方程組進行求解，計算出B樣條曲線的控制頂點Poles。有了節(jié)點矢量，次數(shù)及控制頂點，B樣條就確定下來了：
myCurve =
new Geom_BSplineCurve(poles,
myParameters->Array1(),
mults,
degree) ;
myIsDone = Standard_True ;
OPEN CASCADE提供的插值接口使用還是很簡單的，如對已經(jīng)知點進行插值，其用法如下：
int main(int argc, char* argv[])
{
Handle_TColgp_HArray1OfPnt aPoints = new TColgp_HArray1OfPnt aPoints(1, 3);
Handle_Geom_BSplineCurve aBSplineCurve;
aPoints.SetValue(1, gp_Pnt(0.0, 0.0, 0.0));
aPoints.SetValue(2, gp_Pnt(1.0, 1.0, 0.0));
aPoints.SetValue(3, gp_Pnt(2.0, 6.0, 3.0));
GeomAPI_Interpolate aInterpolater(aPoints, Standard_False, Precision::Approximation());
if (aInterpolater.IsDone())
{
aBSplineCurve = aInterpolater.Curve();
GeomTools::Dump(aBSplineCurve, std::cout);
}
}
4.Conclusion
綜上所述，對給定點進行B樣條插值時，需要考慮參數(shù)值及節(jié)點矢量的選擇。參數(shù)值和節(jié)點矢量確定后，剩下就是利用B樣條基函數(shù)對給定點的參數(shù)計算得到的系數(shù)組成的線性方程進行求解。
在使用OPEN CASCADE的曲線插值類GeomAPI_Interpolate時，需要注間容差Tolerance是用來對插值點進行檢查的，且插值得到的曲線最高只能是三次曲線。
5.Acknowledgments
首先，感謝cnblog提供了一個表現(xiàn)自己的舞臺http://m.shnenglu.com/eryar/，能在網(wǎng)上和世界連通，知道不是一個人在戰(zhàn)斗。
其次，感謝OPEN CASCADE的開源分享，才得以學到幾何造型相關的知識，比起直接啃國內(nèi)教材來，學習效率不可同日而語。正如“Talk is cheap, show me the code”所說，將代碼和書本結合起來學習時，收獲更大。
最后，感謝國內(nèi)外友人對我的肯定和鼓勵，他們自強不息，激情創(chuàng)業(yè)的精神總是讓人興奮。
生活的理想就是為了理想的生活The ideal of life is to live for ideals！人生充滿了起起落落，關鍵在于在頂端時盡情享受，在低谷時不失勇氣。
6.References
1. Hongxin Zhang, Jieqing Feng. B-Spline Interpolation and Approximation. Zhejiang University. 2006-12-18. http://www.cad.zju.edu.cn/home/zhx/GM/009/00-bsia.pdf
2. Les Piegl, Wayne Tiller. The NURBS Book. Springer-Verlag. 1995
3. 趙罡, 穆國旺, 王拉柱譯. 非均勻有理B樣條. 清華大學出版社. 1995
4. 易大義, 陳道琦. 數(shù)值分析引論. 浙江大學出版社. 1998