求B样条曲线各个点坐标及C++实现

功能：根据参数u值和k（大小为阶数值）与节点矢量，计算第i个k次B样条基数

输入参数： u—参数值；k—大小值为阶数；i—第i个k次B样条的支撑区间左端节点的下标；aNode为节点向量。

输出参数：返回函数值。

double GetBaseFunVal(double u, int i, int k, vector <double> m_aNode)

{

	  double Val = 0.0;

	  double Val1 = 0.0;

	   double Val2 = 0.0;

 if (k==0)
{   
    if (u < m_aNode[i] || u > m_aNode[i+1])  
  	  return Val;    
    else   
    {
	    Val = 1.0;
	    return Val;
    }
}

if (k>0)
{ 
	  if (u < m_aNode[i] || u > m_aNode[i+k+1])
    { 
	    return Val;
    }
else
{
    double alpha = 0.0;
    double beta = 0.0;
    double dTemp = 0.0;   
    dTemp = m_aNode[i+k] - m_aNode[i];
    if (dTemp == 0.0)
    {
   		 alpha  = 0;
	}
else
	 alpha = (u - m_aNode[i])/dTemp;

dTemp = m_aNode[i+k+1] - m_aNode[i+1];

if (dTemp == 0.0)
   beta = 0.0;
else
	 beta = (m_aNode[i+k+1] - u)/dTemp;

Val1 = alpha * GetBaseFunVal(u, i, k-1, m_aNode);
Val2 = beta * GetBaseFunVal(u, i+1, k-1, m_aNode);
Val = Val1 + Val2;
}

}
return Val;
}

上述功能模块摘自于计算机辅助几何设计与非均匀有理B样条。已知B样条的n+1控制点坐标，以及相应的节点向量，可求得对应的曲线方程。

先计算各个控制点的基函数

各个基函数的求解可根据上述的功能模块求出。

下面是我的C++实现：曲线是二维的，三维的情况，就Z坐标做同X,Y求解方式相同即可。在求解的过程中，我自己在CAD上画了个样条曲线，然后通过GetBaseFunVal(double u, int i, int k, vector m_aNode)和顶点坐标，及节点向量求各个点的坐标。随着u值的变化，计算各个X,Y,Z值。一个星期的摸爬滚打中，能输出图形，但是与原来的图形对应不上。最终找到的原因在与基函数出问题了。在书本等相关资源中，基函数成员中的k表示的是次数，在我画的样条曲线中，阶数显示为3（为什么是3？CAD的标注里，实体块中的 70 下一行，为3），所以我理所当然的写为了2,。一直有问题。我将它改为3以后，竟然奇迹般的可以用了。而且跟原来图形吻合。这个是我的相关经历，希望对你们能有用。另外，哪位热心人士可以说明下，为什么k改为阶数大小，就可以呢？

#include <iostream>
#include <fstream>
#include <afxtempl.h>

using namespace std;

struct tPoint
{
    double x;
    double y;
    double z;
};


double GetBaseFunVal(double u, int i, int k, vector <double> m_aNode)
{
    double Val = 0.0;
    double Val1 = 0.0;
    double Val2 = 0.0;
    if (k==0)
    {
        if (u < m_aNode[i] || u > m_aNode[i+1])
            return Val;
        else
        {
            Val = 1.0;
            return Val;
        }
    }
    if (k>0)
    {
        if (u < m_aNode[i] || u > m_aNode[i+k+1])
        {
            return Val;
        }
        else
        {
            double alpha = 0.0;
            double beta = 0.0;
            double dTemp = 0.0;
            dTemp = m_aNode[i+k] - m_aNode[i];
            if (dTemp == 0.0)
            {
                alpha  = 0;
            }
            else
                alpha = (u - m_aNode[i])/dTemp;
            dTemp = m_aNode[i+k+1] - m_aNode[i+1];

            if (dTemp == 0.0)
            {
                beta = 0.0;
            }
            else
                beta = (m_aNode[i+k+1] - u)/dTemp;
            Val1 = alpha * GetBaseFunVal(u, i, k-1, m_aNode);
            Val2 = beta * GetBaseFunVal(u, i+1, k-1, m_aNode);
            Val = Val1 + Val2;
        }
    }

    return Val;
}

