随着各种医学成像技术的快速发展, 多模态图像配准在医学诊断和治疗方面发挥着越来越重要的作用. 它可以结合多种成像技术的优点进而为医生所关注的部位提供更加全面、准确的信息^[1]. 例如, 核磁共振 (MR) 图像与同一病人的 CT 图像重新对齐, 以精确定位肿瘤, 进行手术规划. 所谓图像配准, 是指寻找一种空间变换, 使得待配准图像与目标图像的对应点达到空间上的一致. $\Omega$ 是 $\mathbb{R}^2$ 上的一个有界开集, $T$, $F$$: \Omega \to \mathbb{R}$ 分别被称作目标图像和待配准图像. 图像配准的目的是寻找一个空间变换 ${\mathbf {h}}\left( {\mathbf {x}} \right) \buildrel \Delta \over = {\mathbf {x}} + {\mathbf{u}}\left( {\mathbf{x}} \right)$, 其中 $\mathbf x$ 表示恒等部分, ${\mathbf{u}}({\mathbf{x}})$ 为位移场, 使得 $F({\mathbf {h}}( \cdot ))$ 与 $T( \cdot )$ 能够尽可能地相似. 以此为基础, 本文将图像配准问题转化成能量泛函的优化问题, 即寻找最优的位移场 ${\mathbf{u}}({\mathbf{x}})$, 使得以下能量泛函达到最小

其中 $S({\mathbf{u}}) = dis(F({\mathbf{x}} + {\mathbf{u}}({\mathbf{x}})),T({\mathbf{x}}))$ 表示目标图像 $T$ 与浮动图像 $F$ 之间的相似性测度, 即保真项; ${R_0}({\mathbf{u}})$ 表示正则项, 用于对位移场 $\mathbf{u}$ 作先验估计, 以排除不合理的解; $\delta>0$ 用来平衡正则项和保真项.

给定不同的相似性测度和正则项, 产生了许多不同的模型. 常见的相似性度量有图像之间的平方差之和 (sum of squared differences, SSD)^[2]、 互信息 (mutual information, MI)^[3,4,5]、Kullback-Leibler (KL) 散度^[6] 以及其它的一些度量方式^[7]. MI 是多模态图像配准中描述图像对之间相似性最常用的度量方式之一, 而极大化 MI 的方法最大的困难在于需要估计联合概率密度. 这使得计算变得较为复杂, 增大了 CPU 的消耗. 对此, Shi 等^[8] 提出了一种新的度量—— 瑞利度量来简化模型, 并取得了不错的结果. 然而, 文献 [8] 中并没有考虑网格重叠问题. 网格重叠会产生负面积或负体积, 这与物理原理相矛盾. 因此, 在多模态图像配准中引入限制条件来避免网格重叠是非常有必要的. 在此应用需求下, 需要保证形变场 ${\mathbf{h}}:\Omega \to \Omega$ 是一个双射. 根据逆映射定理^[9], 需要保证对 $\forall \mathbf{x} \in \Omega$, 形变场 $\mathbf{h}(\mathbf{x})$ 的雅克比行列式 $\det (\nabla {\mathbf{h}}({\mathbf{x}})) > 0$. 为了避免网格重叠, Zhang 等^[10] 提出了不等式约束来避免网格重叠. 基于他们的工作, Han 等^[11] 提出了一个带有 Cauchy-Riemann 约束条件的松弛的微分同胚图像配准模型, 进而将位移场 $\mathbf{u}$ 限制在一个二维共形集上. 之后, Han 等^[12] 为该微分同胚模型提出了一种快速算法, 并取得了许多令人满意的结果.

然而, 这些模型都是在对位移场 $\mathbf{u}$ 作先验估计的前提下进行的, 即通过极小化 $S({\mathbf{u}}) + \delta {R_0}({\mathbf{u}})$ 来得到一个适当光滑的位移场 $\mathbf{u}$. 但是图像配准的最终目的是寻找 $S(\mathbf{u})$ 具有适当光滑性的全局极小值. 为了能够在没有正则化的情况下直接寻找 $S({\mathbf{u}})$ 的全局极小值点, Modin 等^[13]提出了一种基于大形变微分度量映射 (LDDMM) 的分层图像配准模型, 其本质是一个带有约束的变分问题, 但是成本函数过于复杂, 且没有具体实验证明所提多尺度方法的有效性. 为了简化模型, Han 等^[14]提出了一种二维多尺度图像配准模型, 该方法实现了二维微分同胚图像配准的最优解. 但是在文献 [11,12,14] 中所选取的相似性度量为 SSD, 使用该度量方法的前提是相同组织的灰度值是相同的, 即仅适用于单模态图像配准. 作为上述工作的一个推广, 本文给出如下基于瑞利度量的多尺度多模态微分同胚图像配准模型

