HGMMs的一点点理解～

Robust Ellipse Fitting Using Hierarchical Gaussian Mixture Models

Motivation
Related works
Ellipse fitting based on GMMs
Ellipse fitting based on HGMMs
- Distance-Density-Based Error Metric
- Density-Based Region Growing Method
一点点理解

两篇不错的文章，idea 很新颖
原文链接： HGMMs.
原文链接： GMMs.

Motivation

本文是基于论文"Robust ellipse detection with Gaussian mixture models” (GMMs) 的一种改进方法，主要解决了GMMs在面对outliers (离群点)时不稳定的问题。

Related works

在一个二维欧氏空间中一个圆锥曲线可以表示为一种二阶隐函数方程：
Ax2+2Bxy+Cy2+2Dx+2Ey+F=0.Ax^{2}+2Bxy+Cy^2+2Dx+2Ey+F=0. Ax2+2Bxy+Cy2+2Dx+2Ey+F=0.
为了方便矩阵计算一般标识为向量形式Θv=0\Theta \bm{v}=0Θv=0, 其中Θ=[ABCDEF]\Theta=[A\quad B\quad C\quad D\quad E\quad F]Θ=[ABCDEF],v=[x22xyy22x2y1]\bm{v}=[x^2\quad 2xy\quad y^2\quad 2x\quad 2y\quad 1]v=[x22xyy22x2y1]。
这样根据判别式 (discriminant) B2−AC<0B^2-AC<0B2−AC<0可以判断该双曲线是否为椭圆。这样给定一组数据点C={(xi,yi)∈R2}C= \{(x_{i},y_{i}) \in \mathbb{R}^2\}C={(xi,yi)∈R2},需要处理的问题就是根据这些点拟合出一个合适的椭圆。很多方法都是先假设一个初始椭圆Θ0\Theta_0Θ0通过不断迭代地最优化数据点到椭圆的距离来获得最终的椭圆参数，也就是最优化下面两个式子：
J(Θ)=∑i=1n∥ΘTvi∥2,J(\Theta)=\sum^n_{i=1} \|\Theta^T\bm{v_i}\|^2, J(Θ)=i=1∑n∥ΘTvi∥2,
S=∑i=1n((xi−x^i)2−(yi−y^i)2).S=\sum^n_{i=1}((x_i-\hat{x}_i)^2-(y_i-\hat{y}_i)^2). S=i=1∑n((xi−x^i)2−(yi−y^i)2).
上面一个是代数距离，下面一个是几何距离，其中(xi,yi)(x_i,y_i)(xi,yi)是数据点，(x^,y^)(\hat{x},\hat{y})(x^,y^)是椭圆上离(xi,yi)(x_i,y_i)(xi,yi)最近的一点。

Ellipse fitting based on GMMs

