图像对齐：Parametric Image Alignment Using Enhanced Correlation Coefficient Maximization

Parametric Image Alignment Using Enhanced Correlation Coefficient Maximization

对齐的目的是找到两张图像像素点坐标间的对应关系，对于手持设备，首先假设两张图像间为简单的仿射变换，可以通过求解两张图像间的仿射变换矩阵，实现对齐图像的目的。作者导出ECC这个算法的两个核心思想，一个来源与作者本人的另一篇论文《An Enhanced Correlation-Based Method for Stereo Correspondence with Sub-Pixel Accuracy》中对增强相关系数的与残差之间的等同关系，一个来源与Takeo Kanade在其论文《An Iterative Image Registration Technique with an Application to Stereo Vision》中提出的Forward Additive迭代思想。

Problem Formulation

假设有一对图像，图像间仅存在像素点坐标间的畸变，记其中一帧图像为参考帧$I_r$，对应的像素坐标记为$\mathbf{x}=[x_r,y_r{]}^t$，另一帧图像为畸变帧$I_w$，对应的像素坐标记为$\mathbf{y}=[x_w,y_w{]}^t$，对齐的目的是找到$\mathbf{x}$与$\mathbf{y}$之间的对应关系$?$，即$\mathbf{y}=?(\mathbf{x};p)$，其中$\mathbf{p}=[p_1,......,p_N{]}^T$为$N$维的参数向量，因此可以将对齐问题看作是一个参数估计问题：
$$I_r(\mathbf{x})=I_w(?(\mathbf{x};\mathbf{p}))\text{\hspace{1em}(1)}$$
将参数$\mathbf{p}$的优化问题，定义为参考帧与畸变帧经过畸变矫正的结果之间的差值的最小化问题：
$$\underset{\mathbf{p},\mathbf{\alpha }}{\min }E(\mathbf{p},\mathbf{\alpha })=\underset{\mathbf{p},\mathbf{\alpha }}{\min }\sum _{\mathbf{x}\in \mathcal{T}}{\left|I_r(\mathbf{x})?\Psi \left(I_w(?(\mathbf{x};\mathbf{p})),\mathbf{\alpha }\right)\right|}^p\text{\hspace{1em}(2)}$$
其中$\Psi (I,\alpha )$为简单的亮度转换，目的是降低由图像亮度变化带来的图像间的不相关性。

Proposed Criterion and Main Results

Vectorization

记参考图像与目标图像间的目标对齐区域的中心坐标数为$K$个，对齐区域的亮度值分别构成参考向量$i_w(p)$、畸变向量$i_r$，分别为：
$$\begin{array}{rl}{\mathbf{i}}_r& ={\left[I_r\left({\mathbf{x}}_1\right)I_r\left({\mathbf{x}}_2\right)?I_r\left({\mathbf{x}}_K\right)\right]}^t\\ {\mathbf{i}}_w(\mathbf{p})& ={\left[I_w\left({\mathbf{y}}_1(\mathbf{p})\right)I_w\left({\mathbf{y}}_2(\mathbf{p})\right)?I_w\left({\mathbf{y}}_K(\mathbf{p})\right)\right]}^t\end{array}\text{\hspace{1em}(3)}$$
${\bar{i}}_r$、${\bar{i}}_w(p)$分别为对应的零均值向量(column-zero–mean versions)，用ECC作为参数$p$的性能评价指标：
$$E_{ECC}(p)={\Vert \frac{{\bar{i}}_r}{\Vert {\bar{i}}_r\Vert }?\frac{{\bar{i}}_w(p)}{\Vert {\bar{i}}_w(p)\Vert }\Vert }^2\text{\hspace{1em}(4)}$$
可以看出增益及亮度整体抬升不影响评价指标

Performance Measure Optimization

