RANSAC算法与原理（二）

文章目录

使用RANSAC进行直线拟合总体最小二乘法直线推导根据几何意义推导从解齐次方程的角度推导方法一方法二总结References

之前介绍了RANSAC算法原理的基本原理和算法实现流程，那么我们在具体使用时，应当如何应用呢。由于RANSAC算法的通用性，其适用于各种模型拟合的场景下，如二维（三维）直线拟合、平面拟合、椭圆拟合等等。skimage官方文档中给出了几个使用RANSAC进行模型拟合的例子。

使用RANSAC进行直线拟合

首先先看一个使用RANSAC进行直线拟合的简单例子。参考RANSAC算法详解(附Python拟合直线模型代码) - 知乎 ()，并对迭代次数的更新部分进行修改。

import numpy as npimport matplotlib.pyplot as pltimport randomimport math# 数据量。SIZE = 50# 产生数据。np.linspace 返回一个一维数组，SIZE指定数组长度。# 数组最小值是0，最大值是10。所有元素间隔相等。X = np.linspace(0, 10, SIZE)Y = 3 * X + 10fig = plt.figure()# 画图区域分成1行1列。选择第一块区域。ax1 = fig.add_subplot(1,1, 1)# 标题ax1.set_title("RANSAC")# 让散点图的数据更加随机并且添加一些噪声。random_x = []random_y = []# 添加直线随机噪声for i in range(SIZE):random_x.append(X[i] + random.uniform(-0.5, 0.5)) random_y.append(Y[i] + random.uniform(-0.5, 0.5)) # 添加随机噪声for i in range(SIZE):random_x.append(random.uniform(0,10))random_y.append(random.uniform(10,40))RANDOM_X = np.array(random_x) # 散点图的横轴。RANDOM_Y = np.array(random_y) # 散点图的纵轴。# 画散点图。ax1.scatter(RANDOM_X, RANDOM_Y)# 横轴名称。ax1.set_xlabel("x")# 纵轴名称。ax1.set_ylabel("y")# 使用RANSAC算法估算模型# 迭代最大次数，每次得到更好的估计会优化iters的数值iters = 100000# 数据和模型之间可接受的差值sigma = 0.25# 最好模型的参数估计和内点数目best_a = 0best_b = 0pretotal = 0# 希望的得到正确模型的概率P = 0.99i=0while True:if i<iters:i+=1# 随机在数据中红选出两个点去求解模型sample_index = random.sample(range(SIZE * 2),2)x_1 = RANDOM_X[sample_index[0]]x_2 = RANDOM_X[sample_index[1]]y_1 = RANDOM_Y[sample_index[0]]y_2 = RANDOM_Y[sample_index[1]]# y = ax + b 求解出a，ba = (y_2 - y_1) / (x_2 - x_1)b = y_1 - a * x_1# 算出内点数目total_inlier = 0for index in range(SIZE * 2):y_estimate = a * RANDOM_X[index] + bif abs(y_estimate - RANDOM_Y[index]) < sigma:total_inlier = total_inlier + 1# 判断当前的模型是否比之前估算的模型好if total_inlier > pretotal:iters = math.log(1 - P) / math.log(1 - pow(total_inlier / (SIZE * 2), 2))pretotal = total_inlierbest_a = abest_b = b# 判断是否当前模型已经符合超过一半的点if total_inlier > SIZE:breakelse:break# 用我们得到的最佳估计画图Y = best_a * RANDOM_X + best_b# 直线图ax1.plot(RANDOM_X, Y)text = "best_a = " + str(best_a) + "\nbest_b = " + str(best_b)plt.text(5,10, text,fontdict={'size': 8, 'color': 'r'})plt.show()

上述示例是使用RANSAC进行直线拟合的一个简单实现，我们可以注意到下述这段代码：

for index in range(SIZE * 2):y_estimate = a * RANDOM_X[index] + bif abs(y_estimate - RANDOM_Y[index]) < sigma:total_inlier = total_inlier + 1

在进行内点和外点的判断选择时，是采用循环的方式逐一判断的，实际上，我们可以利用线性代数的知识，利用numpy矩阵运算将代码简化。

上述的代码，在每次进行直线拟合的时候，使用的两个点来进行直线方程的求解。那么当我想使用大于等于2个的数据点进行直线拟合时，最常见的就是使用最小二乘法拟合直线。

