RANSAC是一种常用的剔除数据中异常点的方法。本文以拟合圆为例展示RANSAC的工作方式。首先我们有一组点，假设内点的概率是p。我们要使RANSAC的成功率至少达到${ \eta=99.9\% }$，那么至少需要重复选择多少次样本？首先确定一个圆需要3个点，即每次选择随机选3个点，因此一次选择包含外点的概率是${ 1-p^{3} }$。那么N次选择中每次都包含外点的概率是${ \left( 1-p^{3} \right)^{N} }$。我们需要N次选择至少有一次全是内点的概率超过${ \eta }$，有：
$${ 1-\left( 1-p^{3} \right)^{N} \geqslant \eta }$$
解之得：
$${ N \geqslant \frac{ ln\left( 1-\eta \right) }{ ln\left( 1-p^{3} \right) } }$$
下面给出代码：
class Ransac { public: Ransac(int ifitNumber, float confidence = 0.995f); virtual ~Ransac() = default; int execute(vector<Point2f>& points); protected: virtual void fitModel(const vector<Point2f>& selects, Mat& result) = 0; virtual float evalModel(const Mat& model, const Point2f& pt) = 0; private: int fitNumber; int loopCount; }; Ransac::Ransac(int ifitNumber, float confidence) { fitNumber = ifitNumber; loopCount = logf(1 - confidence) / logf(1 - powf(0.7f, fitNumber)); } int Ransac::execute(vector<Point2f>& points) { std::random_device rd; std::mt19937 e(rd()); int count = (int)points.size(); int bestInner = 0; vector<Point2f> innerDots, outerDots; for (int i = 0; i < loopCount; i++) { for (int j = 0; j < fitNumber; j++) { std::uniform_int_distribution<int> dist(j, count - 1); int r = dist(e); std::swap(points[j], points[r]); } Mat result; fitModel(points, result); int inner = 0; vector<Point2f> dots, outs; dots.reserve(count); outs.reserve(count); for (auto& item : points) { float score = evalModel(result, item); if (score >= 0.5f) { inner++; dots.push_back(item); } else { outs.push_back(item); } } if (inner > bestInner) { bestInner = inner; innerDots = std::move(dots); outerDots = std::move(outs); } } innerDots.insert(innerDots.end(), outerDots.begin(), outerDots.end()); points = std::move(innerDots); return bestInner; } class CircleRansac : public Ransac { public: CircleRansac(float ierror, float confidence = 0.995f); protected: void fitModel(const vector<Point2f>& selects, Mat& result) override; float evalModel(const Mat& model, const Point2f& pt) override; private: float error; }; CircleRansac::CircleRansac(float ierror, float confidence) : Ransac(3, confidence) { error = ierror; } //--------------------------------------------------------------------------------------- // 3点求圆 //--------------------------------------------------------------------------------------- void CircleRansac::fitModel(const vector<Point2f>& selects, Mat& result) { result.create(3, 1, CV_32FC1); float a = selects[0].x - selects[1].x; float b = selects[0].y - selects[1].y; float c = selects[0].x - selects[2].x; float d = selects[0].y - selects[2].y; float e = ((selects[0].x * selects[0].x - selects[1].x * selects[1].x) - (selects[1].y * selects[1].y - selects[0].y * selects[0].y)) * 0.5f; float f = ((selects[0].x * selects[0].x - selects[2].x * selects[2].x) - (selects[2].y * selects[2].y - selects[0].y * selects[0].y)) * 0.5f; float cx = (e * d - b * f) / (a * d - b * c); float cy = (a * f - e * c) / (a * d - b * c); float r = sqrtf((cx - selects[0].x) * (cx - selects[0].x) + (cy - selects[0].y) * (cy - selects[0].y)); result.at<float>(0) = cx; result.at<float>(1) = cy; result.at<float>(2) = r; } float CircleRansac::evalModel(const Mat& model, const Point2f& pt) { float cx = model.at<float>(0); float cy = model.at<float>(1); float r = model.at<float>(2); float dist = fabs(sqrtf((cx - pt.x) * (cx - pt.x) + (cy - pt.y) * (cy - pt.y)) - r); return dist < error ? 1 : 0; }
使用方法如下：
int main() { vector<Point2f> points = { ... }; CircleRansac ransac(2.0f); int innerCount = ransac.execute(points); if (innerCount < 4) { return -1; } points.erase(points.begin() + innerCount, points.end()); // 现在points里都是内点 return 0; }