在本文作者的另外一篇论文《An Enhanced Correlation-Based Method for Stereo Correspondence with Sub-Pixel Accuracy》中证明了对公式(4)的最小化等同于对下式增强相关系数(Enhanced Correction Coefficient)的最大化：
$$\rho (p)=\frac{{\bar{i}}_r^T{\bar{i}}_w(p)}{\Vert {\bar{i}}_r\Vert \Vert {\bar{i}}_w(p)\Vert }={\hat{i}}_r^T\frac{{\bar{i}}_w(p)}{\Vert {\bar{i}}_w(p)\Vert }\text{\hspace{1em}(5)}$$
其中${\hat{i}}_r=\frac{{\bar{i}}_r}{\Vert {\bar{i}}_r\Vert }$，由于不对参考帧做畸变矫正，此向量恒定；由于$\Vert {\bar{i}}_w(p)\Vert$的存在，即使${\bar{i}}_w(p)$相对于参数$p$线性，目标函数$\rho (p)$相当于参数$p$依旧非线性，将参数向量$p$做如下表示：
$$p=\tilde{p}+△p$$
$\tilde{p}：$nominal parameter vector；$△p：$perturbations vector
这里的三个参数向量即是表示对参数向量$p$的近似过程，也表示后续参数向量的迭代求解方式
对图像在x-y坐标系下进行一个简单的一阶泰勒展开，可以将畸变帧做如下近似：
$$I_w(\mathbf{y})\approx I_w(\tilde{\mathbf{y}})+\left[{?}_{\mathbf{y}}{I}_w(\tilde{\mathbf{y}})\right]\frac{??(\mathbf{x};\tilde{\mathbf{p}})}{?\mathbf{p}}\Delta \mathbf{p}\text{\hspace{1em}(6)}$$
其中$[{▽}_Y{I}_w(\tilde{Y})]$为一个$2\times K$的梯度矩阵；$\frac{??(X;\tilde{p})}{?p}$为一个$2\times N$的Jacobian矩阵。同样，将以上线性近似应用于全部目标区域，可得：
$$i_w(p)\approx i_w(\tilde{p})+G(\tilde{p})△p\text{\hspace{1em}(7)}$$
其中$G(\tilde{p})$为$K\times N$的Jacobian矩阵(Jacobian matrix of the warped intensity vector with respect to the parameters)
假设畸变模型为：
$$?(X;p)=\left[\begin{array}{cc}{?}_1(X;p)& {?}_2(X;p{)}^T\end{array}\right]\text{\hspace{1em}(8)}$$
可用以下公式确定$K\times N$的Jacobian矩阵：
$$G(\tilde{p}{)}_{k,n}=\sum _{i=1}^2(\frac{?I_w(Y)}{?Y_i}{|}_{Y=Y_k(\tilde{p})}\times \frac{?{?}_i(X_k;p)}{?p_n}{|}_{p=\tilde{p}})\text{\hspace{1em}(9)}$$
其中$Y$为像素坐标，包含$Y_i、Y_2$，将式(7)带入式(5)，得到新的增强相关系数的表达式：
$$\rho (p)\approx \rho (△p|\tilde{p})={\hat{i}}_r^T\frac{{\bar{i}}_w(\tilde{p})+\bar{G}(\tilde{p})△p}{\Vert {\bar{i}}_w(\tilde{p})+\bar{G}(\tilde{p})△p\Vert }\text{\hspace{1em}(10)}$$
其中$\bar{G}(\tilde{p})$、${\bar{i}}_w(\tilde{p})$分别对应$G(\tilde{p})$、$i_w(\tilde{p})$的零均值向量。为了表达的简洁，之后的推导将省略参数向量$\tilde{p}$，将式(10)直接做如下表达：
$$\rho (\Delta p\mid \tilde{p})=\frac{{\hat{i}}_r^T{\overline{i}}_w+{\hat{i}}_r^T\bar{G}\Delta p}{\sqrt{{\Vert {\overline{i}}_w\Vert }^2+2{\overline{i}}_w^T\bar{G}\Delta p+\Delta p^T{\bar{G}}^T\bar{G}\Delta p}}\text{\hspace{1em}(11)}$$

Closed-Form Expression:

