平滑曲线算法研究

一  样条概述
在绘图术语中样条是通过一组指定点集而生成平滑曲线的柔性带 。
术语 样条曲线 spline curve
绘制样条曲线的方法是给定一组称为控制点的坐标点，可以得到一条样条曲线。
样条曲线分为:
1 插值样条曲线(interpolate)  生成的样条曲线通过这组控制点。（这是我们详细研究的）2 逼近样条曲线(approximate)  生面的样条曲线不通过或通过部分控制点。
   贝塞尔曲线
   B样条曲线
通俗的说构造样条就是 通过给定的一组点（一组 X Y值），比如说 以下四个点 (X在1 到100之间取值)
A(10,15 ) B(28,78) C(68,32) D(76,12)
计算出一个函数Y ＝ F(X)
然后把X 从 1 到100 全代入这个函数，就得到了 100个点，用这100个点画曲线，会平滑许多。
只不过，生成的函数更加的复杂。当控制点很多时，比如有100个，计算函数时，是把控制点分段，比如说 前4个点生成一个Y = F1(X),下四个点生成另一个函数Y = F2(X),如何让这些分段的曲线的连接处保持平滑呢，这就要求各种连续性条件 (continuity condition)
通过在曲线段的公共部分匹配连接的参数导数，从建立参数连续性(parametric continuity)。
二 三次样条插值
   三次样条生成的分段函数是三次多项式，这在灵活性和计算速度之间提供了一个合理的折中方案，与更高层次多项式相比，三次样条只需要较少的计算和存储空间，并且比较稳定，与低次多项式相比，三次样条在模拟任意曲线形状时显得更加灵活。
   三次样条插值
   在数学上,光滑程度的定量描述是：函数(曲线)的k阶导数存在且连续,则称该曲线具有k阶光滑性.自然,阶数越高光滑程度越好.于是,分段线性插值具有零阶光滑性,也就是不光滑；分段三次埃尔米特插值具有一阶光滑性.仅有这些光滑程度,在工程设计和机械加工等实际中是不够的.提高分段函数如多项式函数的次数,可望提高整体曲线的光滑程度.但是,是否存在较低次多项式达到较高阶光滑性的方法？三次样条插值就是一个很好的例子.
    样条曲线本身就来源于飞机、船舶等外形曲线设计中所用的绘图工具.在工程实际中,要求这样的曲线应该具有连续的曲率,也就是连续的二阶导数.值得注意的是分段插值曲线的光滑性关键在于段与段之间的衔接点(节点)处的光滑性.
    三次样条函数记为s(x),它是定义在区间[a,b]上的函数,满足 
    1)s(x)在每一个小区间[Xi-1,Xi]上是一个三次多项式函数；
    2)在整个区间[a,b]上,其二阶导数存在且连续.即在每个节点处的二阶导数连续.
三次样条插值问题的提法:给定函数f(x)在n+1个节点x0,x1,...,xn处的函数值为y0,y1,...,yn,求一个三次样条函数s(x),使其满足
             s(xi)=yi ,       i=0,1,…,n.
    如何确定三次样条函数在每一个小区间上的三次多项式函数的系数呢？这是一个比较复杂的问题,这里只介绍确定系数的思想.
   分段线性插值在每一段的线性函数的两个参数,是由两个方程(两个端点处的函数值为给定值)唯一确定；对于三次样条插值呢,每一个区间上的三次函数有四个参数,而在该区间上由两个端点的函数值为给定值只能够产生两个方程,仅此不足以唯一确定四个参数.注意到三次样条函数对整体光滑性的要求,其二阶导数存在且连续,从全局的角度上考虑参数个数与方程个数的关系如下：
     参数：每个小段上4个,n个小段共计4n个.
     方程：每个小段上由给定函数值得到2个,n个小段共计2n个；光滑性要求每一个内部节点处的一阶、二阶导数连续,得出其左右导数相等,因此,每个节点产生2个方程,共计2(n-1)个.
     现在得到了4n-2个方程,还差两个.为此,常用的方法是对边界节点除函数值外附加要求.这就是所谓的边界条件.需要两个,正好左右两个端点各一个.常用如下三类边界条件.
    m边界条件：s'(X0)=m0,s'(Xn)=mn.即两个边界节点的一阶导数值为给定值：m0,mn.
    M边界条件：s''(x0)=m0, s''(xn)=mn.即两个边界节点的二阶导数值为给定值：m0,mn. 
    特别地,当m0和mn都为零时,称为自然边界条件. 
    周期性边界条件：s'(x0)=s'(xn); s''(x0)=s''(xn).
