圆拟合算法

圆拟合算法

最小二乘法拟合(least square fitting)

使用最小二乘法，对几何特征进行拟合，主要有两大类方式：代数拟合和几何拟合，这两种方式的区别，在于它们对误差距离的定义方式。

代数距离：为要进行拟合的几何特征(直线、圆、椭圆等)写下隐式方程，然后将点的坐标带入此隐式方程就得到了点到该几何特征的距离。这种方法对于直线拟合、圆拟合，返回的就是到其的真实距离。而对于椭圆拟合，它返回的值是与距离有类似的属性，但不是一个真正的距离值。因此这个距离被称为 代数距离。

代数拟合 (algebraic fitting)

代数拟合的几何特征是通过带有参数向量$\vec{a}=(a_1,\cdots ,a_q)$ 的隐式方程 $F(X,a)=0$ 来表示的。其误差距离是由隐式方程与每个给定点的期望值(即零)的偏差定义。但方程的不确定性导致给出的点会出现不位于图形上的现象(即存在一些拟合误差)。

关于圆的最小二乘法拟合，大多数论文给出的其代数距离平方和及其修正公式，如下所示：
$$\sum _{i=1}^m{e}_i^2=\sum _{i=1}^m[(x_i-x_c{)}^2+(y_i-y_c{)}^2-r^2{]}^2\text{\hspace{1em}(1)}$$
尽管代数拟合在实现和计算成本方面具有优势，但在精度方面以及与拟合参数和误差的物理解释相关方面都存在缺陷。其中，已知的代数拟合的缺点，有以下几点：

• 误差距离的定义与测量准则不符；
• 拟合参数的可靠性，难以检验；
• 拟合参数对坐标变换来说是可变的(如，椭圆拟合过程中，对给定点进行一个简单的平移变换，不仅会造成中心点的坐标发生变化，而且还会造成椭圆轴的长度及角度的变化)；
• 预估的拟合参数存在偏差；
• 拟合误差存在不合理的加权；
• 拟合过程有时会以意外的几何特征结束(如，可能会以双曲线代替椭圆的拟合)；

Kåsa Fit

注：网络上查到的大多数最小二乘法求拟合圆的算法，均是该算法的实现。

以下是 Kåsa Fit 的求解过程：

首先，根据圆的公式$(x-A{)}^2+(y-B{)}^2=R^2$ 可得，
$$R^2=x^2-2Ax+A^2+y^2-2By+B^2;$$
令：
$a=-2A,b=-2B,c=A^2+B^2-R^2$
可得圆的方程的另一种形式：
$$x^2+y^2+ax+by+c=0\text{\hspace{1em}(1)}$$

根据公式$(1)$，只要求出参数 a, b, c 即可得到圆心及半径的参数；
$$A=-\frac{a}{2}，B=-\frac{b}{2}，R=\frac{\sqrt{a^2+b^2-4c}}{2}$$

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-HoQMhYw6-1584789138214)(图1.png)]

