波谱学杂志, 2022, 39(3): 291-302 doi: 10.11938/cjmr20212918

研究论文

基于变分推断的磁共振图像群组配准

周勤, 王远军,

上海理工大学 医学影像技术研究所, 上海 200093

Groupwise Registration for Magnetic Resonance Image Based on Variational Inference

ZHOU Qin, WANG Yuan-jun,

Institute of Medical Imaging Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China

通讯作者: 王远军, Tel: 13761603606, E-mail:yjusst@126.com

收稿日期: 2021-05-13  

基金资助: 国家自然科学基金资助项目.  61201067
上海市自然科学基金资助项目.  18ZR1426900

Received: 2021-05-13  

摘要

为解决基于深度学习的成对配准方法精度低和传统配准算法耗时长的问题,本文提出一种基于变分推断的无监督端到端的群组配准以及基于局部归一化互相关(NCC)和先验的配准框架,该框架能够将多个图像配准到公共空间并有效地控制变形场的正则化,且不需要真实的变形场和参考图像.该方法得到的预估变形场可建模为概率生成模型,使用变分推断的方法求解;然后借助空间转换网络和损失函数来实现无监督方式训练.对于公开数据集LPBA40的3D脑磁共振图像配准任务,测试结果表明:本文所提出的方法与基线方法相比,具有较好的Dice得分、运行时间少且产生更好的微分同胚域,同时对噪声具有鲁棒性.

关键词: 深度学习 ; 群组配准 ; 变分推断 ; 可变形配准

Abstract

To address the low precision of pairwise registration method based on the deep learning and the time-consuming nature of traditional registration algorithm, this paper presents a method of unsupervised end-to-end groupwise registration based on variational inference, as well as a registration framework based on normalized cross correlation (NCC) and prior knowledge. The framework can warp all images in the group into a common space and effectively control the deformation field of the regularization, and it doesn't need a real deformation field or a reference image. The estimation of deformation field by this method can be modeled as a probability generation model and solved by variational inference. Then unsupervised training is implemented with the help of spatial transformer network and loss function. The registration results of 3D brain magnetic resonance image from the public data set LPBA40 show that: compared with the baseline method, the proposed method has better Dice score, less running time, better diffeomorphisms domain, and is robust to noise.

Keywords: deep learning ; groupwise registration ; variational inference ; deformable registration

