最小二乘法 直线拟合

/
  功能描述：利用最小二乘法求斜率                                   //
            xyTopX -- 截距                                         //
//方法1：

            xyTopY -- 斜率                                         //
            nCount -- 点数                                         //
/
void  Min2Method(double &xyTopY, double &xyTopX, int  X[], int Y[], int nCount)
{
    int      i;
    double   SumX, SumY, SumXY, SumX2;
   
    SumX = 0;
    SumX2 = 0;
    for( i=0; i<nCount; i++)
    {
        SumX += X[i];
        SumX2 += (X[i]*X[i]);
    }
   
    SumY = 0;
    for( i=0; i<nCount; i++)
    {
        SumY += Y[i];
    }
 
    SumXY = 0;
    for( i=0; i<nCount; i++)
    {
        SumXY += ( X[i]* Y[i]);
    }
 
    xyTopX = ( (SumX2*SumY - SumX*SumXY) / (nCount*SumX2 - SumX * SumX));
    xyTopY = ( (nCount*SumXY - SumX*SumY) / (nCount*SumX2 - SumX * SumX));
}

void CTmfcDlg::OnButton2()
{
 
 int x [100]; int y[100];
 
 
 Graphics graphics(this->m_hWnd); Pen      pen(Color(255, 0, 0, 255));
 
 
 
 for( int i = 0 ; i<100; i++)
 { 
  x[i] = i ;  y[i] = i ;
 
  
  if ( i%3 == 0 )
  {
   x[i] = i +10 ;  y[i] = i ;   
  }
 
  
  graphics.DrawLine(&pen, Point(x[i],y[i]),Point(x[i],y[i]+1)) ;
 }
 

 double k, b;

 Min2Method(k,b, x, y, 100 );         showmessage( k, b );  
 
 
 
 graphics.DrawLine(&pen, 0, b, 100, k*100+b );      // Y=kX+b;
 
}

 

//方法2： (这是网上最多的直线拟合函数)

 

class CFoldPoint 
{
public: 
 CFoldPoint( double x, double y)
 {
  X = x;
  Y = y;
 }
 
 double X;
 double Y; 
};

typedef CTypedPtrList<CPtrList,CFoldPoint*> CFoldPointList;
 
CFoldPointList  FoldList;

BOOL CalculateLineKB(CFoldPointList *m_FoldList,double &k,double &b)
{
 //最小二乘法直线拟合
 //m_FoldList为关键点(x,y)的链表
 //拟合直线方程(Y=kX+b)
 
 if(m_FoldList==NULL)  return FALSE;
 long lCount=m_FoldList->GetCount();
 if(lCount<2)return FALSE;
 CFoldPoint *pFold;
 double mX,mY,mXX,mXY,n;
 mX=mY=mXX=mXY=0;
 n=lCount;
 POSITION pos=m_FoldList->GetHeadPosition();
 while(pos != NULL)
 {
  pFold=m_FoldList->GetNext(pos);
  mX+=pFold->X;
  mY+=pFold->Y;
  mXX+=pFold->X*pFold->X;
  mXY+=pFold->X*pFold->Y;
 }
 if(mX*mX-mXX*n==0)return FALSE;
 k=(mY*mX-mXY*n)/(mX*mX-mXX*n);
 b=(mY-mX*k)/n;
 return TRUE;
}
 

void CTmfcDlg::OnButton3()
{
 // TODO: Add your control notification handler code here
 

 int x [100]; int y[100];

 
 Graphics graphics(this->m_hWnd); Pen      pen(Color(255, 0, 0, 255));
  
 

 for( int i = 0 ; i<100; i++)
 {
  if (i%2 == 0)
  { continue;
   x[i] = i  +22 ;
   y[i] = i;
   
 
   FoldList.AddTail( &CFoldPoint(x[i],y[i]) );

  }
   else
  {
   x[i] = i -22;
   y[i] = i; 
 
   FoldList.AddTail( &CFoldPoint(x[i],y[i]) );

 
  }
 }

 
 double  k,  b ;
 CalculateLineKB( &FoldList, k, b);

 showmessage( k, b );

 
 graphics.DrawLine(&pen, 0, b, 100, k*100+b );
 
}