int main()
{
    tPoint tData;
    vector <tPoint> vtData;

    vtData.clear();
    vector <double> nodeVector;
    nodeVector.push_back(0);
    nodeVector.push_back(0);
    nodeVector.push_back(0);
    nodeVector.push_back(0);
    nodeVector.push_back(1);
    nodeVector.push_back(2);
    nodeVector.push_back(3);
    nodeVector.push_back(4);
    nodeVector.push_back(5);
    nodeVector.push_back(6);
    nodeVector.push_back(6);
    nodeVector.push_back(6);
    nodeVector.push_back(6);

    //节点向量nodeVector， 控制点坐标(0,3),(200,100), (750, 200), k=2

    for (double u = 0; u < 6; u=u+0.01)
    {
        // 样条的数据
        tData.x = (GetBaseFunVal(u, 0, 3, nodeVector)*(-7585) + GetBaseFunVal(u, 1, 3, nodeVector)*(-3427.5) + GetBaseFunVal(u, 2, 3, nodeVector)*46087.5
            + GetBaseFunVal(u, 3, 3, nodeVector)*9220.0 +  GetBaseFunVal(u, 4, 3, nodeVector)*(-14835.0) +  GetBaseFunVal(u, 5, 3, nodeVector)*(-2002.5) + GetBaseFunVal(u, 6, 3, nodeVector)*71975
            +  GetBaseFunVal(u, 7, 3, nodeVector)*45235 +  GetBaseFunVal(u, 8, 3, nodeVector)*83150)/*/(GetBaseFunVal(u, 0, 3, nodeVector) + GetBaseFunVal(u, 1, 3, nodeVector) + GetBaseFunVal(u, 2, 3, nodeVector)
            + GetBaseFunVal(u, 3, 3, nodeVector) +  GetBaseFunVal(u, 4, 3, nodeVector) +  GetBaseFunVal(u, 5, 3, nodeVector))*/;
         tData.y = (GetBaseFunVal(u, 0, 3, nodeVector)*(-3807.5) + GetBaseFunVal(u, 1, 3, nodeVector)*(19850.0) + GetBaseFunVal(u, 2, 3, nodeVector)*14335
             +  GetBaseFunVal(u, 3, 3, nodeVector)*(-17582.5) +  GetBaseFunVal(u, 4, 3, nodeVector)*(-5445.0) +  GetBaseFunVal(u, 5, 3, nodeVector)*(-80735.0) + GetBaseFunVal(u, 6, 3, nodeVector)*(-23817.5)
            +   GetBaseFunVal(u, 7, 3, nodeVector)*5037.5 +  GetBaseFunVal(u, 8, 3, nodeVector)*(-9360))/*/(GetBaseFunVal(u, 0, 3, nodeVector) + GetBaseFunVal(u, 1, 3, nodeVector) + GetBaseFunVal(u, 2, 3, nodeVector)
            + GetBaseFunVal(u, 3, 3, nodeVector) +  GetBaseFunVal(u, 4, 3, nodeVector) +  GetBaseFunVal(u, 5, 3, nodeVector))*/;
        tData.z = 0.0;

        vtData.push_back(tData);
    }

    char *file = "C:/Users/Monkey/Desktop/新建文件夹 (2)/TEST/Last.txt";
    ofstream out(file);
    if (!out)
    {
        cout << "打开文件失败！！！！" << endl;
    }  

    for (int n = 0; n < vtData.size(); n++)
    {
        out << vtData[n].x << "  " << vtData[n].y <<"   " << vtData[n].z << endl;
    }

    out.close();

    return 0;
}