计算机视觉与模式识别（1）—— A4纸边缘提取

转载自：http://blog.csdn.net/simba_scorpio/article/details/50974351

最近跟着老师学习了计算机视觉领域中，一个教科书般的知识点 —— 图形的边缘提取。随后自己写了一个A4纸边缘提取的小程序（200行），觉得蛮有意思的，测试的效果也不错，于是心血来潮，打算把自己整个设计的感想写下来。代码可能不是很规范，因为精力都侧重在理论知识上了，本章的全部难点几乎都在图形处理的理论知识+_+！

        那么首先说说这次的目标是什么：

       输入：普通的A4打印纸，上面可能有手写或打印字体，而且拍照时角度不正。

          

       输出：

        1. 图像的边缘

        2. 计算A4纸边缘的各直线方程

        3. 提取A4纸的4个角点

       我的环境：

        Windows10、C++11、Visual Studio2013

        有很多图形处理的相关库文件，这里我使用了老师推荐的CImg库，可以在官网下载：http://www.cimg.eu/，库的使用十分简单，只需将解压文件中的CImg.h头文件加入工程中即可。更多的有关库的使用这里不过多描述，在后面代码讲解中，我会对关键的函数进行解释。事实上，CImg作为专业的图形库，自身已经实现了类似边缘提取的函数，不过我们这里不会使用，因为这样不就没有意义了吗╮(╯▽╰)╭。

        话不多说，从输入来看要提取出A4纸的边缘绝非易事。从我们最直观的感觉来看，首先我们要区分颜色，然后把白色的那块圈起来就得到我们的边缘了。说实话，我一开始也是这么想的，典型的小白思想。注意上面的第一张图，A4纸中会有手写体，所以并不是全白的，而且由于拍照时的光线原因，纸上还有阴影。除此以外，图像中还存在其他物体，比如桌子、地板、腿等。第二张图中还被另一张A4纸抢镜，我们不能把它也算进来。

        卧cao( ‵o′)，那还怎么做！

       

        莫慌！但不用抱紧我。现在把我的姿势告诉你：

        1. 把彩色图像灰度化。

        2. 高斯模糊一下。

        3. 计算灰度图的梯度。

        4. 把梯度高的坐标映射到Hough图坐标上。（Hough转换，一种Voting算法，关键！）

        5. 提取Hough图峰值

        6. 计算直线方程和角点

        经过这一系列的处理，我们就可以得到这样的结果：

        

        效果还不错吧，没有出现上述任何一项问题，A4纸4条边缘和角点的位置也十分贴近真实情况。

        我们一步步来。

一、彩色图像灰度化处理。

        首先读入BMP格式的原图片：

[cpp]  view plain  copy

 

1. #include "CImg.h"  
2. using namespace cimg_library;  
3.   
4. #define FILE "A4.bmp"  
5.   
6. int main() {  
7.     CImg<double> image(FILE);     // 原始图  
8.     CImg<double> grayImg(image);  // 灰度图  
9.   
10.     image.display();  
11. }  

        CImg类是CImg.h中的主要图形类，存储了图片的所有数据，包括长宽高，色彩通道和像素值等。它的初始化很简单，可以是文件路径或者是自身对象，也可以通过指定长宽高创建一张默认图片。CImg只是默认支持了BMP格式的图像，如果需要读取JPG等其它格式的图片，要配置GraphicsMagick或者ImageMagick，这里不多讲。

        接下来是遍历image的所有像素，将它们变为灰度值。灰度值的计算公式：

        Gray = (R * 299 + G * 587 + B * 114 + 500) / 1000;

[cpp]  view plain  copy

 

1. // Gray scale;  
2. cimg_forXY(grayImg, x, y) {  
3.     double R = grayImg(x, y, 0, 0);  
4.     double G = grayImg(x, y, 0, 1);  
5.     double B = grayImg(x, y, 0, 2);  
6.     double Gray = (R * 299 + G * 587 + B * 114 + 500) / 1000;  
7.     grayImg(x, y, 0, 0) = Gray;  
8.     grayImg(x, y, 0, 1) = Gray;  
9.     grayImg(x, y, 0, 2) = Gray;  
10. }  
11. grayImg.display();  

        cimg_forXY函数等价于：

