opencv学习 霍夫直线 HoughLines和HoughLinesP

霍夫直线原理

先看一下直角坐标系与参数坐标系之间的转化

看y = k1*x + b1

p1 = sqrt(x1*x1 + y1*y1);

sinθ = y1/p1,cosθ = x1/p1;

k1 = tan(Π/2 + θ1 ) = -cotθ = -cosθ/sinθ

由上可以解出 b1 = p1/sinθ;

y = (-cosθ/sinθ)*x + p1/sinθ;

p1 = y*sinθ + x*cosθ;

其中 p1 = y1*sinθ1 + x1*cosθ1 ,因经过(x0,y0)  p1 = y0*sinθ1 + x0*cosθ1; 

 p2 = y2*sinθ2 + x2*sinθ2; 因经过（x0,y0）p2 = y0*sinθ2 + x0*cosθ2; 

其中其他人证明：

 

θ-r坐标系上，是两条当（x1,y1）和（x2,y2）固定不变时，r随着θ的两条变化曲线，也就是说，θ-r坐标系上，每一点都是一条直线。代码中θ-r坐标系，每个点的位置都是一个累计器，进行计数。当达到阈值，就说明是一个条直线。

 

标准霍夫变换本质上是把图像映射到它的参数空间上，它需要计算所有的M个边缘点，这样它的运算量和所需内存空间都会很大。如果在输入图像中只是处理m（m<M）个边缘点，则这m个边缘点的选取是具有一定概率性的，因此该方法被称为概率霍夫变换（Probabilistic Hough Transform）。该方法还有一个重要的特点就是能够检测出线端，即能够检测出图像中直线的两个端点，确切地定位图像中的直线

opencv源码HoughLinesP

实现步骤：

（1） 传入参数   image 单通道8位二值图  rho距离累加器，thera角度累加器，threshold 阈值，累计器用于判断一条直线的阈值  lineGap 直线上两点的最大间隙。 lineLength 直线的最小长度 linesMax 最多要找的点数。 lines 要得到的直线。
（2）根据图像大小、图像数据、建立掩码矩阵mask mask.create(height, width, CV_8UC1); 
         根据图像大小和 距离累加器rho、thera角度累加器  建立霍夫空间矩阵     accum.create(numangle, numrho, CV_32SC1);            我需要解释 图像坐标系和正常二维坐标系的不同

 二维图像都是以左上角为(0,0)，横向为x轴为width，竖向为y轴为hight，横向走的坐标值(1,0)、(2,0).... img.step[0]等于width。可以很明显的观察到，不管怎么旋转，不会有和正常坐标系重合。

（3）先计算好所需的所有正弦和余弦值，其中有一点需要解释，为什么要乘以irho，会在下面解释的

（4）初始化掩码矩阵mask，原图是边缘的位置在对应的mask矩阵中赋予1，同时把对应点加入到CV数据结构CvSeq，不知道CvSeq请看博客：https://blog.csdn.net/NewThinker_wei/article/details/45223589。不是边缘为赋予0。

（5）然后从随机的选择点，对经过的点的直线进行累加器计数，当这条直线超过设定的阈值时，搜索直线的两个端点值，关于搜索方向上（也就是⑥）会在下方做解释，当搜索端点值时，会判断点与点之间的 间隙，一定要小于等于设定的值。搜索到直线的端点值后，会判断直线长度是否大于设定的阈值。

（6）然后再次搜索端点，若满足上述的所有条件，就证明这是我们要找的一条直线，这时更新累加器霍夫矩阵和更新掩码矩阵，更新累加器霍夫矩阵是为了清除那些已经判定是好的直线上的点对应的累加器的值，避免再次利用这些累加值。若不是我们要找的直线，只更新更新掩码矩阵：搜索过的位置，不管是好的直线，还是坏的直线，掩码相应位置都清0，这样下次就不会再重复搜索这些位置了，从而达到减小计算边缘点的目的--->其实这里我是不理解的，谁能给解释一下吗

