eryar

PipeCAD - Plant Piping Design Software.
RvmTranslator - Translate AVEVA RVM to OBJ, glTF, etc.
posts - 603, comments - 590, trackbacks - 0, articles - 0

B样条基函数——B-Spline Basis Functions

Posted on 2013-03-27 10:27 eryar 阅读(16261) 评论(0)  编辑 收藏 引用 所属分类: 2.OpenCASCADE

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样条基函数,其定义为:

wps_clip_image-29792

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(voidconst;

    
/*
    * @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.0f1.0f0.0f0.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.0f0.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样条基函数的结果列出,如下图所示:

wps_clip_image-3302

图1. 一次B样条基函数

wps_clip_image-6258

图2. 二次B样条基函数

本来还想将不同的B样条基函数以不同的颜色显示,试了几次,都没有成功。若以不同的颜色显示,会更直观。若你有设置颜色的方法,欢迎告诉我,eryar@163.com。

三、结论 Conclusion

程序计算结果与书中吻合,效果还不错。

理解了B样条的Cox-deBoor递推定义之后,可以将程序中的递归代码转换为非递归实现,这样就可以深入理解B样条基函数了。

 

PDF Version and Codes: B-Spline Basis Functions


只有注册用户登录后才能发表评论。
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理