[cpp]  view plain  copy

 

1. for (int y = 0; y < grayImg.height(); ++y) {  
2.     for (int x = 0; x < grayImg.width(); ++x) {  
3.         // do something  
4.     }  
5. }  

        由于每个坐标的像素有RGB三个色彩通道，所以需要将所有通道设置成灰度值。

        完成这部分代码，运行后应有如下结果：

        

二、高斯模糊。

        注意到灰度图中有很多噪声点，这不利于待会儿的梯度计算。梯度计算的目的就是把色差较大的像素找出来，通过这种方法找出边缘（注意我没说A4纸的边缘），我们不希望把噪声点也当做边缘。CImg里面有自带的模糊函数，我只需要调用一下就好了，定义一个BLUR宏变量，方便日后修改。

[cpp]  view plain  copy

 

1. #define BLUR 2    

[cpp]  view plain  copy

 

1. grayImg.blur(BLUR);     

三、计算灰度图的梯度。

        接下来就涉及专业的理论知识了，如何计算一个像素点的梯度呢？我要介绍一下Prewitt算子。Prewitt算子是一个3*3的矩阵，它长这样：

      横向模板       纵向模板：

        3*3矩阵的位置描述：9个位置依次是Ipp、Icp、Inp、Ipc、Icc、Inc、Ipn、Icn、Inn。

        （p表示previous，c表示current，n表示next，Ipp = img(x-1,y-1)，Icn=img(x,y+1)，...）

        计算一个点的梯度首先要找出它的邻域矩阵，此时被计算的点位于邻域矩阵的中心，如果是3*3邻域矩阵，则被计算的点位于Icc位置，其它8个点是被计算点的相邻点。通过Prewitt算子，我们可以计算出点的横向平均梯度是(Inc - Ipc)，纵向平均梯度是(Icp - Icn)。

        梯度的值 = sqrt ( (Inc-Ipc)^2 + (Icp-Icn)^2 )

[cpp]  view plain  copy

 

1. #define GRADLIMIT 20  

[cpp]  view plain  copy

 

1. CImg<double> gradnum(image.width(), image.height(), 1, 1, 0);  
2. // 定义3*3领域矩阵I  
3. CImg_3x3(I, double);  
4. // 遍历计算梯度值  
5. cimg_for3x3(grayImg, x, y, 0, 0, I, double) {  
6.     const double ix = Inc - Ipc;  
7.     const double iy = Icp - Icn;  
8.     double grad = std::sqrt(ix*ix + iy*iy);  
9.     // 梯度大于阈值才赋值  
10.     if (grad > GRADLIMIT) {  
11.         gradnum(x, y) = grad;  
12.     }  
13. }  
14. gradnum.display();  

        CImg_3x3是CImg提供的定义领域矩阵的函数，第一个参数是邻域矩阵的名字，第二个参数是邻域矩阵的数据类型。通过将邻域矩阵传入cimg_for3x3循环函数，数据就会自动保存在 Ipp~Inn 九个变量中。我们只打算找出梯度值较大的像素点，因此设置了一个阈值GRADLIMIT，声明为宏变量，只有比阈值大的点才能显示在新创建的图片gradnum里。注意gradnum初始化时只有一个色彩通道，这样就只用修改一个通道值。

        完成这部分代码，运行后应有如下结果：