与上述方法不同的是，Arellano开创性的在2016年提出通过GMMs对观测点和目标椭圆建模，并根据两个GMMs之间的L2范数来迭代的获得椭圆拟合结果 (虽然论文里说的是detection，但本质上还是椭圆拟合的方法)。给定一组观测点C={xj∈R2}j=1nC=\{\bm{x}_j\in\mathbb{R}^2\}_{j=1}^nC={xj∈R2}j=1n，其GMMs表示为
p(x)=1n∑j=1nKh(x−xj),p(\bm{x})=\frac{1}{n}\sum_{j=1}^{n}K_h(\bm{x}-\bm{x_j}), p(x)=n1j=1∑nKh(x−xj),
其中Kh(x−xj)K_h(\bm{x}-\bm{x_j})Kh(x−xj)是核函数, h表示核带宽（也许可以理解成每个高斯核的宽窄？），当它等于下式时p(x)p(\bm{x})p(x)就是高斯混合模型 (GMMs)
N(x;μ,Σ)=12π∥Σ∥1/2exp⁡{−12(x−μ)TΣ−1(x−μ)},N(\bm x;\bm \mu,\bm \Sigma)=\frac{1}{2\pi \|\bm{\Sigma}\|^{1/2}}\exp\{-\frac{1}{2}(\bm x-\bm \mu)^T\bm \Sigma^{-1}(\bm x-\bm \mu) \}, N(x;μ,Σ)=2π∥Σ∥1/21exp{−21(x−μ)TΣ−1(x−μ)},
由此得出观测值所对应的GMMs为
f(x)=1n∑j=1nN(x;,xj,h2I)。f(\bm x)=\frac{1}{n}\sum_{j=1}^{n}N(\bm x;,\bm x_j, h^2\bm I)。 f(x)=n1j=1∑nN(x;,xj,h2I)。
Ok，到这里对观测值的建模已经完成，但是我们需要得到的是观测值所满足的一种分布关系也就是椭圆的参数，因此作者初选一个椭圆参数θ={x0,y0,a,b,γ}\theta=\{x_0,y_0,a,b,\gamma\}θ={x0,y0,a,b,γ}, 在其上均匀取mmm个采样点{ui∈R2}i=1m\{\bm {u}_i\in \mathbb{R}^2\}_{i=1}^m{ui∈R2}i=1m，通过对这些采样点的建模同样可以得到一个GMMs：
g(x∣θ)=∑i=1mωiN(x;μi,Σi),g(\bm x|\bm \theta)=\sum_{i=1}^m\omega_iN(\bm x;\bm \mu_i, \bm \Sigma_i), g(x∣θ)=i=1∑mωiN(x;μi,Σi),
其中μi\bm \mu_iμi是两个采样点之间的中点位置μi=ui+ui+1x\bm \mu_i=\frac{\bm {u}_i + \bm {u}_{i+1}}{x}μi=xui+ui+1, Σi=QiTΛQi\bm \Sigma_i=\bm {Q}^T_i \Lambda\bm Q_iΣi=QiTΛQi是协方差矩阵, 其中Qi=[n1i∣n2i]\bm Q_i=[\bm n_{1i}|\bm {n}_{2i}]Qi=[n1i∣n2i]由在点ui\bm u_iui处两个正交的单位切向量和法向量构成。ωi=hi∑i=1mhi\omega_i=\frac{h_i}{\sum_{i=1}^mh_i}ωi=∑i=1mhihi是GMMs每一个高斯核的权重系数。作者通过找到最优化下式来得到最终的椭圆参数：
θ^=arg⁡min⁡θ∥f(x)−g(x∣θ)∥2,\hat{\bm \theta}=\arg\min_{\bm \theta}\| f(\bm x)-g(\bm x|\bm \theta)\|_2, θ^=argθmin∥f(x)−g(x∣θ)∥2,

利用两个高斯分布之间内积公式上式可以通过梯度下降的方法得到解决（文中为了提高精度加入了模拟退火算法迭代的进行梯度下降，设置hih_ihi在[hmin,hmax][h_{min},h_{max}][hmin,hmax]以0.8位step进行迭代）。
⟨N(μ1,Σ1)∣N(μ2,Σ2)⟩=N(0;μ1−μ2,Σ1+Σ2).\left \langle N(\bm \mu_1,\bm \Sigma_1) | N(\bm \mu_2,\bm \Sigma_2) \right \rangle = N(\bm 0;\bm \mu_1-\bm \mu_2,\bm \Sigma_1+\bm \Sigma_2). ⟨N(μ1,Σ1)∣N(μ2,Σ2)⟩=N(0;μ1−μ2,Σ1+Σ2).

Ellipse fitting based on HGMMs

Distance-Density-Based Error Metric

作者通过实验发现最优的L2L_2L2 error并不能完全代表最优的椭圆参数，因此提出了一种同时结合(1) distance, and (2) density的error metric方法,也就是说需要 (1) 靠近模型的points（inliers）越多越好，（2）inliers对椭圆的损失要越小越好，可以表示为：
E(θ)=∑j=1ndθ(xj)exp⁡{∑jn1(dj2<T2)f(xj)},E(\bm \theta)=\frac{\sum_{j=1}^nd_{\bm \theta}(\bm x_j)}{\exp\{\sum_j^n \mathbf{1}(d_j^2<T^2)f(\bm x_j)\}}, E(θ)=exp{∑jn1(dj2<T2)f(xj)}∑j=1ndθ(xj),
其中dθ(xj)d_{\bm \theta}(\bm x_j)dθ(xj)是一个距离函数，1:R→{0,1}\mathbf{1}:\mathbb{R}\to\{0,1\}1:R→{0,1}是一个映射函数，djd_jdj是一点到椭圆的距离（这一块论文里说的比较清楚）。

