直线拟合原理详解：从极大似然估计到RANSAC算法

直线拟合原理详解：从极大似然估计到RANSAC算法

引用

CSDN

1.

https://blog.csdn.net/weixin_38584764/article/details/143563622

在计算机视觉和图像处理领域，直线拟合是一个常见的问题，特别是在工业检测、目标识别等场景中。本文将详细介绍OpenCV库中fitLine()函数的原理和实现方法，包括极大似然估计、迭代拟合和RANSAC算法等核心概念，并提供具体的C++代码实现。

前言

在工业视觉检测中，直线拟合是一个常见的需求。例如，为了检测玻璃切割区域的宽度，可以先对图像进行直线拟合，然后提取几个检测点，计算平均距离，作为最终的检测宽度。

直线拟合通常分为两个步骤：

1. 在图像上提取用于拟合的特征点
2. 根据特征点计算出直线所对应的参数

拟合的特征点可以通过最大梯度变化的位置来定位。当用于拟合的点与真实直线上的点误差较小时，使用最小二乘法（DIST_L2）可以得到不错的结果。但对于纹理和光照变化较大的场景，选择DIST_HUBER参数会得到更准确的结果。

OpenCV官方文档中关于fitLine()参数的介绍如下：

文档中只给出了简单的原理性说明，本文将深入探讨其详细原理。

一、极大似然估计

fitLine()拟合直线的原理基于M-estimation（最大似然估计），即使用加权的最小二乘算法进行迭代拟合。每次迭代后，权重wi随误差和p(ri)成反比例调整。

最大似然估计（MLE）是一种统计方法，用于估计一个概率模型的参数。其基本思想是：在已知某个参数能使样本出现的概率最大时，选择这个参数作为估计的真实值。

接下来，我们根据这个原理，依次讲解相关的知识点。

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

二、迭代拟合

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

即，给每个拟合点赋予了相同的权重。由于可能可能出现个别拟合点偏离直线距离较大，导致此时得到的结果并不一定是准确的，因此，我们需要给不同的拟合点赋予不同的拟合权重，进行迭代计算。

每个拟合点的权重可以与其残差p(ri)成反比例调整。比如，当 fitLine()拟合参数选择DIST_L2时，残差为

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

据此便可进行第二次的最小二乘计算。

经过数次迭代，当两次迭代结果之差满足设定的阈值时，迭代结束，得到最终结果。

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

三、RANSAC 直线拟合算法

直接使用fitLine()函数进行直线拟合，一般只适用于噪声点较少的场景，当噪声点很多时，如图所示：

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

算法基本思想

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

算法步骤：

1. 随机选择两个点，用于拟合直线方程（一条直线至少需要两个点确定）
2. 通过这两个点，我们可以计算出这两个点所表示的直线方程 y=ax+b
3. 将所有的数据点带入这个直线方程计算误差
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;
}

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