PDF (824KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

周勤, 王远军. 基于变分推断的磁共振图像群组配准. 波谱学杂志[J], 2022, 39(3): 291-302 doi:10.11938/cjmr20212918

ZHOU Qin. Groupwise Registration for Magnetic Resonance Image Based on Variational Inference. Chinese Journal of Magnetic Resonance[J], 2022, 39(3): 291-302 doi:10.11938/cjmr20212918

引 言

各类脑成像中,磁共振成像(magnetic resonance imaging,MRI)由于可以显示血液和其状态的微小变化、无损伤定位大脑的功能活动,帮助病患发现早期脑部疾病,已经成为使用最广泛的脑功能研究手段,并为许多脑部疾病的检测与治疗带来福音.图像配准因具有信息匹配、信息融合的功能而成为在疾病分析中不可或缺的一部分.图像配准是MRI分析中的一项常见任务,也是许多领域中的一个活跃研究课题.它将图像空间对齐到一个共同的解剖空间[1],多幅图像的配准可以作为成对(pairwise,PW)配准或群组(groupwise,GW)配准问题来处理[2].如下图1所示,PW配准指定某一图像作为模板,其他图像与模板配准;GW配准中包括一个联合优化问题,用整个序列的信息创建模板,以避免在后续研究中引入偏差.

图1

图1   (a)成对配准(PW)和(b)群组配准(GW)图像配准示意图.(a)中M代表浮动图像,F为一组图像中任意一幅图像作为模板,φ为变形场,(b)中模板由数据集中N幅图像信息联合构建

Fig.1   Schematic diagrams of image registration for (a) PW and (b) GW. In Fig. (a), M represents a moving image, F represents any image in the group image as a template, and φ is the deformation field. In Fig. (b), the template is jointly constructed by N images in the data set


过去几十年,图像配准算法得到长足的发展.经典的配准算法和基于学习的方法受到很大关注,比较有代表性的算法有elastic模型[3]、B样条[4]和Demons[5].几何变换的微分同胚性是当前配准领域非常看重的优点,比较有代表性的算法有LDDMM[6]、DARTEL[7]和Syn[2].对于以上方法,大多数配准的模板都是任意指定的.但在图像集中随机选择的作为模板的图像往往不能代表图像集的结构变形和复杂性,并可能导致偏差和误导性分析,所以研究GW图像配准具有很重要的意义.近年来GW配准在配准领域越来越受欢迎,是因为它能提供更多有用的信息.Guimond等[8]提出一种建立平均解剖模型的方法,该方法在单个图像(模板图像)中提供平均强度和平均形状,以平均的方式消除了脑形状和强度变化.类似于文献[8]的方法,Seghers等[9]通过选择每个图像作为模板来对齐所有图像,并且使图像与平均变形场非刚性对齐.Wu等[10]提出使用自适应加权策略的SharpMean配准方法,进一步构建基于特征的GW配准方法,该方法在配准过程中实现解剖学上合理的对应.以上提出的几种GW配准算法都是将图像配准到它们的相似图像,其中相似图像称为中间模板[11,12]. Wang等[13]提出一种类似金字塔式的配准框架,该框架可以有效地配准大的图像数据集,配准性能相对较好.随着深度学习技术的兴起,许多研究者也在探索深度神经网络在GW图像配准中的应用.Che等[14,15]提出一个由主成分分析构建的模板图像引导的无偏差的深度GW配准框架,适用于多光谱图像.类似于文献[14,15]的方法,Haase等[16]使用了鲁棒的主成成分分析方法.不同的是,在变分正则化的基础上,他们提出基于一阶原始对偶优化的多级方案来解决由此产生的非参数配准的问题.由于单个模板可能无法捕捉数据集的可变形,Dalca等[17,18]提出一种用于产生条件模板的配准框架和学习策略.类似于文献[17,18]使用可变模板的范例做法,Siebert等[19]提出使用自动编码器实现GW配准的方法,以无监督学习方式学习图像的形状和外观.遵循可变形模板的范例,将其应用到图像集上对齐.最近,He等[20]提出一种无监督的端到端GW框架,该框架具有多步机制来逐步优化输出的变形场,而无需模板.这项工作主要用于二维医学图像配准任务,且预估的变形场不是微分同胚域.Balakrishnan等[21]提出一种基于卷积神经网络(convolutional neural networks,CNN)的无监督PW配准方法,他们使用了一种类似于U-Net的架构,并将其命名为voxelmorph.后来,他们扩展了该方法,提出一个基于无监督学习的推理算法,并将其命名为voxelmorph-diff[22].实验结果表明,他们提出的算法有助于提高Dice系数,性能与ANT和NifTYG相当,但是计算效率是ANT的150倍,是NifTYG的40倍.voxelmorph-diff和voxelmorph是近年来基于深度学习的PW配准领域比较经典的方法,计算效率也高于现有的一般方法.

本文提出一种基于变分推断的无监督端到端的GW图像配准方法.首先,使用平均输入图像的方式构建模板,将输入图像和模板图像在通道处堆叠并送入CNN学习;然后基于变分推断对变形场进行预估,得到的变形场可实现将训练图像变换到模板的公共空间.该方法在训练过程中不需要参考图像和真实的变形场.这种学习策略类似传统的优化迭代方法,但模型和参数由神经网络和其权重代替,可以利用随机梯度下降对模型进行优化.同以往的方法相比,主要创新性如下:(1)不同于传统的GW配准,提出了一个以无监督深度学习方式优化的GW配准框架,它可以同时直接输出所有输入3D图像的位移场,不需要输入参考图像,而且该框架不需要对齐的图像对.(2)为增强所提出框架对噪声和强度变化的鲁棒性,我们使用局部归一化互相关(normalized cross correlation,NCC)来衡量模板和扭曲输入图像之间的相似性损失.(3)与其它基于学习的配准方法相比,所提出的的方法具有网络设计简单,配准精度与Syn算法相当,且比voxelmorph-diff,voxelmorph算法速度更快等特点,并对噪声具有鲁棒性.

1 基于变分推断的无监督配准方法

1.1 相关理论

一般图像配准任务,分为线性配准和可变形配准两步实现.对于大多数可变形图像配准的任务而言,重点在于非线性配准.为遵循这一原则,我们的实验数据均经过线性配准.可变形图像配准是一种非线性配准过程,旨在图像之间建立密集的非线性空间对应关系.假设$M$为浮动图像,$F$为参考图像,$\phi $为变形场,常用的可变配准通过迭代极小化目标函数实现,这一过程可表示为:

$ {\phi ^*} = \arg \mathop {\min }\limits_\phi {L_{{\text{sim}}}}(F,M(\phi )) + \lambda {L_{{\text{smo}}}}(\phi ) $

式中:$ {\phi ^ * } $为最优变形场;$ M(\phi ) $为配准后的图像;$ {L_{{\text{sim}}}}(,) $为相似性度量函数;$ {L_{{\text{smo}}}}({\text{ }}) $为变形场的正则化;$ \lambda $为正则化参数.

1.2 方法流程

本文提出一种基于变分推断的无监督端到端的GW配准框架,该方法得到的预估变形场可建模为概率生成模型,使用变分推断的方法求解;然后借助空间转换网络和损失函数来实现无监督训练.整体流程如图2所示.假设有N幅灰度图像$\left\{ {{I_1}} \right.,{I_2}, \cdots ,\left. {{I_N}} \right\}$,为把这N幅图像配准到公共空间,我们利用N幅输入图像构建模板,模板通过平均配准后的图像来更新.首先,将初始模板和输入图像在通道维度上堆叠,送到CNN网络学习预估变形场,这里我们将预估的变形场建模为概率生成模型,使用变分推断方法求解;然后,将预估的变形场作用于输入图像,得到配准后的图像;随后,更新模板(详细步骤见1.3节).具体来说,对于输入图像和模板图像,我们计算其最佳变换以及由此产生的损失(详细推导见1.4节和1.6节),通过优化损失函数以无监督的学习方式更新模型参数,直至所有图像都与最新的模板对齐.

图2

图2   基于变分推断的磁共振图像群组配准框架

Fig.2   Overview of groupwise registration for magnetic resonance image based on variational inference


1.3 模板更新

GW配准的一个关键点是对组均值图像进行鲁棒和精确的估计.基于微分同胚[2],我们使用如下步骤构建模板.

步骤一:假设我们有N个样本图像$\left\{ {{x_1}} \right.,{x_2}, \cdots ,\left. {{x_N}} \right\}$,取N个样本总和的平均值作为初始模板,即$ {y_0} = \frac{1}{N}\sum\limits_{i = 1}^N {{x_i}} $

步骤二:将初始模板通过CNN配准至N个样本可得N个变形场$ {\phi _i} $,判断所得配准域的均值是否为单位正交阵,即$ \left( {\frac{1}{N}\sum\limits_{i = 1}^N {{\phi _i}} } \right){\left( {\frac{1}{N}\sum\limits_{i = 1}^N {{\phi _i}} } \right)^{ - 1}} $和单位矩阵I是否相等.

步骤三:若不相等,则更新模板$ y = \frac{1}{N}\sum\limits_{i = 1}^N {({x_i} \circ {\phi _i}} ) $,其中$ {x_i} \circ {\phi _i} $表示扭曲第i个输入图像,减少它们的公共偏差.重复第二步,直至生成的变形场的均值为单位正交阵.

1.4 概率生成模型

本文用x表示输入图像,y为构建的模板.假设一个潜在变量zy的强度可由$x \circ {\phi _z}$正态分布表示.概率生成模型图解如图3所示.

图3

图3   概率生成模型图解.白圈表示隐含变量,灰圈表示观测值(输入图像),方块表示参数.虚线矩形框表示里面的变量xyz独立重复N次.配准先验由高斯参数的均值μz和协方差Σz定义

Fig.3   Probability generation model. White circle represents hidden variable, gray circle represents observed value (input image), and squares represent parameters. Dashed rectangle indicates the number of independent repeats N of the x, y, z variables inside. Registration prior is defined by the mean value μz and covariance Σz of Gaussian parameters


通过最大后验概率估计(maximum a posteriori estimation,MAP)获得新图像对(x, y)的最可能变形场$ {\phi _z} $,以及每个体素的速度场方差估计,目的是估计后验配准概率$ p(z|x,y) $.由于后验概率较难计算,用神经网络来近似后验概率$ q(z|x,y) $.基于生成模型,需要最大化边缘似然概率

$ \max L = \sum\limits_{i = 1}^N {\ln p(} {x_i},y) $

根据文献[22],对于每一个样本可知:

$ \ln p(x,y) = KL[{q_\theta }(z|x,y)||p(z|x,y)] + L(\theta |x,y) $

其中$\theta $为神经网络参数;$KL[ \cdot ]$表示KL散度,用来衡量${q_\theta }(z|x,y)$$p(z|x,y)$的接近程度;$L( \cdot )$表示边缘似然概率的下界(evidence lower bound,ELBO).我们将难以求解的边缘似然概率极值问题转化为易于求解的$L( \cdot )$极值问题.因为KL散度非负,则有:

$ \ln p(x,y) \geqslant L(\theta |x,y) = {E_q}[ - \ln {q_\theta }(z|x,y) + \ln p(x,y,z)] $

则有:

$ L(\theta |x,y) = {E_q}[ - \ln {q_\theta }(z|x,y) + \ln p(x,y,z)] $

其中,${E_q}[ \cdot ]$表示对变分分布q的期望.进一步可建模为:

$ \begin{gathered} L(\theta |x,y) = {E_q}[ - \ln {q_\theta }(z|x,y) + \ln p(x,y,z)] \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\text{}} = {E_q}\left[ { - \ln {q_\theta }(z|x,y) + \ln p(z) + \ln \frac{{p(x,y,z)}}{{p(z)}}} \right] \hfill \\ \;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\text{}} = {E_q}[ - \ln {q_\theta }(z|x,y) + \ln p(z)] + {E_q}\left[ {\ln \frac{{p(x,y,z)}}{{p(z)}}} \right] \hfill \\ \;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\text{}} = - KL[{q_\theta }(z|x,y)||p(z)] + {E_q}[\ln p(x|z,y)] \hfill \\ \end{gathered} $

