自适应辛普森迭代求曲线的长度

最近工作都在研究几何相关的算法：需要计算贝塞尔，样条等曲线的弧长。贴一块最近写的递归计算长度的代码。

递归计算弧长的代码：

a,b=[0,1];  //
der1=ArcLenDer(0);//曲线导矢的模
der2= ArcLenDer(0.5);
der3=ArcLenDer(1);
//自适应辛普森计算公式
double simpsonCalLength(double SP,double EP,double eps,double A,der1,der2,der3) //der1,der2端点导矢的模
{
double d1=der1;
double DerMP= der2;  //中间点的导矢
double d2=der3;
  double MP=SP+(EP-SP)/2;//中点

//计算长度
double mid41= SP +( MP - SP)/2;  //上半段中点
double der41=ArcLenDer(mid41); //上半段中点导矢的模
double mid42= MP +( EP - MP)/2;  
double der42=ArcLenDer(mid42);     //下半段中点导矢的模

double L=LengthSimpson(SP,MP,d1, der41，DerMP);
double R=LengthSimpson(Mp,EP,DerMP, der42，d2);

//判断是否满足精度; eps =0.000001
if (abs(L + R - A) < eps)
    return L + R;
else
return simpsonCalLength(SP,MP,eps,L,d1,der41,DerMP)+ simpsonCalLength(MP,EP,eps,R DerMP,der42,d2);
}
//单段辛普森计算公式

double LengthSimpson(double startPoint , double endpoint ,double de1,double de2，double de3) //
{
  return (de1+ 4*de2+ de3) *(endpoint-startPoint)/6;
}

代码解释，计算参数曲线的弧长,参数曲线参数从[0:1]，计算0,0.5,1处的导矢，计算单段曲线长度时，辛普森公式要用到。（注：计算对称曲线，可以计算[0:0.5],然后曲线的长度再乘以2） 然后计算前一半和后一半的长度，此时递归重复利用之前算的0,0.5,1处的导矢（对于计算样条这种效率变得很高通过递归的方式）。满足精度abs(L + R - A) < eps就取消递归，不满足就一层一层的计算下去。实际中算法效率是十分重要的，算法实现是一方面，优化也是重要的方面。 （欢迎交流 ，微博：雪碧可口）