其中 $E({\mathbf{u}}) = \lambda S({\mathbf{u}}) + \delta {R_0}({\mathbf{u}})$, $\lambda > 0$, $S({\mathbf{u}}) = dis(F({\mathbf{x}} + {\mathbf{u}}({\mathbf{x}})),T({\mathbf{x}})) = \left| \Omega \right|(1 - {\rm MCC}(F({\mathbf{x}} + \\ {\mathbf{u}}({\mathbf{x}})),T({\mathbf{x}})){)^2}$ $\left| \Omega \right|$ 表示配准图像区域的面积, ${\rm MCC}(F( \cdot ),T( \cdot ))$ 表示目标图像和待配准图像之间的最大相关系数, 具体定义将在下一节中给出, $\Omega \buildrel \Delta \over = ({a_1},{b_1}) \times ({a_1},{b_1}),{R_0}({\mathbf{u}}) = \int_\Omega {{{\left| {{\nabla ^\alpha }{\mathbf{u}}} \right|}^2}{\rm d}{\mathbf{x}}},$ $\delta > 0,$ $\Psi = \{ {{\mathbf{u}} = {{({u_1},{u_2})}^T} \in {{[H_0^\alpha (\Omega )]}^2}} :\frac{{\partial {u_1}}}{{\partial {x_1}}} = \frac{{\partial {u_2}}}{{\partial {x_2}}}$, $\frac{{\partial {u_1}}}{{\partial {x_2}}} = - \frac{{\partial {u_2}}}{{\partial {x_1}}}\},{{\rm A}_\varepsilon } = \{ {\mathbf{u}} = {({u_1},{u_2})^T} \in \Psi :\det (\nabla ({\mathbf{x}} + {\mathbf{u}}({\mathbf{x}}))) < \varepsilon,\varepsilon > 0\}$ $\alpha > 2$ 目的在于使得 $\mathbf{u}$ 的一阶偏导数有意义, 加入参数 $\lambda$ 是为了引入多尺度方法.

注1.1 根据逆映射定理及 $\det (\nabla {\mathbf{h}}({\mathbf{x}}))$ 的实际意义, $\det (\nabla {\mathbf{h}}({\mathbf{x}}))$ 需要满足条件

${{\rm A}_\varepsilon }$ 的引入是为了排除对 $\forall {\mathbf{x}} \in \Omega$, $\det (\nabla {\mathbf{h}}({\mathbf{x}})) = {(1 + \frac{{\partial {u_1}}}{{\partial {x_1}}})^2} + {(\frac{{\partial {u_1}}}{{\partial {x_2}}})^2} = 0$ 的解, 进而保证 $\mathbf{h}$ 的逆映射存在.

本文拟解决的问题主要分为以下两点: 一、如何在多模态图像配准中避免网格重叠; 二、如何寻找 $S(\mathbf{u})$ 的全局极小值点. 而多尺度方法本质上是多个函数复合的过程, 可以在某个适当的函数空间 $\Psi$ 中逼近变分泛函 $S(\mathbf{u})$ 的下确界, 进而得到一个较为精确的配准结果, 这也是引入多尺度方法的原因.

本文在第2章中介绍了相似性度量的相关概念, 给出了基于多尺度方法的多模态图像配准模型, 并简要证明了模型的解的存在性和收敛性; 在第3章中提出了一种松弛的多尺度模型及其对应的数值计算; 第4章通过实验来论证所提算法的有效性, 在文章的最后对全文进行了总结.

在介绍模型之前, 简要回顾一下本文所用的相似性度量——瑞利度量的相关概念. 在 1959 年, Rényi^[15] 针对任意两个随机变量 $X$, $Y$ 之间的相似性度量 $Q(X,Y)$, 进行了如下假设

(1) 对任意的两个随机变量 $X$ 和 $Y$, $Q(X,Y)$ 都不是概率为 1 的常数;

(2) $Q(X,Y)=Q(Y,X)$;

(3) $0 \le Q(X,Y) \le 1$;

(4) $Q(X,Y)=0$ 当且仅当 $X$ 和 $Y$ 是相互独立的;

(5) $Q(X,Y) = 1$ 当且仅当 $X = g(Y)$ 或 $Y = f(X)$, 其中 $f$ 和 $g$ 都是 Borel 可测函数.

