开心六月综合激情婷婷|欧美精品成人动漫二区|国产中文字幕综合色|亚洲人在线成视频

最小二乘法是一種優(yōu)化算法，最小二乘法名字的緣由有兩個：一是要將誤差最小化，二是將誤差最小化的方法是使誤差的平方和最小化。利用最小二乘法可以簡便地求得未知的數(shù)據(jù)，并使得這些求得的數(shù)據(jù)與實際數(shù)據(jù)之間誤差的平方和為最小。最小二乘法還可用于曲線擬合，所擬合的曲線可以是線性擬合與非線性擬合。

--------------------一元線性函數(shù)--------------------

形式1：借用案例（http://blog.csdn.net/qll125596718/article/details/8248249）首先以一元線性方程參數(shù)估計為例，樣本回歸模型：

殘差平方和：

則通過Q最小確定這條直線，即確定

，以為變量，把它們看作是Q的函數(shù)，就變成了一個求極值的問題，可以通過求導數(shù)得到。求Q對兩個待估參數(shù)的偏導數(shù)：

解得：

#include <iostream>    
#include <algorithm>    
#include <valarray>    
#include <vector>    
using namespace std;
int main()
{
	double x[] = { 1, 2, 3, 4, 5, 6 };
	double y[] = { 3, 5.5, 6.8, 8.8, 11, 12};
	valarray<double> data_x(x, 6);
	valarray<double> data_y(y, 6);
	float A = 0.0;
	float B = 0.0;
	float C = 0.0;
	float D = 0.0;
	A = (data_x*data_x).sum();
	B = data_x.sum();
	C = (data_x*data_y).sum();
	D = data_y.sum();
	float tmp = A*data_x.size() - B*B;
	float k, b;
	k = (C*data_x.size() - B*D) / tmp;
	b = (A*D - C*B) / tmp;
	cout << "y=" << k << "x+" << b << endl;
	return 0;
}

運行結果：

注：valarray類似vector，也是一個模板類，主要用來對一系列元素進行高速的數(shù)字計算，其與vector的主要區(qū)別在于以下兩點：
1、valarray定義了一組在兩個相同長度和相同類型的valarray類對象之間的數(shù)字計算；
2、通過重載operater[]，可以返回valarray的相關信息（valarray其中某個元素的引用、特定下標的值或者其某個子集）。

valarray類構造函數(shù):
valarray( );
explicit valarray(size_t _Count);
valarray( const Type& _Val, size_t _Count);
valarray( const Type* _Ptr, size_t _Count);
valarray( const slice_array<Type>& _SliceArray);
valarray( const gslice_array<Type>& _GsliceArray);
valarray( const mask_array<Type>& _MaskArray);
valarray( const indirect_array<Type>& _IndArray);

valarray 類用法:
1. apply 將 valarray 數(shù)組的每一個值都用 apply 所接受到的函數(shù)進行計算
2. cshift 將 valarray 數(shù)組的數(shù)據(jù)進行循環(huán)移動，參數(shù)為正者左移為負就右移
3. max 返回 valarray 數(shù)組的最大值
4. min 返回 valarray 數(shù)組的最小值
5. resize 重新設置 valarray 數(shù)組大小，并對其進行初始化
6. shift 將 valarray 數(shù)組移動，參數(shù)為正者左移，為負者右移，移動后由 0 填充剩余位
7. size 得到數(shù)組的大小
8. sum 數(shù)組求和

--------------------N元線性函數(shù)--------------------

一元線性方程可以看做多元函數(shù)的特例，現(xiàn)在用矩陣形式表述多元函數(shù)情況下，最小二乘的一般形式：

 設擬合多項式為：

                   

各點到這條曲線的距離之和，即偏差平方和如下：

                

對等式右邊求ai偏導數(shù)，得到： 

              

             

                            .......

             

把這些等式表示成矩陣的形式，就可以得到下面的矩陣：

              

  （3）

進行化簡計算：

           

， 

上面公式（3）可以寫為:

                   

，

              

