krz*_*ych 3 opencv signal-processing image-processing
我想计算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;
}
我究竟做错了什么？
我实际上已经为OpenCV实现了这个方法,但不幸的是它只在这个阶段的SVN主干中.
以下是使用新方法的示例.
如果您想参考我的实现,可以在这里找到.
此外,这是使用它的另一个例子的测试用例.
如果您想使用主干,您可以抓住它来执行此操作:
svn co https://code.ros.org/svn/opencv/trunk/opencv opencv-trunk
这是Linux的CMake构建指南.这是Windows的构建指南.
编辑: 给你补丁:)
我在代码中也发现了一些错误,所以现在应该纠正它们.它现在也应该支持一维相位相关.我还发现cv::sqrt()函数的问题导致一些-nan出现,即使std::sqrt()没有.我猜它是OpenCV的一个bug,或者只是一个准确性问题.尽管如此,还是没有足够的挖掘它.
不幸的是,您不能仅仅使用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在1D情况下,该值始终为0.5,因为它将位于唯一行的中间.
希望对您有所帮助!
| 归档时间： |
|
| 查看次数： |
7119 次 |
| 最近记录： |