试图理解这篇文章，写个笔记。

仅关注方法部分。

相关基础背景

忘光了，记一下

1. 通俗的解释最大似然估计
   也就是说，贝叶斯统计推断和最大似然估计是根据观测值去反推参数的两个不同的路子。贝叶斯统计：未知参数是服从已知的分布的随机变量；最大似然估计：参数虽未知，但它是常数，看看它最大的可能是哪个数。
   前提假设：(1) 样本之间相互独立；(2) 分布确定，即质量密度函数确定

2. EM算法
   最大似然估计和EM算法都是根据实现结果求解概率分布的最佳参数 $\theta$ ，但最大似然估计中知道每个结果对应哪个概率分布（我知道哪个概率分布实现了这个结果），而EM算法面临的问题是：我不知道哪个概率分布实现了该结果。

3. 高斯径向基函数
   函数的形式（$x_c$ 是中心点）：
   $$k\left(x-x_c\right)=\exp \left(-\frac{{\Vert x-x_c\Vert }^2}{2{\sigma }^2}\right)$$
   双竖运算符 $\Vert \cdot \Vert$ 表示Norm运算，即取向量的“度量”。例如，二维下常为距离函数。 二维函数图像：
   

4. 高斯混合模型（Gaussian Mixture Models，GMM）
   GMM可用于无监督学习中的聚类的数据，其方式与$k$-means几乎相同。但是，与$k$-means相比，使用GMM有两个优点：1） $k$-means不考虑方差，GMM考虑方差，并返回数据点属于K个群集中每个群集的概率；2）$k$-means执行硬分类，GMM执行软分类

---

方法框架

该算法的目标是将遥感图像 $I_y$ 进行变换，使其与参考图像 $I_x$ 配准。我们从参考图像中检测一个特征点集 $X$ ，从感知图像中检测一个特征点集 $Y$ 。接下来我们使用基于期望最大化(EM)的过程来获得变换后的 $Y$ 的位置，即 $Z$ 。然后使用 $Y$ 和 $Z$ 来求解用于图像变换的薄板样条函数(thin plate spline，TPS)插值¹。方法的主要过程显示在论文的图一中。

也就是说这个事情的前提是啥呢，是已经有一个训练好的CNN，图像输进去，得到特征图。
$$I_x\to \left(CNN\right)\to X；I_y\to \left(CNN\right)\to Y$$问题来了：这个CNN针对什么训练？分类？目标检测？这里先不管。
貌似在这篇文章里面有写方法：《Discriminative Unsupervised Feature Learning with Exemplar Convolutional Neural Networks》
相关解析：图像配准的前世今生：从人工设计特征到深度学习

相关符号：

• $X_{N\times 2}$,$Y_{M\times 2}$ ——分别从参考图像和感知图像中提取特征点集。
• $Z$ —— $Y$ 的变换位置
• $N$ , $M$ —— 分别代表 $X$ 和 $Y$ 中的点数
• $x_n$ , $y_m$ ——点集 $X$ 中第 $n$ 个元素和点集 $Y$ 中第 $m$ 个元素。

按照后面的描述，$Z=Y+GW$，$G$ 是高斯径向基函数。

特征描述和预匹配

学习特征

使用的网络：VGG-16
输入图片的尺寸：224 $\times$ 224

生成特征描述符

用来生成特征描述符的特征图：pool3，pool4和pool5_1。这些层搜索一组通用的模式，并产生能很好地覆盖不同大小的感受域的特征响应值。

背景知识：Kronecker积

把输入图片分成28 $\times$ 28的格子，每个格子对应一个尺寸为 8 $\times$ 8的切片单元。

pool3出来的特征图是256个通道，也就是说，用一个256维数据来表征一个 8 $\times$ 8的切片的特征。
pool3出来的特征图 $O_{pool3}$ 直接使用，是个28 $\times$ 28 $\times$ 256 的张量，记为 $F_1$，即 $F_1=O_{pool3}$；

pool4出来的特征图 $O_{pool4}$ 是个14 $\times$ 14 $\times$ 512 的张量，为了保持和 $F_1$ 的尺寸一致，做个变形，就是把长和宽都扩大一倍，原来一个点，扩大为4个点，变为 $F_2$， $$F_2=O_{pool4}\otimes I_{2\times 2\times 1}$$

$I$ 表示下标形状的张量，元素全部是1。 也就是说，对于 $F_2$ 每个通道，实际上是对 $O_{pool4}$ 进行了直接的放大，放大2倍，（还没有带消除锯齿的操作）。

pool5_1出来的特征图 $O_{pool5\_1}$ 是个7 $\times$ 7 $\times$ 512 的张量，同上面操作，变为 $F_3$， $$F_3=O_{pool5\_1}\otimes I_{4\times 4\times 1}$$

对 $F_1$，$F_2$ 和 $F_3$ 进行归一化