#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "vector"
using namespace std;
struct point
{
	double x;
	double y;
};
typedef vector<double> doubleVector;
vector<point> getFile(char *File);  //獲取文件數(shù)據(jù)
doubleVector getCoeff(vector<point> sample, int n);   //矩陣方程
void main()
{
	int i, n;
	char *File = "XY.txt";
	vector<point> sample;
	doubleVector  coefficient;
	sample = getFile(File);
	printf("擬合多項式階數(shù)n=");
  	  scanf_s("%d", &n);
	  coefficient = getCoeff(sample, n);
	printf("\n擬合矩陣的系數(shù)為：\n");
	for (i = 0; i < coefficient.size(); i++)
		printf("a%d = %lf\n", i, coefficient[i]);
}
//矩陣方程
doubleVector getCoeff(vector<point> sample, int n)
{
	vector<doubleVector> matFunX;  //公式3左矩陣
	vector<doubleVector> matFunY;  //公式3右矩陣
	doubleVector temp;
	double sum;
	int i, j, k;
	//公式3左矩陣
	for (i = 0; i <= n; i++)
	{
		temp.clear();
		for (j = 0; j <= n; j++)
		{
			sum = 0;
			for (k = 0; k < sample.size(); k++)
				sum += pow(sample[k].x, j + i);
			temp.push_back(sum);
		}
		matFunX.push_back(temp);
	}
	//printf("matFunX.size=%d\n", matFunX.size());
	//printf("matFunX[3][3]=%f\n", matFunX[3][3]);
	//公式3右矩陣
	for (i = 0; i <= n; i++)
	{
		temp.clear();
		sum = 0;
		for (k = 0; k < sample.size(); k++)
			sum += sample[k].y*pow(sample[k].x, i);
		temp.push_back(sum);
		matFunY.push_back(temp);
	}
	printf("matFunY.size=%d\n", matFunY.size());
	//矩陣行列式變換
	double num1, num2, ratio;
	for (i = 0; i < matFunX.size() - 1; i++)
	{
		num1 = matFunX[i][i];
		for (j = i + 1; j < matFunX.size(); j++)
		{
			num2 = matFunX[j][i];
			ratio = num2 / num1;
			for (k = 0; k < matFunX.size(); k++)
				matFunX[j][k] = matFunX[j][k] - matFunX[i][k] * ratio;
			matFunY[j][0] = matFunY[j][0] - matFunY[i][0] * ratio;
		}
	}
	//計算擬合曲線的系數(shù)
	doubleVector coeff(matFunX.size(), 0);
	for (i = matFunX.size() - 1; i >= 0; i--)
	{
		if (i == matFunX.size() - 1)
			coeff[i] = matFunY[i][0] / matFunX[i][i];
		else
		{
			for (j = i + 1; j < matFunX.size(); j++)
				matFunY[i][0] = matFunY[i][0] - coeff[j] * matFunX[i][j];
			coeff[i] = matFunY[i][0] / matFunX[i][i];
		}
	}
	return coeff;
}
//獲取文件數(shù)據(jù)
vector<point> getFile(char *File)
{
	int i = 1;
	vector<point> dst;
	FILE *fp=fopen(File, "r");
	if (fp == NULL)
	{
		printf("Open file error!!!\n");
		exit(0);
	}
	point temp;
	double num;
	while (fscanf(fp, "%lf", &num) != EOF)
	{
		if (i % 2 == 0)
		{
			temp.y = num;
			dst.push_back(temp);
		}
		else
			temp.x = num;
		i++;
	}
	fclose(fp);
	return dst;
}

XY.txt內容：

0    1.0
0.25 1.28
0.5  1.65
0.75 2.12
1    2.72

另外在http://blog.csdn.net/lsh_2013/article/details/46697625里也有相關程序：

#include <iostream>  
#include <vector>  
#include <cmath>  
using namespace std;  
//最小二乘擬合相關函數(shù)定義  
double sum(vector<double> Vnum, int n);  
double MutilSum(vector<double> Vx, vector<double> Vy, int n);  
double RelatePow(vector<double> Vx, int n, int ex);  
double RelateMutiXY(vector<double> Vx, vector<double> Vy, int n, int ex);  
void EMatrix(vector<double> Vx, vector<double> Vy, int n, int ex, double coefficient[]);  
void CalEquation(int exp, double coefficient[]);  
double F(double c[],int l,int m);  
double Em[6][4];  
//主函數(shù)，這里將數(shù)據(jù)擬合成二次曲線  
int main(int argc, char* argv[])  
{  
    double arry1[5]={0,0.25,0,5,0.75};  
    double arry2[5]={1,1.283,1.649,2.212,2.178};  
    double coefficient[5];  
    memset(coefficient,0,sizeof(double)*5);  
    vector<double> vx,vy;  
    for (int i=0; i<5; i++)  
    {  
        vx.push_back(arry1[i]);  
        vy.push_back(arry2[i]);  
    }  
    EMatrix(vx,vy,5,3,coefficient);  
    printf("擬合方程為：y = %lf + %lfx + %lfx^2 \n",coefficient[1],coefficient[2],coefficient[3]);  
    return 0;  
}  
//累加  
double sum(vector<double> Vnum, int n)  
{  
    double dsum=0;  
    for (int i=0; i<n; i++)  
    {  
        dsum+=Vnum[i];  
    }  
    return dsum;  
}  
//乘積和  
double MutilSum(vector<double> Vx, vector<double> Vy, int n)  
{  
    double dMultiSum=0;  
    for (int i=0; i<n; i++)  
    {  
        dMultiSum+=Vx[i]*Vy[i];  
    }  
    return dMultiSum;  
}  
//ex次方和  
double RelatePow(vector<double> Vx, int n, int ex)  
{  
    double ReSum=0;  
    for (int i=0; i<n; i++)  
    {  
        ReSum+=pow(Vx[i],ex);  
    }  
    return ReSum;  
}  
//x的ex次方與y的乘積的累加  
double RelateMutiXY(vector<double> Vx, vector<double> Vy, int n, int ex)  
{  
    double dReMultiSum=0;  
    for (int i=0; i<n; i++)  
    {  
        dReMultiSum+=pow(Vx[i],ex)*Vy[i];  
    }  
    return dReMultiSum;  
}  
//計算方程組的增廣矩陣  
void EMatrix(vector<double> Vx, vector<double> Vy, int n, int ex, double coefficient[])  
{  
    for (int i=1; i<=ex; i++)  
    {  
        for (int j=1; j<=ex; j++)  
        {  
            Em[i][j]=RelatePow(Vx,n,i+j-2);  
        }  
        Em[i][ex+1]=RelateMutiXY(Vx,Vy,n,i-1);  
    }  
    Em[1][1]=n;  
    CalEquation(ex,coefficient);  
}  
//求解方程  
void CalEquation(int exp, double coefficient[])  
{  
    for(int k=1;k<exp;k++) //消元過程  
    {  
        for(int i=k+1;i<exp+1;i++)  
        {  
            double p1=0;  
            if(Em[k][k]!=0)  
                p1=Em[i][k]/Em[k][k];  
            for(int j=k;j<exp+2;j++)   
                Em[i][j]=Em[i][j]-Em[k][j]*p1;  
        }  
    }  
    coefficient[exp]=Em[exp][exp+1]/Em[exp][exp];  
    for(int l=exp-1;l>=1;l--)   //回代求解  
        coefficient[l]=(Em[l][exp+1]-F(coefficient,l+1,exp))/Em[l][l];  
}  
//供CalEquation函數(shù)調用  
double F(double c[],int l,int m)  
{  
    double sum=0;  
    for(int i=l;i<=m;i++)  
        sum+=Em[l-1][i]*c[i];  
    return sum;   
}  

