關於遙感映像的控制點片匹配演算法的一點想法（二）

     上文說到使用OpenCV進行模板匹配的函數matchTemplate，下面就matchTemplate函數的內部處理過程做一個簡單的說明。matchTemplate函數的原始碼在OpenCV的原始碼目錄下的 modules/imgproc/src/templmatch.cpp 檔案中。其核心函數代碼如下（其中的注釋是我添加的）：

void matchTemplate( const Mat& _img, const Mat& _templ, Mat& result, int method ){    CV_Assert( CV_TM_SQDIFF method && method CV_TM_CCOEFF_NORMED );

    //numType用來表示模板匹配的方式，0表示相關匹配法，1表示相關係數匹配法，2表示平方差匹配法    //isNormed表示是否進行歸一化處理，true表示進行歸一化，false表示不進行歸一化處理

    int numType = method == CV_TM_CCORR || method == CV_TM_CCORR_NORMED ? 0 :                  method == CV_TM_CCOEFF || method == CV_TM_CCOEFF_NORMED ? 1 : 2;    bool isNormed = method == CV_TM_CCORR_NORMED ||                    method == CV_TM_SQDIFF_NORMED ||                    method == CV_TM_CCOEFF_NORMED;

    //判斷兩幅映像的大小關係，如果輸入的原始映像比匹配映像要小，則將原始映像作為模板，原來的模板映像作為搜尋圖    Mat img = _img, templ = _templ;    if( img.rows templ.rows || img.cols templ.cols )        std::swap(img, templ);        CV_Assert( (img.depth() == CV_8U || img.depth() == CV_32F) &&               img.type() == templ.type() );
   //crossCorr函數是將輸入映像做了一次DFT變換（離散傅裡葉變換），將空間域的映像轉換到頻率域中來進行處理，並將處理的結果存放在result中    int cn = img.channels();    crossCorr( img, templ, result,               Size(img.cols - templ.cols + 1, img.rows - templ.rows + 1),               CV_32F, Point(0,0), 0, 0);
    //如果是相關匹配方法，此處已經計算完畢，返回
    if( method == CV_TM_CCORR )
        return;
    //將模板看作單位1，計算每一個像元所佔的百分比（也可以理解為整個模板面積為1，計算每個像元的面積）
    double invArea = 1./((double)templ.rows * templ.cols);
    Mat sum, sqsum;
    Scalar templMean, templSdv;
    double *q0 = 0, *q1 = 0, *q2 = 0, *q3 = 0;
    double templNorm = 0, templSum2 = 0;
    //相關係數匹配演算法
    if( method == CV_TM_CCOEFF )
    {
        integral(img, sum, CV_64F);//對原始映像進行求和
        templMean = mean(templ);//計算模板映像的均值向量
    }
    else//其他匹配演算法
    {
        integral(img, sum, sqsum, CV_64F);//計算原始映像的和以及平方和
        meanStdDev( templ, templMean, templSdv );//計算模板映像的均值向量和方差向量
        templNorm = CV_SQR(templSdv[0]) + CV_SQR(templSdv[1]) +
                    CV_SQR(templSdv[2]) + CV_SQR(templSdv[3]);//計算所有通道的方差和
        if( templNorm DBL_EPSILON && method == CV_TM_CCOEFF_NORMED )
        {//如果所有通道的方差的和等於0，並且使用的方法是歸一化相關係數匹配方法，則返回
            result = Scalar::all(1);
            return;
        }
        templSum2 = templNorm +
                     CV_SQR(templMean[0]) + CV_SQR(templMean[1]) +
                     CV_SQR(templMean[2]) + CV_SQR(templMean[3]);//計算所有通道的均值的平方和
        if( numType != 1 )//匹配方式不是相關係數，對模板均值向量和templNorm重新賦值
        {
            templMean = Scalar::all(0);
            templNorm = templSum2;
        }
        templSum2 /= invArea;
        templNorm = sqrt(templNorm);
        templNorm /= sqrt(invArea); // care of accuracy here
        q0 = (double*)sqsum.data;
        q1 = q0 + templ.cols*cn;
        q2 = (double*)(sqsum.data + templ.rows*sqsum.step);
        q3 = q2 + templ.cols*cn;
    }
    //下面就是在結果映像中進行尋找匹配的結果位置，代碼略去，具體可參考OpenCV原始碼
 
    在測試過程中，發現直接用OpenCV開啟超過1G（可能比這個大小還要小）的映像就會導致錯誤，查看OpenCV代碼，發現其一次將所有的映像資料全部載入記憶體，導致記憶體不夠出錯，而且OpenCV只支援普通的常用的映像資料格式，並不能直接開啟Erdas的img格式或者PCI的pix格式，針對這兩個問題，我提出的解決思路是實用GDAL來作為映像資料的讀取庫，然後分塊讀取一部分資料（指定一個大小，我指定的是每次讀64M的資料），如果映像資料超過該大小，則分塊處理，每次處理一塊，OpenCV可以建立一個記憶體資料，就是將映像的RGB值按照特定的順序讀取到記憶體中，然後將該指標交給OpenCV，OpenCV可以通過該記憶體資料來構建一個映像，然後再調用模板匹配演算法，計算出的結果在還原到原來大的映像中的位置。
    在實際的執行過程中，上述方式是可行的，執行效率大概在一個1G的img映像尋找一個128*128的模板需要1分鐘左右，相比按照前一篇部落格中的方式速度提升了很多倍（1G的資料大概計算了4個小時還沒有算完）。對於1min中尋找出結果應該是可以接受的，但是感覺還是有點慢，於是想到首先將模板和原始映像同時按照一定的採樣比例建立金字塔，然後從金字塔的頂層開始進行匹配，根據頂層找到的位置，再到下一層對應的地區中尋找，依次知道尋找到原始映像層級。這樣就可以減少很多的運算。經過測試，使用這種方式，還是上面的1G的資料，尋找只要10S鐘即可完成，而且尋找的結果也是很準確的。
    在編寫程式的過程中，可以使用GDAL的RasteIO函數，直接使用最後面的三個參數，可以直接將映像中的資料讀取成OpenCV中的映像儲存方式 BGRBGRBGR…
   對於OpenCV在進行第一步匹配的時候，是在頻率域中進行處理的，相關的函數也在OpenCV的原始碼目錄下的 modules/imgproc/src/templmatch.cpp 檔案中，關於DFT變換（離散傅裡葉變換）比較複雜，理論知識尤其是需要很好的數學功底才能夠理解，感興趣的童鞋可以google或者參考相關的數位影像處理的書籍，一般的數位影像處理的書籍中都會有的，因為DFT和FFT在數位影像處理以及訊號處理方面使用的非常廣泛。
Technorati 標籤: OpenCV,模板匹配