基于拟共形理论的分数阶多尺度微分同胚图像配准

图1 单元中心分割: $ \Omega_{i,j} $

设 $ \overrightarrow D = \overrightarrow D \big( {PX} \big) \in \mathbb{R}^{{n^2} \times 1} $ 为离散化的目标图像, $ \overrightarrow T = \overrightarrow T \big( {PX + PU} \big) \in \mathbb{R}^{{n^2} \times 1} $ 为离散化的目标图像. 其中 $ P \in \mathbb{R}^{2{n^2} \times 2{{(n+1)}^2}} $ 是从 $ U $ 的节点网络转移到单元中心位置的平均矩阵. 因此, 对于数据项的离散, 可以得到如下式子

(3.3)$\begin{equation} S(\boldsymbol{u}) \approx {h^2}{\big( {\overrightarrow T \big( {PX + PU} \big) - \overrightarrow D \big( {PX} \big)} \big)^T}\big( {\overrightarrow T \big( {PX + PU} \big) - \overrightarrow D \big( {PX} \big)} \big). \end{equation}$

带分数阶正则项 $ R_0(\boldsymbol{u}) $ 的离散化

接下来, 我们对带分数阶导数的正则项 $ {R_0}(\boldsymbol{u}) = \int_\Omega {{{\left| {{\nabla ^\alpha }\boldsymbol{u}} \right|}^2}{\rm d}\boldsymbol{x}} $ 进行离散化处理.

(3.4)$\begin{equation} {R_0}(\boldsymbol{u}) = \int_\Omega {{{\left| {{\nabla ^\alpha }\boldsymbol{u}} \right|}^2}{\rm d}\boldsymbol{x}} \\ = \int_\Omega {{{\big( {\frac{{{\partial ^\alpha }{u_1}}}{{\partial {x_1}^\alpha }}} \big)}^2} + } {\big( {\frac{{{\partial ^\alpha }{u_1}}}{{\partial {x_2}^\alpha }}} \big)^2} + {\big( {\frac{{{\partial ^\alpha }{u_2}}}{{\partial {x_1}^\alpha }}} \big)^2} + {\big( {\frac{{{\partial ^\alpha }{u_2}}}{{\partial {x_2}^\alpha }}} \big)^2}{\rm d}\boldsymbol{x}, \\ \end{equation}$

根据图1(b) 及数值积分中的中点法则, 可以得到

