1 三次Beizer曲线方程介绍

Beizer曲线的一些特性这里不再赘述，大家可以去网上查看一些资料，很详细。最近用到轮廓拟合，所以用三次Beizer曲线效果还可以，有插值和近似拟合（插值就是曲线过点，近似拟合则不过点），就学习了一下。我是做的Beizer曲线插值，插值和近视拟合无非就是控制点选取不一样。
Beizer总方程为
$\sum P_{i}K_{i}^{n}(t)$ (1)，
三次Beizer曲线方程：
$B_{n}(t)=P_{0}(1-t)^{3}+3P_{1}t(1-t)^{2}+3P_{2}t^{2}(1-t)+P_{3}t^{3},t\epsilon [0,1]$ (2)，
这里的$P_1$和$P_2$就是所谓的控制点A、B点。关于控制点AB的求法很多种，我是采用参数设定法构造控制点，思路我是参考百度文库一篇文章，“确定控制点文章链接”。文章写的很详细，对于闭合轮廓的话就想象成一个循环，第0个点的前两个点为n-1和n-2，第n-1个点的后两个点为0和1，文章很好理解。对于参数a和b，通常都是0.25，但我看别人论文中有a=0.125，b=0.05，拟合出来的效果很接近真实曲线，所以我也是采用这个参数。

2 代码实现

上述讲了三次Beizer曲线方程，用在轮廓拟合中怎么实现呢。当然，你先得找到轮廓中的特征点，然后根据两个相邻的特征点拟合成一段三次Beizer曲线，控制点AB是借助周围几个点得到的。思路大概就是这样的。下面讲一下具体步骤。
###2.1 得到控制点AB
控制点AB方程在刚才那篇百度文科文章中有，所以直接根据那个方程来编写代码，下面是我求控制点AB的代码，很简单，大家有需要的可以参考。

//求得控制点AB
void ControlAB(double *Xi,double *Yi, double *Ai_x,double *Ai_y, double *Bi_x,double *Bi_y,int n, double a, double b,int boundType)
{
    if(boundType==1)
    {
        Ai_x[0]=Xi[0]+(Xi[1]-Xi[n-1])*a;
        Ai_y[0]=Yi[0]+(Yi[1]-Yi[n-1])*a;

        Bi_x[n-2]=Xi[n-1]-(Xi[0]-Xi[n-2])*b;
        Bi_y[n-2]=Yi[n-1]-(Yi[0]-Yi[n-2])*b;

        Ai_x[n-1]=Xi[n-1]+(Xi[0]-Xi[n-2])*a;
        Ai_y[n-1]=Yi[n-1]+(Yi[0]-Yi[n-2])*a;
        Bi_x[n-1]=Xi[0]-(Xi[1]-Xi[n-1])*b;
        Bi_y[n-1]=Yi[0]-(Yi[1]-Yi[n-1])*b;
    }
    for(int i=1;i<n-1;i++)
    {
        Ai_x[i]=Xi[i]+(Xi[i+1]-Xi[i-1])*a;
        Ai_y[i]=Yi[i]+(Yi[i+1]-Yi[i-1])*a;
    }
    for(int i=0;i<n-2;i++)
    {
        Bi_x[i]=Xi[i+1]-(Xi[i+2]-Xi[i])*b;
        Bi_y[i]=Yi[i+1]-(Yi[i+2]-Yi[i])*b;
    }
}

代码注意：代码不要瞎贴，这里有些自定义数组，Xi、Yi是点的x、y。

###2.2 拟合曲线生成与绘制
求得AB控制点，可以将$P_1$、$P_2$、A和B代入方程(2)中，得到$P_1$和$P_2$中间的一段三次Beizer曲线方程，就拟合出来了。这时需要看到效果。我用的是C++的MFC画出来的曲线，当然也可以Python，Python一搜有很多例子，网上可以荡到的。我画二维曲线用的是MFC，大家可以参考我的博客“MFC绘制二维曲线”，给个我实验拟合出来的闭合轮廓，其中参数a=0.125，b=0.05。
原轮廓
拟合轮廓