定义公共空间先验的一种方法是对输入图像的潜在变量z进行高斯建模[21],所得概率图谱用作参考.为避免选择固定模板图像带来的偏差并考虑总体信息,对输入图像的潜在变量进行高斯建模,则$ p(z) $表示为:

$ p(z) = N(z;0,{\sum _z}) $

其中$ N(z;0,{\Sigma _z}) $是均值为0、协方差为$ {\sum _z} $的高斯分布函数.我们假设${\sum ^{ - 1}}{\text{ = }}l$,其中像素网格上定义的邻域图拉普拉斯算子$l = {\lambda _a}A - {\lambda _b}B$A为度矩阵,B为像素邻域邻接矩阵.${\lambda _a}$${\lambda _b}$是超参数,它们分别控制由z参数化的域中向量的范数及其位移的范数.z表示速度场.

本文研究单模态变形配准,使用NCC来测量衡量模板和扭曲输入图像之间的相似性损失,以增强其对噪声和强度变化的鲁棒性.任意两幅图像的局部归一化互相关可由以下公式计算:

$ NCC(f,g) = \frac{1}{{\left| \Omega \right|}}\sum\limits_{x \in \Omega } {\frac{{\sum\nolimits_{{x_i}} {[f({x_i}) - } \bar f({x_i})][g({x_i}) - \bar g({x_i})]}}{{\sqrt {\hat f({x_i})\hat g({x_i})} }}} $

