德卡斯特里奥算法——找到Bezier曲线上的一个点

http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/Bezier/de-casteljau.html

随着Bezier曲线的构造，接下来最重要的任务就是通过给定的u值查找Bezier曲线上的点C(u)。一个简单的方法就是将u值插入到每一个基函数中，计算每个基函数的乘积及其相应的控制点，并最终将它们加在一起。尽管它有很好的效果，但这并不是数值稳定的（即在求伯恩斯坦多项式的时候可能会引入数值误差）。

在下文中，我们将只写下控制点的编号。即控制点00当作P0，01当作P1，…，0i当作Pi，…，0n当作PN。其中0表示初始值或第0次迭代，接下来是1、2、3……

德卡斯特里奥算法的基本思想是在向量AB上选择一个点C，使得C分向量AB为u:1-u（即∣AC∣:∣AB∣= u）。我们找到一种方法来确定点C。

向量A到B记做B-A，因为u∈[0,1]，点C在u(B-A)上。以A点坐标作为起始点，点C坐标为A + u(B - A) = (1 - u)A + uB。因此给定一个u值，(1 - u)A + uB使C点分向量AB为u:1-u。

德卡斯特里奥算法的思想如下：
假设我们想要找到点C(u)，u∈[0,1]，从第一条折线开始，00-01-02-03…-0n。利用上面的公式，我们可以在线段0i到0(i+1)上找到一点1i，使得该点分向量0i和0(i+1)位u:1-u。以此类推，我们会找到N个点，10, 11, 12, …., 1(n-1)。它们连在一起，构成了(n-1)长度的新的折线。

上图中，u等于0.4。点10在00和01线段中，点11在01和02线段中，……，点14在04和05线段中。所有这些新加入的点都用蓝色标记。

新点的编号是1i，应用这条新折线我们将会获得由n-1个点组成的第二条折线，20, 21, …, 2(n-2)，长度n-2。通过这条折线，我们又可以获得由n-2个点组成的第三条折线，30, 31, …, 3(n-3)，长度n-3。重复N次这个过程，就会最终产生一个唯一的点N0。德卡斯特里奥证明了这个点就是点C(u)在曲线上的点。

我们继续上面的计算，从图中可以看到，点20在10和11线段中分线段为u:1-u。同样地，点21在11和12线段中，点22在12和13线段中，点23在13和14线段中。通过20、21、22、23四个点连接组成了第三条折线。第三条折线包含4个点和3个线段。持续这么做，我们将会得到由30、31、32三个点组成的第四条折线。通过第四条折线，我们将会得到由40、41两个点组成的第五条折线。再做一次，我们获得了点50。这个点就是C(0.4)在曲线上的点。

以上就是德卡斯特里奥算法的几何说明，最优美的曲线设计之一。

实际计算

通过上面对德卡斯特里奥算法的几何说明，我们提出了一个计算方法，如下图所示。

首先，所有给定的控制点被排成一列，在图中最左边所示。每一对相邻的控制点通过两个箭头指向新的控制点。例如，两个相邻控制点ij和i(j+1)生成新的控制点(i+1)j。指向东南方向的箭头表示在它尾部的点乘以(1-u)，指向东北方向的箭头表示在它尾部的点乘以u。ij和i(j+1)即为新点的总和。

因此，从最初的第0列，我们可以计算出第1列。通过第1列，我们又能获得第2列。以此计算下去，在n此迭代后，最终我们会获得一个单独的点n0，这个点就是曲线上的点。下面的算法总结了我们所讨论的问题。由n+1个点构成的数组P和在0到1范围内的一个u值，返回Bezier曲线上的点C(u)。

Input: array P[0:n] of n+1 points and real number u in [0,1]
Output: point on curve, C(u)
Working: point array Q[0:n] for i := 0 to n do
Q[i] := P[i]; // save input
for k := 1 to n do
for i := 0 to n - k do
Q[i] := (1 - u)Q[i] + u Q[i + 1];
return Q[0];

一个递推关系

上面的计算可以通过递推表示。开始时，设P0,j作为Pj在0到n次循环后的值，即P0,j是第0列第j次。第i列第j次计算公式如下：

更确切地说，Pi,j是(1-u)Pi-1,j和uPi-1,j+1的和。最终结果（即曲线上的点）是Pn,0。有了这个想法，你可能会很快想出以下递归程序：

function deCasteljau(i,j)
begin
if i = 0 then
return P0,j
else
return (1-u)* deCasteljau(i-1,j) + u* deCasteljau(i-1,j+1)
Pn-1,0

这段程序看起来很简短，但是它非常低效。原因就在于，我们调用deCasteljau(n,0)计算Pn,0，这时else部分又调用两次，调用deCasteljau(n-1,0)计算Pn-1,0，调用deCasteljau(n-1,1)计算Pn-1,1。

考虑调用Decasteljau(n-1,0)。它分成两次调用，Decasteljau(n-2,0)计算Pn-2,0同时Decasteljau(n-2,1)计算Pn-2,1。调用Decasteljau(n-1,1)分成两次调用，Decasteljau(n-2,1)计算Pn-2,1同时Decasteljau(n-2,2)计算Pn-2,2。因此，Decasteljau(n-2,1)被调用两次。如果我们继续扩大这些函数的调用，我们会发现，为了计算Pi,j的值，几乎所有的函数都被重复调用过，不仅仅被调用一次，可能是很多次。J是重复的，而不是一次，但许多次。这有多糟糕？事实上，上述的计算方案和以下方式计算第n个Fibonacci数相同：

function Fibonacci(n)
begin
if n = 0 or n = 1 then
return 1
else
return Fibonacci (n-1) + Fibonacci (n-2)
end

这段程序消耗了一个指数量级的函数调用来计算Fibonacci数。因此，德卡斯特里奥算法的递归版本不适合直接实现，虽然它看起来简单和优雅！

一个有趣的性质

德卡斯特里奥算法的三角形计算方案给我们提供了一个有趣的性质。观察下面计算，通过8个控制点00, 01, …, 07定义了7阶Bezier曲线。我们把同一列中一组连续的点作为Bezier曲线的控制点。那么，给定一个u值，u∈[0,1]，我们如何在Bezier曲线上计算出这个对应的点？如果德卡斯特里奥算法使用了这些控制点，那么曲线上的点就是由选定的点组成的等边三角形中与这些控制点相对的顶点。

例如，如果选择02，03，04和05四个点，由这四个控制点确定的曲线上，其对应的u点就是32。观察上图蓝色的三角形。如果选中的点是11，12，13，曲线上的点就是31。观察上图黄色的三角形。如果选中的点是30，31，32，33，34，曲线上的点就是70。

基于同样的原因，70也是Bezier曲线控制点60和61定义的点。它也是由52，51和50定义的曲线上的点，并在40，43，42和41定义的曲线上。总的来说，如果我们选择一个点，画一个等边三角形如上图所示，由选中的这些控制点组成的等边三角形的定点就能被计算出来。