In this chapter, we will cover:

Detecting image contours with the Canny operator

Detecting lines in images with the Hough transform （使用Hough变换检测直线）

Fitting a line to a set of points （用直线拟合一组点）

Extracting the components’ contours（提取连通区域的轮廓）

Computing components’ shape descriptors

使用Canny算子检测轮廓

Canny的使用

　　OpenCV中Canny算法的函数是cv::Canny，该算法需要指定两个阈值，调用如下:

1 2	cv::Mat contours; cv::Canny(image, contours, 125, 350); //125: low threshold, 350: high threshold

Read more here. Canny边缘检测的结果是轮廓用非零元素表示，书上为了方便打印，所以反转了颜色。

theroy

　　参考wikipedia的解释，canny 边缘检测主要分为以下5个步骤:

Apply Gaussian filter to smooth the image in order to remove the noise

Find the intensity gradients of the image

Apply non-maximum suppression to get rid of spurious response to edge detection

Apply double threshold to determine potential edges

Track edge by hysteresis: Finalize the detection of edges by suppressing all the other edges that are weak and not connected to strong edges.

　　先用高斯滤波去噪，然后计算梯度(opencv中使用sobel算子实现)，利用non-maximum suppression(非极大值抑制)的方法得到精细的边缘(仅仅使用sobel算子处理后的结果不精细，会有需要聚集在一起的块，而非边缘，所以使用NMS来取局部最大值，NMS只保留局部最大值，其余的全部使其为0)，这时候得到的结果可能仍然包含噪声的边缘，所以使用double threshold进行修正，梯度值大于high threshold value的我们认为肯定是strong edge pixels，大于low threshold value且小于high threshold value的认为是weak edge pixels，小于low threshold value，就会被抑制，使其为0.但是weak edge pixels里可能仍然有噪点的边缘，所以使用Binary Large Object (BLOB) analysis，看weak edge pixel的八连通区域里是否有strong edge pixel，只要有一个，就认为这个weak edge pixel是边缘，应该保留，否则抑制为0.最后得到的结果就是canny边缘检测算法的结果。此外链接里还讲了canny算法的提高，这里不在叙述，有兴趣请自行查看。

　　其中书上也说了磁滞阈值化（hysteresis thresholding,就是需要两个阈值，一个高阈值，一个低阈值）之前，进行non-maximum suppression，所以canny能得到很薄的边缘，如果还是不理解，把sobel算子的结果跑出来和canny对比一下就会发现，sobel算子的结果比较粗，不精细。

使用Hough变换检测直线

　　Hough Transform是检测直线的经典算法，它最初只用于检测直线，被扩展后能够检测其他的简单结构，如圆等。

Related Work

　　在Hough Transform中，直线常用下列方程表示：
$$p = x \cos \theta + y \sin \theta \tag{1} $$

　　　　　　　　　　　　　　　　
　　　　　　　　　　　　　　　　　　　　　　　　图1

　　如图1所示，参数$p$是直线(直线上的点)到图像原点（左上角）的最短距离（过原点做直线的垂线，该垂线就是最短距离），$\theta$是垂线绕x轴正方向转过的角度。$\theta$取值范围是$[0, \pi)$,$p$的最大值是图像的对角线长度，但也有可能是负值，下面说。
　　直线1的垂直线的角度$\theta$等于0，直线1也同时解释为何我们不用$y=kx+b$来表达直线，而采取极坐标系，因为这时候斜率$k$无穷大，无法表示，所以就换一种方式了；直线5的$\theta$等于$\pi /2$;直线3的$\theta$等于$\pi /4$；直线4的$\theta$约等于$0.7 \pi$。而直线2的$\theta$约等于$0.8 \pi$，大于135度，这时候$p$就是负值，在$\theta$取值范围在$[3\pi / 4, \pi)$时，可以把上述公式这样写成$y = -(\cos \theta / \sin \theta)x + p / \sin \theta $,因为图像的坐标系是向下为正，以左上角为原点，所以自然截矩$p / \sin \theta$是负值，其中$\sin \theta$在取值范围内是正数，所以p是负值。

应用