最大相关系数 ${\rm MCC}(X,Y)$ 很好地符合了上述假设, 其表达式为

其中 $\Upsilon (\mathbb{R})$ 表示 $\mathbb{R}$ 上所有具有有限正方差的 Borel 可测函数组成的空间, ${\rm CC}(f(X),g(Y))$ 表示 $f(X)$ 和 $g(Y)$ 的相关系数, 其表达式为

其中 $Cov( \cdot, \cdot)$ 和 $Var( \cdot )$ 分别表示协方差和方差, 定义为

其中 $E(f(X)) = \frac{{\int_\Omega {f(X)d{\bf{x}}} }}{{\left| \Omega \right|}}$, $E(g(Y)) = \frac{{\int_\Omega {g(Y)d{\bf{x}}} }}{{\left| \Omega \right|}}$.

对每一个 ${\mathbf{x}} \in \Omega$, $F({\mathbf{x}})$ 和 $T({\mathbf{x}})$ 可以看作两个随机变量. 基于假设 (4) 和 (5) 可得, ${(1 - {\rm MCC}(F({\mathbf{x}}),T({\mathbf{x}})))^2}$ 越小, 表明 $F({\mathbf{x}})$ 和 $T({\mathbf{x}})$ 的相关性越强. 因此, 本文的最终目标是极小化 ${(1 - {\rm MCC}(F({\mathbf{x}}),T({\mathbf{x}})))^2}$. 另一方面, 从上述公式中可以发现, 采用极大化瑞利度量, 即极小化 ${(1 - {\rm MCC}(F({\mathbf{x}}),T({\mathbf{x}})))^2}$ 的方式不需要估计联合概率密度函数, 极大地简化了计算.

注2.1 ${\rm MCC}(F({\mathbf{x}}),T({\mathbf{x}}))$ 不同于 ${\rm CC}(F({\mathbf{x}}),T({\mathbf{x}}))$, 后者需要假设 $F(\mathbf{x})$ 和 $T(\mathbf{x})$ 线性相关, 而前者只需要在所有具有有限正方差的 Borel 可测函数空间中找到合适的 $f$ 和 $g$ 即可, 即使 $F(\mathbf{x})$ 和 $T(\mathbf{x})$ 不线性相关, ${\rm MCC}(F({\mathbf{x}}),T({\mathbf{x}}))$ 仍适用.

基于上述瑞利度量的相关知识, 本文采用 $n$ 尺度方法来寻找能量泛函 $S(\mathbf{u})$ 的全局极小值. 该方法表述如下

第 1 步 求解以下变分问题的解

其中 ${\mathbf{h}}({\mathbf{x}}) \buildrel \Delta \over = {\mathbf{x}} + {\mathbf{u}}({\mathbf{x}})$, ${{\mathbf{\tilde h}}_0}({\mathbf{x}}) \buildrel \Delta \over = {\mathbf{x}} + {{\mathbf{u}}^0}({\mathbf{x}})$ 为第 1 步获得的形变, ${\lambda _0} > 0$, ${\varepsilon _0} > 0$.

第 2 步 将 $F \circ {{\mathbf{\tilde h}}_0}( \cdot )$ 作为新的待配准图像, 再与目标图像 $T( \cdot )$ 进行配准. 按照上述同样的方法求解以下变分问题的解

其中 ${{\mathbf{h}}_1}({\mathbf{x}}) \buildrel \Delta \over = {\mathbf{x}} + {{\mathbf{u}}^1}({\mathbf{x}})$ 为第 2 步获得的形变, ${{\mathbf{\tilde h}}_1}({\mathbf{x}}) \buildrel \Delta \over = {{\mathbf{\tilde h}}_0} \circ {{\mathbf{h}}_1}({\mathbf{x}})$ 为前两步复合获得的形变, ${\lambda _1} > 0$, ${\varepsilon _1} > 0$.

第 $n+1$ 步 依次进行, 最后寻找 $\mathbf{u}^n$ 满足以下条件

其中 ${E_n}({\mathbf{u}}) = {\lambda _n}\left| \Omega \right|{(1 - {\rm MCC}(F \circ {{\mathbf{\tilde h}}_{n - 1}}({\mathbf{x}} + {\mathbf{u}}({\mathbf{x}})),T({\mathbf{x}})))^2} + \delta {R_0}({\mathbf{u}})$, ${{\mathbf{h}}_n}({\mathbf{x}}) \buildrel \Delta \over = {\mathbf{x}} + {{\mathbf{u}}^n}({\mathbf{x}})$ 为第 $n+1$ 步获得的形变, ${{\mathbf{\tilde h}}_n}({\mathbf{x}}) \buildrel \Delta \over = {{\mathbf{\tilde h}}_0} \circ {{\mathbf{h}}_1} \circ \cdots \circ {{\mathbf{h}}_n}({\mathbf{x}})$ 为前 $n+1$ 步复合获得的形变, ${\lambda _n} > 0$, ${\varepsilon _n} > 0$.

