曲线数学之贝塞尔曲线Bézier Curves

Tags: ,

本文主要关注的是公式的推导。

在讲贝塞尔曲线之前先复习下组合数学。

组合数学

排列 permutation

注意,排列的英文是permutation,这个词也就是线性代数里的“置换”。联想置换矩阵的概念,就可以近似理解“排列”的意义。

permutation的公式是:

P(n,k)=Pnk=n!(nk)!

含义:从n个数取出某k个数总共有多少种排列。

所以,排列是有顺序的。

组合 combination

组合这个东西可以用“排列”来理解,比如对于某3个不同的数a、b、c,有P33=3!(33)!=6种排列: abc、acb、bac、bca、cab、cba,然而,组合数只有一个,也就是[a、b、c]。

所以,组合是没有顺序的。

combination的公式是:

C(n,k)=Cnk=PnkPkk=n!k!(nk)!,0<k<=n

性质1:

Cnk=n!k!(nk)!

=n!(n(nk))!(nk)!

=n!(nk)!(n(nk))!

=Cnnk

性质2:

Cnk=n!k!(nk)!=0k<=0k>n

(当k<=0或k>n时,右边的式子会出现负数的阶乘。负数阶乘在组合数公式的计算中,可以认为等于0)

combination很重要,比如说二项式定理里的展开式就用到了它:

(a+b)n=r=0nCnranrbr

combination还有一条公式要注意下:

Cnk=Cn1k1+Cn1k,0<k<n

顺便给个简单证明:

Cn1k1+Cn1k=

(n1)!(n1(k1))!(k1)!+(n1)!(n1k)!k!=

(n1)!(nk)!(k1)!+(n1)!(n1k)!k!=

k(n1)!k(nk)!(k1)!+(nk)(n1)!(nk)(n1k)!k!=

k(n1)!(nk)!k!+(nk)(n1)!(nk)!k!=

k(n1)!+(nk)(n1)!(nk)!k!=

n(n1)!(nk)!k!=n!(nk)!k!=Cnk

贝塞尔曲线

定义

给定n个控制点P0,P1,,Pn,贝塞尔曲线的公式如下:

P(t)=i=0nPiBi,n(t),t[0,1]

其中的Bi,n(t)叫Bernstein polynomial(伯思斯坦多项式),定义如下:

Bi,n(t)=Cniti(1t)ni=n!i!(ni)!ti(1t)nii=0,1,,n

Bézier曲线的特性

特性1:改变单个控制点会引起整条曲线的改变。

Bernstein polynomial的特性

递归性

Bi,n(t)=(1t)Bi,n1(t)+tBi1,n1(t)i=0,1,,n

证明:

Bi,n(t)

=Cniti(1t)ni

=(Cn1i1+Cn1i)ti(1t)ni

=Cn1i1ti(1t)ni+Cn1iti(1t)ni

=(1t)Cn1iti(1t)(n1)i+tCn1i1ti1(1t)(n1)(i1)

=(1t)Bi,n1(t)+tBi1,n1(t)

归一性

i=0nBi,n(t)1t(0,1)

证明:

根据二项式定理有:

i=0nBi,n(t) =i=0nCniti(1t)ni [(1t)+t]n1

Partition of Unity

i=0nBi,n(t)=i=0n1Bi,n1(t)=1

证明:

利用递归公式,有:

i=0nBi,n(t)

=i=0n[(1t)Bi,n1(t)+tBi1,n1(t)]

=i=0n[(1t)Bi,n1(t)]+i=0n[tBi1,n1(t)]

=(1t)i=0n[Bi,n1(t)]+ti=0n[Bi1,n1(t)]

=(1t)[i=0n1Bi,n1(t)+Bn,n1(t)]+t[i=1nBi1,n1(t)+B1,n1(t)]

因为:

Bn,n1(t)=Cn1ntn(1t)n1n=0

B1,n1(t)=Cn11t1(1t)(n1)(1)=0

(这里利用了上文提到的组合数公式性质2)

所以可简化为:

i=0nBi,n(t)

=(1t)[i=0n1Bi,n1(t)]+t[i=1nBi1,n1(t)]

=(1t)[i=0n1Bi,n1(t)]+t[i=0n1Bi,n1(t)]

=(1t+t)[i=0n1Bi,n1(t)]

=i=0n1Bi,n1(t)

对称性

Bi,n(1t)=Bni,n(t)

证明:

由定义:

Bi,n(t)=Cniti(1t)ni

有:

Bni,n(t)=Cnnitni(1t)n(ni)

=Cnnitni(1t)i

=Cnitni(1t)i

Bi,n(1t)=Cni(1t)i(1(1t))ni

=Cni(1t)itni

=Cnitni(1t)i

得证。

非负性

当 t = 0 时:

Bi,n(0)=0,i>0

Bi,n(0)=1,i=0

当 t = 1 时:

Bi,n(0)=0,i>0

Bi,n(0)=1,i=0

当 t= (0,1) 时:

Bi,n(t)>0,i=0,1,2,,n1

证明:把数值代入定义公式就可以了。

Bernstein基(Bernstein Basis)到幂基(Power Basis)的转换

由二项式定理:

(a+b)n=r=0nCnranrbr

得到:

(1t)n=r=0nCnr1nr(t)r=r=0nCnr(t)r

所以:

Bi,n(t)=Cniti(1t)ni

=Cnitik=0niCnik(t)k

=Cnitik=0niCnik(1)ktk

=k=0niCnitiCnik(1)ktk

