图像对齐（ECC，Enhanced Correlation Coefficient）

该方法收录在OpenCV中：cv2.findTransformECC
论文名为：《Parametric Image Alignment Using Enhanced Correlation Coefficient Maximization》（2008年）
教程在此：https://github.com/spmallick/learnopencv/tree/master/ImageAlignment
更多教程在此：https://github.com/spmallick/learnopencv

---

2 问题定义

假设我们有一对图像轮廓（强度）： $I_r(\text{x})$ 和 $I_w(\text{y})$。
第一个是 reference（参考）图像，或者叫做 template（模板）图像。
第二个是 warped（扭曲的）图像。
$\text{x}=[x_1,x_2{]}^t，\text{y}=[y_1,y_2{]}^t$ 表示坐标。 （注意这里小写的 ${}^t$ 表示转置，全文都是）

再假设我们有一组 reference 图像上的坐标 T = { x k , k = 1 , … , K } \mathcal{T} = \{ \mathbf{x}_k, k=1, \dots,K \} T={xk​,k=1,…,K} ，这个叫做 target area（目标区域）。

alignment problem （对齐问题）就是在 warped 图像上找到对应的坐标点。 当然我们感兴趣的不是任意对应关系，而是那些结构化的、能够用 well-defined 的向量映射关系 $\mathbf{y}=\mathbf{\varphi }(\mathbf{x};\mathbf{p})$ 来建模的那些点。这里 $\mathbf{p}=[p_1,\cdots ,p_N{]}^t$ 是一个未知参数的向量。

在实践中经常会面临这种对齐问题，最常见的场景就是对一个序列的图像做运动估计（motion estimation）。在这种应用里，由于场景和相机之间存在相对运动，整个（目标）区域在时间上表现得不同。

假设给定了变换模型（并且亮度恒定的假设有效），对齐问题就简化为 —— 估计参数 $\mathbf{p}$ 使得
$$I_r(\mathbf{x})=I_w(\mathbf{\varphi }(\mathbf{x};\mathbf{p})),\forall \mathbf{x}\in \mathcal{T}\text{\hspace{1em}(1)}$$
为了能够获取唯一解，未知参数量 $N$ 必须小于目标坐标数 $K$。
当然，在实际中往往 $K$ 远大于 $N$，这表明公式（1）是一个 over-determined（超定）的（非线性）方程组。

大多数现有算法尝试通过最小化两组轮廓之间的 difference 或者 dissimilarity 来求解参数向量 $\mathbf{p}$。
Dissimilarity 通过目标函数 $E(\mathbf{p})$ 来表示，它涉及到两个图像的强度差异的 $l_p$ 范数。
而在实际中，由于视角方向不同，并且/或者光照条件不同，违反了亮度恒定的假设。
有必要引入额外的光度变换 $\Psi (I,\mathbf{\alpha })$ 来考虑光照的变化，它用一组未知参数 $\mathbf{\alpha }$ 来表示。

典型优化问题具有以下形式：
$$\underset{\mathbf{p},\mathbf{\alpha }}{\min }E(\mathbf{p},\mathbf{\alpha })=\underset{\mathbf{p},\mathbf{\alpha }}{\min }\sum _{\mathbf{x}\in \mathcal{T}}|I_r(\mathbf{x})-\Psi (I_w(\mathbf{\varphi }(\mathbf{x};\mathbf{p})),\mathbf{\alpha }){|}^p\text{\hspace{1em}(2)}$$
我们需要提一下，这种形式的优化问题通常是不适定的，为了获得可接受的解，通常需要施加额外的规则（平滑性）条件。

由于对应性（correspondence）的部分涉及到非线性，这个最优化问题的求解显然不简单。现有方案的计算复杂度和估计质量取决于特定的 $l_p$ 范数，以及“ warping(扭曲) 和 photometric distortion(光度失真)”所采用的模型。就范数等级 $p$ 而言，大多数方法用的是 $p=2$（欧式范数）。我们也是用这个的。

3 提出的标准与主要结果

在 warping transformation(扭曲变换) $\mathbf{\varphi }(\mathbf{x};\mathbf{p})$ 下，目标区域 $\mathcal{T}$ 里面的坐标 ${\mathbf{x}}_k，k=1,\dots ,K$ 被映射到坐标 ${\mathbf{y}}_k(\mathbf{p})=\mathbf{\varphi }({\mathbf{x}}_k;\mathbf{p})，k=1,\dots ,K$ 上。

我们来定义 reference vector（参考向量）${\mathbf{i}}_r$ 以及对应的 warped vector（扭曲向量）${\mathbf{i}}_w(\mathbf{p})$ ：

$$\begin{array}{rl}{\mathbf{i}}_r& =[I_r({\mathbf{x}}_1)I_r({\mathbf{x}}_2)\cdots I_r({\mathbf{x}}_K){]}^t\\ {\mathbf{i}}_w(\mathbf{p})& =[I_w({\mathbf{y}}_1(\mathbf{p}))I_w({\mathbf{y}}_2(\mathbf{p}))\cdots I_w({\mathbf{y}}_K(\mathbf{p})){]}^t\end{array}\text{\hspace{1em}(3)}$$

