高斯拟合

高斯拟合(GaussianFitting)即使用形如：Gi(x)=Ai*exp((x-Bi)^2/Ci^2)的高斯函数对数据点集进行函数逼近的拟合方法。其实可以跟多项式拟合类比起来，不同的是多项式拟合是用幂函数系，而高斯拟合是用高斯函数系。使用高斯函数来进行拟合，优点在于计算积分十分简单快捷。这一点在很多领域都有应用，特别是计算化学。著名的化学软件Gaussian98就是建立在高斯基函数拟合的数学基础上的。高斯拟合算法解决思路发布于：2014-02-2523:11:14浏览：0次0高斯拟合算法请问，如何用c++实现高斯拟合算法？------解决方案--------------------#includeiostream.h#includemath.h#includestdlib.h#includewindows.hdoublef(intn,doublex){//f(n,x)用来返回x的n次方doubley=1.0;if(n==0)return1.0;else{for(inti=0;in;i++)y*=x;returny;}}intxianxingfangchengzu(double**a,intn,double*b,double*p,doubledt)//用高斯列主元法来求解法方程组{inti,j,k,l;doublec,t;for(k=1;k=n;k++){c=0.0;for(i=k;i=n;i++)if(fabs(a[i-1][k-1])fabs(c)){c=a[i-1][k-1];l=i;}if(fabs(c)=dt)return(0);if(l!=k){for(j=k;j=n;j++){t=a[k-1][j-1];a[k-1][j-1]=a[l-1][j-1];a[l-1][j-1]=t;}t=b[k-1];b[k-1]=b[l-1];b[l-1]=t;}c=1/c;for(j=k+1;j=n;j++){a[k-1][j-1]=a[k-1][j-1]*c;for(i=k+1;i=n;i++)a[i-1][j-1]-=a[i-1][k-1]*a[k-1][j-1];}b[k-1]*=c;for(i=k+1;i=n;i++)b[i-1]-=b[k-1]*a[i-1][k-1];}for(i=n;i=1;i--)for(j=i+1;j=n;j++)b[i-1]-=b[j-1]*a[i-1][j-1];cout.precision(12);for(i=0;in;i++)p[i]=b[i];}double**create(inta,intb)//动态生成数组{double**P=newdouble*[a];for(inti=0;ib;i++)P[i]=newdouble[b];returnP;}voidzuixiaoerchengnihe(doublex[],doubley[],intn,doublea[],intm){inti,j,k,l;double**A,*B;A=create(m,m);B=newdouble[m];for(i=0;im;i++)for(j=0;jm;j++)A[i][j]=0.0;for(k=0;km;k++)for(l=0;lm;l++)for(j=0;jn;j++)A[k][l]+=f(k,x[j])*f(l,x[j]);//计算法方程组系数矩阵A[k][l]cout法方程组的系数矩阵为:endl;for(i=0;im;i++)for(j=0,k=1;jm;j++,k++){coutA[i][j]'\t';if(k&&k%m==0)coutendl;}for(i=0;im;i++)B[i]=0.0;for(i=0;im;i++)for(j=0;jn;j++)B[i]+=y[j]*f(i,x[j]);for(i=0;im;i++)coutB[i]=B[i]endl;//记录B[n]xianxingfangchengzu(A,m,B,a,1e-6);delete[]A;deleteB;}doublepingfangwucha(doublex[],doubley[],intn,doublea[],intm)//计算最小二乘解的平方误差{doubledeta,q=0.0,r=0.0;inti,j;double*B;B=newdouble[m];for(i=0;im;i++)B[i]=0.0;for(i=0;im;i++)for(j=0;jn;j++)B[i]+=y[j]*f(i,x[j]);for(i=0;in;i++)q+=y[i]*y[i];for(j=0;jm;j++)r+=a[j]*B[j];deta=fabs(q-r);returndeta;deleteB;}voidmain(void){inti,n,m;double*x,*y,*a;charch='y';do{system(cls);cout请输入所给拟合数据点的个数n=;cinn;cout请输入所要拟合多项式的项数m=;cinm;while(n=m){cout你所输入的数据点无法确定拟合项数，请重新输入endl;Sleep(1000);system(cls);cout请输入所给拟合数据点的个数n=;cinn;cout请输入所要拟合多项式的项数m=;cinm;}x=newdouble[n];//存放数据点xy=newdouble[n];//存放数据点ya=newdouble[m];//存放拟合多项式的系数cout请输入所给定的n个数据xendl;for(i=0;in;i++){coutx[i+1]=;cinx[i];}cout请输入所给定的n个数据yendl;for(i=0;in;i++){couty[i+1]=;ciny[i];}zuixiaoerchengnihe(x,y,n,a,m+1);coutendl;cout拟合多项式的系数为:endl;for(i=0;i=m;i++)couta[i]=a[i]'\t';coutendl;cout平方误差为：pingfangwucha(x,y,n,a,m+1)endl;deletex;deletey;cout按y继续，按其他字符退出endl;cinch;}while(ch=='y'------解决方案--------------------（1）二维高斯去曲面拟合推导一个二维高斯方程可以写成如下形式：其中，G为高斯分布的幅值，,为x,y方向上的标准差，对式（1）两边取对数，并展开平方项，整理后为：假如参与拟合的数据点有N个，则将这个N个数据点写成矩阵的形式为：A=BC，其中：A为N*1的向量，其元素为：B为N*5的矩阵：C为一个由高斯参数组成的向量：（2）求解二维高斯曲线拟合N个数据点误差的列向量为：E=A-BC，用最小二乘法拟合，使其N个数据点的均方差最小，即：在图像数据处理时，数据量比较大，为减小计算量，将矩阵B进行QR分解，即：B=QR，分解后Q为一个N*N的正交矩阵，R为一个N*5的上三角矩阵，对E=A-BC进行如下推导：由于Q为正交矩阵，可以得到：令：其中：S为一个5维列向量；T为一个N-5维列向量；R1为一个5*5的上三角方阵，则上式中，当S=R1C时取得最小值，因此只需解出：即可求出：中的这些参数，这里先给出：这里：（3）C++代码实现，算法的实现过程中由于涉及大量的矩阵运算，所以采用了第三方的开源矩阵算法Eigen，这里真正用于高斯拟合的函数是boolGetCentrePoint(float&x0,float&y0)关于Eigen的用法请参考：[cpp]viewplaincopy1.#includestdafx.h2.#includeGaussAlgorithm.h3.4.VectorXfm_Vector_A;5.MatrixXfm_matrix_B;6.intm_iN=-1;7.8.boolInitData(intpSrc[100][100],intiWidth,intiHeight)9.{10.if(NULL==pSrc||iWidth=0||iHeight=0)11.returnfalse;12.m_iN=iWidth*iHeight;13.VectorXftmp_A(m_iN);14.MatrixXftmp_B(m_iN,5);15.inti=0,j=0,iPos=0;16.while(iiWidth)17.{18.j=0;19.while(jiHeight)20.{21.tmp_A(iPos)=pSrc[i][j]*log((float)pSrc[i][j]);22.23.tmp_B(iPos,0)=pSrc[i][j];24.tmp_B(iPos,1)=pSrc[i][j]*i;25.tmp_B(iPos,2)=pSrc[i][j]*j;26.tmp_B(iPos,3)=pSrc[i][j]*i*i;27.tmp_B(iPos,4)=pSrc[i][j]*j*j;28.++iPos;29.++j;30.}31.++i;32.}33.m_Vector_A=tmp_A;34.m_matrix_B=tmp_B;35.36.}37.boolGetCentrePoint(float&x0,float&y0)38.{39.if(m_iN=0)40.returnfalse;41.//QR分解42.HouseholderQRMatrixXfqr;43.qr.compute(m_matrix_B);44.MatrixXfR=qr.matrixQR().triangularViewUpper();45.MatrixXfQ=qr.householderQ();46.47.//块操作，获取向量或矩阵的局部48.VectorXfS;49.S=(Q.transpose()*m_Vector_A).head(5);50.MatrixXfR1;51.R1=R.block(0,0,5,5);52.53.VectorXfC;54.C=R1.inverse()*S;55.56.x0=-0.5*C[1]/C[3];57.y0=-0.5*C[2]/C[4];58.59.returntrue;60.}其中：函数boolInitData(intpSrc[100][100],intiWidth,intiHeight)主要用于将数组转换成相应的矩阵，因为我的所有数据都在一个整形的二维数组中存着，所以需要做这种转换。函数boolGetCentrePoint(float&x0,float&y0)主要用于对数据点进行二维高斯曲面拟合，并返回拟合的光点中心。