$\begin{equation*} \left\{\begin{aligned} &\int_{\Omega_{i,j}^{x_1}}{{\big( {\frac{{{\partial ^\alpha }{u_l}}}{{\partial {x_1}^\alpha }}} \big)}^2}{\rm d}\boldsymbol{x} \approx h^2{{\big( {\partial _{{x_1},\alpha }^{i + \frac{1}{2},j}{u_l}} \big)}^2} \quad 1\leqslant j\leqslant n-1, & l=1,2,\\ &\int_{\Omega_{i,j}^{x_1}}{{\big( {\frac{{{\partial ^\alpha }{u_l}}}{{\partial {x_1}^\alpha }}} \big)}^2}{\rm d}\boldsymbol{x}\approx \frac{h^2}{2}{{\big( {\partial _{{x_1},\alpha }^{i + \frac{1}{2},j}{u_l}} \big)}^2} \quad j=0,n, & l=1,2.\\ \end{aligned}\right. \end{equation*}$

同理, 对于 $ \int_{\Omega_{i,j}^{x_2}}{{\big( {\frac{{{\partial ^\alpha }{u_l}}}{{\partial {x_2}^\alpha }}} \big)}^2}{\rm d}\boldsymbol{x} $, $ l=1,2 $, 可以得到类似的结果.

另外, 对于 $ {\partial _{{x_1},\alpha }{u_l}} $ 和 $ {\partial _{{x_2},\alpha }{u_l}} $, $ l=1,2 $, 用

$\begin{equation*} \partial _{x{}_1,\alpha }^{i + \frac{1}{2},j}{u_l} \approx \frac{{\partial _{{x_{{}_1}},\alpha }^{i,j}{u_l} + \partial _{{x_{{}_1}},\alpha }^{i + 1,j}{u_l}}}{2}\quad \partial _{x{}_2,\alpha }^{i,j + \frac{1}{2}}{u_l} \approx \frac{{\partial _{{x_{{}_2}},\alpha }^{i,j}{u_l} + \partial _{{x_{{}_2}},\alpha }^{i,j + 1}{u_l}}}{2} \end{equation*}$

来估计.

因此, 可以得到

$\begin{equation*} \begin{aligned} &{R_0}(\boldsymbol{u}) = \int_\Omega {{{\left| {{\nabla ^\alpha }\boldsymbol{u}} \right|}^2}{\rm d}\boldsymbol{x}} \\ &\approx\sum\limits_{l = 1}^2 {\sum\limits_{i = 0}^{n - 1}{{\sum\limits_{j = 1}^{n - 1}}}{{h^2}{{\big( {\partial _{{x_1},\alpha }^{i + \frac{1}{2},j}{u_l}} \big)}^2}} +\sum\limits_{l = 1}^2{\sum\limits_{i = 0}^{n - 1} {\frac{{{h^2}}}{2}{{\big( {\partial _{{x_{{}_1}},\alpha }^{i + \frac{1}{2},0}{u_l}} \big)}^2}}}+ \sum\limits_{l = 1}^2{\frac{{{h^2}}}{2}{{\big( {\partial _{{x_{{}_1}},\alpha }^{i + \frac{1}{2},n}{u_l}} \big)}^2}}}\\ &+\sum\limits_{l = 1}^2 {\sum\limits_{j = 0}^{n - 1} {\sum\limits_{i = 1}^{n - 1} {{h^2}{{\big( {\partial _{{x_{{}_2}},\alpha }^{i,j + \frac{1}{2}}{u_l}} \big)}^2}} + \sum\limits_{l = 1}^2{\sum\limits_{j = 0}^{n - 1}}{\frac{{{h^2}}}{2}{{\big( {\partial _{{x_{{}_2}},\alpha }^{0,j + \frac{1}{2}}{u_l}} \big)}^2}} +\sum\limits_{l = 1}^2{\frac{{{h^2}}}{2}{{\big( {\partial _{{x_{{}_2}},\alpha }^{n,j + \frac{1}{2}}{u_l}} \big)}^2}}}}\\ &\approx \frac{1}{4}{h^2}\left[ {\sum\limits_{i = 0}^{n - 1} {\sum\limits_{j = 1}^{n - 1} {{{\big( {\partial _{{x_1},\alpha }^{i,j}{u_1}\!+\! \partial _{{x_1},\alpha }^{i + 1,j}{u_1}} \big)}^2}\! +\! \sum\limits_{i = 0}^{n - 1} {\frac{1}{2}{{\big( {\partial _{{x_{{}_1}},\alpha }^{i,0}{u_1}\! +\! \partial _{{x_{{}_1}},\alpha }^{i + 1,0}{u_1}} \big)}^2}\! +\! \frac{1}{2}{{\big( {\partial _{{x_{{}_1}},\alpha }^{i,n}{u_1}\! +\! \partial _{{x_{{}_1}},\alpha }^{i + 1,n}{u_1}} \big)}^2}} } } } \right] \\ &+ \frac{1}{4}{h^2}\left[ {\sum\limits_{i = 0}^{n - 1} {\sum\limits_{j = 1}^{n - 1} {{{\big( {\partial _{{x_1},\alpha }^{i,j}{u_2}\! +\! \partial _{{x_1},\alpha }^{i + 1,j}{u_2}} \big)}^2}\! +\! \sum\limits_{i = 0}^{n - 1} {\frac{1}{2}{{\big( {\partial _{{x_{{}_1}},\alpha }^{i,0}{u_2}\! +\! \partial _{{x_{{}_1}},\alpha }^{i + 1,0}{u_2}} \big)}^2}\! +\! \frac{1}{2}{{\big( {\partial _{{x_{{}_1}},\alpha }^{i,n}{u_2}\! +\! \partial _{{x_{{}_1}},\alpha }^{i + 1,n}{u_2}} \big)}^2}} } } } \right] \\ &+\frac{1}{4}{h^2}{\left[ {\sum\limits_{j = 0}^{n - 1} {\sum\limits_{i = 1}^{n - 1} {{{\big( {\partial _{{x_{{}_2}},\alpha }^{i,j}{u_{\text{1}}}\! +\! \partial _{{x_{{}_2}},\alpha }^{i,j + 1}{u_{\text{1}}}} \big)}^2}\! +\! \sum\limits_{j = 0}^{n - 1} {\frac{1}{2}{{\big( {\partial _{{x_{{}_2}},\alpha }^{0,j}{u_{\text{1}}}\! +\! \partial _{{x_{{}_2}},\alpha }^{0,j + 1}{u_{\text{1}}}} \big)}^2}\! +\! \frac{1}{2}{{\big( {\partial _{{x_2},\alpha }^{n,j}{u_{\text{1}}}\! +\! \partial _{{x_2},\alpha }^{n,j + 1}{u_{\text{1}}}} \big)}^2}} } } } \right]} \\ &+\frac{1}{4}{h^2}\left[ {\sum\limits_{j = 0}^{n - 1} {\sum\limits_{i = 1}^{n - 1} {{{\big( {\partial _{{x_{{}_2}},\alpha }^{i,j}{u_{\text{2}}}\! +\! \partial _{{x_{{}_2}},\alpha }^{i,j + 1}{u_{\text{2}}}} \big)}^2}\! +\! \sum\limits_{j = 0}^{n - 1} {\frac{1}{2}{{\big( {\partial _{{x_{{}_2}},\alpha }^{0,j}{u_{\text{2}}}\! +\! \partial _{{x_{{}_2}},\alpha }^{0,j + 1}{u_{\text{2}}}} \big)}^2}\! +\! \frac{1}{2}{{\big( {\partial _{{x_{{}_2}},\alpha }^{n,j}{u_{\text{2}}}\! +\! \partial _{{x_{{}_2}},\alpha }^{n,j + 1}{u_{\text{2}}}} \big)}^2}} } } } \right] \\ &=\frac{{{h^2}}}{4}{{\big( {{\partial _\alpha }U} \big)^T}{A^T}GA\big( {{\partial _\alpha }U} \big)}, \\ \end{aligned} \end{equation*}$

其中有关矩阵 $ A $, $ G $ 的详细信息, 请参见附录 A.

根据文献 [30,31] 中的 Grunwald 近似法, 对于 $ \frac{{{\partial ^\alpha }{u_l}(x)}}{{\partial {x_p}^\alpha }}(p = 1,2;l = 1,2) $ 可以按以下方式离散化

(3.5)$\begin{equation} \frac{{{\partial ^\alpha }{u_l}\big( {{x_{i,j}}} \big)}}{{\partial {x_p}^\alpha }} = \delta _{p - }^\alpha {u_l}\big( {{x_{i,j}}} \big) + O(h), \end{equation}$

其中 $ \delta _{1 - }^\alpha {u_l}\big( {{x^{_{i,j}}}} \big) = \frac{1}{{{h^\alpha }}}\sum\limits_{q = 1}^i {\rho _q^{(\alpha)}{u_l}\big( {{x^{_{i - q + 1,j}}}} \big)} $, $ \delta _{2 - }^\alpha {u_l}\big( {{x^{_{i,j}}}} \big) = \frac{1}{{{h^\alpha }}}\sum\limits_{q = 1}^{j + 1} {\rho _q^{(\alpha)}{u_l}\big( {{x^{_{i,j - q + 1}}}} \big)} $, $ \rho _0^{(\alpha)} = 1 $, $ \rho _q^{(\alpha)} = \big( {1 - \frac{{q + \alpha }}{q}} \big)\rho _{q - 1}^{(\alpha)} $.

记 $ {U_j} = {(u_{_1}^{0,j},u_1^{1,j},\cdots,u_{_1}^{n,j})^T} \in {\mathbb{R}^{({n + 1}) \times 1}},u_{_l}^{0,j} \triangleq {u_l}\big( {{x^{0,j}}} \big),j = 0,1, \cdots, n $, 根据 (3.5) 式则有

(3.6)$\begin{equation} \frac{{{\partial ^\alpha }{U_j}}}{{\partial {x_1}^\alpha }} = {\Big( {\frac{{{\partial ^\alpha }u_1^{0,j}}}{{{{\partial {x_1}^\alpha }} }},\frac{{{\partial ^\alpha }u_1^{1,j}}}{{{{\partial {x_1}^\alpha }} }},\cdots,\frac{{{\partial ^\alpha }u_{_1}^{n,j}}}{{{{\partial {x_1}^\alpha }} }}} \Big)^T} \approx {B_{n + 1,\alpha }^h} {U_j}, \end{equation}$

其中

(3.7)$\begin{equation} {B_{n + 1,\alpha }^h} =\frac{1}{{{h^\alpha }}}{B_{n + 1,\alpha }}=\frac{1}{{{h^\alpha }}}\left( {\begin{array}{*{20}{c}} \begin{gathered} \rho _1^{(\alpha)} \hfill \\ \rho _2^{(\alpha)} \hfill \\ \end{gathered} &\begin{gathered} \rho _0^{(\alpha)} \hfill \\ \rho _1^{(\alpha)} \hfill \\ \end{gathered} &\begin{gathered} 0 \hfill \\ \rho _0^{(\alpha)} \hfill \\ \end{gathered} \\ \vdots & \vdots & \ddots \\ \begin{gathered} \rho _n^{(\alpha)} \hfill \\ \rho _{n + 1}^{(\alpha)} \hfill \\ \end{gathered} &\begin{gathered} \rho _{n - 1}^{(\alpha)} \hfill \\ \rho _n^{(\alpha)} \hfill \\ \end{gathered} &\begin{gathered} \cdots \hfill \\ \cdots \hfill \\ \end{gathered} \end{array}\begin{array}{*{20}{c}} \begin{gathered} \cdots \hfill \\ \cdots \hfill \\ \ddots \hfill \\ \rho _1^{(\alpha)} \hfill \\ \end{gathered} &\begin{gathered} 0 \hfill \\ 0 \hfill \\ \vdots \hfill \\ \rho _0^{(\alpha)} \hfill \\ \end{gathered} \\ {\rho _2^{(\alpha)}}&{\rho _1^{(\alpha)}} \end{array}} \right) \in {\mathbb{R}^{(n + 1) \times (n + 1)}}. \end{equation}$

由以上分数阶离散理论可知 $ {\partial _\alpha }U = \frac{1}{h^{\alpha}}CU $. 有关矩阵 $ C $ 的详细信息,请参见附录 A.

因此, 对带分数阶导数的正则项 $ {R_0}\big( u \big) $ 的离散, 我们可以得出以下结论

(3.8)$\begin{equation} {R_0}\big( u \big) \approx \frac{{{h^{2-2\alpha}}}}{4}{U^T}{C^T}{A^T}GACU. \end{equation}$

正则项 $ R_1(\boldsymbol{u}) $ 的离散化

最后, 是对正则项 $ {R_1}\big( \boldsymbol{u} \big) $ 的离散. 根据文献 [11] 可以得到

(3.9)$\begin{equation} {R_1}\big( \boldsymbol{u} \big) \approx \frac{{{h^2}}}{4}\phi \big( {\vec r\big( U \big)} \big){{\rm e}^T}, \end{equation}$

其中 $ \phi \big( {\vec r\big( U \big)} \big) = \big( {\phi \big( {\vec r{{\big( U \big)}_1}} \big), \cdots, \phi \big( {\vec r{{\big( U \big)}_{4{n^2}}}} \big)} \big) $, $ e = \big( {1,1, \cdots 1} \big) \in {\mathbb{R}}^{4{n^2}} $. $ \vec r\big( U \big) $ 为离散化的形变场 $ \boldsymbol{h} $ 的 Beltrami 系数的模的平方, 详见文献 [11] 中的附录 B.

最终, 结合 (3.7), (3.8) 和 (3.9) 式, 我们可以得到模型 (3.1) 的离散化形式如下

(3.10)$\begin{equation} \begin{aligned} {U^n} &\triangleq \arg \mathop {\inf }\limits_{\boldsymbol{u} \in {{\left[ {H_0^\alpha ( \Omega )} \right]}^2}} J\big( U \big) \\ &= \arg \mathop {\inf }\limits_{\boldsymbol{u} \in {{\left[ {H_0^\alpha ( \Omega )} \right]}^2}} {\lambda _n}{h^2}{\big( {\overrightarrow T \big( {PX + PU} \big) - \overrightarrow D \big( {PX} \big)} \big)^T}\big( {\overrightarrow T \big( {PX + PU} \big) - \overrightarrow D \big( {PX} \big)} \big) \\ & + {\gamma _n}\frac{{{h^{2-2\alpha}}}}{4}{U^T}{C^T}{A^T}GACU + {\Theta _n}\frac{{{h^2}}}{4}\phi \big( {\vec r\big( U \big)} \big){{\rm e}^T}. \end{aligned} \end{equation}$

3.2.2 离散化问题 (3.10) 的最优化方法

在数值计算中, 我们选择线性搜索法求解无约束优化问题 (3.10). 为了保证搜索方向是一个下降方向, 我们采用 Gauss-Newton 方向. 另外, 使用 Gauss-Newton 法时, 我们不需要计算二阶项, 可以加速 CPU 时间消耗.

$ J( U ):{\mathbb{R}}^{2{{( {n + 1} )}^2}} \to \mathbb{R} $ 是两次连续可微分的, $ {U^i} \in {{\mathbb{R}}^{2{{( {n + 1} )}^2}}} $, 且近似 Hessian 矩阵 $ H( {{U^i}} ) $ 是正定的. 用二次逼近的方法将 $ J( U) $ 在 $ {U^i} $ 处展开

(3.11)$\begin{equation} J\big( {{U^i} + s} \big) \approx {q^i}(s) = J\big( {{U^i}} \big) + {d_J}{\big( {{U^i}} \big)^T}s + \frac{1}{2}{s^T}H{\big( {{U^i}} \big)^T}s, \end{equation}$

其中 $ s = U - {U^i} $, $ {d_J}\big( {{U^i}} \big) = \nabla J\big( {{U^i}} \big) $. 同时, 由 (3.11) 式可以得到

(3.12)$\begin{equation} {U^{i + 1}} = {U^i} - {\left[ {H\big( {{U^i}} \big)} \right]^{ - 1}}{d_J}\big( {{U^i}} \big). \end{equation}$

为了保证 Gauss-Newton 法的全局收敛性, 我们采用了线性搜索法. 其迭代过程如下

(3.13)$\begin{equation} {U^{i + 1}} = {U^i} - {\theta _i}{\left[ {H\big( {{U^i}} \big)} \right]^{ - 1}}{d_J}\big( {{U^i}} \big), \end{equation}$

其中 $ {\theta _i} $ 表示步长.

接下来, 我们将研究近似的 Hessian 矩阵 $ H\big( {{U^i}} \big) $、搜索方向, 步长及算法终止标准.

Hessian 矩阵 $ H $ 估计

我们分别考虑 (3.10) 式中 $J (U)$ 的三个项.

首先, 是离散化的 SSD 项

(3.14)$\begin{equation} {h^2}{\big( {\overrightarrow T \big( {PX + PU} \big) - \overrightarrow D \big( {PX} \big)} \big)^T}\big( {\overrightarrow T \big( {PX + PU} \big) - \overrightarrow D \big( {PX} \big)} \big), \end{equation}$

其梯度和 Hessian 矩阵分别为

(3.15)$\begin{equation} \left\{ \begin{aligned} &{d_1} = 2{h^2}{P^T}\overrightarrow T _{\widetilde U}^T{\big( {\overrightarrow T \big( {\widetilde U} \big) - \overrightarrow D } \big)^T} \in {\mathbb{R}^{2{{(n+1)}^2} \times 1}}, \\ &{H_1} = 2{h^2}{P^T}\big( {\overrightarrow T _{\widetilde U}^T{{\overrightarrow T }_{\widetilde U}} + \sum\limits_{l = 1}^{{n^2}} {{{\big( {\overrightarrow T \big( {\widetilde U} \big) - \overrightarrow D } \big)}_l}{\nabla ^2}{{\big( {\overrightarrow T \big( {\widetilde U} \big) - \overrightarrow D } \big)}_l}} } \big)P, \end{aligned} \right. \end{equation}$

其中 $ \widetilde {U} = PX + PU $, $ \overrightarrow T _{\widetilde {U}} = \frac{{\partial \overrightarrow T ( \widetilde {U} )}}{{\partial \widetilde {U}}} $.

对于 $ {H_1} $, 我们无法确保它是半正定的. 如果它不是正定的, 我们可能得不到下降方向. 因此, 我们考虑省略 $ {H_1} $ 的二阶项, 得到 (3.14) 式的近似 Hessian 矩阵

(3.16)$\begin{equation} \widehat {{H_1}} = 2{h^2}{P^T} \big( {\overrightarrow T _{\widetilde {U}}^T}{{\overrightarrow T }_{\widetilde {U}}} \big)P. \end{equation}$

另外, 对于离散带分数阶导数的正则项 $ \frac{{{h^{2-2\alpha}}}}{4}{U^T}{C^T}{A^T}GACU $, 其梯度和 Hessian 矩阵分别为

(3.17)$\begin{equation} \left\{ \begin{aligned} &{d_2} = \frac{h^{2-2\alpha}}{2}{C^T}{A^T}GACU \in \mathbb{R}^{2{{(n+1)}^2} \times 1}, \\ &{H_2} =\frac{h^{2-2\alpha}}{2}{C^T}{A^T}GAC \in {\mathbb{R}}^{2{{(n+1)}^2} \times 2{{(n+1)}^2}}. \end{aligned} \right. \end{equation}$

最后, 对于 $ J\big( U \big) $ 的第三项 $ \frac{{{h^2}}}{4}\phi \big( {\vec r\big( U \big)} \big){{\rm e}^T} $, 其梯度和 Hessian 矩阵 $ H $ 分别为

(3.18)$\begin{equation} \left\{ \begin{aligned} &{d_3} = \frac{h^2}{4}{\rm d}{{\vec r}^T}{\rm d}\phi \big( {\vec r} \big) \in {\mathbb{R}}^{2{{(n+1)}^2} \times 1}, \\ &{H_3} =\frac{h^2}{4}\big( {{\rm d}{{\vec r}^T}{d^2}\phi \big( {\vec r} \big){\rm d}\vec r + \sum\limits_{l = 1}^{4{n^2}} {{{\left[ {{\rm d}\phi \big( {\vec r} \big)} \right]}_l}{\nabla ^2}{{\vec r}_l}} } \big) \in {\mathbb{R}^{2{{(n+1)}^2} \times 2{{(n+1)}^2}}}. \end{aligned} \right. \end{equation}$

为了从 (3.18) 式中提取半正定部分, 我们省略其二阶项, 得到近似的 Hessian 矩阵为

(3.19)$\begin{equation} \widehat {{H_3}} = \frac{h^2}{4}{h^2}{\rm d}{\vec r^T}{d^2}\phi \big( {\vec r} \big){\rm d}\vec r. \end{equation}$

因此, 对于 (3.10) 式中的函数 $ J\big( U \big) $, 我们可以得到它的梯度矩阵为

(3.20)$\begin{equation} {d_J}^n = \lambda^n{d_1} + \gamma^n{d_2} + \Theta^n{d_3}, \end{equation}$

以及近似的 Hessian 矩阵

(3.21)$\begin{equation} H^n = \lambda^n \widehat {{H_1}} + \gamma^{n}{H_2} + \Theta^n\widehat {{H_3}}. \end{equation}$

搜索方向

每次迭代时, 利用式 (3.20) 和 (3.21) 式, 我们需要求解高斯-牛顿系统, 找到 (3.10) 式的搜索方向

(3.22)$\begin{equation} H\delta U = - {d_J}, \end{equation}$

其中 $ \delta U $ 是搜索方向. 在我们的实验过程中, 我们使用带有对角线预处理的极小化残差算法 (MINRES)^[32,33]来求解这个线性系统.

步长和终止标准

我们使用标准化的 Armijo 线搜索策略, 通过反向追踪来找到合适的步长 $ \theta $. 在算法过程中, 我们还需要确保 $ \vec r\big( U \big) $ 小于 1. 我们分别选择 $ T0=10^{-8} $ 和 $ Tol=10^{-12} $ 作为步长 $ \theta $ 和 $ \theta\left\| \delta U \right\| $ 的下限^{[34⇓⇓-37]}. 因此给出以下算法 1

算法 1

输入: 初始化 ID=0, $ \theta=1 $, $ T0=10^{-8} $, $ Tol=10^{-12} $, $ \eta=10^{-4} $, $ U $, $ \delta U $;

a) 计算 $ J(U) $ 和 $ d_J $;

$\text { b) while } \theta\|\delta U\| \geqslant T o l \text {, 则 } U^{n e w}=U+\theta\|\delta U\| \text {; }$

$if \vec{r}\left(U^{n e w}\right) \leqslant 1, 则 $

$if \theta \geqslant T 0, 退出 while 循环, 执行 c );$

$else \ if \theta＜T 0, 执行 d );$

