Reputation: 2176
Phase correlation with Open Cv
I want to calculate phase correlation of function in OpenCV. I've found http://nashruddin.com/phase-correlation-function-in-opencv.html but there is used some other library and I want to do it with OpenCV.
My problem is differen cause I want to calculate phase correlation of two arrays, both 360 elements. I tried to evaluate from documentation how to do it, but I don't know if my method is good. My R matrix isn't normalized, how to normalize it and is it necessary? If you give me some examples I would be grateful. My function to do this task:
void calcPhaseCorrelation(int* x1, int* x2){
Mat array1 = Mat(1,360,DataType<float>::type);
Mat array2 = Mat(1,360,DataType<float>::type);
uchar* begin = array1.data;
uchar* end = begin + (array1.step.p[0]/ sizeof(float)) * array1.size().height;
uchar *ptr = begin;
int ctr1 = 0, ctr2 = 0; //control in loops
while(ptr<end)
{
*ptr = (float)x1[ctr1];
ctr1++;
ptr++;
}
begin = array2.data;
end = begin + (array2.step.p[0]/ sizeof(float)) * array2.size().height;
ptr = begin;
while(ptr<end)
{
*ptr = x2[ctr2];
ctr2++;
ptr++;
}
Mat outputArray;
outputArray.create(abs(array1.rows - array2.rows)+1,
abs(array1.cols - array2.cols)+1, array1.type());
Size dftSize;
dftSize.width = getOptimalDFTSize(array1.cols + array2.cols - 1);
dftSize.height = getOptimalDFTSize(array1.rows + array2.rows - 1);
Mat resultA(dftSize, array1.type(), Scalar::all(0));
Mat resultB(dftSize, array2.type(), Scalar::all(0));
dft(array1,resultA);
dft(array2,resultB);
Mat R;
mulSpectrums(resultA,resultB,R,DFT_ROWS,true);
Mat NormR;
normalize(R,NormR);
idft(NormR,outputArray);
double minVal, maxVal;
Point minLoc, maxLoc;
minMaxLoc(outputArray,&minVal,&maxVal,&minLoc,&maxLoc);
std::cout<<"Min value: "<<minVal<<", max value: "<<maxVal<<std::endl;
std::cout<<"Min loc x: "<<minLoc.x<<", min loc y: "<<minLoc.y<<std::endl;
std::cout<<"Max loc x: "<<maxLoc.x<<", max loc y: "<<maxLoc.y<<std::endl;
}
I know that code is not clear etc., but it's only fast test. I want to know if method is correct. But also there every advice is appreciated.
//Edit: I used mevatron code and get #QNAN also while debugging seen that in function that it does not find a peak, value is (-1,-1). Function I'm using with new OpenCV correlation function:
void phaseCorrelationOpenCvTrunk(int* array1, int* array2)
{
Mat hann;
vector<double> arr1, arr2;
for(int i = 0; i < 360; i++)
{
arr1.push_back((double)array1[i]);
arr2.push_back((double)array2[i]);
}
Mat firstArray = Mat(arr1);
Mat secondArray = Mat(arr2);
std::cout<<"Type control: "<<firstArray.type()<<std::endl;
createHanningWindow(hann, firstArray.size(), CV_64F);
Point2d shift = phaseCorrelate(firstArray, secondArray,hann);
std::cout<<"shift: "<<shift.x<<";"<<shift.y<<std::endl;
}
What am I doing wrong?
Upvotes: 3
Views: 8685
Answers (1)
Reputation: 14011
I have actually implemented this method for OpenCV, but unfortunately it is only in the SVN trunk at this stage.
Here is a sample using the new method.
If you would like to reference my implementation of it, you can find that here.
Also, here is the test case for an additional example of using it.
If you would like to use the trunk you can get a hold of it doing this:
svn co https://code.ros.org/svn/opencv/trunk/opencv opencv-trunk
Here is the CMake build guide for Linux. Here is the build guide for Windows.
EDIT : Got a patch for you :)
I found a few bugs in my code as well, so they should be corrected now. It also should support 1D phase-correlation now. I also found an issue with the cv::sqrt() function causing some -nan's to show up even though std::sqrt() did not. I'm guessing it's either a bug with OpenCV, or just an accuracy issue. Haven't dug into it enough to find out though.
Unfortunately, you can't just use svn update to get my latest changes because I have to wait for the OpenCV developers to apply this patch. So, instead of you having to wait on that here is a patch you can apply to the $(OPENCV_SRC)/modules/imgproc/src/phasecorr.cpp file. Save this file as something like phasecorr.patch, and place it in the root of the OpenCV source directory. Here is a short TortoiseSVN guide for creating/applying patches from/to source trees.
Index: modules/imgproc/src/phasecorr.cpp
===================================================================
--- modules/imgproc/src/phasecorr.cpp (revision 6971)
+++ modules/imgproc/src/phasecorr.cpp (working copy)
@@ -83,8 +83,8 @@
for( j = 1; j <= rows - 2; j += 2 )
{
- dataDst[j*stepDst] = (float)((double)dataSrc[j*stepSrc]*dataSrc[j*stepSrc] +
- (double)dataSrc[(j+1)*stepSrc]*dataSrc[(j+1)*stepSrc]);
+ dataDst[j*stepDst] = (float)std::sqrt((double)dataSrc[j*stepSrc]*dataSrc[j*stepSrc] +
+ (double)dataSrc[(j+1)*stepSrc]*dataSrc[(j+1)*stepSrc]);
}
if( k == 1 )
@@ -103,7 +103,7 @@
for( j = j0; j < j1; j += 2 )
{
- dataDst[j] = (float)((double)dataSrc[j]*dataSrc[j] + (double)dataSrc[j+1]*dataSrc[j+1]);
+ dataDst[j] = (float)std::sqrt((double)dataSrc[j]*dataSrc[j] + (double)dataSrc[j+1]*dataSrc[j+1]);
}
}
}
@@ -127,8 +127,8 @@
for( j = 1; j <= rows - 2; j += 2 )
{
- dataDst[j*stepDst] = dataSrc[j*stepSrc]*dataSrc[j*stepSrc] +
- dataSrc[(j+1)*stepSrc]*dataSrc[(j+1)*stepSrc];
+ dataDst[j*stepDst] = std::sqrt(dataSrc[j*stepSrc]*dataSrc[j*stepSrc] +
+ dataSrc[(j+1)*stepSrc]*dataSrc[(j+1)*stepSrc]);
}
if( k == 1 )
@@ -147,13 +147,10 @@
for( j = j0; j < j1; j += 2 )
{
- dataDst[j] = dataSrc[j]*dataSrc[j] + dataSrc[j+1]*dataSrc[j+1];
+ dataDst[j] = std::sqrt(dataSrc[j]*dataSrc[j] + dataSrc[j+1]*dataSrc[j+1]);
}
}
}
-
- // do batch sqrt to use SSE optimizations...
- cv::sqrt(dst, dst);
}
static void divSpectrums( InputArray _srcA, InputArray _srcB, OutputArray _dst, int flags, bool conjB)
@@ -196,9 +193,9 @@
{
if( k == 1 )
dataA += cols - 1, dataB += cols - 1, dataC += cols - 1;
- dataC[0] = dataA[0] / dataB[0];
+ dataC[0] = dataA[0] / (dataB[0] + eps);
if( rows % 2 == 0 )
- dataC[(rows-1)*stepC] = dataA[(rows-1)*stepA] / dataB[(rows-1)*stepB];
+ dataC[(rows-1)*stepC] = dataA[(rows-1)*stepA] / (dataB[(rows-1)*stepB] + eps);
if( !conjB )
for( j = 1; j <= rows - 2; j += 2 )
{
@@ -239,9 +236,9 @@
{
if( is_1d && cn == 1 )
{
- dataC[0] = dataA[0] / dataB[0];
+ dataC[0] = dataA[0] / (dataB[0] + eps);
if( cols % 2 == 0 )
- dataC[j1] = dataA[j1] / dataB[j1];
+ dataC[j1] = dataA[j1] / (dataB[j1] + eps);
}
if( !conjB )
@@ -281,9 +278,9 @@
{
if( k == 1 )
dataA += cols - 1, dataB += cols - 1, dataC += cols - 1;
- dataC[0] = dataA[0] / dataB[0];
+ dataC[0] = dataA[0] / (dataB[0] + eps);
if( rows % 2 == 0 )
- dataC[(rows-1)*stepC] = dataA[(rows-1)*stepA] / dataB[(rows-1)*stepB];
+ dataC[(rows-1)*stepC] = dataA[(rows-1)*stepA] / (dataB[(rows-1)*stepB] + eps);
if( !conjB )
for( j = 1; j <= rows - 2; j += 2 )
{
@@ -323,9 +320,9 @@
{
if( is_1d && cn == 1 )
{
- dataC[0] = dataA[0] / dataB[0];
+ dataC[0] = dataA[0] / (dataB[0] + eps);
if( cols % 2 == 0 )
- dataC[j1] = dataA[j1] / dataB[j1];
+ dataC[j1] = dataA[j1] / (dataB[j1] + eps);
}
if( !conjB )
@@ -354,31 +351,57 @@
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;
-
- for(size_t i = 0; i < planes.size(); i++)
+
+ bool is_1d = xMid == 0 || yMid == 0;
+
+ if(is_1d)
{
- // 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);
+ 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);
}
@@ -548,38 +571,67 @@
int rows = dst.rows;
int cols = dst.cols;
+ bool is_1d = rows == 1 || cols == 1;
+
+ if(is_1d)
+ {
+ cols = cols + rows - 1;
+ }
+
if(dst.depth() == CV_32F)
{
float* dstData = (float*)dst.data;
- for(int i = 0; i < rows; i++)
+ if(is_1d)
{
- double wr = 0.5 * (1.0f - cos(2.0f * CV_PI * (double)i / (double)(rows - 1)));
- for(int j = 0; j < cols; j++)
+ for(int i = 0; i < cols; i++)
{
- double wc = 0.5 * (1.0f - cos(2.0f * CV_PI * (double)j / (double)(cols - 1)));
- dstData[i*cols + j] = (float)(wr * wc);
+ double w = 0.5 * (1.0f - cos(2.0f * CV_PI * (double)i / (double)(cols - 1)));
+ dstData[i] = (float)w;
}
}
+ else
+ {
+ for(int i = 0; i < rows; i++)
+ {
+ double wr = 0.5 * (1.0f - cos(2.0f * CV_PI * (double)i / (double)(rows - 1)));
+ for(int j = 0; j < cols; j++)
+ {
+ double wc = 0.5 * (1.0f - cos(2.0f * CV_PI * (double)j / (double)(cols - 1)));
+ dstData[i*cols + j] = (float)(wr * wc);
+ }
+ }
- // perform batch sqrt for SSE performance gains
- cv::sqrt(dst, dst);
+ // perform batch sqrt for SSE performance gains
+ cv::sqrt(dst, dst);
+ }
}
else
{
double* dstData = (double*)dst.data;
- for(int i = 0; i < rows; i++)
+ if(is_1d)
{
- double wr = 0.5 * (1.0 - cos(2.0 * CV_PI * (double)i / (double)(rows - 1)));
- for(int j = 0; j < cols; j++)
+ for(int i = 0; i < cols; i++)
{
- double wc = 0.5 * (1.0 - cos(2.0 * CV_PI * (double)j / (double)(cols - 1)));
- dstData[i*cols + j] = wr * wc;
+ double w = 0.5 * (1.0f - cos(2.0f * CV_PI * (double)i / (double)(cols - 1)));
+ dstData[i] = w;
}
}
+ else
+ {
+ for(int i = 0; i < rows; i++)
+ {
+ double wr = 0.5 * (1.0 - cos(2.0 * CV_PI * (double)i / (double)(rows - 1)));
+ for(int j = 0; j < cols; j++)
+ {
+ double wc = 0.5 * (1.0 - cos(2.0 * CV_PI * (double)j / (double)(cols - 1)));
+ dstData[i*cols + j] = wr * wc;
+ }
+ }
- // perform batch sqrt for SSE performance gains
- cv::sqrt(dst, dst);
+ // perform batch sqrt for SSE performance gains
+ cv::sqrt(dst, dst);
+ }
}
}
Finally, here is a code sample of using the 1D phase correlation I was using.
int main(int argc, char* argv[])
{
Mat firstArray = Mat::zeros(Size(360, 1), CV_64F);
Mat secondArray = Mat::zeros(Size(360, 1), CV_64F);
for(int i = 0; i < firstArray.cols; i++)
{
if(i < 8)
{
firstArray.at<double>(0, i) = 1;
}
if(i < 6)
{
secondArray.at<double>(0, i) = 1;
}
}
Mat hann;
createHanningWindow(hann, firstArray.size(), CV_64F);
Point2d shift = phaseCorrelate(firstArray, secondArray, hann);
std::cout<< "shift: " << shift.x << ";" << shift.y << std::endl;
return 0;
}
You should see an output of (-2, 0.5). The y value will always be 0.5 in the 1D case because that would be midway through the only row.
Hope that is helpful to you!
Upvotes: 8
Related Questions
- Is the OpenCV phase() function working correctly in this edge case?
- What is the difference between phase correlation and template matching in OpenCV?
- How to find phase spectrum of an image in opencv python
- Get rotational shift using phase correlation and log polar transform
- Interpreting results of the cv2 method: phaseCorrelate? | python opencv
- Calculate phase angle in Optical flow
- How to do inverse DFT using magnitude and phase of a image in opencv?
- phase correlation using FFTW
- phaseCorrelate function in OpenCV Python
- phase correlation 2D (Stitch 2D in ImageJ) for image stitching