其中$\Omega $为任一非空集合,$f,g:\Omega \to R$$\Omega \subset {R^3}$,表示非空3D图像集中的图像;$ \bar f(x) = \sum\nolimits_{{x_i}} {f({x_i})/{n^3}} $$ \hat f(x){\text{ = }}\sum\nolimits_{{x_i}} {[f({x_i}) - \bar f} (x){]^2} $分别表示局部均值和方差图像,${x_i}$表示在体素$x$处尺寸大小为$ {n^3} $的体积上循环遍历.类似文献[20],这里我们使用$n = 5$.因为NCC对强度变化是鲁棒的,本文将$x$的概率建模为归一化互相关的指数分布,如下:

$ p(x|z,y) \propto \exp \left\{ {\frac{1}{N}\sum\limits_i {NCC({x_i} \circ {\phi _{{z_i}}},y} )} \right\} $

1.5 神经网络模型

GW方法旨在找到GW图像到模板空间的最佳变换,以及模板图像的计算.本文的网络架构包括解码网络、CNN、空间转换器和损失函数.在文献[22]中,Balakrishnan等明确的选择模板图像和运动图像来进行PW配准.相比之下,在我们的方法中,在输入到CNN网络之前,先将输入图像和生成的模板图像在通道堆叠并送入CNN网络,学习群组图像到模板图像的公共空间变换.

CNN模型由卷积、缩小/放大和跳跃连接组成,详细结构如图4所示该架构由编码器和解码器两部分组成.网络采用通过将xy连接成2通道3D图像而形成单个输入.在本文实验中,输入大小为160*160*192*2.编码阶段和解码阶段均使用内核大小为3的3D卷积,步长为2. 在每个卷积操作后紧跟一个参数为0.2的LeakyReLU激励函数,增加网络的非线性,且步长均为1. 编码阶段使用4个步长为2、大小为3*3*3的卷积核进行下采样(最大值池化),使空间维度减半.编码阶段最终输出一个10*10*12*32的特征图.解码阶段使用4个上采样层和1*1*1卷积核.最终输出两个160*160*192*3大小的特征图分别对应近似后验概率的均值和协方差.由于上采样过程易造成信息丢失,我们在编码阶段和解码阶段之间加入跳跃连接(skip connection).其中,跳跃连接层可以将编码阶段阶段学到的信息对应传到解码阶段阶段,起到了提供更精细空间尺度的作用,从而补充信息,使得精确配准成为可能.

图4

图4   本文使用的CNN网络框架

Fig.4   CNN network frame in this research


总体结构与大多数医学图像配准网络使用的U-net相同.但是,为满足GW配准的需要,做了以下几处改变:

(1)在原始的U-net中,下尺度层和上尺度层是通过最大池化和转置卷积实现的,取而代之的是一个更为简单的群组网格.

(2)批次数量将始终为一,因为在优化过程中,只有一组图像被送入网络.因此,批处理规范化被实例规范化所取代.

(3)卷积-归一化-激活操作的两个连续集合被减少到一个.这一变化提高了效率,但不影响性能.使用泄漏整流激活层(LeakyReLU)来代替原始整流线性激活(ReLU).

(4)由于内存大小限制,输入图像在输入CNN之前被缩小到原始分辨率的一半,然后将输出位移场提升到原始分辨率,用以变换输入图像.

1.6 损失函数

制定的损失函数旨在优化模型和减少公共空间生成的偏差.通过不断优化损失函数从而指导模型的学习.本文使用的损失函数为

$\begin{aligned}-L(\theta \mid x, y) &=-E_q[\ln p(x \mid z, y)]+K I\left[q_\theta(z \mid x, y) \| p(z)\right] \\&=-\frac{1}{N} \sum_i N C C\left(x_i \circ \phi_{z i}, y\right)+\frac{1}{2}\left\{t r\left[\left(\lambda_a A-\lambda_b B\right) \Sigma_z-\ln \Sigma_z\right]+\mu_z^T \Lambda^{-1} \mu_z^T\right\}+\text { const }\end{aligned}$