$结束 if 判断;$

更新 $ \theta=\frac{\theta}{2} $;

结束 while 循环;

c) while $ \theta\left\| \delta U \right\|\geqslant Tol $,则计算 $ J(U^{new}) $;

if $ J(U^{new})\leqslant J(U)+\theta\eta {d_J}^T\delta U $, 则

if $ \theta\geqslant T0 $, 结束算法, 输出 $ U=U^{new} $;

else if $ \theta < T0 $, 执行 d);

结束 if 判断;

更新 $ \theta=\frac{\theta}{2} $, $ U^{new}=U+\theta\delta U $;

结束 while 循环;

d) ID=1, $ U=U^{new} $;

输出: ID, $ U $.

根据文献 [38], 给定

$\begin{equation*} \begin{aligned} &{\rm (i.a)} \left\| J(U^{i+1})-J(U^{i}) \right\|\leqslant\tau_J(1+\left\| J(U^0) \right\|), \\ &{\rm (i.b)} \left\| {\boldsymbol{h}}^{i+1}-{\boldsymbol{h}}^{i}\right\|\leqslant\tau_W\left\| (1+{\boldsymbol{h}}^0) \right\|,\\ &{\rm (i.c)} \left\| d_J \right\|\leqslant\tau_G(1+\left\| J(U^0) \right\|), \\ &{\rm (ii)} d_J\leqslant {\rm eps},\\ &{\rm (iii)} i\geqslant {\rm MaxIter}. \\ \end{aligned} \end{equation*}$