可以看到，拟合效果还是不错的。图中的小黑点是我给的特征点，MFC画出来就是这个样子。

3 总结

三次Beizer曲线拟合算法还是很简单的，主要是控制点的选取，还有最关键的特征点选取。特征点是最基本的。当然Beizer也有一些不足之处，对于非闭合轮廓的情况拟合情况不是很理想，这就需要其他的曲线拟合，我会再发一篇B样条曲线拟合算法的。

关键代码： 
Bezier曲线：
void CalcBZPoints(){
       float a0,a1,a2,a3,b0,b1,b2,b3;
       a0=pt[0].x;
       a1=-3*pt[0].x+3*pt[1].x;
       a2=3*pt[0].x-6*pt[1].x+3*pt[2].x;
       a3=-pt[0].x+3*pt[1].x-3*pt[2].x+pt[3].x;
       b0=pt[0].y;
       b1=-3*pt[0].y+3*pt[1].y;
       b2=3*pt[0].y-6*pt[1].y+3*pt[2].y;
       b3=-pt[0].y+3*pt[1].y-3*pt[2].y+pt[3].y;
       float t = 0;
       float dt = 0.01;
       for(int i = 0; t<1.1; t+=0.1, i++){
              bz[i].x = a0+a1*t+a2*t*t+a3*t*t*t;
              bz[i].y = b0+b1*t+b2*t*t+b3*t*t*t;           
       }
}
画点函数：
void ControlPoint(vector<Point> vpt){
       glPointSize(2);//点的粗细
       for(int i=0; i<vpt.size(); i++){
              glBegin (GL_POINTS);	  
              glColor3f (1.0f, 1.0f, 0.5f);//控制点颜色
			  glVertex2i (vpt[i].x,vpt[i].y);
              glEnd ();
       }
}
连线函数：
void PolylineGL(Point *pt, int num) {
       glBegin (GL_LINE_STRIP);
       for(int i=0;i<num;i++){           
              glColor3f (0.0f, 0.0f, 0.0f);//线的颜色
              glVertex2i (pt[i].x,pt[i].y);           
       }
       glEnd ();
}
画点画线：
void myDisplay(){
       glColor3f (0.0f, 0.0f, 0.0f);//线的颜色为黑色
	   glPointSize(5);//点的宽度
       if (vpt.size() > 0) {
              ControlPoint(vpt);
       }
       if(bDraw){ 
              CalcBZPoints();
              PolylineGL(bz, 11);//画曲线
       }
       glFlush();
}
初始：
void Init(){
	glClearColor(1.0, 1.0, 1.0, 0.0);//背景颜色
	glClear(GL_COLOR_BUFFER_BIT);//清除当前帧的颜色-黑色
	glShadeModel(GL_SMOOTH);
	printf("Please Click left button of mouse to input control point of Bezier Curve!\n");
}
鼠标控制四个控制点：
void mouse(int button, int state, int x, int y){
switch (button){
	case GLUT_LEFT_BUTTON:
	if (state == GLUT_DOWN){
	if (nInput == 0){
	pt[0].x = x;
	pt[0].y = 480 - y;
	printf("(%.3f,%.3f)",pt[0].x/640.0,pt[0].y/480.0);//输出控制顶点坐标0
	nInput = 1;
	vpt.clear();
	vpt.push_back(pt[0]);
	bDraw = false;
	glutPostRedisplay();
	}
	else if (nInput == 1){
	pt[1].x = x;
	pt[1].y = 480 - y;
	printf("(%.3f,%.3f)",pt[1].x/640.0,pt[1].y/480.0);//输出控制顶点坐标1
	vpt.push_back(pt[1]);
	nInput = 2;
	glutPostRedisplay();//
	}
	else if (nInput == 2){
	pt[2].x = x;
	pt[2].y = 480 - y;
	printf("(%.3f,%.3f)",pt[2].x/640.0,pt[2].y/480.0);//输出控制顶点坐标2
	vpt.push_back(pt[2]);
	nInput = 3;
	glutPostRedisplay();//
	}
	else if (nInput == 3){
	pt[3].x = x;
	pt[3].y = 480 - y;
	printf("(%.3f,%.3f)\n",pt[3].x/640.0,pt[3].y/480.0);//输出控制顶点坐标3
	bDraw = true;
	vpt.push_back(pt[3]);
	nInput = 0;
	glutPostRedisplay();
	}
	}
	break;
	default:
	break;
	}
}
void KeyboardFunc(unsigned char Key,int x,int y){
if(Key == 'w' || Key == 'W'){ //按键清屏
        glClear(GL_COLOR_BUFFER_BIT);//清屏
		printf("清屏一次\n");
    }
	if(Key == 'x' || Key == 'X'){ //按键清除上次所画区域——将所画曲线变为背景色
        glColor3f (1.0f, 1.0f, 1.0f);
		if (vpt.size() > 0) {
		ControlPoint(vpt);
		}
		if(bDraw){
			CalcBZPoints();
			glBegin (GL_LINE_STRIP);
			for(int i=0;i<101;i++){
				glColor3f (1.0f, 1.0f, 1.0f);//线颜色
				glVertex2i (bz[i].x,bz[i].y);
			}
			glEnd ();
		}
		glFlush();
		printf("取消上次所画内容\n");
    }
}

