RANSAC估计——以直线拟合为例

最新推荐文章于 2022-06-19 22:31:48 发布

XDWX

最新推荐文章于 2022-06-19 22:31:48 发布

阅读量2.2k

点赞数 4

分类专栏： RANSAC

RANSAC 专栏收录该内容

9 篇文章 0 订阅

订阅专栏

RANSAC(RANdom SAmple Consensus)，即随机采样一致性。该方法最早是由Fischler和Bolles提出的一种鲁棒估计方法，最早用于计算机视觉中位姿估计问题，现在已广泛应用于已知模型的参数估计问题中。对于数据中存在大比例外点（错误数据、野值点）时，RANSAC方法十分有效。下面就以直线估计的例子来说明RANSAC的基本思想。

直线拟合RANSAC估计基本原理

RANSAC的思想比较简单，主要有以下几步：

随机选择两点(确定一条直线所需要的最小点集)；由这两个点确定一条线l；
使用几何距离，求最大一致集的最佳拟合直线，作为数据点的最佳匹配直线。

如下图所示，首先在点集中随机选择两个点，求解它们构成的直线参数，再计算点集中剩余点到该直线的距离，距离小于设置的距离阈值的点为内点，统计内点个数；接下来再随机选取两个点，同样统计内点个数，以此类推；其中内点个数最多的点集即为最大一致集，最后将该最大一致集里面的点利用最小二乘拟合出一条直线。
这里写图片描述

参数设置

根据上面RANSAC基本原理的介绍，在这算法流程中存在两个重要的参数需要设置，采样次数和距离阈值。

为了保证算法的效率，当数据量比较大时，我们不可能采集所有的样本，将所有的可能情况都计算一遍。如何在抽样当中保证至少有一个好样本（不包含外点）是我们需要考虑的问题。Fischler和Bolles在原文中通过概率统计的方法得出了采样次数与数据中外点比例和得到一个好样本概率之间的数学关系，在这里我直接给出结论，有兴趣推导该关系的博友可以自己阅读原文或者由吴福朝编写的《计算机视觉中的数学方法》。

$K = log ( 1 - z ) log ( 1 - w n )$

其中，K为模型参数估计需要的最小点个数，直线拟合最少需要2个点。

距离阈值一般需要根据经验来设定，在直线拟合实验中，我设置的参数为2.99。当观测误差符合0均值和sigma。

RANSAC直线拟合代码

结合OpenCV库（2.0版本以上）实现了RANSAC直线拟合，在主函数中给出了一个例子，仅供大家参考。

#include <opencv2\opencv.hpp>
#include <iostream>
#include <ctime>

using namespace std;
using namespace cv;

//生成[0,1]之间符合均匀分布的数
double uniformRandom(void)
{
    return (double)rand() / (double)RAND_MAX;
}

//生成[0,1]之间符合高斯分布的数
double gaussianRandom(void)
{
    /* This Gaussian routine is stolen from Numerical Recipes and is their
    copyright. */
    static int next_gaussian = 0;
    static double saved_gaussian_value;

    double fac, rsq, v1, v2;

    if (next_gaussian == 0) {
        do {
            v1 = 2 * uniformRandom() - 1;
            v2 = 2 * uniformRandom() - 1;
            rsq = v1*v1 + v2*v2;
        } while (rsq >= 1.0 || rsq == 0.0);
        fac = sqrt(-2 * log(rsq) / rsq);
        saved_gaussian_value = v1*fac;
        next_gaussian = 1;
        return v2*fac;
    }
    else {
        next_gaussian = 0;
        return saved_gaussian_value;
    }
}

//根据点集拟合直线ax+by+c=0，res为残差
void calcLinePara(vector<Point2d> pts, double &a, double &b, double &c, double &res)
{
    res = 0;
    Vec4f line;
    vector<Point2f> ptsF;
    for (unsigned int i = 0; i < pts.size(); i++)
        ptsF.push_back(pts[i]);

    fitLine(ptsF, line, CV_DIST_L2, 0, 1e-2, 1e-2);
    a = line[1];
    b = -line[0];
    c = line[0] * line[3] - line[1] * line[2];

    for (unsigned int i = 0; i < pts.size(); i++)
    {
        double resid_ = fabs(pts[i].x * a + pts[i].y * b + c);
        res += resid_;
    }
    res /= pts.size();
}

//得到直线拟合样本，即在直线采样点集上随机选2个点
bool getSample(vector<int> set, vector<int> &sset)
{
    int i[2];
    if (set.size() > 2)
    {
        do
        {
            for (int n = 0; n < 2; n++)
                i[n] = int(uniformRandom() * (set.size() - 1));
        } while (!(i[1] != i[0]));
        for (int n = 0; n < 2; n++)
        {
            sset.push_back(i[n]);
        }
    }
    else
    {
        return false;
    }
    return true;
}

