功能:根据参数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 <double> 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;
}