RANSAC LMedS 详细分析

RANSAC的经典步骤是

1. 随机采样

2. 计算模型参数

3. 阈值分类inlier outlier

4. 迭代2. 3.步直至获得符合阈值的最优解

5. 精确优化模型参数（LM算法迭代优化）

RANSAC的函数接口 参照opencv来说主要需要3-4个参数（第四个不是必须的）

1. 误差阈值ransacThreshold：区分inlier和outliner的依据

2. 置信度confidence：设置之后代表RANSAC采样n次过程中会出现（至少一次）采样点数据集中的点都为内点的概率

这个值设置的太大，会增加采样次数。太小，会使结果不太理想。一般取0.95-0.99

3. 最大采样迭代次数maxIters：为了防止一直在采样计算

4.*最大精细迭代次数refineIters：在采样之后，选取最优解。可以增加精确优化，比如使用LM算法获得更优解。

以上4个参数，有三个是经验值。其中最大采样迭代次数maxIters是可以有数值解的。

在RANSAC计算中，在已经采样计算之后，可以根据inlier num 和 outlier num求得一个更准确的maxIters

这里先列出一个等式：

1 − p = (1 − wⁿ)^{k                           （1）}

其中

p 表示置信度confidence

w 表示数据集中outlier占的比例

n 表示采样点数

k 表示需要迭代采样的最少次数

1 − wⁿ 表示采样一次，n个点中至少有一个outlier的概率

(1 − wⁿ)^k 表示采样k次，n个点中至少有一个outlier的概率

因为p为采样k次，能有至少一次n个点都是inlier的概率

所以(1 − p)和(1 − wⁿ)^k相等时，k 为需要迭代采样的最少次数

下面是k的解析解：

                                （2）

关于这个部分的实现 opencv中代码如下

        for( iter = 0; iter < niters; iter++ )
        {
            int i, nmodels;
            if( count > modelPoints )
            {
                bool found = getSubset( m1, m2, ms1, ms2, rng, 10000 );
                if( !found )
                {
                    if( iter == 0 )
                        return false;
                    break;
                }
            }

            nmodels = cb->runKernel( ms1, ms2, model );
            if( nmodels <= 0 )
                continue;
            CV_Assert( model.rows % nmodels == 0 );
            Size modelSize(model.cols, model.rows/nmodels);

            for( i = 0; i < nmodels; i++ )
            {
                Mat model_i = model.rowRange( i*modelSize.height, (i+1)*modelSize.height );
                int goodCount = findInliers( m1, m2, model_i, err, mask, threshold );

                if( goodCount > MAX(maxGoodCount, modelPoints-1) )
                {
                    std::swap(mask, bestMask);
                    model_i.copyTo(bestModel);
                    maxGoodCount = goodCount;
                    niters = RANSACUpdateNumIters( confidence, (double)(count - goodCount)/count, modelPoints, niters );
                }
            }
        }

可以看到在采样迭代的环节中，最大迭代次数niters是动态更新的。在循环内，会被修改。

int RANSACUpdateNumIters( double p, double ep, int modelPoints, int maxIters )
{
    if( modelPoints <= 0 )
        CV_Error( Error::StsOutOfRange, "the number of model points should be positive" );

    p = MAX(p, 0.);
    p = MIN(p, 1.);
    ep = MAX(ep, 0.);
    ep = MIN(ep, 1.);

    // avoid inf's & nan's
    double num = MAX(1. - p, DBL_MIN);
    double denom = 1. - std::pow(1. - ep, modelPoints);
    if( denom < DBL_MIN )
        return 0;

    num = std::log(num);
    denom = std::log(denom);

    return denom >= 0 || -num >= maxIters*(-denom) ? maxIters : cvRound(num/denom);
}

这里的实现也是按照上面的公式（2）书写的。

末尾增加了一个功能，当求得的迭代次数比之前的大，就更新迭代次数。

LMedS的经典步骤是

1. 随机采样

2. 计算模型参数

3. 计算相对模型的点集偏差err，并求出偏差中值Med(err)

4. 迭代2. 3.步直至获得符合阈值的最优解：Med(err)最小

5. 精确优化模型参数（LM算法迭代优化）

LMedS的函数接口 参照opencv来说主要需要2-3个参数（第三个不是必须的）

1. 置信度confidence：设置之后代表RANSAC采样n次过程中会出现（至少一次）采样点数据集中的点都为内点的概率

这个值设置的太大，会增加采样次数。太小，会使结果不太理想。一般取0.95-0.99

2. 最大采样迭代次数maxIters：为了防止一直在采样计算

3.*最大精细迭代次数refineIters：在采样之后，选取最优解。可以增加精确优化，比如使用LM算法获得更优解。

注意： 相对于RANSAC，LMedS有一个优点：不需要指定 - 误差阈值ransacThreshold：区分inlier和outliner的依据

RANSAC的阈值在具有物理意义或者几何意义的时候比较容易确定，

但是当阈值不具有这些特征的时候，就成了一个不太好调整的参数了。

这时LMedS可以自适应迭代获得最优解。

此外，LMedS也能自适应获得inlier和outliner

这里要引用 Robust Regression Methods for Computer Vision_A Review【1】

这里先列出一个等式：

其中

n 表示点集的个数

p 表示计算模型一次采样的点个数

ri²  表示误差

med(ri² ) 表示误差中值

筛选条件为

关于这个部分的实现 opencv中代码如下

            double sigma = 2.5*1.4826*(1 + 5./(count - modelPoints))*std::sqrt(minMedian);
            sigma = MAX( sigma, 0.001 );

            count = findInliers( m1, m2, bestModel, err, mask, sigma );

findInliers的计算如下

    int findInliers( const Mat& m1, const Mat& m2, const Mat& model, Mat& err, Mat& mask, double thresh ) const
    {
        cb->computeError( m1, m2, model, err );
        mask.create(err.size(), CV_8U);

        CV_Assert( err.isContinuous() && err.type() == CV_32F && mask.isContinuous() && mask.type() == CV_8U);
        const float* errptr = err.ptr<float>();
        uchar* maskptr = mask.ptr<uchar>();
        float t = (float)(thresh*thresh);
        int i, n = (int)err.total(), nz = 0;
        for( i = 0; i < n; i++ )
        {
            int f = errptr[i] <= t;
            maskptr[i] = (uchar)f;
            nz += f;
        }
        return nz;
    }

可以看到 里面用的是二阶范数，sigma输入到thresh参数入口之后会和err的计算一样，使用二阶范数来衡量

由于LMedS会需要对整个点集的err求中值，当点集很大的时候，求中值的过程会很消耗时间

文【1】中提到可以使用Monte Carlo方法，利用局部采样来模拟全局。

这里的式子和上面的式（1）是一样的 ：1 − p = (1 − wⁿ)^k   

使用多次采样，得到至少一次采样全是inlier的最少采样次数k。详细过程可以看论文。

【1】Robust Regression Methods for Computer Vision_A Review