代码
先给出代码，再详细解释一下过程：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 |
#include <opencv2/core/core.hpp> #include <opencv2/imgproc/imgproc.hpp> #include <opencv2/highgui/highgui.hpp> #include <iostream> using namespace cv; using namespace std; #define GRAY_THRESH 150 #define HOUGH_VOTE 100 //#define DEGREE 27 int main(int argc, char **argv) { //Read a single-channel image const char* filename = "imageText.jpg"; Mat srcImg = imread(filename, CV_LOAD_IMAGE_GRAYSCALE); if(srcImg.empty()) return -1; imshow("source", srcImg); Point center(srcImg.cols/2, srcImg.rows/2); #ifdef DEGREE //Rotate source image Mat rotMatS = getRotationMatrix2D(center, DEGREE, 1.0); warpAffine(srcImg, srcImg, rotMatS, srcImg.size(), 1, 0, Scalar(255,255,255)); imshow("RotatedSrc", srcImg); //imwrite("imageText_R.jpg",srcImg); #endif //Expand image to an optimal size, for faster processing speed //Set widths of borders in four directions //If borderType==BORDER_CONSTANT, fill the borders with (0,0,0) Mat padded; int opWidth = getOptimalDFTSize(srcImg.rows); int opHeight = getOptimalDFTSize(srcImg.cols); copyMakeBorder(srcImg, padded, 0, opWidth-srcImg.rows, 0, opHeight-srcImg.cols, BORDER_CONSTANT, Scalar::all(0)); Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)}; Mat comImg; //Merge into a double-channel image merge(planes,2,comImg); //Use the same image as input and output, //so that the results can fit in Mat well dft(comImg, comImg); //Compute the magnitude //planes[0]=Re(DFT(I)), planes[1]=Im(DFT(I)) //magnitude=sqrt(Re^2+Im^2) split(comImg, planes); magnitude(planes[0], planes[1], planes[0]); //Switch to logarithmic scale, for better visual results //M2=log(1+M1) Mat magMat = planes[0]; magMat += Scalar::all(1); log(magMat, magMat); //Crop the spectrum //Width and height of magMat should be even, so that they can be divided by 2 //-2 is 11111110 in binary system, operator & make sure width and height are always even magMat = magMat(Rect(0, 0, magMat.cols & -2, magMat.rows & -2)); //Rearrange the quadrants of Fourier image, //so that the origin is at the center of image, //and move the high frequency to the corners int cx = magMat.cols/2; int cy = magMat.rows/2; Mat q0(magMat, Rect(0, 0, cx, cy)); Mat q1(magMat, Rect(0, cy, cx, cy)); Mat q2(magMat, Rect(cx, cy, cx, cy)); Mat q3(magMat, Rect(cx, 0, cx, cy)); Mat tmp; q0.copyTo(tmp); q2.copyTo(q0); tmp.copyTo(q2); q1.copyTo(tmp); q3.copyTo(q1); tmp.copyTo(q3); //Normalize the magnitude to [0,1], then to[0,255] normalize(magMat, magMat, 0, 1, CV_MINMAX); Mat magImg(magMat.size(), CV_8UC1); magMat.convertTo(magImg,CV_8UC1,255,0); imshow("magnitude", magImg); //imwrite("imageText_mag.jpg",magImg); //Turn into binary image threshold(magImg,magImg,GRAY_THRESH,255,CV_THRESH_BINARY); imshow("mag_binary", magImg); //imwrite("imageText_bin.jpg",magImg); //Find lines with Hough Transformation vector<Vec2f> lines; float pi180 = (float)CV_PI/180; Mat linImg(magImg.size(),CV_8UC3); HoughLines(magImg,lines,1,pi180,HOUGH_VOTE,0,0); int numLines = lines.size(); for(int l=0; l<numLines; l++) { float rho = lines[l][0], theta = lines[l][1]; Point pt1, pt2; double a = cos(theta), b = sin(theta); double x0 = a*rho, y0 = b*rho; pt1.x = cvRound(x0 + 1000*(-b)); pt1.y = cvRound(y0 + 1000*(a)); pt2.x = cvRound(x0 - 1000*(-b)); pt2.y = cvRound(y0 - 1000*(a)); line(linImg,pt1,pt2,Scalar(255,0,0),3,8,0); } imshow("lines",linImg); //imwrite("imageText_line.jpg",linImg); if(lines.size() == 3){ cout << "found three angels:" << endl; cout << lines[0][1]*180/CV_PI << endl << lines[1][1]*180/CV_PI << endl << lines[2][1]*180/CV_PI << endl << endl; } //Find the proper angel from the three found angels float angel=0; float piThresh = (float)CV_PI/90; float pi2 = CV_PI/2; for(int l=0; l<numLines; l++) { float theta = lines[l][1]; if(abs(theta) < piThresh || abs(theta-pi2) < piThresh) continue; else{ angel = theta; break; } } //Calculate the rotation angel //The image has to be square, //so that the rotation angel can be calculate right angel = angel<pi2 ? angel : angel-CV_PI; if(angel != pi2){ float angelT = srcImg.rows*tan(angel)/srcImg.cols; angel = atan(angelT); } float angelD = angel*180/(float)CV_PI; cout << "the rotation angel to be applied:" << endl << angelD << endl << endl; //Rotate the image to recover Mat rotMat = getRotationMatrix2D(center,angelD,1.0); Mat dstImg = Mat::ones(srcImg.size(),CV_8UC3); warpAffine(srcImg,dstImg,rotMat,srcImg.size(),1,0,Scalar(255,255,255)); imshow("result",dstImg); //imwrite("imageText_D.jpg",dstImg); waitKey(0); return 0; } |
过程
读取图片
|
1 2 3 |
Mat srcImg = imread(filename, CV_LOAD_IMAGE_GRAYSCALE); if(srcImg.empty()) return -1; |
srcImg.empty()用来判断是否成功读进图像，如果srcImg中没有数据，在后面的步骤会产生内存错误。
由于处理的是文本，彩色信息不会提供额外帮助，所以要用CV_LOAD_IMAGE_GRAYSCALE表明以灰度形式读进图像。
假定读取的图像如下：
旋转原图像(可选)
|
1 2 3 4 5 6 7 8 9 |
Point center(srcImg.cols/2, srcImg.rows/2); #ifdef DEGREE //Rotate source image Mat rotMatS = getRotationMatrix2D(center, DEGREE, 1.0); warpAffine(srcImg, srcImg, rotMatS, srcImg.size(), 1, 0, Scalar(255,255,255)); imshow("RotatedSrc", srcImg); //imwrite("H:\\imageText_02_R.jpg",srcImg); #endif |
如果手头没有这样的倾斜图像，可以选择一张正放的文本图像，再把第12行#define DEGREE那行前的注释符号去掉。然后这部分代码就会把所给的图像旋转你规定的角度，再交给后面处理。
图像延扩
|
1 2 3 4 |
Mat padded; int opWidth = getOptimalDFTSize(srcImg.rows); int opHeight = getOptimalDFTSize(srcImg.cols); copyMakeBorder(srcImg, padded, 0, opWidth-srcImg.rows, 0, opHeight-srcImg.cols, BORDER_CONSTANT, Scalar::all(0)); |
OpenCV中的DFT采用的是快速算法，这种算法要求图像的尺寸是2、3和5的倍数时处理速度最快。所以需要用getOptimalDFTSize()找到最适合的尺寸，然后用copyMakeBorder()填充多余的部分。这里是让原图像和扩大的图像左上角对齐。填充的颜色如果是纯色对变换结果的影响不会很大，后面寻找倾斜线的过程又会完全忽略这一点影响。
DFT
|
1 2 3 4 5 6 |
Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)}; Mat comImg; merge(planes,2,comImg); dft(comImg, comImg); |
DFT要分别计算实部和虚部，把要处理的图像作为输入的实部、一个全零的图像作为输入的虚部。dft()输入和输出应该分别为单张图像，所以要先用merge()把实虚部图像合并，分别处于图像comImg的两个通道内。计算得到的实虚部仍然保存在comImg的两个通道内。
获得DFT图像
|
1 2 3 4 5 6 |
split(comImg, planes); magnitude(planes[0], planes[1], planes[0]); Mat magMat = planes[0]; magMat += Scalar::all(1); log(magMat, magMat); |
一般都会用幅度图像来表示图像傅里叶的变换结果（傅里叶谱）。
幅度的计算公式：magnitude = sqrt(Re(DFT)^2 + Im(DFT)^2)。
由于幅度的变化范围很大，而一般图像亮度范围只有[0,255]，容易造成一大片漆黑，只有几个点很亮。所以要用log函数把数值的范围缩小。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
magMat = magMat(Rect(0, 0, magMat.cols & -2, magMat.rows & -2)); int cx = magMat.cols/2; int cy = magMat.rows/2; Mat q0(magMat, Rect(0, 0, cx, cy)); Mat q1(magMat, Rect(0, cy, cx, cy)); Mat q2(magMat, Rect(cx, cy, cx, cy)); Mat q3(magMat, Rect(cx, 0, cx, cy)); Mat tmp; q0.copyTo(tmp); q2.copyTo(q0); tmp.copyTo(q2); q1.copyTo(tmp); q3.copyTo(q1); tmp.copyTo(q3); normalize(magMat, magMat, 0, 1, CV_MINMAX); Mat magImg(magMat.size(), CV_8UC1); magMat.convertTo(magImg,CV_8UC1,255,0); |
dft()直接获得的结果中，低频部分位于四角，高频部分位于中间。习惯上会把图像做四等份，互相对调，使低频部分位于图像中心，也就是让频域原点位于中心。
虽然用log()缩小了数据范围，但仍然不能保证数值都落在[0,255]之内，所以要先用normalize()规范化到[0,1]内，再用convertTo()把小数映射到[0,255]内的整数。结果保存在一幅单通道图像内：
Hough直线检测
从傅里叶谱可以明显地看到一条过中心点的倾斜直线。要想求出这个倾斜角，首先要在图像上找出这条直线。
一个很方便的方法是采用霍夫（Hough）变换检测直线。
|
1 |
threshold(magImg,magImg,GRAY_THRESH,255,CV_THRESH_BINARY); |
Hough变换要求输入图像是二值的，所以要用threshold()把图像二值化。
二值化的一种结果：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
vector<Vec2f> lines; float pi180 = (float)CV_PI/180; Mat linImg(magImg.size(),CV_8UC3); HoughLines(magImg,lines,1,pi180,HOUGH_VOTE,0,0); int numLines = lines.size(); for(int l=0; l<numLines; l++) { float rho = lines[l][0], theta = lines[l][1]; Point pt1, pt2; double a = cos(theta), b = sin(theta); double x0 = a*rho, y0 = b*rho; pt1.x = cvRound(x0 + 1000*(-b)); pt1.y = cvRound(y0 + 1000*(a)); pt2.x = cvRound(x0 - 1000*(-b)); pt2.y = cvRound(y0 - 1000*(a)); line(linImg,pt1,pt2,Scalar(255,0,0),3,8,0); } |
这一部分用HoughLines()检测图像中可能存在的直线，并把直线参数保存在向量组lines中，然后绘制出找到的直线。
两个参数GRAY_THRESH和HOUGH_VOTE需要手动指定，不同的图像需要设置不同的参数，同一段文本旋转不同的角度也需要不同的参数。GRAY_THRESH越大，二值化的阈值就越高；HOUGH_VOTE越大，霍夫检测的投票数就越高（需要更多的共线点来确定一条直线）。说白了，如果发现二值化图像中直线附近有很多散点，就要适当提高GRAY_THRESH；如果发现从二值图像的一条直线上检测到了几条角度相差很小的直线，就需要适当提高HOUGH_VOTE。我们希望得到的结果时刚好检测到三条直线（有时只能检测到一条直线，后面会给出一个例子）。
检测到的直线：
计算倾斜角
上面得到了三个角度，一个是0度，一个是90度，另一个就是我们所需要的倾斜角。要把这个角找出来，而且要考虑误差。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
float angel=0; float piThresh = (float)CV_PI/90; float pi2 = CV_PI/2; for(int l=0; l<numLines; l++) { float theta = lines[l][1]; if(abs(theta) < piThresh || abs(theta-pi2) < piThresh) continue; else{ angel = theta; break; } } angel = angel<pi2 ? angel : angel-CV_PI; if(angel != pi2){ float angelT = srcImg.rows*tan(angel)/srcImg.cols; angel = atan(angelT); } float angelD = angel*180/(float)CV_PI; |
由于DFT的特点，只有输入图像是正方形时，检测到的角才是文本真正旋转的角度。但我们的输入图像不一定是正方形的，所以要根据图像的长宽比改变这个角度。
还有一个需要注意的细节，虽然HoughLines()输出的倾斜角在[0,180)之间，但在[0,90]和(90,180)之间这个角的含义是不同的。请看图示：
当倾斜角大于90度时，(180-倾斜角)才是直线相对竖直方向的偏离角度。在OpenCV中，逆时针旋转，角度为正。要把图像转回去，这个角度就变成了(倾斜角-180)。
校正图像
最后一步，当然是把图像转回去~
|
1 2 3 |
Mat rotMat = getRotationMatrix2D(center,angelD,1.0); Mat dstImg = Mat::ones(srcImg.size(),CV_8UC3); warpAffine(srcImg,dstImg,rotMat,srcImg.size(),1,0,Scalar(255,255,255)); |
先用getRotationMatrix2D()获得一个2*3的仿射变换矩阵，再把这个矩阵输入warpAffine()，做一个单纯旋转的仿射变换。warpAffine()的最后一个参数Scalar(255,255,255)是把由于旋转产生的空白用白色填充。
校正的结果：
一个检测单条直线的例子
原始图像：
傅里叶谱：
只有一条明显的直线。还好仅有的这条直线正是我们所需要的。
检测直线：
校正结果：
对中文的效果
我们来试试看这段程序对中文的校正效果。
输入图像：
傅里叶谱：
可以发现有许多条平行的亮线，其中过频域原点的那条长度最长，最容易检测出来。
检测直线：
校正结果：
虽然中文和英文在文字上有很大的不同，但字母（或者文字）的高度比较一致，使得行与行之间的分隔很明显。所以它们的频域特征是相似的。
对其他语言文字的效果
我从IMDB.com摘取影片《教父》的英文介绍，然后用谷歌翻译成其他文字进行测试。
阿拉伯语
一枚反例
老挝语：
傅里叶谱：
一种二值化的结果：
直线检测：
这种文字的很多字母的上下方多了很多“笔画”（我不知道该怎么称呼那些小曲线），让行与行之间的分离变得不明显，使得频域特征变得不明显。
虽然用肉眼可以看出傅里叶谱中存在一条倾斜的直线，但它的亮度太低，二值化过程很难排除噪声，导致直线检测会首先检出噪声产生的直线。这也是我的程序目前受限之处。需要增加一个过滤散点噪声的步骤以增加程序的适用范围。
参考：Discrete Fourier Transform — OpenCV 2.4.7.0 documentation
代码还可以在这里下载：https://github.com/johnhany/textRotCorrect
2014.1.3更新：
由于文章内的图片右下角存在水印，若直接使用文章内的图片进行处理会使频域原点附近增加一团亮点，妨碍直线的检出。而且为了节省空间，图片是经过缩小的，使得字母的边缘变得模糊，频域特征也减弱。为此我提供了十幅没有水印的图片，供想要亲手实验的朋友使用。下载链接
博主，你好，我是用lena图，按你说的设置了旋转角度为27，但是通过Hough变换检测直线计算出的角度是63,校正回去就表现为更远离正图了。hough变换那里我检测出的有四条线，我尝试把HOUGH_VOTE增大，就完全检测不出来直线了
很想试验一下opencv中离散傅里叶变换进行图像倾斜校正，正苦于找不到倾斜的图像，没想到博主还提供图片，在此先感谢博主，我先去试试，有问题再来请教！
博主，好像这个是基于大段文字的检测哦，对于就只有两个字符，比如：4m，这样的短字符检测有什么方法改进的么？不然识别不了啊!
对的，基于频谱的方法对一两个字符是无效的，因为无论从哪个方向上看亮度的变化都是有限的。只有一个字符且已知字典的话，可能需要用带旋转的模板匹配等传统的方法来做。两个以上字符时，完全可以找出字符区域的长宽结合字符匹配来校正。
博主有90度和180度图片旋转的解决方案吗？
有java版本的实现么？
你好，博主，这个地方不理解，可以解释一下吗？
由于DFT的特点，只有输入图像是正方形时，检测到的角才是文本真正旋转的角度。但我们的输入图像不一定是正方形的，所以要根据图像的长宽比改变这个角度。
因为DFT之前的原图像在x y方向上表示空间坐标，DFT是经过x y方向上的傅里叶变换来统计像素在这两个方向上不同频率的分布情况，所以DFT得到的图像在x y方向上不再表示空间上的长度，而是频率。opencv的DFT函数设计为输出（频谱）图像与输入（像素）图像尺寸一致。如果输入图像是正方形，那么输出图像在x y方向上的（频率）坐标间隔是相同的；如果输入图像不是正方形，就会造成输出图像在x或y方向上的缩放，即两轴的坐标间隔不一致。正方形被拉伸为长方形，其对角线的斜率也势必产生变化。所以在检测直线之后，还要根据图像的原始尺寸调整这个对角线的倾角。
楼主，你好。我尝试用你的代码做倾斜矫正，但是对于同一张不同倾斜角度的图片，需要调整灰度化参数值和霍夫变换投票率才能达到比较好的效果。想请教一下有没有什么好的自适应方法。
你好，我想问问，为什么频谱图象进行中心转移之后，计算得到的倾角还是原来的那个角呢？不会变嘛？
Android下使用傅里叶变换检测出linenum=0是怎么回事，谢谢解惑！
您好，请问一下您使用的这个方法是有没有发过相应的文章，因为我的实验里用到了您的方法，论文的时候需要写这个方法的来源，所以想问一下您这个方法来自您自己发过的文章里的还是其他书本或期刊上的。谢谢您。
我的方法并不是完全原创，参考的资料也只有文章末尾的那个官方文档页面。那个页面的最后一节只是提到水平与倾斜的文本频谱图像是有差异的，也没有提及相关的文献，我也是自己用代码实验了一下，所以也没有发过相应的文章。某些杂志是允许把网址作为参考文献的，可以按照格式要求直接把这篇文章的链接写上。
谢谢。
您好，看了您的方法后，我用C语言试用了一下，但是不知道为什么在HoughLines后，呈现出来的图象是全黑的。为了方便，我贴出这部分的代码，不知道哪里出了问题，想向您请教一下。谢谢。
CvMemStorage *storage = cvCreateMemStorage(0);
CvSeq* lines = NULL;
float pi180 = (float)CV_PI/180;
lines = cvHoughLines2(Image,storage,CV_HOUGH_STANDARD,1,pi180,HOUGH_VOTE,0,0); //Image是之前的二值化图象
IplImage *lineImg;
lineImg = cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1);
int numLines = lines->total;
int i=0;
for(i=0;i<numLines;i++)
{
float *line = (float*) cvGetSeqElem(lines,i);
float rho = line[0];
float theta = line[1];
CvPoint pt1,pt2;
double a = cos(theta), b = sin(theta);
double x0 = a*rho, y0 = b*rho;
pt1.x = cvRound(x0+1000*(-b));
pt1.y = cvRound(y0+1000*(a));
pt2.x = cvRound(x0-1000*(-b));
pt2.y = cvRound(y0-1000*(a));
cvLine(lineImg,pt1,pt2,CV_RGB(255,0,0),3,8,0);
}
cvNamedWindow("line",0);
cvShowImage("line",lineImg);
cvWaitKey(10000);
抱歉，最近事情比较多，现在问题解决了吗？
已经解决了，抱歉，最近都没有看信息。
您好，请问houghline之后呈现图像是黑的这个问题要怎么解决？我是一枚编程小白，最近在研究畸变矫正的相关知识，希望得到您的帮助
你指的是直线图像linImg全黑还是原图像magImg经过HoughLines()函数之后发生了变化？如果是前者的话，有可能由于参数的设置没有正确检出直线（numLines=0）；后者是由于HoughLines()函数在计算的过程中改变了输入图像的值，所以最好先复制一下输入图像再输入该函数。
您好，我是指直线图像lineImg全黑，在此之前显示的图像是正确的，关闭图像后会出现如图所示的错误（在这次实验中所用的代码&图像都是您在上文中使用的内容，我没有更改QAQ由于在HOUGH线获取环节出错，之后的输出：“find three angles:”也没有显示出来。）
抱歉，图片显示错误，下面是调试时出错的代码段
// MEMBER FUNCTIONS FOR _Container_base12
inline void _Container_base12::_Orphan_all()
{ // orphan all iterators
#if _ITERATOR_DEBUG_LEVEL == 2
if (_Myproxy != 0)
{ // proxy allocated, drain it
_Lockit _Lock(_LOCK_DEBUG);
for (_Iterator_base12 **_Pnext = &_Myproxy->_Myfirstiter;
*_Pnext != 0; *_Pnext = (*_Pnext)->_Mynextiter)
(*_Pnext)->_Myproxy = 0;（在这一行出现了错误）
_Myproxy->_Myfirstiter = 0;
}
#endif /* _ITERATOR_DEBUG_LEVEL == 2 */
}
其中局部变量的名称 值 类型 如下：
_Pnext 0x0000000000473918 {0x3fc90fdb43700000 {_Myproxy=??? _Mynextiter=??? }} std::_Iterator_base12 * *
_Lock {_Locktype=3 } std::_Lockit
this 0x000000000026f278 {_Myproxy=0x0000000000473910 {_Mycont=0x0000000043f00000 {_Myproxy=??? } _Myfirstiter=…} } std::_Container_base12 *
麻烦您了
很好用的代码，十分感谢！
多谢支持！
作者您好，我想做的是一幅图像经过傅里叶变换之后，分别得到实部和虚部的数据，看了您的文章，深受启发，但是遇到个问题，下面三步算是对原有图像填充的尺寸，但是填充值后傅里叶变换得到的虚部和实部尺寸也已改变，怎么让实部虚部数据都还原到填充前的尺寸呢？非常感谢！
int opWidth = getOptimalDFTSize(srcImg.rows);
int opHeight = getOptimalDFTSize(srcImg.cols);
copyMakeBorder(srcImg, padded, 0, opWidth-srcImg.rows, 0, opHeight-srcImg.cols, BORDER_CONST
原有图像的纯色填充对频域的影响应该是可以忽略的。况且这三个函数改变的是图像的空间域，经过FFT之后得到的频域图像，其尺寸也失去了与原来空间域尺寸比较的意义。虽然OpenCV库函数默认产生的频谱图像本身的大小与输入图像相同，但实际是对频域横纵坐标缩放的结果，其本质上仍然是沿坐标轴对称的正方形区域。如果直接对频谱图像切割，会造成原始图像频率上信息的丢失（高频或是低频，取决于是否将四象限对调）。
关于频域对应空间域的关系，可以参考这个问题的第一个答案，个人觉得解释得很精湛：http://dsp.stackexchange.com/questions/1637/what-does-frequency-domain-denote-in-case-of-images
感谢您的详细解答，可能我太笨啦，有个问题还是想不太明白，希望您解惑。问题描述：其实我是在把MATLAB代码变成c++时遇到的提取图像经过fft变换后的实部和虚部数据这个问题，下面是MATLAB代码： … 阅读更多 »
为方便您的理解我贴出函数PCE的MATLAB代码：
shift_range = [0,0];
Cinrange = rel(end-shift_range(1):end,end-shift_range(2):end);
[max_cc, imax] = max(Cinrange(:));
[ypeak, xpeak] = ind2sub(size(Cinrange),imax(1));
peakheight = Cinrange(ypeak,xpeak);
C_without_peak = RemoveNeighborhood(C,[ypeak, xpeak],squaresize);//RemoveNeighborhood是一个函数，下面有代码
correl = rel(end,end); clear rel
PCE_energy = mean(C_without_peak.*C_without_peak);
detection = peakheight.^2/PCE_energy * sign(peakheight);
//RemoveNeighborhood函数
function Y = RemoveNeighborhood(X,x,ssize)
% Remove a 2-D neighborhood around x=[x1,x2] from matrix X and output a 1-D vector Y
% ssize square neighborhood has size (ssize x ssize) square
[M,N] = size(X);
radius = (ssize-1)/2;
X = circshift(X,[radius-x(1)+1,radius-x(2)+1]);
Y = X(ssize+1:end,1:ssize);
Y = Y(:);
Y = [Y;X(M*ssize+1:end)'];
非常感谢！
在fft前不进行速度优化，直接对图像进行操作即可
是把高频部分移到中心位置吧？
象限换位的目的就是把频域原点移到图像中心，所以移位后低频是位于中心区域的。频谱图像中心的小亮点F(0,0)表示图像的平均灰度，也能说明这一点。
感兴趣的话，可以试着把一幅图像作高斯模糊，再比较一下前后的频谱图像，坐标轴附近的点亮度是会增加的。
试了一下 英文很好用，中文就不行了，葡萄牙语也不行，西班牙语还能凑合着用。总之，非常感谢博主分享。另有什么改善的方法么？
多谢支持！就方法上的改进可以用调整参数，过滤散点噪声等手段，但从理论上我目前还没有更好的改进思路。
非常感谢如此详细的文章！
多谢支持！
能不能判断出一个图片中的三角形倾斜的角度？
抱歉回复的很晚。本文的方法是不能用来计算单个三角形倾斜的角度的，只能考虑其他方法