查看“︁德卡斯特里奥算法”︁的源代码

{{noteTA|G1=Math}}

[[数学]]子领域[[数值分析]]中的'''德卡斯特里奥算法'''（{{lang-en|De Casteljau's algorithm}}），以发明者保尔·德·卡斯特里奥命名，是计算[[伯恩斯坦形式]]的多项式或[[貝茲曲線]]的[[递归]]方法。

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

== 定义 ==

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

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

其中''b''为{{link-en|伯恩施坦多项式|Bernstein polynomial|伯恩施坦基本多项式}}

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

曲线在''t''<sub>0</sub>点上可以用[[遞迴關係式]]运算

:<math>\beta_i^{(0)} := \beta_i \mbox{ , } i=0,\ldots,n</math>
:<math>\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</math>

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

:<math>B(t_0)=\beta_0^{(n)}.</math>

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

== 注意事项 ==

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

:<math>
\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}
</math>

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

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

把它变成

:<math>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]</math>

以及

:<math>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]</math>

== 例子 ==

我们要计算2次波恩斯坦多项式，其伯恩斯坦系数为
:<math>\beta_0^{(0)} = \beta_0</math>
:<math>\beta_1^{(0)} = \beta_1</math>
:<math>\beta_2^{(0)} = \beta_2</math>
在''t''<sub>0</sub>点计算。

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

递归的第二次重复结束于
:<math> 
\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}
</math>

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

== 贝塞尔曲線 ==

在计算带''n''+1个控制点'''P'''<sub>''i''</sub>的三维空间中的''n''次[[贝塞尔曲線]] (Bézier curve) 时

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

其中

:<math>\mathbf{P}_i := 
\begin{pmatrix} 
x_i \\ 
y_i \\
z_i 
\end{pmatrix}
</math>.

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

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

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

== 伪代码例子 ==

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

<syntaxhighlight lang="c">
    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
</syntaxhighlight>

== 代码实现 ==

=== [[Haskell]] ===

<syntaxhighlight lang="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] ]
</syntaxhighlight>

=== [[Python]] ===

（该代码用到[http://www.pythonware.com/products/pil/ ''Python图像库'']{{Wayback|url=http://www.pythonware.com/products/pil/ |date=20120415045020 }}）
<syntaxhighlight lang="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."
</syntaxhighlight>

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

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

[[Category:算法]]
[[Category:数值分析]]
[[Category:样条]]