memset相關介紹：

http://baike.baidu.com/link?url=p7JreiRCj9yPs3r3WAfsXgynjvtGrWoQ_exF9tFGK6fsVP7V6tdm-_13QhCZxqPrfRi0wH0EihhRL_-qVvrewq

http://c.biancheng.net/cpp/html/157.html

--------------------擬合圓的方程--------------------

/*
最小二乘法擬合圓，擬合出的圓以圓心坐標和半徑的形式表示
 */
typedef complex<int> POINT;
bool FitCircle(const std::vector<POINT> &points, double &er_x, double &er_y, double &radius)
{
     cent_x = 0.0f;
     cent_y = 0.0f;
     radius = 0.0f;
     if (points.size() < 3)
     {
         return false;
     }
     double sum_x = 0.0f, sum_y = 0.0f;
     double sum_x2 = 0.0f, sum_y2 = 0.0f;
     double sum_x3 = 0.0f, sum_y3 = 0.0f;
     double sum_xy = 0.0f, sum_x1y2 = 0.0f, sum_x2y1 = 0.0f;
     int N = points.size();
     for (int i = 0; i < N; i++)
     {
         double x = points[i].real();
         double y = points[i].imag();
         double x2 = x * x;
         double y2 = y * y;
         sum_x += x;
         sum_y += y;
         sum_x2 += x2;
         sum_y2 += y2;
         sum_x3 += x2 * x;
         sum_y3 += y2 * y;
         sum_xy += x * y;
         sum_x1y2 += x * y2;
         sum_x2y1 += x2 * y;
     }
     double C, D, E, G, H;
     double a, b, c;
     C = N * sum_x2 - sum_x * sum_x;
     D = N * sum_xy - sum_x * sum_y;
     E = N * sum_x3 + N * sum_x1y2 - (sum_x2 + sum_y2) * sum_x;
     G = N * sum_y2 - sum_y * sum_y;
     H = N * sum_x2y1 + N * sum_y3 - (sum_x2 + sum_y2) * sum_y;
     a = (H * D - E * G) / (C * G - D * D);
     b = (H * C - E * D) / (D * D - G * C);
     c = -(a * sum_x + b * sum_y + sum_x2 + sum_y2) / N;
     cent_x = a / (-2);
     cent_y = b / (-2);
     radius = sqrt(a * a + b * b - 4 * c) / 2;
     return true;
}

--------------------擬合橢圓方程--------------------

//LSEllipse.h

/*************************************************************************
	功能說明： 對平面上的一些列點給出最小二乘的橢圓擬合，利用奇異值分解法
	           解得最小二乘解作為橢圓參數(shù)。
	調用形式： cvFitEllipse2f（arrayx,arrayy,box）；    
	參數(shù)說明： arrayx: arrayx[n],每個值為x軸一個點
	           arrayx: arrayy[n],每個值為y軸一個點
			   n     : 點的個數(shù)
			   box   : box[5],橢圓的五個參數(shù)，分別為center.x,center.y,2a,2b,xtheta
			   esp: 解精度，通常取1e-6,這個是解方程用的說
***************************************************************************/
#include<cstdlib>
#include<float.h>
#include<vector>
using namespace std;
class LSEllipse
{
public:
	LSEllipse(void);
	~LSEllipse(void);
	vector<double> getEllipseparGauss(vector<CPoint> vec_point);
	void cvFitEllipse2f( int *arrayx, int *arrayy,int n,float *box );
private:
	int SVD(float *a,int m,int n,float b[],float x[],float esp);
	int gmiv(float a[],int m,int n,float b[],float x[],float aa[],float eps,float u[],float v[],int ka);
    int ginv(float a[],int m,int n,float aa[],float eps,float u[],float v[],int ka);
	int muav(float a[],int m,int n,float u[],float v[],float eps,int ka);
};

//LSEllipse.cpp

