论文阅读——椭圆检测 2016 Robust ellipse detection with Gaussian mixture models
这篇文章是16年发表的椭圆检测文章，论文题目为：《Robust ellipse detection with Gaussian mixture models》，发表在《Pattern Recognition》（2区SCI）上。这里最为新颖的地方就是使用高斯混合模型GMM算法进行椭圆检测。下面我就对这篇文章进行分析。
注：
① 2019-1-24 更新：
在明阳师弟的辛苦努力下，终于联系到作者，原版的程序因为作者已经毕业，无法获得，作者自己花了2周时间进行了简单的复现，这个程序主要展示拟合的优点，与论文的方法有一定不同，暂时并不能直接在图像数据集上使用，后期我会在此基础上进行修改，后续我会进一步的进行改进。
GitHub下载地址：https://github.com/clarella/L2-Ellipse-Fitting ，
网盘下载地址：链接：https://pan.baidu.com/s/15aIlNyHHC3wYNdWFRmTrmQ 提取码：ze3v ，在文件夹“Robust ellipse detection with Gaussian mixture models”中。
CSDN下载：https://download.csdn.net/download/zhaoxi_li/10935361 。
文章目录
- 〇 摘要部分
- 一 介绍
- 二 使用L2\mathcal{L}_2L2的椭圆检测
- 2.1 对椭圆gθg_\thetagθ进行建模
- 2.2 对观测点fff进行建模
- 2.3 密度函数的简化表示
- 三 椭圆检测算法
- 3.1 退火算法用于椭圆检测
- 3.2 多个椭圆的检测
- 四 GMM模型的维数增广
- 五 实验部分
- 六 个人总结
〇 摘要部分
高斯混合之间的欧式距离已经被证明在进行点集配准上是鲁棒的。作者采用这个思想用于匹配一系列形状（椭圆）。使用退火策略执行优化，并且多次重复寻找来检测感兴趣中的多个实例1。
一 介绍
介绍这部分第一段只是泛泛介绍想基于配准的算法的思想，没有阅读过这篇文章的肯定读的一头雾水，所以直接先看创新点是啥。相关研究这里也不说了，同类文章大同小异，具体的可以参考前面介绍的椭圆检测文章。这里只想了解是如何将GMM应用在椭圆检测上。创新点部分一共有三点，如下所示：
-
扩展了基于L2\mathcal{L}_2L2用于估计一系列曲线参数的框架。
-
提出了一种用于密度函数的多维模型，以便包含其他可用的信息，而且具有很少的计算。
-
作者提出了一个检测多个椭圆的方法。这个方法用于在2D点云和图片中中进行椭圆检测。
二 使用L2\mathcal{L}_2L2的椭圆检测
高斯混合模型用于表示观测点（定义为fff）和椭圆模型（定义为gθg_{\theta}gθ）之间的形状，其中θ\thetaθ就是椭圆的参数(γ,a,b,xo,yo)(\gamma,a,b,x_o,y_o)(γ,a,b,xo,yo)，γ\gammaγ是椭圆旋转角，a,ba,ba,b是椭圆半短轴和半长轴，xo,yox_o,y_oxo,yo就是椭圆中心点。椭圆的参数可以通过最小化两个密度函数fff和gθg_{\theta}gθ之间的欧式距离来估计。接下来介绍如何对这两个密度函数进行建模，为了方便理解，我会对公式进行详细的推导与分析。
θ^=argminθ{f−gθ}\hat{\theta}=arg\min_{\theta}\left\{f-g_{\theta}\right\}θ^=argθmin{f−gθ}
2.1 对椭圆gθg_\thetagθ进行建模
参数为θ\thetaθ的椭圆上的一点u(j)u^{(j)}u(j)可以由吐下等式表示。其中τj∈[0,2π]\tau_j\in[0,2\pi]τj∈[0,2π]。
u(j)=(cosγ−sinγsinγcosγ)(acosτjbsinτj)+(xoyo)u^{(j)} = \left(\begin{array}{cc} cos\gamma & -sin\gamma \\ sin\gamma & cos\gamma \end{array}\right)\left(\begin{array}{c} acos\tau_j \\ bsin\tau_j \end{array}\right)+\left(\begin{array}{c} x_o \\ y_o \end{array}\right)u(j)=(cosγsinγ−sinγcosγ)(acosτjbsinτj)+(xoyo)
模型gθg_\thetagθ使用均匀分布的N=20N=20N=20个点u(j)u^{(j)}u(j)来定义一个非等方向（non-isotropic）的GMM。
gθ(x)=∑j=1NwjN(x;μj,Σj)g_\theta(x) = \sum_{j=1}^Nw_j\mathcal{N}(x;\mu_j,\Sigma_j)gθ(x)=j=1∑NwjN(x;μj,Σj)
其中N(x;μj,Σj)\mathcal{N}(x;\mu_j,\Sigma_j)N(x;μj,Σj)正态密度函数，x∈R2x\in R^2x∈R2是一个随机变量。其中μj\mu_jμj定义为两个连续定点的均值u(j),u(j+1)u^{(j)},u^{(j+1)}u(j),u(j+1)。协方差矩阵Σj=QjTΛjQj\Sigma_j=Q_j^T\Lambda_jQ_jΣj=QjTΛjQj，其中QjQ_jQj是一个旋转矩阵Qj=[n⃗1j,n⃗2j]Q_j=[\vec{n}_{1j},\vec{n}_{2j}]Qj=[n1j,n2j]，n⃗2j\vec{n}_{2j}n2j与u(j),u(j+1)u^{(j)},u^{(j+1)}u(j),u(j+1)方向相同，是一个单位向量，n⃗1j\vec{n}_{1j}n1j就是n⃗2j\vec{n}_{2j}n2j对应的正交单位向量。Λj\Lambda_jΛj的定义方式如下式，其中htj=∣∣u(j)−u(j+1)∣∣h_{tj} = ||u^{(j)}-u^{(j+1)}||htj=∣∣u(j)−u(j+1)∣∣，参数hhh控制轮廓发现方向的模糊性，这个参数由用户设定（也就是需要调参的参数）。权重参数wjw_jwj与htjhh_{tj}hhtjh成正比，且权重的和为1，利用这个条件可以计算出wjw_jwj实际计算方法，wj=htjh∑i=1Nhtihw_j = \dfrac{h_{tj}h}{\sum_{i=1}^Nh_{ti}h}wj=∑i=1Nhtihhtjh。
Λj=(h200htj2)\Lambda_j=\left(\begin{array}{cc}h^2 & 0 \\ 0 & h_{tj}^2\end{array}\right)Λj=(h200htj2)
下图是不同的采样点个数NNN和不同的模糊阈值hhh对概率的影响，个人认为，边缘点有时候会受到误差影响，有一定的便宜，这个两个参数就是控制这个参数的接受程度吧。NNN过小就形成不了一个严格的山脊，就像图（e）一样，也就是保证椭圆上的区域概率要高。最后选择，N=20，h=1N=20，h=1N=20，h=1。
2.2 对观测点fff进行建模
GMM模型fff对观测点进行建模，依赖可使用的信息和观测点的结构。例如，假设，我们得到一组观测点{v(i)}i=1,...,n\{v^{(i)}\}_{i=1,...,n}{v(i)}i=1,...,n。在这种情况下，没有关于这些点是如何相连的信息。因此，我们使用各向同性协方差矩阵来定义GMM，每个高斯的均值就是其本身。下式就是对应的GMM公式。
f(x)=1c∑i=1nN(x;vi,h2,I)f(x) = \dfrac{1}{c}\sum_{i=1}^n\mathcal{N}(x;v_i,h^2,I)f(x)=c1i=1∑nN(x;vi,h2,I)
我们可以使用边缘点{v(i)}i=1,...,n\{v^{(i)}\}_{i=1,...,n}{v(i)}i=1,...,n作为观测点。除此之外，每个观测点的梯度信息是已知的，是可以用于前面说的GMM模型的。具体的模型定义如下所示。像素梯度信息根据上一小节就可以用于计算Σi\Sigma_iΣi，因为每个边缘点之间的距离很近，因此原来的htjh_{tj}htj在这里就是ht=1h_t=1ht=1。权重wiw_iwi的计算方式同之前的介绍。
f(x)=∑i=1nwiN(x;vi,Σi)f(x) = \sum_{i=1}^nw_i\mathcal{N}(x;v_i,\Sigma_i)f(x)=i=1∑nwiN(x;vi,Σi)
2.3 密度函数的简化表示
两个密度函数的目标函数可以化简为如下所示。
∣∣f−gθ∣∣2=∣∣f∣∣2+∣∣gθ∣∣2−2<f∣gθ>||f-g_\theta||^2 = ||f||^2+||g_\theta||^2-2<f|g_\theta>∣∣f−gθ∣∣2=∣∣f∣∣2+∣∣gθ∣∣2−2<f∣gθ>
∣∣f∣∣2||f||^2∣∣f∣∣2项没有涉及到未知参数θ\thetaθ就不需要进行计算了。其他项的计算应用到了如下的等式。函数的内积就是这两个函数的乘积在其定义域上的积分，下面这个公式是作者参考的那篇文章的引用（等我手头工作完事，就详细去证明这个等式）。
<N(μ1,Σ1)∣N(μ2,Σ2)>=N(0;μ1−μ2,Σ1+Σ2)<\mathcal{N}(\mu_1,\Sigma_1)|\mathcal{N}(\mu_2,\Sigma_2)>=\mathcal{N}(0; \mu_1-\mu_2,\Sigma_1+\Sigma_2)<N(μ1,Σ1)∣N(μ2,Σ2)>=N(0;μ1−μ2,Σ1+Σ2)
之后就可以得到目标函数的两个项的表达形式。这里面的未知参数只有θ\thetaθ。
{∣∣gθ∣∣2=∑j=1N∑k=1NwjwkN(0;μj−μk,Σj+Σk)<f∣gθ>=∑i=1n∑k=1NwiwkN(0;μi−μk,Σi+Σk)\left\{ \begin{array}{l} ||g_{\theta}||^2 = \sum_{j=1}^{N}\sum_{k=1}^{N}w_jw_k\mathcal{N}(0; \mu_j-\mu_k,\Sigma_j+\Sigma_k)\\ <f|g_{\theta}> = \sum_{i=1}^{n}\sum_{k=1}^{N}w_iw_k\mathcal{N}(0; \mu_i-\mu_k,\Sigma_i+\Sigma_k) \end{array}\right.{∣∣gθ∣∣2=∑j=1N∑k=1NwjwkN(0;μj−μk,Σj+Σk)<f∣gθ>=∑i=1n∑k=1NwiwkN(0;μi−μk,Σi+Σk)
三 椭圆检测算法
我们首先展示了一个算法从观测数据中拟合数据。与作者参考的算法相似，在模拟退火框架中进行优化来避免局部解，并限制起始猜测点对优化结果的影响。然后作者有提出了一种用于从观测数据中检测出多个椭圆的策略2。
3.1 退火算法用于椭圆检测
在计算L2\mathcal{L}_2L2时唯一自由的参数是GMM中的fff和gθg_{\theta}gθ。该参数在所提出的估计框架中起着两个重要的作用。首先，它影响形状的描述。它控制曲线在法线方向上的模糊性，其次它影响了代价函数的凸性。使用基于梯度的优化算法（GOA）来执行优化3，这个算法依赖于起始猜测点θ(0)\theta^{(0)}θ(0)和正交带hhh的选择。
θ^=GOA(C(θ),θ(0),h)\hat{\theta} = GOA(\mathcal{C}(\theta),\theta^{(0)},h)θ^=GOA(C(θ),θ(0),h)
hhh的选择越大，损失函数就越平滑。因此，未来使我们的方法对初始点θ(0)\theta^{(0)}θ(0)不明娜，我们使用了一个模拟退火算法框架，其中正交带hhh是随着几何速率的降低而降低（由参数β\betaβ控制）的温度。hhh从hmaxh_{max}hmax开始到小于hminh_{min}hmin停止。模拟退火的使用有助于收敛到全局解。
这个小节说白了，就是使用模糊退火算法进行优化求解，最终得到目标椭圆参数θ\thetaθ。
3.2 多个椭圆的检测
之前的算法都在介绍如何得到一个椭圆，这个小节的目标是从观测点中检测出一组椭圆。作者提出了一个迭代算法来处理这个问题。
-
全局估计： 给定一组观测点，通过最小化损失函数估计出形状参数θ^\hat{\theta}θ^，方法就是前面介绍的拟合出一个椭圆的方法。
-
观测模型fff的更新： 一旦得到一个形状实例θ^\hat{\theta}θ^，所有关联这个形状的观测点将会被移除。剩下的观测点被用于检测其他的椭圆。更新集SfS_fSf应该具有如下的性质，其中Hellipse=1N∑i=1Ngθ^(μj)H_{ellipse}=\dfrac{1}{N}\sum_{i=1}^{N}g_{\hat{\theta}}(\mu_j)Hellipse=N1∑i=1Ngθ^(μj)，t1=0.3t_1=0.3t1=0.3。通俗点讲，HellipseH_{ellipse}Hellipse就是所有观测点对应的形状概率图。小于概率值的就被认为是不关联的。得到SfS_fSf。
vi∈Sf,ifgθ^(vi)<t1×Hellipsev_i\in S_f, ~if~ g_{\hat{\theta}}(v_i)<t_1\times H_{ellipse}vi∈Sf, if gθ^(vi)<t1×Hellipse
关于HellipseH_{ellipse}Hellipse的理解，可以参考下图的介绍。参与构成一个椭圆的观测点在GMM上具有较大的概率，也就是下图中凸起的部分，这个凸起部分的平均值代表整体概率，其他观测点在这个模型上不应该具有较大的概率。这样我们可以得到一堆其他观测点，这些观测点肯定不在这个得到的椭圆上，然后重复检测，重复利用这个想法，直到剩下10%的点结束。
四 GMM模型的维数增广
之前的模型直接采用的每个图像的像素点，所以每个GMM模型都是2维的，作者又利用了梯度角的信息，将GMM模型增光到3维，更新了均值与协方差矩阵的定义，其他的计算方法都不变，目的就是提高检测精度。思想方法都很容易理解。
五 实验部分
作者将自己的方法与一些拟合算法进行了对比，验证了其在噪声较高的地方具有高鲁棒性。但是，我觉得这个问题在椭圆检测领域是不重要的，根据提供的测试图可以看出来，参与拟合的数据都有一些错误值，这些值肯定不属于椭圆。作者提供的这个算法可以理解为一种剔除错误值的算法，所以肯定会比那些没有剔除错误值的方法要好。但是，实际中，椭圆拟合主要就是弧段组合，很少会出现这种有异常值的点，每个点都是通过Canny检测出来的，会有1~3个像素的失真，但像测试图这种我个人认为出现的情况太少了。
从下面的这个仿真图测试来看，其效果不是十分理想，很多看起来很容易检测出来的椭圆没有检测出来。个人认为，观测点找到之后最好进行一次直接拟合，不要进行优化求解，毕竟优化求解不一定能很好的找到好的拟合解。效果不好，不代表这个方法不好，作者是新的尝试，是一种大胆的创新，作者已经证明了其可行性，在未来的工作中如果解决这些问题，一定会比边缘链接的要容易处理，效果更好。
下图是检测一个图片多个椭圆的过程，其方法是有效的，而且大部分是交给算法进行自动识别，很符合机器学习的思想，结果看起来也十分有趣。
六 个人总结
这个文章创新性很好，尽管结果差强人意，但是这个想法非常值得我们去借鉴。但是，从这个算法，仍然有很多地方需要分析研究。
-
简单场景，这个算法会十分实用，而且对于手写的方法一定会有良好的匹配效果，匹配的思想绝对是优于边缘链接的。真实场景的Canny边缘图相当复杂，这个算法的耗时是否增加的非常严重。
-
基于边缘的方法如果边缘被分割的十分严重，效果肯定不好。但这个方法存在大量的错误检测，而且观测点正确，缺拟合出了错误的结果，这后期也要深度研究。
-
观测点的组合，我觉得会有效的提高检测精度。
总而言之，这个算法非常值得我们去阅读，如果我们想采用机器学习的思想去进行椭圆检测，不建议直接去与基于边缘的去对比。个人建议，制作一个手绘椭圆数据集，与本文和UpWrite算法去对比。这样绝对会是很好的创新。
简而言之就是参考了点集配准的思想，从图片上配准椭圆这个形状，并采用退火算法进行优化求解。但是需要注意的是，我之前阅读过CPD点击配准算法，配准的是两个点云，而在椭圆检测中，椭圆像素点占据很小的部分，如何剔除大量无用的点对配准造成的影响非常重要。 ↩︎
这才是关键，如果只有一个椭圆直接拟合就可以了，椭圆拟合算法肯定比这种匹配更加直接有效。 ↩︎
论文中基于梯度的优化算法的缩写为GA。我们知道遗传算法的缩写也是GA，为了避免这种非常严重的模糊，我将GA缩写改为GOA。 ↩︎