关于变分模型 (2.3) 解的存在性, 有以下结论成立

定理2.1 假设 ess ${\sup\limits _{\mathbf x \in \Omega }}\left | {F(\mathbf x)} \right| \le M$, ess ${\sup\limits _{\mathbf x \in \Omega }}\left | {T(\mathbf x)} \right| \le M(M < + \infty )$, 对 $\forall {\mathbf{x}} \in \Omega$, $F(\mathbf{x})$, $T(\mathbf x)$ 不恒为常数, 且 $F$ 在 $\mathbf{x}$ 处的不连续点构成一个零测度集, 则 (2.3) 式的解存在.

证 定理 2.1 中解的存在性证明是变分法中的一个标准过程, 该过程主要分为两步: 第一步, 证明 $E_n(\mathbf{u})$ 的弱下半连续性; 第二步, 根据文献 [16] 中第 2.1.2 节, 即可得 (2.3) 式的解的存在性. 具体证明过程与文献 [11] 中的定理 2.2 的证明类似, 因此不再予以详细证明.

由 (2.3) 式得, 对 $\forall n \in {N^ + }$, ${E_n}({{\mathbf{u}}^n}) \le {E_n}({\mathbf{0}})$, 即

则

进而可得 $\{ {(1 - {\rm MCC}(F \circ {{\mathbf{\tilde h}}_n}( \cdot ),T( \cdot )))^2}\}$ 是一个单调递减有下界的序列.

不妨令 ${\varphi _0} = \mathop {\lim }\limits_{n \to + \infty } {(1 - {\rm MCC}(F \circ {{\mathbf{\tilde h}}_n}( \cdot ),T( \cdot )))^2}$ 及 ${\gamma _0} = {\inf\limits _{{\mathbf{u}} \in \Psi }}(1 - {\rm MCC}(F({\mathbf{x}} + {\mathbf{u}}({\mathbf{x}})), T({\mathbf{x}})){)^2}$. 关于在集合 $\Psi$ 上是否能找到 $S(\mathbf{u})$ 的全局极小值点问题, 有以下结论成立

定理2.2 ${{\mathbf{h}}_n}$, ${{\mathbf{\tilde h}}_n}$ 是由多尺度方法 (2.1)-(2.3) 产生的形变场, $B = B(\Omega )$, $M$, $\lambda_n$ 是三个充分大的数, 其中 $B$ 是依赖于 $\Omega$ 的数. 若 $\mathop {\lim }\limits_{n \to + \infty } \frac{{{B^{4n - 3}}{M^{{4^n}}}}}{{{\lambda _n}}} = 0,\mathop {\lim }\limits_{n \to + \infty } {\varepsilon _n} = 0$, 则有 ${\varphi _0} = {\gamma _0}$.

证 由文献 [14] 中引理 2.6-引理 2.10 可得, 微分同胚映射的复合仍为微分同胚映射, 微分同胚映射的逆映射也为微分同胚映射. 在此基础上, 可以将定理 2.2 的证明概述为以下步骤: 第一步, 证明 ${\varphi _0} \ge {\gamma _0}$. 基于 ${\gamma _0}$ 的定义, ${\varphi _0} \ge {\gamma _0}$ 显然成立; 第二步, 证明 ${\varphi _0} \le {\gamma _0}$, 可以直接通过反证法和相关的泛函分析理论推得矛盾. 该证明思想与文献 [14] 中定理 2.3 证明本质上相同, 因此本文不再重述.

注2.2 $\lambda _n$ 是满足定理 2.2 中条件的单调递增的序列, 本文令 $\lambda _n = \lambda _0 \times {a^{n - 1}}(n \in {\mathbb{N}^ + })$. 随着 ${\lambda _n}$ 的增大, 保真项 $S(\mathbf u)$ 逐渐占据主导地位, 进而逼近 $S(\mathbf u)$ 的下确界, 实现“贪婪配准”.

注2.3 由定理 2.2 可得, 多尺度方法 (2.1)-(2.3) 是收敛的, 当 $n$ 充分大时, 多尺度方法 (2.1)-(2.3) 与以下变分问题等价