#include "LSEllipse.h"
#include <cmath>
LSEllipse::LSEllipse(void)
{
}
LSEllipse::~LSEllipse(void)
{
}
//列主元高斯消去法
//A為系數(shù)矩陣，x為解向量，若成功，返回true，否則返回false，并將x清空。
bool RGauss(const vector<vector<double> > & A, vector<double> & x)
{
    x.clear();
    int n = A.size();
    int m = A[0].size();
    x.resize(n);
    //復制系數(shù)矩陣，防止修改原矩陣
    vector<vector<double> > Atemp(n);
    for (int i = 0; i < n; i++)
    {
        vector<double> temp(m);
        for (int j = 0; j < m; j++)
        {
            temp[j] = A[i][j];
        }
        Atemp[i] = temp;
        temp.clear();
    }
    for (int k = 0; k < n; k++)
    {
        //選主元
        double max = -1;
        int l = -1;
        for (int i = k; i < n; i++)
        {
            if (abs(Atemp[i][k]) > max)
            {
                max = abs(Atemp[i][k]);
                l = i;
            }
        }
        if (l != k)
        {
            //交換系數(shù)矩陣的l行和k行
            for (int i = 0; i < m; i++)
            {
                double temp = Atemp[l][i];
                Atemp[l][i] = Atemp[k][i];
                Atemp[k][i] = temp;
            }
        }
        //消元
        for (int i = k+1; i < n; i++)
        {
            double l = Atemp[i][k]/Atemp[k][k];
            for (int j = k; j < m; j++)
            {
                Atemp[i][j] = Atemp[i][j] - l*Atemp[k][j];
            }
        }
    }
    //回代
    x[n-1] = Atemp[n-1][m-1]/Atemp[n-1][m-2];
    for (int k = n-2; k >= 0; k--)
    {
        double s = 0.0;
        for (int j = k+1; j < n; j++)
        {
            s += Atemp[k][j]*x[j];
        }
        x[k] = (Atemp[k][m-1] - s)/Atemp[k][k];
    }
    return true;
}
vector<double>  LSEllipse::getEllipseparGauss(vector<CPoint> vec_point)
{
	vector<double> vec_result;
	double x3y1 = 0,x1y3= 0,x2y2= 0,yyy4= 0, xxx3= 0,xxx2= 0,x2y1= 0,yyy3= 0,x1y2= 0 ,yyy2= 0,x1y1= 0,xxx1= 0,yyy1= 0;
	int N = vec_point.size();
	for (int m_i = 0;m_i < N ;++m_i )
	{
		double xi = vec_point[m_i].x ;
		double yi = vec_point[m_i].y;
		x3y1 += xi*xi*xi*yi ;
		x1y3 += xi*yi*yi*yi;
		x2y2 += xi*xi*yi*yi; ;
		yyy4 +=yi*yi*yi*yi;
		xxx3 += xi*xi*xi ;
		xxx2 += xi*xi ;
		x2y1 += xi*xi*yi;
		x1y2 += xi*yi*yi;
		yyy2 += yi*yi;
		x1y1 += xi*yi;
		xxx1 += xi;
		yyy1 += yi;
		yyy3 += yi*yi*yi;
	}
	double resul[5];
	resul[0] = -(x3y1);
	resul[1] = -(x2y2);
	resul[2] = -(xxx3);
	resul[3] = -(x2y1);
	resul[4] = -(xxx2);
	long double Bb[5],Cc[5],Dd[5],Ee[5],Aa[5];
	Bb[0] = x1y3, Cc[0] = x2y1, Dd[0] = x1y2, Ee[0] = x1y1, Aa[0] = x2y2;
	Bb[1] = yyy4, Cc[1] = x1y2, Dd[1] = yyy3, Ee[1] = yyy2, Aa[1] = x1y3;
	Bb[2] = x1y2, Cc[2] = xxx2, Dd[2] = x1y1, Ee[2] = xxx1, Aa[2] = x2y1;
	Bb[3] = yyy3, Cc[3]= x1y1, Dd[3] = yyy2, Ee[3] = yyy1, Aa[3] = x1y2;
	Bb[4]= yyy2, Cc[4]= xxx1, Dd[4] = yyy1, Ee[4] = N, Aa[4]= x1y1;
    vector<vector<double>>Ma(5);
	vector<double>Md(5);
	for(int i=0;i<5;i++)
	{
		Ma[i].push_back(Aa[i]);
		Ma[i].push_back(Bb[i]);
        Ma[i].push_back(Cc[i]);
        Ma[i].push_back(Dd[i]);
		Ma[i].push_back(Ee[i]);
		Ma[i].push_back(resul[i]);
	}
	RGauss(Ma,Md);
	long double A=Md[0];
	long double B=Md[1];
	long double C=Md[2];
	long double D=Md[3];
	long double E=Md[4];
	double XC=(2*B*C-A*D)/(A*A-4*B);
	double YC=(2*D-A*C)/(A*A-4*B);
	long double fenzi=2*(A*C*D-B*C*C-D*D+4*E*B-A*A*E);
	long double fenmu=(A*A-4*B)*(B-sqrt(A*A+(1-B)*(1-B))+1);
	long double fenmu2=(A*A-4*B)*(B+sqrt(A*A+(1-B)*(1-B))+1);
	double XA=sqrt(fabs(fenzi/fenmu));
	double XB=sqrt(fabs(fenzi/fenmu2));
	double Xtheta=0.5*atan(A/(1-B))*180/3.1415926;
	if(B<1)
		Xtheta+=90;
	vec_result.push_back(XC);
	vec_result.push_back(YC);
	vec_result.push_back(XA);
	vec_result.push_back(XB);
	vec_result.push_back(Xtheta);
	return vec_result;
}
void  LSEllipse::cvFitEllipse2f(  int *arrayx, int *arrayy,int n,float *box )
{   
	float cx=0,cy=0;
    double rp[5], t;
	float *A1=new float[n*5];
	float *A2=new float[2*2];
	float *A3=new float[n*3];
	float *B1=new float[n],*B2=new float[2],*B3=new float[n];
    const double min_eps = 1e-6;
    int i;
    for( i = 0; i < n; i++ )
    {
        cx += arrayx[i]*1.0;
        cy += arrayy[i]*1.0;
    }
    cx /= n;
    cy /= n;
    for( i = 0; i < n; i++ )
    {
		int step=i*5;
		float px,py;
        px = arrayx[i]*1.0;
        py = arrayy[i]*1.0;
        px -= cx;
        py -= cy;
        B1[i] = 10000.0;
        A1[step] = -px * px;
        A1[step + 1] = -py * py;
        A1[step + 2] = -px * py;
        A1[step + 3] = px;
        A1[step + 4] = py;
    }
	float *x1=new float[5];
	//解出Ax^2+By^2+Cxy+Dx+Ey=10000的最小二乘解！
	SVD(A1,n,5,B1,x1,min_eps);
	A2[0]=2*x1[0],A2[1]=A2[2]=x1[2],A2[3]=2*x1[1];
	B2[0]=x1[3],B2[1]=x1[4];
	float *x2=new float[2];
	//標準化，將一次項消掉，求出center.x和center.y;
	SVD(A2,2,2,B2,x2,min_eps);
    rp[0]=x2[0],rp[1]=x2[1];
    for( i = 0; i < n; i++ )
    {
        float px,py;
        px = arrayx[i]*1.0;
        py = arrayy[i]*1.0;
        px -= cx;
        py -= cy;
        B3[i] = 1.0;
		int step=i*3;
        A3[step] = (px - rp[0]) * (px - rp[0]);
        A3[step+1] = (py - rp[1]) * (py - rp[1]);
        A3[step+2] = (px - rp[0]) * (py - rp[1]);
    }
	//求出A(x-center.x)^2+B(y-center.y)^2+C(x-center.x)(y-center.y)的最小二乘解
	SVD(A3,n,3,B3,x1,min_eps);
    rp[4] = -0.5 * atan2(x1[2], x1[1] - x1[0]);
    t = sin(-2.0 * rp[4]);
    if( fabs(t) > fabs(x1[2])*min_eps )
        t = x1[2]/t;
    else
        t = x1[1] - x1[0];
    rp[2] = fabs(x1[0] + x1[1] - t);
    if( rp[2] > min_eps )
        rp[2] = sqrt(2.0 / rp[2]);
    rp[3] = fabs(x1[0] + x1[1] + t);
    if( rp[3] > min_eps )
        rp[3] = sqrt(2.0 / rp[3]);
    box[0] = (float)rp[0] + cx;
    box[1]= (float)rp[1] + cy;
    box[2]= (float)(rp[2]*2);
    box[3] = (float)(rp[3]*2);
    if( box[2] > box[3] )
    {
		double tmp=box[2];
		box[2]=box[3];
		box[3]=tmp;
    }
	box[4] = (float)(90 + rp[4]*180/3.1415926);
    if( box[4] < -180 )
        box[4] += 360;
    if( box[4] > 360 )
        box[4] -= 360;
	delete []A1;
	delete []A2;
	delete []A3;
    delete []B1;
	delete []B2;
	delete []B3;
	delete []x1;
	delete []x2;
}
int LSEllipse::SVD(float *a,int m,int n,float b[],float x[],float esp)
{  
	float *aa;
	float *u;
	float *v;
    aa=new float[n*m];
	u=new  float[m*m];
    v=new  float[n*n];
   int ka;
   int  flag;
   if(m>n)
   { 
	ka=m+1;
   }else
   {
	   ka=n+1;
   }
   flag=gmiv(a,m,n,b,x,aa,esp,u,v,ka);
	delete []aa;
	delete []u;
	delete []v;
	return(flag);
}
int LSEllipse::gmiv( float a[],int m,int n,float b[],float x[],float aa[],float eps,float u[],float v[],int ka)  
{ 
	int i,j;
    i=ginv(a,m,n,aa,eps,u,v,ka);
    if (i<0) return(-1);
    for (i=0; i<=n-1; i++)
      { x[i]=0.0;
        for (j=0; j<=m-1; j++)
          x[i]=x[i]+aa[i*m+j]*b[j];
      }
    return(1);
  }