//直线样本中两随机点位置不能太近
bool verifyComposition(const vector<Point2d> pts)
{
    cv::Point2d pt1 = pts[0];
    cv::Point2d pt2 = pts[1];
    if (abs(pt1.x - pt2.x) < 5 && abs(pt1.y - pt2.y) < 5)
        return false;

    return true;
}

//RANSAC直线拟合
void fitLineRANSAC(vector<Point2d> ptSet, double &a, double &b, double &c, vector<bool> &inlierFlag)
{
    double residual_error = 2.99; //内点阈值

    bool stop_loop = false;
    int maximum = 0;  //最大内点数

    //最终内点标识及其残差
    inlierFlag = vector<bool>(ptSet.size(), false);
    vector<double> resids_(ptSet.size(), 3);
    int sample_count = 0;
    int N = 500;

    double res = 0;

    // RANSAC
    srand((unsigned int)time(NULL)); //设置随机数种子
    vector<int> ptsID;
    for (unsigned int i = 0; i < ptSet.size(); i++)
        ptsID.push_back(i);
    while (N > sample_count && !stop_loop)
    {
        vector<bool> inlierstemp;
        vector<double> residualstemp;
        vector<int> ptss;
        int inlier_count = 0;
        if (!getSample(ptsID, ptss))
        {
            stop_loop = true;
            continue;
        }

        vector<Point2d> pt_sam;
        pt_sam.push_back(ptSet[ptss[0]]);
        pt_sam.push_back(ptSet[ptss[1]]);

        if (!verifyComposition(pt_sam))
        {
            ++sample_count;
            continue;
        }

        // 计算直线方程
        calcLinePara(pt_sam, a, b, c, res);
        //内点检验
        for (unsigned int i = 0; i < ptSet.size(); i++)
        {
            Point2d pt = ptSet[i];
            double resid_ = fabs(pt.x * a + pt.y * b + c);
            residualstemp.push_back(resid_);
            inlierstemp.push_back(false);
            if (resid_ < residual_error)
            {
                ++inlier_count;
                inlierstemp[i] = true;
            }
        }
        // 找到最佳拟合直线
        if (inlier_count >= maximum)
        {
            maximum = inlier_count;
            resids_ = residualstemp;
            inlierFlag = inlierstemp;
        }
        // 更新RANSAC迭代次数，以及内点概率
        if (inlier_count == 0)
        {
            N = 500;
        }
        else
        {
            double epsilon = 1.0 - double(inlier_count) / (double)ptSet.size(); //野值点比例
            double p = 0.99; //所有样本中存在1个好样本的概率
            double s = 2.0;
            N = int(log(1.0 - p) / log(1.0 - pow((1.0 - epsilon), s)));
        }
        ++sample_count;
    }

    //利用所有内点重新拟合直线
    vector<Point2d> pset;
    for (unsigned int i = 0; i < ptSet.size(); i++)
    {
        if (inlierFlag[i])
            pset.push_back(ptSet[i]);
    }

    calcLinePara(pset, a, b, c, res);
}

void main()
{
    //demo
    int width = 640;
    int height = 320;

    //直线参数
    double a = 1, b = 2, c = -640;

    //随机获取直线上20个点
    int ninliers = 0;
    vector<Point2d> ptSet;
    srand((unsigned int)time(NULL)); //设置随机数种子
    while (true)
    {
        double x = uniformRandom()*(width - 1);
        double y = -(a*x + c) / b;
        //加0.5高斯噪声
        x += gaussianRandom()*0.5;
        y += gaussianRandom()*0.5;
        if (x >= 640 && y >= 320)
            continue;
        Point2d pt(x, y);
        ptSet.push_back(pt);
        ninliers++;
        if (ninliers == 20)
            break;
    }

    int nOutliers = 0;
    //随机获取10个野值点
    while (true)
    {
        double x = uniformRandom() * (width - 1);
        double y = uniformRandom() * (height - 1);

        if (fabs(a*x + b*y + c) < 10)  //野值点到直线距离不小于10个像素
            continue;

        Point2d pt(x, y);
        ptSet.push_back(pt);
        nOutliers++;
        if (nOutliers == 10)
            break;
    }

    //绘制点集中所有点
    Mat img(321, 641, CV_8UC3, Scalar(255, 255, 255));
    for (unsigned int i = 0; i < ptSet.size(); i++)
        circle(img, ptSet[i], 3, Scalar(255, 0, 0), 3, 8);

    double A, B, C;
    vector<bool> inliers;
    fitLineRANSAC(ptSet, A, B, C, inliers);

    B = B / A;
    C = C / A;
    A = A / A;

    //绘制直线
    Point2d ptStart, ptEnd;
    ptStart.x = 0;
    ptStart.y = -(A*ptStart.x + C) / B;
    ptEnd.x = -(B*ptEnd.y + C) / A;
    ptEnd.y = 0;
    line(img, ptStart, ptEnd, Scalar(0, 255, 255), 2, 8);
    cout << "A:" << A << " " << "B:" << B << " " << "C:" << C << " " << endl;

    imshow("line fitting", img);
    waitKey();
}
 
 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258