德卡斯特里奥算法

维基百科,自由的百科全书
(重定向自De Casteljau算法
跳转至: 导航搜索

数学子领域数值分析中的德卡斯特里奥算法英语De Casteljau's algorithm),以发明者保尔·德·卡斯特里奥命名,是计算伯恩斯坦形式的多项式或貝茲曲線递归方法。

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

定义[编辑]

n次波恩斯坦形式给定一个多项式B

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

即多項式

\mathbf{b}_{i,n}(t) = {n\choose i} t^i (1-t)^{n-i},\quad i=0,\ldots n

定義00 = 1。


在点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)=\beta_0^{(n)}.

注意事项[编辑]

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


\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) \qquad \mbox{ , } \in [0,1]

把它变成

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

以及

B_2(t) = \sum_{i=0}^n \beta_{n-i}^{(i)} b_{i,n}\left(\frac{t-t_0}{1-t_0}\right) \qquad \mbox{ , } \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 + \beta_1 2t_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算法分别计算。

伪代码例子[编辑]

这是一个递归的画出一条从点P1P4,弯向P2P3的曲线的伪代码例子。级数参数是递归的次数。该过程用增加了的级数参数来递归的调用它自己。当级别达到最大级别这个全局变量时,在P1P4之间就画上直线。函数中点(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

参看[编辑]