四、梯度图映射到Hough图（Voting算法）

        细心的你可能已经发现，经过梯度处理，我们得到的图片已经得到了A4纸的大概边缘，并且消去了很多字体的边缘和其它物体的边缘。但是，仍有一些顽强的噪声点幸存下来了，它们可能是纸上的浓墨重彩，也可能是另一些边缘明显的物体，比如另一张A4纸的一角。因此我们不能利用最小二乘法或者线性回归方法来拟合A4纸的边缘，这些方法受干扰点的影响很大。

        伟大的Paul Hough提出了一种解决的办法，称为Hough变换，中文叫霍夫变换。它的主要思路是通过投票的方式，让原始数据点们自己选出一个可以代表他们的函数表达式，因此是一种投票算法。说得这么玄，不如举个例子来说明：

        假设我们现在有这么一些数据点：

        我们知道直线的表达式形式是：y = mx + b。通常我们把x和y看成自变量，把m和b看成是常数。

        现在我们做一次变换，把x和y看成常数，把m和b看成变量：b = -xm + y

        由于我们已经获得了所有的坐标 （x，y） ，每个坐标（x，y）代进表达式 b = -xm + y 都表示一条关于m和b的直线，因此可以在mb坐标系上画出一簇直线：

        这些直线会近乎相交于一点（注意是近乎不是完全），这一点的坐标（1，1）就是原始点们选举出来的值，这个值能最好地表示原始点的直线表达式，也就是：y = x + 1。

        为什么原始图中，右下方的噪声点没有影响到直线表达式呢？我们可以看到，噪声点在mb坐标系上表示的直线偏离了大众交点，它固然也会和其它直线相交，但是那些交点的值（每有一条直线经过就加一），远远比不上大众交点的值，因此被毫不留情地淘汰掉了，这不正体现了投票的精髓吗——少数服从多数。

        回到我们的程序，现实是，我们不可能定义一个无限长的m轴和b轴，因此需要把直角坐标系转换为极坐标系。这样就有：p = x*cos(a)+y*sin(a)

        极坐标系的好处是，角度a的定义域是 [0, 360] 度，距离p的定义域是 [0, sqrt(x^2+y^2)]。

        为了完成hough变换，首先创建一个新的CImg对象hough，用来记录点的值：

[cpp]  view plain  copy

 

1. int maxDistance = distance(image.width(), image.height());  
2. CImg<double> hough(360, maxDistance, 1, 1, 0);  

        distance的函数定义如下：

[cpp]  view plain  copy

 

1. // 计算两点间距离  
2. double distance(double x, double y) {  
3.     return sqrt(x*x + y*y);  
4. }  

        修改计算梯度的代码，因为我们只对梯度值大于阈值的像素点感兴趣。找出该点后在hough图上对经过 p = x*cos(a)+y*sin(a) 的点值加1：

[cpp]  view plain  copy

 

1. // 定义3*3领域矩阵I  
2. CImg_3x3(I, double);  
3. // 遍历计算梯度值  
4. cimg_for3x3(grayImg, x, y, 0, 0, I, double) {  
5.     const double ix = Inc - Ipc;  
6.     const double iy = Icp - Icn;  
7.     double grad = std::sqrt(ix*ix + iy*iy);  
8.     // 梯度大于阈值才赋值  
9.     if (grad > GRADLIMIT) {  
10.         gradnum(x, y) = grad;  
11.         cimg_forX(hough, angle) {      // 遍历x轴计算y，[x,y]+=1  
12.             double rangle = (double)angle*PI / 180.0;  
13.             int polar = (int)(x*cos(rangle) + y*sin(rangle));  
14.             if (polar >= 0 && polar < hough.height()) {      // 确定y在图像高度内  
15.                 hough(angle, polar) += 1;  
16.             }  
17.         }  
18.     }  
19. }  

        显示hough图：

[cpp]  view plain  copy

 

1. hough.display();  

        完成这部分代码，运行后应有如下结果：（右边是局部放大）

        

五、提取Hough峰值。

        在hough图上，我们可以找到4个亮点，那就是我们希望得到的4条直线。但是每个亮点放大后看，都不是规规矩矩的一个点，而是一簇相邻的亮点，所以要进行提取。Hough峰值的提取有很多种办法，比较理想的是质心算法。我使用了一种较为简单的算法，只是单纯的找出范围内最亮的点，这个范围通过一个阈值DIFF表示：

[cpp]  view plain  copy

 

1. #define DIFF 200  

        我希望把找到的峰值点信息都保存起来，所以使用了STL的Vector类，并且定义了一个新结构体Dot：

[cpp]  view plain  copy

 

1. // 2D Dot(x, y) [value]    
2. struct Dot {  
3.     int x, y, value;  
4.     Dot(int _x, int _y, int _value)  
5.         :x(_x), y(_y), value(_value){}  
6. };  

        另外要注意，我们的峰值所表示的关于x和y的直线方程 并不一定合法 ，一般是由极坐标转换回直线坐标的正负误差等造成的。因此首先要进行一番筛选，判断直线是否经过我们的图像。我的判断方法是，如果直线与图像任意一条边有交点，则直线是合法的。该方法涉及两个函数CrossX和CrossY，分别计算直线在指定x或者y上的交点。先看峰值提取代码：