　　OpenCV中提供两种Hough变换的实现,基础版本是cv::HoughLines, 它要求的输入图像是八位单通道的二进制图像，但是即使你按照要求输入也行，函数会帮你改，非0的变为1；此外他的输出是一个向量$(p, \theta)$，通常保存在vector里面，具体见下面代码；参数rho就是移动步长，一般设为1；参数theta是转动角度，一般都是直线每转一度就算一次，看看该直线经过多少点，这样在区间$[0, \pi)$内我们都计算了一下，避免遗漏；参数threshold是阈值，我们前面在每个rho和theta下面就会计算该直线经过多少点，通过阈值来去除一些不满足要求的点，只有该直线上经过的点大于阈值所定数目时，我们才会认为他是直线，阈值怎么定，看实际情况吧。See more here.

代码如下：

// [20151225]Opencv_chapter7.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"

#include <opencv2/opencv.hpp>

const int PI = 3.14;

void DrawHoughLines(cv::Mat &iamge, std::vector<cv::Vec2f> lines)
{
	printf("%d", lines.size());
	for (size_t i = 0; i < lines.size(); ++i)
	{
		float rho = lines[i][0], theta = lines[i][1];
		cv::Point pt1, pt2;
		double	a = cv::cos(theta), b = cv::sin(theta);
		double x0 = a * rho, y0 = b *rho;
		printf("1 is %f, 2 is %f\n", 1000 * (-b), 1000 * (a));
		pt1.x = cvRound(x0 + 1000 * (-b));
		pt1.y = cvRound(y0 + 1000 * (a));
		pt2.x = cvRound(x0 - 1000 * (-b));
		pt2.y = cvRound(y0 - 1000 * (a));
		cv::line(iamge, pt1, pt2, cv::Scalar(0, 0, 255), 1, CV_AA);
	}
}

void DrawHoughLines2(cv::Mat &image, std::vector<cv::Vec2f> lines)
{
	for (size_t i = 0; i < lines.size(); ++i)
	{
		float rho = lines[i][0], theta = lines[i][1];
		cv::Point pt1(0, 0), pt2(0, 0);
		pt1.x = rho / cv::cos(theta);	pt1.y = 0;
		pt2.x = (rho - image.rows * cv::sin(theta)) / cv::cos(theta);	pt2.y = image.rows;
// 		if (theta < CV_PI / 4 || theta > 3 * CV_PI / 4)
// 		{
// 			pt1.x = rho / cv::cos(theta);	pt1.y = 0;
// 			pt2.x = (rho - image.rows * cv::sin(theta)) / cv::cos(theta);	pt2.y = image.rows;
// 		}
// 		else
// 		{
// 			pt1.x = 0;		pt1.y = rho / cv::sin(theta);
// 			pt2.x = image.cols;		pt2.y = (rho - image.cols * cv::cos(theta)) / cv::sin(theta);
// 		}
		cv::line(image, pt1, pt2, cv::Scalar(0, 0, 255), 1, CV_AA);
	}
}

int _tmain(int argc, _TCHAR* argv[])
{
	cv::Mat image = cv::imread("../../resource2/road.jpg");
	cv::Mat imageCopy = image.clone();

	cv::Mat edge;
	cv::Canny(imageCopy, edge, 125, 350);

	/*
	cv::Mat dx, dy;
	cv::Sobel(imageCopy, dx, imageCopy.depth(), 1, 0, 3);
	cv::Sobel(imageCopy, dy, imageCopy.depth(), 0, 1, 3);
	cv::threshold(dx, dx, 128, 255, CV_THRESH_BINARY_INV);
	cv::threshold(dy, dy, 128, 255, CV_THRESH_BINARY_INV);
*/
//	cv::imshow("dx", dx);
//	cv::imshow("dy", dy);

	cv::Mat result;
	std::vector<cv::Vec2f> lines;
	cv::HoughLines(edge, lines, 1, CV_PI / 180, 80);
	DrawHoughLines(imageCopy, lines);
	DrawHoughLines2(imageCopy, lines);

	cv::Mat test;
	test = cv::Mat::zeros(500, 500, CV_8UC3);
	cv::Point pt1(0, 552), pt2(500, 203);
	cv::line(test, pt1, pt2, cv::Scalar(0, 0, 255), 3, CV_AA);

//	cv::imshow("result", result);
	cv::imshow("edge", edge);
	cv::imshow("iamgeCopy", imageCopy);
	cv::imshow("test", test);

	cv::waitKey(0);

	return 0;
}

