贝塞尔曲线应用（贝塞尔插值）

贝塞尔曲线应用（贝塞尔插值）
所谓贝塞尔插值就是指给定n个顶点，要求把这n个顶点连接成为平滑的曲线。那肯定得在这些顶点之间插值了，但这些插值的点怎么找到，可不能随便插值，否则整体上未必是平滑曲线，所以必须找到一个曲线方程，根据这个曲线方程来找到这些插值的点，而且要求这条曲线方程过原来条件中规定的n个顶点。由于贝塞尔曲线可以由几个控制点决定，所以想到用一条贝塞尔曲线作为所求的曲线方程，这就是所谓的贝塞尔插值【个人理解哈】。

下面是一个典型的例子，他根据每两个顶点作为一条贝塞尔曲线的端点（即起始点和终止点），并由这两个端点结合相邻的其它两个顶点求得和这两个端点对应的贝塞尔曲线的控制点，然后根据端点和控制点就绘制一条过两个顶点的贝塞尔曲线了。由于平均每两个顶点决定一条贝塞尔曲线，如果画的多条贝塞尔曲线在相邻交接处也是平滑的，那么画出来的整个曲线就是平滑的。

原文地址：http://www.antigrain.com/research/
bezier_interpolation/index.html#PAGE_BEZIER_INTERPOLATION

Interpolation with Bezier Curves 意思是：贝塞尔插值

A very simple method of smoothing polygons 意思是：一种非常简单的多边形平滑方法

上面的文章是英文的，英文不好的请看下面的翻译

翻译：唐风翻译转自：http://liyiwen.iteye.com/blog/705489

文章内容：

之前 comp.graphic.algorithms 上有一个讨论，是关于怎么样使用曲线对多边形进行插值处理，使得最终产生的曲线是光滑的而且能通过所有的顶点。Gernot Hoffmann 建议说使用著名的 B-Spline 来进行插值。这里有他当时的文章。B-Spline 在这里效果很好，它看起来就像是一个固定在多边形顶点上的橡皮尺（elastic ruler）。