其中 eps 是机器精度, MaxIter 是外部迭代的最大次数. 我们设置 $ \tau_J=10^{-3} $, $ \tau_W=10^{-2} $, $ \tau_G=10^{-2} $, MaxIter=500. 如果 (i)(ii) 或 (iii) 中的任何一项得到满足, 迭代就会终止. 因此利用 Armijo 搜索得到的 Gauss-Newton 算法可由下述算法 2 给出

算法 2

输入: 目标图像 $ D(\cdot) $, 待配准图像 $ T(\cdot) $, $ U^0 $, $ \Theta $, $ \gamma $, $ \alpha $;

a) $ i=0 $, $ U^i=U^0 $;

b) 根据 (3.10) 式计算 $ J(U^i) $, 根据 (3.20) 和 (3.21) 式计算 $ {d_J}^i $ 和 $ H^i $;

c) while (i)--(iii) 均不满足, 则根据 $ H^i\delta U^i=-{d_J}^i $ 得到 $ U^i $ 和 $ \delta U^i $, 执行算法 1, 得到 $ U^{i+1} $ 和 ID;

if ID=1, 则结束该算法, 输出 $ U=U^{i+1} $ 和 ID;

else $ i=i+1 $, 计算 $ J(U^i) $, $ {d_J}^i $ 和 $ H^i $;

