高斯拟合Word文档下载推荐.docx
《高斯拟合Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《高斯拟合Word文档下载推荐.docx（11页珍藏版）》请在冰豆网上搜索。
for（int
i=0;
i<
n;
i++）y*=x;
return
y;
}
}
int
xianxingfangchengzu（double
**a,int
*b,double
*p,double
dt）//用高斯列主元法来求解法方程组
{
i,j,k,l;
c,t;
for（k=1;
k<
=n;
k++）
c=0.0;
for（i=k;
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<
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;
a[k-1][j-1]=a[k-1][j-1]*c;
for（i=k+1;
a[i-1][j-1]-=a[i-1][k-1]*a[k-1][j-1];
b[k-1]*=c;
b[i-1]-=b[k-1]*a[i-1][k-1];
for（i=n;
i>
=1;
i--）
for（j=i+1;
b[i-1]-=b[j-1]*a[i-1][j-1];
cout.precision（12）;
for（i=0;
i++）p[i]=b[i];
double**
create（int
a,int
b）//动态生成数组
**P=new
*[a];
for（int
b;
P[i]=new
double[b];
return
P;
void
zuixiaoerchengnihe（double
x[],double
y[],int
a[],int
m）
int
**A,*B;
A=create（m,m）;
B=new
double[m];
for（i=0;
m;
for（j=0;
j++）A[i][j]=0.0;
for（k=0;
for（l=0;
l<
l++）
for（j=0;
j++）A[k][l]+=f（k,x[j]）*f（l,x[j]）;
//计算法方程组系数矩阵A[k][l]
cout<
"
法方程组的系数矩阵为:
endl;
for（j=0,k=1;
j++,k++）{
cout<
A[i][j]<
'
\t'
;
if（k&
&
k%m==0）cout<
}
i++）B[i]=0.0;
for（j=0;
j++）B[i]+=y[j]*f（i,x[j]）;
i++）cout<
B["
]="
B[i]<
//记录B[n]
xianxingfangchengzu（A,m,B,a,1e-6）;
delete[]A;
delete
B;
pingfangwucha（double
m）//计算最小二乘解的平方误差
deta,q=0.0,r=0.0;
i,j;
*B;
i++）q+=y[i]*y[i];
j++）r+=a[j]*B[j];
deta=fabs（q-r）;
deta;
}
main（void）{
i,n,m;
*x,*y,*a;
char
ch='
y'
do{
system（"
cls"
）;
请输入所给拟合数据点的个数n="
cin>
>
请输入所要拟合多项式的项数m="
while（n<
=m）{
cout<
你所输入的数据点无法确定拟合项数，请重新输入"
Sleep（1000）;
system（"
cin>
x=new
double[n];
//存放数据点x
y=new
//存放数据点y
a=new
//存放拟合多项式的系数
请输入所给定的"
n<
个数据x"
{
cout<
x["
i+1<
cin>
x[i];
个数据y"
y["
y[i];
zuixiaoerchengnihe（x,y,n,a,m+1）;
拟合多项式的系数为:
=m;
a["
a[i]<
平方误差为：
pingfangwucha（x,y,n,a,m+1）<
x;
delete
按y继续，按其他字符退出"
ch;
}while（ch=='
------解决方案--------------------
（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]
viewplaincopy
1.#include
stdafx.h"
2.#include
GaussAlgorithm.h"
3.
4.VectorXf
m_Vector_A;
5.MatrixXf
m_matrix_B;
6.int
m_iN
=-1;
7.
8.bool
InitData（int
pSrc[100][100],
iWidth,
iHeight）
9.{
10.
if
（NULL
==
pSrc
||
iWidth
=0
iHeight
=
0）
11.
return
false;
12.
iWidth*iHeight;
13.
VectorXf
tmp_A（m_iN）;
14.
MatrixXf
tmp_B（m_iN,
5）;
15.
i
=0,
j=0,
iPos
=0;
16.
while（i<
iWidth）
17.
{
18.
j=0;
19.
while（j<
20.
21.
tmp_A（iPos）
pSrc[i][j]
*
log（（float）pSrc[i][j]）;
22.
23.
tmp_B（iPos,0）
24.
tmp_B（iPos,1）
i;
25.
tmp_B（iPos,2）
j;
26.
tmp_B（iPos,3）
27.
tmp_B（iPos,4）
j
28.
++iPos;
29.
++j;
30.
}
31.
++i;
32.
33.
m_Vector_A
tmp_A;
34.
m_matrix_B
tmp_B;
35.
36.}
37.bool
GetCentrePoint（float&
x0,float&
y0）
38.{
39.
（m_iN<
=0）
40.
41.
//QR分解
42.
HouseholderQR<
MatrixXf>
qr;
43.
pute（m_matrix_B）;
44.
R
qr.matrixQR（）.triangularView<
Upper>
（）;
45.
Q
qr.householderQ（）;
46.
47.
//块操作，获取向量或矩阵的局部
48.
S;
49.
S
（Q.transpose（）*
m_Vector_A）.head（5）;
50.
R1;
51.
R1
R.block（0,0,5,5）;
52.
53.
C;
54.
C
R1.inverse（）
55.
56.
x0
-0.5
C[1]
/
C[3];
57.
y0
C[2]
C[4];
58.
59.
true;
60.}
函数
boolInitData（intpSrc[100][100],intiWidth,intiHeight）主要用于将数组转换成相应的矩阵，因为我的所有数据都在一个整形的二维数组中存着，所以需要做这种转换。
函数boolGetCentrePoint（float&
y0）主要用于对数据点进行二维高斯曲面拟合，并返回拟合的光点中心。
升级会员 