与 (1.2) 式相比, (2.4) 式与参数 $\lambda$, $\alpha$, $\delta$ 无关, 说明了多尺度方法具有良好的鲁棒性这一优势. 进一步通过极小化 (2.4) 式可以直接寻找 $S(\mathbf{u})$ 的全局极小值, 这也是本文引入多尺度方法的主要原因.

为了进一步简化模型和计算, 本文将模型 (2.3) 转化成以下无约束变分问题

其中 $\Theta$ 是一个充分大的数, $R_1(\mathbf{u}) =\int_\Omega {[(\frac{{\partial {u_1}}}{{\partial {x_1}}}} - \frac{{\partial {u_2}}}{{\partial {x_2}}}{)^2} + {(\frac{{\partial {u_1}}}{{\partial {x_2}}} + \frac{{\partial {u_2}}}{{\partial {x_1}}})^2}]{\rm d}{\mathbf{x}}$, $E_n(\bf u)$ 的定义在模型 (2.3) 中已给出.

注3.1 关于模型 (3.1) 解的存在性可以采用定理 2.1 的方法进行证明. 同样地, 令 $R({\mathbf{u}}) = \delta {R_0}({\mathbf{u}}) + \Theta {R_1}({\mathbf{u}})$, 可以采用定理 2.2 同样的方法来证明该松弛多尺度模型的收敛性. 进一步地, 基于该收敛性, 可以得出当 $n$ 充分大时, 模型 (3.1) 也等价于 (2.4) 式, 故仍与参数无关, 且能得到泛函 $S(\mathbf{u})$ 的全局极小值. 值得注意的是, $\Theta$ 对于整体和局部都具有十分重要的意义: 对于整个多尺度方法而言, $\Theta$ 充分大的目的在于将 $\mathbf{u}$ 近似地控制到集合 $\Psi$ 中; 而对于每个尺度而言, $\Theta$ 充分大可以消除网格重叠. 若满足此条件, 则 $\Theta$ 对最终配准结果没有影响.

注3.2 一般情况下, 在图像配准中极小化 ${\lambda _n}S\left( {\bf{u}} \right) + \Theta {R_1}({\mathbf{u}})$ 也能得到较为光滑的形变场. 因此, 为了进一步简化模型和计算, 本文令 $\delta = 0$, 这足以得到较为光滑的解. 另外, 在一些极端情况下, 如非常大的形变等, 设置 $\delta >0$ 有助于得到更平滑的配准结果.

根据瑞利度量的定义, 本文的任务是需要在整个具有有限正方差的 Borel 可测函数空间中寻找合适的 $f$ 和 $g$, 这并不是一件容易的事. 对此, Rényi 等人^[8] 将对 $f$ 和 $g$ 的寻找限制在更小的函数空间内, 即再生核希尔伯特空间 (RKHS). 在这个更小的空间内, $f(F( \cdot ))$ 和 $g(T( \cdot ))$ 的相关系数的上确界仍然保持不变. 利用再生核希尔伯特空间理论, 在该空间的任意一个函数都可以表示成有限个核函数的线性组合, 且这个核函数是连续、对称、正定的, 在无穷远处值为 0. 在本文中, 采用高斯核函数作为光滑函数对 RKHS 中的函数进行逼近, 其定义为

其中 $\delta$ 表示带宽, 控制核函数的作用范围, $x$ 表示 $\Omega$ 中任意一点, $y$ 表示核函数的中心. 则 RKHS 中任意函数 $f(x)$ 和 $g(y)$ 都可以近似地表示为以下形式

其中 ${\alpha _s}$, ${\beta _{q}}$ 表示系数, $d$, $l$ 为有限正整数, ${\varepsilon _s}$ 和 ${\eta _{q}}$ 表示灰度值的量值. 进而将寻找 $f$ 和 $g$ 的问题转化为寻找有限个系数 ${\alpha _s}$ 和 ${\beta _{q}}$, 其个数取决于 ${\varepsilon _s}$ 和 ${\eta _{q}}$^[1].

基于上述描述, 模型 (3.1) 进一步地可以改写成

其中

