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样条基函数了。