int LSEllipse::ginv(float a[],int m,int n,float aa[],float eps,float u[],float v[],int ka)
  { 
 //  int muav(float a[],int m,int n,float u[],float v[],float eps,int ka);
	int i,j,k,l,t,p,q,f;
    i=muav(a,m,n,u,v,eps,ka);
    if (i<0) return(-1);
    j=n;
    if (m<n) j=m;
    j=j-1;
    k=0;
    while ((k<=j)&&(a[k*n+k]!=0.0)) k=k+1;
    k=k-1;
    for (i=0; i<=n-1; i++)
    for (j=0; j<=m-1; j++)
      { t=i*m+j; aa[t]=0.0;
        for (l=0; l<=k; l++)
          { f=l*n+i; p=j*m+l; q=l*n+l;
            aa[t]=aa[t]+v[f]*u[p]/a[q];
          }
      }
    return(1);
  }
int LSEllipse::muav(float a[],int m,int n,float u[],float v[],float eps,int ka)
  { int i,j,k,l,it,ll,kk,ix,iy,mm,nn,iz,m1,ks;
    float d,dd,t,sm,sm1,em1,sk,ek,b,c,shh,fg[2],cs[2];
    float *s,*e,*w;
    //void ppp();
   // void sss();
     void ppp(float a[],float e[],float s[],float v[],int m,int n);
     void sss(float fg[],float cs[]);
	s=(float *) malloc(ka*sizeof(float));
    e=(float *) malloc(ka*sizeof(float));
    w=(float *) malloc(ka*sizeof(float));
    it=60; k=n;
    if (m-1<n) k=m-1;
    l=m;
    if (n-2<m) l=n-2;
    if (l<0) l=0;
    ll=k;
    if (l>k) ll=l;
    if (ll>=1)
      { for (kk=1; kk<=ll; kk++)
          { if (kk<=k)
              { d=0.0;
                for (i=kk; i<=m; i++)
                  { ix=(i-1)*n+kk-1; d=d+a[ix]*a[ix];}
                s[kk-1]=(float)sqrt(d);
                if (s[kk-1]!=0.0)
                  { ix=(kk-1)*n+kk-1;
                    if (a[ix]!=0.0)
                      { s[kk-1]=(float)fabs(s[kk-1]);
                        if (a[ix]<0.0) s[kk-1]=-s[kk-1];
                      }
                    for (i=kk; i<=m; i++)
                      { iy=(i-1)*n+kk-1;
                        a[iy]=a[iy]/s[kk-1];
                      }
                    a[ix]=1.0f+a[ix];
                  }
                s[kk-1]=-s[kk-1];
              }
            if (n>=kk+1)
              { for (j=kk+1; j<=n; j++)
                  { if ((kk<=k)&&(s[kk-1]!=0.0))
                      { d=0.0;
                        for (i=kk; i<=m; i++)
                          { ix=(i-1)*n+kk-1;
                            iy=(i-1)*n+j-1;
                            d=d+a[ix]*a[iy];
                          }
                        d=-d/a[(kk-1)*n+kk-1];
                        for (i=kk; i<=m; i++)
                          { ix=(i-1)*n+j-1;
                            iy=(i-1)*n+kk-1;
                            a[ix]=a[ix]+d*a[iy];
                          }
                      }
                    e[j-1]=a[(kk-1)*n+j-1];
                  }
              }
            if (kk<=k)
              { for (i=kk; i<=m; i++)
                  { ix=(i-1)*m+kk-1; iy=(i-1)*n+kk-1;
                    u[ix]=a[iy];
                  }
              }
            if (kk<=l)
              { d=0.0;
                for (i=kk+1; i<=n; i++)
                  d=d+e[i-1]*e[i-1];
                e[kk-1]=(float)sqrt(d);
                if (e[kk-1]!=0.0)
                  { if (e[kk]!=0.0)
                      { e[kk-1]=(float)fabs(e[kk-1]);
                        if (e[kk]<0.0) e[kk-1]=-e[kk-1];
                      }
                    for (i=kk+1; i<=n; i++)
                      e[i-1]=e[i-1]/e[kk-1];
                    e[kk]=1.0f+e[kk];
                  }
                e[kk-1]=-e[kk-1];
                if ((kk+1<=m)&&(e[kk-1]!=0.0))
                  { for (i=kk+1; i<=m; i++) w[i-1]=0.0;
                    for (j=kk+1; j<=n; j++)
                      for (i=kk+1; i<=m; i++)
                        w[i-1]=w[i-1]+e[j-1]*a[(i-1)*n+j-1];
                    for (j=kk+1; j<=n; j++)
                      for (i=kk+1; i<=m; i++)
                        { ix=(i-1)*n+j-1;
                          a[ix]=a[ix]-w[i-1]*e[j-1]/e[kk];
                        }
                  }
                for (i=kk+1; i<=n; i++)
                  v[(i-1)*n+kk-1]=e[i-1];
              }
          }
      }
    mm=n;
    if (m+1<n) mm=m+1;
    if (k<n) s[k]=a[k*n+k];
    if (m<mm) s[mm-1]=0.0;
    if (l+1<mm) e[l]=a[l*n+mm-1];
    e[mm-1]=0.0;
    nn=m;
    if (m>n) nn=n;
    if (nn>=k+1)
      { for (j=k+1; j<=nn; j++)
          { for (i=1; i<=m; i++)
              u[(i-1)*m+j-1]=0.0;
            u[(j-1)*m+j-1]=1.0;
          }
      }
    if (k>=1)
      { for (ll=1; ll<=k; ll++)
          { kk=k-ll+1; iz=(kk-1)*m+kk-1;
            if (s[kk-1]!=0.0)
              { if (nn>=kk+1)
                  for (j=kk+1; j<=nn; j++)
                    { d=0.0;
                      for (i=kk; i<=m; i++)
                        { ix=(i-1)*m+kk-1;
                          iy=(i-1)*m+j-1;
                          d=d+u[ix]*u[iy]/u[iz];
                        }
                      d=-d;
                      for (i=kk; i<=m; i++)
                        { ix=(i-1)*m+j-1;
                          iy=(i-1)*m+kk-1;
                          u[ix]=u[ix]+d*u[iy];
                        }
                    }
                  for (i=kk; i<=m; i++)
                    { ix=(i-1)*m+kk-1; u[ix]=-u[ix];}
                  u[iz]=1.0f+u[iz];
                  if (kk-1>=1)
                    for (i=1; i<=kk-1; i++)
                      u[(i-1)*m+kk-1]=0.0;
              }
            else
              { for (i=1; i<=m; i++)
                  u[(i-1)*m+kk-1]=0.0;
                u[(kk-1)*m+kk-1]=1.0;
              }
          }
      }
    for (ll=1; ll<=n; ll++)
      { kk=n-ll+1; iz=kk*n+kk-1;
        if ((kk<=l)&&(e[kk-1]!=0.0))
          { for (j=kk+1; j<=n; j++)
              { d=0.0;
                for (i=kk+1; i<=n; i++)
                  { ix=(i-1)*n+kk-1; iy=(i-1)*n+j-1;
                    d=d+v[ix]*v[iy]/v[iz];
                  }
                d=-d;
                for (i=kk+1; i<=n; i++)
                  { ix=(i-1)*n+j-1; iy=(i-1)*n+kk-1;
                    v[ix]=v[ix]+d*v[iy];
                  }
              }
          }
        for (i=1; i<=n; i++)
          v[(i-1)*n+kk-1]=0.0;
        v[iz-n]=1.0;
      }
    for (i=1; i<=m; i++)
    for (j=1; j<=n; j++)
      a[(i-1)*n+j-1]=0.0;
    m1=mm; it=60;
    while (1==1)
      { if (mm==0)
          { ppp(a,e,s,v,m,n);
            free(s); free(e); free(w); return(1);
          }
        if (it==0)
          { ppp(a,e,s,v,m,n);
            free(s); free(e); free(w); return(-1);
          }
        kk=mm-1;
	while ((kk!=0)&&(fabs(e[kk-1])!=0.0))
          { d=(float)(fabs(s[kk-1])+fabs(s[kk]));
            dd=(float)fabs(e[kk-1]);
            if (dd>eps*d) kk=kk-1;
            else e[kk-1]=0.0;
          }
        if (kk==mm-1)
          { kk=kk+1;
            if (s[kk-1]<0.0)
              { s[kk-1]=-s[kk-1];
                for (i=1; i<=n; i++)
                  { ix=(i-1)*n+kk-1; v[ix]=-v[ix];}
              }
            while ((kk!=m1)&&(s[kk-1]<s[kk]))
              { d=s[kk-1]; s[kk-1]=s[kk]; s[kk]=d;
                if (kk<n)
                  for (i=1; i<=n; i++)
                    { ix=(i-1)*n+kk-1; iy=(i-1)*n+kk;
                      d=v[ix]; v[ix]=v[iy]; v[iy]=d;
                    }
                if (kk<m)
                  for (i=1; i<=m; i++)
                    { ix=(i-1)*m+kk-1; iy=(i-1)*m+kk;
                      d=u[ix]; u[ix]=u[iy]; u[iy]=d;
                    }
                kk=kk+1;
              }
            it=60;
            mm=mm-1;
          }
        else
          { ks=mm;
            while ((ks>kk)&&(fabs(s[ks-1])!=0.0))
              { d=0.0;
                if (ks!=mm) d=d+(float)fabs(e[ks-1]);
                if (ks!=kk+1) d=d+(float)fabs(e[ks-2]);
                dd=(float)fabs(s[ks-1]);
                if (dd>eps*d) ks=ks-1;
                else s[ks-1]=0.0;
              }
            if (ks==kk)
              { kk=kk+1;
                d=(float)fabs(s[mm-1]);
                t=(float)fabs(s[mm-2]);
                if (t>d) d=t;
                t=(float)fabs(e[mm-2]);
                if (t>d) d=t;
                t=(float)fabs(s[kk-1]);
                if (t>d) d=t;
                t=(float)fabs(e[kk-1]);
                if (t>d) d=t;
                sm=s[mm-1]/d; sm1=s[mm-2]/d;
                em1=e[mm-2]/d;
                sk=s[kk-1]/d; ek=e[kk-1]/d;
                b=((sm1+sm)*(sm1-sm)+em1*em1)/2.0f;
                c=sm*em1; c=c*c; shh=0.0;
                if ((b!=0.0)||(c!=0.0))
                  { shh=(float)sqrt(b*b+c);
                    if (b<0.0) shh=-shh;
                    shh=c/(b+shh);
                  }
                fg[0]=(sk+sm)*(sk-sm)-shh;
                fg[1]=sk*ek;
                for (i=kk; i<=mm-1; i++)
                  { sss(fg,cs);
                    if (i!=kk) e[i-2]=fg[0];
                    fg[0]=cs[0]*s[i-1]+cs[1]*e[i-1];
                    e[i-1]=cs[0]*e[i-1]-cs[1]*s[i-1];
                    fg[1]=cs[1]*s[i];
                    s[i]=cs[0]*s[i];
                    if ((cs[0]!=1.0)||(cs[1]!=0.0))
                      for (j=1; j<=n; j++)
                        { ix=(j-1)*n+i-1;
                          iy=(j-1)*n+i;
                          d=cs[0]*v[ix]+cs[1]*v[iy];
                          v[iy]=-cs[1]*v[ix]+cs[0]*v[iy];
                          v[ix]=d;
                        }
                    sss(fg,cs);
                    s[i-1]=fg[0];
                    fg[0]=cs[0]*e[i-1]+cs[1]*s[i];
                    s[i]=-cs[1]*e[i-1]+cs[0]*s[i];
                    fg[1]=cs[1]*e[i];
                    e[i]=cs[0]*e[i];
                    if (i<m)
                      if ((cs[0]!=1.0)||(cs[1]!=0.0))
                        for (j=1; j<=m; j++)
                          { ix=(j-1)*m+i-1;
                            iy=(j-1)*m+i;
                            d=cs[0]*u[ix]+cs[1]*u[iy];
                            u[iy]=-cs[1]*u[ix]+cs[0]*u[iy];
                            u[ix]=d;
                          }
                  }
                e[mm-2]=fg[0];
                it=it-1;
              }
            else
              { if (ks==mm)
                  { kk=kk+1;
                    fg[1]=e[mm-2]; e[mm-2]=0.0;
                    for (ll=kk; ll<=mm-1; ll++)
                      { i=mm+kk-ll-1;
                        fg[0]=s[i-1];
                        sss(fg,cs);
                        s[i-1]=fg[0];
                        if (i!=kk)
                          { fg[1]=-cs[1]*e[i-2];
                            e[i-2]=cs[0]*e[i-2];
                          }
                        if ((cs[0]!=1.0)||(cs[1]!=0.0))
                          for (j=1; j<=n; j++)
                            { ix=(j-1)*n+i-1;
                              iy=(j-1)*n+mm-1;
                              d=cs[0]*v[ix]+cs[1]*v[iy];
                              v[iy]=-cs[1]*v[ix]+cs[0]*v[iy];
                              v[ix]=d;
                            }
                      }
                  }
                else
                  { kk=ks+1;
                    fg[1]=e[kk-2];
                    e[kk-2]=0.0;
                    for (i=kk; i<=mm; i++)
                      { fg[0]=s[i-1];
                        sss(fg,cs);
                        s[i-1]=fg[0];
                        fg[1]=-cs[1]*e[i-1];
                        e[i-1]=cs[0]*e[i-1];
                        if ((cs[0]!=1.0)||(cs[1]!=0.0))
                          for (j=1; j<=m; j++)
                            { ix=(j-1)*m+i-1;
                              iy=(j-1)*m+kk-2;
                              d=cs[0]*u[ix]+cs[1]*u[iy];
                              u[iy]=-cs[1]*u[ix]+cs[0]*u[iy];
                              u[ix]=d;
                            }
                      }
                  }
              }
          }
      }
   	free(s);free(e);free(w); 
	  return(1);
  }