以上分析说明,理论上三次样条插值函数是确定的,具体如何操作,可以查阅有关资料.
    例题，观测数据为
         x=[0 1 2 3   4   5   6   7   8 9 10];
         y=[0 2 0 -4   0   4   0 -2   0 3 1];
待求的三次多项式函数s(x)在[0 10]上有连续的一阶，二阶导数.我们通过简单的讨论来认识问题。在第一区间[0 1]、第二区间[1 2]上考虑两个三次多项式
        s(x)=s1*x^3+s2*x^2+s3*x+s4
        r(x)=r1*x^3+r2*x^2+r3*x+r4
示意图: 
可以得到
          s(0)=s1*0^3+s2*0^2+s3*0+s4=0     (1)
          s(1)=s1*1^3+s2*1^2+s3*1+s4=2     (2)
          r(1)=r1*2^3+r2*2^2+r3*2+r4=2     (3)
          r(2)=r1*2^3+r2*2^2+r3*2+r4=0     (4)
一阶导函数
          s'(x)=3*s1*x^2+2*s2*x+s3
          r'(x)=3*r1*x^2+2*r2*x+r3
由一阶导数的连续性且在1点处相等，有
          3*s1*1^2+2*s2*1+s3=3*r1*1^2+2*r2*1+r3   (5)
二阶导函数
          s''(x)=6*s1*x+2*s2
          r''(x)=6*r1*x+2*r2
由二阶导数的连续性且在1点处相等，有
          6*s1*1+2*s2=6*r1*1+2*r2          (6)
由m边界条件
         s'(0)=1.6,r'(2)=0.3
有     
         3*s1*0^2+2*s2*0+s3=1.6            (7)
         3*r1*2^2+2*r2*2+r3=0.3            (8)
M边界条件
         s''(0)=-1,r''(2)=1