但我有个大胆的推测，我觉得肯定还存在更简单的方法。比如，使用三次贝塞曲线(cubic Bezier)进行近似。贝塞尔曲线有两个固定点（起点和终点），另加两个决定曲线形状的控制点(CP）。关于贝塞尔曲线的更多知识可以在搜索引擎中找到，比如，你可以参考 Paul Bourke 的站点。现在给贝塞尔曲线的锚点（固定点），也就是多边形的某一对顶点，那么问题是，我们怎么计算控制点的位置？我运行 Xara X 然后画出了右边这个图形，这很简单，所以我决定尝试下计算出它们的坐标。很显然，多边形两条相邻边的两个控制点与这两个控制点之间的顶点应该在一条直线上，只有这样，两条相邻的插值曲线才能平滑地连接在一起。所以，这两个控制点应该是相对于顶点是对称的，不过，也不全是……，因为真正的对称就要求它们与中心点的距离应该是相等的，但对于我们的情况中并不完全是这样的。一开始，我试着先算出多边形两条边的角平分线，然后把控制点放在这条角平分线的垂直线上。但从图上可以看到，控制点的连线并不会总是垂直于角平分线的。

最终，我找到一个非常简单的办法，不需要任何复杂的数学计算。首先，我们计算出多边形所有边线的中点，Ai。

然后连接起相邻边中点，得到很多线段，记为 Ci 。并用图记的方法（下图中Step2下面的文字）计算出 Bi 点。

最后一步，只需要简单地将 Ci 进行平移，平移的路径就是每条线段上 Bi 到对应顶点的路径。就这样，我们计算出了贝塞尔曲线的控制点，平滑的结果看起来也很棒。

这里还可以做一点小小的改进，因为我们已经得到了一条决定控制点的直线，所以，我们可以根据需要，使控制点在这条直线上移动，这样可以改变插值曲线的状态。我使用了一个与控制点和顶点初始距离相关的系数 K ，用来沿直线移动控制点。控制点离顶点越远，图形看起来就越锐利。

下面是用原始形状和系统K=1.0的贝塞尔插值两种方法来描画的 SVG 的狮子。

下面是放大图

这个方法对于自相关的多边形也适用，下面的例子可以看到，结果非常有意思：

这个方法只是探索和经验式的，如果从严格的数学模型的角度看它可能是错误的。但在实际使用中的效果已经足够好了，而且这个方法只需要最小的计算量。下面的代码就是用来画出上面狮子图像的。这些代码并没有进行优化，只是用来演示的。里面有些变量计算了两次，在实际程序中，如果连续的步骤中都用到同一个变量值，我们可以先缓存变量值进行复用（以避免重复的计算）。
```
// Assume we need to calculate the control
// points between (x1,y1) and (x2,y2).
// Then x0,y0 - the previous vertex,
//      x3,y3 - the next one. 
double xc1 = (x0 + x1) / 2.0;
double yc1 = (y0 + y1) / 2.0;
double xc2 = (x1 + x2) / 2.0;
double yc2 = (y1 + y2) / 2.0;
double xc3 = (x2 + x3) / 2.0;
double yc3 = (y2 + y3) / 2.0; 
double len1 = sqrt((x1-x0) * (x1-x0) + (y1-y0) * (y1-y0));
double len2 = sqrt((x2-x1) * (x2-x1) + (y2-y1) * (y2-y1));
double len3 = sqrt((x3-x2) * (x3-x2) + (y3-y2) * (y3-y2)); 
double k1 = len1 / (len1 + len2);
double k2 = len2 / (len2 + len3); 
double xm1 = xc1 + (xc2 - xc1) * k1;
double ym1 = yc1 + (yc2 - yc1) * k1; 
double xm2 = xc2 + (xc3 - xc2) * k2;
double ym2 = yc2 + (yc3 - yc2) * k2; 
// Resulting control points. Here smooth_value is mentioned
// above coefficient K whose value should be in range [0...1].
ctrl1_x = xm1 + (xc2 - xm1) * smooth_value + x1 - xm1;
ctrl1_y = ym1 + (yc2 - ym1) * smooth_value + y1 - ym1; 
ctrl2_x = xm2 + (xc2 - xm2) * smooth_value + x2 - xm2;
ctrl2_y = ym2 + (yc2 - ym2) * smooth_value + y2 - ym2; 
```
使用三次贝塞尔近似的代码：
```
// Number of intermediate points between two source ones,
// Actually, this value should be calculated in some way,
// Obviously, depending on the real length of the curve.
// But I don't know any elegant and fast solution for this
// problem.
#define NUM_STEPS 20 
void curve4(Polygon* p,
            double x1, double y1,   //Anchor1
            double x2, double y2,   //Control1
            double x3, double y3,   //Control2
            double x4, double y4)   //Anchor2
{
    double dx1 = x2 - x1;
    double dy1 = y2 - y1;
    double dx2 = x3 - x2;
    double dy2 = y3 - y2;
    double dx3 = x4 - x3;
    double dy3 = y4 - y3; 
    double subdiv_step  = 1.0 / (NUM_STEPS + 1);
    double subdiv_step2 = subdiv_step*subdiv_step;
    double subdiv_step3 = subdiv_step*subdiv_step*subdiv_step; 
    double pre1 = 3.0 * subdiv_step;
    double pre2 = 3.0 * subdiv_step2;
    double pre4 = 6.0 * subdiv_step2;
    double pre5 = 6.0 * subdiv_step3; 

    double tmp1x = x1 - x2 * 2.0 + x3;
    double tmp1y = y1 - y2 * 2.0 + y3; 
    double tmp2x = (x2 - x3)*3.0 - x1 + x4;
    double tmp2y = (y2 - y3)*3.0 - y1 + y4; 
    double fx = x1;
    double fy = y1; 
    double dfx = (x2 - x1)*pre1 + tmp1x*pre2 + tmp2x*subdiv_step3;
    double dfy = (y2 - y1)*pre1 + tmp1y*pre2 + tmp2y*subdiv_step3; 
    double ddfx = tmp1x*pre4 + tmp2x*pre5;
    double ddfy = tmp1y*pre4 + tmp2y*pre5; 
    double dddfx = tmp2x*pre5;
    double dddfy = tmp2y*pre5; 
    int step = NUM_STEPS; 
    // Suppose, we have some abstract object Polygon which
    // has method AddVertex(x, y), similar to LineTo in
    // many graphical APIs.
    // Note, that the loop has only operation add!
    while(step--)
    {
        fx   += dfx;
        fy   += dfy;
        dfx  += ddfx;
        dfy  += ddfy;
        ddfx += dddfx;
        ddfy += dddfy;
        p->AddVertex(fx, fy);
    }
    p->AddVertex(x4, y4); // Last step must go exactly to x4, y4
}
```
你可以下载一个能运行的画狮子的例子，对它进行旋转和缩放，也可以生成一些随机的多边形。点左键并拖动它可以围绕中心点旋转和缩放图像。点右键并从左向右拖动，可以改变系统数K。 K＝1时大约是距窗口左边100像素处。每次双击会产生一个随机的多边形，对于这些多边形，也可以进行旋转、缩放以及改变K值的操作。

运用上述原理，某位网友给出了C++封装好的代码，如下：

为了把一串点连成光滑的曲线，运用贝塞尔曲线的光滑性来穿过这些点。

大致思路就是先算出相邻原始点的中点，在把相邻中点连成的线段平移到对应的原始点，以平移后的中点作为控制点，相邻原始点为起始点画贝塞尔曲线，这样就保证了连接处的光滑。而贝塞尔曲线本身是光滑的，所以就把这些原始点用光滑曲线连起来了。

我封装了一个函数，留着以后用。

（c++版，其它语言只要把数组和可变数组稍微变一下就能用）
[cpp] view plain copy print
1. void createCurve(CvPoint *originPoint,int originCount,vector<CvPoint> &curvePoint){
2. //控制点收缩系数，经调试0.6较好，CvPoint是opencv的，可自行定义结构体(x,y)
3. float scale = 0.6;
4. CvPoint midpoints[originCount];
5. //生成中点
6. for(int i = 0 ;i < originCount ; i++){
7. int nexti = (i + 1) % originCount;
8. midpoints[i].x = (originPoint[i].x + originPoint[nexti].x)/2.0;
9. midpoints[i].y = (originPoint[i].y + originPoint[nexti].y)/2.0;
10. }
12. //平移中点
13. CvPoint extrapoints[2 * originCount];
14. for(int i = 0 ;i < originCount ; i++){
15. int nexti = (i + 1) % originCount;
16. int backi = (i + originCount - 1) % originCount;
17. CvPoint midinmid;
18. midinmid.x = (midpoints[i].x + midpoints[backi].x)/2.0;
19. midinmid.y = (midpoints[i].y + midpoints[backi].y)/2.0;
20. int offsetx = originPoint[i].x - midinmid.x;
21. int offsety = originPoint[i].y - midinmid.y;
22. int extraindex = 2 * i;
23. extrapoints[extraindex].x = midpoints[backi].x + offsetx;
24. extrapoints[extraindex].y = midpoints[backi].y + offsety;
25. //朝 originPoint[i]方向收缩
26. int addx = (extrapoints[extraindex].x - originPoint[i].x) * scale;
27. int addy = (extrapoints[extraindex].y - originPoint[i].y) * scale;
28. extrapoints[extraindex].x = originPoint[i].x + addx;
29. extrapoints[extraindex].y = originPoint[i].y + addy;
31. int extranexti = (extraindex + 1)%(2 * originCount);
32. extrapoints[extranexti].x = midpoints[i].x + offsetx;
33. extrapoints[extranexti].y = midpoints[i].y + offsety;
34. //朝 originPoint[i]方向收缩
35. addx = (extrapoints[extranexti].x - originPoint[i].x) * scale;
36. addy = (extrapoints[extranexti].y - originPoint[i].y) * scale;
37. extrapoints[extranexti].x = originPoint[i].x + addx;
38. extrapoints[extranexti].y = originPoint[i].y + addy;
40. }
42. CvPoint controlPoint[4];
43. //生成4控制点，产生贝塞尔曲线
44. for(int i = 0 ;i < originCount ; i++){
45. controlPoint[0] = originPoint[i];
46. int extraindex = 2 * i;
47. controlPoint[1] = extrapoints[extraindex + 1];
48. int extranexti = (extraindex + 2) % (2 * originCount);
49. controlPoint[2] = extrapoints[extranexti];
50. int nexti = (i + 1) % originCount;
51. controlPoint[3] = originPoint[nexti];
52. float u = 1;
53. while(u >= 0){
54. int px = bezier3funcX(u,controlPoint);
55. int py = bezier3funcY(u,controlPoint);
56. //u的步长决定曲线的疏密
57. u -= 0.005;
58. CvPoint tempP = cvPoint(px,py);
59. //存入曲线点
60. curvePoint.push_back(tempP);
61. }
62. }
63. }
64. //三次贝塞尔曲线
65. float bezier3funcX(float uu,CvPoint *controlP){
66. float part0 = controlP[0].x * uu * uu * uu;
67. float part1 = 3 * controlP[1].x * uu * uu * (1 - uu);
68. float part2 = 3 * controlP[2].x * uu * (1 - uu) * (1 - uu);
69. float part3 = controlP[3].x * (1 - uu) * (1 - uu) * (1 - uu);
70. return part0 + part1 + part2 + part3;
71. }
72. float bezier3funcY(float uu,CvPoint *controlP){
73. float part0 = controlP[0].y * uu * uu * uu;
74. float part1 = 3 * controlP[1].y * uu * uu * (1 - uu);
75. float part2 = 3 * controlP[2].y * uu * (1 - uu) * (1 - uu);
76. float part3 = controlP[3].y * (1 - uu) * (1 - uu) * (1 - uu);
77. return part0 + part1 + part2 + part3;
78. }
void createCurve(CvPoint *originPoint,int originCount,vector<CvPoint> &curvePoint){ //控制点收缩系数，经调试0.6较好，CvPoint是opencv的，可自行定义结构体(x,y) float scale = 0.6; CvPoint midpoints[originCount]; //生成中点 for(int i = 0 ;i < originCount ; i++){ int nexti = (i + 1) % originCount; midpoints[i].x = (originPoint[i].x + originPoint[nexti].x)/2.0; midpoints[i].y = (originPoint[i].y + originPoint[nexti].y)/2.0; } //平移中点 CvPoint extrapoints[2 * originCount]; for(int i = 0 ;i < originCount ; i++){ int nexti = (i + 1) % originCount; int backi = (i + originCount - 1) % originCount; CvPoint midinmid; midinmid.x = (midpoints[i].x + midpoints[backi].x)/2.0; midinmid.y = (midpoints[i].y + midpoints[backi].y)/2.0; int offsetx = originPoint[i].x - midinmid.x; int offsety = originPoint[i].y - midinmid.y; int extraindex = 2 * i; extrapoints[extraindex].x = midpoints[backi].x + offsetx; extrapoints[extraindex].y = midpoints[backi].y + offsety; //朝 originPoint[i]方向收缩 int addx = (extrapoints[extraindex].x - originPoint[i].x) * scale; int addy = (extrapoints[extraindex].y - originPoint[i].y) * scale; extrapoints[extraindex].x = originPoint[i].x + addx; extrapoints[extraindex].y = originPoint[i].y + addy; int extranexti = (extraindex + 1)%(2 * originCount); extrapoints[extranexti].x = midpoints[i].x + offsetx; extrapoints[extranexti].y = midpoints[i].y + offsety; //朝 originPoint[i]方向收缩 addx = (extrapoints[extranexti].x - originPoint[i].x) * scale; addy = (extrapoints[extranexti].y - originPoint[i].y) * scale; extrapoints[extranexti].x = originPoint[i].x + addx; extrapoints[extranexti].y = originPoint[i].y + addy; } CvPoint controlPoint[4]; //生成4控制点，产生贝塞尔曲线 for(int i = 0 ;i < originCount ; i++){ controlPoint[0] = originPoint[i]; int extraindex = 2 * i; controlPoint[1] = extrapoints[extraindex + 1]; int extranexti = (extraindex + 2) % (2 * originCount); controlPoint[2] = extrapoints[extranexti]; int nexti = (i + 1) % originCount; controlPoint[3] = originPoint[nexti]; float u = 1; while(u >= 0){ int px = bezier3funcX(u,controlPoint); int py = bezier3funcY(u,controlPoint); //u的步长决定曲线的疏密 u -= 0.005; CvPoint tempP = cvPoint(px,py); //存入曲线点 curvePoint.push_back(tempP); } } } //三次贝塞尔曲线 float bezier3funcX(float uu,CvPoint *controlP){ float part0 = controlP[0].x * uu * uu * uu; float part1 = 3 * controlP[1].x * uu * uu * (1 - uu); float part2 = 3 * controlP[2].x * uu * (1 - uu) * (1 - uu); float part3 = controlP[3].x * (1 - uu) * (1 - uu) * (1 - uu); return part0 + part1 + part2 + part3; } float bezier3funcY(float uu,CvPoint *controlP){ float part0 = controlP[0].y * uu * uu * uu; float part1 = 3 * controlP[1].y * uu * uu * (1 - uu); float part2 = 3 * controlP[2].y * uu * (1 - uu) * (1 - uu); float part3 = controlP[3].y * (1 - uu) * (1 - uu) * (1 - uu); return part0 + part1 + part2 + part3; }
相关阅读:
【Proteogenomis】工具评述
 【6】蛋白质组学鉴定定量软件之MaxQuant
【5】蛋白质组学鉴定定量软件之PD
【4】蛋白质组学鉴定软件之MSGFPlus
堆叠式降噪自动编码器（SDA）
机器学习
 决策优先级的算法
 软件的安装
 一定要读高水平的论文，论文的档次越高，收益越深
 五种世界的顶级定律
原文地址：https://www.cnblogs.com/hnfxs/p/3148743.html