C++编写任意次clampedB样条曲线（曲线分别与第一个控制点和最后一个控制点的第一边和最后一边相切）

我们知道，B样条曲线根据其起始点与控制点的关系可分为：openB样条曲线，clampedB样条曲线，closedB样条曲线。OpenB样条曲线不会与第一个控制点和最后一个控制点的第一条边和最后一条边相切，closedB样条曲线首尾相连，一般工程运用中，更多的是我们希望B样条曲线的起始点分别于第一个控制点和最后一个控制点的第一条边和最后一条边相切，这种B样条就是clampedB样条。

本文中，生成这种B样条曲线的方法是：假设B样条的次数为p，控制顶点数为n+1，那么节点向量为{u₀，u₁，u₂，，，，u_m}，m=p+n+1.我们只需要将这节点向量的前p+1个节点和最后p+1个节点分别相等，就能生成clampedB样条曲线。即：u₀=u₁=u₂=，，，=u_p，u_m=u_m-1=u_m-2=，，，=u_m-p。代码如下：

1.生成节点向量

std::vector<double> pcxBSpline::get_NoteVector(int p, int N)
{
	std::vector<double> result;
	int n = p + N + 1;
	result.reserve(n);
	double dt = 1.0 / double(n);
	for (int i = 0; i <= n; i++)
	{
		if (i <= p)result.push_back(0.0);
		else if (i < N + 1)
		{
			result.push_back(double(i)*dt);
		}
		else result.push_back(1.0);
	}
	return result;
}

2.生成基函数权值

std::pair<std::vector<double>, int> pcxBSpline::get_BaseWeight(std::vector<double>& notevector, int N, int p, double u)
{
	std::vector<double> result;
	result.resize(N+p+1,0.0);
	std::vector<double>& uu = notevector;
	int k;

	//排除特例;
	if (u == uu[0])
	{
		result[0] = 1.0;
		return std::make_pair(result,0);
	}
	else if (u == uu[uu.size() - 1])
	{
		result[N] = 1.0;
		return std::make_pair(result,N-p);
	}
	else
	{
		for (int i = 0; i < uu.size(); i++)
		{
			if ((uu[i] <= u) && (uu[i + 1]>u))
			{
				k = i;
				break;
			}
		}
	}

	//仅仅计算k-p到k的基函数权值，因为其他控制点的基函数权值为0;
	result[k] = 1.0;
	for (int i = 1; i <= p; i++)
	{
		if (uu[k + 1] - uu[k + 1 - i] == 0.0)result[k - p] = 0.0;
		else result[k - p] = (uu[k + 1] - u) / (uu[k + 1] - uu[k + 1-i])*result[k - p+1];
		for (int j = k - p + 1; j < k; j++)
		{
			double n1, n2;
			if ((uu[j + i] - uu[j]) == 0.0)n1 = 0.0;
			else n1 = (u - uu[j]) / (uu[j + i] - uu[j])*result[j];
			if ((uu[i + j + 1] - uu[j + 1]) == 0.0)n2 =0.0;
			else n2 = (uu[j + i + 1] - u) / (uu[i + j + 1] - uu[j + 1])*result[j + 1];
			result[j] = n1 + n2;
				
		}
		if ((uu[k + i] - uu[k]) == 0.0)result[k] = 0.0;
		else result[k] = (u - uu[k]) / (uu[k + i] - uu[k])*result[k];
	}
	return std::make_pair(result,k-p);
}

3.计算参数化顶点的坐标

CPoint pcxBSpline::get_point(std::vector<CPoint>& contr_points, std::vector<double>& notevector, int N, int p, double u)
{
	std::pair<std::vector<double>,int> baseweight = this->get_BaseWeight(notevector, N, p, u);
	int k = baseweight.second;
	std::vector<double> &W = baseweight.first;
	CPoint result(0.0, 0.0);
	for (int i = 0; i <= p; i++)
	{
		result += (CPoint(contr_points[k+i].x*W[k+i], contr_points[k+i].y*W[k+i]));
	}
	return result;
}

4.运行结果如下图所示