C++实现二次、三次B样条曲线
原文：Bezier曲线、B样条和NURBS的基本概念
下面是一个有四个控制点的Bezier曲线：
可以通过改变一个控制点的位置来改变曲线的形状，比如将上图曲线中左边第二个控制点往上移，就可以得到下面的曲线:
可以看到，这种曲线生成方式比较直观和灵活，我只需要放置控制点，然后调整控制点的位置来得到想要的曲线，这就避免了和复杂的数学方程打交道，岂不快哉？
Bezier曲线、B样条和NURBS都是根据控制点来生成曲线的，那么他们有什么区别了？简单来说，就是:
Bezier曲线中的每个控制点都会影响整个曲线的形状，而B样条中的控制点只会影响整个曲线的一部分，显然B样条提供了更多的灵活性；
Bezier和B样条都是多项式参数曲线，不能表示一些基本的曲线，比如圆，所以引入了NURBS，即非均匀有理B样条来解决这个问题；
B样条克服了Bezier曲线的一些缺点，Bezier曲线的每个控制点对整条曲线都有影响，也就是说，改变一个控制点的位置，整条曲线的形状都会发生变化，而B样条中的每个控制点只会影响曲线的一段参数范围，从而实现了局部修改；
以下内容原文：二次与三次B样条曲线c++实现
B样条曲线从Bezier曲线演变而来，了解B样条曲线首先得了解Bezier曲线。对于平面上的三个点P0，P1，P2，其坐标分别是(x0,y0)、(x1,y1)、(x2,y2)。二次Bezier曲线用一条抛物线进行拟合，其参数方程是
二次Bezier曲线有如下要求：（1）t=0时，曲线经过P0，并与P0P1相切；（2）t=1时，曲线经过P2，并与P1P2相切。即有
其中
，称为二次Bezier曲线基函数。
每3个离散点形成一条Bezier曲线，由于每条Bezier都经过端点，导致分段的Bezier曲线在端点处难以平滑过渡。二次B样条曲线克服了这个缺点，把端点移到线段中点（如下图所示），这样就能保证各段曲线在连接处能够一阶导数连续。
二次B样条曲线实现，只需将曲线参数t划分成k等分，t从0开始取值，间隔dt，直至k*dt。如果想让整条曲线两端与起始点P0和终止点Pn重合，只需以P0和Pn为中点，构造新点PP1=2*P0-P1，与PP2=2*Pn-Pn-1，替换掉P0与Pn即可。
三次B样条曲线通过样条基函数可以得到如下形式：
三次B样条具有以下物理意义，曲线起点S位于三角形P0P1P2的中线P1M1上，距离P1点1/3倍P1M1处；曲线中点E位于三角形P1P2P3的中线P2M2上，距离1/3倍P2M2处；曲线起点切线平行于P0P2，终点切线平行于P1P3.
三次B样条曲线要想让曲线两端与起始端P0与Pn重合，只需构造新点PP1=2*P0-P1，与PP2=2*Pn-Pn-1，分别加到P0之前与Pn之后即可。由此，参与计算点的增加了2个（注意，二次B样条是替换不是增加）。
程序实现了二次与三次B样条曲线，封装成了BSpline类
BSpine.h
1 #pragma once 2 //#include "position.h" 3 4 typedef struct tagPosition 5 { 6 double x; 7 double y; 8 tagPosition(double _x,double _y) { x=_x; y=_y;} 9 tagPosition() {}; 10 bool operator==(const tagPosition & pt) { return (x==pt.x && y==pt.y);} 11 } CPosition; 12 13 class CBSpline 14 { 15 public: 16 CBSpline(void); 17 ~CBSpline(void); 18 19 void TwoOrderBSplineSmooth(CPosition *pt,int Num); 20 void TwoOrderBSplineInterpolatePt(CPosition *&pt,int &Num,int *InsertNum); 21 double F02(double t); 22 double F12(double t); 23 double F22(double t); 24 25 void ThreeOrderBSplineSmooth(CPosition *pt,int Num); 26 void ThreeOrderBSplineInterpolatePt(CPosition *&pt,int &Num,int *InsertNum); 27 double F03(double t); 28 double F13(double t); 29 double F23(double t); 30 double F33(double t); 31 };
BSpine.cpp
1 //**************************** BSpline.cpp *********************************** 2 // 包含功能：二次B样条平滑，三次B样条平滑；二次B样条平滑后节点插值 3 // 4 // 作者： 蒋锦朋 1034378054@qq.com 5 // 单位： 中国地质大学（武汉） 地球物理与空间信息学院 6 // 日期： 2014/12/03 7 //************************************************************************************* 8 #include "StdAfx.h" 9 #include "BSpline.h" 10 11 12 CBSpline::CBSpline(void) 13 { 14 } 15 16 17 CBSpline::~CBSpline(void) 18 { 19 } 20 //====================================================================== 21 // 函数功能： 二次B样条平滑，把给定的点，平滑到B样条曲线上，不增加点的数目 22 // 输入参数： *pt ：给定点序列，执行完成后，会被替换成新的平滑点 23 // Num：点个数 24 // 返回值： 无返回值 25 26 // 编辑日期： 2014/12/03 27 //====================================================================== 28 void CBSpline::TwoOrderBSplineSmooth(CPosition *pt,int Num) 29 { 30 CPosition *temp=new CPosition[Num]; 31 for(int i=0;i<Num;i++) 32 temp[i]=pt[i]; 33 34 temp[0].x=2*temp[0].x-temp[1].x; // 将折线两端点换成延长线上两点 35 temp[0].y=2*temp[0].y-temp[1].y; 36 37 temp[Num-1].x=2*temp[Num-1].x-temp[Num-2].x; 38 temp[Num-1].y=2*temp[Num-1].y-temp[Num-2].y; 39 40 CPosition NodePt1,NodePt2,NodePt3; 41 double t; 42 for(int i=0;i<Num-2;i++) 43 { 44 NodePt1=temp[i]; NodePt2=temp[i+1]; NodePt3=temp[i+2]; 45 if(i==0) // 第一段取t=0和t=0.5点 46 { 47 t=0; 48 pt[i].x=F02(t)*NodePt1.x+F12(t)*NodePt2.x+F22(t)*NodePt3.x; 49 pt[i].y=F02(t)*NodePt1.y+F12(t)*NodePt2.y+F22(t)*NodePt3.y; 50 t=0.5; 51 pt[i+1].x=F02(t)*NodePt1.x+F12(t)*NodePt2.x+F22(t)*NodePt3.x; 52 pt[i+1].y=F02(t)*NodePt1.y+F12(t)*NodePt2.y+F22(t)*NodePt3.y; 53 }else if(i==Num-3) // 最后一段取t=0.5和t=1点 54 { 55 t=0.5; 56 pt[i+1].x=F02(t)*NodePt1.x+F12(t)*NodePt2.x+F22(t)*NodePt3.x; 57 pt[i+1].y=F02(t)*NodePt1.y+F12(t)*NodePt2.y+F22(t)*NodePt3.y; 58 t=1; 59 pt[i+2].x=F02(t)*NodePt1.x+F12(t)*NodePt2.x+F22(t)*NodePt3.x; 60 pt[i+2].y=F02(t)*NodePt1.y+F12(t)*NodePt2.y+F22(t)*NodePt3.y; 61 }else // 中间段取t=0.5点 62 { 63 t=0.5; 64 pt[i+1].x=F02(t)*NodePt1.x+F12(t)*NodePt2.x+F22(t)*NodePt3.x; 65 pt[i+1].y=F02(t)*NodePt1.y+F12(t)*NodePt2.y+F22(t)*NodePt3.y; 66 } 67 } 68 delete []temp; 69 } 70 71 //================================================================ 72 // 函数功能： 二次B样条拟合,在节点之间均匀插入指定个数点 73 // 输入参数： *pt ：给定点序列，执行完成后，会被替换成新的数据点 74 // Num：节点点个数 75 // *InsertNum: 节点之间需要插入的点个数指针 76 // 返回值： 无返回值 77 // 78 // 编辑日期： 2014/12/07 79 //================================================================= 80 void CBSpline::TwoOrderBSplineInterpolatePt(CPosition *&pt,int &Num,int *InsertNum) 81 { 82 if(pt==NULL || InsertNum==NULL) return; 83 84 int InsertNumSum=0; // 计算需要插入的点总数 85 for(int i=0;i<Num-1;i++) InsertNumSum+=InsertNum[i]; 86 87 CPosition *temp=new CPosition[Num]; // 二次B样条不需要增加点数，需要将首尾点替换掉 88 for(int i=0;i<Num;i++) 89 temp[i]=pt[i]; 90 91 temp[0].x=2*temp[0].x-temp[1].x; // 将折线两端点换成延长线上两点 92 temp[0].y=2*temp[0].y-temp[1].y; 93 94 temp[Num-1].x=2*temp[Num-1].x-temp[Num-2].x; 95 temp[Num-1].y=2*temp[Num-1].y-temp[Num-2].y; 96 97 delete []pt; // 点数由原来的Num个增加到Num+InsertNumSum个，删除旧的存储空间，开辟新的存储空间 98 99 pt=new CPosition[Num+InsertNumSum]; 100 101 CPosition NodePt1,NodePt2,NodePt3,NodePt4; // 两节点间均匀插入点，需要相邻两段样条曲线,因此需要四个节点 102 103 double t; 104 int totalnum=0; 105 for(int i=0;i<Num-1;i++) // 每条线段均匀插入点 106 { 107 if(i==0) // 第一段只需计算第一条样条曲线，无NodePt1 108 { 109 NodePt2=temp[i]; NodePt3=temp[i+1]; NodePt4=temp[i+2]; 110 111 double dt=0.5/(InsertNum[i]+1); 112 for(int j=0;j<InsertNum[i]+1;j++) 113 { 114 t=0+dt*j; 115 pt[totalnum].x=F02(t)*NodePt2.x+F12(t)*NodePt3.x+F22(t)*NodePt4.x; 116 pt[totalnum].y=F02(t)*NodePt2.y+F12(t)*NodePt3.y+F22(t)*NodePt4.y; 117 totalnum++; 118 } 119 }else if(i==Num-2) // 最后一段只需计算最后一条样条曲线，无NodePt4 120 { 121 NodePt1=temp[i-1]; NodePt2=temp[i]; NodePt3=temp[i+1]; 122 123 double dt=0.5/(InsertNum[i]+1); 124 for(int j=0;j<InsertNum[i]+2;j++) 125 { 126 t=0.5+dt*j; 127 pt[totalnum].x=F02(t)*NodePt1.x+F12(t)*NodePt2.x+F22(t)*NodePt3.x; 128 pt[totalnum].y=F02(t)*NodePt1.y+F12(t)*NodePt2.y+F22(t)*NodePt3.y; 129 totalnum++; 130 } 131 }else 132 { 133 NodePt1=temp[i-1],NodePt2=temp[i]; NodePt3=temp[i+1]; NodePt4=temp[i+2]; // NodePt1,2,3计算第一条曲线，NodePt2,3,4计算第二条曲线 134 135 int LeftInsertNum,RightInsertNum; // 计算线段间左右曲线段上分别需要插入的点数 136 double rightoffset=0; // 左边曲线段从t=0.5开始，又边曲线段从t=rightoffset开始 137 double Leftdt=0,Rightdt=0; // 左右曲线取点t步长 138 if(InsertNum[i]==0 ) 139 { 140 LeftInsertNum=0; 141 RightInsertNum=0; 142 }else if(InsertNum[i]%2==1) // 插入点数为奇数，左边曲线段插入点个数比右边多一点 143 { 144 RightInsertNum=InsertNum[i]/2; 145 LeftInsertNum=RightInsertNum+1; 146 Leftdt=0.5/(LeftInsertNum); 147 Rightdt=0.5/(RightInsertNum+1); 148 rightoffset=Rightdt; 149 }else // 插入点数为偶数，左右边曲线段插入个数相同 150 { 151 RightInsertNum=InsertNum[i]/2; 152 LeftInsertNum=RightInsertNum; 153 Leftdt=0.5/(LeftInsertNum+0.5); 154 Rightdt=0.5/(RightInsertNum+0.5); 155 rightoffset=Rightdt/2; 156 } 157 158 for(int j=0;j<LeftInsertNum+1;j++) 159 { 160 t=0.5+Leftdt*j; 161 pt[totalnum].x=F02(t)*NodePt1.x+F12(t)*NodePt2.x+F22(t)*NodePt3.x; 162 pt[totalnum].y=F02(t)*NodePt1.y+F12(t)*NodePt2.y+F22(t)*NodePt3.y; 163 totalnum++; 164 } 165 166 for(int j=0;j<RightInsertNum;j++) 167 { 168 t=rightoffset+Rightdt*j; 169 pt[totalnum].x=F02(t)*NodePt2.x+F12(t)*NodePt3.x+F22(t)*NodePt4.x; 170 pt[totalnum].y=F02(t)*NodePt2.y+F12(t)*NodePt3.y+F22(t)*NodePt4.y; 171 totalnum++; 172 } 173 } 174 } 175 delete []temp; 176 Num=Num+InsertNumSum; 177 178 } 179 //================================================================ 180 // 函数功能： 二次样条基函数 181 // 182 // 编辑日期： 2014/12/03 183 //================================================================ 184 double CBSpline::F02(double t) 185 { 186 return 0.5*(t-1)*(t-1); 187 } 188 double CBSpline::F12(double t) 189 { 190 return 0.5*(-2*t*t+2*t+1); 191 } 192 double CBSpline::F22(double t) 193 { 194 return 0.5*t*t; 195 } 196 //======================================================================== 197 // 函数功能： 三次B样条平滑，把给定的点，平滑到B样条曲线上，不增加点的数目 198 // 输入参数： *pt ：给定点序列，执行完成后，会被替换成新的平滑点 199 // Num：点个数 200 // 返回值： 无返回值 201 // 202 // 编辑日期： 2014/12/03 203 //======================================================================== 204 void CBSpline::ThreeOrderBSplineSmooth(CPosition *pt,int Num) 205 { 206 CPosition *temp=new CPosition[Num+2]; 207 for(int i=0;i<Num;i++) 208 temp[i+1]=pt[i]; 209 210 temp[0].x=2*temp[1].x-temp[2].x; // 将折线延长线上两点加入作为首点和尾点 211 temp[0].y=2*temp[1].y-temp[2].y; 212 213 temp[Num+1].x=2*temp[Num].x-temp[Num-1].x; 214 temp[Num+1].y=2*temp[Num].y-temp[Num-1].y; 215 216 CPosition NodePt1,NodePt2,NodePt3,NodePt4; 217 double t; 218 for(int i=0;i<Num-1;i++) 219 { 220 NodePt1=temp[i]; NodePt2=temp[i+1]; NodePt3=temp[i+2]; NodePt4=temp[i+3]; 221 222 if(i==Num-4) // 最后一段取t=0.5和t=1点 223 { 224 t=0; 225 pt[i].x=F03(t)*NodePt1.x+F13(t)*NodePt2.x+F23(t)*NodePt3.x+F33(t)*NodePt4.x; 226 pt[i].y=F03(t)*NodePt1.y+F13(t)*NodePt2.y+F23(t)*NodePt3.y+F33(t)*NodePt4.y; 227 t=1; 228 pt[i+1].x=F03(t)*NodePt1.x+F13(t)*NodePt2.x+F23(t)*NodePt3.x+F33(t)*NodePt4.x; 229 pt[i+1].y=F03(t)*NodePt1.y+F13(t)*NodePt2.y+F23(t)*NodePt3.y+F33(t)*NodePt4.y; 230 }else // 中间段取t=0.5点 231 { 232 t=0; 233 pt[i].x=F03(t)*NodePt1.x+F13(t)*NodePt2.x+F23(t)*NodePt3.x+F33(t)*NodePt4.x; 234 pt[i].y=F03(t)*NodePt1.y+F13(t)*NodePt2.y+F23(t)*NodePt3.y+F33(t)*NodePt4.y; 235 } 236 } 237 delete []temp; 238 } 239 240 //================================================================ 241 // 函数功能： 三次B样条拟合,在节点之间均匀插入指定个数点 242 // 输入参数： *pt ：给定点序列，执行完成后，会被替换成新的数据点 243 // Num：节点点个数 244 // *InsertNum: 节点之间需要插入的点个数指针 245 // 返回值： 无返回值 246 // 247 // 编辑日期： 2014/12/07 248 //================================================================= 249 void CBSpline::ThreeOrderBSplineInterpolatePt(CPosition *&pt,int &Num,int *InsertNum) 250 { 251 if(pt==NULL || InsertNum==NULL) return; 252 253 int InsertNumSum=0; // 计算需要插入的点总数 254 for(int i=0;i<Num-1;i++) InsertNumSum+=InsertNum[i]; 255 256 CPosition *temp=new CPosition[Num+2]; 257 for(int i=0;i<Num;i++) 258 temp[i+1]=pt[i]; 259 260 temp[0].x=2*temp[1].x-temp[2].x; // 将折线延长线上两点加入作为首点和尾点 261 temp[0].y=2*temp[1].y-temp[2].y; 262 263 temp[Num+1].x=2*temp[Num].x-temp[Num-1].x; 264 temp[Num+1].y=2*temp[Num].y-temp[Num-1].y; 265 266 CPosition NodePt1,NodePt2,NodePt3,NodePt4; 267 double t; 268 269 delete []pt; // 点数由原来的Num个增加到Num+InsertNumSum个，删除旧的存储空间，开辟新的存储空间 270 271 pt=new CPosition[Num+InsertNumSum]; 272 273 int totalnum=0; 274 for(int i=0;i<Num-1;i++) // 每条线段均匀插入点 275 { 276 NodePt1=temp[i]; NodePt2=temp[i+1]; NodePt3=temp[i+2]; NodePt4=temp[i+3]; 277 double dt=1.0/(InsertNum[i]+1); 278 279 for(int j=0;j<InsertNum[i]+1;j++) 280 { 281 t=dt*j; 282 pt[totalnum].x=F03(t)*NodePt1.x+F13(t)*NodePt2.x+F23(t)*NodePt3.x+F33(t)*NodePt4.x; 283 pt[totalnum].y=F03(t)*NodePt1.y+F13(t)*NodePt2.y+F23(t)*NodePt3.y+F33(t)*NodePt4.y; 284 totalnum++; 285 } 286 287 if(i==Num-2){ // 最后一个尾点 288 t=1; 289 pt[totalnum].x=F03(t)*NodePt1.x+F13(t)*NodePt2.x+F23(t)*NodePt3.x+F33(t)*NodePt4.x; 290 pt[totalnum].y=F03(t)*NodePt1.y+F13(t)*NodePt2.y+F23(t)*NodePt3.y+F33(t)*NodePt4.y; 291 totalnum++; 292 } 293 } 294 295 delete []temp; 296 Num=Num+InsertNumSum; 297 298 } 299 300 //================================================================ 301 // 函数功能： 三次样条基函数 302 // 303 // 编辑日期： 2014/12/03 304 //================================================================ 305 double CBSpline::F03(double t) 306 { 307 return 1.0/6*(-t*t*t+3*t*t-3*t+1); 308 } 309 double CBSpline::F13(double t) 310 { 311 return 1.0/6*(3*t*t*t-6*t*t+4); 312 } 313 double CBSpline::F23(double t) 314 { 315 return 1.0/6*(-3*t*t*t+3*t*t+3*t+1); 316 } 317 double CBSpline::F33(double t) 318 { 319 return 1.0/6*t*t*t; 320 }
程序调用
1 #include "stdafx.h" 2 #include "math.h" 3 #include "BSpline.h" 4 5 int _tmain(int argc, _TCHAR* argv[]) 6 { 7 int num=8; 8 double x[8]={9.59,60.81,105.57,161.59,120.5,100.1,50.0,10.0}; 9 double y[8]={61.97,107.13,56.56,105.27,120.5,150.0,110.0,180.0}; 10 11 CPosition *testpt=new CPosition[num]; 12 for(int i=0;i<num;i++) testpt[i]=CPosition(x[i],y[i]); 13 14 int *Intnum=new int[num-1]; 15 for(int i=0;i<num-1;i++){ 16 Intnum[i]=10; // 每一个样条曲线内插入10个点 17 } 18 19 int num2=num; 20 CBSpline bspline; 21 bspline.TwoOrderBSplineInterpolatePt(testpt,num2,Intnum); // 二次B样条曲线 22 //bspline.ThreeOrderBSplineInterpolatePt(testpt,num2,Intnum); // 三次B样条曲线 23 24 //bspline.TwoOrderBSplineSmooth(testpt,num2); // 二次B样条平滑 25 //bspline.ThreeOrderBSplineSmooth(testpt,num2); // 三次B样条平滑 26 delete Intnum; 27 28 29 FILE *fp_m_x = fopen("Bspline_test_x.txt", "wt"); 30 FILE *fp_m_y = fopen("Bspline_test_y.txt", "wt"); 31 for (int i = 0; i < num2; i++){ 32 fprintf(fp_m_x, "%lf\n", testpt[i].x); 33 fprintf(fp_m_y, "%lf\n", testpt[i].y); 34 } 35 fclose(fp_m_x); 36 fclose(fp_m_y); 37 38 return 0; 39 }