${E_n}({{\mathbf{u}}^n},{{ {\alpha }}^n},{{ {\beta }}^n}) \buildrel \Delta \over = {\lambda _n}S({{\mathbf{u}}^n},{{ {\alpha }}^n},{{ {\beta }}^n}) +$ $\Theta {R_1}({{\mathbf{u}}^n}),$ $S({{\mathbf{u}}^n},{{ {\alpha }}^n},{{ {\beta }}^n}) = \left| \Omega \right|(1 -{\rm MCC}(F({\mathbf{x}} +$ ${{\mathbf{u}}^n}({\mathbf{x}})), T({\mathbf{x}})){)^2},$ ${\rm MCC}(F({\mathbf{x}} + {{\mathbf{u}}^n}{\mathbf{(x}})),T({\mathbf{x}})) =$ $\sup {\rm CC}(f(F({\mathbf{x}} + {{\mathbf{u}}^n}({\mathbf{x}}))),g(T({\mathbf{x}}))) =$ $\mathop {\sup }\limits_{\alpha _{_s}^n,\beta _{_{q}}^n} {\rm CC}(S({\mathbf{x}}), D({\mathbf{x}})),$ $S({\mathbf{x}}) \buildrel \Delta \over =$ $f(F({\mathbf{x}} + {{\mathbf{u}}^n}({\mathbf{x}}))) =$ $\sum\limits_{1 \le s \le d} {{\alpha _s}} K(F({\mathbf{x}} + {{\mathbf{u}}^n}({\mathbf{x}})),{\varepsilon _s}),$ $D({\mathbf{x}}) \buildrel \Delta \over = g(T({\mathbf{x}})) =$ $\sum\limits_{1 \le q \le l} {{\beta _{q}}K(T({\mathbf{x}}),{\eta _{q}})}, \quad {{ {\alpha }}^n} =$ ${(\alpha _s^n)_{1 \times d}}, {{ {\beta }}^n} = {(\beta _{_{q}}^n)_{1 \times l}}.$

通过交替极小化方法, 可以将模型 (3.2) 转化成以下三个子问题

采用与文献 [12] 中定理 2.1 相同的方法, 根据变分基本引理, 可以计算出上述三个子问题对应的 Euler-Lagrange 方程分别为

其中 $G({{\mathbf{u}}^n}) = [(D - E(D))Var(S) - (S - E(S))Cov(S,D)]$ $\cdot \sum\limits_s {{\alpha _s}\nabla K(F({\mathbf{x}}$ $ + {{\mathbf{u}}^n}({\mathbf{x}})),T({\mathbf{x}}))}$ $\cdot Va{r^{ - \frac{3}{2}}}(S)Va{r^{ - \frac{1}{2}}}(D),$ ${M_s} = K(F({\mathbf{x}} + {\mathbf{u}}({\mathbf{x}})), {\varepsilon _s}),$ ${N_{q}} = K(T({\mathbf{x}}),{\eta _{q}})$.

由梯度下降法, 引入时间变量 $t$, 令 ${{\mathbf{u}}^n} = {{\mathbf{u}}^n}({\mathbf{x}},t)$, 则方程 (3.3) 的稳态解可以转化成以下带有初始条件和 Dirichlet 边界条件的演化方程的解

类似地, 由梯度下降法, 可以得到 ${\alpha _s}$ 和 ${\beta _{q}}$ 的演化方程分别为

在本章的最后, 分别给出图像配准算法 1 和对应的二维微分同胚多尺度算法 2

图像配准算法 1

输入: 目标图像 $T( \cdot )$ 和待配准图像 $F( \cdot )$, 初始值 ${{\mathbf{u}}^0= {0}}$, ${{ {\alpha }}^0={c_1}(1,1,\cdots,1 )_{1 \times d}}$, ${{ {\beta }}^0}={c_2}(1,1,\cdots,1 )_{1 \times l}$, $n=1$, 初始误差 E 和最大迭代次数 K.

a) 根据初始值计算对应的方差 $Var( \cdot )$, 协方差 $Cov( \cdot, \cdot)$ 和最大相关系数 ${\rm CC}(S,D)$;

b) 由 (3.7) 式更新系数 ${{ {\alpha }}^n}$;

c) 由 ${{ {\alpha }}^n}$, ${ {\beta}}^{n-1}$ 及 ${\bf{u}}^{n-1}$ 更新 $Var(S)$, $Cov(S,\cdot)$, ${\rm CC}(S,D)$;

d) 根据新得到的 ${{ {\alpha }}^n}$ 和初始值 ${{\mathbf{u}}^{n-1}}$, ${{ {\beta }}^{n-1}}$ 及 (3.8) 式更新系数 ${{ {\beta }}^n}$;

e) 根据新得到的 ${{ {\alpha }}^n}$, ${{ {\beta }}^n}$, 初始值 ${{\mathbf{u}}^{n-1}}$ 及 (3.6) 式, 利用经典的多重网格算法^[12]--前磨光、限制、精确求解、延拓、后磨光五步来更新 $\mathbf{u}^n$;

