B-Spline(六):给定参数求点(de Boor 算法)

给定参数u,计算参数曲线上对应点的运算称为求值(Evalute)操作。反之，称之为求解(Solve)。这一节的目标是介绍德布尔(de Boor)在1972年提出的数值稳定的求值方法。它是德利斯特里奥(De Casteljau’s)算法的一般化形式。

这一节我的安排与Dr.Shene略有不同，首先是在章节顺序上，Dr.Shene在B-Spline性质说明后先讲解了B-Spline的求导，然后是在节点插入算法章节中介绍了De Boor算法。我打算先行介绍。另外，Dr.Shene 是从节点插入的角度来讲解的De Boor算法，我因为希望“强调”De Boor与De Casteljau 算法的对应关系，计划参考wikipedia对de Boor算法的讲解来介绍De Boor算法。实质上，Dr.shene对De Boor算法的理解当然更加深入，而我打算给出它的“Naive”实现方法。

给定一个参数\(u \in[u_p,u_{m-p}]\),求解曲线上一点的步骤如下:

1）首先确定 u 所在的区间\([u_k,u_{k+1})\)。因为节点区间为有序非递减数组，因此可以使用std::upper_bound()-1  来确定k。

2）确定在区间\([u_k,u_{k+1})\)的p+1个构成曲线凸包（即对应基函数不为0）的控制点\(P_{k-p},…,P_k\)。

3）对控制点凸包进行p次“切角”。迭代式为 $$\begin{align*} & P_{i,r}=(1-\alpha _{i,r})P_{i-1,r-1}+\alpha _{i,r}P_{i,r-1} \\ & \alpha _{i,r} = \frac{u-u_i}{u_{i+p+1-r}-u_i}  \\ \end{align*}$$需要注意的是对\(\alpha _{i,r}\)的处理，因为有可能存在重复节点，所以规定\(\frac{0}{0}=0\)。

我的C++版本的实现如下：

Point BSpline::Evaluate(double u)
{
	const std::vector<double>& knots = m_vecKnots;
	int p = m_nDegree;

	//查找 u 所在区间
	int k = std::distance(knots.begin(),std::upper_bound(knots.begin(), knots.end(), u))-1;
	
	//申请 p+1个point,分别存放 P_{k-p},...,P_k
	Point *pArray = new Point[p+1];
	std::copy(m_vecCVs.begin() + k - p, m_vecCVs.begin() + k + 1, pArray);

	//迭代 p次
	for (int r = 1;r <= p;++r)
	{
		//i 从 k 到 k-p+1
		for (int i = k,j=p;i >= k - p+r;--i,--j)
		{
			double alpha = u - knots[i];
			double dev = (knots[i+p+1-r]-knots[i]);
			alpha = (dev !=0)? alpha / dev : 0;

			pArray[j] = (1.0 - alpha)*pArray[j - 1] + alpha * pArray[j];
		}
	}

	//所求点是最后一个
	Point pt  = pArray[p];
	delete[]pArray;
	return pt;
}

与德卡斯特利奥算法的对比：

1）德卡斯特利奥算法对控制多边形的“切角”的比例始终是 1-u : u 不变，而德布尔算法的比例与参数区间相关。

2）德卡斯特利奥作用在全部 p+1个控制点上，de Boor算法也是作用在p+1个控制点上，但是不是B样条的全部控制点。这就是我们说的Bezier曲线是全局的，而B样条具有“局部”性。

3）这个是纯从编程角度提醒大家注意的，就是de Casteljau’s 方法给的递推式是 \(P_{i,r} = (1-\alpha)P_{i,r-1}+\alpha P_{i+1,r-1}\)方式的，所以递推时i是递增的，而de-Boor的递推式是\(P_{i,r} = (1-\alpha)P_{i-1,r-1}+\alpha P_{i,r-1}\)方式的，递推时i要递减，否则就会写错。

德布尔因为在对数学研究得卓著贡献，获得2003年美国国家科学奖章（National Medal of Science）