Фазовая корреляция с 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.

Здесь - образец с использованием нового метода.

Если вы хотите сослаться на мою реализацию, вы можете найти, что здесь .

Также здесь - это тестовый пример для дополнительного примера его использования.

Если вы хотите использовать сундук, вы можете получить его, выполнив следующее:

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

Здесь - руководство по сборке CMake для Linux. Здесь - руководство по сборке для Windows.

РЕДАКТИРОВАТЬ: Получил для вас патч :)

Я также нашел несколько ошибок в своем коде, так что теперь они должны быть исправлены. Теперь он также должен поддерживать одномерную фазовую корреляцию. Я также обнаружил проблему с функцией cv::sqrt(), из-за которой некоторые -nan отображались, хотя std::sqrt() этого не делал. Я предполагаю, что это либо ошибка в OpenCV, либо просто проблема точности. Хотя не вникнул в это достаточно, чтобы выяснить это.

К сожалению, вы не можете просто использовать svn update, чтобы получить мои последние изменения, потому что мне нужно подождать, пока разработчики OpenCV установят этот патч. Таким образом, вместо того, чтобы ждать, вот патч, который вы можете применить к файлу $(OPENCV_SRC)/modules/imgproc/src/phasecorr.cpp. Сохраните этот файл как phasecorr.patch и поместите его в корень исходного каталога OpenCV. Здесь - это краткое руководство 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);
+        }
     }
 }

Наконец, вот пример кода использования 1D фазовой корреляции, которую я использовал.

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). Значение y всегда будет 0,5 в одномерном случае, потому что это будет на полпути через единственный ряд.

Надеюсь, что это полезно для вас!