1. 程式人生 > >相位相關法——Opencv原始碼分析——phaseCorrelate

相位相關法——Opencv原始碼分析——phaseCorrelate

相位相關法(phase correlate)可以用於檢測兩幅內容相同的影象之間的相對位移量。它是基於傅立葉變換的位移定理:一個平移過的函式的傅立葉變換僅僅是未平移函式的傅立葉變換與一個具有線性相位的指數因子的乘積,即空間域中的平移會造成頻域中頻譜的相移。它的公式定義為:設二維函式(影象)f(x,y)的傅立葉變換為F(u,v),即DFT[f(x,y)]=F(u,v),如果f(x,y)平移(a,b),則平移後的傅立葉變換為:

        (1)

因此,當兩幅函式f1(x,y)和f2(x,y)僅僅有位移的差異,即f2(x,y)= f1(x-a,y-b),則它們的傅立葉變換F1(u,v)和F2(u,v)有如下關係:

         (2)

由上式很容易得到f1(x,y)和f2(x,y)的互功率譜為(這裡還用到了f1(x,y)和f2(x,y)的頻譜的模相等):

         (3)

式中F*表示F的共軛,上式表示平移定理保證了互功率譜的相位等於兩幅影象之間的相移。

Opencv的文件給出了詳細的用相位相關法求解位移量的過程,

1、對待處理的兩幅影象src1和src2應用窗函式去除影象的邊界效應,文件中推薦使用漢寧窗,它可用createHanningWindow函式生成;

2、求傅立葉變換:Ga=DFT[scr1]和Ga=DFT[scr1];

3、計算互功率譜:

4、對互功率譜求傅立葉逆變換:r=DFT-1[R];

5、對r計算最大值的位置,並在以該位置為中心的5×5的窗體內應用下列公式獲得亞畫素級的精度位置:

                (4)

最終(a,b)為兩個影象之間的位移量。

 

phaseCorrelate函式的原型為:

Point2d phaseCorrelate(InputArray src1,InputArray src2, InputArray window=noArray())

src1和src2為兩個要比較的影象,window為步驟1中所要用到的窗函式,該引數可選,預設情況下不使用任何窗函式,函式返回是具有亞畫素級精度的位移量。

下面就給出具體的原始碼分析,該函式在source/modules/imgproc/src/phasecorr.cpp檔案內,

cv::Point2d cv::phaseCorrelate(InputArray _src1, InputArray _src2, InputArray _window)
{
    return phaseCorrelateRes(_src1, _src2, _window, 0);
}

phaseCorrelateRes函式的第4個引數表示的是最大響應值,也就是公式4中分母部分的歸一化後的值,在這裡我們沒有用到它。

