B-Spline Basis Functions
eryar@163.com
摘要Abstract：直接根據B樣條的Cox-deBoor遞推定義寫出計算B樣條基函數的程序，并將計算結果在OpenSceneGraph中顯示。
關鍵字Key Words：B Spline Basis Functions、OpenSceneGraph
一、概述Overview
有很多方法可以用來定義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) 可微性；
二、程序 Codes
直接根據B樣條基函數的Cox-deBoor遞推定義，寫出計算B樣條基函數的程序如下：
頭文件BSplineBasisFunction.h：
/**//*
* Copyright (c) 2013 eryar All Rights Reserved.
*
* File : BSplineBasisFunction.h
* Author : eryar@163.com
* Date : 2013-03-23 22:13
* Version : V1.0
*
* Description : Use Cox-deBoor formula to implemente the
* B-Spline Basis functions.
*
*/
#ifndef _BSPLINEBASISFUNCTION_H_
#define _BSPLINEBASISFUNCTION_H_
#include <vector>
class BSplineBasisFunction
{
public:
BSplineBasisFunction(const std::vector<double>& U);
~BSplineBasisFunction(void);
public:
/**//*
* @brief Binary search of the knot vector.
*/
int FindSpan(double u);
/**//*
* @brief
* @param [in] i: span of the parameter u;
* [in] p: degree;
* [in] u: parameter;
*/
double EvalBasis(int i, int p, double u);
/**//*
* @breif Get knot vector size.
*/
int GetKnotVectorSize(void) const;
/**//*
* @breif Get the knot value of the given index.
*/
double GetKnot(int i) const;
private:
std::vector<double> mKnotVector;
};
#endif // _BSPLINEBASISFUNCTION_H_
實現文件BSplineBasisFunction.cpp：
/**//*
* Copyright (c) 2013 eryar All Rights Reserved.
*
* File : BSplineBasisFunction.cpp
* Author : eryar@163.com
* Date : 2013-03-23 22:14
* Version : V1.0
*
* Description : Use Cox-deBoor formula to implemente the
* B-Spline Basis functions.
*
*/
#include "BSplineBasisFunction.h"
BSplineBasisFunction::BSplineBasisFunction( const std::vector<double>& U )
:mKnotVector(U)
{
}
BSplineBasisFunction::~BSplineBasisFunction(void)
{
}
int BSplineBasisFunction::GetKnotVectorSize( void ) const
{
return static_cast<int> (mKnotVector.size());
}
double BSplineBasisFunction::GetKnot( int i ) const
{
return mKnotVector[i];
}
/**//*
* @brief Binary search of the knot vector.
*/
int BSplineBasisFunction::FindSpan( double u )
{
int iSize = static_cast<int> (mKnotVector.size());
if (u >= mKnotVector[iSize-1])
{
return iSize;
}
int iLow = 0;
int iHigh = iSize;
int iMiddle = (iLow + iHigh) / 2;
while (u < mKnotVector[iMiddle] || u > mKnotVector[iMiddle+1])
{
if (u < mKnotVector[iMiddle])
{
iHigh = iMiddle;
}
else
{
iLow = iMiddle;
}
iMiddle = (iLow + iHigh) / 2;
}
return iMiddle;
}
double BSplineBasisFunction::EvalBasis( int i, int p, double u )
{
if ((i+p+1) >= GetKnotVectorSize())
{
return 0;
}
if (0 == p)
{
if (u >= mKnotVector[i] && u < mKnotVector[i+1])
{
return 1;
}
else
{
return 0;
}
}
double dLeftUpper = u - mKnotVector[i];
double dLeftLower = mKnotVector[i+p] - mKnotVector[i];
double dLeftValue = 0;
double dRightUpper = mKnotVector[i+p+1] - u;
double dRightLower = mKnotVector[i+p+1] - mKnotVector[i+1];
double dRightValue = 0;
if (dLeftUpper != 0 && dLeftLower != 0)
{
dLeftValue = (dLeftUpper / dLeftLower) * EvalBasis(i, p-1, u);
}
if (dRightUpper != 0 && dRightLower != 0)
{
dRightValue = (dRightUpper / dRightLower) * EvalBasis(i+1, p-1, u);
}
return (dLeftValue + dRightValue);
}
主函數：
/**//*
* Copyright (c) 2013 eryar All Rights Reserved.
*
* File : Main.cpp
* Author : eryar@163.com
* Date : 2013-03-23 22:11
* Version : V1.0
*
* Description : Use Cox-deBoor formula to implemente the
* B-Spline Basis functions.
*
*/
#include <osgDB/ReadFile>
#include <osgViewer/Viewer>
#include <osgGA/StateSetManipulator>
#include <osgViewer/ViewerEventHandlers>
#include "BSplineBasisFunction.h"
#pragma comment(lib, "osgd.lib")
#pragma comment(lib, "osgDBd.lib")
#pragma comment(lib, "osgGAd.lib")
#pragma comment(lib, "osgViewerd.lib")
osg::Node* MakeBasisFuncLine(BSplineBasisFunction& bf, int i, int p)
{
// The property basis functions.
int iLen = bf.GetKnotVectorSize();
int iStep = 800;
double dStart = bf.GetKnot(0);
double dEnd = bf.GetKnot(iLen-1);
double dDelta = (dEnd - dStart) / iStep;
double u = 0;
double v = 0;
// Create the Geode (Geometry Node) to contain all our osg::Geometry objects.
osg::Geode* geode = new osg::Geode;
// Create Geometry object to store all the vertices and lines primitive.
osg::ref_ptr<osg::Geometry> linesGeom = new osg::Geometry;
// Set the vertex array to the points geometry object.
osg::ref_ptr<osg::Vec3Array> pointsVec = new osg::Vec3Array;
for (int s = 0; s <= iStep; s++)
{
u = s * dDelta;
v = bf.EvalBasis(i, p, u);
if (v != 0)
{
pointsVec->push_back(osg::Vec3(u, 0, v));
}
}
linesGeom->setVertexArray(pointsVec);
// Set the colors.
osg::ref_ptr<osg::Vec4Array> colors = new osg::Vec4Array;
colors->push_back(osg::Vec4(1.0f, 1.0f, 0.0f, 0.0f));
linesGeom->setColorArray(colors.get());
linesGeom->setColorBinding(osg::Geometry::BIND_OVERALL);
// Set the normal in the same way of color.
osg::ref_ptr<osg::Vec3Array> normals = new osg::Vec3Array;
normals->push_back(osg::Vec3(0.0f, -1.0f, 0.0f));
linesGeom->setNormalArray(normals.get());
linesGeom->setNormalBinding(osg::Geometry::BIND_OVERALL);
// Add the points geometry to the geode.
linesGeom->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::LINE_STRIP, 0, pointsVec->size()));
geode->addDrawable(linesGeom.get());
return geode;
}
osg::Node* CreateScene(void)
{
osg::Group* root = new osg::Group;
// Knot vector: U={0,0,0,1,2,3,4,4,5,5,5}.
std::vector<double> knotVector;
knotVector.push_back(0);
knotVector.push_back(0);
knotVector.push_back(0);
knotVector.push_back(1);
knotVector.push_back(2);
knotVector.push_back(3);
knotVector.push_back(4);
knotVector.push_back(4);
knotVector.push_back(5);
knotVector.push_back(5);
knotVector.push_back(5);
BSplineBasisFunction basisFunc(knotVector);
for (int i = 0; i < basisFunc.GetKnotVectorSize(); i++)
{
//
//root->addChild(MakeBasisFuncLine(basisFunc, i, 1));
//
root->addChild(MakeBasisFuncLine(basisFunc, i, 2));
}
return root;
}
int main(int argc, char* argv[])
{
osgViewer::Viewer viewer;
viewer.setSceneData(CreateScene());
viewer.addEventHandler(new osgGA::StateSetManipulator(viewer.getCamera()->getOrCreateStateSet()));
viewer.addEventHandler(new osgViewer::StatsHandler);
viewer.addEventHandler(new osgViewer::WindowSizeHandler);
return viewer.run();
}
若想顯示出所有次數的B樣條基函數，只需要在CreateScene中添加就好了。
以《The NURBS Book》中的例子2.2，節點矢量U={0, 0, 0, 1, 2, 3, 4, 4, 5, 5, 5}，次數p=2，分別將程序計算的一次、二次B樣條基函數的結果列出，如下圖所示：
圖1. 一次B樣條基函數
圖2. 二次B樣條基函數
本來還想將不同的B樣條基函數以不同的顏色顯示，試了幾次，都沒有成功。若以不同的顏色顯示，會更直觀。若你有設置顏色的方法，歡迎告訴我，eryar@163.com。
三、結論 Conclusion
程序計算結果與書中吻合，效果還不錯。
理解了B樣條的Cox-deBoor遞推定義之后，可以將程序中的遞歸代碼轉換為非遞歸實現，這樣就可以深入理解B樣條基函數了。