然后用 ${\bar{\mathbf{i}}}_r$ 和 ${\bar{\mathbf{i}}}_w(\mathbf{p})$ 表示各自均值为 $0$ 的版本，通过对每个向量减掉对应的算术平均值获得。
然后我们提出了下面的方法来度量 warping transformation(扭曲变换) 的性能：
$$E_{ECC}(\mathbf{p})={\Vert \frac{{\bar{\mathbf{i}}}_r}{\Vert {\bar{\mathbf{i}}}_r\Vert }-\frac{{\bar{\mathbf{i}}}_w(\mathbf{p})}{\Vert {\bar{\mathbf{i}}}_w(\mathbf{p})\Vert }\Vert }^2\text{\hspace{1em}(4)}$$
其中 $\Vert \cdot \Vert$ 表示常规的欧式范数。

从（4）中可以看出我们的 criterion (准则) 对 bias (偏差) 和 gain (增益) 是 invariant (不变) 的。这也表明我们的度量对于任何 brightness (亮度) 和/或 constrast (对比度) 中的 photometric distortions (光度失真) 也是不变的。因此对于第一近似，我们可以完全忽略光度变换，而只关注几何。值得一提的是，我们的度量对于异常值表现出统计的鲁棒性。所有这些良好特性表明我们提出的 criterion (准则) 是对齐问题的合适的目标函数。

3.1 性能度量优化

指定了性能度量之后，我们继续进行最小化，以得到最佳参数值。
显然，最小化 $E_{ECC}(\mathbf{p})$ 等价于最大化下面的 enhanced correlation coefficient (ECC，增强的相关系数)：
$$\rho (\mathbf{p})=\frac{{\bar{\mathbf{i}}}_r^t{\bar{\mathbf{i}}}_w(\mathbf{p})}{\Vert {\bar{\mathbf{i}}}_r\Vert \Vert {\bar{\mathbf{i}}}_w(\mathbf{p})\Vert }={\hat{\mathbf{i}}}_r^t\frac{{\bar{\mathbf{i}}}_w(\mathbf{p})}{\Vert {\bar{\mathbf{i}}}_w(\mathbf{p})\Vert }\text{\hspace{1em}(5)}$$
其中，为了简单起见，我们记 ${\hat{\mathbf{i}}}_r=\frac{{\bar{\mathbf{i}}}_r}{\Vert {\bar{\mathbf{i}}}_r\Vert }$ 为均值为零的参考向量（${\bar{\mathbf{i}}}_r$）的归一化，这是 constant（常数）。

因为 reference（参考）图像是已知的？即 $\mathbf{x}$ 都是已知的，所以是 constant。

注意，尽管 ${\bar{\mathbf{i}}}_w(\mathbf{p})$ 线性地依赖参数向量 $\mathbf{p}$，由于扭曲向量的归一化，得到的目标函数相对于 $\mathbf{p}$ 仍然是非线性的。当然，这表明最大化需要非线性优化技术。

正如在第一章讲的，最大化 $\rho (\mathbf{p})$ 可以用 直接搜索 或者基于 梯度 的方法。这里我们准备用后者。按照迭代技术中的惯例，我们准备用一系列二次优化来代替原始优化问题。每个二次优化依赖于前一次优化的结果，从而产生一系列参数估计，这些参数估计有望收敛到我们想要的优化向量。每次迭代中，我们不必优化目标函数，而是对它进行 近似 。当然，必须选择近似值，以便优化器容易计算。下面我们来介绍对目标函数的近似，推导使其最大化的解。

假设 $\mathbf{p}$“接近”某个 nominal (标称) 参数向量 $\tilde{\mathbf{p}}$，记作 $\mathbf{p}=\tilde{\mathbf{p}}+\Delta \mathbf{p}$，其中 $\Delta \mathbf{p}$ 表示 perturbations (扰动) 向量。
令 $\tilde{\mathbf{y}}=\mathbf{\varphi }(\mathbf{x};\tilde{\mathbf{p}})$ 为标称参数向量下的 warped (扭曲的) 坐标， $\mathbf{y}=\mathbf{\varphi }(\mathbf{x};\mathbf{p})$ 是 perturbed (扰动的) 坐标。

考虑在 warped (扭曲的) 图像上坐标 $\mathbf{y}$ 处的 intensity (强度) ，对参数做一阶泰勒展开，有：
$$I_w(\mathbf{y})\approx I_w(\tilde{\mathbf{y}})+[{\nabla }_{\mathbf{y}}{I}_w(\tilde{\mathbf{y}}){]}^t\frac{\partial \mathbf{\varphi }(\mathbf{x};\tilde{\mathbf{p}})}{\partial \mathbf{p}}\Delta \mathbf{p}\text{\hspace{1em}(6)}$$

