基于B样条仿射和GANs的医学图像配准

宋 枭，朱家明，徐婷宜

(扬州大学 信息工程学院，江苏 扬州 225000)

0 引言

图像配准是图像研究领域的一个典型问题和技术难点，其目的在于比较或融合针对同一对象在不同条件下获取的图像。随着医学成像设备的不断进步，可以采集更加精准的器官解剖图像，但仅凭医务工作者的空间想象和主观经验观察病情的发展十分困难。精确快速的医学图像配准是医学深入分析的基础，是分析病灶、观察器官变化、引导临床手术的关键。为此，广大学者提出了许多医学图像配准方法，从传统的基于特征和灰度的方法到深度学习方法，图像配准的性能得到了质的飞跃。

传统的配准大多数在变分框架中实现[1]，即迭代最小化变形空间中的能量问题，计算量大且受局部极值影响，鲁棒性差。深度学习方法直接通过深层网络提取不同层次的灰度、结构特征，最大化配准图像之间的相似度。Maes等人[2]提出的最大化互信息准则的配准方法，比较灰度概率分布来度量图像配准程度，利用powell算法进行参数优化，当图像对齐时互信息最大，该方法易收敛至局部极值，且采样及插值算法对配准结果影响较大；Panigrahi等人[3]提出的基于仿射变换的刚性配准，通过优化多个空间形变多项式的参数，消除了因成像角度、时间造成的不规则失真，但对于局部形变，不能达到良好的配准效果；Fan等人[4]提出的基于生成对抗网(Generative Adversarial Network，GAN)的图像配准模型，通过最小化对抗损失，直接输出配准图像，并未对配准图像进行空间变换以及相似性度量，仅用概率表示配准效果，难以达到精准的配准效果。为了解决上述问题，本文将B样条、仿射变换与GAN结合，提出了二次判别机制，对B样条仿射变换结果进行初次判别，以弱化全局及局部形变影响，利用图像小波变换进行插值得到粗配准结果，最后利用GAN最大化相似性度量输出配准结果。

1 B样条仿射变换模型

1.1 B样条仿射变换

仿射变换(AF)源于缩放、旋转、剪切和平移变换的组合的全局线性变换，能保证方向、比例的一致性[5]。但仿射变换的全局性不能解决图像局部形变带来的配准问题，为此引入B样条，对形变图像进行粗到精配准。

B样条仿射变换的目标是找到全局形变场Tglobal以及局部形变场Tlocal，将浮动的图像Ifloat映射到固定(参考)图像Ifixed，表达为：

Ifixed=Tlocal(TglobalIfloat)。

(1)

由于B样条的引入，仿射变换的粗配准只需6个参数来完成图像的平移、旋转和缩放，不需要额外多的控制点加入。引入齐次坐标，增广维度，全局形变场Tglobal可写成：

(2)

式中，a11，a12，a13，a21，a22，a23是仿射变换的约束参数，仿射变换的过程就是最优化上述参数的过程。

为了引入局部控制，将浮动图像网格化，利用控制点的位置更新实现局部配准。假设二维图像(x,y)，利用m×n网格控制像素点在x轴和y轴上的移动，可以表达为：

(3)

式中，i，j表示控制点编号；u，v表示与x，y的相对位置。3次B样条B0，B1，B2，B3表达式为：

(4)

1.2 模型代价函数

C=Csimilarity(T(Ifloat),Ifixed)+λCsmooth(T)，

(5)

式中，λCsmooth为平滑约束项；Csimilarity为相似性度量，该度量由归一化误差均值以及内点比构成，表示为：

Csimilarity(T(Ifloat),Ifixed)=E(T(Ifloat),Ifixed)-

ninliers/(ninliers+noutliers)。

(6)

1.3 插值实现

变换后的像素不一定落在整数上，需对非整数坐标上的点进行灰度插值，通过内插点附近的像素的灰度级来估计插值灰度[8]。由于线性插值在插值过程中采用插值内核而不考虑像素点所处位置，会产生方块效应从而导致边缘模糊、纹理特征损失，显然不符合医学图像配准的要求，故本文采用基于小波变换(DWT)的非线性插值方法[9]。该方法往往具有多分辨率分析和逐渐局部细化等优点，其基本思想是将信号分解到不同的尺度或者分辨率层上，因此可以在不同的尺度上独立对信号进行研究。