　　进行边缘检测后使用Hough变换做直线检测，最后调用画线函数将其画出来，这里给了两种画线的方式。
　　第一种是OpenCV自带的demo里的方式，可能刚开始看到后有点不理解，为何搞个1000，如果不理解可以看这里。
　　第二种方法是书上的，是分开角度来算，我测试了一下，发现完全没有必要，你把if里面的或者else里面的任意一个单独拎出来都可以计算出结果，刚开始我还比较疑惑，比如对于一个500*500的图片，如果给一个点$(0,525)$,一个点$(550,203)$,那这样岂不是越界了，点都不在图像范围内，后来发现OpenCV提供的画线函数cv::line可以hold住，只要你提供两个点就ok，正事因为这个函数的强大，所以就不用像书上那样麻烦，其实书上搞的这个东西说解决近似水平或者近似竖直的直线，我觉得这个好像是有点…….因为如果一条直线经过图像的最后一行，最后一列，按照书上的意思他的这个if..else是hold不住的，但是其实没问题，没问题是因为cv::line函数的强大，和他的这个if..else是没有关系的。当然这是我的观点，难免会有边界条件等没考虑到，如果你发现错误，请予以指正，不胜感激。也许你可能不理解这句代码pt2.x = (rho - image.rows * cv::sin(theta)) / cv::cos(theta);,但是如果这样写或许就会明白了，pt2.x = rho/cv::cos(theta) - image.rows * cv::sin(theta)/cv::cos(theta);如图2所示：

　　　　　　　　　　　　　　　　　　　　　　　　　　　　图2

另外那个同理，可以自己画画看，推倒一下。

　　然而由于基本的Hough Transform可能由于意外的像素排列或者多条直线经过同一组像素,有误检，所以提出概率Hough变换(Probabilistic Hough Transform),对应的函数是cv::HoughLinesP,大部分参数一致，不同的是输出为$(x_1, y_1, x_2, y_2)$,$(x_1, y_1)$和$(x_2, y_2)$是线的两个端点坐标；此外添加了两个阈值，一个minLineLength,表示检测到的直线长度小于该阈值则被out，一个maxLineGap，表示在同一直线上点之间的距离不能超过这个值，这样一条直线同样是经过十个点，如果这十个点之间的距离大于maxLineGap我们就认为这并不是一条直线，如果小于则认为是直线。
代码如下：

std::vector<cv::Vec4i> lines;
cv::HoughLinesP(edge, lines, 1, CV_PI / 180, 80, 30, 10);
for (size_t i = 0; i < lines.size(); ++i)
{
	cv::line(imageCopy, cv::Point(lines[i][0], lines[i][1]), cv::Point(lines[i][2], lines[i][3]), cv::Scalar(0, 0, 255), 1, 8);
}
printf("%d", lines.size());

Theory

Hough Transform用来检测直线，这里可以参考wikipedia里的资料，可以直接看Example部分，我将链接里的图引用过来如图3所示：