泰勒展开的公式为： $f(x)=f(x_0)+\frac{f^{\prime }(x_0)}{1!}(x-x_0)+\frac{f^{\prime \prime }(x_0)}{2!}(x-x_0{)}^2+\cdots +\frac{f^{(n)}(x_0)}{n!}(x-x_0{)}^n+R_n$
按照我的理解，这里从 $\bar{\mathbf{y}}$ 处展开（相当于是公式里的 $x_0$）， $\mathbf{y}$ 相当于公式里的 $x$，所以 $x-x_0=\Delta \mathbf{p}$。

其中 ${\nabla }_{\mathbf{y}}{I}_w(\tilde{\mathbf{y}})$ 表示 warped (扭曲) 图像的 intensity (强度) 函数 $I_w(\mathbf{y})$ 的长度为 $2$ 的梯度向量，这是在 nominal (标称) warped (扭曲的) 坐标 $\tilde{\mathbf{y}}$ 处 evaluated (评估) 的。

由于 $\mathbf{\varphi }(\mathbf{x};\mathbf{p})$ 是一个长度为 $2$ 的向量转换（为了生成 warped 坐标），$\frac{\partial \mathbf{\varphi }(\mathbf{x};\tilde{\mathbf{p}})}{\partial \mathbf{p}}$ 表示这个转换的尺寸为 $2\times N$ 的雅克比矩阵，这是在 nominal (标称) 参数值上评估的。

注意我们已经假设了intensity (强度) 函数$I_w$ 以及 warping transformation (扭曲变换) $\mathbf{\varphi }$ 都具有足够的光滑性，以允许偏导数的存在。

我们可以对目标区域 $\mathcal{T}$ 的所有坐标 ${\mathbf{x}}_k,k=1,\dots ,K$ 应用（6）。这样就产生了 warped vector (扭曲向量) （参数为 $\mathbf{p}$）的线性版本：
$${\mathbf{i}}_w(\mathbf{p})\approx {\mathbf{i}}_w(\tilde{\mathbf{p}})+G(\tilde{\mathbf{p}})\Delta (\mathbf{p})\text{\hspace{1em}(7)}$$
其中 $G(\tilde{\mathbf{p}})$ 表示 warped intensity vector (扭曲强度向量) 关于其参数的 $K\times N$ 雅克比矩阵，这是在 nominal (标称) 参数值 $\tilde{\mathbf{p}}$ 上评估的。为了准确地表述这个矩阵，我们假设 warping transformation (扭曲变换) 的形式为：
$$\mathbf{\varphi }(\mathbf{x};\mathbf{p})=[{\varphi }_1(\mathbf{x};\mathbf{p}),{\varphi }_2(\mathbf{x};\mathbf{p}){]}^t\text{\hspace{1em}(8)}$$
其中 ${\varphi }_1,{\varphi }_2$ 为标量函数。矩阵 $G$ 的 $(k,n)$ 元素可以写为：
$$G(\tilde{\mathbf{p}}{)}_{k,n}=\sum _{i=1}^2\left({\frac{\partial I_w(\mathbf{y})}{\partial y_i}|}_{\mathbf{y}={\mathbf{y}}_k(\tilde{\mathbf{p}})}\times {\frac{\partial {\varphi }_i({\mathbf{x}}_k;\mathbf{p})}{\partial p_n}|}_{\mathbf{p}=\tilde{\mathbf{p}}}\right)\text{\hspace{1em}(9)}$$
其中 $k=1,\dots ,K;n=1,\dots ,N$，我们回想一下 $\mathbf{y}=[y_1,y_2{]}^t$ 为 warped image 上的坐标。

我们现在需要计算 warped vector (扭曲向量) 零均值版本。在（7）的帮助下，我们得到（5）里面定义的目标函数 $\rho (\mathbf{p})$ 的近似：
$$\rho (\mathbf{p})\approx \rho (\Delta \mathbf{p}|\bar{\mathbf{p}})={\hat{\mathbf{i}}}_r^t\frac{{\bar{\mathbf{i}}}_w(\tilde{\mathbf{p}})+\bar{G}(\tilde{\mathbf{p}})\Delta \mathbf{p}}{\Vert {\bar{\mathbf{i}}}_w(\tilde{\mathbf{p}})+\bar{G}(\tilde{\mathbf{p}})\Delta \mathbf{p}\Vert }\text{\hspace{1em}(10)}$$
其中 $\bar{G}(\tilde{\mathbf{p}})$ 和 $\bar{\mathbf{i}}(\tilde{\mathbf{p}})$ 分别是 $G(\tilde{\mathbf{p}})$ 和 $\mathbf{i}(\tilde{\mathbf{p}})$ 的 column-zero-mean (列-零均值) 版本。

现在开始让我们简化一下符号的写法，去除 $\mathbf{p}$ 上对 warped vector (扭曲向量) 的依赖（意思就是省掉括号里的$\bar{\mathbf{p}}$），我们可以把上面的近似写成：