张建勋等人[10]对原始图像进行小波分解，用邻接子带的相似性计算高频分量，最后对原始图像低频分量做插值运算。小波变换的插值实现过程主要由小波分解和小波逆变换组成，正交小波分解不仅可将图像的高低频信息很好地分离，而且分解后各层子带之间具有相似性[11]。分解后的低频子带信号ML中包含了图像的绝大部分能量；高频子带信号MH、MV、MD则对应图像的边缘信息。通过小波变换，将图像的高低频信息分离后，可以单独对高频信息进行处理。低分辨率中的高频分量MH2、MV2、MD2相似于高分辨率中的MH1、MV1、MD1。利用重构理论，将得到的高频与原有的低频相叠加，离散小波逆变换就可以得到分辨率更高的图像[12]，如图1所示。

图1 小波变换Fig.1 Wavelet transform

2 构建配准模型

2.1 生成对抗网络

GANs基于两人零和对策模型，单独交替迭代训练生成模型和判别模型，自动学习原始真实样本集的数据分布，假样本在训练过程中真假变换，以达到纳什平衡[13]。

本文引入二次判别机制，由仿射变换生成变形场，经DWT插值得到中间配准图像，由一次判别器判别，根据RANSAC算法有限迭代修正参数，直到一次判别为真，再输入到生成模型，直到二次判别输出为真。整个网络由B样条仿射变换和GAN两部分组成，B样条仿射变换利用RANSAC算法迭代优化数学模型参数生成粗变形场，再由GAN网络完成学习。模型结构如图2所示。

图2 本文配准模型Fig.2 Proposed registration model

2.2 构造生成和判别网络

生成器由卷积层、密集残差块和反卷积层组成。张桂梅等人[14]在GAN中加入5层密集残差块提取深层特征，提高了配准精度。由于仿射变换的加入，采用一层残差块以减小网络的学习成本。由4层卷积层组成的残差块进行下采样，每层均由64个3×3步长为1的卷积核、激活函数和归一化层组成。由6组反卷积层进行上采样，最后经卷积层输出。生成网络结构如图3所示。

图3 生成网络结构Fig.3 Structure of generative network

生成器输出作为判别器输入，经7个卷积层、2个全连层，最终由分类器输出结果。其中卷积层由卷积核、归一化以及激活函数组成，第1个卷积层使用32个步长为1的3×3卷积核，第2个和第3个使用64个步长为1的3×3卷积核，第4个和第5个使用128个步长为1的3×3卷积核，第6个和第7个使用256个步长为1的3×3卷积核。判别网络结构如图4所示。

图4 判别网络结构Fig.4 Structure of discriminate network

2.3 构造损失函数

GAN的判别器可以认为是一种自主学习的损失函数[15]，区别于其他固定的损失函数(如MAE、MSE等)，判别器自适应地度量Itrans，Ifixed之间连续的概率分布。判别网络从真实数据分布pdata、生成数据分布pG中采样，用于判别网络的损失计算：

θG,θD=argminmaxV(G,D)=Ex～pdata[lg(D(x))]+

Eg～pG[lg(1-D(G(z)))]，

(7)

式中，G，D分别表示生成网络和判别网络；x，g分别为真实数据与生成数据；pdata为真实数据分布；pG为生成数据分布，且x=G(z);z表示生成网络输入。对于GAN，生成网络的期望是最小化损失函数，而判别网络的期望是最大化损失函数。

2.3.1 判别器损失函数构造

本文二次判别机制需要构造2个判别损失函数LD1，LD2，第一判别器主要目的是优化仿射变换的6个参数以估计待配准图像对之间的形变大小，即最大化inliers数量，写成归一化形式：

(8)

由于判别器目标是最大化真实数据pdata的期望lg(D(x))，最小化生成数据PG的期望lg(D(G(z)))，因此第二判别器的期望是最大化损失：

LD2=argmaxV(G,D)=Ex～pdata[lg(D(x))]+

Eg～pG[lg(1-D(G(z)))]。

(9)

2.3.2 生成器损失函数构造

GAN的生成器目的在于使输出无限接近真实数据，用p(x;θG)表示生成器输出G(z)与真实数据分布pdata相似的概率，因此需要p(x;θG)越大越好，即G的损失越小越好：

(10)

式中，损失LG包含内容损失Lcon以及对抗损失Ladv，保证了Itrans与Ifixed具有相似的内容和结构特征。

互信息NMI匹配2个图像之间的联合灰度概率分布，广泛用作多模态可变形配准的代价函数[17]。结构相似性指标SSIM基于灰度值，能够准确量化不同图像之间的对应特征关系。因此内容损失函数表达为：

Lcon=-NMI(Itrans,Ifixed)-SSIM(Itrans,Ifixed)。

(11)

对抗损失约束了配准图像之间的灰度和结构信息，包括正向、反向映射转换的对抗损失，直接引用张桂梅等在文献[11]中的公式：