样本集$(X_i,Y_i),i\in (1,2,3,\cdots ,n$ 中，点到圆心的距离为 $d_i$：
$$d_i^2=(X_i-A{)}^2+(Y_i-B{)}^2$$
令，点$(X_i,Y_i)$到圆边缘的距离的平方与半径的平方差为:
$$\delta =d_i^2-R^2=(X_i-A{)}^2+(Y_i-B{)}^2-R^2=X_i^2+Y_i^2+a{X}_i+b{Y}_i+c$$

令，$Q(a,b,c)$ 为 ${\delta }^2$ 的平方和:
$$Q(a,b,c)=\sum _{i=1}^n{\delta }_i^2=\sum _{i=1}^n[(X_i^2+Y_i^2+a{X}_i+b{Y}_i+c){]}^2$$
求参数$(a,b,c)$ 使得 $Q(a,b,c)$ 的值最小。

解：平方差 $Q(a,b,c)$ 大于等于 $0$，因此，函数存在大于或等于 $0$ 的极小值，其最大值为无穷大；
$F(a,b,c)$ 对 $(a,b,c)$，令偏导等于 $0$，则得到极值点，比较所有的极值点函数，即可得到极小值；
$$\frac{\partial Q(a,b,c)}{\partial a}=\partial (\sum _{i=1}^n[(X_i^2+Y_i^2+a{X}_i+b{Y}_i+c){]}^2)$$
令，常数 $D=(X_i^2+Y_i^2+b{Y}_i+c)$，则上式可转化为：
$$\frac{\partial Q(a,b,c)}{\partial a}=\partial (\sum _{i=1}^n[(D+a{X}_i){]}^2)$$
则，
$$\frac{\partial Q(a,b,c)}{\partial a}=\partial (\sum _{i=1}^n[(D^2+2Da{X}_i+X_i^2{a}^2)])=\sum _{i=1}^n(2X_i^{2}a+2D{X}_i)=\sum _{i=1}^{n}2(a{X}_i+D)X_i$$
将 $D$代入上式，得:
$$\frac{\partial Q(a,b,c)}{\partial a}=\sum _{i=1}^{n}2(X_i^2+Y_i^2+a{X}_i+b{Y}_i+c)X_i=0\text{\hspace{1em}(2)}$$
同理，可得 $b,c$ 的偏导数：
$$\frac{\partial Q(a,b,c)}{\partial b}=\sum _{i=1}^{n}2(X_i^2+Y_i^2+a{X}_i+b{Y}_i+c)Y_i=0\text{\hspace{1em}(3)}$$

$$\frac{\partial Q(a,b,c)}{\partial c}=\sum _{i=1}^{n}2(X_i^2+Y_i^2+a{X}_i+b{Y}_i+c)=0\text{\hspace{1em}(4)}$$

将上述等式，改写成方程的形式：
$$\{\begin{array}{l}\frac{\partial Q(a,b,c)}{\partial a}={\sum }_{i=1}^{n}2(X_i^2+Y_i^2+a{X}_i+b{Y}_i+c)X_i=0\\ \frac{\partial Q(a,b,c)}{\partial b}={\sum }_{i=1}^{n}2(X_i^2+Y_i^2+a{X}_i+b{Y}_i+c)Y_i=0\\ \frac{\partial Q(a,b,c)}{\partial c}={\sum }_{i=1}^{n}2(X_i^2+Y_i^2+a{X}_i+b{Y}_i+c)=0\end{array}$$

根据 $\sum$ 运算法则，将常数提出，并且因为等式等于$0$，所以，可以直接省略，得下式；

$$\{\begin{array}{l}\frac{\partial Q(a,b,c)}{\partial a}={\sum }_{i=1}^n(X_i^3+X_i{Y}_i^2+a{X}_i^2+b{X}_i{Y}_i+c{X}_i)=0\\ \frac{\partial Q(a,b,c)}{\partial b}={\sum }_{i=1}^n(X_i^2{Y}_i+Y_i^3+a{X}_i{Y}_i+b{Y}_i^2+c{Y}_i)=0\\ \frac{\partial Q(a,b,c)}{\partial c}={\sum }_{i=1}^n(X_i^2+Y_i^2+a{X}_i+b{Y}_i+c)=0\end{array}$$

将上式展开，成下列形式：

$$\{\begin{array}{l}{\sum }_{i=1}^n{X}_i^3+{\sum }_{i=1}^n{X}_i{Y}_i^2+a{\sum }_{i=1}^n{X}_i^2+b{\sum }_{i=1}^n{X}_i{Y}_i+c{\sum }_{i=1}^n{X}_i=0\\ {\sum }_{i=1}^n{X}_i^2{Y}_i+{\sum }_{i=1}^n{Y}_i^3+a{\sum }_{i=1}^n{X}_i{Y}_i+b{\sum }_{i=1}^n{Y}_i^2+c{\sum }_{i=1}^n{Y}_i=0\\ {\sum }_{i=1}^n{X}_i^2+{\sum }_{i=1}^n{Y}_i^2+a{\sum }_{i=1}^n{X}_i+b{\sum }_{i=1}^n{Y}_i+{\sum }_{i=1}^{n}c=0\end{array}$$

根据求和公式，得${\sum }_{i=1}^{n}c=Nc$，代入方程得：
$$\{\begin{array}{l}{\sum }_{i=1}^n{X}_i^3+{\sum }_{i=1}^n{X}_i{Y}_i^2+a{\sum }_{i=1}^n{X}_i^2+b{\sum }_{i=1}^n{X}_i{Y}_i+c{\sum }_{i=1}^n{X}_i=0\\ {\sum }_{i=1}^n{X}_i^2{Y}_i+{\sum }_{i=1}^n{Y}_i^3+a{\sum }_{i=1}^n{X}_i{Y}_i+b{\sum }_{i=1}^n{Y}_i^2+c{\sum }_{i=1}^n{Y}_i=0\\ {\sum }_{i=1}^n{X}_i^2+{\sum }_{i=1}^n{Y}_i^2+a{\sum }_{i=1}^n{X}_i+b{\sum }_{i=1}^n{Y}_i+Nc=0\end{array}$$

所以，解方程组，消去 $c$，式$(2)\ast N-(4)\ast {\sum }_{i=1}^n{X}_i$ 得
$$(N\sum _{i=1}^n{X_i}^2-\sum _{i=1}^n{X}_i\sum _{i=1}^n{X}_i)a+N(\sum _{i=1}^n{X}_i{Y}_i-\sum _{i=1}^n{X}_i\sum _{i=1}^n{Y}_i)b+$$
$$N\sum _{i=1}^n{X}^3+N\sum _{i=1}^n{X}_i{Y}_i^2-\sum _{i=1}^n(X_i^2+Y_i^2)\sum _{i=1}^n{X}_i=0\text{\hspace{1em}(5)}$$

$(3)\ast N-(4)\ast {\sum }_{i=1}^n{Y}_i$，得
$$N(\sum _{i=1}^n{X}_i{Y}_i-\sum _{i=1}^n{X}_i\sum _{i=1}^n{Y}_i)a+(N\sum _{i=1}^n{Y_i}^2-\sum _{i=1}^n{Y}_i\sum _{i=1}^n{Y}_i)b+$$
$$N\sum _{i=1}^n{X}_i^2{Y}_i+N\sum _{i=1}{Y}_i^3-\sum _{i=1}^n(X_i^2+Y_i^2)\sum _{i=1}^n{Y}_i=0\text{\hspace{1em}(6)}$$

令：
$$C=(N\sum _{i=1}^n{X_i}^2-\sum _{i=1}^n{X}_i\sum _{i=1}^n{X}_i),$$
$$D=N(\sum _{i=1}^n{X}_i{Y}_i-\sum _{i=1}^n{X}_i\sum _{i=1}^n{Y}_i),$$
$$E=N\sum _{i=1}^n{X}_i^3+N\sum _{i=1}^n{X}_i{Y}_i^2-\sum _{i=1}^n(X_i^2+Y_i^2)\sum _{i=1}^n{X}_i$$
$$G=(N\sum _{i=1}^n{Y_i}^2-\sum _{i=1}^n{Y}_i\sum _{i=1}^n{Y}_i)$$
$$H=N\sum _{i=1}^n{X}_i^2{Y}_i+N\sum _{i=1}{Y}_i^3-\sum _{i=1}^n(X_i^2+Y_i^2)\sum _{i=1}^n{Y}_i$$
则，得方程:
$$\{\begin{array}{l}Ca+Db+E=0\\ Da+Gb+H=0\end{array}$$

解得：
$$a=\frac{HD-EG}{CG-D^2}$$
$$b=\frac{HC-ED}{D^2-CG}$$
$$c=-\frac{{\sum }_{i=1}^n(X_i^2+Y_i^2)+a{\sum }_{i=1}^n{X}_i+b{\sum }_{i=1}^n{Y}_i}{N}$$

得到 $A、B、R$ 的估计值，
$$A=-\frac{a}{2},B=-\frac{b}{2}，R=\frac{\sqrt{a^2+b^2-4c}}{2}$$

注：该方法常用于求迭代方式拟合圆的初始值；

使用 Kåsa 的好处之一就是，如果点$(x_i,y_i)$ 碰巧在圆上，则该方法会找到正确的圆；但当样本点位于小圆弧上时， Kåsa 会偏离很大地估计出一个更小的圆。

以下两种，仅作了解，不再给出详细推导过程

Pratt Fit

kasa方法圆拟合作为最常见的圆拟合方法，虽然计算方法简单，效率快，但是拟合结果存在较大偏差(Heavy bias)。通过将 $D=\sqrt{(x-a{)}^2+(y-b{)}^2}$ 与半径R的差值 $D-R$ 转换为 $D^2-R^2$，将非线性问题转换为线性问题。但是因为 $D^2-R^2=d(d+2R)$(令 $d=D-R$)，当偏离值 $d$ 较大时，Kåsa 方法导致R明显变小。Pratt 通过将较大时，Kåsa 方法的目标除以 $(2R{)}^2$ 的方式，取得更准确的结果。
$$f=\frac{{\sum }_{i=1}^n((x_i-a{)}^2+(y_i-b{)}^2-R^2{)}^2}{(2R{)}^2}$$

类似于Kåsa method，pratt 方法将圆方程描述为 $A(x^2+y^2)+Bx+Cy+D=0$。这样圆心 $(a,b)=(-\frac{B}{2A},-\frac{C}{2A}$，半径 $R=\frac{\sqrt{B^2+X^2-4AD}}{2A}$。注意，A有可能为0，此时圆拟合方法得到的是一条直线，或者说一个曲率为0的圆。但是在这里，因为A,B,C,D的值乘以一个标量也不会改变圆的方程，可以让A=1,所以目标方程可以改写为：
$$f=(\frac{\sum A(x^2+Y^2)+Bx+Cy+D}{B^2+C^2-4AD}{)}^2$$

Pratt 将这个方程转化为求:$g=(\sum A(x_i^2+y_i^2)+B{x}_i+C{y}_i+D{)}^2$ 的最小值，并在约束条件:$B^2+C^2-4AD=1$ 的条件下。
惯例上将上式写成矩阵的形式：
$$g=A^{T}MA(\vec{A}=[A,B,C,D];M=\left[\begin{array}{cccc}\sum z^2& \sum xz& \sum yz& \sum z\\ \sum xz& \sum x^2& \sum xy& \sum x\\ \sum yz& \sum xy& \sum y^2& \sum y\\ \sum z& \sum x& \sum y& n\end{array}\right];z_i=x_i^2+y_i^2)$$
$$A^{T}BA=\left[\begin{array}{cccc}A& B& C& D\end{array}\right]\left[\begin{array}{cccc}0& 0& 0& -2\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ -2& 0& 0& 0\end{array}\right]\left[\begin{array}{c}A\\ B\\ C\\ D\end{array}\right]=1$$

通过引入拉格朗日乘子解决最小值问题：
$$g(A,\eta )=ATMA–\eta (ATBA-1)$$

对 $\vec{A}$ 求偏导: $MA–\eta BA=0$。因为$B$ 可逆，可以变换为 $B^{-1}MA=\eta A$。可知 $\eta$为矩阵 $B^{-1}M$ 的特征值，A为其特征向量。

由 $MA–\eta BA=0$可得$ATMA–\eta ATBA=ATMA-\eta =0$,即$ATMA$的极值为$\eta$，那么矩阵$B^{-1}M$有四个特征值，最小的特征值就是$ATMA$的最小值？不然。如果 $\eta <0$ ,那么拟合出的圆将没有意义，因为$ATMA\ge 0$。（Nikolai Chernov 教授经过复杂的证明，$B^{-1}M$ 存在一个负特征值。）因此，$A$的值应该是对应矩阵$B^{-1}M$ 最小非负特征值的特征向量。

算法分析

pratt 方法有两个明显优势，一个是消除了半径的影响，使得曲率半径非常大、噪音大的情况下算法也能正确运行。

一个是在方程式上融合了直线方程，当系数A=0时，求得的曲线是直线，在实际应用中，不清楚待拟合曲线为直线还是曲线的情况下非常有效。

还有一点就是 pratt 方法的效率也比较高，因此是一种非常实用的方法。

Taubin Fit

注:暂时只作简单了解

Taubin Fit 的最小拟合函数为：
$$F_T=\frac{[\sum (x_i-a{)}^2+(y_i-b{)}^2-R^2{]}^2}{4n^{-1}\sum [(x_i-a{)}^2+(y_i-b{)}^2]}$$
转化为下列形式：
$$F_T=\frac{[A{z}_i^2+B{x}_i+C{y}_i+D{]}^2}{n^{-1}\sum 4A^2{z}_i+4AB{x}_i+4AC{y}_i+B^2+C^2}$$
将上述公式，求最小值，等效为：
$$4A^2\overline{z}+4AB\overline{x}+4AC\overline{y}+B^2+C^2=1$$
其中，$\overline{x}=\frac{1}{n}\sum x_i$。

迭代重加权最小二乘法(IRLS, iteratively reweighted least-squares)

加权最小二乘法是对原模型进行加权，使之成为一个新的不存在异方差性的模型，然后采用普通最小二乘法估计其参数的一种数学优化技术。加权迭代最小二乘 IRLS属于 稳健回归方法，常用于样本较少的回归，IRLS 属于 M-Estimators。

迭代加权最小二乘法在图形拟合中的应用

最小二乘法拟合图形，对于离群值是不足够可靠的，由于采用的是平方距离，所以，与图形距离远的那些点，在最优化过程中会拥有非常大的权重值。为了减轻这些影响，我们可以为每个点引入权重 $w_i$。对于远离线的那些点，其权重应该 $\ll 1$。

权重值 $w_i$ 该如何定义？因为我们想要给出距离远的那些点较小的权重，所以权重必须是基于点到图形的距离。不幸的是，在我们还没有拟合出图形的方程时，我们并不能得到这些距离，所以这看起来是个不可能的要求。解决的方法是用多次的迭代来拟合直线。在第一次迭代中使用 $w_i=1$，即执行一个标准的直线拟合计算出距离${\delta }_i$。通过用一个权重函数 $w(\delta )$ 可用已经计算出的距离来定义权重，这些权重用于后续的迭代处理中，这个方法称为迭代重加权最小二乘法(IRLS)。

在实践中，可以使用以下两个权重函数中的任何一个，这两个权重函数的应用效果都很好。第一个权重函数是由 Huber 提出的，其定义为：
$$w(\delta )=\{\begin{array}{l}1,|\delta |\le \tau \\ \frac{\tau }{|\delta |}|\delta |>\tau \end{array}$$

上式中，参数 $\tau$ 是削波因数，它定义哪些点应该被视为离群值。请记住所有距离 $\le \tau$ 的点对应的权重都是 $1$。这意味着，对于小的距离的点，在极小化处理中就直接使用平方距离。另一方面，对于距离 $>\tau$ 的点，将获得一个更小些的权重。事实上，在最优化计算中，此权重函数为那些距离远的点选定了其距离值而不是平方距离值参加运算。有时，选定的权重值不足够的小，所以不能完抑制离群值。这种情况下，可以使用 Tukey 权重函数，其定义为：
$$w(\delta )=\{\begin{array}{l}{(1-(\frac{\delta }{\tau }{)}^2)}^2,|\delta |\le \tau \\ 0,|\delta |>\tau \end{array}$$

上式中的 $\tau$ 也是削波因数。注意此权重函数完全忽略那些距离 $>\tau$ 的点。而对于 $\le \tau$ 的点，其权重值在 $1$ 到 $0$ 之间平滑变化。

削波因数的计算

在上面的两个权重函数中，削波因数指明哪些点应该被视为离群值。因为削波因数代表的是一个距离，所以它可以被手动设置。但是，这将忽视噪声的分布和数据中的离群值，所以不得不针对每个 应用来进行调整。从数据自身来推出削波因数会更方便。这通常可以基于到图形的这些距离值的标准偏差来实现。因为我们想得到的是数据中离群值，所以我们不能使用正规的标准偏差，而n必须使用对于离群值是可靠的标准偏差。一般情况下，下式用来计算这个可靠的标准偏差：
$${\sigma }_{\delta }=\frac{median|{\delta }_i|}{0.6745}$$

对于正态分布的距离值,分母上的常数量被选定为适用于基于正态分布计算的标准偏差。削波因子则被设置为 ${\sigma }_{\delta }$ 的一个小倍数，如 $\tau =2\sigma \delta$。

Huber 函数

在统计学习角度，Huber损失函数是一种使用鲁棒性回归的损失函数，它相比均方误差来说，它对异常值不敏感。常常被用于分类问题上。

Huber 函数的定义:

只给出了定义，没有给出原因，且论文中的文献为相关理论的图书，若有需要，可详细阅读

$${\rho }_H=\{\begin{array}{l}\frac{1}{2}{e}^2,|e|\le k,\\ k|e|-\frac{1}{2}{k}^2,|k|>e\end{array}$$

Huber 损失函数常常用于鲁棒性系统分析，M元估计和适应性建模。

下面给出最小二乘法(LS)、Huber 和 Tukey Bisquare 的函数及其权重函数对应表

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-AMhJ2wbD-1584789138217)(图3.png)]

上图中，所示的函数图像、导函数图像及权重函数图像，如下图所示：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-HYfMpaQY-1584789138218)(图4.png)]

如上图所示，从左到，分别为对应的函数图像、导函数图像及权重函数图像。

M 估计

M 估计是基于最小二乘估计发展起来的一种抗差估计(Robust Estimation)方法。

M 估计的发展

理论研究和实践经验表明高斯一马尔可夫模型的最小二乘估计,存在不足不稳健,也就是粗差发生时,参数的最小二乘估计不可靠,与其真值偏离太远。误差分布不是正态分布时,最小二乘估计不是最优估计。针对问题,最小二乘法的创立者之一,法国科学家勒让得指出,在使用最小二乘法时,要注意那些大的偏离值粗差可见,粗差问题不是现代才注意到,但其系统研究却是始于世纪年代,并被定义为估计,年发表的论文“位置参数的稳健估计”是其标志。在年,他将此方法用于一般的多元线性模型,自那以后,这个领域不仅受到统计学家的重视,更受到应用工作者的重视和欢迎,在测绘科学也不例外,几十年来,理论和应用上产生了一批重要成果并且对问题当误差不是正态分布时,在估计中可以找到优于最小二乘法的估计。

尽管估计的初衷是为了解决最小二乘法的不稳健问题,但从它的发展来看永远超越了它最初的目的,实际上它是包括了稳健估计、最小二乘估计在内的一个广义的估计类。

M 估计的极值法定义

对独立同分布等精度模型：
$$L=AX-\Delta$$
$$D(\Delta )={\delta }^{2}I$$

选定一个定义于一维欧氏空间 $R$ 的实函数 $\rho$,令
$$\sum _{i=1}^n\rho (a_i^{T}X-L)$$

式中 $a_i^T$ 是设计矩阵的行向量：$X$ 又是极值解这种估计就称为M估计,是 Huber 于1964年在发表的论文“位置参数的稳健估计”中引进的位置参数模型就是上式中只有一个未知参数、设计矩阵,A=（1,1…1）的情况。1973年,Huber 进一步将这种估计拓展到独立同分布等精度的线性函数模型。应用中,我们常常面临独立同分布不等精度模型,由此拓展上式为：
$$\sum _{i=1}^n\rho (P_i^{\frac{1}{2}}(a_i^{T}X-L))=min\sum _{i=1}^n\rho (P_i^{\frac{1}{2}}(a_i^{T}X-L))$$

式中 $P$ 是相应观测量或观测误差的权。

Huber引入估计的动机是追求参数估计的稳健性,这是最小二乘估计所缺乏的。但要注意到,稳健性只是参数估计的一个性质,这个性质与最优性不同,我们不能说愈稳健就愈好。对基于正态分布的参数估计,过于强调稳健性会导致效率损失。

几何拟合(geometric fitting)

几何拟合(也称为最佳拟合)，其误差距离定义为从给定点到要拟合的几何特征的正交距离或最短距离。几何拟合可能是上述所有代数拟合问题的唯一解决方案。各种误差量度的拟合优度测试已将几何误差量度用作基准。

圆的几何拟合是一个非线性问题，必须用迭代法来进行求解。

雅克比矩阵

在向量微积分中，雅可比矩阵是一阶偏导数以一定方式排列成的矩阵，其行列式称为雅可比行列式。雅可比矩阵的重要性在于它体现了一个可微方程与给出点的最优线性逼近。因此，雅可比矩阵类似于多元函数的导数。

雅克比矩阵的定义

假设 $F:{\mathbb{R}}_n\to {\mathbb{R}}_m$ ，是一个从 $n$ 维欧氏空间映射到到 $m$ 维欧氏空间的函数。这个函数由 $m$ 个实数组成：
$$\{\begin{array}{l}{y}_1=f(x_1,x_2,\cdots ,x_n),\\ y_2=f(x_1,x_2,\cdots ,x_n),\\ \cdots \\ y_m=f(x_1,x_2,\cdots ,x_n)\end{array}$$
这些函数的偏导数(如果存在)可以组成一个 $m$ 行 $n$ 列的矩阵，这个矩阵就是所谓的雅可比矩阵：

$$J=\left[\begin{array}{cccc}\frac{\partial f}{\partial x_1}& \frac{\partial f}{\partial x_2}& \cdots & \frac{\partial f}{\partial x_n}\end{array}\right]=\left[\begin{array}{ccc}\frac{\partial y_1}{\partial x_1}& \cdots & \frac{\partial y_1}{\partial x_n}\\ ⋮& \ddots & ⋮\\ \frac{\partial y_m}{\partial x_1}& \cdots & \frac{\partial y_m}{\partial x_n}\end{array}\right]$$

向量微积分

把向量中的对应位置的元素进行求导，但是，这些元素的组织方式有两种，分别是分子布局和分母布局，二者并无本质上的差别，只是结果相差个转置。

例如求$\frac{\partial y}{\partial x}$，其中 $y$ 是 $n$ 维列向量，$x$ 是标量。这个求导就是把 $y$里每个元素分别对 $y$ 求导，但求导后是得到列向量还是行向量呢？

对于分子布局：

$$\frac{\partial y}{\partial x}=\left[\begin{array}{c}\frac{\partial y_1}{\partial x}\\ \\ \frac{\partial y_2}{\partial x}\\ ⋮\\ \frac{\partial y_n}{\partial x}\end{array}\right]$$

对于分子布局：

$$\frac{\partial y}{\partial x}=\left[\begin{array}{cccc}\frac{\partial y_1}{\partial x}& \frac{\partial y_2}{\partial x}& \cdots & \frac{\partial y_n}{\partial x}\end{array}\right]$$

其中输入向量 $\vec{x}=(x_1,\cdots ,x_n)$，输出向量 $\vec{y}=(y_1,\cdots ,y_m)$，$x \in \mathbb{R_n} $，$ y = f(x) \in \mathbb{R_m}$。

如果 $p$ 是中的一个点，$F$ 在 $p$ 点可微分，根据数学分析，雅可比矩阵是在这点的导数。在此情况下，这个线性映射，即 $F$ 在点 $p$ 附近的最优线性逼近，也就是说当 $x$ 足够靠近点 $p$ 时，我们有:
$$F(x)\approx F(p)+J_F(p)\cdot (x-p)$$

高斯 - 牛顿迭代法

高斯- 牛顿迭代法的定义

高斯一牛顿迭代法(Gauss-Newton iteration method)是非线性回归模型中求回归参数进行最小二乘的一种迭代方法，该法使用泰勒级数展开式去近似地代替非线性回归模型，然后通过多次迭代，多次修正回归系数，使回归系数不断逼近非线性回归模型的最佳回归系数，最后使原模型的残差平方和达到最小。其直观思想是先选取一个参数向量的参数值 $\beta$，若函数 $f_t(X_t,\beta )$ 在 ${\beta }_0$ 附近有连续二阶偏导数，则在 ${\beta }_0$ 的邻域内可近似地将 $f_t(X_t,\beta )$看作是线性，因而可近似地用线性最小二乘法求解。

高斯-牛顿迭代的求解步骤

高斯-牛顿法的一般步骤,如下所示:

1. 初始值的选择

   初始值的选择有三种方法：

   • 根据以往的经验选定初始值；
   • 用分段法求出初始值；
   • 对于可线性化的非线性回归模型，通过线性变换，然后实行最小平方法求出初始值。

2. 泰勒级数展开式

   设非线性回归模型为
   $$\hat{y}=f(x_i,r)+{\epsilon }_i,(i=1,2,\cdots n)\text{\hspace{1em}(1)}$$
   其中，$r$ 为待估计回归系数，误差项 ${\epsilon }_i\sim N(0,{\sigma }^2)$。

   设 $g^0(g_0^{(0)},g_1^{(0)},g_2^{(0)},\cdots ,g_{p-1}^{(0)}{)}^T$ 为待估回归系数 $r=(r_0,r_1,r_2,\cdots ,r_{p-1}{)}^T$ 的初始值，将式$(1)f(x_i,r)$ 在 $g$ 点附近的泰勒展开，并略去非线性回归模型的二阶及二阶以上的偏导数项，得:
   $$f(x_i,r)=f(x_i,g^{(0)})+\sum _{k=0}^{(p-1)}[\frac{\partial f(x_i,r)}{\partial r_k}{]}_{r=g^{(0)}}(r_k-g_k^{(0)})\text{\hspace{1em}(2)}$$
   将式 $(2)$ 代人式 $(1)$，则
   $$y_i\approx f(x_i,g^{(0)})+\sum _{k=0}^{(p-1)}[\frac{\partial f(x_i,r)}{\partial r_k}{]}_{r=g^{(0)}}(r_k-g_k^{(0)})+{\epsilon }_i$$

   移项，有
   $$y_i-f(x_i,g^{(0)})\approx \sum _{k=0}^{(p-1)}[\frac{\partial f(x_i,r)}{\partial r_k}{]}_{r=g^{(0)}}(r_k-g_k^{(0)})+{\epsilon }_i$$
   令， $y_i^0=y_i-f(x_i,g^{(0)}),D_{ik}^{(0)}=[\frac{\partial f(x_i,r)}{\partial r_k}{]}_{r=g^{(0)}}，b_k^{(0)}=r_k-g_k^{(0)}$ 则,
   $$y_i^0\approx \sum _{k=0}^{(p-1)}{D}_{(ik)}^{(0)}{b}_k^{(0)}+{\epsilon }_i(i=1,2,\cdots ,n)$$

   用矩阵形式表示，则上式为：
   $$Y^{(0)}\approx D^{(0)}{b}^{(0)}+\epsilon ,\text{\hspace{1em}(3)}$$
   其中，
   $$Y_{n\times p}^{(0)}=\left[\begin{array}{c}{y}_1-f(x_1,g^{(0)})\\ \cdots \\ y_n-f(x_n,g^{(0)})\end{array}\right],D_{n\times p}^0=\left[\begin{array}{ccc}{D}_{10}^{(0)}& \cdots & D_{1(p-1)}^{(0)}\\ ⋮& & ⋮\\ D_{(n0)}^{(0)}& \cdots & D_{n(p-1)}^{(0)}\end{array}\right],D_{p\times 1}^0=\left[\begin{array}{c}{b}_0^{(0)}\\ ⋮\\ b_{p-1}^{(0)}\end{array}\right]$$

3. 估计修正因子

   用最小二乘法对式 $(3)$ 估计修正则
   $$b^{(0)}=(D^{(0)T}{D}^{(0)}{)}^{-1}{D}^{(0)T}{Y}^{(0)},\text{\hspace{1em}(4)}$$
   设 $g^{(0)}$为第一个迭代值，则$g^{(1)}=g^{(0)}+b^{(0)}$。

4. 精度的检验
   设残差平方和为：
   $$SS{R}^{(s)}=\sum _{i=0}^n[y_i-f(x_i,g^{(s)}){]}^2$$
   其中，$s$ 内重夏迭代次数。对于给定的允许误差率$K$，当$|\frac{SS{R}^{(s)}-SS{R}^{(s-1)}}{SS{R}^{(s)}}|\le K$ 时，终止迭代，否则，对式$(4)$做下一次迭代。

5. 重复迭代
   重复式$(4)$，当重复迭代 $s$ 次时，则有修正因子 $b^{(s)}=(D^{(s)T}{D}^{(s)}{)}^{-1}{D}^{(s)T}{Y}^{(s)}$，第 $(s+1)$ 次迭代值$g^{(s+1)}=g^{(s)}+b^{(s)}$。

泰勒公式

泰勒公式，应用于数学、物理领域，是一个用函数在某点的信息描述其附近取值的公式。如果函数足够平滑的话，在已知函数在某一点的各阶导数值的情况之下，泰勒公式可以用这些导数值做系数构建一个多项式来近似函数在这一点的邻域中的值。泰勒公式还给出了这个多项式和实际的函数值之间的偏差。

泰勒公式的定义

泰勒公式是将一个在 $x=x_0$ 处具有 $n$ 阶导数的函数 $f(x)$ 利用关于 $(x-x_0)$ 的 $n$ 次多项式来逼近函数的方法。

若函数 $f(x)$ 在包含 $x_0$ 的某个闭区间 $[a,b]$ 上具有 $n$ 阶导数，且在开区间 $(a,b)$ 上具有 $(n+1)$ 阶导数，则对闭区间 $[a,b]$ 上任意一点x，成立下式：
$$f(x)=\frac{x_0}{0!}+\frac{f^{\prime }(x_0)}{1!}(x-x_0)+\frac{f^{\prime \prime }(x_0)}{2!}(x-x_0)+\cdots +\frac{f^{(n)}(x_0)}{n!}+R_n(x)$$

其中, $f^{(n)}(x)$ 表示 $f(x)$ 的 $n$ 阶导数，等号后的多项式称为函数 $f(x)$ 在$x_0$ 处的泰勒展开式，剩余的 $R_n(x)$ 是泰勒公式的余项，是 $(x-x_0{)}^n$ 的高阶无穷小。

泰勒公式的余项 $R_n(x) $，可以写成以下几种形式：

1. 佩亚诺(Peano)余项：
   $$R_n(x)=o[(x-x_n{)}^n]$$

   这里只需要n阶导数存在。

2. 施勒米尔希-罗什(Schlomilch-Roche)余项：
   $$R_n(x)=f^{(n+1)}[x_0+\theta (x-x_0)]\frac{(1-\theta {)}^{n+1-p}(x-x_0{)}^{n+1}}{n!p}$$

   其中, $\theta \in (0,1)$，$p$ 为任意正实数。(注意到 $p=n+1$ 与 $p=1$ 分别对应拉格朗日余项与柯西余项)

3. 拉格朗日(Lagrange)余项：
   $$R_n(x)=f^{(n+1)}[x_0+\theta (x-x_0)]\frac{(x-x_0{)}^{n+1}}{(n+1)!}$$

   其中, $\theta \in (0,1)$。

4. 柯西(Cauchy)余项：
   $$R_n(x)=f^{(n+1)}[x_0+\theta (x-x_0)]\frac{(x-x_0{)}^{n+1}}{n!}$$
   其中, $\theta \in (0,1)$。

5. 积分余项
   $$R_n(x)=\frac{(-1{)}^n}{n!}{\int }_a^x(t-x{)}^n{f}^{(n+1)}(t)dt$$

   其中，以上诸多余项事实上很多是等价的。

泰勒级数

按照上述泰勒公式，如果f(x)在x0处无限可导，那么泰勒公式则变为
$$f(x)=\sum _{n=0}^{\infty }\frac{f^{(n)}}{n!}(x-x_0{)}^n+R_{\infty }$$

其中幂级数(Power Series)
$$\sum _{n=0}^{\infty }\frac{f^{(n)}}{n!}(x-x_0{)}^n$$
称为 $f(x)$ 在点 $x_0$ 处的泰勒级数。

算法原理：用高斯-牛顿迭代法，计算非线性方程组中雅克比矩阵的求解。

非线性最小二乘拟合

TODO:完善公式的推导过程

假设向量$a$的参数 $q$ 与 $p(p\ge q)$ 的在 $X$ 中的测量有关，
$$X=F(\hat{\vec{a}})+\vec{e}\text{\hspace{1em}(1)}$$

其中，$F$ 是 $a$ 的非线性连续函数，$e$ 表示误差的零平均值。
针对 $a$，给定的 $X$ 的非线性最小二乘估计必须最小化性能指标为：
$${\delta }_0^2=[X-F(\hat{\vec{a}}){]}^T[X-F(\hat{\vec{a}})]\text{\hspace{1em}(2)}$$
为了方便起见，将权重矩阵(weight matrix)(或噪声协方差矩阵(noise covariance matrix))选用单位矩阵(单位矩阵是一个主对角线元素都为 1 的方阵)；

根据高斯-牛顿迭代，将初始值设置为$\overrightarrow{a_0}$，步长参数设置为 $\lambda$ 则
$$\frac{\partial F}{\partial \vec{a}}{|}_{a_k}\Delta \vec{a}=X-F(\vec{a})\text{\hspace{1em}(3)}$$

$${\vec{a}}_{k+1}={\vec{a}}_k+\lambda \Delta \vec{a}\text{\hspace{1em}(4)}$$

其中，等式(3)的雅克比矩阵，如下所示：
$$J_{ij}=\frac{\partial F_i}{\partial a_j}\text{\hspace{1em}(5)}$$

对于上述类型的数值解，必须提供向量 $\vec{F}$ 的函数值和雅克比矩阵 $J$。在提供了向量 $\vec{F}$ 的函数值和雅克比矩阵 $J$后，可使用 $SVD$ 对雅克比矩阵进行分解。
$$\frac{\partial \vec{F}}{\partial \vec{a}}{|}_{a_k}=UW{V}^T$$
$$U^{T}U=V^{T}V=I,W=(diag(w_1,w_2,\cdots ,w_q))\text{\hspace{1em}(6)}$$

根据上述公式，完成迭代后，即可得到公式(2)中的${\delta }^2$。
协方差矩阵的公式，如下所示：
$$Cov(\vec{a})=(J^{T}J{)}^{-1}=(VW{U}^{T}UW{V}^T{)}^{-1}=V{W}^{-2}{V}^T$$
$$Cov(a_j,a_k)=\sum _{i=1}^q\left(\frac{V_{ji}{V}_{ki}}{w_i^2}\right),i=1,2,\cdots ,q.\text{\hspace{1em}(7)}$$
因此，预估参数的协方差可以写成以下形式：
$${\delta }^2(a_j)=\frac{{\delta }^2}{p-q}Cov(a_j,a_j),j=1,2,\cdots ,q\text{\hspace{1em}(8)}$$
相关系数为：
$$\rho (a_j,a_k)=\frac{Cov(a_j,a_k)}{\sqrt{Cov(a_j,a_j)Cov(a_k,a_k)}},j=1,2,\cdots ,q,k=1,2,\cdots ,q.$$

圆的几何拟合

TODO: 完成公式的理论推导过程

圆心在$X_c$ 处, 且半径为$R$ 圆或球，在 $n(n\ge 2)$ 维空间，可表示为以下形式：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-v0NUhZ5W-1584789138219)(图2.png)]

根据公式，$R^2=||X-X_c|{|}^2$,其中，$X$ 为拟合圆上的点，则点 $X_i(\ne X_c)$ 对应的最近的点的坐标:
$$X_i^{\prime }=X_c+R\frac{X_i-X_c}{||X_i-X_c||}(i=1,2,\cdots ,m)\text{\hspace{1em}(9)}$$
$X_i^{\prime }$ 的误差距离的正交向量为：
$$X_i^{\prime \prime }=X_i-X_i^{\prime }=[||X_i-X_c||-R]\frac{X_i-X_c}{||X_i-X_c||}\text{\hspace{1em}(10)}$$
根据公式(9),可以推导出，雅克比矩阵 $J$ 的高斯-牛顿迭代公式为：
$$J_{X_i^{\prime },\vec{a}}=\frac{\partial }{\partial \vec{a}}\left[X_c+R\frac{X_i-X_c}{||X_i-X_c||}\right]$$
$$=\frac{\partial X_c}{\partial \vec{a}}+\frac{X_i-X_c}{||X_i-X_c||}\frac{\partial R}{\partial \vec{a}}-\frac{R}{||X_i-X_c||}$$
$$\times \left[I-\frac{(X_i-X_c)(X_i-X_c{)}^T}{||X_i-X_c|{|}^2}\right]\frac{\partial X_c}{\partial \vec{a}}$$

注：论文中，以球体的拟合为例：

我们定义参数向量 $\vec{a}$, 则：
$${\vec{a}}^T=(R,X_c^T)=(R,X_c,Y_c,Z_c),\text{\hspace{1em}(12)}$$
对上式中的 $R$ 与 $X_c$ 求偏导数，得：
$$\frac{\partial R}{\partial \vec{a}}=(1|{\vec{0}}^T)和\frac{\partial X_c}{\partial \vec{a}}=(0|\vec{I})\text{\hspace{1em}(13)}$$

根据，雅克比矩阵 $J_{X_i^{\prime },a}$ 和误差距离向量 $X_i^{\prime \prime }$,我们可以构建出 $m$ 个点在 $n$ 维空间中的线性关系：$p=m\cdot n$, 因为是球体的拟合，所以,$(n=3,p=3m,q=4)$,因此，该线性等式，如下所示：
$$\left(\begin{array}{cccc}{J}_{X_1^{\prime },R}& J_{X_1^{\prime },X_c}& J_{X_1^{\prime },Y_c}& J_{X_1^{\prime },Z_c}\\ J_{Y_1^{\prime },R}& J_{Y_1^{\prime },X_c}& J_{Y_1^{\prime },Y_c}& J_{Y_1^{\prime },Z_c}\\ J_{Z_1^{\prime },R}& J_{Z_1^{\prime },X_c}& J_{Z_1^{\prime },Y_c}& J_{Z_1^{\prime },Z_c}\\ ⋮& ⋮& ⋮& ⋮\\ J_{X_m^{\prime },R}& J_{X_m^{\prime },X_c}& J_{X_m^{\prime },Y_c}& J_{X_m^{\prime },Z_c}\\ J_{Y_m^{\prime },R}& J_{Y_m^{\prime },X_c}& J_{Y_m^{\prime },Y_c}& J_{Y_m^{\prime },Z_c}\\ J_{Z_m^{\prime },R}& J_{Z_m^{\prime },X_c}& J_{Z_m^{\prime },Y_c}& J_{Z_m^{\prime },Z_c}\end{array}\right)\left(\begin{array}{c}\Delta R\\ \Delta X_c\\ \Delta Y_c\\ \Delta Z_c\end{array}\right)=\left(\begin{array}{c}{X}_1^{\prime \prime }\\ Y_1^{\prime \prime }\\ Z_1^{\prime \prime }\\ ⋮\\ X_m^{\prime \prime }\\ Y_m^{\prime \prime }\\ Z_m^{\prime \prime }\end{array}\right)\text{\hspace{1em}(14)}$$

上述过程中，高斯-牛顿迭代初始参数向量，可以由代数拟合圆提供，也可以通过重心和RMS(Root Mean Square) 获取，计算公式如下所示：
$$X_c=\overline{X}=\frac{1}{m}\sum _{i=1}^m{X}_i,R=\sqrt{\frac{1}{m}\sum _{i=1}^m|X_i-\overline{X}{|}^2}\text{\hspace{1em}(15)}$$

RANSAC (Random Sample Consensus)拟合

随机抽样一致算法（RANdom SAmple Consensus,RANSAC）,采用迭代的方式从一组包含离群的被观测数据中估算出数学模型的参数。RANSAC算法假设数据中包含正确数据和异常数据（或称为噪声）。正确数据记为内点（inliers），异常数据记为外点（outliers）。同时RANSAC也假设，给定一组正确的数据，存在可以计算出符合这些数据的模型参数的方法。该算法核心思想就是随机性和假设性，随机性是根据正确数据出现概率去随机选取抽样数据，根据大数定律，随机性模拟可以近似得到正确结果。假设性是假设选取出的抽样数据都是正确数据，然后用这些正确数据通过问题满足的模型，去计算其他点，然后对这次结果进行一个评分。

RANSAC算法被广泛应用在计算机视觉领域和数学领域，例如直线拟合、平面拟合、计算图像或点云间的变换矩阵、计算基础矩阵等方面，使用的非常多。

RANSAC的基本假设是：
（1）数据由“局内点”组成，例如：数据的分布可以用一些模型参数来解释；
（2）“局外点”是不能适应该模型的数据；
（3）除此之外的数据属于噪声。
局外点产生的原因有：噪声的极值；错误的测量方法；对数据的错误假设。

RANSAC 算法的流程

RANSAC 的算法流程，如下所述:

1. 随机地从数据集 $S$ 中，选取 $s$ 个数据点，一般取 $s$ 刚好能够组成一个最小的子集，以 $s$ 的数据点构成待求模型，比如直线拟合的实验中，最小的子集就是两个点，因为，只要取两个点就可以构成一条直线。

2. 根据上一步得到模型，我们可以遍历 $S$ 集合中的每一个数据点，判断它们是否满足这个模型，也就是说，是否在距离阈值 $t$ 内。如果在，则这个点就是内点，将所有这些内点集合，称为一致集 $S_i$。

3. 如果 $S_i$ 的大小(比如，内点的数目) 大于某个阈值 $T$，就用 $S_i$ 中所有点重新估计模型并结束。

4. 如果 $S_i$ 小于阈值 $T$, 就选择一个新的最小子集 $s$，并重复上述过程。

5. 经过 $N$ 次实验选择一个最大一致集 $S_i$， 并用 $S_i$ 的所有点重新估计模型。

RANSAC 与最小二乘的区别

最小二乘法尽量去适应包括局外点的内在的所有点。相反，RANSAC 能得出一个仅仅用局内点计算出的模型，并且概率还足够高。但是，RANSAC 并不能保证结果一定正确，为了保证算法有足够高的合理概率，必须小心的选择算法的参数(参数配置)。经实验验证，对于包含 80% 误差的数据集，RANSAC 的效果，远优于直接的最小二乘法。

验证思路

RANSAC算法的输入是一组观测数据，一个可以解释或者适应于观测数据的参数化模型，一些可信的参数。RANSAC通过反复选择数据中的一组随机子集来达成目标。被选取的子集被假设为局内点，并用下述方法进行验证：

1.有一个模型适应于假设的局内点，即所有的未知参数都能从假设的局内点计算得出。

2.用1中得到的模型去测试所有的其它数据，如果某个点适用于估计的模型，认为它也是局内点。

3.如果有足够多的点被归类为假设的局内点，那么估计的模型就足够合理。

4.然后，用所有假设的局内点去重新估计模型，因为它仅仅被初始的假设局内点估计过。

5.最后，通过估计局内点与模型的错误率来评估模型。

这个过程被重复执行固定的次数，每次产生的模型要么因为局内点太少而被舍弃，要么因为比现有的模型更好而被选用。

参数选择

根据特定的问题和数据集通过实验来确定参数 $t$ 和 $d$。然而参数 $k$（迭代次数）可以从理论结果推断。当估计模型参数时，用 $p$ 表示一些迭代过程中从数据集内随机选取出的点均为局内点的概率；此时，结果模型很可能有用，因此, $p$ 也表征了算法产生有用结果的概率。用 w表示每次从数据集中选取一个局内点的概率，如下式所示：
w = 局内点的数目 / 数据集的数目
通常情况下，我们事先并不知道w的值，但是可以给出一些鲁棒的值。假设估计模型需要选定 $n$个点，$w^n$ 是所有 $n$ 个点均为局内点的概率；$1-w^n$是 $n$ 个点中至少有一个点为局外点的概率，此时表明我们从数据集中估计出了一个不好的模型。 $(1-w^n{)}^k$ 表示算法永远都不会选择到 $n$ 个点均为局内点的概率，它和 $1-p$ 相同。因此，
$$1-p=(1-w^n{)}^k$$

我们对上式的两边取对数，得出
$$k=\frac{\log (1-p)}{log(1-w^n)}$$

值得注意的是，这个结果假设 $n$ 个点都是独立选择的；也就是说，某个点被选定之后，它可能会被后续的迭代过程重复选定到。这种方法通常都不合理，由此推导出的k值被看作是选取不重复点的上限。例如，要从数据集寻找适合的直线，RANSAC算法通常在每次迭代时选取2个点，计算通过这两点的直线 maybe_model，要求这两点必须唯一。
为了得到更可信的参数，标准偏差或它的乘积可以被加到 $k$ 上。$k$的标准偏差定义为：
$$SD(k)=\frac{\sqrt{1-w^n}}{w^n}$$

RANSAC 的数学原理

RANSAC的数学原理如下：对于数量为N的样本集，其中正确数量为 $M$，则每一次获取到正确数据的概率 $p=\frac{M}{N}$，计算模型参数需要 $K$ 个数据，则一次性得到 $K$ 个正确数据的概率为 $p^K$（概率很低）；反之至少存在一个错误数据的概率为$(1-p^K)$。若重复抽取 $L$ 次，则 $L$ 组数据中都存在错误数据的概率为$(1-p^K{)}^L$，注意这里当$L$足够大时，这个概率就会变得很小，也就意味着$L$组数据中至少有一组正确数据的概率为$1-(1-p^K{)}^L$，会变得足够大。直观解释就是抽取次数越多，抽到一组正确数据的概率就会越大。

一般RANSAC的实现代码中有一个输入参数称为置信度 $P$，$P=1-（1-p^K{）}^L$，我们可以通过P反算得到抽取次数$$，$L=\frac{\log(1-P)}{\log(1-p^K)}$。

优点与缺点

优点

RANSAC 的优点是它能鲁棒的估计出模型参数，对噪声具有一定容忍度；例如，它能从包含大量局外点的数据集中估计出高精度的参数。

缺点

1. RANSAC 计算参数的迭代次数没有上限；如果设置迭代次数的上限，得到的结果可能不是最优的结果，甚至可能得到错误的结果。

2. RANSAC只有一定的概率得到可信的模型，概率与迭代次数成正比。

3. RANSAC的另一个缺点是它要求设置跟问题相关的阀值。

4. RANSAC只能从特定的数据集中估计出一个模型，如果存在两个（或多个）模型，RANSAC不能找到别的模型。

应用

计算机视觉中RANSAC经常用于多视图间计算基本矩阵或者单应矩阵的计算。如由于镜头的限制，往往需要多张照片才能拍下那种巨幅的风景。在多幅图像合成时，事先会在待合成的图片中提取一些关键的特征点。计算机视觉的研究表明，不同视角下物体往往可以通过一个透视矩（3X3或2X2）阵的变换而得到。RANSAC被用于拟合这个模型的参数（矩阵各行列的值），由此便可识别出不同照片中的同一物体。$$
$K$ 个正确数据的概率为 $p^K$（概率很低）；反之至少存在一个错误数据的概率为$(1-p^K)$。若重复抽取 $L$ 次，则 $L$ 组数据中都存在错误数据的概率为$(1-p^K{)}^L$，注意这里当$L$足够大时，这个概率就会变得很小，也就意味着$L$组数据中至少有一组正确数据的概率为$1-(1-p^K{)}^L$，会变得足够大。直观解释就是抽取次数越多，抽到一组正确数据的概率就会越大。

一般RANSAC的实现代码中有一个输入参数称为置信度 $P$，$P=1-（1-p^K{）}^L$，我们可以通过P反算得到抽取次数，$L=\frac{\log (1-P)}{\log (1-p^K)}$。

优点与缺点

优点

RANSAC 的优点是它能鲁棒的估计出模型参数，对噪声具有一定容忍度；例如，它能从包含大量局外点的数据集中估计出高精度的参数。

缺点

1. RANSAC 计算参数的迭代次数没有上限；如果设置迭代次数的上限，得到的结果可能不是最优的结果，甚至可能得到错误的结果。

2. RANSAC只有一定的概率得到可信的模型，概率与迭代次数成正比。

3. RANSAC的另一个缺点是它要求设置跟问题相关的阀值。

4. RANSAC只能从特定的数据集中估计出一个模型，如果存在两个（或多个）模型，RANSAC不能找到别的模型。

应用

计算机视觉中RANSAC经常用于多视图间计算基本矩阵或者单应矩阵的计算。如由于镜头的限制，往往需要多张照片才能拍下那种巨幅的风景。在多幅图像合成时，事先会在待合成的图片中提取一些关键的特征点。计算机视觉的研究表明，不同视角下物体往往可以通过一个透视矩（3X3或2X2）阵的变换而得到。RANSAC被用于拟合这个模型的参数（矩阵各行列的值），由此便可识别出不同照片中的同一物体。