直线拟合原理详解

前言

在平时的工作中，经常使用Opencv中的fitLine()函数，解决一些工业视觉检测的问题。如图所示：为了检测出玻璃切割上下磨损区域的宽度，可以先对其进行直线拟合，再提取几个检测点，计算平均距离，作为最终的检测宽度。

对于这种宽度的磨损区宽度的检测，通常分为一下两个步骤：

1. 如何在图像上提取出用于拟合的特征点；
2. 根据特征点计算出直线所对应的参数。
   拟合的特征点我们可以通过最大梯度变化的位置，进行定位。当用于拟合的点与真实直线上的点误差较小时，使用DIST_L2参数，既最小二乘法，即可得到不错的结果。但对于上图场景来说，磨损区域的纹理和拍照的阴影区域的变化较大，有可能会导致出现一些误差较大的定位点。这时，选择DIST_HUBER参数，则能得到较DIST_L2更为准确的结果。Opencv官方文档中，关于fitLine()参数的介绍如下图所示：
   
   文档中只给出了简单的原理性说明，对于喜欢搞清详细原理的博主来说，还是远远不够的，因此便做了一点研究，写下此文，便于总结分享。

一、极大似然估计

fitLine()拟合直线原理，其文档中原理说明如下：

该算法是基于M-estimation(最大似然估计)实现的，即使用加权的最小二乘算法进行迭代拟合。每次迭代后，权重wi随误差和p(ri)成反比例调整。
是不是感觉跟神经网络中每个感知机的权重调节方式是类似啊呀！
接下来，我们根据这个原理，依次讲解相关的知识点

最大似然估计（MLE），是一种统计方法，用于估计一个概率模型的参数。其基本思想是：在已知某个参数能使样本出现的概率最大时，选择这个参数作为估计的真实值。直观地说，如果一个参数使得观测到的数据出现的概率最大，那么这个参数就是最可能的真实值

