The following code based on openCV's phaseCorrelateRes()
will do correlation in 2 dimensions.
static void fftShift(InputOutputArray _out)
{
Mat out = _out.getMat();
if(out.rows == 1 && out.cols == 1)
{
// trivially shifted.
return;
}
vector<Mat> planes;
split(out, planes);
int xMid = out.cols >> 1;
int yMid = out.rows >> 1;
bool is_1d = xMid == 0 || yMid == 0;
if(is_1d)
{
xMid = xMid + yMid;
for(size_t i = 0; i < planes.size(); i++)
{
Mat tmp;
Mat half0(planes[i], Rect(0, 0, xMid, 1));
Mat half1(planes[i], Rect(xMid, 0, xMid, 1));
half0.copyTo(tmp);
half1.copyTo(half0);
tmp.copyTo(half1);
}
}
else
{
for(size_t i = 0; i < planes.size(); i++)
{
// perform quadrant swaps...
Mat tmp;
Mat q0(planes[i], Rect(0, 0, xMid, yMid));
Mat q1(planes[i], Rect(xMid, 0, xMid, yMid));
Mat q2(planes[i], Rect(0, yMid, xMid, yMid));
Mat q3(planes[i], Rect(xMid, yMid, xMid, yMid));
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
}
}
merge(planes, out);
}
void Correlate2d(
const cv::Mat& src1,
const cv::Mat& src2,
cv::Mat& dst,
double* response)
{
CV_Assert( src1.type() == src2.type());
CV_Assert( src1.type() == CV_32FC1 || src1.type() == CV_64FC1 );
CV_Assert( src1.size == src2.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));
}
else
{
padded1 = src1;
padded2 = src2;
}
Mat FFT1, FFT2, P, Pm, C;
// correlation equation
// Reference: http://en.wikipedia.org/wiki/Phase_correlation
dft(padded1, FFT1, DFT_REAL_OUTPUT);
dft(padded2, FFT2, DFT_REAL_OUTPUT);
mulSpectrums(FFT1, FFT2, dst, 0, true);
idft(dst, dst, DFT_SCALE); // gives us the correlation result...
fftShift(dst); // shift the energy to the center of the frame.
// locate the highest peak
Point peakLoc;
minMaxLoc(dst, NULL, NULL, NULL, &peakLoc);
// max response is scaled
if( response )
*response = dst.at<float>(peakLoc);
}
You can find the code in \opencv\sources\modules\imgproc\src\phasecorr.cpp
In order to change the code to convolution simply change this line:
mulSpectrums(FFT1, FFT2, dst, 0, true);
to
mulSpectrums(FFT1, FFT2, dst, 0, false);
This is equivalent to doing in matlab:
dst = fftshift(ifft2(fft2(src1).*conj(fft2(src2))))