德卡斯特里奧算法

數學子領域數值分析中的德卡斯特里奧算法（英語：De Casteljau's algorithm），以發明者保爾·德·卡斯特里奧命名，是計算伯恩斯坦形式的多項式或貝茲曲線的遞歸方法。

雖然對於大部分的體系結構，該算法和直接方法相比較慢，但它在數值上更為穩定。

定義

貝茲曲線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}$.

曲線在t₀點上可以用遞迴關係式運算

$\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}}$

當選擇一點t₀來計算波恩斯坦多項式時，我們可以用三角形形式的兩個對角線來構造多項式的分段表示。

$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}$

在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{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個控制點P_i的三維空間中的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<-[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圖像庫（頁面存檔備份，存於網際網路檔案館））

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