$$\rho (\Delta \mathbf{p}|\tilde{\mathbf{p}})=\frac{{\hat{\mathbf{i}}}_r^t{\bar{\mathbf{i}}}_w+{\hat{\mathbf{i}}}_r^t\bar{G}\Delta \mathbf{p}}{\sqrt{\Vert {\bar{\mathbf{i}}}_w{\Vert }^2+2{\bar{\mathbf{i}}}_w^t\bar{G}\Delta \mathbf{p}+\Delta {\mathbf{p}}^t{\bar{G}}^t\bar{G}\Delta \mathbf{p}}}\text{\hspace{1em}(11)}$$

尽管 $\rho (\Delta \mathbf{p}|\tilde{\mathbf{p}})$ 在 $\Delta \mathbf{p}$ 是非线性的，它的最大化很简单，结果是一个闭合形式的表达式。这是下一个定理的结论，它提供了必要的结果。

---

定理1

考虑标量函数：
$$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$ 是方形的，对称的，正定的，大小为 $N$ 的矩阵；
$v，\mathbf{v}，Q$ 满足：
$$v>{\mathbf{v}}^t{Q}^{-1}\mathbf{v}\text{\hspace{1em}(13)}$$
然后，就 $f(\mathbf{x})$ 的最大值而言，我们区分了以下两种情况：

$\bullet$ $u>{\mathbf{u}}^t{Q}^{-1}\mathbf{v}$ 时，我们有一个最大值，具体地：
$$\underset{\mathbf{x}}{\max }f(\mathbf{x})=\sqrt{\frac{(u-{\mathbf{u}}^t{Q}^{-1}\mathbf{v}{)}^2}{v-{\mathbf{v}}^t{Q}^{-1}\mathbf{v}}+{\mathbf{u}}^t{Q}^{-1}\mathbf{u}}\text{\hspace{1em}(14)}$$

它可以由这得到：
$$\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}(15)}$$

$\bullet$ $u\le {\mathbf{u}}^t{Q}^{-1}\mathbf{v}$ 时，我们有一个上界，它等于：
$$\underset{\mathbf{x}}{\sup }f(\mathbf{x})=\sqrt{{\mathbf{u}}^t{Q}^{-1}\mathbf{u}}\text{\hspace{1em}(16)}$$
这可以它来无限接近：
x = Q − 1 { λ u − v } (17) \mathbf{x} = Q^{-1} \left\{ \lambda \mathbf{u} - \mathbf{v} \right\} \tag{17} x=Q−1{λu−v}(17)
$\lambda$ 是足够大的正数（更准确地说，我们的意思是，对于所有 $ϵ>0$，存在足够大的标量 ${\lambda }_{ϵ}$，使得 $f(\mathbf{x})$ 为 $ϵ$ 接近上限）

证明见附录。

---

让我们来检验一下定理1是否可以用来最大化（11）里定义的 $\rho (\Delta \mathbf{p}|\tilde{\mathbf{p}})$。
为此，我们需要验证（13）的有效性。对于感兴趣的问题，可以转成如下不等式：$\Vert {\bar{\mathbf{i}}}_w{\Vert }^2>{\bar{\mathbf{i}}}_w^t{P}_G{\bar{\mathbf{i}}}_w$ ，其中 $P_G=\bar{G}({\bar{G}}^t\bar{G}{)}^{-1}{\bar{G}}^t$。这个关系式很容易满足，因为 $P_G$ 是一个正交投影操作（比方说 $P_G^2=P_G$以及$P_G^t=P_G$），于是我们可以写：
$$\Vert {\bar{\mathbf{i}}}_w{\Vert }^2=\Vert P_G{\bar{\mathbf{i}}}_w{\Vert }^2+\Vert [I-P_G]{\bar{\mathbf{i}}}_w{\Vert }^2\ge \Vert P_G{\bar{\mathbf{i}}}_w{\Vert }^2={\bar{\mathbf{i}}}_w^t{P}_G{\bar{\mathbf{i}}}_w\text{\hspace{1em}(18)}$$

其中 $I$ 表示单位矩阵。当 $[I-P_G]{\bar{\mathbf{i}}}_w=0$ 的时候我们得到一个等式，当 ${\bar{\mathbf{i}}}_w$ 是 $\bar{G}$ 的线性组合时为真。显然这种情况发生的可能性为零，特别是在有噪声的情况下。因此对于所有实际情况都是不等的。

由于我们可以应用定理1，根据（15），优化 perturbation (扰动) 等价于：
$$\Delta \mathbf{p}=({\bar{G}}^t\bar{G}{)}^{-1}{\bar{G}}^t\left\{\frac{\Vert {\bar{\mathbf{i}}}_w{\Vert }^2-{\bar{\mathbf{i}}}_w^t{P}_G{\bar{\mathbf{i}}}_w}{{\hat{\mathbf{i}}}_r^t{\bar{\mathbf{i}}}_w-{\hat{\mathbf{i}}}_r^t{P}_G{\bar{\mathbf{i}}}_w}{\hat{\mathbf{i}}}_r-{\bar{\mathbf{i}}}_w\right\}\text{\hspace{1em}(19)}$$