f) $n=n+1$, 重复步骤 b) 至步骤 e), 直至达到误差限制或最大迭代次数;

输出: 配准结果 $F({\mathbf{x}} + {{\mathbf{u}}^n}({\mathbf{x}}))$.

基于上述图像配准算法 1, 将二维微分同胚多尺度多模态图像配准算法 2 表述如下

算法 2

输入: 最大尺度数 ${K_M}$.

for $scale = 1:{K_M}$

a) 令 ${\mathbf{u}} = {\mathbf{0}}$, 更新 ${\lambda _{scale}} = {\lambda _0} \times {a^{scale - 1}}$ 和 ${\varepsilon _{scale}} = {\varepsilon _0}^{scale - 1}(0 < {\varepsilon _0} < 1)$ 以满足定理 2.2 中的条件;

b) 由算法 1 得到 ${{\mathbf{u}}^{scale}}$;

c) 更新 ${{\mathbf{h}}_{^{scale}}} = {\mathbf{x}} + {{\mathbf{u}}^{scale}}$, ${{\mathbf{\tilde h}}_{^{scale}}} = {{\mathbf{\tilde h}}_{^{scale - 1}}} \circ {{\mathbf{h}}_{^{scale}}}$, 计算 ${F_{scale}}( \cdot ) = {F_{scale - 1}} \circ {{\mathbf{\tilde h}}_{^{scale}}}( \cdot )$;

end

输出: 配准结果 $F \circ {{\mathbf{\tilde h}}_{^{{K_M}}}}(\cdot)$.

本实验在 Microsoft Windows 10 操作系统下, 内存为 12.00GB, 显卡为 GeForce 920M, 处理器为 Intel(R) Core(TM) i5-6200U CPU@2.30GHz(2400 MHz) 的配置上使用 MATLAB R2018a 进行的图像配准实验.

本文为了验证算法在单模态和多模态图像配准中均有效, 主要采用以下两种数据

(1) 单模态图像数据

合成图像: A-A, A-R;

对比算法: Multi-grid algorithm of Diffeomorphic Image Registration (A4)^[12], Diffeomorphic Log Demons Image Registration algorithm (DLDIR)^[17,18], TV model^[19] 和 the Fourth order model^[20].

(2) 多模态图像数据

合成图像: A-A1, A-R1;

自然图像: Apple, Pig;

医学图像: CT-MR 肺部图像^[8], CT-synthesized 肺部图像^[21], T1-T2 脑切片图像^[22].

对比算法: the Multi-modality Non-rigid Demon (MND), Imregister and Deformable Multi-modal Image Registration (DFMIR), 其中 MND 的代码是开放的, 可以直接从文献 [23,24] 中下载, Imregister 是 MATLAB 自带的多模态图像配准函数, DFMIR 是 Shi 等^[18] 基于瑞利度量提出的多模态图像配准算法. 为了更好地与本文所提出的算法进行对比, 本文将文献 [8] 中原为多分辨率的求解器换成与本文相同的求解器--多重网格, 再进行相应的结果对比.

1) 在单模态图像配准中, 主要选取残差平方和 (Re-SSD) 和网格重叠个数 (NMF) 作为评价指标, 其表达式分别为

其中 SSD$(F,T) = \frac{1}{2}\sum\nolimits_{i,j} {{{(F(i,j) - T(i,j))}^2}}$.

其中 $\det J(\mathbf{u}) = (1 + \frac{{\partial {u_1}}}{{\partial {x_1}}})(1 + \frac{{\partial {u_2}}}{{\partial {x_2}}}) - \frac{{\partial {u_1}}}{{\partial {x_2}}}\frac{{\partial {u_2}}}{{\partial {x_1}}}$, $\# (\Psi )$ 表示集合 $\Psi$ 中元素个数.

2) 在多模态图像配准中, 除了以 NMF 作为指标来记录网格重叠个数以外, 还以 CC$\left( {F,T} \right)$ 和 Dice^[25]作为指标判断配准结果的好坏, 其表达式分别为

其中 $\left| {T \cap F({\mathbf{x}} + {\mathbf{u}}({\mathbf{x}}))} \right|$ 表示正确配准的像素点个数, $\left| T \right| + \left| {F({\mathbf{x}} + {\mathbf{u}}({\mathbf{x}}))} \right|$ 表示所有配准图像像素点的个数. Dice 值越接近于 1, 说明配准效果越好.

