优化了的过关键点的光滑曲线拟合算法

这个是我一个数学老师(教授,数学高手,经常自己做算法)给我的例子,用于多个离散点拟合光滑曲线的,他优化了追赶法,这个例子适用于闭合和不闭合两种情况。当时由于工程情况,写的急,代码不好看,但是很好用。为了方便传递参数,我做了一个链表,用时候根据自己情况可以修改,核心算法不动即可。

class CFoldPoint

{public:

    double X;    double Y;

};

typedef CTypedPtrList CFoldPointList;

typedef CArray CDoubleArray;

三个函数,SPLine 调用另外两个。用时候直接调用SPLine函数,入口pList是已知离散点链表,pDestList是生成的点的链表。SM是在两个点中间插入点的数目,continue=0是采样点无规律,要求生成闭合曲线。1是采样点x坐标连续 2是y连续

void ZG(CDoubleArray *A,CDoubleArray *B,CDoubleArray *C,CDoubleArray *G,int &LOGI)
{
 //追赶法
 register long I;
 int N;
 N=A->GetSize();
 if(LOGI==0)
 {
  (*C)[0]=(*C)[0]/(*B)[0];
   for(I=1;I    {
    (*B)[I]=(*B)[I]-(*A)[I]*(*C)[I-1];
   (*C)[I]=(*C)[I]/(*B)[I];
  }
  (*A)[0]=0.;
   (*C)[N-1]=0.;
   LOGI=1;
  }
  (*G)[0]=(*G)[0]/(*B)[0];
  for(I=1;I   {
   (*G)[I]=((*G)[I]-(*A)[I]*(*G)[I-1])/(*B)[I];
  }
  for(I=N-2;I>-1;I--)//DO 30 I=N-1,1,-1
  {
   (*G)[I]=(*G)[I]-(*C)[I]*(*G)[I+1];
  }
  return;
}

void SPLine4(CDoubleArray *X,CDoubleArray *Y,double &XI,double&YI,CDoubleArray *A,CDoubleArray *B,CDoubleArray *C,CDoubleArray *G,int &LOGI,int MD)
{
 
 register long I;
 double W1,W2,H;
 int N=X->GetSize();
 
 if(LOGI==0)
 {
  for(I=1;I   {
   (*B)[I]=(*X)[I]-(*X)[I-1];
   (*C)[I]=((*Y)[I]-(*Y)[I-1])/(*B)[I];
  }
  for(I=1;I   {
   (*A)[I]=(*B)[I]+(*B)[I+1];
   (*G)[I]=6.*((*C)[I+1]-(*C)[I])/(*A)[I];
   (*A)[I]=(*B)[I]/(*A)[I];
  }
  for(I=1;I   {
   (*C)[I]=1.-(*A)[I];
   (*B)[I]=2.;
  }
  (*B)[0]=2.;
  (*B)[N-1]=2.;
  if(MD==3)
  {
   (*C)[0]=-1.;
   (*A)[N-1]=-1.;
   (*A)[0]=0.;
   (*C)[N-1]=0.;
  }
  ZG(A,B,C,G,LOGI);
 }
 for(I=1;I  {
  if(XI>=(*X)[I-1] && XI<=(*X)[I])//GE LE
  {
   H=(*X)[I]-(*X)[I-1];
   W1=(*X)[I]-XI;
   W2=XI-(*X)[I-1];
   YI=W1*W1*W1*(*G)[I-1]/6./H;
   YI=YI+W2*W2*W2*(*G)[I]/6./H;
   YI=YI+W1*((*Y)[I-1]-(*G)[I-1]*H*H/6.)/H;
   YI=YI+W2*((*Y)[I]-(*G)[I]*H*H/6.)/H;
  }
 }
}
void SPLine(CFoldPointList *pList,CFoldPointList *pDestList,int SM,int Continue=0)
{
 CFoldPoint *pFoldHead,*pFoldTail;
 POSITION pos;
 CDoubleArray A,B,C,G,X,Y,T;
 double XI,YI,XX,YY;
 register long i;
 long N;
 int LOGI;
 long RealSM;
 long Bei,Yu;
 CFoldPoint *pFold;
  file://赋初值
 N=pList->GetCount();
 A.SetSize(N);
 B.SetSize(N);
 C.SetSize(N);
 G.SetSize(N);
 X.SetSize(N);
 Y.SetSize(N);
 T.SetSize(N);
 RealSM=(N-1)*SM+N;
 pos=pList->GetHeadPosition();
 for(i=0;i
 {
  pFold=pList->GetNext(pos);
  X[i]=pFold->X;
  Y[i]=pFold->Y;
 }
 
 pFoldHead=pList->GetHead();
 pFoldTail=pList->GetTail();
 if(Continue==0)//pFoldHead->X==pFoldTail->X && pFoldHead->Y==pFoldTail->Y)
 {  file://闭合
  T[0]=0;
  for(i=0;i
  {
   T[i+1]=T[i]+CalculateDistance(X[i],Y[i],X[i+1],Y[i+1])+0.000000001;
  }
  LOGI=0;
  YI=0;
  for(i=0;i
  {
   Bei=i/(SM+1);
   Yu=i%(SM+1);
   if(Yu!=0)
   {
    XI=T[Bei]+(T[Bei+1]-T[Bei])/(SM+1)*Yu;
    SPLine4(&T,&Y,XI,YI,&A,&B,&C,&G,LOGI,3);
    YY=YI;//+Y[Bei];
   }
   else
   {
    YY=Y[Bei];
   }
   pFold=new CFoldPoint;
   pFold->Y=YY;
   pDestList->AddTail(pFold);
  }
  LOGI=0;
  YI=0;
  pos=pDestList->GetHeadPosition();
  for(i=0;i
  {
   Bei=i/(SM+1);
   Yu=i%(SM+1);
   if(Yu!=0)
   {
    XI=T[Bei]+(T[Bei+1]-T[Bei])/(SM+1)*Yu;
    SPLine4(&T,&X,XI,YI,&A,&B,&C,&G,LOGI,3);
    YY=YI;//+X[Bei];
   }
   else
   {
    YY=X[Bei];
   }
   pFold=pDestList->GetNext(pos);
   pFold->X=YY;
  }
 }
 else if(Continue==1)
 {
   file://x连续
  LOGI=0;
  YI=0;
  for(i=0;i
  {
   Bei=i/(SM+1);
   Yu=i%(SM+1);
   if(Yu!=0)
   {
    XI=X[Bei]+(X[Bei+1]-X[Bei])/(SM+1)*Yu;
    SPLine4(&X,&Y,XI,YI,&A,&B,&C,&G,LOGI,3);
    XX=XI;
    YY=YI;
   }
   else
   {
    XX=X[Bei];
    YY=Y[Bei];
   }
   pFold=new CFoldPoint;
   pFold->X=XX;
   pFold->Y=YY;
   pDestList->AddTail(pFold);
  }
 }
 else
 {
   file://y连续
  LOGI=0;
  YI=0;
  for(i=0;i
  {
   Bei=i/(SM+1);
   Yu=i%(SM+1);
   if(Yu!=0)
   {
    XI=Y[Bei]+(Y[Bei+1]-Y[Bei])/(SM+1)*Yu;
    SPLine4(&Y,&X,XI,YI,&A,&B,&C,&G,LOGI,3);
    XX=YI;
    YY=XI;
   }
   else
   {
    XX=X[Bei];
    YY=Y[Bei];
   }
   pFold=new CFoldPoint;
   pFold->X=XX;
   pFold->Y=YY;
   pDestList->AddTail(pFold);
  }
 }
 return;
}

  • 1
    点赞
  • 2
    评论
  • 3
    收藏
  • 一键三连
    一键三连
  • 扫一扫,分享海报

相关推荐
matlab离散连成曲线-intersections.m 本帖最后由 kastin 于 2012-12-29 11:47 编辑 引言     曾经思考曲面求交,结果发现是学术界一个难题,并且也想出一个当前广泛使用方法原理一样近似解法(追踪法)。当然网上也有很多方法,只不那些方法非常粗糙,无非就是meshgrid出离散网格,比较两曲面在某位置坐标是否在某一精度范围内,然后标记显示之。这个方法仅仅当离散网格非常细时候才比较精确。除此之外,还有个非常严重问题:上面“精度范围”不是你随心所欲给,而且也没规律寻找,当给得不恰当时候,在格处两曲面作比较,会出很多个符合要求,或者一个也没有。这样就会使得交线非常曲折,甚至断裂等,严重影响精确度。 ———————————————————分割线————————————————————————     当然,既然有曲面求交,那么也有曲线求交,其基本结构就是两曲线求交。只是曲线求交问题,事先得澄清一些注意:     1. 数学分析层面求两曲线,其实就是方程组求解;     2. “曲线”概念包括“直线”(处处曲率半径为无穷大);     3. Matlab离散 矩阵运算,因此所有运算都是基于离散,因而这里曲线并不是绝对光滑。     4. 近似试探与未知函数表达式。 对于1,我想说是,如果你想要求得两曲线精确交,并且一个不漏,那就直接求解方程组,不用看本帖下文; 对于2,直线在Matlab里面是两个确定,因此交如果是一段线(无穷个情况,可能只是显示两端为交; 对于3,很简单例子,参数方程 x=cos,y=sin 在数学分析(即连续空间)层面上是个圆,但是如果你在离散t时候,间距比较大,那么最后Matlab绘制图像不是圆,而是正多边形。因此,此时我们讨论曲线是这个离散连线图形与其他图形,而非圆与其他交。这也是我在标题中加离散连成”修饰词,防止被误会。 对于4,既然是求曲线,那么本方法可以作为求方程组近似解。当然,如果离散够多,解精确度可以保证,不不能保证一个不漏。另外就是,对于一组离散构成曲线,很难知道它们解析表达式,因此想通非线性方程组求解方法来求交,就不大可能(不你可以用曲线拟合出函数解析式),因此,本帖方法将会是一个较为有效求交方法。     废话那么多,下面就说说曲线求交方法吧。除求解方程组,很多人想到方法就是“离散 判断距离是否足够接近”,这个方法原理跟引言中曲面求交方法是一样。因此缺也是一样——太粗糙。网上这种方法代码也很多,这里就不上。 下面将阐述我方法以及给出例子代码。     我有两种思路,一种是高级绘图层面(不涉及到底层操作),一种是底层。我只给出第一种代码,因为我不会底层操作。     思路一:既然matlab曲线绘图是通有序离散依次连线形成,也就是说,通“以直代曲”程,那么曲线无非就是离散(结)或者两线段交。这比上面直接用交附近替代交方法要精确得多。而两直线交很容易求,只要知道四个坐标,那么交精确坐标自然可以表示出来。这就是求交原理。只是还有一些细节处理和要注意地方,我会留到后面再详细说。     思路二:仔细观察两曲线特性,很容易发现,其实交就是操作系统底层绘图重叠那些像素。因此,只要给要绘制像素做个标记,将那些重合突出显示(比如换个颜色),那么就相当于显示出交。这种方法由于是本质性,因此不会遗漏任何交,而且精确度极高,适用范围广。Matlab提供plot plot3 surf等绘图函数都属于高级绘图,底层绘图(或称低级绘图)只有line surface以及patch等少数函数。但是,这里“底层”并非真正底层,因为它还是经封装,而C MFC里面直接用刷子绘图,那才是依靠操作系统完成真正“底层”绘图操作(包括所有窗口都是操作系统绘制)。这里扯远,想要说明就是底层绘图概念而已。只是我不会用matlab实现这些底层绘图。     上面说思路,下面就详细说说一些注意和需要处理细节。     为算法健壮性,就必须考虑各种奇异情况,防止bug。我们要考虑曲线有分支(很多代数曲线是这样,代数几何里面研究东西)、间断跳跃(有绝对值函数或者存在渐近线情况)、首尾是交、在切相交,等等这些情况。而且对于定位交处附近四个最近端也是个问题(因为这里存在一个情况,如果曲线1上一条线段与曲线2上两条
©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值