结束 if 循环;

结束 while 循环;

输出: ID, $ U $.

3.2.3 多尺度图像配准算法

在本章的最后, 给出对应的二维多尺度微分同胚算法 3

算法 3

输入: 最大尺度数 $ K_M $, $ \lambda_{0} $, $ \gamma_0 $, $ \Theta_0 $, $ \alpha $;

for $ scale = 1:K_M $

a) 令 $ \boldsymbol{u}=\boldsymbol{0} $, 更新 $ \lambda_{scale}= \lambda_0\times a^{scale-1} $, $ \gamma=\frac{\gamma_0}{\lambda_{scale}} $, $ \Theta=\frac{\Theta_0}{\lambda_{scale}} $;

b) 由算法 2 得到 $ \boldsymbol{u}^{scale} $;

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

end

输出: 配准结果 $ T\circ\overline{{\boldsymbol{h}}_{K_M}}(\cdot) $.

4 数值实验

4.1 实验数据及对比算法

为验证算法的有效性, 本文主要采用以下几类图像数据及对比算法

实验数据

(i) 合成图像数据: A-A, A-R;

(ii) 自然图像数据: Watermelon, Pineapple-Pepper;

(iii) 医学图像数据: Liver MRI, Brain MRI, Hand CT.

对比算法: TV model^[39], Demons, the Multimodality Nonrigid Demon (MND)^[40,41].

