# De Casteljau算法

## 定义

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$

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

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

## 例子

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

## 貝茲曲線

$\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_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]$

## 伪代码例子

    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


## 示范实现

用线性插值计算P和Q之间的一点R，插值参数为t

P = 代表一个点的表
Q = 代表一个点的表
t = 线性插值的参数值, t<-[0..1]

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

t = 线性插值的参数值, t<-[0..1]
b = 控制点的表

>	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)

t = 线性插值的参数值, t<-[0..1]
b = 控制点的表

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

n = 要计算的点的个数
b = Bezier控制点列表

>	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