其中 ${\hat{\mathbf{i}}}_r^t{\bar{\mathbf{i}}}_w>{\hat{\mathbf{i}}}_r^t{P}_G{\bar{\mathbf{i}}}_w$；或者根据（17）有：
Δ p = ( G ˉ t G ˉ ) − 1 G ˉ t { λ i ^ r − i ˉ w } (20) \Delta\mathbf{p} = (\={G}^t \={G})^{-1} \={G}^t \{ \lambda \hat{\mathbf{i}}_r - \={\mathbf{i}}_w \} \tag{20} Δp=(GˉtGˉ)−1Gˉt{λi^r​−iˉw​}(20)
其中 ${\hat{\mathbf{i}}}_r^t{\bar{\mathbf{i}}}_w\le {\hat{\mathbf{i}}}_r^t{P}_G{\bar{\mathbf{i}}}_w$，其中 $\lambda$ 的选择要使 $\rho (\Delta \mathbf{p}|\tilde{\mathbf{p}})$ 满足 $\rho (\Delta \mathbf{p}|\tilde{\mathbf{p}})>\rho (0|\tilde{\mathbf{p}})$。换句话说，我们想要选择一个 perturbation (扰动) 可以增加相关性，并且使其非负。下面的引理给出了 $\lambda$ 可能的值。

---

引理1

令 ${\bar{\mathbf{i}}}_r^t{\bar{\mathbf{i}}}_w\le {\hat{\mathbf{i}}}_r^t{P}_G{\bar{\mathbf{i}}}_w$，并且给 $\lambda$ 定义两个值：
$${\lambda }_1=\sqrt{\frac{{\bar{\mathbf{i}}}_w^t{P}_G{\bar{\mathbf{i}}}_w}{{\hat{\mathbf{i}}}_r^t{P}_G{\hat{\mathbf{i}}}_r}}，{\lambda }_2=\frac{{\hat{\mathbf{i}}}_r^t{P}_G{\bar{\mathbf{i}}}_w-{\hat{\mathbf{i}}}_r^t{\bar{\mathbf{i}}}_w}{{\hat{\mathbf{i}}}_r^t{P}_G{\hat{\mathbf{i}}}_r}\text{\hspace{1em}(21)}$$
然后对于 $\lambda \ge {\lambda }_1$，我们有 $\rho (\Delta \mathbf{p}|\tilde{\mathbf{p}})>\rho (0|\tilde{\mathbf{p}})$；
对于 $\lambda \ge {\lambda }_2$，有 $\rho (\Delta \mathbf{p}|\tilde{\mathbf{p}})\ge 0$；
最后，对于 λ ≥ max ⁡ { λ 1 , λ 2 } \lambda \geq \max \{ \lambda_1, \lambda_2\} λ≥max{λ1​,λ2​}，这两个不等式都成立。

证明
通过把（20）里的 $\Delta \mathbf{p}$ 代入到（11），目标函数就变成下面这个关于 $\lambda$ 的式子：

$$f(\lambda )=\frac{\left({\hat{\mathbf{i}}}_r^t{\bar{\mathbf{i}}}_w-{\hat{\mathbf{i}}}_r^t{P}_G{\bar{\mathbf{i}}}_w\right)+\lambda {\hat{\mathbf{i}}}_r^t{P}_G{\hat{\mathbf{i}}}_r}{\sqrt{\left(\Vert {\bar{\mathbf{i}}}_w{\Vert }^2-{\hat{\mathbf{i}}}_w^t{P}_G{\bar{\mathbf{i}}}_w\right)+{\lambda }^2{\hat{\mathbf{i}}}_r^t{P}_G{\hat{\mathbf{i}}}_r}}\text{\hspace{1em}(22)}$$

容易验证 $f(\lambda )$ 的导数是非负的；因此 $f(\lambda )$ 是关于 $\lambda$ 递增的。这表明对于 $\lambda \ge {\lambda }_2$ 我们有 $f(\lambda )\ge 0$。
注意现在对于 $\lambda ={\lambda }_1$ 我们可以写：
$$f({\lambda }_1)=\frac{{\hat{\mathbf{i}}}_r^t{\bar{\mathbf{i}}}_w-{\hat{\mathbf{i}}}_r^t{P}_G{\bar{\mathbf{i}}}_w+\sqrt{\left({\bar{\mathbf{i}}}_w^t{P}_G{\bar{\mathbf{i}}}_w\right)\left({\hat{\mathbf{i}}}_r^t{P}_G{\hat{\mathbf{i}}}_r\right)}}{\Vert {\bar{\mathbf{i}}}_w\Vert }\ge \rho (0|\tilde{\mathbf{p}})\text{\hspace{1em}(23)}$$
这条不等式是在 ${\hat{\mathbf{i}}}_r^t{P}_G{\bar{\mathbf{i}}}_w$ 上应用 Schwartz inequality (瓦茨不等式) 的结果，回想 $P_G$ 是一个 orthogonal projection (正交投影) 操作。

