医学图像配准是医学图像分析处理任务中最基本的步骤, 是指通过对一幅图像经过一系列的几何变换, 使变换后图像与参考图像处于相同的坐标空间, 以达到两幅图像中相同的组织结构在空间位置上一一对齐. 配准算法的配准精度和配准效率将直接影响后续的医疗诊断分析和研究结果. 因此, 研究快速、精准的医学图像配准算法具有重要的理论意义和临床实用价值.
传统的配准算法框架主要包含4个部分[1]: 特征空间、搜索空间、搜索策略和相似性度量. 特征空间是指从参考图像和浮动图像中提取可用于配准的特征; 搜索空间定义了图像变换的范围和方式, 即线性变换或非线性变换; 搜索策略是指采用合适最优化算法在搜索空间中寻找最优的配准参数; 相似性度量用于评估每次变换结果的优劣, 为搜索策略的搜索方向提供依据. 在此基础上, 已有大量的配准算法被提出并被广泛用于实际临床应用, 如SyN[2]、Demons[3]、FSL[4]等. 尽管这些配准算法都能取得较好的配准效果, 但是都面临着一个共同的问题: 对于每一对待配准图像对, 传统配准算法均需要从零开始迭代优化目标函数, 严重影响了配准效率[5], 难以满足临床实时的配准需求.
近年来, 深度学习技术被广泛应用于医学图像配准. 目前, 基于深度学习的配准技术主要被分为两大类: 有监督学习配准算法和无监督学习配准算法. 有监督学习配准算法是利用已知几何变换的待配准图像对训练网络, 并通过计算网络估计的几何变换和已知几何变换的均方误差 (mean squared error, MSE) 更新网络的权重参数. 如Miao等人[6]在ConvNet的最后一层设计一个回归器去预测二维X光和三维CT图像的刚性配准参数, 解决了现有基于灰度的二维/三维配准算法速度慢、捕获范围小等两个缺陷. 不同于回归的方法, Yang 等人[7]采用已知几何变换的待配准图像对训练深度编解码网络, 使其能够精准预测大变形微分同胚度量映射 (large deformation diffeomorphic metric mapping, LDDMM) 配准模型[8]中的动量, 随后通过贝叶斯概率网络求得DVF. Sokooti等人[9]则直接使用多尺度ConvNet去估计待配准图像对的DVF, 而无需将几何变换约束至具体的配准模型. 一些作者也尝试利用ConvNet估计薄板样条模型[10]的参数, 如Cao等人[11]和Eppenhof等人[12], 他们分别使用模板样条模型实现了大脑MRI图像和胸部CT图像的可变形配准. 与传统配准算法相比, 有监督学习配准算法极大地提高了配准效率, 但这类方法的配准性能十分依赖带有标签的训练数据. 然而, 获得带有标签的医学图像数据成本非常高, 并且标记信息通常利用传统配准算法获得或者人工合成, 并不精准, 因此, 限制了有监督配准算法的发展.
基于无监督学习的配准算法不需要任何标签, 其通过最小化给定的损失函数约束神经网络预测最优的DVF. Jaderberg等人[13]提出的空域变换网络 (spatial transformer network, STN) 因其能对图像重采样并且可微而极大推动了无监督配准算法的发展. 自STN提出后, 大量无监督配准算法涌现出来. 如Balakrishnan等人[14]提出了一个优秀的无监督配准框架VoxelMorph, 该框架使用一个类似Unet的架构估计稠密的DVF, 然后在训练过程中借助STN模块对浮动图像重采样, 并通过计算重采样后的图像和参考图像的相似性度量指导网络训练, 最终在Dice分数评价指标上取得了与SyN配准算法相当的精度. Zhao等人[15]和Tang等人[16]基于VoxelMorph分别设计了无监督配准框架VTN (volume tweening network) 和 ADMIR (affine and deformable medical image Registration), 他们均在可变形网络框架中新增仿射预配准模型, 使得网络能够同时执行仿射配准和可变形配准, 真正地做到端到端配准. 但是他们也存在区别, VTN框架能够级联多个可变形配准的子网络, 递进式提高配准性能. 实验结果证明了级联策略能够处理大位移形变, 提高配准的精度. 除了模型上的创新, 一些学者也致力于探索几何变换映射的微分同胚性质. 如Zhang[17]提出使用逆一致正则化项来约束相应逆映射的差异, 并在损失函数中引入逆一致性损失和防折叠损失来保证DVF的平滑属性. Mok等人[18]使用类似SyN算法的思想, 利用ConvNet输出一对微分同胚映射, 用于将两幅图像从两条测地线映射到两幅图像的中间地带, 实验结果表明, 该方法能够获得平滑且拓扑保持的DVF.
尽管无监督学习配准算法不需要标签数据, 并且在多个数据集上也取得了十分优秀的配准性能, 但仍存在以下问题: 1)现存的大多数无监督配准算法仅能配准单模态图像, 在多模态图像配准上, 性能低下, 且针对QSM图像和T1加权图像的配准算法仍鲜见报道. 2)现存的配准算法大多采用全局相似度量如MSE、NCC (normalized cross-correlation)、NMI (normalized mutual information)[19] 来驱动网络学习, 但这些相似性度量并不适用于衡量QSM图像和T1加权图像的相似性, 难以驱动网络学习.
为了解决以上问题, 本文设计了RF-RegNet用于QSM图像和T1加权图像配准, 通过编码金字塔提取QSM图像和T1加权图像特征, 两个解码器分别用于估计图像对之间的局部和全局几何形变, 最后利用局部特征相似损失和DVF的正则化共同指导网络模型训练, 完成高效精确地配准.
2 RF-RegNet配准方法分别定义F, M为在n维空间
$ \begin{split} \hat \varphi =& \arg \mathop {\min }\limits_\varphi L(F, M, \varphi ) \\ =& \arg \mathop {\min }\limits_\varphi {L_{\rm sim}}(F, M \circ \varphi ) + \lambda {L_{\rm smooth}}(\varphi ) \end{split} $ | (1) |
其中,
$ \hat \varphi = {f_\theta }(F, M) $ | (2) |
其中,
2.1 编码金字塔
编码金字塔结构类似于拉普拉斯金字塔, 它是由一系列残差图像
构建高斯金字塔是构建拉普拉斯金字塔的前提, 高斯金字塔是通过对原始图像逐级高斯滤波再下采样得到的, 图2中
$ {G_i} = \left\{ {\begin{array}{*{20}{c}} {D(\Phi ({G_{i - 1}})),\quad i \geqslant 1} \\ {{G_0},\quad \quad \quad \quad \;i = 0} \end{array}} \right. $ | (3) |
其中,
$ {L_i} = {G_i} - U(\Phi ({G_{i + 1}})), 0 \leqslant i \leqslant n - 1 $ | (4) |
其中,
本文使用可学习的卷积核代替编码金字塔中的高斯核, 让本文所提出的编码金字塔不仅能够捕获待配准图像对的纹理细节差异, 还能够捕获大脑微小结构的特征差异. 编码金字塔由语义编码和残差编码两部分组成, 其中, 语义编码部分将参考图像和浮动图像堆叠形成一个二通道图像作为输入, 并通过一个步长为1、4个卷积核为
医学图像配准的目标旨在寻找最优的像素级几何形变参数, 需要将编码金字塔学习到的特征映射至像素空间, 因此本文使用与语义编码器和残差编码器对应的语义解码器和残差解码器去完成这一过程. 语义解码器和残差解码器均是通过步长为2、卷积核为
本文使用上下文自相似性度量 (self-similarity context, SSC) 衡量重采样后的QSM图像与参考图像T1图像的相似性从而驱动卷积神经网络学习全局共享参数. 自相似性是用于描述一幅图像内的两个图像块之间的距离度量, 对于三维图像
则在像素点
$ {\textit{SSC}}({x_p}) = \left( \begin{gathered} {{d_{\rm gauss}}(I, x_p^1, x_p^2)} \\ {{d_{\rm gauss}}(I, x_p^1, x_p^3)} \\ {{d_{\rm gauss}}(I, x_p^1, x_p^5)} \\ {d_{\rm gauss}}(I, x_p^1, x_p^6) \\ {d_{\rm gauss}}(I, x_p^2, x_p^3) \\ {d_{\rm gauss}}(I, x_p^2, x_p^4) \\ {d_{\rm gauss}}(I, x_p^2, x_p^6) \\ {d_{\rm gauss}}(I, x_p^3, x_p^4) \\ {d_{\rm gauss}}(I, x_p^3, x_p^5) \\ {d_{\rm gauss}}(I, x_p^4, x_p^5) \\ {d_{\rm gauss}}(I, x_p^4, x_p^6) \\ {d_{\rm gauss}}(I, x_p^5, x_p^6) \\ \end{gathered} \right) $ | (5) |
其中, 对于以图像块
$ {d_{\rm gauss}}(x_p^i, x_p^j) = \exp \left( { - \frac{{{D_p}(x_p^i, y_p^j)}}{{{\sigma ^2}}}} \right) $ | (6) |
其中,
$ {\sigma }^{2}=\frac{1}{12}{\displaystyle \sum _{j=1}^{6}{\displaystyle \sum _{i=1}^{6}{D}_{p}({x}_{i}, {x}_{j}), j\ne i, 且 j\ne i/2\times2+(1-i{\text{%}}2)}} $ | (7) |
提取了两幅图像的自相似特征后, 则
$ {l_{\textit{SSC}}}(F, M \circ \varphi ) = \frac{1}{N}\sum {|{\textit{SSC}}(F) - {\textit{SSC}}(M \circ \varphi )|} $ | (8) |
其中,
2.3.2 正则化约束
图像配准不仅要使得配准结果与参考图像尽可能相似, 同时还需要使得几何变换尽可能平滑. 不平滑的几何变换不满足微分同胚性质, 表明配准后的组织结构被破坏. 为了解决这一问题, 本文采用正则项
$ R(\varphi ) = \sum\limits_{p \in \Omega } {||\nabla \varphi (p)|{|^2}} $ | (9) |
其中,
$ \begin{split} \nabla \varphi (p) =& \left(\frac{{\partial \varphi (p)}}{{\partial x}}, \frac{{\partial \varphi (p)}}{{\partial y}}, \frac{{\partial \varphi (p)}}{{\partial {\textit{z}}}}\right) \\ =& \left( {\begin{array}{*{20}{c}} {\dfrac{{\partial {\varphi _x}(p)}}{{\partial x}}}&{\dfrac{{\partial {\varphi _x}(p)}}{{\partial y}}}&{\dfrac{{\partial {\varphi _x}(p)}}{{\partial {\textit{z}}}}} \\ {\dfrac{{\partial {\varphi _y}(p)}}{{\partial x}}}&{\dfrac{{\partial {\varphi _y}(p)}}{{\partial y}}}&{\dfrac{{\partial {\varphi _y}(p)}}{{\partial {\textit{z}}}}} \\ {\dfrac{{\partial {\varphi _{\textit{z}}}(p)}}{{\partial x}}}&{\dfrac{{\partial {\varphi _{\textit{z}}}(p)}}{{\partial y}}}&{\dfrac{{\partial {\varphi _{\textit{z}}}(p)}}{{\partial {\textit{z}}}}} \end{array}} \right) \end{split} $ | (10) |
因此, 网络模型所使用的总损失
$ {l_{\rm total}} = {l_{\textit{SSC}}}(F, M \circ \varphi ) + \lambda R(\varphi ) $ | (11) |
其中,
本文研究所采用的模板图像为MNI152的T1加权图像, 分辨率为1 mm×1 mm×1 mm. QSM图像数据采集自贵州省人民医院, 共计104名受试者. 所有采集数据均在3T MRI (GE 750w)扫描, 使用16通道磁头线圈. 其采集参数如下: 重复采集时间/回波时间/翻转角度
本文将采集自贵州省人民医院的90幅QSM图像作为训练集, 剩余的14幅QSM图像作为测试集. 训练集用于网络模型预测, 测试集用于评估模型的性能和泛化能力. 本文所提出的方法均在Keras中实现, 使用Adam优化器, 在32 GB的NVIDIA Tesla V100 GPU上进行加速训练, 共迭代5 000次, 训练约18小时左右, 初始学习率设置为0.000 1, BatchSize设置为1, 训练的总轮数Epoch设置为50, 每轮迭代次数Iteration设置为100. 训练RF-RegNet模型的步骤如算法1.
算法1. RF-RegNet网络模型训练算法
初始化:
参考图像
训练数据集
学习率
批次大小
迭代次数
训练轮数
损失函数:
训练开始:
1. For epoch=1 to 500 do
2. 随机打乱训练数据集
3. For iteration=1 to 100 do
随机选择数据集中的一幅图像作为浮动图像M;
4. 利用网络模型预测位移矢量场DVF;
5. 依据DVF对浮动图像
6. 计算损失值并更新网络参数(每一个
7. End for each
8. End for epoch
3.3 评价标准为了有效评估QSM图像和T1加权图像配准的性能, 本文采用目标配准误差 (target registration error, TRE)[23]、Dice分数[24]、Hausdorff距离 (Hausdorff distance, HD) 和平均对称表面距离 (average symmetric surface distance, ASD)[25] 来评估图像之间的表面相似性. 其中, TRE表示同一个标注点在参考图像和配准后图像的差异, 其计算方法如式 (12) 所示:
$ TRE = \sum\nolimits_{i = 1}^n {\sqrt {{{({x_i} - {x_i}')}^2} + {{({y_i} - {y_i}')}^2} + {{({{\textit{z}}_i} - {{\textit{z}}_i}')}^2}} } $ | (12) |
其中,
HD用于反映两个区域的最大差异, 其定义如下:
$ HD(M_W^ * , {F^ * }) = \max \left( {\mathop {\max }\limits_{m \in M_W^ * } \mathop {\min }\limits_{f \in {F^ * }} (m - f), \mathop {\max }\limits_{f \in {F^ * }} \mathop {\min }\limits_{m \in M_W^ * } (m - f)} \right) $ | (13) |
其中,
$ Dice = 2\frac{{|M_W^* \cap {F^ * }|}}{{|M_W^* \cup {F^*}|}} $ | (14) |
Dice分数的取值范围为0–1, 越高的Dice分数值意味着更好配准性能.
最后, 本文也利用ASD来评估图像之间的表面相似性, 其计算方法如式(15)所示:
$ {\textit{ASD}} = \frac{1}{{M_W^* + {F^*}}}\left( {\sum\limits_{m \in M_W^*} {d(m, {F^*})} + \sum\limits_{f \in {F^*}} {d(f, M_W^*)} } \right) $ | (15) |
其中,
$ d(x, Y) = \mathop {\min }\limits_{y \in Y} x - y $ | (16) |
由于QSM是功能图像, 难以分割出图像局部解剖结构, 因此本文仅在全局脑部结构上进行定量评估配准的性能.
4 实验结果与定量分析为了评估本文模型的性能, 本文与VoxelMorph、ADMIR、VTN、Cascaded[26] 4种先进的深度学习配准算法进行了比较. 在对比实验中, 由于对比算法均是针对单模态图像配准的, 因此, 本文采用对比算法中的网络结构, 而损失函数和插值方式与本文提出的算法相同. 不同网络模型的训练过程迭代曲线如图4所示.
从图4可以看出, 相比与其他4种方法, RF-RegNet损失下降速度最快且最低. Cascade、VoxelMorph、ADMIR的损失收敛曲线有大量重叠部分, VTN的损失收敛速度最慢且难以收敛至最优值.
4.1 可视化结果图5为5种配准算法的可视化结果, 其中第1、2列分别为参考图像和浮动图像, 后续5列分别为5种方法的配准结果并使用绿色椭圆框标记出本文方法明显改善的部分解剖区域. 从图中可以直观地看出, 5种方法在全局脑部结构上均取得了较好的配准效果, 但各个方法之间也存在明显差异. VTN、Cascaded两种方法的配准结果更为平滑, 丢失了许多的高频细节信息. VoxelMorph、ADMIR两种方法相比于VTN和Cascaded, 配准效果有一定的提升, 但依然丢失较多的细节信息. 本文所提出的方法取得了最好的视觉效果.
4.2 定量分析
为了进一步验证本文所提出的模型RF-RegNet的性能, 本文使用TRE、Dice分数、HD和ASD四种评价指标在测试集上对全局的配准结果进行了定量分析. 定量结果如表1所示.
从表1中可以直观地观察到, 本文的方法取得了最佳的平均TRE、Dice分数、HD, 其中, TRE相较于其他4种方法分别下降了35.5%、20%、5.4%、6.4%; Dice分数较其他方法分别提高了1.7%、1.1%、0.5%、0.4%; HD较其他几种方法分别降低了2.9%、3%、1.9%、1.6%; ASD较VTN、Cascaded、Voxel-Morph降低了6.2%、4.7%、1.2%, 较ADMIR方法升高了1.9%. 越高的Dice分数和越低的TRE、HD、ASD意味着更好的配准效果. 因此, 本文所提出的模型显示了最好的配准性能. 从配准所需的时间上看, VoxelMorph的网络结构最为简单, 仅包含一个编码器和一个解码器, 因此取得了最快的运行时间. VTN、ADMIR和Cascaded的设计思想均是通过堆叠多个网络模型来实现由粗到细配准, 模型参数量大, 因此需要耗费更多的时间, 本文所提出的方法所需配准时间略高于VoxelMorph, 低于其他方法, 满足临床实时的配准需求.
图6展示了5种方法分别在Dice分数、ASD、HD、TRE四个评价指标上的箱线图, 可以直观地看出, 本文提出的方法相比于其他方法的性能更为稳定, 进一步证明了RF-RegNet的优越性.
5 结论本文针对QSM纹理结构特点, 借鉴拉普拉斯金字塔的思想, 设计了残差编码器和语义编码器, 用于提取和学习配准图像对之间的局部纹理和边缘信息的匹配变形参数. 最后, 通过重采样层重建以得到最后的配准结果. 在网络优化的过程中, 本文采用SSC的特征相似性作为损失函数驱动网络模型学习. 实验结果表明, 本文提出的方法显著提升了QSM图像和T1加权图像的配准精度. 但仍存在些许不足, 主要体现在: 在配准之前需要使用外部工具包进行预对齐, 降低了配准效率. 其次, 本文使用SSC特征相似性作为损失函数, 虽然能够较好衡量T1加权图像和QSM图像的相似性, 但是计算复杂, 致使网络训练速度缓慢.
[1] |
陈显毅. 图像配准技术及其MATLAB编程实现. 北京: 电子工业出版社, 2009.
|
[2] |
Avants BB, Epstein CL, Grossman M, et al. Symmetric diffeomorphic image registration with cross-correlation: Evaluating automated labeling of elderly and neurodegenerative brain. Medical Image Analysis, 2008, 12(1): 26-41. DOI:10.1016/j.media.2007.06.004 |
[3] |
Vercauteren T, Pennec X, Perchant A, et al. Diffeomorphic demons: Efficient non-parametric image registration. NeuroImage, 2009, 45(1S): S61-S72. |
[4] |
Smith SM, Jenkinson M, Woolrich MW, et al. Advances in functional and structural MR image analysis and implementation as FSL. NeuroImage, 2004, 23 Suppl 1: S208–S219.
|
[5] |
邹茂扬, 杨昊, 潘光晖, 等. 深度学习在医学图像配准上的研究进展与挑战. 生物医学工程学杂志, 2019, 36(4): 677-683. |
[6] |
Miao S, Wang ZJ, Liao R. A CNN regression approach for real-time 2D/3D registration. IEEE Transactions on Medical Imaging, 2016, 35(5): 1352-1363. DOI:10.1109/TMI.2016.2521800 |
[7] |
Yang X, Kwitt R, Niethammer M. Fast predictive image registration. In: Carneiro G, Mateus D, Peter L, et al. eds. Deep Learning and Data Labeling for Medical Applications. Cham: Springer, 2016. 48–57.
|
[8] |
Cao Y, Miller MI, Winslow RL, et al. Large deformation diffeomorphic metric mapping of vector fields. IEEE Transactions on Medical Imaging, 2005, 24(9): 1216-1230. DOI:10.1109/TMI.2005.853923 |
[9] |
Sokooti H, de Vos B, Berendsen F, et al. Nonrigid image registration using multi-scale 3D convolutional neural networks. In: Descoteaux M, Maier-Hein L, Franz A, et al. eds. Medical Image Computing and Computer Assisted Intervention—MICCAI 2017. Cham: Springer, 2017. 232–239.
|
[10] |
Rueckert D, Sonoda LI, Hayes C, et al. Nonrigid registration using free-form deformations: Application to breast MR images. IEEE Transactions on Medical Imaging, 1999, 18(8): 712-721. DOI:10.1109/42.796284 |
[11] |
Cao XH, Yang JH, Zhang J, et al. Deformable image registration based on similarity-steered CNN regression. Proceedings of the 20th International Conference on Medical Image Computing and Computer-assisted Intervention. Quebec City: Springer, 2017. 300–308.
|
[12] |
Eppenhof KAJ, Lafarge MW, Moeskops P, et al. Deformable image registration using convolutional neural networks. Proceedings of the SPIE 10574, Medical Imaging 2018: Image Processing. Houston: SPIE, 2018. 192–197.
|
[13] |
Jaderberg M, Simonyan K, Zisserman A, et al. Spatial transformer networks. Proceedings of the 28th International Conference on Neural Information Processing Systems, Volume 2. Montreal: ACM, 2015. 2017–2025.
|
[14] |
Balakrishnan G, Zhao A, Sabuncu MR, et al. VoxelMorph: A learning framework for deformable medical image registration. IEEE Transactions on Medical Imaging, 2019, 38(8): 1788-1800. DOI:10.1109/TMI.2019.2897538 |
[15] |
Zhao SY, Lau T, Luo J, et al. Unsupervised 3D end-to-end medical image registration with volume tweening network. IEEE Journal of Biomedical and Health Informatics, 2020, 24(5): 1394-1404. DOI:10.1109/JBHI.2019.2951024 |
[16] |
Tang K, Li Z, Tian LL, et al. ADMIR—Affine and deform-able medical image registration for drug-addicted brain images. IEEE Access, 2020, 8: 70960-70968. DOI:10.1109/ACCESS.2020.2986829 |
[17] |
Zhang J. Inverse-consistent deep networks for unsupervised deformable image registration. arXiv: 1809.03443, 2018.
|
[18] |
Mok TCW, Chung ACS. Fast symmetric diffeomorphic image registration with convolutional neural networks. Proceedings of 2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition. Seattle: IEEE, 2020. 4643–4652.
|
[19] |
Fu YB, Lei Y, Wang TH, et al. Deep learning in medical image registration: A review. Physics in Medicine & Biology, 2020, 65(20): 20TR01. |
[20] |
Schofield MA, Zhu YM. Fast phase unwrapping algorithm for interferometric applications. Optics Letters, 2003, 28(14): 1194-1196. DOI:10.1364/OL.28.001194 |
[21] |
Kan H, Arai N, Takizawa M, et al. Background field removal technique based on non-regularized variable kernels sophisticated harmonic artifact reduction for phase data for quantitative susceptibility mapping. Magnetic Resonance Imaging, 2018, 52: 94-101. DOI:10.1016/j.mri.2018.06.006 |
[22] |
Paige CC, Saunders MA. LSQR: An algorithm for sparse linear equations and sparse least Squares. ACM Transactions on Mathematical Software, 1982, 8(1): 43-71. DOI:10.1145/355984.355989 |
[23] |
Khamene A, Bloch P, Wein W, et al. Automatic registration of portal images and volumetric CT for patient positioning in radiation therapy. Medical Image Analysis, 2006, 10(1): 96-112. DOI:10.1016/j.media.2005.06.002 |
[24] |
Dice LR. Measures of the amount of ecologic association between species. Ecology, 1945, 26(3): 297-302. DOI:10.2307/1932409 |
[25] |
de Vos BD, Berendsen FF, Viergever MA, et al. A deep learning framework for unsupervised affine and deformable image registration. Medical Image Analysis, 2019, 52: 128-143. DOI:10.1016/j.media.2018.11.010 |
[26] |
Zhao SY, Dong Y, Chang EIC, et al. Recursive cascaded networks for unsupervised medical image registration. Proceedings of 2019 IEEE/CVF International Conference on Computer Vision. Seoul: IEEE, 2019. 10599–10609.
|