（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)

#include "stdafx.h"  
#include "GaussAlgorithm.h"  
  
VectorXf m_Vector_A;  
MatrixXf m_matrix_B;  
int m_iN =-1;  
  
bool InitData(int pSrc[100][100], int iWidth, int iHeight)  
{  
    if (NULL == pSrc || iWidth <=0 || iHeight <= 0)  
        return false;  
    m_iN = iWidth*iHeight;  
    VectorXf tmp_A(m_iN);  
    MatrixXf tmp_B(m_iN, 5);  
    int i =0, j=0, iPos =0;  
    while(i<iWidth)  
    {  
         j=0;  
        while(j<iHeight)  
        {  
            tmp_A(iPos) = pSrc[i][j] * log((float)pSrc[i][j]);  
  
            tmp_B(iPos,0) = pSrc[i][j] ;  
            tmp_B(iPos,1) = pSrc[i][j] * i;  
            tmp_B(iPos,2) = pSrc[i][j] * j;  
            tmp_B(iPos,3) = pSrc[i][j] * i * i;  
            tmp_B(iPos,4) = pSrc[i][j] * j * j;  
            ++iPos;  
            ++j;  
        }  
        ++i;  
    }  
    m_Vector_A = tmp_A;  
    m_matrix_B = tmp_B;  
  
}  
bool GetCentrePoint(float& x0,float& y0)  
{  
    if (m_iN<=0)  
        return false;  
    //QR分解  
    HouseholderQR<MatrixXf> qr;  
    qr.compute(m_matrix_B);  
    MatrixXf R = qr.matrixQR().triangularView<Upper>();  
    MatrixXf Q =  qr.householderQ();  
  
    //塊操作，獲取向量或矩陣的區域性  
    VectorXf S;  
    S = (Q.transpose()* m_Vector_A).head(5);  
    MatrixXf R1;  
    R1 = R.block(0,0,5,5);  
  
    VectorXf C;  
    C = R1.inverse() * S;  
  
    x0 = -0.5 * C[1] / C[3];  
    y0 = -0.5 * C[2] / C[4];  
  
    return true;  
} 
 

其中：

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

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