图3
在图3中我们可以看出，它在每一个点(在该图中一共三个，在二值图像中，计算每一个点)上每次转动30度来计算距离(我们在上面的代码中中计算是以每次转动1度来计算距离，防止有遗漏，这里仅作示例，所以转动30度)，这样的话，我们就可以用一个矩阵来存储这些数据，如果有相同的角度和距离就对其加1，这样最后看矩阵里的值，我们定义一个阈值，超过该阈值就认为这是一条直线即可。图3中可以看出在60度，距离为80左右时，累加次数我们可以算做3(不要钻牛角尖，81.2,80.5,80.6都可以近似为80了，毕竟其他数据差距比这大多了)，那我们就认为该直线经过这三个点，链接里下面的那个图(三条曲线的交点)也说明这一点。
引用书上的话说：Hough Transform的目的是找到二值图像中经过足够数量的点的所有直线，它分析每一个点（注意不是分析每一个像素点，从代码中我们可以看出输入的是边缘图像，即二值图像，我们可以理解成我们要检测这些边缘中有没有直线边缘，所以这里的每一个点指的是边缘上的点，而非图像中每一个像素点，所以需要逐行扫描看哪些是边缘点，一般情况下边缘点都是用非零像素标记，而HoughLines函数也是默认非零像素点才是边缘点，如果使用threshold函数把像素值颜色反转，就会造成结果很大差别，因为本来的零像素值点变成了非零像素了，有兴趣的试试，并识别出所有可能经过它的直线。当同一条直线穿过许多点时，就可以认为检测到这条直线了。Hough transform使用二维的累加器来统计特定的直线被识别多少次，累加器的大小由$(p, \theta)$决定，比如定义步长$p=1$和$\theta = \pi/180$，在一个图像的对角线为200的情况下，那么这个累加器的尺寸就是$200*180$（对角线是距离原点最远的了，所以是200），然后我们就统计每一个点转动180次的结果，然后在累加器的对应位置加1，在统计完很多次以后，累加器中肯定有的位置的值比较大，有的比较小，值比较大的位置，他的坐标就是$(p, \theta)$,也就是说我们只要在这累加器中找到比规定阈值大的，就能得到相应的$(p, \theta)$，那自然直线就得到了。这就是Hough变换的原理。
此外书上也给出了阈值比较小时，就会检测到更多的直线，然而这些都是误检，所以看出阈值很重要。
概率Hough Transform在上面基本的算法中增加了少许改动。首先不再是逐行扫描图像而是随机挑选像素点。一旦累加器中某一项达到给定的最小值(参数threshold)，就扫描沿着对应直线的像素并移出所有经过的点（即便他们并没有被计算过）。这次扫描决定了接受线段的长度。为此，定义了两个额外的参数，之前已经介绍过了。额外的算法步骤增加了复杂度，但是由于参与投票的点减少得到部分的补偿。

Circles Hough Transform(CHT)

Theory

　　可以参考Circles Hough Transform.在二维空间里，圆的方程表示为如下形式：
$$(x-a)^2 + (y-b)^2 = r^2 \tag{2} $$
　　$(a,b)$是圆心位置，$r$是半径，不像直线方程，在二维参数空间下就可以表示(即二维矩阵就能hold住，$(r, \theta)$),这里未知的参数是3个，所以需要三维的参数空间才行，然而三维参数空间想找局部最大值很难，不像二维那么好找（找局部最大值就是为了找到圆，对应局部最大值的三个坐标就能确定一个圆，类似前面直线检测里对应局部最大值的两个坐标就能确定一条直线）。那既然三维空间不好找，大家就折中一下，换个方式找，怎么弄的，先不管待会说，先说怎么在半径确定的情况下怎么找圆心。
　　根据数学知识，如果同一平面上的三个点共圆的画，假设该圆半径为$r$,然后以这三个点为圆心画半径为$r$的圆，则画出来的三个圆的交点即为这个三个点的圆心。那这样就好办了，比如一个平面上有10个点，其中九个点共圆，一个点在其他位置，那么在我们已知半径为$r$的情况下，我们就在每个点上画半径为$r$的圆，在敲代码时，这个平面肯定是有限大的，可以用个二维数组（矩阵）来存储，然后在每个点上画圆时，圆上没经过一处，我们就在那个位置累加1，这个最后画完之后，角点出的累计值就为9，那么最后我们在这个二维数组里找局部最大值（超过某个阈值）还是可以比较快的找到的（这个过程和直线检测还是蛮相似的）。
　　那么新的问题来了，上面所说的都是已知半径，那么如果不知道半径怎么办，那就一个个试，比如图像是$600*600$的，那么圆的半径肯定不会超过300了，那就可以从$[1, 300]$这个范围里找，实际OpenCV里的函数提供了参数，让你自己定义可能的半径范围。

应用