$$F_i\leftarrow \frac{F_i}{\sigma \left(F_i\right)},i=1,2,3$$
分别将点 $x$ 的pool3、pool4和pool5_1的描述符记为 $D_1\left(x\right)$, $D_2\left(x\right)$ 和 $D_3\left(x\right)$

特征的预匹配

$x$, $y$ 两点之间的特征距离
$$d(x,y)=\sqrt{2}{d}_1\left(x,y\right)+d_2\left(x,y\right)+d_3\left(x,y\right)$$
为什么要加这个 $\sqrt{2}$ 呢，是因为$D_1$是256通道的，$D_2$、$D_3$都是512通道的。

$d_i\left(x,y\right)是$ $x$ 和 $y$ 两点的欧拉距离。

$x$ 和 $y$ 二者匹配的条件：

1. $d\left(\cdot ,y\right)$之中，$d\left(x,y\right)$最小。
2. 不存在使 $d\left(z,y\right)<\theta \cdot d\left(x,y\right)$ 的 $d\left(z,y\right)$ 值。

$\theta$ 取大于1的值，被称为匹配阈值。

这种匹配方法不能保证双射（bijection）

个人理解：遥感的图片有的点参考图片里没有，参考图片里有的点遥感图片里没有。

动态的内层（inlier ）选择

inlier = 内层？没相关解释。

我们的特征点生成在正方形图像块的中心。$Y$ 变形为 $Z$ 的情况下，$Z(i)$ 可能和 $X(i)$ 对应特征点的图像斑块可能部分或完全重叠。因此，为了实现更准确的配准，重叠比例较大的特征点应具有更好的对准程度，而部分重叠的图块的中心距离应较小。对齐程度是使用我们的动态内层选择来确定的。

在使用EM算法迭代求解 $Z$ ($Y$ 在每一次迭代中变换的位置)的同时，我们在每 $k$ 次迭代中更新内层的选择。选取的点作为内线引导点位置的运动，而离群点的运动则是一致的。在特征预匹配阶段，利用低阈值 ${\theta }_0$ 选择大量的特征点来过滤掉不相关的点。然后我们指定一个大的起始阈值 $\hat{\theta }$ ，只有可信的内层点(有重叠斑块的特征点)满足。在其余的配准过程中，阈值 $\theta$ 在每 $k$ 次迭代中被步长 $\delta$ 减去，允许更多的特征点影响变换。这种做法可以使强匹配的特征点决定整体的变换，而其他特征点优化配准精度。

内层选择过程产生一个维度为 $M\times N$ 的先验概率矩阵 $P_R$，然后由我们的基于高斯混合模型(Gaussian Mixed Model，GMM)的变换求解器得到。矩阵元素$P_R\left[m,n\right]$是$x_n$和$y_m$对应的假定概率（putative probability）。假设$x_n$对应于$y_m$，则得到$P_R\left[m,n\right]$较大的假定概率。而且大的概率会进一步导致$y_m$上的一个明显的转换，通过这个转换对应的对可以对齐。

假定概率使用卷积特征和几何结构信息来确定。先验概率矩阵$P_R$通过以下步骤得到：

1）准备 $M\times N$的卷积特征代价矩阵 $C_{\theta }^{conv}$：