（7） 若是好的直线就加入到Lines中。

源码：

HoughLinesP函数是在sources/modules/imgproc/src/hough.cpp文件中被定义的：

void cv::HoughLinesP( InputArray _image, OutputArray _lines,
                      double rho, double theta, int threshold,
                      double minLineLength, double maxGap )
{
    Ptr<CvMemStorage> storage = cvCreateMemStorage(STORAGE_SIZE);
    Mat image = _image.getMat();
    CvMat c_image = image;
    CvSeq* seq = cvHoughLines2( &c_image, storage, CV_HOUGH_PROBABILISTIC,
                    rho, theta, threshold, minLineLength, maxGap );
    seqToMat(seq, _lines);
}

从HoughLinesP函数可以看出，该函数会调用cvHoughLines2函数。它通过参数CV_HOUGH_PROBABILISTIC，最终调用了icvHoughLinesProbabilistic函数： 

icvHoughLinesProbabilistic(CvMat* image,
	float rho, float theta, int threshold,
	int lineLength, int lineGap,
	CvSeq *lines, int linesMax)       // image 单通道8位二值图  rho距离累加器，thera角度累加器，threshold 阈值，累计器用于判断一条直线的阈值。
{                                     // lineGap 直线上两点的最大间隙。 lineLength 直线的最小长度 linesMax 最多要找的点数。
	//accum为累加器矩阵，mask为掩码矩阵
	cv::Mat accum, mask;
	std::vector<float> trigtab;       //用于存储事先计算好的正弦和余弦值
								  
	cv::MemStorage storage(cvCreateMemStorage(0));     //开辟一段内存空间
	//用于存储特征点坐标，即边缘像素的位置
	CvSeq* seq;
	CvSeqWriter writer;   
	int width, height;         //图像的宽和高
	int numangle, numrho;      //角度和距离的离散数量
	float ang;
	int r, n, count;
	CvPoint pt;
	float irho = 1 / rho;      //距离分辨率的倒数
	CvRNG rng = cvRNG(-1);     //随机数
	const float* ttab;         //向量trigtab的地址指针
	uchar* mdata0;             //矩阵mask的地址指针
	//确保输入图像的正确性 图像存在且为单通道 8位 二值图
	CV_Assert(CV_IS_MAT(image) && CV_MAT_TYPE(image->type) == CV_8UC1);      

	width = image->cols;      //提取出输入图像的宽
	height = image->rows;     //提取出输入图像的高

	//由角度和距离分辨率，得到角度和距离的离散数量 cvRound 四舍五入，近似值
	numangle = cvRound(CV_PI / theta);   
	numrho = cvRound(((width + height) * 2 + 1) / rho);
	//创建累加器矩阵，即霍夫空间
	accum.create(numangle, numrho, CV_32SC1);        //heigth：numangle width：numrho
	//创建掩码矩阵，大小与输入图像相同
	mask.create(height, width, CV_8UC1);
	//定义trigtab的大小，因为要存储正弦和余弦值，所以长度为角度离散数的2倍
	trigtab.resize(numangle * 2);
	//累加器矩阵清零
	accum = cv::Scalar(0);
	//避免重复计算，事先计算好所需的所有正弦和余弦值
	for (ang = 0, n = 0; n < numangle; ang += theta, n++)
	{
		trigtab[n * 2] = (float)(cos(ang) * irho);                  // ① 为什么要乘以irho。  
		trigtab[n * 2 + 1] = (float)(sin(ang) * irho);
	}
	//赋值首地址
	ttab = &trigtab[0];       // trigtab 存每个角度的正余弦值。
	mdata0 = mask.data;
	//开始写入序列
	cvStartWriteSeq(CV_32SC2, sizeof(CvSeq), sizeof(CvPoint), storage, &writer);

	// stage 1. collect non-zero image points
	//收集图像中的所有非零点，因为输入图像是边缘图像，所以非零点就是边缘点
	for (pt.y = 0, count = 0; pt.y < height; pt.y++)               //② mask.create(height, width, CV_8UC1); 图像坐标系的转化。 
	{
		//提取出输入图像和掩码矩阵的每行地址指针
		const uchar* data = image->data.ptr + pt.y*image->step;
		uchar* mdata = mdata0 + pt.y*width;
		for (pt.x = 0; pt.x < width; pt.x++)
		{
			if (data[pt.x])    //是边缘点
			{	
				//掩码的相应位置置1
				mdata[pt.x] = (uchar)1;         
				//把该坐标位置写入序列
				CV_WRITE_SEQ_ELEM(pt, writer);       
			}
			else    //不是边缘点
				mdata[pt.x] = 0;    //掩码的相应位置清0
		}
	}
	//终止写序列，seq为所有边缘点坐标位置的序列
	seq = cvEndWriteSeq(&writer);                                     // ③ CVseq 动态开辟空间里面的方法使用
	count = seq->total;     //得到边缘点的数量

						   // stage 2. process all the points in random order
						   //随机处理所有的边缘点
	for (; count > 0; count--)
	{
		// choose random point out of the remaining ones
		//步骤1，在剩下的边缘点中随机选择一个点，idx为不大于count的随机数
		int idx = cvRandInt(&rng) % count;                            // ④.1 为什么在剩下的点中随机取点，并且用最后一个覆盖去除
		//max_val为累加器的最大值，max_n为最大值所对应的角度
		int max_val = threshold - 1, max_n = 0;
		//由随机数idx在序列中提取出所对应的坐标点
		CvPoint* point = (CvPoint*)cvGetSeqElem(seq, idx);          //    取出点的坐标。
		//定义直线的两个端点
		CvPoint line_end[2] = { { 0,0 },{ 0,0 } };
		float a, b;
		//累加器的地址指针，也就是霍夫空间的地址指针
		int* adata = (int*)accum.data;
		int i, j, k, x0, y0, dx0, dy0, xflag;
		int good_line;
		const int shift = 16;
		//提取出坐标点的横、纵坐标
		i = point->y;
		j = point->x;

		// "remove" it by overriding it with the last element
		//用序列中的最后一个元素覆盖掉随机内存空间中的数，也是进行了删除操作了
		*point = *(CvPoint*)cvGetSeqElem(seq, count - 1);                 //④.2 覆盖去除操作。

		// check if it has been excluded already (i.e. belongs to some other line)
		//检测这个坐标点是否已经计算过，也就是它已经属于其他直线
		//因为计算过的坐标点会在掩码矩阵mask的相对应位置清零
		// 掩码矩阵，为1则为边缘点，为0则为不为边缘点或被处理过的点。

		if (!mdata0[i*width + j])    //该坐标点被处理过                     
			continue;    //不做任何处理，继续主循环

		// update accumulator, find the most probable line
		//步骤2，更新累加器矩阵，找到最有可能的直线
		// numrho： 霍夫空间的每一行的大小，也就是霍夫矩阵宽。
		for (n = 0; n < numangle; n++, adata += numrho)               
		{
			//由角度计算距离
			r = cvRound(j * ttab[n * 2] + i * ttab[n * 2 + 1]);      // 固定点不变，变换角度，计算ρ
			r += (numrho - 1) / 2;                                   // 为什么这样安排呢?  可能是因为θ、ρ正负问题

			//在累加器矩阵的相应位置上数值加1，并赋值给val
			int val = ++adata[r];
			//更新最大值，并得到它的角度
			if (max_val < val)
			{
				max_val = val;
				max_n = n;                     //记录角 找到点最多的直线。
			}
		}

		// if it is too "weak" candidate, continue with another point
		//步骤3，如果上面得到的最大值小于阈值，则放弃该点，继续下一个点的计算
		if (max_val < threshold)                                   //  ⑤当前点最多的直线的点数若大于阈值了，再做进一步判断。       
			continue;

		// from the current point walk in each direction
		// along the found line and extract the line segment
		//步骤4，从当前点出发，沿着它所在直线的方向前进，直到达到端点为止
		a = -ttab[max_n * 2 + 1];    //a=-sinθ
		b = ttab[max_n * 2];         //b=cosθ
								     //当前点的横、纵坐标值
		x0 = j;
		y0 = i;
		//θ为极坐标系中的角度，真正的直线斜率为 tan(π/2+θ) = -cotθ = - cosθ/sinθ
		//当θ在 45度~135度和 225度~315度 之间时，才满足fabs(-sinθ) > fabs(cosθ) 而真正的直线角度 也就是在0～45或135度～180度之间
		if (fabs(a) > fabs(b))                          //直线在0～45或135度～180度之间
		{
			xflag = 1;                                 //置标识位，标识直线的粗略方向
			//确定横、纵坐标的位移量			                                         
			dx0 = a > 0 ? 1 : -1; 
			// 先扩大2^16倍，无精度损失的计算，后面又减少2^16倍，恢复原值。
			dy0 = cvRound(b*(1 << shift) / fabs(a));                 // ⑥ 为什么x，y位移量要这样设置呢。
			y0 = (y0 << shift) + (1 << (shift - 1));
		}
		else                                            //直线在45度~135度
		{
			xflag = 0;   //清标识位
						 //确定横、纵坐标的位移量
			dy0 = b > 0 ? 1 : -1;
			dx0 = cvRound(a*(1 << shift) / fabs(b));
			//确定横坐标
			x0 = (x0 << shift) + (1 << (shift - 1));
		}
		//搜索直线的两个端点
		for (k = 0; k < 2; k++)
		{
			//gap表示两条直线的间隙，x和y为搜索位置，dx和dy为位移量
			int gap = 0, x = x0, y = y0, dx = dx0, dy = dy0;
			//搜索第二个端点的时候，反方向位移
			if (k > 0)
				dx = -dx, dy = -dy;

			// walk along the line using fixed-point arithmetics,
			// stop at the image border or in case of too big gap
			//沿着直线的方向位移，直到到达图像的边界或大的间隙为止
			for (;; x += dx, y += dy)
			{
				uchar* mdata;
				int i1, j1;
				//确定新的位移后的坐标位置
				if (xflag)
				{
					j1 = x;
					i1 = y >> shift;
				}
				else
				{
					j1 = x >> shift;
					i1 = y;
				}
				//如果到达了图像的边界，停止位移，退出循环
				if (j1 < 0 || j1 >= width || i1 < 0 || i1 >= height)
					break;
				//定位位移后掩码矩阵位置
				mdata = mdata0 + i1*width + j1;
				// for each non-zero point:
				//    update line end,
				//    clear the mask element
				//    reset the gap
				//该掩码不为0，说明该点可能是在直线上
				if (*mdata)
				{
					gap = 0;    //设置间隙为0
								//更新直线的端点位置
					line_end[k].y = i1;
					line_end[k].x = j1;
				}
				//掩码为0，说明不是直线，但仍继续位移，直到间隙大于所设置的阈值为止
				//⑦ 点与点之间的 间隙，一定要小于等于设定的值。
				else if (++gap > lineGap)    //间隙加1               
					break;
			}
		}
		//步骤5，由检测到的直线的两个端点粗略计算直线的长度
		//当直线长度大于所设置的阈值时，good_line为1，否则为0
		//⑧ 同时直线长度也要满足设定的阈值。
		good_line = abs(line_end[1].x - line_end[0].x) >= lineLength ||
			abs(line_end[1].y - line_end[0].y) >= lineLength;
		//再次搜索端点，目的是更新累加器矩阵和更新掩码矩阵，以备下一次循环使用
		for (k = 0; k < 2; k++)
		{
			int x = x0, y = y0, dx = dx0, dy = dy0;

			if (k > 0)
				dx = -dx, dy = -dy;

			// walk along the line using fixed-point arithmetics,
			// stop at the image border or in case of too big gap
			for (;; x += dx, y += dy)
			{
				uchar* mdata;
				int i1, j1;

				if (xflag)
				{
					j1 = x;
					i1 = y >> shift;               
				}
				else
				{
					j1 = x >> shift;
					i1 = y;
				}

				mdata = mdata0 + i1*width + j1;
				// for each non-zero point:
				//    update line end,
				//    clear the mask element
				//    reset the gap
				if (*mdata)
				{
					//⑨ if语句的作用是清除那些已经判定是好的直线上的点对应的累加器的值，避免再次利用这些累加值
					if (good_line)    //在第一次搜索中已经确定是好的直线
					{
						//得到累加器矩阵地址指针
						adata = (int*)accum.data;
						for (n = 0; n < numangle; n++, adata += numrho)
						{
							r = cvRound(j1 * ttab[n * 2] + i1 * ttab[n * 2 + 1]);
							r += (numrho - 1) / 2;
							adata[r]--;    //相应的累加器减1
						}
					}
					//⑩ 搜索过的位置，不管是好的直线，还是坏的直线，掩码相应位置都清0，这样下次就不会再重复搜索这些位置了，从而达到减小计算边缘点的目的
					*mdata = 0;
				}
				//如果已经到达了直线的端点，则退出循环
				if (i1 == line_end[k].y && j1 == line_end[k].x)
					break;
			}
		}
		//如果是好的直线
		if (good_line)
		{
			CvRect lr = { line_end[0].x, line_end[0].y, line_end[1].x, line_end[1].y };
			//把两个端点压入序列中
			cvSeqPush(lines, &lr);
			//如果检测到的直线数量大于阈值，则退出该函数
			if (lines->total >= linesMax)
				return;
		}
	}
}