考虑如下标量函数：
$$f(\mathbf{x})=\frac{u+{\mathbf{u}}^t\mathbf{x}}{\sqrt{v+2{\mathbf{v}}^t\mathbf{x}+{\mathbf{x}}^{t}Q\mathbf{x}}}\text{\hspace{1em}(12)}$$
其中$u,v$为标量；$\mathbf{u},\mathbf{v}$是长度为$N$的向量；$Q$是正定阵，如果满足$v>{\mathbf{v}}^t{Q}^{?1}\mathbf{v}$，则可以分情况确定函数$f(x)$是极大值。
1.当$u>{\mathbf{u}}^t{Q}^{?1}\mathbf{v}$时，$f(x)$的最大值如下:
$$\underset{\mathbf{x}}{\max }f(\mathbf{x})=\sqrt{\frac{{\left(u?{\mathbf{u}}^t{Q}^{?1}\mathbf{v}\right)}^2}{v?{\mathbf{v}}^t{Q}^{?1}\mathbf{v}}+{\mathbf{u}}^t{Q}^{?1}\mathbf{u}};\mathbf{x}=Q^{?1}\left\{\frac{v?{\mathbf{v}}^t{Q}^{?1}\mathbf{v}}{u?{\mathbf{u}}^t{Q}^{?1}\mathbf{v}}\mathbf{u}?\mathbf{v}\right\}\text{\hspace{1em}(14、15)}$$
2.当$u\le {\mathbf{u}}^t{Q}^{?1}\mathbf{v}$时，$f(x)$的上确界如下：
sup?xf(x)=utQ?1u;x=Q?1{λu?v};λ>0(16、17)\sup _{\mathbf{x}} f(\mathbf{x})=\sqrt{\mathbf{u}^{t} Q^{-1} \mathbf{u}};\mathbf{x}=Q^{-1}\{\lambda \mathbf{u}-\mathbf{v}\};\lambda>0 \tag{16、17} xsup?f(x)=utQ?1u ?;x=Q?1{ λu?v};λ>0(16、17)

Derive Forward Additive ECC

首先给出正交投影矩阵$P_G=\bar{G}{\left({\bar{G}}^t\bar{G}\right)}^{?1}{\bar{G}}^t$，因此，以下不等式成立：
$${\Vert {\overline{\mathbf{i}}}_w\Vert }^2={\Vert P_G{\overline{\mathbf{i}}}_w\Vert }^2+{\Vert \left[I?P_G\right]{\overline{\mathbf{i}}}_w\Vert }^2\ge {\Vert P_G{\overline{\mathbf{i}}}_w\Vert }^2={\overline{\mathbf{i}}}_w^t{P}_G{\overline{\mathbf{i}}}_w\text{\hspace{1em}(18)}$$
其中$I$为单位阵，当且仅当$\left[I?P_G\right]{\overline{\mathbf{i}}}_w=0$时等式成立，而$\left[I?P_G\right]{\overline{\mathbf{i}}}_w?=0$，所以不等式(18)严格成立。又有式(18)形如式(13)，式(11)形如式(12)，即将式(15)带入式(11)可解得$△p$：
$$\Delta \mathbf{p}={\left({\bar{G}}^t\bar{G}\right)}^{?1}{\bar{G}}^t\left\{\frac{{\Vert {\overline{\mathbf{i}}}_w\Vert }^2?{\overline{\mathbf{i}}}_w^t{P}_G{\overline{\mathbf{i}}}_w}{{\hat{\mathbf{i}}}_r^t{\overline{\mathbf{i}}}_w?{\hat{\mathbf{i}}}_r^t{P}_G{\overline{\mathbf{i}}}_w}?{\overline{\mathbf{i}}}_w\right\}\text{\hspace{1em}(19)}$$
当${\hat{\mathbf{i}}}_r^t{\overline{\mathbf{i}}}_w>{\hat{\mathbf{i}}}_r^t{P}_G{\overline{\mathbf{i}}}_w$时，根据式(17)，$\lambda ＞0$，既有：
$$\Delta \mathbf{p}={\left({\bar{G}}^t\bar{G}\right)}^{?1}{\bar{G}}^t\left\{\lambda {\hat{\mathbf{i}}}_r?{\overline{\mathbf{i}}}_w\right\}\text{\hspace{1em}(20)}$$
当${\hat{\mathbf{i}}}_r^t{\overline{\mathbf{i}}}_w\le {\hat{\mathbf{i}}}_r^t{P}_G{\overline{\mathbf{i}}}_w$时，为保证$\rho (\Delta \mathbf{p}\mid \tilde{\mathbf{p}})>\rho (0\mid \tilde{\mathbf{p}})$根据文中Lemma 1，确定$\lambda$的取值：
λ1=iwtPGi?wi^rtPGi^r,λ2=i^rtPGi?w?i^rti?wi^rtPGi^r,λ≥max?{λ1,λ2}(21)\lambda_{1}=\sqrt{\frac{\mathbf{i}_{w}^{t} P_{G} \overline{\mathbf{i}}_{w}}{\hat{\mathbf{i}}_{r}^{t} P_{G} \hat{\mathbf{i}}_{r}}}, \quad \lambda_{2}=\frac{\hat{\mathbf{i}}_{r}^{t} P_{G} \overline{\mathbf{i}}_{w}-\hat{\mathbf{i}}_{r}^{t} \overline{\mathbf{i}}_{w}}{\hat{\mathbf{i}}_{r}^{t} P_{G} \hat{\mathbf{i}}_{r}},\quad \lambda \geq \max \left\{\lambda_{1}, \lambda_{2}\right\}\tag{21} λ1?=i^rt?PG?i^r?iwt?PG?iw?? ?,λ2?=i^rt?PG?i^r?i^rt?PG?iw??i^rt?iw??,λ≥max{ λ1?,λ2?}(21)
将式(20)带入是(11)，得到关于$\lambda$的表达式：
$$\rho (\Delta p\mid \tilde{p})=f(\lambda )=\frac{\left({\hat{\mathbf{i}}}_r^t{\overline{\mathbf{i}}}_w?{\hat{\mathbf{i}}}_r^t{P}_G{\overline{\mathbf{i}}}_w\right)+\lambda {\hat{\mathbf{i}}}_r^t{P}_G{\hat{\mathbf{i}}}_r}{\sqrt{\left({\Vert {\overline{\mathbf{i}}}_w\Vert }^2?{\overline{\mathbf{i}}}_w^t{P}_G{\overline{\mathbf{i}}}_w\right)+{\lambda }^2{\hat{\mathbf{i}}}_r^t{P}_G{\hat{\mathbf{i}}}_r}}\text{\hspace{1em}(22)}$$
简单计算式(22)关于$\lambda$的导数，可知$f(\lambda )$单调增
1.当$\lambda \ge {\lambda }_2$时，$f(\lambda)\geq$0

