德卡斯特里奥算法

数学子领域数值分析中的德卡斯特里奥算法（），以发明者保尔·德·卡斯特里奥命名，是计算伯恩斯坦形式的多项式或貝茲曲線的递归方法。

虽然对于大部分的体系结构，该算法和直接方法相比较慢，但它在数值上更为稳定。

定义
贝兹曲线B（角度为n，控制点\beta_0, \ldots, \beta_n）可用以下方式运用德卡斯特里奥算法

:B(t) = \sum_{i=0}^{n}\beta_{i}b_{i,n}(t) ,

其中b为

: b_{i,n}(t) = {n \choose i}(1-t)^{n-i}t^i.

曲线在t0点上可以用遞迴關係式运算

:\beta_i^{(0)} := \beta_i \mbox{ , } i=0,\ldots,n
:\beta_i^{(j)} := \beta_i^{(j-1)} (1-t_0) + \beta_{i+1}^{(j-1)} t_0 \mbox{ , } i = 0,\ldots,n-j \mbox{ , } j= 1,\ldots,n

然后，B在t_0点上的计算可以此算法的n步计算。B(t_0)的结果为：

:B(t_0)=\beta_0^{(n)}.

再者，贝兹曲线B可在t_0分成带有各种控制点的两段曲线：
:\beta_0^{(0)},\beta_0^{(1)},\ldots,\beta_0^{(n)}
:\beta_0^{(n)},\beta_1^{(n-1)},\ldots,\beta_n^{(0)}

注意事项
进行手算时把系数写成三角形形式很有用。

:
\begin{matrix}
\beta_0 & = \beta_0^{(0)} & & & \\
& & \beta_0^{(1)} & & \\
\beta_1 & = \beta_1^{(0)} & & & \\
& & & \ddots & \\
\vdots & & \vdots & & \beta_0^{(n)} \\
& & & & \\
\beta_{n-1} & = \beta_{n-1}^{(0)} & & & \\
& & \beta_{n-1}^{(1)} & & \\
\beta_n & = \beta_n^{(0)} & & & \\
\end{matrix}

当选择一点t0来计算波恩斯坦多项式时，我们可以用三角形形式的两个对角线来构造多项式的分段表示。

:B(t) = \sum_{i=0}^n \beta_i^{(0)} b_{i,n}(t) \mbox{ , } \qquad t \in [0,1]

把它变成

:B_1(t) = \sum_{i=0}^n \beta_0^{(i)} b_{i,n}\left(\frac{t}{t_0}\right) \mbox{ , } \qquad t \in [0,t_0]

以及

:B_2(t) = \sum_{i=0}^n \beta_i^{(n-i)} b_{i,n}\left(\frac{t-t_0}{1-t_0}\right) \mbox{ , } \qquad t \in [t_0,1]

例子
我们要计算2次波恩斯坦多项式，其伯恩斯坦系数为
:\beta_0^{(0)} = \beta_0
:\beta_1^{(0)} = \beta_1
:\beta_2^{(0)} = \beta_2
在t0点计算。

我们有下式开始递归
:\beta_0^{(1)} = \beta_0^{(0)} (1-t_0) + \beta_1^{(0)}t = \beta_0(1-t_0) + \beta_1 t_0
:\beta_1^{(1)} = \beta_1^{(0)} (1-t_0) + \beta_2^{(0)}t = \beta_1(1-t_0) + \beta_2 t_0

递归的第二次重复结束于
:
\begin{matrix}
\beta_0^{(2)} & = & \beta_0^{(1)} (1-t_0) + \beta_1^{(1)} t_0 \\
\ & = & \beta_0(1-t_0) (1-t_0) + \beta_1 t_0 (1-t_0) + \beta_1(1-t_0)t_0 + \beta_2 t_0 t_0 \\
\ & = & \beta_0 (1-t_0)^2 + 2\beta_1 t_0(1-t_0) + \beta_2 t_0^2 \\
\end{matrix}

这就是我们所预料的n阶伯恩斯坦多项式。

贝塞尔曲線
在计算带n+1个控制点Pi的三维空间中的n次贝塞尔曲線 (Bézier curve) 时

:\mathbf{B}(t) = \sum_{i=0}^{n} \mathbf{P}_i b_{i,n}(t) \mbox{ , } t \in [0,1]