其中$ \mu _z^T{\Lambda ^{ - 1}}\mu _z^T{\text{ = }}\frac{{{\lambda _a}}}{2}\Sigma {\Sigma _{j \in N(q)}}{(\mu [q] - \mu [j])^2} $$ N(q) $为体素q的邻域;$ \mu [q] $$ \mu [j] $分别代表体素q的均值和体素j的均值;const表示常数,为方便计算,实验中我们使用const=0.第一项控制模板和扭曲输入图像之间的相似性损失,第二项中$ \frac{1}{2}tr[({\lambda _a}A - {\lambda _b}B){\Sigma _z} - \ln {\Sigma _z}] $是使后验损失尽可能的靠近先验损失,第二项中$ \frac{1}{2}\mu _z^T{\Lambda ^{ - 1}}\mu _z^T $是控制空间平滑损失.我们使用神经网络实现我们的概率生成模型,并使用随机梯度下降算法来优化模型参数.

2 实验部分

2.1 实验描述

为验证本文提出算法的性能,使用公开数据集LONI LPBA40[23]进行实验.对数据集的所有3D磁共振图像都使用FSL和FreeSurfer软件进行偏移场矫正、大脑提取及线性配准、体素重采样至 1 mm*1 mm*1 mm和尺寸裁剪为160*160*192等预处理.并将含有40个3D磁共振图像的LPBA40数据集划分为训练集(30)、测试集(10).每个3D图像均有对应的标签(54个脑区标记),这些标签可用于评价配准的精度.

2.2 深度学习模型训练与测试

本文构建的深度模型使用无监督的学习方式进行训练.将初始模板和训练数据在通道维度上堆叠,然后送到CNN网络学习得到变形场,将得到的变形场作用于训练图像,输出配准后的图像,随后,更新模板直到所有图像都与最新的模板图像对齐.同训练一样,测试时将测试样本的总和取平均输入到训练好的模型进行多次迭代,计算迭代后的平均变形场.当迭代的变形场接近正交网格时,迭代终止,得到的最新模板作为测试时的模板图像.将得到的模板图像和测试图像输入到训练好的模型,将得到的变形场作用于测试图像标签,得到形变后的标签图,通过这种方式计算不同标签的Dice值.

我们使用Tensorflow在6-core Intel i7-8700K CPU和6 GB NVIDIA GeForce RTX 2060 GPU机器上实现模型的训练与测试.实验环境为DUDA10.0并行计算架构,操作系统为Win10,软件为PyCharm.并设计了以下两个不同的实验:(1)LPBA40数据集配准实验,以评估算法的配准精度;(2)在LPBA40数据集上叠加不同参数的噪声,以分析算法对噪声的鲁棒性.

本文还选取配准精度较高的Syn、voxelmorph-diff和voxelmorph三种算法进行了对比,其中后面两种方法是基于深度学习的经典的可变形PW配准算法.Syn使用ANTsPy实现,使用互相关(cross correlation,CC)系数作为其度量标准.voxelmorph-diff和voxelmorph两种方法使用与文献[22]相同的网络参数进行训练.

2.3 评估方法

参考我们以前的工作[1,24,25],通过Dice、$E{M_1}$$E{M_2}$、变形场的雅克比行列式和时间五个度量标准来对GW配准结果进行衡量.

Dice表示模型的预测结果与真实值的相似程度,其值越高,图像配准效果越好.$ S_{M(\phi )}^K $$ S_F^K $分别表示$ M(\phi ) $F的体素,K为对应标签的序号,计算公式为

$ {\text{Dice(}}S_{M(\phi )}^K,S_F^K) = 2*\frac{{\left| {S_{M(\phi )}^K \cap S_F^K} \right|}}{{\left| {S_{M(\phi )}^K} \right| \cup \left| {S_F^K} \right|}} $

$E{M_1}$为平均偏差指数,表示参考模板F和总体之间的偏差[1]N为样本图像个数,D为欧几里得距离测度.$E{M_1}$值越小,偏置越小;反之,偏置越大.$E{M_1}$计算如下:

$ E{M_1} = \frac{1}{N}\sum\limits_{i = 1}^N {D(F,{x_i}} ) $

$E{M_2}$为不同模板间的锐度差,可用平均梯度法计算.Q为模板F在点p处的体素总数.$E{M_2}$越小代表模板越清晰.$E{M_2}$计算如下:

$ E{M_2} = \frac{1}{Q}\sum\limits_p {\left\| {\nabla F(p)} \right\|} $

变形场的雅克比行列式中负值的个数一般用来评价变形场的重叠程度.变形场雅克比矩阵$ {J_\phi }(p) = \nabla \phi (p) $,我们计算雅克比行列式中体素小于0的点,即$\det ({J_\phi }(p)) \leqslant 0$.雅克比行列式的负值越少,表示具有更好的微分同胚变形场.

3 结果与讨论

3.1 LPBA40数据集配准实验