Ladv=EIfloat～pdata(pfor(Ifloat,Ifixed))2+

EIfixed～pdata(1-pfor(Ifloat,Gfor(Ifixed)))2+

EIfixed～pdata(prev(Ifloat,Ifixed))2+

EIfloat～pdata(1-prev(Ifixed,Grev(Ifloat)))2，

(12)

式中，pfor表示前向映射中浮动图像与参考图像之间相似的概率；prev表示反向映射中浮动图像与参考图像之间相似的概率；E表示相应的期望值。

3 实验分析

3.1 实验数据及配置

Vanderbilt大学主持的医学图像配准评估项目RIRE作为目前较完善的配准项目，提供了大量的病人实例数据，包含了19位患有脑部肿瘤、颅内出血、肺部肿瘤的患者的医学数字影像。

每位病人实例数据中包括一套CT及6套MR数据，本文挑选了480对结构纹理清晰、灰度均匀的CT-MR作为数据集。用该项目提供的测试集进行测试实验硬件及软件环境如表1所示。

表1 实验环境Tab.1 Experimental environment

3.2 实验分析

训练前，将所有样本均预处理为256 pixel×256 pixel的图像。用Adam算法，对判别器和生成器进行交替训练，利用RIRE项目的测试集评估模型性能。设置初始学习率为0.000 1，batch-size设置为16，随着迭代次数的增加，生成器损失函数及判别器损失函数趋势如图5所示。

图5 损失函数趋势Fig.5 Tendency chart of loss function

由图5可以看出，生成损失随着迭代次数的不断增加而减小，最终收敛到0.25左右，并且在整个迭代优化的过程中没有明显的突变，而对抗损失稳定上升，说明本文模型在整个训练过程中稳定可靠。

医学图像包括CT、MR等通过不同灰度等级反映人体组织的密度，因此医学图像的灰度值可以看成是概率测度[18]，当配准效果越好时，对应像素的灰度互信息(NMI)最大，结构相似性(SSIM)最大，偏差(RMSE)越小，表达式分别为：

(13)

(14)

(15)

式中，H(Gtrans)，H(Gfixed)表示配准后的图像与参考图像的香农熵；H(Gtrans,Gfixed)表示2幅图像之间的联合熵；μItrans，μIfixed分别表示图像Itrans，Ifixed的均值；σItrans，σIfixed分别表示图像Itrans，Ifixed的标准差；σItrans，Ifixed表示Itrans和Ifixed的协方差。不同配准模型的配准结果如表2所示。

表2 不同模型图像配准指标Tab.2 Registration indexes of different models

表2对比了基于当下几种常见的配准算法的相关数据，从表中可以看出，本文的算法与传统的Powell相比，NMI提高了30.48%，RANSAC提高了35.41%，SSIM提高了4.35%，RMSE减小了47.93%；与其他基于GANs的配准模型相比NMI提高了5.72%，1.39%，RANSAC提高了2.77%，9.59%，SSIM提高了3.85%，2.28%，RMSE减小了17.57%，5.93%，因此本文所提出的配准模型能够取得较好的配准效果。不同模型配准效果如图6所示。

(a) 参考图像

由图6中可以很明显地看出，基于Powell模型的效果图中高亮部分导致参考CT图像中灰度信息的缺失，而且边缘信息损失严重。cycleGANs和wGAN模型效果图灰度分布不均，软组织以及骨骼组织信息丢失严重，其中cycleGANs头骨边缘模糊，软组织边缘呈锯齿状。而本文模型相较于其他模型不仅在轮廓及拓扑保持上有巨大的优势，而且灰度分布更加均匀，配准更加精准。

为了验证本文模型对形变图像的配准能力，利用Matlab将浮动图像进行不同程度的扭曲变形后输入配准模型，配准效果如图7所示。

(a) 形变图1

由图7可以看出，形变图1形变量小，效果图清晰，细节丰富；形变图2形变量中等，画面模糊，对比度低；形变图3因形变量大亮部信息扭曲严重，配准效果图画面模糊严重。由于医学图像成像设备姿态稳定，信道抗干扰强，成像过程中很难发生大形变，因此本文模型在小形变医学图像配准中能够取得较好的配准效果。

4 结束语

本文构建了基于仿射变换和生成对抗网络的新型配准模型，将传统配准方法与深度学习结合，采用基于随机采样一致算法的仿射变换模型，迭代估计模型参数，不断优化变形场。并在生成对抗网络原有的判别机制上引入了二次判别机制，实现了配准变形场的最优化。通过实验，对比多种配准模型的性能指标，证明了该模型能够取得较好的配准效果。