在使用最小二乘法进行直线拟合时，需要注意传统的最小二乘法（least squares)和总体最小二乘法(total least squares)的区别，详细说明请参考wiki链接，分别为最小二乘法 - 维基百科，自由的百科全书 ()和[Total least squares - Wikipedia](/wiki/Total_least_squares#:~:text=In applied statistics%2C total least squares is a,dependent and independent variables are taken into account.)。

两种方法的示意图如下图所示：

知乎上也有作者对这两个问题进行了说明和推导，参见线性拟合1-最小二乘 - 知乎 ()和

线性拟合2-正交回归 - 知乎 ()。

下面我将对总体最小二乘法的推导进行详细说明。

总体最小二乘法直线推导

根据几何意义推导

下面的推导主要是根据线性拟合2-正交回归 - 知乎 ()进行整理得到。

首先需要说明的是，下述推导是二维直线的点法式推导（三维直线没有点法式方程）。点法式的方程相比斜截式的方程，能更完整地的表示所有的直线情况。

横坐标和纵坐标的误差定义如下：

ηi^=xi−xi^\hat{\eta _{i}} =x_{i}-\hat{x_{i}} ηi​^​=xi​−xi​^​

ϵi^=yi−yi^\hat{\epsilon _{i}} =y_{i}-\hat{y_{i}} ϵi​^​=yi​−yi​^​

综合考虑横纵坐标的误差，目标函数的形式如下：

J2=∑in[(ηi^)2+(ϵi^)2]=∑in[(xi−xi^)2+(yi−yi^)2]J_{2}=\sum_{i}^{n} [(\hat{\eta _{i}})^{2}+(\hat{\epsilon _{i}})^{2}]=\sum_{i}^{n} [(x_{i}-\hat{x_{i}})^{2}+(y_{i}-\hat{y_{i}})^{2}] J2​=i∑n​[(ηi​^​)2+(ϵi​^​)2]=i∑n​[(xi​−xi​^​)2+(yi​−yi​^​)2]

因为要求目标函数的最小值，所以就是求观测点(xi,yi)(x_{i},y_{i})(xi​,yi​)到估计点(xi^,yi^)(\hat{x_{i}},\hat{y_{i}})(xi​^​,yi​^​)的最小值，因为点(xi,yi)(x_{i},y_{i})(xi​,yi​)我们是已知的，也就是求点(xi,yi)(x_{i},y_{i})(xi​,yi​)到拟合直线的距离的和的最小值。估计点(xi^,yi^)(\hat{x_{i}},\hat{y_{i}})(xi​^​,yi​^​)是点(xi,yi)(x_{i},y_{i})(xi​,yi​)在拟合直线上的正交投影点。所以目标函数可以写成：

J2=∑indi2J_{2}=\sum_{i}^{n}d_{i}^{2} J2​=i∑n​di2​

did_{i}di​为第iii个点(xi,yi)(x_{i},y_{i})(xi​,yi​)到拟合直线的距离。

用点法式拟合直线的表示形式为

a(x−x0)+b(y−y0)=0a(x-x_0)+b(y-y_0)=0 a(x−x0​)+b(y−y0​)=0

其中，(x0,y0)(x_0,y_0)(x0​,y0​)为经过直线的一个点，(a,b)(a,b)(a,b)表示直线的法向量。(a,b)(a,b)(a,b)仅仅表示一个方向，为了计算方便，我们取单位向量，即满足下式：

a2+b2=1a^2+b^2=1 a2+b2=1

第iii个点到直线的距离，可以表示为向量(xi−x0,yi−y0)(x_i-x_0,y_i-y_0)(xi​−x0​,yi​−y0​)在法向量方向(a,b)(a,b)(a,b)投影的长度，目标函数可以写成：

J2=∑indi2=∑in((xi−x0,yi−y0)⋅(a,b))2a2+b2=∑in[a(xi−x0)+b(yi−y0)]2J_2=\sum_{i}^{n}d_i^{2}=\sum_{i}^{n}\frac{((x_i-x_0,y_i-y_0)\cdot(a,b))^2}{a^2+b^2}=\sum_{i}^{n}[a(x_i-x_0)+b(y_i-y_0)]^2 J2​=i∑n​di2​=i∑n​a2+b2((xi​−x0​,yi​−y0​)⋅(a,b))2​=i∑n​[a(xi​−x0​)+b(yi​−y0​)]2

由极值点和可导函数之间的关系得，极值点可以推出一阶导数为0，那么一阶导数为0就是极值点的必要条件。那么可以推出：一阶导数不为0，不是极值点。

将目标函数J2J_2J2​对x0x_0x0​和y0y_0y0​求偏导，并令其等于0，得：

∂J2∂x0=−2a∑in[a(xi−x0)+b(yi−y0)]=0\frac{\partial J_2}{\partial x_0}=-2a\sum_{i}^{n}[a(x_i-x_0)+b(y_i-y_0)]=0 ∂x0​∂J2​​=−2ai∑n​[a(xi​−x0​)+b(yi​−y0​)]=0

∂J2∂y0=−2b∑in[a(xi−x0)+b(yi−y0)]=0\frac{\partial J_2}{\partial y_0}=-2b\sum_{i}^{n}[a(x_i-x_0)+b(y_i-y_0)]=0 ∂y0​∂J2​​=−2bi∑n​[a(xi​−x0​)+b(yi​−y0​)]=0

上式等号两边同时除以n得

a(xˉ−x0)+b(yˉ−y0)=0a(\bar{x}-x_0)+b(\bar{y}-y_0)=0 a(xˉ−x0​)+b(yˉ​−y0​)=0

可以看到点(xˉ,yˉ)(\bar{x},\bar{y})(xˉ,yˉ​)满足直线方程，所以令x0=xˉ,y0=yˉx_0=\bar{x},y_0=\bar{y}x0​=xˉ,y0​=yˉ​,此时目标函数变为：

J2=∑in[a(xi−xˉ)+b(yi−yˉ)]2=[ab][∑in(xi−xˉ)2∑in(xi−xˉ)(yi−yˉ)∑in(xi−xˉ)(yi−yˉ)∑in(yi−yˉ)2][ab]J_2=\sum_{i}^{n}[a(x_i-\bar{x})+b(y_i-\bar{y})]^2=\\ \begin{bmatrix} a & b \end{bmatrix} \begin{bmatrix} \sum_{i}^{n}(x_i-\bar{x})^2 &\sum_{i}^{n}(x_i-\bar{x})(y_i-\bar{y}) \\ \sum_{i}^{n}(x_i-\bar{x})(y_i-\bar{y}) &\sum_{i}^{n}(y_i-\bar{y})2 \end{bmatrix} \begin{bmatrix} a\\ b \end{bmatrix} J2​=i∑n​[a(xi​−xˉ)+b(yi​−yˉ​)]2=[a​b​][∑in​(xi​−xˉ)2∑in​(xi​−xˉ)(yi​−yˉ​)​∑in​(xi​−xˉ)(yi​−yˉ​)∑in​(yi​−yˉ​)2​][ab​]

对目标函数J2J_2J2​除以n可得：

J2n=[ab][SxxSxySxySyy][ab]=vTSv\frac{J_2}{n}= \begin{bmatrix} a & b \end{bmatrix} \begin{bmatrix} S_{xx} & S_{xy}\\ S_{xy} & S_{yy} \end{bmatrix} \begin{bmatrix} a\\ b \end{bmatrix}=v^TSv nJ2​​=[a​b​][Sxx​Sxy​​Sxy​Syy​​][ab​]=vTSv

其中SxxS_{xx}Sxx​和SyyS_{yy}Syy​为xxx和yyy的方差,SxyS_{xy}Sxy​为xxx和yyy的协方差。

很明显，这是一个二次型求最小值的问题。相关定理证明也可参考线性代数及其应用 (豆瓣) ()

因为SSS是实对称矩阵，所以可以将其进行正交对角化分解：

S=[q1q2][λ100λ2][q1Tq2T]=QΛQTS= \begin{bmatrix} q_1 & q_2 \end{bmatrix} \begin{bmatrix} \lambda_1 &0 \\ 0 &\lambda_2 \end{bmatrix} \begin{bmatrix} q_{1}^{T}\\ q_{2}^{T} \end{bmatrix}=Q\Lambda Q^T S=[q1​​q2​​][λ1​0​0λ2​​][q1T​q2T​​]=QΛQT

有

J2n=vTSv=(vTQ)Λ(vTQ)T\frac{J_2}{n}=v^TSv=(v^TQ)\Lambda(v^TQ)^T nJ2​​=vTSv=(vTQ)Λ(vTQ)T

令u1=vTq1,u2=vTq2,u=[u1u2]u_1=v^Tq_1,u_2=v^Tq_2,u=\begin{bmatrix} u_1\\ u_2 \end{bmatrix}u1​=vTq1​,u2​=vTq2​,u=[u1​u2​​]则：

J2n=uTΛu=λ1u12+λ2u22\frac{J_2}{n}=u^T\Lambda u=\lambda_1u_1^2+\lambda_2u_2^2 nJ2​​=uTΛu=λ1​u12​+λ2​u22​

因为：

uTu=vTQQTv=vTv=1u^Tu=v^TQQ^Tv=v^Tv=1 uTu=vTQQTv=vTv=1

所以uuu为单位矩阵，即u12+u22=1u_1^2+u_2^2=1u12​+u22​=1

不妨设λ1≤λ2\lambda_1\le\lambda_2λ1​≤λ2​,即当u1=1,u2=0u_1=1,u_2=0u1​=1,u2​=0时，J2J_2J2​取得最小值λ1\lambda_1λ1​，即v=q1v=q_1v=q1​

所以最终的直线拟合结果是法向量vvv对应矩阵SSS的最小特征值的特征向量。

结果整理：

x0=xˉ,y0=yˉx_0=\bar{x},y_0=\bar{y} x0​=xˉ,y0​=yˉ​

拟合直线法向量：

v=[ab]=q1v=\begin{bmatrix} a\\ b \end{bmatrix}=q_1 v=[ab​]=q1​

那么如何求矩阵SSS的最小特征值对应的特征向量呢

下面介绍使用奇异值分解求解矩阵SSS的最小特征值对应的特征向量，参考奇异值分解(SVD)原理与在降维中的应用 - 刘建平Pinard - 博客园 ()

假设现在有nnn个二维点(x1,y1),(x2,y2)...(xn,yn)(x_1,y_1),(x_2,y_2)...(x_n,y_n)(x1​,y1​),(x2​,y2​)...(xn​,yn​),那么我们可以直接得到

A=[x1−xˉy1−yˉx2−xˉy2−yˉ......xn−xˉyn−yˉ]A=\begin{bmatrix} x_1-\bar{x}&y_1-\bar{y} \\ x_2-\bar{x}&y_2-\bar{y} \\ ...&... \\ x_n-\bar{x}&y_n-\bar{y} \end{bmatrix} A=⎣⎢⎢⎡​x1​−xˉx2​−xˉ...xn​−xˉ​y1​−yˉ​y2​−yˉ​...yn​−yˉ​​⎦⎥⎥⎤​

那么我们可以观察得到，矩阵ATAA^TAATA就是矩阵SSS，那么根据奇异值分解的定义

A=UΣVTA=U\Sigma V^T A=UΣVT

上式中的VVV就是矩阵ATAA^TAATA也就是SSS的特征值对应的特征向量构成的正交矩阵，又因为在奇异值分解时，矩阵VVV中的特征值已经根据大小从大到小重新排列了，所以我们需要的就是VVV的最后一列，也就是VTV^TVT的最后一行。

在python中，我们可以很方便的使用numpy中的np.linalg.svd求得相应的特征矩阵和特征向量。

从解齐次方程的角度推导

参考SVD之最小二乘【推导与证明】 - Brad_Lucas - 博客园 ()和奇异值分解（SVD）方法求解最小二乘问题的原理 - 一抹烟霞 - 博客园 ()

求解直线拟合的问题，可以转化为求解齐次方程Ax=0Ax=0Ax=0的问题，其中

A=[x1−xˉy1−yˉx2−xˉy2−yˉ......xn−xˉyn−yˉ]x=[ab]A=\begin{bmatrix} x_1-\bar{x}&y_1-\bar{y} \\ x_2-\bar{x}&y_2-\bar{y} \\ ...&... \\ x_n-\bar{x}&y_n-\bar{y} \end{bmatrix}\\x= \begin{bmatrix} a\\b \end{bmatrix} A=⎣⎢⎢⎡​x1​−xˉx2​−xˉ...xn​−xˉ​y1​−yˉ​y2​−yˉ​...yn​−yˉ​​⎦⎥⎥⎤​x=[ab​]

可以想到，当直线拟合的越好时，∣∣Ax∣∣||Ax||∣∣Ax∣∣越小

方法一

求min∣∣Ax∣∣||Ax||∣∣Ax∣∣，根据正交矩阵的性质：

∣∣Ax∣∣=∣∣UΣVTx∣∣=∣∣ΣVTx∣∣||Ax||=||U\Sigma V^Tx||=||\Sigma V^Tx|| ∣∣Ax∣∣=∣∣UΣVTx∣∣=∣∣ΣVTx∣∣

令y=VTxy=V^Txy=VTx，则：

∣∣ΣVTx∣∣=∣∣Σy∣∣||\Sigma V^Tx||=||\Sigma y|| ∣∣ΣVTx∣∣=∣∣Σy∣∣

进一步处理，求min(yTDTDy)min(y^TD^TDy)min(yTDTDy)

yTΣTΣy=σ12y12+⋯+σn2yn2σ1≥σ2≥…σn≥0y^T\Sigma^T\Sigma y=\sigma_1^2y_1^2+\dots+\sigma_n^2y_n^2\\ \sigma_1\ge\sigma_2\ge\dots\sigma_n\ge0 yTΣTΣy=σ12​y12​+⋯+σn2​yn2​σ1​≥σ2​≥…σn​≥0

约束条件为∣∣y∣∣=1||y||=1∣∣y∣∣=1,如上所示，y=[00…1]y=\begin{bmatrix} 0&0&\dots&1\end{bmatrix}y=[0​0​…​1​]是最小解，即：

x=Vy=[V1V2…Vn][000…1]=Vnx=Vy=\begin{bmatrix}V_1&V_2&\dots&V_n\end{bmatrix}\begin{bmatrix}0\\0\\0\\ \dots\\1\end{bmatrix}=V_n x=Vy=[V1​​V2​​…​Vn​​]⎣⎢⎢⎢⎢⎡​000…1​⎦⎥⎥⎥⎥⎤​=Vn​

即最优解是VVV的最后一列

方法二

求min∣∣Ax∣∣=min(xTATAx)min||Ax||=min(x^TA^TAx)min∣∣Ax∣∣=min(xTATAx),其中A=UΣVT,UTU=I,VTV=IA=U\Sigma V^T,U^TU=I,V^TV=IA=UΣVT,UTU=I,VTV=I,,,Σ\SigmaΣ为对角阵。

ATA=VΣTUTUΣVT=VΣTΣVT=V[σmax2⋱σmin2]VTA^TA=V\Sigma^TU^TU\Sigma V^T=V\Sigma^T\Sigma V^T\\=V\begin{bmatrix}\sigma_{max}^2&&&\\ &\ddots&&\\&&&\sigma_{min}^2 \end{bmatrix}V^T\\ ATA=VΣTUTUΣVT=VΣTΣVT=V⎣⎡​σmax2​​⋱​​σmin2​​⎦⎤​VT

令V=[V0V1…Vn]V=\begin{bmatrix}V_0&V_1&\dots&V_n\end{bmatrix}V=[V0​​V1​​…​Vn​​],得：

ATA=σmax2V0V0T+⋯+σmin2VnVnTA^TA=\sigma_{max}^2V_0V_0^T+\dots+\sigma_{min}^2V_nV_n^T ATA=σmax2​V0​V0T​+⋯+σmin2​Vn​VnT​

由正交矩阵的性质：

VTV=IViTVj=0.i≠j∣∣Vi∣∣=1V^TV=I\\ V_i^TV_j=0.i\ne j\\ ||V_i||=1 VTV=IViT​Vj​=0.i​=j∣∣Vi​∣∣=1

令x=∑0nkiVix=\sum_0^nk_iV_ix=∑0n​ki​Vi​，其中ViV_iVi​为基向量

xTATAx=σmax2k02V0TV0V0TV0+⋯+σmin2kn2VnTVnVnTVn=k02σmax2+⋯+kn2σmin2x^TA^TAx=\sigma_{max}^2k_0^2V_0^TV_0V_0^TV_0+\dots+\sigma_{min}^2k_n^2V_n^TV_nV_n^TV_n\\=k_0^2\sigma_{max}^2+\dots+k_n^2\sigma_{min}^2 xTATAx=σmax2​k02​V0T​V0​V0T​V0​+⋯+σmin2​kn2​VnT​Vn​VnT​Vn​=k02​σmax2​+⋯+kn2​σmin2​

由于∣∣x∣∣=1||x||=1∣∣x∣∣=1，则

∑0nki2=1\sum_0^nk_i^2=1 0∑n​ki2​=1

上式取最小的情况为kn=1,ki=0(i≠n)k_n=1,k_i=0(i\ne n)kn​=1,ki​=0(i​=n)

此时x=Vnx=V_nx=Vn​，为σmin\sigma_{min}σmin​对应的列向量。

总结

我们发现，如果把上述的推导过程变为三维，把二维直线的点法式方程变为空间平面的点法式方程，上述推导依旧是成立的，所以我们已经得到了拟合空间平面的理论基础了。

References

RANSAC算法详解(附Python拟合直线模型代码) - 知乎 ()

线性拟合1-最小二乘 - 知乎 ()

线性拟合2-正交回归 - 知乎 ()

线性代数及其应用 (豆瓣) ()

奇异值分解(SVD)原理与在降维中的应用 - 刘建平Pinard - 博客园 ()

奇异值分解（SVD）方法求解最小二乘问题的原理 - 一抹烟霞 - 博客园 ()

SVD之最小二乘【推导与证明】 - Brad_Lucas - 博客园 ()