有
         6*s1*0+2*s2=-1                   (7')
         6*r1*2+2*r2=1                    (8')
由周期性边界条件
         s'(0)=r'(2)
         s''(0)=r''(2)
有
         3*s1*0^2+2*s2*0+s3=-1            (7'')
         3*r1*2^2+2*r2*2+r3=1             (8'')
这样，对于两个多项式的8个未知量，我们给出了8个方程。三次样条曲线的难点在于，我们不能分段去求解方程，完成绘图。
 CubicSplineInterpolation.h
view plaincopy to clipboardprint?
 // CubicSplineInterpolation.h: interface for the CCubicSplineInterpolation class.  
//  
/  
#if !defined(AFX_CUBICSPLINEINTERPOLATION_H__58C0304E_719D_4738_B86C_26BFA529B203__INCLUDED_)  
#define AFX_CUBICSPLINEINTERPOLATION_H__58C0304E_719D_4738_B86C_26BFA529B203__INCLUDED_  
 
#if _MSC_VER > 1000  
#pragma once  
#endif // _MSC_VER > 1000  
 
class CCubicSplineInterpolation    
{  
public:  
   
 //  
 // Construction/Destruction  
 //  
 //ctrlPtX 控制点X数组  
 //ctrlPtY 控制点Y数组  
 //nCtrlPtCount 控制点数目，控制点要大于等于3.  
 //  
    //  
 CCubicSplineInterpolation(double *ctrlPtX,double *ctrlPtY,int nCtrlPtCount);  
 virtual ~CCubicSplineInterpolation();  
 
public:  
   
 //  
 //outPtCount 想要输出的插值点数目,输出的点数组要大于1  
 //outPtX     已经分配好内存的X值的数组。  
 //outPtY     已经分配好内存的Y值的数组。  
 //  
 //调用此函数，获得插值点数组  
 //  
 //计算成功返回true，计算失败返回false  
    //  
 bool GetInterpolationPts(int outPtCount,double* outPtX,double* outPtY);  
 
 //  
 //根据X值 计算Y值  
 //dbInX   x自变量值，输入  
 //dbOutY  计算得到的Y值 输出  
 //  
 bool GetYByX(const double &dbInX, double &dbOutY);  
 
protected:  
 void ReleaseMem();  
 void InitParam();  
 bool InterPolation();  
 bool Spline();  
 
protected:  
 bool m_bCreate; //类是否创建成功，即控制点是否有效  
 
 int N;   //输入控制点数量  
 int M;   //输出的插入点数量  
 
 typedef double* pDouble;  
 
 pDouble X,Y; //输入的控制点数组  
 pDouble Z,F; //输出的控制点数组  
 
 pDouble H,A,B,C,D; //间距，缓存运算中间结果。  
 
};  
 
#endif // !defined(AFX_CUBICSPLINEINTERPOLATION_H__58C0304E_719D_4738_B86C_26BFA529B203__INCLUDED_) 
 // CubicSplineInterpolation.h: interface for the CCubicSplineInterpolation class.
//
/
#if !defined(AFX_CUBICSPLINEINTERPOLATION_H__58C0304E_719D_4738_B86C_26BFA529B203__INCLUDED_)
#define AFX_CUBICSPLINEINTERPOLATION_H__58C0304E_719D_4738_B86C_26BFA529B203__INCLUDED_
#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000
class CCubicSplineInterpolation 
{
public:
 
 //
 // Construction/Destruction
 //
 //ctrlPtX 控制点X数组
 //ctrlPtY 控制点Y数组
 //nCtrlPtCount 控制点数目，控制点要大于等于3.
 //
    //
 CCubicSplineInterpolation(double *ctrlPtX,double *ctrlPtY,int nCtrlPtCount);
 virtual ~CCubicSplineInterpolation();
public:
 
 //
 //outPtCount 想要输出的插值点数目,输出的点数组要大于1
 //outPtX     已经分配好内存的X值的数组。
 //outPtY     已经分配好内存的Y值的数组。
 //
 //调用此函数，获得插值点数组
 //
 //计算成功返回true，计算失败返回false
    //
 bool GetInterpolationPts(int outPtCount,double* outPtX,double* outPtY);
 //
 //根据X值 计算Y值
 //dbInX   x自变量值，输入
 //dbOutY  计算得到的Y值 输出
 //
 bool GetYByX(const double &dbInX, double &dbOutY);
protected:
 void ReleaseMem();
 void InitParam();
 bool InterPolation();
 bool Spline();
protected:
 bool m_bCreate; //类是否创建成功，即控制点是否有效
 int N;   //输入控制点数量
 int M;   //输出的插入点数量
 typedef double* pDouble;
 pDouble X,Y; //输入的控制点数组
 pDouble Z,F; //输出的控制点数组
 pDouble H,A,B,C,D; //间距，缓存运算中间结果。
};
#endif // !defined(AFX_CUBICSPLINEINTERPOLATION_H__58C0304E_719D_4738_B86C_26BFA529B203__INCLUDED_)
 
 
CubicSplineInterpolation.cpp
view plaincopy to clipboardprint?
// CubicSplineInterpolation.cpp: implementation of the CCubicSplineInterpolation class.  
//  
//  
 
#include "stdafx.h"  
#include "CubicSplineInterpolation.h"  
 
#ifdef _DEBUG  
#undef THIS_FILE  
static char THIS_FILE[]=__FILE__;  
#define new DEBUG_NEW  
#endif  
 
 
#include <math.h>  
//  
// Construction/Destruction  
//  
//ctrlPtX 控制点X数组  
//ctrlPtY 控制点Y数组  
//nCtrlPtCount 控制点数目，控制点要大于等于3.  
//  
//  
 
CCubicSplineInterpolation::CCubicSplineInterpolation(double *ctrlPtX,double *ctrlPtY,int nCtrlPtCount)  
{  
    InitParam();  
      
    if (NULL == ctrlPtX || NULL == ctrlPtY || nCtrlPtCount < 3 )  
    {  
        m_bCreate = FALSE;  
    }  
    else 
    {  
        N = nCtrlPtCount - 1;  
          
        int nDataCount = N + 1;  
        X = new double[nDataCount];  
        Y = new double[nDataCount];  
          
        A = new double[nDataCount];  
        B = new double[nDataCount];  
        C = new double[nDataCount];  
        D = new double[nDataCount];  
        H = new double[nDataCount];  
          
        memcpy(X,ctrlPtX,nDataCount*sizeof(double));  
        memcpy(Y,ctrlPtY,nDataCount*sizeof(double));  
          
        m_bCreate = Spline();         
    }     
}  
 
//  
//outPtCount 想要输出的插值点数目,输出的点数组要大于1  
//outPtX     已经分配好内存的X值的数组。  
//outPtY     已经分配好内存的Y值的数组。  
//  
//调用此函数，获得插值点数组  
//  
//计算成功返回true，计算失败返回false  
//  
bool CCubicSplineInterpolation::GetInterpolationPts(int outPtCount, double *outPtX, double *outPtY)  
{  
    if (!m_bCreate)  
    {  
        return m_bCreate;  
    }  
 
    M = outPtCount - 1;  
 
    if (M == 0)  
    {  
        return false;  
    }  
 
    Z = outPtX;  
    F = outPtY;  
 
 
 
    return InterPolation();  
 
      
}  
 
CCubicSplineInterpolation::~CCubicSplineInterpolation()  
{  
     ReleaseMem();  
}  
 
void CCubicSplineInterpolation::InitParam()  
{  
    X = Y = Z = F = A = B = C = D = H = NULL;  
 
    N = 0;  
    M = 0;  
}  
 
void CCubicSplineInterpolation::ReleaseMem()  
{  
    delete [] X;  
    delete [] Y;  
//  delete [] Z;  
//  delete [] F;  
    delete [] A;  
    delete [] B;  
    delete [] C;  
    delete [] D;  
    delete [] H;  
      
    InitParam();  
}  
 
 
bool CCubicSplineInterpolation::Spline()  
{  
    int i,P,L;  
      
    for (i=1;i<=N;i++)  
    {  
        H[i-1]=X[i]-X[i-1];  
    }  
      
    L=N-1;  
    for(i=1;i<=L;i++)  
    {  
        A[i]=H[i-1]/(H[i-1]+H[i]);  
        B[i]=3*((1-A[i])*(Y[i]-Y[i-1])/H[i-1]+A[i]*(Y[i+1]-Y[i])/H[i]);  
    }  
    A[0]=1;  
    A[N]=0;  
    B[0]=3*(Y[1]-Y[0])/H[0];  
    B[N]=3*(Y[N]-Y[N-1])/H[N-1];  
      
    for(i=0;i<=N;i++)  
    {  
        D[i]=2;  
    }  
      
    for(i=0;i<=N;i++)  
    {  
        C[i]=1-A[i];  
    }  
    P=N;  
    for(i=1;i<=P;i++)  
    {  
        if (  fabs(D[i]) <= 0.000001 )                                 
        {  
            return false;  
            //    MessageBox(0,"无解","提示,MB_OK);  
            //break;  
        }  
        A[i-1]=A[i-1]/D[i-1];  
        B[i-1]=B[i-1]/D[i-1];  
        D[i]=A[i-1]*(-C[i])+D[i];  
        B[i]=-C[i]*B[i-1]+B[i];  
    }  
    B[P]=B[P]/D[P];  
    for(i=1;i<=P;i++)  
    {  
        B[P-i]=B[P-i]-A[P-i]*B[P-i+1];  
    }  
    return true;  
}  
 
bool CCubicSplineInterpolation::InterPolation()  
{  
 
    double dbStep = (X[N] - X[0])/(M);  
      
    for (int i = 0;i <= M ;++i)  
    {  
        Z[i] = X[0] + dbStep*i;  
    }  
 
    for(i=1;i<=M;i++)  
    {  
        if (!GetYByX(Z[i],F[i]))  
        {  
            return false;  
        }  
 
    }  
 
    F[0] = Y[0];  
 
    return true;  
}  
 
 
 
bool CCubicSplineInterpolation::GetYByX(const double &dbInX, double &dbOutY)  
{  
 
    if (!m_bCreate)  
    {  
        return m_bCreate;  
    }  
 
    double E,E1,K,K1,H1;  
    int j ;   
    if(dbInX<X[0])  
    {  
        j = 0;  
          
    }  
    else if (dbInX > X[N])  
    {  
        j = N-1;  
    }  
    else 
    {  
        for (j=1;j<=N;j++)  
        {  
            if(dbInX<=X[j])  
            {  
                j=j-1;  
                  
                break;  
            }  
        }  
          
    }  
      
    //  
    E=X[j+1]-dbInX;  
    E1=E*E;  
    K=dbInX-X[j];  
    K1=K*K;  
    H1=H[j]*H[j];  
      
    dbOutY=(3*E1-2*E1*E/H[j])*Y[j]+(3*K1-2*K1*K/H[j])*Y[j+1];  
    dbOutY=dbOutY+(H[j]*E1-E1*E)*B[j]-(H[j]*K1-K1*K)*B[j+1];  
    dbOutY=dbOutY/H1;  
 
    return true;  
} 
// CubicSplineInterpolation.cpp: implementation of the CCubicSplineInterpolation class.
//
//
#include "stdafx.h"
#include "CubicSplineInterpolation.h"
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif


#include <math.h>
//
// Construction/Destruction
//
//ctrlPtX 控制点X数组
//ctrlPtY 控制点Y数组
//nCtrlPtCount 控制点数目，控制点要大于等于3.
//
//
CCubicSplineInterpolation::CCubicSplineInterpolation(double *ctrlPtX,double *ctrlPtY,int nCtrlPtCount)
{
 InitParam();
 
 if (NULL == ctrlPtX || NULL == ctrlPtY || nCtrlPtCount < 3 )
 {
  m_bCreate = FALSE;
 }
 else
 {
  N = nCtrlPtCount - 1;
  
  int nDataCount = N + 1;
  X = new double[nDataCount];
  Y = new double[nDataCount];
  
  A = new double[nDataCount];
  B = new double[nDataCount];
  C = new double[nDataCount];
  D = new double[nDataCount];
  H = new double[nDataCount];
  
  memcpy(X,ctrlPtX,nDataCount*sizeof(double));
  memcpy(Y,ctrlPtY,nDataCount*sizeof(double));
  
  m_bCreate = Spline();  
 }  
}
//
//outPtCount 想要输出的插值点数目,输出的点数组要大于1
//outPtX     已经分配好内存的X值的数组。
//outPtY     已经分配好内存的Y值的数组。
//
//调用此函数，获得插值点数组
//
//计算成功返回true，计算失败返回false
//
bool CCubicSplineInterpolation::GetInterpolationPts(int outPtCount, double *outPtX, double *outPtY)
{
 if (!m_bCreate)
 {
  return m_bCreate;
 }
 M = outPtCount - 1;
 if (M == 0)
 {
  return false;
 }
 Z = outPtX;
 F = outPtY;
 
 return InterPolation();
 
}
CCubicSplineInterpolation::~CCubicSplineInterpolation()
{
     ReleaseMem();
}
void CCubicSplineInterpolation::InitParam()
{
 X = Y = Z = F = A = B = C = D = H = NULL;
 N = 0;
 M = 0;
}
void CCubicSplineInterpolation::ReleaseMem()
{
 delete [] X;
 delete [] Y;
//  delete [] Z;
//  delete [] F;
 delete [] A;
 delete [] B;
 delete [] C;
 delete [] D;
 delete [] H;
   
 InitParam();
}


bool CCubicSplineInterpolation::Spline()
{
 int i,P,L;
 
 for (i=1;i<=N;i++)
 {
  H[i-1]=X[i]-X[i-1];
 }
 
 L=N-1;
 for(i=1;i<=L;i++)
 {
  A[i]=H[i-1]/(H[i-1]+H[i]);
  B[i]=3*((1-A[i])*(Y[i]-Y[i-1])/H[i-1]+A[i]*(Y[i+1]-Y[i])/H[i]);
 }
 A[0]=1;
 A[N]=0;
 B[0]=3*(Y[1]-Y[0])/H[0];
 B[N]=3*(Y[N]-Y[N-1])/H[N-1];
 
 for(i=0;i<=N;i++)
 {
  D[i]=2;
 }
 
 for(i=0;i<=N;i++)
 {
  C[i]=1-A[i];
 }
 
 P=N;
 for(i=1;i<=P;i++)
 {
  if (  fabs(D[i]) <= 0.000001 )                              
  {
   return false;
   //    MessageBox(0,"无解","提示,MB_OK);
   //break;
  }
  A[i-1]=A[i-1]/D[i-1];
  B[i-1]=B[i-1]/D[i-1];
  D[i]=A[i-1]*(-C[i])+D[i];
  B[i]=-C[i]*B[i-1]+B[i];
 }
 B[P]=B[P]/D[P];
 for(i=1;i<=P;i++)
 {
  B[P-i]=B[P-i]-A[P-i]*B[P-i+1];
 }
 return true;
}
bool CCubicSplineInterpolation::InterPolation()
{
 double dbStep = (X[N] - X[0])/(M);
 
 for (int i = 0;i <= M ;++i)
 {
  Z[i] = X[0] + dbStep*i;
 }
 for(i=1;i<=M;i++)
 {
        if (!GetYByX(Z[i],F[i]))
        {
            return false;
        }
 }
 F[0] = Y[0];
 return true;
}
 
bool CCubicSplineInterpolation::GetYByX(const double &dbInX, double &dbOutY)
{
 if (!m_bCreate)
 {
  return m_bCreate;
 }
 double E,E1,K,K1,H1;
 int j ; 
 if(dbInX<X[0])
 {
  j = 0;
  
 }
 else if (dbInX > X[N])
 {
  j = N-1;
 }
 else
 {
  for (j=1;j<=N;j++)
  {
   if(dbInX<=X[j])
   {
    j=j-1;
    
    break;
   }
  }
  
 }
 
 //
 E=X[j+1]-dbInX;
 E1=E*E;
 K=dbInX-X[j];
 K1=K*K;
 H1=H[j]*H[j];
 
 dbOutY=(3*E1-2*E1*E/H[j])*Y[j]+(3*K1-2*K1*K/H[j])*Y[j+1];
 dbOutY=dbOutY+(H[j]*E1-E1*E)*B[j]-(H[j]*K1-K1*K)*B[j+1];
 dbOutY=dbOutY/H1;
 return true;

}

原文链接http://blog.sina.com.cn/s/blog_9de219560100yxqr.html