由于作业要求，数据会每点击一下，输出一个坐标点。

取消上次所画内容是为了画起来方便点，不过只能取消一次。其实就是给覆盖了。看着网上有说用异或会比较好，但是现在还不会用，以后会试着改进。

int main(int argc, char *argv[])
{
	printf("清屏请按w或W\n取消上次所画内容x或X（仅限于清除线）\n注意不能连续取消代码！\n");
	glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_RGB | GLUT_SINGLE);
	glutInitWindowPosition(100, 100);
	glutInitWindowSize(640, 480);
	glutCreateWindow("Hello World!");
	Init();
    glutKeyboardFunc(KeyboardFunc);//数字、字母键的按键检测的回调函数
	glutDisplayFunc(myDisplay);
	glutReshapeFunc(Reshape);
	glutMouseFunc(mouse);
	glutMainLoop();
	return 0;
}

    其中glutKeyboardFunc(KeyboardFunc);需要写在glutDisplayFunc(myDisplay);前面，不能写最后！

所画如图：这个图是照着个海报画的，水平一般尝试做画

周围会看到一些点。有颜色的话做图比较方便但美观性可能差一些

上一讲讲解了伯恩斯坦多项式，现在就开始对Bezier曲线进行研究。Bezier曲线采用伯恩斯坦多项式作为基函数。

首先，我们定义Bezier曲线的表达式:
$$\begin{array}{r}\mathcal{C}(t)=\sum _{k=0}^n{\mathbf{P}}_k{B}_k^n(t)\end{array}\text{\hspace{1em}(1)}$$

其中，${\left(\mathbf{P}\right)}_{0\le k\le n}$是Bezier的控制顶点。虽然该表达式看起来很简单，但还是不那么直观，下面我们使用代码实现一下Bezier曲线。

计算伯恩斯坦多项式的值

要绘制出Bezier曲线，首先需要先求出$n$次伯恩斯坦多项式在$t$处的值，这里我们使用伯恩斯坦多项式的递归定义
$$B_{k,n}(t)=(1-t)B_{k,n-1}(t)+t{B}_{k-1,n-1}(t)$$
代码实现如下

def all_bernstein(n, t):
    b = np.zeros(n+1)
    b[0] = 1.
    t1 = 1.-t
    for j in range(1, n+1):
        saved = 0.
        for i in range(0, j):
            tmp = b[i]
            b[i] = saved + t1*tmp
            saved = t*tmp
        b[j] = saved
    return b

计算Bezier曲线的值

我们根据式(1)，可以直接求得Bezier曲线在某一点处的值。伯恩斯坦多项式的次数根据控制顶点$\mathbf{P}$的个数确定。