cv::Point2d cv::phaseCorrelateRes(InputArray _src1, InputArray _src2, InputArray _window, double* response)
{
    //分別得到兩幅輸入影象和窗函式的矩陣形式
    Mat src1 = _src1.getMat();
    Mat src2 = _src2.getMat();
    Mat window = _window.getMat();
    //輸入影象的型別和大小的判斷,必須一致,而且型別必須是32位或64位浮點灰度影象
    CV_Assert( src1.type() == src2.type());
    CV_Assert( src1.type() == CV_32FC1 || src1.type() == CV_64FC1 );
    CV_Assert( src1.size == src2.size);
    //如果使用窗函式,則窗函式的大小和型別必須與輸入影象的一致
    if(!window.empty())
    {
        CV_Assert( src1.type() == window.type());
        CV_Assert( src1.size == window.size);
    }
    //因為要進行離散傅立葉變換,所以為了提高效率,就要得到最佳的影象尺寸
    int M = getOptimalDFTSize(src1.rows);
    int N = getOptimalDFTSize(src1.cols);
 
    Mat padded1, padded2, paddedWin;
    //生成尺寸修改以後的矩陣
    if(M != src1.rows || N != src1.cols)   //最佳尺寸不是原影象的尺寸
    {
        //通過補零的方式填充多出來的畫素
        copyMakeBorder(src1, padded1, 0, M - src1.rows, 0, N - src1.cols, BORDER_CONSTANT, Scalar::all(0));
        copyMakeBorder(src2, padded2, 0, M - src2.rows, 0, N - src2.cols, BORDER_CONSTANT, Scalar::all(0));
 
        if(!window.empty())
        {
            copyMakeBorder(window, paddedWin, 0, M - window.rows, 0, N - window.cols, BORDER_CONSTANT, Scalar::all(0));
        }
    }
    else    //最佳尺寸與原影象的尺寸一致
    {
        padded1 = src1;
        padded2 = src2;
        paddedWin = window;
    }
 
    Mat FFT1, FFT2, P, Pm, C;
 
    // perform window multiplication if available
    //執行步驟1,兩幅輸入影象分別與窗函式逐點相乘
    if(!paddedWin.empty())
    {
        // apply window to both images before proceeding...
        multiply(paddedWin, padded1, padded1);
        multiply(paddedWin, padded2, padded2);
    }
 
    // execute phase correlation equation
    // Reference: http://en.wikipedia.org/wiki/Phase_correlation
    //執行步驟2,分別對兩幅影象取傅立葉變換
    dft(padded1, FFT1, DFT_REAL_OUTPUT);
    dft(padded2, FFT2, DFT_REAL_OUTPUT);
    //執行步驟3
    //計算互功率譜的分子部分,即公式3中的分子,其中P為輸出結果,true表示的是對FF2取共軛,所以得到的結果為:P=FFT1×FFT2*,mulSpectrums函式為通用函式
    mulSpectrums(FFT1, FFT2, P, 0, true);
    //計算互功率譜的分母部分,即公式3中的分母,結果為:Pm=|P|,magSpectrums函式就是在phasecorr.cpp檔案內給出的,它的作用是對複數取模。
    magSpectrums(P, Pm);
    //計算互功率譜,即公式3,結果為:C=P / Pm,divSpectrums函式也是在phasecorr.cpp檔案內給出的,它仿照mulSpectrums函式的寫法,其中引數false表示不取共軛
    divSpectrums(P, Pm, C, 0, false); // FF* / |FF*| (phase correlation equation completed here...)
    //執行步驟4,傅立葉逆變換
    idft(C, C); // gives us the nice peak shift location...
    /*平移處理,fftShift函式也是在phasecorr.cpp檔案內給出的,它的作用是把影象平均分割成——左上、左下、右上、右下,把左上和右下對調,把右上和左下對調。它的目的是把能量調整到影象的中心,也就是影象的中心對應於兩幅影象相頻差為零的地方,即沒有發生位移的地方。*/
    fftShift(C); // shift the energy to the center of the frame.
 
    //執行步驟5
    // locate the highest peak
    //找到最大點處的畫素位置,minMaxLoc為通用函式
    Point peakLoc;
    minMaxLoc(C, NULL, NULL, NULL, &peakLoc);
 
    // get the phase shift with sub-pixel accuracy, 5x5 window seems about right here...
    //在5×5的窗體內確定亞畫素精度的座標位置
    Point2d t;
    // weightedCentroid也是在phasecorr.cpp檔案內給出的,它是利用公式4來計算更精確的座標位置
    t = weightedCentroid(C, peakLoc, Size(5, 5), response);
 
    // max response is M*N (not exactly, might be slightly larger due to rounding errors)
    //求最大響應值
    if(response)
        *response /= M*N;
 
    // adjust shift relative to image center...
    //最終確定位移量
    //先找到影象中點,然後用中點減去由步驟5得到的座標位置
    Point2d center((double)padded1.cols / 2.0, (double)padded1.rows / 2.0);
 
    return (center - t);
}

phaseCorrelate函式的原始碼可讀性較強,而且它的原理在文件中說明得也十分清晰,因此對相位相關的理解難度應該不大。下面就給出具體的實現程式。我是對彩色影象進行測試,因為phaseCorrelate函式只能處理32位或64位浮點型的單通道影象,因此如果不是處理這類影象,是需要進行轉換的。
 

#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>
using namespace cv;
using namespace std;
int main( int argc, char** argv )
{
	Mat src1, src2;
	Mat dst1, dst2;
	if( argc != 3 || !(src1=imread(argv[1], 1)).data || !(src2=imread(argv[2], 1)).data)
		return -1;
 
	cvtColor( src1, src1, CV_BGR2GRAY );     //轉換為灰度影象
        src1.convertTo(dst1,CV_32FC1);       //轉換為32位浮點型
	cvtColor( src2, src2, CV_BGR2GRAY );
        src2.convertTo(dst2,CV_32FC1);
 
	Point2d phase_shift;
        phase_shift = phaseCorrelate(dst1,dst2);
        cout<<endl<<"warp :"<<endl<<"\tX shift : "<<phase_shift.x<<"\tY shift : "<<phase_shift.y<<endl;
 
        waitKey(0);
 
        return 0;
}

最終得到的結果是:

warp:

      X  shift :  3.803         Y  shift :  21.8661


我們來驗證一下上述結果是否正確。我用的測試影象(如下所示)的原始大小為1761×1761,兩幅影象的相對位移量為(16,96),然後我把兩幅影象的大小調整為400×400,這樣經過計算,相對位移量應該是(3.6343,21.8058)。可以看出兩者之間還是存在一定的誤差的,但我認為誤差還是能夠在可接受的範圍內的。