Open Cv相位相關

Question

我想計算OpenCV中函數的相位相關性。我發現http://nashruddin.com/phase-correlation-function-in-opencv.html，但有一些其他庫，我想用OpenCV來做。

我的問題是不同的，因爲我想計算兩個數組的相位相關性，包括360個元素。 我試圖從文檔評估如何做到這一點，但我不知道我的方法是否好。 我的R矩陣未標準化，如何標準化它，是否有必要？ 如果你給我一些例子，我將不勝感激。 我的函數來完成這個任務：

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; 
} 

我知道，代碼是不清晰等優點，但它只是快速測試。我想知道方法是否正確。但也有每一個建議表示讚賞。

//編輯： 我使用mevatron代碼，並得到#QNAN也調試時發現，它沒有找到峯值的函數，值是（-1，-1）。 功能我使用的是新的OpenCV的相關函數：

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; 
} 

我在做什麼錯？

Answer 1

我實際上已經實現了OpenCV的這種方法，但不幸的是它只在SVN主幹中處於這個階段。

Here是一個使用新方法的示例。

如果你想引用我的實現它，你可以找到那here。

此外，here是另一個使用它的例子的測試案例。

如果你想使用主幹，你可以得到它保持這樣：

svn co https://code.ros.org/svn/opencv/trunk/opencv opencv-trunk 

Here是CMake的製作指南的Linux版本。 Here是Windows的編譯指南。

編輯： 有一個補丁，你:)

我發現在我的代碼中的一些錯誤，以及，所以他們現在應該糾正。它現在也應該支持一維相位關聯。我還發現一個cv::sqrt()函數的問題，即使std::sqrt()沒有。我猜這是OpenCV的一個bug，或者只是一個準確性問題。儘管如此，還是沒有深入研究。

不幸的是，您不能僅僅使用svn update來獲得我的最新更改，因爲我必須等待OpenCV開發人員應用此修補程序。所以，而不是你不得不等待，這是一個補丁，你可以申請$(OPENCV_SRC)/modules/imgproc/src/phasecorr.cpp文件。將此文件保存爲類似phasecorr.patch的文件，並將其放在OpenCV源目錄的根目錄下。 Here是一個簡短的TortoiseSVN指南，用於創建/應用源/目錄樹中的補丁。

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); 
+  } 
    } 
} 

最後，這是一個使用我正在使用的一維相位相關的代碼示例。

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; 
} 

您應該看到（-2,0.5）的輸出。在1D情況下，y的值將始終爲0.5，因爲這將在唯一行的中間。

希望對你有幫助！