其实代码中，已经注释好多了，但是我现在还是要对一些地方进行解释：

 ①的地方： 为什么要乘以irho。

irho 等于距离累加器rho的导数，可以看到就那个地方用到了irho了，当存所有角度的sin 和 cos值的时候：

trigtab[n * 2] = (float)(cos(ang) * irho);     trigtab[n * 2 + 1] = (float)(sin(ang) * irho);

还有计算ρ的时候，r = cvRound(j * ttab[n * 2] + i * ttab[n * 2 + 1]); 

其实为了霍夫空间图像坐标(θ,ρ),中ρ和实际坐标的转化，因为ρ的取值范围根据 numrho = cvRound(((width + height) * 2 + 1) / rho); 这个确定的，除以rho是把ρ和实际坐标之间缩小或放大多少倍。

当计算 r = cvRound(j * ttab[n * 2] + i * ttab[n * 2 + 1]);  在这个公式中(i,j) 是 实际矩阵中的坐标，r是霍夫空间中的ρ，在这转化时，就用到了*irho。

④随机取点：
⑥ r = cvRound(j * ttab[n * 2] + i * ttab[n * 2 + 1]);  当我们看到这个公式时，是不是会想到(i,j)点图像那个位置，直线(θ,ρ)  θ角这条直线又是在哪，当时我想了很长时间也没有相同，其实换一种思路，把图像坐标系中的点，对应到正常X-Y轴坐标系上的点就可以了。

θ为极坐标系中的角度，真正的直线斜率为 tan(π/2+θ) = -cotθ = - cosθ/sinθ

当θ在 45度~135度和 225度~315度 之间时，才满足fabs(-sinθ) > fabs(cosθ) 而真正的直线角度 也就是在0～45或135度～180度之间

关于 直线角度 在0～45或135度～180度之间，为什么是  dx0 = a > 0 ? 1 : -1;    dy0 = cvRound(b*(1 << shift) / fabs(a)); 

如图：当直线角度为0~45度时，x变化，要比y变化快，因为直线斜率等于 tan(π/2+θ) = -cotθ = - cosθ/sinθ = △y/△x，当△x等于1时，△y一定小于1，等于 - cosθ/sinθ，也是代码所示，关于正负号，大家自己想吧。同理，直线角度在45度~135度也是同理分析的。