由于训练集数据较少,我们在不发生图像形变的前提下,使用左右翻转和上下翻转的方式进行数据扩增[26],以大大增强模型的泛化能力,减少模型过拟合的情况.算法在Tensorflow框架中实现,使用学习率为1e-5的Adam优化器优化训练,迭代次数为15 000次.超参数$ {\lambda _a} $${\lambda _b}$根据网络收敛性调参,最终分别设置为1和0.04.由于实验环境内存大小的限制,仅将批处理大小设为1.测试时,使用LPBA40数据集中的10个样本构建模板图像,并将其作为CNN网络配准的模板图像.我们对平均输入图像做了2次迭代和2次平均变形场,如图5所示.因为第2次配准域接近正交网格,停止迭代,并将其作为测试的模板图像.

图5

图5   测试结果中心切片:(a)第一迭代的变形场;(b)第二次迭代的变形场;(c)浮动图像;(d)模板图像;(e)配准后的图像;(f)变形场

Fig.5   Central slice of test results: (a) Deformation field of the first iteration; (b) Deformation field of the second iteration; (c) Moving image; (d) Template image; (e) Image after registration; (f) Deformation field


表1给出不同算法在测试集上的配准结果.从表1中,可以看出本文提出的方法,在大多指标上性能表现优异(粗体表示最优结果).同Syn、voxelmorph-diff和voxelmorph相比,它具有较高的Dice得分、GPU运行时间少,产生更好的微分同胚变形场(具有更低的非负雅可比位置数对),生成的平均模板更清晰($E{M_2}$相对较小).综合来看,voxelmorph表现较差,可能是训练时随机选取的参考图像的外观离公共的解剖空间较远.Syn的$E{M_1}$略小于我们的方法,这可能是我们的方法使用了过多的平滑和正则化,累计较多的近似误差造成的.图6给出本文算法的部分结果图.

表1   不同算法在LPBA40测试集上配准结果的多指标分析

Table 1  Groupwise registration results of multiple metrics on LPBA40 test set by different algorithms

方法Avg. Dice*GPU运行时间/s*Percentages (%) of ${(\det ({J_\phi }(p)) \leqslant 0)^*}$EM1EM2
ANTs (Syn)0.689 (0.036)45.6 (19)930.986.54
voxelmorph-diff0.698 (0.030)0.51 (0.03)49.6 (11)960.676.29
voxelmorph0.672 (0.031)0.49 (0.02)48.9 (26)969.467.08
本文方法0.701 (0.025)0.40 (0.02)45.0 (18)940.895.07

*:括号内指方差. Avg. Dice为54个脑区配准后Dice的平均值.

新窗口打开| 下载CSV


图6

图6   基于本文算法得到的10幅脑磁共振图像冠状位中心切片测试集群组配准结果.(a)原图像;(b)配准后的图像;(c)变形场的彩色图像;(d)变形场的网格图像

Fig.6   Results of groupwise registration of coronal section central slices of 10 brain magnetic resonance images by the proposed algorithms. (a) The original image; (b) The registered image (c) The color image of the deformation field; (d) The grid image of the deformation field


GW配准的重点是将每幅图像变形到的公共空间,因此获得一幅具有代表性的模板图像是组配准的关键.其核心是找到一个无偏的模板图像,将其他图像配准到图像的公共空间.但很难验证模板的无偏向,我们取相关算法配准后的组均值图像,如图7所示,结果表明,我们的方法生成的模板图像比voxelmorph和voxelmorph-diff方法更清晰,其结构与Syn方法的结果大致一致,这意味着我们的方法至少与基线一样无偏.

图7

图7   不同算法得到的扭曲图像的均值图像.(a)矢状位中心切片;(b)冠状位中心切片;(c)水平位中心切片.每幅小图中第一行从左到右依次为voxelmorph算法的均值图像,voxelmorph-diff算法的均值图像.第二行从左到右依次为Syn算法的均值图像和本文方法的均值图像

Fig.7   Mean images of distorted images by different algorithms. (a) Sagittal section central slice; (b) Coronal section central slice; (c) Horizontal section central slice. In each sub-figure, the first row from left to right is the mean image by voxelmorph algorithm, and voxelmorph-diff algorithm, the second row from left to right is the mean image by Syn algorithm and the method proposed in this work


3.2 本文方法对噪声的鲁棒性分析

为验证本文所提算法对噪声的鲁棒性,我们设计噪声干扰实验.对LPBA40数据集训练集添加均值为0,方差分别为0.001、0.002、0.003、0.004的高斯噪声,训练集扩增至共150幅.对测试集做同样的处理,扩增至50幅,部分结果如图8所示.算法在Tensorflow框架中实现,使用学习率为1e-5的Adam优化器优化训练,迭代次数为15 000次.超参数${\lambda _a}$${\lambda _b}$分别设置为1和0.04.由于Syn迭代,voxelmorph和voxelmorph-diff模型训练需要大量的时间,短时间内工作站计算能力不足,因此仅对本文提出的算法做了纵向对比.配准结果指标计算如表2所示.高斯噪声均值为0,方差为0.002的配准结果如图9所示.从表2可以看出,随着噪声方差由0.001增加至0.002时,我们得到Dice和得到的微分同胚场的指标较好且保持稳定.事实上,由图8可以看出,当方差增加至0.003时,图像的质量已经很差.因添加较多的高斯噪声,图像灰度差值分布更广泛,噪声会影响NCC的分布,进而影响基于NCC作为相似性度量的配准方法的性能.从图9中可以看出,与浮动图像相比,配准后的图像显示更少的噪声.综合分析来看,我们提出的算法对噪声具有一定的鲁棒性.