void ppp(float a[],float e[],float s[],float v[],int m,int n) 
{ int i,j,p,q;
    float d;
    if (m>=n) i=n;
    else i=m;
    for (j=1; j<=i-1; j++)
      { a[(j-1)*n+j-1]=s[j-1];
        a[(j-1)*n+j]=e[j-1];
      }
    a[(i-1)*n+i-1]=s[i-1];
    if (m<n) a[(i-1)*n+i]=e[i-1];
    for (i=1; i<=n-1; i++)
    for (j=i+1; j<=n; j++)
      { p=(i-1)*n+j-1; q=(j-1)*n+i-1;
        d=v[p]; v[p]=v[q]; v[q]=d;
      }
    return;
  }
  void sss(float fg[],float cs[])
 { float r,d;
    if ((fabs(fg[0])+fabs(fg[1]))==0.0)
      { cs[0]=1.0; cs[1]=0.0; d=0.0;}
    else 
      { d=(float)sqrt(fg[0]*fg[0]+fg[1]*fg[1]);
        if (fabs(fg[0])>fabs(fg[1]))
          { d=(float)fabs(d);
            if (fg[0]<0.0) d=-d;
          }
        if (fabs(fg[1])>=fabs(fg[0]))
          { d=(float)fabs(d);
            if (fg[1]<0.0) d=-d;
          }
        cs[0]=fg[0]/d; cs[1]=fg[1]/d;
      }
    r=1.0;
    if (fabs(fg[0])>fabs(fg[1])) r=cs[1];
    else
      if (cs[0]!=0.0) r=1.0f/cs[0];
    fg[0]=d; fg[1]=r;
    return;
  }

參考:

線性，非線性多項式：

http://blog.csdn.net/qll125596718/article/details/8248249

http://blog.csdn.net/poxiaozhuimeng/article/details/41117947

http://blog.csdn.net/ouyangying123/article/details/53996403

http://blog.csdn.net/jairuschan/article/details/7517773/

http://blog.csdn.net/zang141588761/article/details/50523036

http://www.cnblogs.com/gnuhpc/archive/2012/12/09/2809997.html

http://download.csdn.net/download/biaobiao11/9755119

圓擬合：

http://blog.sina.com.cn/s/blog_b27f71160101gxun.html  

http://www.cnblogs.com/dotLive/archive/2007/04/06/524633.html

http://blog.csdn.net/andylao62/article/details/24522365

http://blog.csdn.net/liyuanbhu/article/details/50889951

http://blog.csdn.net/liyuanbhu/article/details/50890587

橢圓擬合：

http://doc.okbase.net/u013708970/archive/121532.html