其中

:\mathbf{P}_i :=
\begin{pmatrix}
x_i \\
y_i \\
z_i
\end{pmatrix}
.

我们把Bézier曲线分成三个分立的方程

:B_1(t) = \sum_{i=0}^{n} x_i b_{i,n}(t) \mbox{ , } t \in [0,1]
:B_2(t) = \sum_{i=0}^{n} y_i b_{i,n}(t) \mbox{ , } t \in [0,1]
:B_3(t) = \sum_{i=0}^{n} z_i b_{i,n}(t) \mbox{ , } t \in [0,1]

然后我们用de Casteljau算法分别计算。

伪代码例子
这是一个递归的画出一条从点P1到P4，弯向P2和P3的曲线的伪代码例子。级数参数是递归的次数。该过程用增加了的级数参数来递归的调用它自己。当级别达到最大级别这个全局变量时，在P1和P4之间就画上直线。函数中点（midpoint）去两个点，并返回这两点间的线段的中点。

global max_level = 5
procedure draw_curve(P1, P2, P3, P4, level)
if (level > max_level)
draw_line(P1, P4)
else
L1 = P1
L2 = midpoint(P1, P2)
H = midpoint(P2, P3)
R3 = midpoint(P3, P4)
R4 = P4
L3 = midpoint(L2, H)
R2 = midpoint(R3, H)
L4 = midpoint(L3, R2)
R1 = L4
draw_curve(L1, L2, L3, L4, level + 1)
draw_curve(R1, R2, R3, R4, level + 1)
end procedure draw_curve

代码实现
Haskell
用线性插值计算P和Q之间的一点R，插值参数为t
用法：linearInterp P Q t
P = 代表一个点的表
Q = 代表一个点的表
t = 线性插值的参数值, t linearInterp :: [Float]->[Float]->Float->[Float]

linearInterp [] [] _ = [] linearInterp (p:ps) (q:qs) t = (1-t)p + tq : linearInterp ps qs t

计算一对控制点间的线性插值的中间结果
用法：eval t b
t = 线性插值的参数值, t eval :: Float->Float->Float

eval t（bi:bj:[]）= [linearInterp bi bj t] eval t (bi:bj:bs) = (linearInterp bi bj t) : eval t (bj:bs)

用de Casteljau算法计算Bezier曲线上一点
用法：deCas t b
t = 线性插值的参数值, t deCas :: Float->Float->[Float]

deCas t（bi:[]）= bi deCas t bs = deCas t (eval t bs)

用de Casteljau算法计算沿着Bezier曲线的一系列点。点用一个列表返回。
用法：bezierCurve n b
n = 要计算的点的个数
b = Bezier控制点列表
返回：Bezier曲线上n＋1个点
例子：bezierCurve 50 [[0,0]，[1,1]，[0,1]，[1,0]]

bezierCurve :: Int->Float->Float bezierCurve n b = [deCas (fromIntegral x / fromIntegral n) b | x

Python
（该代码用到[http://www.pythonware.com/products/pil/ Python图像库]）

import Image
import ImageDraw

SIZE=128
image = Image.new("RGB", (SIZE, SIZE))
d = ImageDraw.Draw(image)

def midpoint((x1, y1), (x2, y2)):
return ((x1+x2)/2, (y1+y2)/2)

MAX_LEVEL = 5
def draw_curve(P1, P2, P3, P4, level=1):
if level == MAX_LEVEL:
d.line((P1, P4))
else:
L1 = P1
L2 = midpoint(P1, P2)
H = midpoint(P2, P3)
R3 = midpoint(P3, P4)
R4 = P4
L3 = midpoint(L2, H)
R2 = midpoint(R3, H)
L4 = midpoint(L3, R2)
R1 = L4
draw_curve(L1, L2, L3, L4, level+1)
draw_curve(R1, R2, R3, R4, level+1)

draw_curve((10,10),(100,100),(100,10),(100,100))

image.save(r"c:\DeCasteljau.png", "PNG")

print "ok."

参考
Farin, Gerald & Hansford, Dianne (2000). The Essentials of CAGD*. Natic, MA: A K Peters, Ltd. ISBN 1-56881-123-3

参看
*Horner法：计算单项式形式多项式
*Clenshaw算法：计算切比雪夫形式多项式