　　OpenCV里提供的函数是cv::HoughCircles, see more here in detail.
　　输入和cv::HoughLines一样，都是边缘检测后的边缘图像（一般多用cv::Canny）；输出是一个三维向量$(x,y,r)$,即圆心坐标和半径；Method参数一般多用CV_HOUGH_GRADIENT,dp参数我暂时不太理解，但是在Stackoverflow上写了个提问，等待大神回答，大神已回答，可以参考，底下有我一点自己的理解，如有错误，恳请指正。minDist参数就是你检测的圆的圆心之间的最小距离，所有检测出来的圆，圆心小于该距离的会被reject。param1 参数就是之前Canny边缘检测里面较大的那个参数，他们两值一样。param2 参数就是阈值，检测到圆的投票数，低于该阈值的被reject。minRadius 和maxRadius 就是圆的最小半径和最大半径，按照前面的解释，在这个范围里找。
　　关于dp参数的理解，根据上面链接里大神的讲解，dp参数的作用是控制累加器的size，dp越大，累加器的size越小，即dp=2，累加器的宽高只有图像的一半，dp=3，累加器的宽高只有图像的三分之一。我之前一直不理解当累加器的大小与图像不一样时是怎么在hough space下投票的，根据大神说的，相近的像素被合并，这样我就懂了。前面说了，我们在边缘点上画圆，然后圆上经过的点都会被vote，那时候刚好一个像素点对应矩阵（hough space）里的一个位置，我懂只要在对应位置加1即可，但是累加器缩小以后，怎么办呢？其实可以这样理解，还是在跟图像一样大小的累加器里投票，然后把这个累加器按照dp参数的大小来缩小，这样比如本来某几个邻近位置的投票数就会被合并到一起，这样导致得不到非常精确的圆，但是这个投票是robust，就是说虽然得不到非常精确的，但是我能保证我检测到的这个圆很大很大可能性是存在的，而且和真实的圆很接近。大神还补充了OpenCV里的实现是用梯度来算的，不是用标准hough变换来做的，如何使用梯度来算，可以参考这个文章，讲的挺好。补充一下，在这个链接里可以看到使用梯度计算时，有时因为误差，导致同一个圆上的点映射的圆心不在同一个位置，但是比较接近，这时候如果按照前面说的dp参数的作用，把dp设为2，或者更大的数，那么累加器中相近的位置被合并，那么是不是就可以确定一个唯一的圆心了，是不是，如果不合并即dp=1，那么可能vote数量没有超过阈值而导致检测不出来圆或者一下子出来好几个圆。

代码我是直接把文档里的拷贝过来了，如下：

#include <cv.h>
#include <highgui.h>
#include <math.h>

using namespace cv;

int main(int argc, char** argv)
{
    Mat img, gray;
    if( argc != 2 && !(img=imread(argv[1], 1)).data)
        return -1;
    cvtColor(img, gray, CV_BGR2GRAY);
    // smooth it, otherwise a lot of false circles may be detected
    GaussianBlur( gray, gray, Size(9, 9), 2, 2 );		//先降噪，减少噪点的影响
    vector<Vec3f> circles;				//注意，这里必须是Vec3f
    HoughCircles(gray, circles, CV_HOUGH_GRADIENT,
                 2, gray->rows/4, 200, 100 );
    for( size_t i = 0; i < circles.size(); i++ )
    {
         Point center(cvRound(circles[i][0]), cvRound(circles[i][1]));
         int radius = cvRound(circles[i][2]);
         // draw the circle center
         circle( img, center, 3, Scalar(0,255,0), -1, 8, 0 );
         // draw the circle outline
         circle( img, center, radius, Scalar(0,0,255), 3, 8, 0 );
    }
    namedWindow( "circles", 1 );
    imshow( "circles", img );
    return 0;
}

广义hough变换

伯克利的讲义关于Hough 变换的，其中广义hough变换(Generalised Hough transform,GHT)说的比较浅显易懂，wikipedia里的讲的也不错.

Fitting a line to a set of points

本节主要讲给定一组点，如何找出那条最贴近这组点的直线（有点像最小二乘法）。

实现方法

//points可以是一组二维或三维的点集，例std::Vector<cv::Point>points;
cv::Vec4f line;
cv::fitLine(cv::Mat(points),line,
CV_DIST_L2, // distance type
0, // not used with L2 distance
0.01,0.01); // accuracy