[cpp]  view plain  copy

 

1. // Find peaks  
2.     std::vector<Dot*> peaks;  
3.     cimg_forXY(hough, angle, polar) {  
4.         if (hough(angle, polar) > THRESHOLD) {          // 是否是峰值  
5.             bool flag = false;  
6.             const int ymin = 0;  
7.             const int ymax = image.height() - 1;  
8.             const int x0 = CrossY(angle, polar, ymin);  
9.             const int x1 = CrossY(angle, polar, ymax);  
10.   
11.             const int xmin = 0;  
12.             const int xmax = image.width() - 1;  
13.             const int y0 = CrossX(angle, polar, xmin);  
14.             const int y1 = CrossX(angle, polar, xmax);  
15.   
16.             if (x0 >= 0 && x0 <= xmax ||              // 表示的直线是否在图像内  
17.                 x1 >= 0 && x1 <= xmax ||  
18.                 y0 >= 0 && y0 <= ymax ||  
19.                 y1 >= 0 && y1 <= ymax) {  
20.                 for (int i = 0; i < peaks.size(); ++i) {    // 遍历数组，找相邻峰值  
21.                     if (distance(peaks[i]->x - angle, peaks[i]->y - polar) < DIFF) {   // 存在相邻峰值  
22.                         flag = true;  
23.                         if (peaks[i]->value < hough(angle, polar)) {  // 如果比相邻峰值还大  
24.                             peaks[i] = new Dot(angle, polar, hough(angle, polar));   // 替换为当前峰值  
25.                         }  
26.                     }  
27.                 }  
28.                 if (flag == false) {        // 当前峰值无相邻峰值  
29.                     peaks.push_back(new Dot(angle, polar, hough(angle, polar)));   // 加入新峰值  
30.                 }  
31.             }  
32.         }  
33.     }  

        CrossX和CrossY的定义如下：

[cpp]  view plain  copy

 

1. // Polar coordinate intersection at x  
2. const int CrossX(int theta, int distance, int x) {  
3.     double angle = (double)theta*PI / 180.0;  
4.     double m = -cos(angle) / sin(angle);  
5.     double b = (double)distance / sin(angle);  
6.     return m*x + b;  
7. }  
8.   
9. // Polar coordinate intersection at y  
10. const int CrossY(int theta, int distance, int y) {  
11.     double angle = (double)theta*PI / 180.0;  
12.     double m = -cos(angle) / sin(angle);  
13.     double b = (double)distance / sin(angle);  
14.     return ((double)(y - b) / m);  
15. }  

        如果你的阈值设置得当，那么你会得到一个含有4个值的峰值数组。

六、计算直线方程和角点坐标。

        有了峰值数组，实际上我们已经得到了A4纸4条边缘的极坐标直线方程，我们要把它们转换成直角坐标方程，并且计算直线的交点：

[cpp]  view plain  copy

 

1. // Transform polar coordinates to rectangular coordinates  
2.     std::vector<Line*> lines;  
3.     for (int i = 0; i < peaks.size(); ++i) {  
4.         double angle = (double)peaks[i]->x*PI / 180.0;  
5.         double m = -cos(angle) / sin(angle);  
6.         double b = (double)peaks[i]->y / sin(angle);  
7.         lines.push_back(new Line(m, b));  
8.         std::cout << "y = (" << m << ") x + (" << b << ")" << std::endl;  
9.     }  
10.   
11.     std::cout << std::endl << "intersections: " << std::endl;  
12.   
13.     // Calculate line intersections  
14.     std::vector<Dot*> intersections;  
15.     for (int i = 0; i < lines.size(); ++i) {  
16.         for (int j = i + 1; j < lines.size(); ++j) {  
17.             double m0 = lines[i]->m;  
18.             double m1 = lines[j]->m;  
19.             double b0 = lines[i]->b;  
20.             double b1 = lines[j]->b;  
21.             double x = (b1 - b0) / (m0 - m1);  
22.             double y = (m0*b1 - m1*b0) / (m0 - m1);  
23.             if (x >= 0 && x < result.width() && y >= 0 && y < result.height()) {  
24.                 intersections.push_back(new Dot(x, y, 0));  
25.                 std::cout << "(" << x << ", " << y << ")" << std::endl;  
26.             }  
27.         }  
28.     }  

        其中Line是新定义的结构体：

