RANSAC算法——过程及代码

一、RANSAC算法概述

RANSAC(RANdom SAmple Consensus)算法是一种可以去除噪声影响，估计模型的算法。它比类似的最小二乘法更robust。
RANSAC使用数据的子集(inlier)来估计模型，然后通过现有模型参数，尽可能扩大模型在样本集合中的影响，拟合较多的点。

二、算法过程

1. RANSAC算法的过程如下：

对于$N$个点构成的集合$P$，我们假定最少通过$n$个点可以拟合出正确的模型，也就是说样本集合中有$N-n$个噪声点。具体的执行下列操作：

1. 从数据集中随机选择一组样本构建成内点集合，拟合基础模型$H$
2. 使用剩余数据对$H$进行测试，将落在预定误差($\sigma$)范围内的点划到内点集合中。
3. 使用全部的内点集合拟合模型$H$
4. 使用内点集合来估计模型的误差
5. 如果模型的性能达到了设定的阈值，或者达到了预定的迭代次数，则算法终止，否则跳转到第2步继续执行。

2. 参数说明

a. 初始内点的确定

算法第1步中选择内点时，通常选择能确定模型的最少点，比如两点确定一条直线，3点确定一个平面。在求解单应性矩阵的时候，就是4个点对8个点确定一个矩阵模型。

b. 迭代次数的确定

通常来说$N$的值比较大，那么随机选择点的时候，点与点的组合就很多，如果不对迭代次数进行限制，运算量就会很大。因为算法第2步要用剩余点对模型进行测试，如果点的数量较多，这一步的运算量很大。包括后面的第3，4步。其实我们只要保证我们估计模型的时候使用的都是内点即可。
假设是内点的概率为$t$：
$$t=\frac{n_i}{n_i+n_o}$$
$n_i$为内点数量，$n_o$为外点数量。那么我们每次计算模型使用$n$个点的情况下，选取的点集中至少有一个外点的概率就是：$1-t^n$。那么在迭代$k$次的情况下，$(1-t^n{)}^k$就是$k$次迭代都至少有个为外点的概率。
那么能够得到$n$个正确的点来估计模型的概率就是：
$$P=1-(1-t^n{)}^k$$
两边取对数，就可以得到最少迭代次数：
$$k=\frac{log(1-P)}{log(1-t^n)}$$
这里$P$是希望通过算法得到正确模型的最少概率，$k$是关于$P$单调递增的，$k$就是在$P$确定情况下的最少迭代次数。$t$通常未知，那么我们可以采用自适应的办法，在开始的时候设置一个较大的迭代次数$k$，然后通过每次迭代中更新$t$来更新迭代次数。

三、代码示例

这里使用RANSAC算法来拟合直线，代码如下：

import numpy as np
import matplotlib.pyplot as plt
import random
import math

class ransacMatchingTest(object):

    def __init__(self, X, Y, sigma=0.2, prob=0.95):

        self.X = X
        self.Y = Y
        self.sigma = sigma
        self.prob = prob

    def runMatch(self):

        preInlier = 0
        iters = 10000
        best_a = 0
        best_b = 0
        for i in range(iters):
            # sample points
            sampleIndex = random.sample(range(self.X.shape[0]), 2)
            x1 = self.X[sampleIndex[0]]
            y1 = self.Y[sampleIndex[0]]

            x2 = self.X[sampleIndex[1]]
            y2 = self.Y[sampleIndex[1]]

            # Compute mode
            a = (y2 - y1) / (x2 - x1)
            b = y1 - a * x1

            # Get number of inlier
            inlier = 0
            for index in range(self.X.shape[0]):
                y_estimate = a * self.X[index] + b
                if abs(y_estimate - self.Y[index]) < self.sigma:
                    inlier = inlier + 1

            if inlier > preInlier:
                # Update iteration times
                iters = math.log(1-self.prob) / math.log(1 - pow(inlier / self.X.shape[0], 2))
                preInlier = inlier
                best_a = a
                best_b = b
            # While inliers over half of input set then break
            if inlier > (self.X.shape[0] / 2):
                break

        return best_a, best_b


if __name__ == '__main__':
    X = np.linspace(0, 10, 50)
    Y = 4 * X + 6

    randomX = []
    randomY = []
    # Add mode noise
    for i in range(50):
        randomX.append(X[i] + random.uniform(-0.5, 0.5))
        randomY.append(Y[i] + random.uniform(-0.5, 0.5))
    # Add random noise
    for _ in range(50):
        randomX.append(random.uniform(0, 6))
        randomY.append(random.uniform(6, 30))

    RANDOMX = np.array(randomX)
    RANDOMY = np.array(randomY)

    ransac = ransacMatchingTest(RANDOMX, RANDOMY, sigma=0.2, prob=0.95).runMatch()

    preY = ransac[0] * RANDOMX + ransac[1]

    fig = plt.figure()
    ax1 = fig.add_subplot(1, 1, 1)
    ax1.scatter(RANDOMX, RANDOMY)

    ax1.plot(RANDOMX, preY)
    plt.show()