數學子領域數值分析中的德卡斯特里奧算法(英語:De Casteljau's algorithm),以發明者保爾·德·卡斯特里奧命名,是計算伯恩斯坦形式的多項式或貝茲曲線的遞歸方法。
雖然對於大部分的體系結構,該算法和直接方法相比較慢,但它在數值上更為穩定。
定義
貝茲曲線B(角度為n,控制點
)可用以下方式運用德卡斯特里奧算法
,
其中b為伯恩施坦基本多項式
.
曲線在t0點上可以用遞歸關係式運算
![{\displaystyle \beta _{i}^{(0)}:=\beta _{i}{\mbox{ , }}i=0,\ldots ,n}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2fc8fa1843a2b22addc72880c7ccf2d68a19b1dc)
![{\displaystyle \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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2d5dd3746c03b87584a6a379595332ff0f30bf0c)
然後,
在
點上的計算可以此算法的
步計算。
的結果為:
![{\displaystyle B(t_{0})=\beta _{0}^{(n)}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1e58a9b6e472c671520531b405d5aedb10e51797)
再者,貝茲曲線
可在
分成帶有各種控制點的兩段曲線:
![{\displaystyle \beta _{0}^{(0)},\beta _{0}^{(1)},\ldots ,\beta _{0}^{(n)}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e80d473388ce635945b4cc751ed704789e8e2966)
![{\displaystyle \beta _{0}^{(n)},\beta _{1}^{(n-1)},\ldots ,\beta _{n}^{(0)}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f09ff951b277c404ca6d4f40ac34187c9cc775d8)
注意事項
進行手算時把系數寫成三角形形式很有用。
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/58fd085f9048ffc54657fc80dfae9bab07d4f783)
當選擇一點t0來計算波恩斯坦多項式時,我們可以用三角形形式的兩個對角線來構造多項式的分段表示。
![{\displaystyle B(t)=\sum _{i=0}^{n}\beta _{i}^{(0)}b_{i,n}(t){\mbox{ , }}\qquad t\in [0,1]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/eb84d69bb346f499c657d7e0737dda3ee72535fa)
把它變成
![{\displaystyle 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}]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ddb872182151fb849a91f142f00b0a7ecb6c0eed)
以及
![{\displaystyle 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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/38ea4239f7425782fc773820c1c01f6450856354)
例子
我們要計算2次波恩斯坦多項式,其伯恩斯坦系數為
![{\displaystyle \beta _{0}^{(0)}=\beta _{0}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b906ff14964ea77ba3a8d1e1484c835a53c54c65)
![{\displaystyle \beta _{1}^{(0)}=\beta _{1}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0c0afa61238e8ecc64244d7f5f80046ae9fb6e52)
![{\displaystyle \beta _{2}^{(0)}=\beta _{2}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ef05314d830a915317409fbbf9e79c43bb2fa87c)
在t0點計算。
我們有下式開始遞歸
![{\displaystyle \beta _{0}^{(1)}=\beta _{0}^{(0)}(1-t_{0})+\beta _{1}^{(0)}t=\beta _{0}(1-t_{0})+\beta _{1}t_{0}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/da05473f0decd8d51d2050a89931402346a7ddc8)
![{\displaystyle \beta _{1}^{(1)}=\beta _{1}^{(0)}(1-t_{0})+\beta _{2}^{(0)}t=\beta _{1}(1-t_{0})+\beta _{2}t_{0}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/120f37276c5210b7ce79149582d873df8ff994f8)
遞歸的第二次重複結束於
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b9bcb28605f63867570a0712edccf9c04a9fa491)
這就是我們所預料的n階伯恩斯坦多項式。
貝塞爾曲線
在計算帶n+1個控制點Pi的三維空間中的n次貝塞爾曲線 (Bézier curve) 時
![{\displaystyle \mathbf {B} (t)=\sum _{i=0}^{n}\mathbf {P} _{i}b_{i,n}(t){\mbox{ , }}t\in [0,1]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7841aeeff9c1b50c687890c1a35194fc036309a8)
其中
.
我們把Bézier曲線分成三個分立的方程
![{\displaystyle B_{1}(t)=\sum _{i=0}^{n}x_{i}b_{i,n}(t){\mbox{ , }}t\in [0,1]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/5fc49c852d0ddc38093382f5562bfd61cbfeb0fb)
![{\displaystyle B_{2}(t)=\sum _{i=0}^{n}y_{i}b_{i,n}(t){\mbox{ , }}t\in [0,1]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b169ad3d55321750016c7bf014e46319f229c1f8)
![{\displaystyle B_{3}(t)=\sum _{i=0}^{n}z_{i}b_{i,n}(t){\mbox{ , }}t\in [0,1]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/464a9fad7d1512dd399ade0d86684e1ba47e9fb2)
然後我們用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
代碼實現
用线性插值计算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圖像庫(頁面存檔備份,存於互聯網檔案館))
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
參看