[cpp]  view plain  copy

 

1. // 2D Straight Line(x, y) [y = m*x + b]  
2. struct Line {  
3.     double m, b;  
4.     Line(double _m, double _b)  
5.         :m(_m), b(_b){}  
6. };  

        完成这部分代码，运行后可以在控制台看到直线方程和角点坐标：

        Oops！看来之前只找出了3条直线，不过不要紧，只需改一下阈值就好，后面会详细讲解如何调阈值。

七、绘图。

        有了直线方程和角点坐标，我想把他们绘制在原图上，这样可以有直观的感受，知道是否出现了问题。直线和点的绘制可以直接使用CImg的draw_line和draw_circle函数：

[cpp]  view plain  copy

 

1. CImg<double> result(image);  
2.   
3. // Draw lines  
4. for (int i = 0; i < lines.size(); ++i) {  
5.     const int ymin = 0;  
6.     const int ymax = result.height() - 1;  
7.     const int x0 = (double)(ymin - lines[i]->b) / lines[i]->m;  
8.     const int x1 = (double)(ymax - lines[i]->b) / lines[i]->m;  
9.   
10.     const int xmin = 0;  
11.     const int xmax = result.width() - 1;  
12.     const int y0 = xmin*lines[i]->m + lines[i]->b;  
13.     const int y1 = xmax*lines[i]->m + lines[i]->b;  
14.   
15.     const double color[] = { 255, 0, 255 };  
16.   
17.     if (abs(lines[i]->m) > SLOPE_FLAG) {  
18.         result.draw_line(x0, ymin, x1, ymax, color);  
19.     }  
20.     else {  
21.         result.draw_line(xmin, y0, xmax, y1, color);  
22.     }  
23. }  
24.   
25. // Draw intersections  
26. for (int i = 0; i < intersections.size(); ++i) {  
27.     const double color[] = { 255, 100, 255 };  
28.     result.draw_circle(intersections[i]->x, intersections[i]->y, 50, color);  
29. }  
30.   
31. result.display();  

        注意特殊情况，要考虑直线是水平或者垂直情况下的绘制。

        完成绘制后，我们终于可以看见结果了：

        正如我们第6步结束后看到的，我们只找到了A4纸的3条边缘，左侧的边缘似乎被遗弃掉了。这是由于阈值设置不得当引起的，如果足够细心的话，会发现原图A4纸左侧是有阴影的，因此得到的梯度图中，左侧的点比较少，所以在hough图上亮点就不明显。

        目前为止，整个程序的阈值有：

[cpp]  view plain  copy

 

1. #define BLUR 3  
2. #define GRADLIMIT 20  
3. #define THRESHOLD 800  
4. #define DIFF 200  
5. #define SLOPE_FLAG 1  

        其中SLOPE_FLAG是用来判断直线是否水平的，因此可以不用管它。我们来讨论剩下的4个阈值。
        

        BLUR：模糊参数，通常在原图噪声点较多的时候适当增大。

        GRADLIMIT：梯度阈值，如果边缘与周围物体的区分不明显，则应适当减小，反之则适当增大。

        THRESHOLD：Hough峰值阈值，这里应通过调整阈值来获取到刚好4个峰值。

        DIFF：峰值提取时判断相邻峰值的距离阈值，如果一条边提取后有两个斜率十分相近的直线方程，那么应增大该阈值。

        这里，我们应该减小THRESHOLD阈值，从原先的800降到650，这样就能得到正确的结果：

        直线的方程和角点坐标显示在命令行窗口：

        结果图的左侧边缘似乎并没有完全贴合我们的A4纸，只是由于我的峰值提取算法导致的，如果采用质心算法效果会更好。

        Hough的强大之处是，如果我们把A4纸的一角盖住了，我们也能预测出A4纸原来的样子，仔细想想，是不是这样？

        好了，今天给出的是满满的干货，要吸收还得多打打代码。我写得太急，代码一点都不优雅，大家千万别学我，最好把各个功能都结构化，方便日后重用。

        最后贴出我的完整代码吧：