图8

图8   不同噪声强度的浮动图像:(a)~(d)噪声均值均为0,方差分别为0.001、0.002、0.003、0.004

Fig.8   Moving images with different noise intensities: (a)~(d) The mean values of noise are 0, and the variances are 0.001, 0.002, 0.003, and 0.004, respectively


表2   使用本文方法对包含不同强度噪声的LPBA40测试集进行配准的多指标分析

Table 2  Groupwise registration results of multiple metrics on LPBA40 testing set with different noise intensities using the proposed method

噪声方差(均值为0)Avg. Dice*GPU运行时间/s*Percentages (%) of ${(\det ({J_\phi }(p)) \leqslant 0)^*}$EM1EM2
0.0010.686 (0.026)0.50 (0.02)48.8 (21)1330.9789.54
0.0020.683 (0.030)0.60 (0.04)49.3 (27)1360.67100.29
0.0030.646 (0.031)0.58 (0.02)51.7 (19)1469.4699.08
0.0040.656 (0.020)0.62 (0.03)54.5 (28)1480.89113.07

*:括号内指方差. Avg. Dice为54个脑区配准后Dice的平均值.

新窗口打开| 下载CSV


图9

图9   使用本文算法对含均值为0、方差为0.002的噪声图像进行配准的结果:(a)浮动的图像;(b)构建的模板;(c)配准后的图像

Fig.9   Groupwise registration result of LPBA 40 testing set with noise of mean value=0 and variance=0.002 using the proposed method: (a) Moving image; (b) Constructed template; (c) Image after registration


4 结论

为解决传统配准算法配准时间长和精度低的问题,本文提出了一种基于变分推断的无监督GW配准方法.该方法将整个序列的信息融入到配准过程中,对预估的变形场使用概率生成模型进行建模,并用变分推断的方法求解.此外,我们设计一种基于NCC和先验的配准方法,以有效地控制变形场的正则化.利用LPBA40数据集配准实验以及方法相对噪声的鲁棒性分析实验,验证本文配准方法的性能.实验结果表明,我们的方法在多数评价指标上表现良好且具对噪声具有一定的鲁棒性.而且我们的GW配准模型可以同时配准多幅图像,其性能类似于Syn和基于学习的voxelmorph-diff,接近实时的GW图像配准.不足的是,尽管使用了微分同胚变换,但雅克比行列式也有负值且损失函数使用了过多的平滑项和正则化项,这可能增加近似误差的累积,并导致配准的速度下降.同时,本文对算法鲁棒性实验设计简单,后续可以考虑使用多线圈图像进行测试,并考虑将算法应用于多模态图像间的配准.

利益冲突


参考文献

WANG Y, JIANG F, LIU Y

Reference-free brain template construction with population symmetric registration

[J]. Med Biol Eng Comput, 2020,58 (9): 2083- 2093.

DOI:10.1007/s11517-020-02226-5      [本文引用: 3]

MARTÍN-GONZÁLEZ E, SEVILLA T, REVILLA-ORODEA A, et al

Groupwise non-rigid registration with deep learning: an affordable solution applied to 2D cardiac cine MRI reconstruction

[J]. Entropy, 2020,22 (6): 687.

DOI:10.3390/e22060687      [本文引用: 3]

GEE J C, REIVICH M, BILANIUK L, et al

Evaluation of multiresolution elastic matching using mri data

[J]. Proc Spie, 1991,1445,226- 234.

DOI:10.1117/12.45220      [本文引用: 1]

RUECKERT D, SONODA L I

Nonrigid registration using free-form deformations: application to breast mr images

[J]. IEEE T Med Imaging, 1999,18 (8): 712- 721.

DOI:10.1109/42.796284      [本文引用: 1]

THIRION J P

Image matching as a diffusion process: an analogy with Maxwell's demons

[J]. Med Image Anal, 2011,2 (3): 243- 260.

[本文引用: 1]

ZHONG J, PHUA D, QIU A

Quantitative evaluation of lddmm, freesurfer, and caret for cortical surface mapping

[J]. Neuroimage, 2010,52 (1): 131- 141.

DOI:10.1016/j.neuroimage.2010.03.085      [本文引用: 1]

ASAMI T, BOUIX S, WHITFORD T J, et al

Longitudinal loss of gray matter volume in patients with first-episode schizophrenia: dartel automated analysis and roi validation

[J]. Neuroimage, 2012,59 (2): 986- 996.

DOI:10.1016/j.neuroimage.2011.08.066      [本文引用: 1]

GUIMOND A, MEUNIER J, THIRION J P

Average brain models: a convergence study

[J]. Comput Vis Image Und, 2000,77 (2): 192- 210.

DOI:10.1006/cviu.1999.0815      [本文引用: 2]