附注
当 ${\bar{\mathbf{i}}}_w$ 接近 ${\bar{\mathbf{i}}}_r$ 时，主要用（19），于是对于 ${\bar{\mathbf{i}}}_w\approx {\bar{\mathbf{i}}}_r$，我们有 ${\hat{\mathbf{i}}}_r^t{\bar{\mathbf{i}}}_w\approx {\hat{\mathbf{i}}}_r^t{\bar{\mathbf{i}}}_r>{\hat{\mathbf{i}}}_r^t{P}_G{\bar{\mathbf{i}}}_r\approx {\hat{\mathbf{i}}}_r^t{P}_G{\bar{\mathbf{i}}}_w$，这很有意思。
但是，注意的是，如果一直坚持用（19），每当 ${\hat{\mathbf{i}}}_r^t{\bar{\mathbf{i}}}_w\le {\hat{\mathbf{i}}}_r^t{P}_G{\bar{\mathbf{i}}}_w$ 成立，我们得到一个负的相关 $\rho (\Delta \mathbf{p}|\tilde{\mathbf{p}})$（即使 $\rho (0|\tilde{\mathbf{p}})>0$ 时也为是的），它总是小于 $\rho (0|\tilde{\mathbf{p}})$。
换句话说，这种情况下我们不是增加了 correlation coefficient (相关系数) （这是我们所期望的），而是减少了它。这清晰地表明应该优选使用（20），把（21）中的 $\lambda$ 的值代进去。

---

3.2 Forward Additive ECC Iterative Algorithm（前向可加的ECC迭代算法）

现在我们把上述结果转为迭代的方案，以获得原始非线性优化问题的解。
假设参数向量的估计 ${\mathbf{p}}_{j-1}$ 可以从 $j-1$ 次迭代获得，我们可以计算 ${\bar{\mathbf{i}}}_w({\mathbf{p}}_{j-1})$ 和 $\bar{G}({\mathbf{p}}_{j-1})$；然后我们可以在 $\rho (\Delta {\mathbf{p}}_j|{\mathbf{p}}_{j-1})$ 的帮助下根据（10）来近似 $\rho (\mathbf{p})$，并且优化这个关于 $\Delta {\mathbf{p}}_j$ 的近似。这样会有参数更新规则 ${\mathbf{p}}_j={\mathbf{p}}_{j-1}+\Delta {\mathbf{p}}_j$。
如步骤 $S_4$，当更新向量 $\Delta {\mathbf{p}}_j$ 的范数小于某个预定的阈值 $T$ 时，停止迭代。迭代步骤如表1所示。我们把相应的算法称为 Forward Additive ECC (FA-ECC)

给定目标区域 $\mathcal{T}$ 中的 $K$ 个像素，以及长度为 $N$ 的 parameter vector (参数向量) 估计 ${\mathbf{p}}_{j-1}$，可以简单地估计出所提出的方案的每次迭代的复杂度。考虑到通常 $K\gg N$，我们发现步骤 $S_3$ 中最费计算量的部分是：用（19）或（20）中的式子计算 $\Delta {\mathbf{p}}_j$，我们可以看到这里我们需要计算 ${\bar{G}}^t\bar{G}$，这需要 $O(K{N}^2)$ 次操作。这使我们的算法变得复杂，而其它步骤每次迭代大概需要 $O(KN)$。

3.3 Inverse Compositional ECC Iterative Algorithm

当对齐问题限制为指定类型的参数模型时，可以设计计算更高效的版本，因为算法的某些部分可以离线计算。如果我们采用[19]中提出的方法，我们可以得到我们算法的 Inverse Compositional ECC (IC-ECC) (转置组成的ECC) 版本，在每次迭代中可以显著减少 $O(KN)$ 的计算量。

我们简要提一下，这个方法可以在 [3] [13] 中找到，它基于交换 ${\mathbf{i}}_w$ 和 ${\mathbf{i}}_r$ 的作用。这样矩阵 $G$ 就变成了 reference intensity vector (参考强度向量) 的雅可比矩阵，由于这个向量的 warping function (扭曲函数) 是 identity (恒等式) ，矩阵 $G$ 是常量，${\bar{G}}^t\bar{G}$ 可以离线计算。后者是计算复杂度可以减少一个量级的原因。

通过适当地修改我们的 FA-ECC 版本，我们可以很容易从表1获得我们可替代的算法版本 IC-ECC 的概览。

关于 inverse 算法（additive 和 compositional ）以及 forward compositional 算法，要指出的是它们只能被用于指定的 warps 类型。并且 inverse 算法比对应的 forward 版更容易收到噪声的影响。这个缺点限制了这类算法在实际中的应用。

3.4 Relation to Existing SSD-Based Measures