$$C_{\theta }^{conv}\left[m,n\right]=\{\begin{array}{rl}\frac{d\left(y_m,x_n\right)}{d_{\theta }^{max}},& cond1\\ 1,& otherwise\end{array}$$

Cond 1的定义：$x_n$ 和 $y_m$ 在阈值 $\theta$ 的条件下是匹配的。 $d\left(\cdot ,\cdot \right)$是之前定义的特征距离。$d_{\theta }^{max}$ 是在阈值 $\theta$ 下所有匹配的特征点对之中最大的距离。

2）用论文《Shape matching and object recognition using shape contexts》提出的方法来计算几何结构的代价矩阵 $C^{geo}$，它是一种基于直方图的描述符，描述点的邻域结构。描述符将轮廓点置于极坐标系统的中心，并记录落在弧形箱中的点的数量。$C^{geo}$ 通过执行卡方检验 ${\mathcal{X}}^2$ 来计算

$$C^{geo}\left[m,n\right]=\frac{1}{2}\sum _{b=1}^B\frac{{\left[h_m^y\left(b\right)-h_n^x\left(b\right)\right]}^2}{h_m^y\left(b\right)+h_n^x\left(b\right)}$$
$h_m^y\left(b\right)$ 和 $h_n^x\left(b\right)$ 分别表示落在 $y_m$ 和 $x_n$ 附近的第 $b$ 个单元柱中的数量。

3） $C^{geo}$ 和 $C_{\theta }^{conv}$ 的取值范围都是 $\left[0,1\right]$，我们使用Hadamard积来计算整合的代价矩阵 $C$

$$C=C_{\theta }^{conv}\odot C^{geo}$$

4） 应用Jonker-Volgenant算法（《A shortest augmenting path algorithm for
dense and sparse linear assignment problems》）求解了代价矩阵 $C$ 的线性分配问题。指定的点对被认为是假定对应的。最后，我们用下面的公式计算先验概率矩阵

$$P_R\left[m,n\right]=\{\begin{array}{rl}1,& \text{ if }{y}_m\text{ and }{x}_n\text{ are corresponding}\\ \frac{1-ϵ}{N},& otherwise\end{array}$$

$ϵ$ 是在$\left[0,1\right]$之间取值的超参数，它应该根据我们对内层选择的置信值来指定，以确保准确。先验概率需要进行归一化：

$$P_R\left[m,n\right]:=\frac{P_R\left[m,n\right]}{{\sum }_{k=1}^N{P}_R\left[m,k\right]}$$

我们考虑点集 $Y$ 作为高斯混合模型（GMM）质心。GMM概率密度函数定义为：

$$p\left(x\right)=\omega \frac{1}{N}+\left(1-\omega \right)\sum _{m=1}^M{g}_m\left(x\right)$$

$g_m\left(x\right)$ 是个正态分布密度函数

$$g_m\left(x\right)=\frac{1}{2\pi {\sigma }^2}\exp \left(-\frac{1}{2{\sigma }^2}{\Vert x-y_m\Vert }^2\right)$$

该模型对混合变量中的每个高斯质心使用各向同性方差 ${\sigma }^2$。一个额外的均匀分布项 $\frac{1}{N}$被加入到一个权重参数 $\omega$ 中来考虑异常值（outliers），$0<\omega <1$ 。

然后，我们使用期望最大化(EM)算法寻找最优的转换参数$\left(W,{\sigma }^2,\omega \right)$。这种方法的目标是最大化似然函数，或等效地最小化负对数似然函数。

$$L\left(W,{\sigma }^2,\omega \right)=-\sum _{n=1}^N\log \sum _{m=1}^{M+1}{P}_R\left[m,n\right]g_m\left(x_n\right)$$

从这个公式中，由于存在不能被直接观测的变量 $m$，不能直接计算梯度。取而代之，EM算法使负对数似然函数的期望最小：

$$Q=-\sum _{n=1}^N\sum _{m=1}^{M+1}{P}^{old}\left(m|x_n\right)\log \left(P_R\left[m,n\right]g_m\left(x_n\right)\right)$$

$P^{old}\left(m|x_n\right)$ 表示用上次迭代的参数计算的后验概率项。将方程展开，省略导数冗余项后，方程可改写为：

$$\begin{array}{rl}Q\left(W,{\sigma }^2,\omega \right)=& \frac{1}{2{\sigma }^2}\sum _{n=1}^N\sum _{m=1}^M{P}^{old}\left(m|x_n\right){\Vert x_n-\tau \left(y_m,W\right)\Vert }^2\\ -& \frac{1}{2}{N}_P\log \left(\frac{{\sigma }^2\omega }{1-\omega }\right)-N\log \left(\omega \right)\end{array}$$

$N_P={\sum }_{n=1}^N{\sum }_{m=1}^M{P}^{old}\left(m|x_n\right)$，$\tau \left(y_m,W\right)$ 表示变换后 $y_m$ 的位置。

非刚性变换定义为：
$$Z=Y+GW$$

式子中 $G$ 是高斯径向基函数（Gaussian radial basis function ，GRBF），$W$ 包括要学习的变换参数。

$$G\left[i,j\right]=\exp \left(-\frac{{\Vert x_j-y_i\Vert }^2}{2{\beta }^2}\right)$$

加入一项基于运动一致性理论的正规项（见《A mathematical analysis of the motion coherence theory》）（Motion Coherence Theory，MCT），得到

$$Q_r=Q+\frac{\lambda }{2}tr\left(W^{T}GW\right)$$

$tr(\cdot )$ 表示迹操作（trace operation）。迹等于对角线之和等于特征值之和。

EM算法迭代计算期望和最小梯度，直到收敛。

E-Step：从最后迭代的得到的参数中计算后验概率矩阵 $P_O$

$$P_O\left[m,n\right]=P^{old}\left(m|x_n\right)=\frac{P_R\left[m,n\right]g_m\left(x_n\right)}{p\left(x_n\right)}$$

M-Step：求导和更新参数。

$$W:={\left(G+\lambda {\sigma }^2{P}_d^{-1}\right)}^{-1}\cdot \left(P_d^{-1}PX-Y\right)$$

$${\sigma }^2:=\frac{1}{2N_P}\left(tr\left(X^T{P}_{d}X\right)-2tr\left(X^T{P}^{T}Z\right)+tr\left(Z^T{P}_{d}Z\right)\right)$$

$$\omega :=1-\frac{N_P}{N}$$

$P_d=diag\left(P1\right)$， $1$ 表示元素都是1的列向量。

---

1. 样条函数其实就是分段函数，表示在相邻点 $x_k$ 和 $x_{k+1}$ 之间用低阶多项式函数 $S_k(x)$ 进行插值。分段线性插值和三次样条插值都属于样条插值。 ↩︎