Density-Based Region Growing Method

考虑到GMMs随意选取初始椭圆参数并且粗暴的直接将所有点进行建模和最优化，效率很低且容易受到outliers（离群点）的影响。因此，作者提出了一个包含两阶段的递进式椭圆拟合方法（内层和外层）。给定一组points C={xi∈R2}i=1nC=\{\bm x_i\in\mathbb{R}^2\}_{i=1}^nC={xi∈R2}i=1n, 外层HGMMs的遵循以下步骤：

计算每一个数据点的概率密度值ρ(xi)=f(xi)\rho(\bm x_i)=f(\bm x_i)ρ(xi)=f(xi), 然后选定最大密度值点ρ0=max⁡i=1nρ(xi)\rho_0=\max_{i=1}^n\rho(\bm x_i)ρ0=maxi=1nρ(xi)作为种子点开始区域增长，如果一个点在种子点半径为ϵ\epsilonϵ圆内并且改点的密度值也就是ρ(xi)\rho(\bm x_i)ρ(xi)不能比 ρ0\rho_0ρ0 小太多 (ρ0\rho_0ρ0 已经是最大密度点)，那么该点就被加入到增长区域 RRR 中，也就是满足下式：
NEps(x0)={xi∈C∣∥xi−x0∥≤ϵ∩ρ(xi)ρ0≥τ}.{\rm NEps}(\bm x_0)=\{\bm x_i\in C|\|\bm x_i - \bm x_0\|\le\epsilon\cap \frac{\rho(\bm x_i)}{\rho_0}\ge\tau \}. NEps(x0)={xi∈C∣∥xi−x0∥≤ϵ∩ρ0ρ(xi)≥τ}.
当 RRR 无法加入新的 xi\bm x_ixi 时，选择 RRR 中离 x0\bm x_0x0 距离最远的一个点作为新的种子点开始新的区域增长（更新同一个 RRR），当 ∣R∣|R|∣R∣ （也就是 RRR 中点的个数）超过给定阈值 NNN 时，区域增长结束输出点集 RRR (如果无法更新 RRR时也没达到 NNN 就将这一组点认为是 outliers 暂时从 CCC 中删除)。
根据步骤1中输出的 inliers 集合 RRR，利用GMMs算法拟合出一个椭圆θ0\bm \theta_0θ0.

内层HGMMs以外层输出 θ0\bm \theta_0θ0 作为初始椭圆，以 CCC 中所有数据点作为inputs 点集，再次进行GMMs拟合出一个最终的椭圆。

一点点理解

GMMs中的参数 mmm (也就是下图的 ‘NNN’) 是控制椭圆模型 g(x∣θ)g(\bm x|\bm \theta)g(x∣θ) 中采样点 ui\bm u_iui 的个数，也就是GMMs中高斯核的个数，过小的高斯核过于孤立不能很好的表达一个椭圆的轮廓（如图e），过多的高斯核又有太多的冗余信息影响算法效率，论文中经过实验选择为了20。hhh 是控制是控制单个高斯核的“胖瘦”（协方差矩阵），胖一点的高斯核可以接收到更多input points的反馈但也容易受离群点的影响，瘦一点的高斯核更精确但容易过拟合，论文中用模拟退火算法来得到最佳 hhh。

HGMMs主要贡献是两个：（1）在GMMs外部套了一层寻找最优初始椭圆 θ0\bm \theta_0θ0 的预处理步骤，从论文的实验来看这一步的效果很好，最优化问题的初始值选取一直是一个十分重要的问题（CNN网络的权值初始化也是这个道理）；（2）优化目标不仅仅是考虑 density 的 L2L_2L2 loss 而是同时考虑了distance 和 density，论文中也用实验证明了最优的 L2L_2L2 loss 并不能映射为最优的椭圆参数 θ\bm \thetaθ.