4.2 评价指标

在本文中, 主要选取相对残差平方和 ($ Re-SSD $) 和网格重叠个数 (MFN) 两个指标来判断配准结果的好坏, 其表达式分别为

(4.1)$\begin{equation} Re-SSD=\frac{SSD(T(\boldsymbol{x}+\boldsymbol{u}(\boldsymbol{x})),D)}{SSD(T,D)}, \end{equation}$

其中 $ SSD(T,D)=\frac{1}{2}\sum_{i,j}(T_{i,j}-D_{i,j})^2 $;

(4.2)$\begin{equation} MFN(\boldsymbol{u})=\#(\det J(\boldsymbol{u})\leqslant0), \end{equation}$

其中 $ \det(J(\boldsymbol{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}}} $, $ \#(A) $ 表示集合 $ A $ 中元素个数.

4.3 数值实验

4.3.1 合成图像配准

为了验证所提算法在图像配准中的有效性, 在此部分以 A-A 和 A-R 两对典型合成数据为例, 以 $ Re-SSD $ 和 $ MFN $ 为指标来判断所提算法在合成图像中配准结果的好坏, 并给出了对应的配准结果示意图, 形变场 $ \boldsymbol{h} $ 及对应的 $ Re-SSD $ 随尺度的变化趋势, 结果如图2-5 所示. 从图中可以看出, 随着尺度的增加, 误差 $ Re-SSD $ 在逐渐减小, 最终误差分别为 3.86 $ \% $ 和 6.86 $ \% $, 这说明配准结果 $ T\circ{\overline{\boldsymbol{h}_K}}{(\cdot)} $ 与目标图像 $ D(\cdot) $ 较为接近. 另外, 如表1 所示, 通过与其它三种算法进行对比, 结果表明, 所提出的算法误差较小并且能够有效地避免网格重叠.

图2

图2 A-A配准结果

图3

图3 A-A 在不同尺度下的配准结果

图4

图4 A-R 在不同尺度下的配准结果

图5

图5 A-R在不同尺度下的配准结果

表1 合成图像配准中四种不同配准算法的配准结果定量比较

新窗口打开| 下载CSV

4.3.2 自然图像配准

在本节中, 我们对两个不同的自然图像对“西瓜、菠萝-辣椒”进行了配准.

使用所提出的算法对其进行配准, 定量评估指标和最终配准结果见图6-9. 为了体现所提算法的竞争力, 还采用了另外三种算法同时对上述自然图像进行配准, 并将配准结果与所提算法的最终结果进行比较; 比较结果如表2 所示. 表2 显示, 与其他三种算法相比, 所提出的算法具有很强的竞争力.

图6

图6 Watermelon 配准结果

图7

图7 Watermelon 在不同尺度下的配准结果

图8

图8 Pineapple 配准结果

图9

图9 Pineapple 在不同尺度下的配准结果

表2 自然图像配准中四种不同配准算法的配准结果定量比较

新窗口打开| 下载CSV

4.3.3 医学图像配准

在这个部分, 我们对三个不同的医学图像对“肝脏、大脑和手部”进行了配准.

使用所提算法对上述医学图像对进行配准, 定量评价指标及最终配准结果如图10-15 所示. 另外, 为体现所提算法的竞争性, 同时用其他三种算法对上述医学图像进行配准, 将得到的配准结果与所提算法的配准结果进行对比, 比较结果如表3 所示. 由表3 可知, 本文所提算法与其他 3 种算法相比, 具有较强的竞争力.

图10

图10 Liver 配准结果

图11

图11 Liver 在不同尺度下的配准结果

图12

图12 Brain 配准结果

图13

图13 Brain 在不同尺度下的配准结果

图14

图14 Hand 配准结果

图15

图15 Hand 在不同尺度下的配准结果

表3 医学图像配准中四种不同配准算法的配准结果定量比较

新窗口打开| 下载CSV

注 4.1 TV 模型是算法 3 中 $ \alpha=1 $ 的特例. 从表1-3 的对比可以看出, $ \alpha $ 取整数虽然能加速算法的 CPU 时间消耗, 但计算精度要低很多. 本文更关心于精度方法的提升, 这也是为什么引入分数阶正则项 (非局部算子) 的重要原因.

5 结论

本文提出了一种基于拟共形理论的分数阶微分同胚多尺度方法, 该方法有效地避免了网格重叠, 同时在没有先验正则项的情况下找到了微分同胚图像配准的最优解, 并且通过数值实验验证了所提算法的有效性.

附录 A 3.2.1 节中矩阵 $ A, G, C $ 的计算

令矩阵 $ {A_1}={I_2} \otimes {I_{n + 1}} \otimes \partial _n^1 \in {R^{2n(n+1) \times 2{{(n+1)}^2}}} $, $ {A_2} = {I_2} \otimes \partial _n^1 \otimes {I_{n + 1}} \in {R^{2n(n+1) \times 2{{(n+1)}^2}}} $,

其中 $ \partial _n^1 = \left[ {\begin{array}{*{20}{c}} 1&1&{}&{}&{} \\ {}&1&1&{}&{} \\ {}& \cdots & \cdots & \cdots &{} \\ {}&{}&1&1&{} \\ {}&{}&{}&1&1 \end{array}} \right] \in {R^{n,n + 1}} $, $ \otimes $ 表示克罗内克积,则

$A = \left[ {\begin{array}{*{20}{c}} {A_1}&{{0_{2n(n+1) \times 2{{(n+1)}^2}}}} \\ {{0_{2n(n+1) \times 2{{(n+1)}^2}}}}&{A_2} \end{array}} \right] \in {R^{4n(n+1) \times 4{{(n+1)}^2}}}.$