$\bullet$敏感性测试 本文以数据 A-A 为例, 将算法对参数 $\lambda_n$ 及参数 $\Theta$ 的敏感性进行了测试. 通过给出不同的 $\Theta$ 记录每个尺度对应的误差 Re-SSD, 并将结果绘制在图1(a) 中, 从图中可以看出, 每个 $\Theta$ 得到了几乎相同的 Re-SSD, 说明不同的参数影响了算法的收敛速度, 但不改变算法的收敛结果. 类似地, 本文还给出了算法对参数 $\lambda_n$ 的敏感性测试, 如图1(b) 所示. 数值结果表明, 当参数满足定理 2.2 的条件时, 最终配准结果可以逐步逼近 $S(\mathbf u)$ 的下确界.

$\bullet$单模态图像配准实验 为了验证所提算法在单模态图像配准中的有效性, 本文以 A-A 和 A-R 两个典型数据为例, 以 Re-SSD 和 NMF 为指标来判断所提算法在单模态图像中配准结果的好坏, 并给出了对应的配准结果示意图, 形变场 $\mathbf{h}$ 及对应的误差 Re-SSD 随尺度的变化趋势, 结果如图2-图5 所示, 随着尺度的增加, 误差在逐渐减小, 最终结果分别达到 $2.01\%$ 和 $2.74\%$, 说明配准结果 $F \circ {{\mathbf{\tilde{h}}}_{K_M}}( \cdot )$ 与目标图像 $T( \cdot )$ 较为接近. 通过与其它四种算法--A4、DLDIR、TV、Fourth order 进行对比, 突出了本文所提出的多尺度微分同胚算法能很好地避免网格重叠, 且计算速度较快.

$\bullet$多模态图像配准实验 该实验分为两部分: 有效性检验和各种多模态图像配准实验, 通过多种类型的数据来突出算法在多模态图像配准中的有效性.

算法有效性检验 在进行多模态图像配准对比实验之前, 先用一组已知形变场的图像对来进行有效性测试. 如图6 所示, $O(\cdot)$ 为原始图像, ${{\mathbf{h}}_0}({\cdot})$ 为已知的形变场, 将 ${{\mathbf{h}}_0}({\cdot})$ 作用到 $O(\cdot)$ 上, 得到一个新的图像 $F(\cdot)$, 作为该测试中的待配准图像. 进一步地, 给 $O(\cdot)$ 一定的灰度转换, 又得到一个新的图像 $T(\cdot)$, 作为该测试中的目标图像. 则 $F(\cdot)$ 和 $T(\cdot)$ 可以看作一组多模态图像对. 通过图像配准得到结果 $F_\mathbf{u}(\cdot)$, $\mathbf{h}(\cdot)$ 为新得到的形变场. 如图6(h) 所示, 原始图像与配准后的图像几乎完全重合, 进一步说明了算法的有效性.

多模态图像配准 本文分别选用了合成图像、自然图像和医学图像来研究算法的有效性. 配准结果如图7-图20 所示, 并给出了对应的 MCC 随尺度变化的趋势图. 通过与其它三种不同的算法—— MND、Imregister、DFMIR 进行比较, 结果如表2-表4 所示, 以 MCC 和 NMF 作为主要研究指标, Dice 作为辅助指标, 进一步说明了所提算法有效地避免了网格重叠, 且配准结果也优于其它三种算法.

合成图像配准

从图7-图10 中可以看出, A-A1 和 A-R1 对应的 MCC 都在随着尺度的增加逐渐增大, 配准后的图像与目标图像越来越相似, 且通过与其它三种算法配准结果的定量比较, 进一步突出了本文所提算法的优势.

自然图像配准

如图11-图14 及表3 结果所示, 通过与其它三种不同的算法结果进行对比, 发现 Imregister 有较多的网格重叠, 与物理原理相矛盾, 而算法 MND 和 DFMIR 的 MCC 和 Dice 对应的数值没有算法 2 好, 综合网格重叠个数 NMF 和 MCC 等指标, 进一步说明了算法 2 的配准效果更好.

医学图像配准

从表4 三组医学图像的配准结果中可以看出, DFMIR 算法虽然配准结果较好, 但是有明显的网格重叠, 算法 2 有效地避免了网格重叠且配准效果较好, 说明了算法 2 在多模态医学图像配准中的有效性.

本文在多模态图像配准中提出了一种基于瑞利度量的二维微分同胚多尺度算法, 有效地避免了网格重叠, 且实现了在没有先验正则项的情况下寻找二维微分多模态图像配准的最优解, 并通过具体的实验证明了算法在单模态和多模态图像配准中的有效性.