SEGHERS D, D'AGOSTINO E, MAES F, et al. Construction of a brain template from mr images using state-of-the-art registration and segmentation techniques[C]// Medical Image Computing and Computer-Assisted Intervention--MICCAI 2004, 7th International Conference Saint-Malo, France, September 26-29, 2004, Proceedings, Part I. 2004.

[本文引用: 1]

WU G, JIA H, WANG Q, et al

Sharpmean: groupwise registration guided by sharp mean image and tree-based registration

[J]. Neuroimage, 2011,56 (4): 1968- 1981.

DOI:10.1016/j.neuroimage.2011.03.050      [本文引用: 1]

WU G, WANG Q, JIA H, et al

Feature-based groupwise registration by hierarchical anatomical correspondence detection

[J]. Hum Brain Mapp, 2012,33 (2): 253- 271.

DOI:10.1002/hbm.21209      [本文引用: 1]

YANOVSKY I, THOMPSON P M, OSHER S, et al. Topology preserving log-unbiased nonlinear image registration: theory and implementation[C]// IEEE Conference on Computer Vision & Pattern Recognition. Minneapolis, Minnesota, USA: IEEE, 2007.

[本文引用: 1]

WANG Q, CHEN L Y YAP P T, et al

Groupwise registration based on hierarchical image clustering and atlas synthesis

[J]. Hum Brain Mapp, 2010,31,1128- 1140.

[本文引用: 1]

CHE T, ZHENG Y, CONG J, et al

Deep group-wise registration for multi-spectural images from fundus images

[J]. IEEE Access, 2019,7,27650- 27661.

DOI:10.1109/ACCESS.2019.2901580      [本文引用: 2]

CHE T, ZHENG Y, SUI X, et al. Dgr-net: deep groupwise registration of multispectral images[C]// Information Processing in Medical Imaging - 26th International Conference , Hong Kong, china: IPMI, 2019: 706-717.

[本文引用: 2]

HAASE R, HELDMANN S, LELLMANN J. Deformable groupwise image registration using low-rank and sparse decomposition[EB/OL]. [2020-06-10].https://arxiv.org/abs/2001.03509.

[本文引用: 1]

DALCA A V, RAKIC M, GUTTAG J, et al. Learning conditional deformable templates with convolutional networks[C]//Neural Information Processing Systems 2019(NeurIPS 2019), Vancouver, BC, Canada, 2019: 804-816.

[本文引用: 2]

YU E M, DALCA A V, SABUNCU M R Learning conditional deformable shape templates for brain anatomy[M]. Lima: Machine Learning in Medical Imaging, 2020, 353- 352.

[本文引用: 2]

SIEBERT H, HEINRICH M P Deep groupwise registration of mri using deforming autoencoders[M]. Berlin: Springer, 2020, 236- 241.

[本文引用: 1]

HE Z Y, CHUNG A C S. Unsupervised end-to-end groupwise registration framework without generating templates[C]// 2020 IEEE International Conference on Image Processing (ICIP). IEEE, 2020: 375-379.

[本文引用: 2]

DALCA A V, BALAKRISHNAN G, GUTTAG J, et al

Unsupervised learning of probabilistic diffeomorphic registration for images and surfaces

[J]. Med Image Anal, 2019,57,226- 236.

DOI:10.1016/j.media.2019.07.006      [本文引用: 2]

BALAKRISHNAN G, ZHAO A, SABUNCU M R, et al

Voxelmorph: a learning framework for deformable medical image registration

[J]. IEEE T Med Imaging, 2019,1788- 1800.

[本文引用: 4]

SHATTUCK D W, MIRZA M, ADISETIYO V, et al

Construction of a 3d probabilistic atlas of human cortical structures.

[J]. Neuroimage, 2008,39 (3): 1064- 1080.

DOI:10.1016/j.neuroimage.2007.09.031      [本文引用: 1]

王远军, 刘玉

基于图像集拓扑中心的群体配准方法

[J]. 波谱学杂志, 2018,35 (4): 60- 67.

URL     [本文引用: 1]

WANG Y J, LIU Y

Group registration method based on topological center of images

[J]. Chinese J Magn Reson, 2018,35 (4): 60- 67.

URL     [本文引用: 1]

蔡文琴, 王远军

基于磁共振成像的人脑图谱构建方法研究进展

[J]. 波谱学杂志, 2020,37 (2): 241- 253.

URL     [本文引用: 1]

CAI W Q, WANG Y J

Advances in the construction of the human brain map based on magnetic resonance imaging

[J]. Chinese J Magn Reson, 2020,37 (2): 241- 253.

URL     [本文引用: 1]

刘可文, 刘紫龙, 汪香玉, 等

基于级联卷积神经网络的前列腺磁共振图像分类

[J]. 波谱学杂志, 2020,37 (2): 152- 161.

URL     [本文引用: 1]

LIU K W, LIU Z L, WANG X Y, et al

Prostate magnetic resonance image classification based on cascading convolutional neural networks

[J]. Chinese J Magn Reso, 2020,37 (2): 152- 161.

URL     [本文引用: 1]

/