=k=0niCniCnik(1)ktk+i

这时设g = k + i ,则有 k = g - i,i = g - k,上式变成:

=gi=0n(gk)CniCnigi(1)gitg

=g=ing+kCniCnigi(1)gitg

把g换成k,上式变成:

=k=ink+kCniCniki(1)kitk

=k=inCniCniki(1)kitk

其中CniCniki可以进一步简化:

CniCniki =n!i!(ni)!(ni)!(ki)!((ni)(ki))!

=n!i!(ni)!(ni)!(ki)!((nk)!

=n!(ni)!i!(ni)!(ki)!((nk)!

=n!i!(ki)!((nk)!

=n!k!i!(ki)!((nk)!k!

=n!k!(nk)!k!i!((ki)!

=CnkCki

所以:

k=inCniCniki(1)kitk=k=inCnkCki(1)kitk

综上:

Bi,n(t)=k=inCnkCki(1)kitk

bk,i=CnkCki(1)ki

则上式变成:

Bi,n(t)=k=inbk,itk

展开后:

Bi,n(t)=bi,iti+bi+1,iti+1++bn,itn

Bézier曲线的递推形式(de Casteljau算法)

前面讲的是Bézier曲线的曲线方程定义,现在介绍一个简单实用的算法:de Casteljau's Algorithm。

先分享我找到的一些演示程序:

http://myst729.github.io/bezier-curve/

https://www.jasondavies.com/animated-bezier/

以及油管上的:https://www.youtube.com/watch?v=YATikPP2q70

递推公式如下:

Pik={Pik=0(1t)Pik1+tPi+1k1k=1,2,,n,i=0,1,,nk

2.png

以上图为例演示下这条公式:

因为有P0,P1,P2,P34个控制点,所以n的值是3(要减1)。

然后求该贝塞尔曲线在 t = 1/2时的坐标点B(1/2)的步骤如下:

k = 0时,i=0,1,,nk=0,1,2,3,所以有:

Pik=0=Pi

P0k=0=P0

P1k=0=P1

P2k=0=P2

P3k=0=P3

k = 1时,i=0,1,,nk=0,1,2,所以有:

Pik=1=(1t)Pi11+tPi+111=(1t)Pi0+tPi+10

P0k=1=(1t)P00+tP10=0.5P0+0.5P1=m0

P1k=1=(1t)P10+tP20=0.5P1+0.5P2=m1

P2k=1=(1t)P20+tP30=0.5P2+0.5P3=m2

k = 2时,i=0,1,,nk=0,1,所以有:

Pik=2=(1t)Pi1+tPi+11

P0k=2=(1t)P01+tP11=0.5m0+0.5m1=q0

P1k=2=(1t)P11+tP21=0.5m1+0.5m2=q1

k = 3时,i=0,1,,nk=0,所以有:

Pik=3=(1t)Pi2+tPi+12

P0k=3=(1t)P02+tP12=0.5q0+0.5q1=B(1/2)

Bézier曲线的矩阵形式

由上上一节推导出来的2条式子:

bk,i=CnkCki(1)ki

Bi,n(t)=k=inbk,itk=bi,iti+bi+1,iti+1++bn,itn

可以推导出矩阵:

[B0,n(t)B1,n(t)Bn,n(t)]

=[k=0nbk,0tkk=1nbk,1tkk=nnbk,ntk]

=[t0t1t2tn][b0,0000b1,0b1,100b2,0b2,1b2,20bn,0bn,1bn,2bn,n]

再由Bézier曲线的公式:

P(t)=i=0nPiBi,n(t),t[0,1]

有:

P(t)=[t0t1t2tn][b0,0000b1,0b1,100b2,0b2,1b2,20bn,0bn,1bn,2bn,n][P0P1P2Pn]

注意,中间的矩阵B是常量(取决于阶数):

n = 2时:

bk,i=C2kCki(1)ki

b0,0=C20C00(1)00=11=1

b1,0=C21C10(1)10=21(1)=2

b2,0=C22C20(1)20=11=1

b1,1=C21C11(1)11=21=2

b2,1=C22C21(1)21=12(1)=2

b2,2=C22C22(1)22=11=1

B=[100220121]

n = 3时:

B=[1000330036301331]

测试一下正确性

测试代码基于我正在开发中的renderer

#include "transform.hpp"
#include "geometry.hpp"

using namespace renderer;

int main(int argc, char ** argv){
    Matrix3x3 P = {
        20.f, 20.f,     0,  //P0
        770.f, 30.f,    0,  //P1
        400.f, 780.f,   0,  //P2
    };

    Matrix3x3 B = {
        1,  0,  0,
        -2, 2,  0,
        1,  -2, 1,
    };

    Matrix3x3 BP = B * P;

    typedef Matrix<MxN<float, 1, 3>> Matrix1x3;

    cil::CImg<unsigned char> img(800, 800, 1, 3);

    img.atXYZC(P[0], P[1], 0, 1) = 255;
    img.atXYZC(P[3], P[4], 0, 1) = 255;
    img.atXYZC(P[6], P[7], 0, 1) = 255;

    for (float t = 0; t < 1; t += 0.0001f) {
        Matrix1x3 T = { 1, t, t * t };
        Matrix1x3 TBP = T * BP;
        int x = int(TBP[0]);
        int y = int(TBP[1]);
        img.atXYZC(x, y, 0, 0) = 255;
    }
    img.display("");
    return 0;
}

1.png

(未经授权禁止转载)
Written on December 26, 2015

写作不易,您的支持是我写作的动力!