[cpp]  view plain  copy

 

1. #include "CImg.h"  
2. using namespace cimg_library;  
3. #include <vector>  
4. #include <iostream>  
5.   
6. #define FILE "A4.bmp"  
7.   
8. #define BLUR 3  
9. #define GRADLIMIT 20  
10. #define THRESHOLD 650   
11. #define DIFF 200  
12. #define SLOPE_FLAG 1  
13. #define PI 3.14159265358979323846  
14.   
15. // 2D Dot(x, y) [value]    
16. struct Dot {  
17.     int x, y, value;  
18.     Dot(int _x, int _y, int _value)  
19.         :x(_x), y(_y), value(_value){}  
20. };  
21.   
22. // 2D Straight Line(x, y) [y = m*x + b]  
23. struct Line {  
24.     double m, b;  
25.     Line(double _m, double _b)  
26.         :m(_m), b(_b){}  
27. };  
28.   
29. // Calculate the distance  
30. double distance(double x, double y) {  
31.     return sqrt(x*x + y*y);  
32. }  
33.   
34. // Polar coordinate intersection at x  
35. const int CrossX(int theta, int distance, int x) {  
36.     double angle = (double)theta*PI / 180.0;  
37.     double m = -cos(angle) / sin(angle);  
38.     double b = (double)distance / sin(angle);  
39.     return m*x + b;  
40. }  
41.   
42. // Polar coordinate intersection at y  
43. const int CrossY(int theta, int distance, int y) {  
44.     double angle = (double)theta*PI / 180.0;  
45.     double m = -cos(angle) / sin(angle);  
46.     double b = (double)distance / sin(angle);  
47.     return ((double)(y - b) / m);  
48. }  
49.   
50. int main() {  
51.     CImg<double> image(FILE);  
52.     CImg<double> grayImg(image);  
53.     CImg<double> result(image);  
54.   
55.     // Gray scale;  
56.     cimg_forXY(grayImg, x, y) {  
57.         double R = grayImg(x, y, 0, 0);  
58.         double G = grayImg(x, y, 0, 1);  
59.         double B = grayImg(x, y, 0, 2);  
60.         double Gray = (R * 299 + G * 587 + B * 114 + 500) / 1000;  
61.         grayImg(x, y, 0, 0) = Gray;  
62.         grayImg(x, y, 0, 1) = Gray;  
63.         grayImg(x, y, 0, 2) = Gray;  
64.     }  
65.   
66.     // Blur  
67.     grayImg.blur(BLUR);  
68.   
69.     CImg<double> gradnum(image.width(), image.height(), 1, 1, 0);  
70.     int maxDistance = distance(image.width(), image.height());  
71.     CImg<double> hough(360, maxDistance, 1, 1, 0);  
72.     // 定义3*3领域矩阵I  
73.     CImg_3x3(I, double);  
74.     // 遍历计算梯度值  
75.     cimg_for3x3(grayImg, x, y, 0, 0, I, double) {  
76.         const double ix = Inc - Ipc;  
77.         const double iy = Icp - Icn;  
78.         double grad = std::sqrt(ix*ix + iy*iy);  
79.         // 梯度大于阈值才赋值  
80.         if (grad > GRADLIMIT) {  
81.             gradnum(x, y) = grad;  
82.             cimg_forX(hough, angle) {  
83.                 double rangle = (double)angle*PI / 180.0;  
84.                 int polar = (int)(x*cos(rangle) + y*sin(rangle));  
85.                 if (polar >= 0 && polar < hough.height()) {  
86.                     hough(angle, polar) += 1;  
87.                 }  
88.             }  
89.         }  
90.     }  
91.   
92.     // Find peaks  
93.     std::vector<Dot*> peaks;  
94.     cimg_forXY(hough, angle, polar) {  
95.         if (hough(angle, polar) > THRESHOLD) {          // 是否是峰值  
96.             bool flag = false;  
97.             const int ymin = 0;  
98.             const int ymax = image.height() - 1;  
99.             const int x0 = CrossY(angle, polar, ymin);  
100.             const int x1 = CrossY(angle, polar, ymax);  
101.   
102.             const int xmin = 0;  
103.             const int xmax = image.width() - 1;  
104.             const int y0 = CrossX(angle, polar, xmin);  
105.             const int y1 = CrossX(angle, polar, xmax);  
106.   
107.             if (x0 >= 0 && x0 <= xmax ||              // 表示的直线是否在图像内  
108.                 x1 >= 0 && x1 <= xmax ||  
109.                 y0 >= 0 && y0 <= ymax ||  
110.                 y1 >= 0 && y1 <= ymax) {  
111.                 for (int i = 0; i < peaks.size(); ++i) {    // 遍历数组，找相邻峰值  
112.                     if (distance(peaks[i]->x - angle, peaks[i]->y - polar) < DIFF) {   // 存在相邻峰值  
113.                         flag = true;  
114.                         if (peaks[i]->value < hough(angle, polar)) {  // 如果比相邻峰值还大  
115.                             peaks[i] = new Dot(angle, polar, hough(angle, polar));   // 替换为当前峰值  
116.                         }  
117.                     }  
118.                 }  
119.                 if (flag == false) {        // 当前峰值无相邻峰值  
120.                     peaks.push_back(new Dot(angle, polar, hough(angle, polar)));   // 加入新峰值  
121.                 }  
122.             }  
123.         }  
124.     }  
125.     std::cout << "lines: " << std::endl;  
126.   
127.     // Transform polar coordinates to rectangular coordinates  
128.     std::vector<Line*> lines;  
129.     for (int i = 0; i < peaks.size(); ++i) {  
130.         double angle = (double)peaks[i]->x*PI / 180.0;  
131.         double m = -cos(angle) / sin(angle);  
132.         double b = (double)peaks[i]->y / sin(angle);  
133.         lines.push_back(new Line(m, b));  
134.         std::cout << "y = (" << m << ") x + (" << b << ")" << std::endl;  
135.     }  
136.   
137.     std::cout << std::endl << "intersections: " << std::endl;  
138.   
139.     // Calculate line intersections  
140.     std::vector<Dot*> intersections;  
141.     for (int i = 0; i < lines.size(); ++i) {  
142.         for (int j = i + 1; j < lines.size(); ++j) {  
143.             double m0 = lines[i]->m;  
144.             double m1 = lines[j]->m;  
145.             double b0 = lines[i]->b;  
146.             double b1 = lines[j]->b;  
147.             double x = (b1 - b0) / (m0 - m1);  
148.             double y = (m0*b1 - m1*b0) / (m0 - m1);  
149.             if (x >= 0 && x < result.width() && y >= 0 && y < result.height()) {  
150.                 intersections.push_back(new Dot(x, y, 0));  
151.                 std::cout << "(" << x << ", " << y << ")" << std::endl;  
152.             }  
153.         }  
154.     }  
155.   
156.     std::cout << std::endl;  
157.   
158.     // Draw lines  
159.     for (int i = 0; i < lines.size(); ++i) {  
160.         const int ymin = 0;  
161.         const int ymax = result.height() - 1;  
162.         const int x0 = (double)(ymin - lines[i]->b) / lines[i]->m;  
163.         const int x1 = (double)(ymax - lines[i]->b) / lines[i]->m;  
164.   
165.         const int xmin = 0;  
166.         const int xmax = result.width() - 1;  
167.         const int y0 = xmin*lines[i]->m + lines[i]->b;  
168.         const int y1 = xmax*lines[i]->m + lines[i]->b;  
169.   
170.         const double color[] = { 255, 0, 255 };  
171.   
172.         if (abs(lines[i]->m) > SLOPE_FLAG) {  
173.             result.draw_line(x0, ymin, x1, ymax, color);  
174.         }  
175.         else {  
176.             result.draw_line(xmin, y0, xmax, y1, color);  
177.         }  
178.     }  
179.   
180.     // Draw intersections  
181.     for (int i = 0; i < intersections.size(); ++i) {  
182.         const double color[] = { 255, 100, 255 };  
183.         result.draw_circle(intersections[i]->x, intersections[i]->y, 50, color);  
184.     }  
185.   
186.     result.display();  
187.   
188.     // Display  
189.     //grayimg.display();  
190.     //gradnum.display();  
191.     //hough.display();  
192.     result.display();  
193.   
194. }