设矩阵 $ G_1 $ 满足 $ {G_{{1_{i + 1 + jn,i + 1 + jn}}}} = \left\{ \begin{gathered} 1,{\text{ }}0 \leqslant i \leqslant n - 1,{\text{ }}1 \leqslant j \leqslant n - 1 \hfill \\ \frac{1}{2},{\text{ }}0 \leqslant i \leqslant n - 1,{\text{ }}j = 0,n \hfill \\ \end{gathered} \right.,{\text{ }}{G_{{1_{i,j}}}} = 0,{\text{ }}i \ne j $, 矩阵 $ G_2 $ 满足 $ {G_{{2_{i + 1 + j(n+1),i + 1 + j(n+1)}}}} = \left\{ \begin{gathered} 1,{\text{ }}1 \leqslant i \leqslant n - 1,{\text{ }}0 \leqslant j \leqslant n - 1 \hfill \\ \frac{1}{2},{\text{ }}i = 0,n,{\text{ }}0 \leqslant j \leqslant n - 1 \hfill \\ \end{gathered} \right.,{\text{ }}{G_{{2_{i,j}}}} = 0,{\text{ }}i \ne j $, 则

$G = \left[ {\begin{array}{*{20}{c}} {{G_1}}&0&0&0 \\ 0&{{G_1}}&0&0 \\ 0&0&{{G_2}}&0 \\ 0&0&0&{{G_2}} \end{array}} \right] \in {R^{4n(n+1) \times 4n(n+1)}}.$

令矩阵 $ C_1={{I_2} \otimes {I_{n + 1}} \otimes {B_{n + 1,\alpha }}}\in {R^{ 2(n+1)^2\times 2(n+1)^2}} $, $ {C_2} = \left[ {{I_2} \otimes {B_{n + 1,\alpha }} \otimes {I_{n + 1}}} \right](list,:) $ $ \in {R^{ 2(n+1)^2\times 2(n+1)^2}} $,其中 $ list = \big( 1 + 0(n+1), \cdots, 1 + n(n+1), \cdots, n + 1 + 0(n+1), \cdots, n + 1 + n(n+1) \big)^T\in{R^{(n+1)^2\times 1}} $,则

$C=\Bigg[ {\begin{array}{*{20}{c}} {{C_1}} \\ {{C_2}} \end{array}} \Bigg]\in {R^{4(n+1)^2\times 2(n+1)^2}}.$

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Brown

L G

A survey of image registration techniques

ACM Computing Surveys, 1992, 24(4): 325-376

[2]

Lester

, Arridge

S R

A survey of hierarchical non-linear medical image registration

Pattern Recognition, 1998, 32(1): 129-149

[3]

Maintz

J B

, Viergever

M A

A survey of medical image registration

Medical Image Analysis, 1998, 2(1): 1-36

DOI:10.1016/s1361-8415(01)80026-8 PMID:10638851 [本文引用: 1]

The purpose of this paper is to present a survey of recent (published in 1993 or later) publications concerning medical image registration techniques. These publications will be classified according to a model based on nine salient criteria, the main dichotomy of which is extrinsic versus intrinsic methods. The statistics of the classification show definite trends in the evolving registration techniques, which will be discussed. At this moment, the bulk of interesting intrinsic methods is based on either segmented points or surfaces, or on techniques endeavouring to use the full information content of the images involved.

[4]

Mohamed

, Zacharaki

E I

, Shen

et al.

Deformable registration of brain tumor images via a statistical model of tumor-induced deformation

Medical Image Analysis, 2006, 10(5): 752-763

PMID:16860588 [本文引用: 1]

An approach to the deformable registration of three-dimensional brain tumor images to a normal brain atlas is presented. The approach involves the integration of three components: a biomechanical model of tumor mass-effect, a statistical approach to estimate the model's parameters, and a deformable image registration method. Statistical properties of the sought deformation map from the atlas to the image of a tumor patient are first obtained through tumor mass-effect simulations on normal brain images. This map is decomposed into the sum of two components in orthogonal subspaces, one representing inter-individual differences in brain shape, and the other representing tumor-induced deformation. For a new tumor case, a partial observation of the sought deformation map is obtained via deformable image registration and is decomposed into the aforementioned spaces in order to estimate the mass-effect model parameters. Using this estimate, a simulation of tumor mass-effect is performed on the atlas image in order to generate an image that is similar to tumor patient's image, thereby facilitating the atlas registration process. Results for a real tumor case and a number of simulated tumor cases indicate significant reduction in the registration error due to the presented approach as compared to the direct use of deformable image registration.

[5]

Sotiras

, Davatzikos

, Paragios

Deformable medical image registration: a survey

IEEE Transactions on Medical Imaging, 2013, 32(7): 1153-1190

DOI:10.1109/TMI.2013.2265603 PMID:23739795 [本文引用: 1]

Deformable image registration is a fundamental task in medical image processing. Among its most important applications, one may cite: 1) multi-modality fusion, where information acquired by different imaging devices or protocols is fused to facilitate diagnosis and treatment planning; 2) longitudinal studies, where temporal structural or anatomical changes are investigated; and 3) population modeling and statistical atlases used to study normal anatomical variability. In this paper, we attempt to give an overview of deformable registration methods, putting emphasis on the most recent advances in the domain. Additional emphasis has been given to techniques applied to medical images. In order to study image registration methods in depth, their main components are identified and studied independently. The most recent techniques are presented in a systematic fashion. The contribution of this paper is to provide an extensive account of registration techniques in a systematic manner.

[6]

Zitova

, Flusser

Image registration methods: a survey

Image and Vision Computing, 2003, 21(11): 977-1000

[7]

Budhiraja

, Dupuis

, Maroulas

Large deviations for stochastic flows of diffeomorphisms

Bernoulli, 2010, 16(1): 234-257

[8]

Chen

Y M

, Shi

J L

, Rao

et al.

Deformable multi-modal image registration by maximizing rényi's statistical dependence measure

