De Casteljau算法

题目 在Bezier曲线上找到点:De Casteljau算法

翻译文章来自http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/Bezier/de-casteljau.html

在贝塞尔曲线绘制完成后，下一个重要的任务是选定特定u的情况下在曲线上找到点C(u)。一个简单的方法是把u插到每一个基函数上，计算每个其与基函数的乘积以及其相应的控制顶点，最后将其相加。虽然这种方法很好，但是缺乏数值稳定性（例如在计算Bernstein多项式的时候可能引进数值（numerical errors）误差）。

在下文中，我们仅仅记下控制点的数字。例如00对应控制点P0,01对应P1,...,0i对应Pi,...,0n对应Pn。因此在这些数字中0s表示初始的或者是第0次迭代。在之后，0被其他数字1,2,3代替表示第1,2,3次迭代，等等。

De Casteljau算法的基本观点是选择在AB中的一个点C，C将AB分为u:1-u(A到C的距离与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是在A和B之间的点C，将AB分为u:1-u的两段

De Casteljau算法的想法如下。假设我们想要找到C(u)，u在[0,1]中。由第一个多段线00-01-02-03...-0n开始，利用上面的法则找到在线段上的点1i，1i在0i到0(i+1)的连线上并且将这段线分为u:1-u的两部分。依次地，我们可以得到n个点10,11,12,...,1(n-1)，他们定义了一个新的多段线（polyline），一共有n-1段



在上图中，u是0.4,10是在00和01的线段（leg），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，Casteljau已经证明在曲线上的点C(u)对应u

以上图举例，20是在10和11上将这段线分为u:1-u的点，类似地选取21在11和12上，22在12和13上，23在13和14上.第三个多段线是20,21,22,23.第三个多段线有四个点和三条边。继续做得到30,31,32组成的多段线，这是第四个多段线。继续得到有40和41组成的线段，再做一次得到50，为点C(0.4)。

这是de Casteljau算法的几何解释，是在曲线设计中最漂亮（elegant）的结果之一。

实际计算

给定了de Casteljau的几何解释之后，我们现在展示一个计算方法，由下图



首先，在如图中的最左边一列是给定的控制点。对于每一对临近的控制点，可以画出一条右上方和右下方的箭头，并且在两个箭头的交点处写下一个新点。例如相邻的两个点分别为ij 和i(j+1)，新点是(i+1)j，右下方（相对应的左下方）的箭头表示将其尾数ij（相对应的为i(j+1)）乘以1-u（相对应的乘以u），新的点是两个的和。

因此，从初始的第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];

一个递归关系

上面的计算过程可以用递归的方法表示，对于j=0,1,...,n用P0,j表示Pj，也就是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)
end

这个过程看起来简单，表达上也很简练，但是效率却不高。因为我们想要对于Pn,0计算deCasteljau(n,0) ,我们要分别对Pn-1,0计算deCasteljau(n-1,0),0，对Pn-1,1计算deCasteljau(n-1,1)



再考虑计算deCasteljau(n-1,0)。分为对于Pn-2,0计算deCasteljau(n-2,0),和对于Pn-2,1计算deCasteljau(n-2,1)。对于deCasteljau(n-1,1)也是分为两部分如上图所示。因此deCasteljau(n-2,1)被计算了两遍，如果我们扩展这些函数计算，我们会发现对于Pi,j的计算重复了很多遍，而不是仅仅计算了一遍。事实上，上面的计算和计算Fibonacci数列的第n项是类似的

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

为了计算Fibonacci(n)函数被调用了指数次，因此上面的de Casteljau算法对于直接的实现来说是不合适的，虽说代码看上去比较简洁。

一个有趣的观察

de Casteljau算法的三角形计算方法提供了一种有趣的现象，举一个例子，八个控制点00,01,...,07组成的7阶Bezier曲线。选择控制点列中的连续序列，对于给定的u，怎样计算其在Bezier曲线上对应的点呢？事实上，如果利用de Casteljau算法，发现在曲线上的点就是等边三角形的顶点！其底是选择的连续点列。



举例来说，如果我们选择02,03,04,05，其对应于u这四个控制顶点的曲线上的点就是32，就是上图中的蓝色三角形，如果选择11,12,13，在曲线上的点就是31，如黄色三角形。如果选择30,31,32,33和34，那么在曲线上的点就是70.

也是同样的原因，70在由60和61组成的Bezier曲线上。同样地，也是50,51和52上的点，也是40,41,42,43.一般地，如果我们选择一个点并且如上所示画一个等边三角形，其底上的控制点就是被选择用来计算的点。