这一节我们用一种不同的方式来导出我们的性能度量。这也将帮助我们将它与文献中目前最流行的两种 SSD 方法联系起来。对于我们的分析，我们将假设 photometric distorition (光度失真) 仅限于全局亮度(brightness) 和对比度 (contrast) 的变化，在这种简单的光度变化下，我们给我们的参数对齐问题定义以下的性能度量：
$$E(\mathbf{p},\alpha )=\Vert {\alpha }_1{\mathbf{i}}_w(\mathbf{p})+{\alpha }_2-{\mathbf{i}}_r{\Vert }^2\text{\hspace{1em}(24)}$$

其中 $\alpha =[{\alpha }_1,{\alpha }_2{]}^t$ 是 photometric transformation (光度转换) 的参数向量。我们的目标是最小化关于所有参数的目标函数。对于第一个光度参数，我们必须指出 ${\alpha }_1$ 的负值会产生 inversion (倒置?) 的影响，其颜色会反转。因此，如果存在一种先验知识：不能发生这种颜色反转，那么将 ${\alpha }_1$ 的值限制为正，是合乎逻辑的。现在，如果我们先最小化关于 ${\alpha }_1,{\alpha }_2$ 的目标函数，我们得到下面有趣的结果：
E ( p ) = min ⁡ α 1 ≥ 0 , α 2 E ( p , α ) = ∥ i ˉ r ∥ 2 { 1 − [ max ⁡ { ρ ( p ) , 0 } ] 2 } (25) E(\mathbf{p}) = \min_{\alpha_1 \geq 0, \alpha_2} E(\mathbf{p}, \alpha) = \| \={\mathbf{i}}_r \|^2 \left\{1 - [\max \{ \rho(\mathbf{p}), 0\}]^2 \right\} \tag{25} E(p)=α1​≥0,α2​min​E(p,α)=∥iˉr​∥2{1−[max{ρ(p),0}]2}(25)

其中 $\rho (\mathbf{p})$ 是（5）中定义的 correlation function (相关函数) 。注意由于 reference image (参考图像) 是恒定的，所以范数 $\Vert {\bar{\mathbf{i}}}_r{\Vert }^2$ 从上一次关联中得到，于是接下来最小化 $\mathbf{p}$ 等价于最小化 1 − [ max ⁡ { ρ ( p ) , 0 } ] 2 1 - [\max \{ \rho(\mathbf{p}), 0\}]^2 1−[max{ρ(p),0}]2。但是这个式子在 $\rho (\mathbf{p})$ 中正在减小，因此我们可以等价地最大化相关函数 $\rho (\mathbf{p})$，这又符合我们的准则了。最后这个优化问题很有意义。确实，注意由于 $\rho (\mathbf{p})$ 是光度失真无关的（则立我们考虑简单的情况），并且基于没有颜色反转的认知，寻找最正相关的关系是很合理的。

如果我们去掉约束 ${\alpha }_1\ge 0$，则（25）中目标函数的最小化问题就变成了 Fuh 和 Maragos 提出的最优化问题[6]。通过首先优化 ${\alpha }_1$ 和 ${\alpha }_2$ 得到：
E FM ( p ) = min ⁡ α 1 , α 2 E ( p , α ) = ∥ i ˉ r ∥ 2 { 1 − ρ 2 ( p ) } (26) E_{\text{FM}}(\mathbf{p}) = \min_{\alpha_1, \alpha_2} E(\mathbf{p}, \alpha) = \|\={\mathbf{i}}_r\|^2 \{ 1 - \rho^2 ({\mathbf{p}}) \} \tag{26} EFM​(p)=α1​,α2​min​E(p,α)=∥iˉr​∥2{1−ρ2(p)}(26)
注意现在得到的度量是关于 $|\rho (\mathbf{p})|$ 的减函数；因此接下来关于 $\mathbf{p}$ 的最小化问题等价于最大化 correlation function (相关函数) 的绝对值 $|\rho (\mathbf{p})|$。很明显，这个优化问题没有考虑到 “没有颜色反转”(no color inversion ) 的先验知识。在 [6] 中是通过在 $N$-$\text{D}$ 量化参数空间中 ($N$-$\text{D}$ 就是 $N$ 维的意思) 通过穷举搜索的方法实现最大化的。显然，在没有颜色反转的情况下，这样的搜索将会产生正确的最大正相关（当然，前提是扭曲图像不包含目标区域中的负相关部分）。然而正如我们在第1节提到的，穷举搜索方法的特点是计算复杂度高，当我们要求亚像素的精度时，这就变得非常高了。

我们还可以用迭代的方法，类似于我们提出的。然而，如果我们试图通过和（10）相同的近似方法来最大化 $|\rho (\mathbf{p})|$ ，然后可以证明最优扰动（perturbation） $\Delta \mathbf{p}$ 总是由（19）给出。正如我们附注（引理1后面）指出的，用这种策略可能导致 $\rho (\mathbf{p})$ 的局部最小值的负相关，而不是所期望的最大值。换句话说，迭代的方法更容易被锁定在错误的局部极值中。