def point_on_bezier_curve(P,t):
    n = len(P) - 1
    b = all_bernstein(n, t)
    c = 0.
    for k in range(0, n+1):
        c += b[k]*P[k]
    return c

一些例子

下面例子中的ts表示在参数域的采样点

def example_1():    
    nt = 200
    ts = np.linspace(0., 1., nt)

    P = np.zeros((2, 2))
    P[:, 0] = [0., 1.]
    P[:, 1] = [1., 0.]

    Q = np.zeros((nt, 2))
    for i,t in enumerate(ts):
        Q[i,:] = point_on_bezier_curve(P,t)

    plt.plot(Q[:,0], Q[:,1], '-b')
    plt.plot(P[:,0], P[:,1], '--ok', linewidth=0.7)
    
    for i in range(0, 2):
        x,y = P[i,:]
        plt.text(x+0.05,y+0.05,'$\mathbf{P}_{' + str(i) + '}$')

    plt.axis([-0.2, 1.2, -0.2, 1.2])

    #ax = plt.axes()
    #ax.set_aspect('equal', 'box')

    
example_1()

def example_3a():    
    nt = 200
    ts = np.linspace(0., 1., nt)

    P = np.zeros((4, 2))
    P[:, 0] = [1., 1., 0., -1.]
    P[:, 1] = [0., 1., 1.,  0.]

    Q = np.zeros((nt, 2))
    for i,t in enumerate(ts):
        Q[i,:] = point_on_bezier_curve(P,t)

    plt.plot(Q[:,0], Q[:,1], '-b')
    plt.plot(P[:,0], P[:,1], '--ok', linewidth=0.7)
    
    for i in range(0, 3):
        x,y = P[i,:]
        plt.text(x+0.05,y+0.05,'$\mathbf{P}_{' + str(i) + '}$')
    i = 3
    x,y = P[i,:]
    plt.text(x-0.05,y+0.05,'$\mathbf{P}_{' + str(i) + '}$')

    plt.axis([-1.2, 1.2, -0.2, 1.2])

    #ax = plt.axes()
    #ax.set_aspect('equal', 'box')

    
example_3a()

复合Bezier曲线

def plot_composite_bezier_curve():
    n = 2

    nt = 200
    ts = np.linspace(0., 1., nt)

    # ...
    P = np.zeros((3, 2))
    P[:, 0] = [1., 1., 0.]
    P[:, 1] = [0., 1., 1.]

    Q = np.zeros((nt, 2))
    for i,t in enumerate(ts):
        Q[i,:] = point_on_bezier_curve(P,t)

    plt.plot(Q[:,0], Q[:,1], '-r')
    plt.plot(P[:,0], P[:,1], '--or', linewidth=0.7)
    # ...

    # ...
    P = np.zeros((3, 2))
    P[:, 0] = [0., -1., 0.]
    P[:, 1] = [1., -1., -1.]

    Q = np.zeros((nt, 2))
    for i,t in enumerate(ts):
        Q[i,:] = point_on_bezier_curve(P,t)

    plt.plot(Q[:,0], Q[:,1], '-b')
    plt.plot(P[:,0], P[:,1], '--ob', linewidth=0.7)
    # ...

    # ...
    P = np.zeros((3, 2))
    P[:, 0] = [1., 1., 0.]
    P[:, 1] = [0., -1., -1.]

    Q = np.zeros((nt, 2))
    for i,t in enumerate(ts):
        Q[i,:] = point_on_bezier_curve(P,t)

    plt.plot(Q[:,0], Q[:,1], '-g')
    plt.plot(P[:,0], P[:,1], '--og', linewidth=0.7)
    # ...

    #plt.axis('equal')
    #plt.title('Composite Bézier curve')

plot_composite_bezier_curve()

一些Beizer曲线的性质

• 伯恩斯坦基函数次数=控制顶点数-1
• 控制顶点数=伯恩斯坦基函数的个数
• 伯恩斯坦基函数的个数=伯恩斯坦基函数次数+1
• 物理域的取值只和控制顶点的取值有关