德卡斯特里奥算法

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

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

定义

贝兹曲线B（角度为n，控制点${\beta }_0,\dots ,{\beta }_n$）可用以下方式运用德卡斯特里奥算法

$B(t)=\sum _{i=0}^n{\beta }_i{b}_{i,n}(t)$,

其中b为Template:Link-en

$b_{i,n}(t)=(\genfrac{}{}{0ex}{}{n}{i})(1-t{)}^{n-i}{t}^i$.

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

${\beta }_i^{(0)}:={\beta }_i\text{ , }i=0,\dots ,n$

${\beta }_i^{(j)}:={\beta }_i^{(j-1)}(1-t_0)+{\beta }_{i+1}^{(j-1)}{t}_0\text{ , }i=0,\dots ,n-j\text{ , }j=1,\dots ,n$

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

$B(t_0)={\beta }_0^{(n)}.$

再者，贝兹曲线$B$可在$t_0$分成带有各种控制点的两段曲线：

${\beta }_0^{(0)},{\beta }_0^{(1)},\dots ,{\beta }_0^{(n)}$

${\beta }_0^{(n)},{\beta }_1^{(n-1)},\dots ,{\beta }_n^{(0)}$

注意事项

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

$\begin{array}{ccccc}{\beta }_0& ={\beta }_0^{(0)}& & & \\ & & {\beta }_0^{(1)}& & \\ {\beta }_1& ={\beta }_1^{(0)}& & & \\ & & & \ddots & \\ ⋮& & ⋮& & {\beta }_0^{(n)}\\ & & & & \\ {\beta }_{n-1}& ={\beta }_{n-1}^{(0)}& & & \\ & & {\beta }_{n-1}^{(1)}& & \\ {\beta }_n& ={\beta }_n^{(0)}& & & \end{array}$

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

$B(t)=\sum _{i=0}^n{\beta }_i^{(0)}{b}_{i,n}(t)\text{ , }t\in [0,1]$

把它变成

$B_1(t)=\sum _{i=0}^n{\beta }_0^{(i)}{b}_{i,n}\left(\frac{t}{t_0}\right)\text{ , }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)\text{ , }t\in [t_0,1]$

例子

我们要计算2次波恩斯坦多项式，其伯恩斯坦系数为

${\beta }_0^{(0)}={\beta }_0$

${\beta }_1^{(0)}={\beta }_1$

${\beta }_2^{(0)}={\beta }_2$

在t₀点计算。

我们有下式开始递归

${\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{array}{ccc}{\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{array}$

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

贝塞尔曲線

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

$𝐁(t)=\sum _{i=0}^n{𝐏}_i{b}_{i,n}(t)\text{ , }t\in [0,1]$

其中

${𝐏}_i:=\left(\begin{array}{c}{x}_i\\ y_i\\ z_i\end{array}\right)$.

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

$B_1(t)=\sum _{i=0}^n{x}_i{b}_{i,n}(t)\text{ , }t\in [0,1]$

$B_2(t)=\sum _{i=0}^n{y}_i{b}_{i,n}(t)\text{ , }t\in [0,1]$

$B_3(t)=\sum _{i=0}^n{z}_i{b}_{i,n}(t)\text{ , }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<-[0..1]
返回：代表点（1-t）P + tQ的表

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

计算一对控制点间的线性插值的中间结果
用法：eval t b
          t = 线性插值的参数值, t<-[0..1]
          b = 控制点的表
返回：对n个控制点，返回n-1个插值点的表

>	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<-[0..1]
          b = 控制点的表
返回：代表Bezier曲线上一个点的列表

>	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 <nowiki>[[0,0]，[1,1]，[0,1]，[1,0]]</nowiki>

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

Python

（该代码用到Python图像库Template:Wayback）

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算法：计算切比雪夫形式多项式