如果在（24）中我们交换 ${\mathbf{i}}_w$ 和 ${\mathbf{i}}_r$ ，则会出现另一种度量方法，即：
$$E_{\mathbf{p},\mathbf{\alpha }}=\Vert {\alpha }_1{\mathbf{i}}_r+{\alpha }_2-{\mathbf{i}}_w(\mathbf{p}){\Vert }^2\text{\hspace{1em}(27)}$$
这个是 Lucas 和 Kanade 提出的方法[10]（L-K光流法），总所周知，它可以生成许多在实际中被广泛应用的算法变体。按照前面两种情况类似的步骤，我们首先最小化两个光度参数，有：
E LK ( p ) = min ⁡ α 1 , α 2 E ( p , α ) = ∥ i ˉ w ( p ) ∥ 2 { 1 − ρ 2 ( p ) } (28) E_{\text{LK}} (\mathbf{p}) = \min_{\alpha_1, \alpha_2} E(\mathbf{p}, \boldsymbol{\alpha}) = \| \={\mathbf{i}}_w (\mathbf{p}) \|^2 \{ 1 - \rho^2 ({\mathbf{p}}) \} \tag{28} ELK​(p)=α1​,α2​min​E(p,α)=∥iˉw​(p)∥2{1−ρ2(p)}(28)
我们从当前 criterion(标准) 结果中观察到，这个有两项是依赖于参数 $\mathbf{p}$ 的，即很熟悉的部分 { 1 − ρ 2 ( p ) } \{ 1 - \rho^2 ({\mathbf{p}}) \} {1−ρ2(p)}，以及 warped image(扭曲图像) 的 magnitude(量级) $\Vert {\bar{\mathbf{i}}}_w(\mathbf{p}){\Vert }^2$（这个不是 constant 的）。因此最小化 $E_{\text{LK}}(\mathbf{p})$ 就涉及到这两项的组合的最小化。
第一个观察是， 这个 criterion(标准) 不一定会产生和我们的 measure(度量) 相同的解。
第二，由于这个 $\Vert {\bar{\mathbf{i}}}_w(\mathbf{p}){\Vert }^2$，显然迭代算法会陷入 $\Vert {\bar{\mathbf{i}}}_w(\mathbf{p}){\Vert }^2\approx 0$ 时的解（例如具有均匀强度(intensity)的区域）。
第三，由于 ${\rho }^2(\mathbf{p})$，算法会陷入负相关(negative correlations)。

尽管有上述观察，但是 Lucas-Kanade performance measure 还是产生了图像对齐问题的最流行的迭代算法。因此我们将把它作为参考，与我们的方案进行对比。因此我们来详细讲一下 Forward Additive LK（FA-LK） 的更新版本。代入（28）中的 ${\bar{\mathbf{i}}}_w(\mathbf{p})$ 的线性近似，然后相对于 $\Delta \mathbf{p}$ 进行最小化，我们得到以下最优更新扰动 (optimum updating perturbation)：
$$\Delta {\mathbf{p}}_{\text{LK}}=({\bar{G}}^t\bar{G}{)}^{-1}{\bar{G}}^t\left\{\frac{{\hat{\mathbf{i}}}_r^t{\bar{\mathbf{i}}}_w-{\hat{\mathbf{i}}}_r^t{P}_G{\bar{\mathbf{i}}}_w}{1-{\hat{\mathbf{i}}}_r^t{P}_G{\hat{\mathbf{i}}}_r}{\hat{\mathbf{i}}}_r-{\bar{\mathbf{i}}}_w\right\}\text{\hspace{1em}(29)}$$

这个在任何时候都是适用的。将（19）与（29）进行比较，我们意识到区别仅仅是向量 ${\hat{\mathbf{i}}}_r$ 前面的标量。正如我们将看到的，这个看似微小的变化，结合了（20）之后，将会带来显著的性能提升。

对于 Lucas-Kanade 方法，可以定义一种特殊的基于 SSD 的度量 (measure)，该度量可以处理任意的线性外观变化。为了最小化它，Hager 和 Belhumeur 在 [3] 中提出了一种利用逆加性到更新规则(inverse additive update rule) 的算法。基于同样的 SSD 度量 (measure)，Baker [19] 等人通过采用逆合成方法(inverse compositional approach) 提出了几个 Hager-Belhumeur 算法的变种。在这些替代算法方案中，根据报告，SIC 算法性能最好 [19]。因此下面也会测试该算法。

---

Reference

[3] Efficient Region Tracking with Parametric Models of Geometry and Illumination.
[6] Motion Displacement Estimation Using an Affine Model for Image Matching.
[10] An Iterative Image Registration Technique with an Application to Stereo Vision.
[13] Lucas-Kanade 20 Years On: A Unifying Framework: Part 1. The Quantity Approximated, the Warp Update Rule, and the Gradient Descent Approximation.
[19] Lucas-Kanade 20 Years On: A Unifying Framework: Part 3.