下文中 对数似然函数 和 最小二乘法 的讲解，均引用自该博文：(https://blog.youkuaiyun.com/afgc223/article/details/133681873)


使用最小二乘法求解的方法网上资源很多，这里就不进行说明了。

二、迭代拟合

第一次最小二乘计算，我们使对数似然函数最大的公式为：

即，给每个拟合点赋予了相同的权重。由于可能可能出现个别拟合点偏离直线距离较大，导致此时得到的结果并不一定是准确的，因此，我们需要给不同的拟合点赋予不同的拟合权重，进行迭代计算。
每个拟合点的权重可以与其残差p(ri)成反比例调整。比如，当 fitLine()拟合参数选择DIST_L2时，残差为

r 为拟合点到本次计算出的直线的距离。根据残差，计算出每个拟合点的权重a1,a2,…,an，得到新的使对数似然函数最大的公式：

据此便可进行第二次的最小二乘计算。
经过数次迭代，当两次迭代结果之差满足设定的阈值时，迭代结束，得到最终结果

至此，当拟合点中有很多误差较大的点时，我们对fitLine()函数选择拟合参数DIST_HUBER 优于 DIST_L2 的原因，大概可以有点感觉了。因为使用的残差计算方法不同，导致权重的分配也不相同。残差计算方法的原因还需进一步详细了解

三、RANSAN 直线拟合算法

参考文章 ： https://blog.youkuaiyun.com/chinaswin/article/details/68066295
直接使用fitLine()函数进行直线拟合，一般只适用于噪声点较少的场景，当噪声点很多时，如图所示：

此时，就需要使用 RANSAC（RANdom SAmple Consensus,RANSAC）算法，即 随机抽样一致算法。该算法假设数据中包含正确数据和异常数据（或称为噪声）。正确数据记为内点（inliers），异常数据记为外点（outliers）。得到一系列内点后，再进行 直线拟合，即可减小过多噪声的影响。

算法基本思想

如下图所示，存在很多离散的点，而我们认为这些点构成了一条直线。当然，人眼能很清晰地拟合出这条直线，找到外点。但要让计算机找到这条直线，在很久之前是很难的，RACSAC的出现是通过数学之美解决这一难题的重要发明。

算法步骤：

（1）随机选择两个点，用于拟合直线方程（一条直线至少需要两个点确定）
（2）通过这两个点，我们可以计算出这两个点所表示的直线方程 y=ax+b
（3）将所有的数据点带入这个直线方程计算误差
（4）找到所有满足误差阈值的点（第4幅图中可以看到，只有很少的点支持这个模型）
（5）重复（1）~（4）这个过程，直到达到一定迭代次数后，我们选出那个被支持的最多的模型，作为问题的解。如下图所示：

C++代码实现

cv::Vec4f AlgKernel::fitLineRansac(const std::vector<cv::Point2d>& vecPts, DistanceTypes fitlinedistype, std::vector<cv::Point2d>& inlierPts, const double & residErr)
{
	Vec4f line;
	int numPts = (int)vecPts.size();
	if (numPts < 2) return Vec4f();

	int iterCnt = 50;
	const int minCnt = 2;								//拟合直线最少抽样个数
	int nIter = 0;										//统计迭代次数
	int nInlierMax = 0;									//保存内点数量最大值

	double modelRes = 0;								//模型残差初始化
	vector<bool> inlierFlag = vector<bool>(numPts, false);	//内点标识初始化

	//step1 随机生成点的下标
	static default_random_engine rng;
	static uniform_int_distribution<unsigned> uniform;
	decltype(uniform.param()) new_range(0, numPts - 1);

	uniform.param(new_range);					//设定随机数范围
	rng.seed(1);								//初始化
	std::set<unsigned int> selectIndexs;		//选择点的下标，自动排序
	vector<Point2f> selectPts;					//选择点
	double A = 0.0, B = 0.0, C = 0.0;

	//step2 使用Ransac提取内点
	while (nIter < iterCnt) {	//执行迭代
		selectIndexs.clear();
		selectPts.clear();

		int inlierCnt = 0;						//统计当前迭代中内点个数
		vector<bool> inliersTemp;				//临时内点标识

		//1.随机选取minCnt个点
		while (1) {
			unsigned int index = uniform(rng);
			selectIndexs.insert(index);
			if (selectIndexs.size() == minCnt) { break; }
		}

		//2.获取抽样点
		std::set<unsigned int>::iterator selectItc = selectIndexs.begin();
		while (selectItc != selectIndexs.end()) {
			unsigned int index = *selectItc;
			selectPts.push_back(vecPts[index]);
			selectItc++;
		}

		//3.根据残差统计内点个数
		calcLinePara(selectPts, A, B, C, modelRes);
		for (unsigned int i = 0; i < numPts; i++) {
			Point2f pt = vecPts[i];
			double resid_ = fabs(pt.x * A + pt.y * B + C);
			//residualsTemp.push_back(resid_);
			inliersTemp.push_back(false);
			if (resid_ < residErr) {		//阈值范围内
				++inlierCnt;
				inliersTemp[i] = true;
			}
		}

		//4.根据内点个数更新结果 
		if (inlierCnt >= nInlierMax) {
			nInlierMax = inlierCnt;
			//resids_ = residualsTemp;
			inlierFlag = inliersTemp;
		}
		else if (inlierCnt == 0) {			//此时迭代次数趋于无穷大
			iterCnt = 500;
		}
		else {	//内点数量减少，更新迭代次数
			double inlierRatio = double(inlierCnt) / numPts;
			double p = 0.99;
			double s = 2.0;
			iterCnt = int(std::log(1.0 - p) / std::log(1.0 - pow(inlierRatio, s))); //计算迭代次数
		}
		++nIter;
	}//end_while

	 //step3 对内点进行最小二乘法拟合
	inlierPts.clear();
	for (int i = 0; i < numPts; ++i)
	{
		if (inlierFlag[i])
			inlierPts.push_back(vecPts[i]);
	}

	if (inlierPts.size() < 2) return Vec4f();
	cv::fitLine(inlierPts, line, fitlinedistype, 0, 0.01, 0.01);

	return line;
}

本文主要是博主的自我总结与分享，如有不当之处，还望帮忙指出，万分感谢！