2.当$\lambda ={\lambda }_1$时，$f\left({\lambda }_1\right)=\frac{{\hat{\mathbf{i}}}_r^t{\overline{\mathbf{i}}}_w?{\hat{\mathbf{i}}}_r^t{P}_G{\overline{\mathbf{i}}}_w+\sqrt{\left({\overline{\mathbf{i}}}_w^t{P}_G{\overline{\mathbf{i}}}_w\right)({\hat{\mathbf{i}}}_r^t{P}_G{\hat{\mathbf{i}}}_r})}{\Vert {\overline{\mathbf{i}}}_w\Vert }\ge \rho (0\mid \tilde{\mathbf{p}})$

即当λ≥max?{λ1,λ2}\lambda \geq \max \left\{\lambda_{1}, \lambda_{2}\right\}λ≥max{ λ1?,λ2?}时$f(\lambda )＞0$，$\Delta p$总会带来相关系数的增大，即是$\rho (\tilde{p}+\Delta p)＞\rho (\tilde{p})$，可以通过迭代求解$\Delta P$实现对参数$p$的优化

Forward Additive ECC Iterative Algorithm

S0.用参考帧$I_r$计算${\hat{i}}_r$；初始化参数向量$p_0$；给定迭代次数$T$；给定$\varepsilon $；$j=1$

S1.用$?(X;p_{j?1})$对畸变帧$I_w$做畸变矫正，并用矫正后的结果计算${\hat{i}}_w(p_{j?1})$

S2.用公式(9)计算Jacobin矩阵

S3.根据式(19、20、21)计算$\Delta p_j$

S4.更新参数向量：$p_j=p_{j?1}+\Delta p_j，j=j+1$

S5.当$\Vert \Delta p_j\Vert \ge \epsilon$ 时，停止迭代

Inverse Compositional ECC Iterative Algorithm

FA_ECC需要计算到$G^{t}G$，算法复杂度为$O(K{N}^2)$，除此之外的算法复杂度为$O(KN)$。可以通过变换$I_r，I_{}w$的位置来将$G^{t}G$的计算次数降为1，用参考帧$I_r$计算Jacobin矩阵，此矩阵恒定，不需要再次计算