（1）二維高斯去曲面擬合推導

一個二維高斯方程可以寫成如下形式：

其中，G為高斯分佈的幅值，,為x,y方向上的標準差，對式（1）兩邊取對數，並展開平方項，整理後為：

假如參與擬合的資料點有N個，則將這個N個數據點寫成矩陣的形式為：A = B C，

其中：

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 = R1C時取得最小值，因此只需解出：

即可求出：

中的

這些引數，這裡先給出：

這裡：

（3）C++程式碼實現，演算法的實現過程中由於涉及大量的矩陣運算，所以採用了第三方的開源矩陣演算法Eigen，這裡真正用於高斯擬合的函式是

bool GetCentrePoint(float& x0,float& y0)

關於Eigen的用法請參考：http://blog.csdn.NET/hjx_1000/article/details/8490941

1. #include "stdafx.h"
2. #include "GaussAlgorithm.h"
3. VectorXf m_Vector_A;  
4. MatrixXf m_matrix_B;  
5. int m_iN =-1;  
6. bool InitData(int pSrc[100][100], int iWidth, int iHeight)  
7. {  
8.     if (NULL == pSrc || iWidth <=0 || iHeight <= 0)  
9.         returnfalse;  
10.     m_iN = iWidth*iHeight;  
11.     VectorXf tmp_A(m_iN);  
12.     MatrixXf tmp_B(m_iN, 5);  
13.     int i =0, j=0, iPos =0;  
14.     while(i<iWidth)  
15.     {  
16.          j=0;  
17.         while(j<iHeight)  
18.         {  
19.             tmp_A(iPos) = pSrc[i][j] * log((float)pSrc[i][j]);  
20.             tmp_B(iPos,0) = pSrc[i][j] ;  
21.             tmp_B(iPos,1) = pSrc[i][j] * i;  
22.             tmp_B(iPos,2) = pSrc[i][j] * j;  
23.             tmp_B(iPos,3) = pSrc[i][j] * i * i;  
24.             tmp_B(iPos,4) = pSrc[i][j] * j * j;  
25.             ++iPos;  
26.             ++j;  
27.         }  
28.         ++i;  
29.     }  
30.     m_Vector_A = tmp_A;  
31.     m_matrix_B = tmp_B;  
32. }  
33. bool GetCentrePoint(float& x0,float& y0)  
34. {  
35.     if (m_iN<=0)  
36.         returnfalse;  
37.     //QR分解
38.     HouseholderQR<MatrixXf> qr;  
39.     qr.compute(m_matrix_B);  
40.     MatrixXf R = qr.matrixQR().triangularView<Upper>();  
41.     MatrixXf Q =  qr.householderQ();  
42.     //塊操作，獲取向量或矩陣的區域性
43.     VectorXf S;  
44.     S = (Q.transpose()* m_Vector_A).head(5);  
45.     MatrixXf R1;  
46.     R1 = R.block(0,0,5,5);  
47.     VectorXf C;  
48.     C = R1.inverse() * S;  
49.     x0 = -0.5 * C[1] / C[3];  
50.     y0 = -0.5 * C[2] / C[4];  
51.     returntrue;  
52. }  

其中：

函式 bool InitData(int pSrc[100][100], int iWidth, int iHeight)主要用於將陣列轉換成相應的矩陣，因為我的所有資料都在一個整形的二維陣列中存著，所以需要做這種轉換。

函式bool GetCentrePoint(float& x0,float& y0)主要用於對資料點進行二維高斯曲面擬合，並返回擬合的光點中心。

原文地址：http://blog.csdn.net/houjixin/article/details/8490653/