给定参数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)
你好,请问为什么这里的alpha的表达式和求基函数的时候不一样?这两者不是一个东西吗?
同问??