Inverse Problems and Imaging, 2015, 9(1): Article 79

[9]

Chen

Y M

, Ye

X J

Inverse consistent deformable image registration

The Legacy of Auadi Ramakrishnan in the Mathematical Sciences, 2010: 419-440

[10]

Yang

X M

, Shen

C M

, Li

et al.

A combination of the total variation filter and a fourth-order filter for image registration

Mathematical Problems in Engineering, 2015, 1: 235134

[11]

Zhang

D P

, Chen

A novel diffeomorphic model for image registration and its algorithm

Journal of Mathematical Imaging and Vision, 2018, 60: 1261-1283

[本文引用: 4]

[12]

Zhang

J P

, Chen

Variational image registration by a total fractional-order variation model

Journal of Computational Physics, 2015, 293: 442-461

[本文引用: 4]

[13]

Zhang

, Chen

, Yu

An improved discontinuity-preserving image registration model and its fast algorithm

Applied Mathematical Modelling, 2016, 40(23/24): 10740-10759

[14]

Zhang

, Chen

, Yu

A novel high-order functional based image registration model with inequality constraint

Computers and Mathematics with Applications, 2016, 72(12): 2887-2899

[15]

Machado

J T

, Kiryakova

, Mainardi

Recent history of fractional calculus

Communications in Nonlinear Science and Numerical Simulation, 2011, 16(3): 1140-1153

[16]

Abdeljawad

On conformable fractional calculus

Journal of Computational and Applied Mathematics, 2015, 279: 57-66

[17]

Han

A variational model with fractional-order regularization term arising in registration of diffusion tensor image

Inverse Problems and Imaging, 2018, 12(6): 1263-1291

[本文引用: 3]

[18]

Metzler

, Klafter

The random walk's guide to anomalous diffusion: a fractional dynamics approach

Physics Reports, 2000, 339(1): 1-77

[19]

Han

, Wang

Z Q

An alternating direction implicit scheme of a fractional-order diffusion tensor image registration model

Applied Mathematics and Computation, 2019, 356: 105-118

[20]

Haber

, Modersitzki

Numerical methods for volume preserving image registration

Inverse Problems, 2004, 20(5): 1621-1638

[21]

Lui

L M

, Thiruvenkadam

, Wang

Y L

et al.

Optimized conformal surface registration with shape-based landmark matching

SIAM Journal on Imaging Sciences, 2010, 3(1): 52-78

[22]

Lui

L M

, Wong

T W

, Zeng

et al.

Optimization of surface registrations using Beltrami holomorphic flow

Journal of Scientific Computing, 2012, 50(3): 557-585

[23]

Lam

K C

, Lui

L M

Landmark and intensity-based registration with large deformations via quasi-conformal maps

SIAM Journal on Imaging Sciences, 2014, 7(4): 2364-2392

[24]

Han

, Wang

Z P

, Zhang

Y M

Multi-scale approach for two-dimensional diffeomorphic image registration

Multiscale Modeling and Simulation, 2021, 19(4): 1538-1572

[本文引用: 3]

[25]

Lehto

, Virtanen

K I

. Quasiconformal Mappings in the Plane. New York: Springer, 1973

[26]

Gardiner

F P

, Lakic

Quasiconformal Teichmuller Theory

Providence RI: American Mathematical Society, 2000

[27]

Ahlfors

L V

Lectures on Quasiconformal Mappings

Providence RI: American Mathematical Society, 2006

[28]

Aubert

, Kornprobst

, Aubert

. Mathematical Problems in Image Processing:Partial Differential Equations and the Calculus of Variations. New York: Springer, 2006

[29]

Han

, Wang

Z P

A diffeomorphic image registration model with fractional-order regularization and Cauchy-Riemann constraint

SIAM Journal on Imaging Sciences, 2020, 13(3): 1240-1271

[30]

Hou

T L

, Tang

, Yang

Numerical analysis of fully discretized Crank-Nicolson scheme for fractional-in-space Allen-Cahn equations

Journal of Scientific Computing, 2017, 72: 1214-1231

[31]

Tian

W Y

, Zhou

, Deng

W H

A class of second order difference approximations for solving space fractional diffusion equations

Mathematics of Computation, 2015, 84(294): 1703-1727

[32]

Barrett

, Berry

, Chan

T F

et al.

Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods

Philadelphia: Society for Industrial and Applied Mathematics, 1994

[33]

Paige

C C

, Saunders

M A

Solution of sparse indefinite systems of linear equations

SIAM journal on numerical analysis, 1975, 12(4): 617-629

[34]

Burger

, Modersitzki

, Ruthotto

A hyperelastic regularization energy for image registration

SIAM Journal on Scientific Computing, 2013, 35(1): 132-148

[35]

Kelley

C T

Iterative Methods for Optimization

Philadelphia: Society for Industrial and Applied Mathematics, 1999

[36]

Nocedal

, Wright

S J

. Numerical Optimization. Berlin: Springer, 2006

[37]

Sun

, Yuan

Y X

. Optimization Theory and Methods:Nonlinear Programming. Berlin: Springer, 2006

[38]

Modersitzki

. FAIR: Flexible Algorithms for Image Registration. Philadelphia: Society for Industrial and Applied Mathematics, 2009

[39]

Pock

, Urschler

, Zach

et al. A duality based algorithm for TV-L1-optical-flow image registration//International Conference on Medical Image Computing and Computer-Assisted Intervention. Berlin: Springer, 2007: 511-518

[40]

Kroon

D J

Multimodality Non-rigid Demon Algorithm Image Registration

MatlabCentral. http://www.mathworks.com/matlabCentral/fileexchange/21451-multimodality-non-rigid-demon-algorithm-image-registration, 2008

URL [本文引用: 1]

[41]

Kroon

D J

, Slump

C H

MRI modalitiy transformation in demon registration//IEEE International Symposium on Biomedical Imaging: From Nano to Macro

Boston: IEEE, 2009: 963-966