cv::fitLine的各个参数含义见文档.其中输出line为(vx, vy, x0, y0)（二维）或(vx, vy, vz, x0, y0, z0)（三维），(vx, vy)是一组与拟合直线平行的向量，相当于提供斜率（即vy/vk）不过他们的绝对值应该都是小于1的，(x0, y0)就是拟合直线上的一点。
想要将该直线画出来，引用书上的代码如下：

int x0= line[2]; // a point on the line
int y0= line[3];
int x1= x0-200*line[0]; // add a vector of length 200
int y1= y0-200*line[1]; // using the unit vector
image= cv::imread("../road.jpg",0);
cv::line(image,cv::Point(x0,y0),cv::Point(x1,y1),
cv::Scalar(0),3);

画出来的方法很多，书上举例，将直线的长度定成200。

实现原理

直线拟合在fitLine函数中，有多个距离函数可供选用，最快的就是欧式距离（Euclidean distance），即参数CV_DIST_L2。该参数对应的是标准的最小二乘法。可以参考文档链接或者M-estimator

椭圆拟合

1 2	cv::RotatedRect rrect= cv::fitEllipse(cv::Mat(points)); cv::ellipse(image,rrect,cv::Scalar(0));

cv::ellipse就能画出这个拟合后的椭圆了。其中返回的矩形的内切椭圆就是拟合的椭圆。

Extracting the components’ contours

实现方法

如何提取连通域，OpenCV给出的实现函数cv::findContours，可参考文档:

std::vector<std::vector<cv::Point> > contours;
cv::findContours(image, contours, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE);
std::cout << "contours.size()" << contours.size() << std::endl;
cv::Mat result(image.size(), CV_8U, cv::Scalar(255));
cv::drawContours(result, contours, -1, cv::Scalar(0), 2);

移出过长或过短的轮廓代码：

int cmin = 100;
int cmax = 1000;
std::vector<std::vector<cv::Point> >::const_iterator itc = contours.begin();
while (itc != contours.end())
{
	if (itc->size() < cmin || itc->size() > cmax)
	{
		itc = contours.erase(itc);
	}
	else
	{
		++itc;
	}
}

将这段代码加在cv::drawContours前面即可。这是在物体检测中，我们知道先验知识（比如知道物体的大小），通过去掉一些较小的或者较大的轮廓达到识别出感兴趣的物体。
cv::findContours的参数mode：

CV_RETR_EXTERNAL 只取外轮廓，忽略孔洞。文档里的hierarchy[i][2]=hierarchy[i][3]=-1，代表没有父轮廓和子轮廓，因为孔洞的轮廓相当于外轮廓的儿子，现在只取外轮廓，自然就没有子轮廓，没有子轮廓，也就没有父轮廓了。
CV_RETR_LIST 遍历所有轮廓，忽略轮廓之间的关系（如父子关系，前后关系（前后关系我翻译的不恰当，英文是previous和next））
CV_RETR_CCOMP 遍历所有的轮廓，但是所有轮廓的关系只有两层（种）。顶层关系中，包含的是外轮廓，第二层关系中，是外轮廓中的孔洞。如果孔洞中还有孔洞，将其加入顶层中。依次类推，孔洞中孔洞中孔洞，应该就是第二层中。（个人理解）
CV_RETR_TREE 用树结构来遍历所有轮廓，外轮廓是根节点，孔洞是其子节点，孔洞的孔洞是其子节点的子节点，层次关系不像上面只有2层。
cv::findContours的参数method：
CV_CHAIN_APPROX_NONE 存储轮廓上所有的点，其中任意两点(x1,y1) and (x2,y2)，他们遵循 max(abs(x1-x2),abs(y2-y1))==1这个规则.
CV_CHAIN_APPROX_SIMPLE 该参数适用与存储特殊图像的轮廓点，比如矩形（其中两条边和x轴垂直的那种矩形，而不是斜着的矩形，个人理解，原文有个up-right），只需要存储四个顶点即可，而不用像上面那样存储所有的点，节省空间
CV_CHAIN_APPROX_TC89_L1,CV_CHAIN_APPROX_TC89_KCOS 使用teh-Chinl chain 近似算法。See TehChin89 for details.