张 纠 刘晓芳* 杨 兵

1(中国计量大学计算机应用与技术研究所 浙江 杭州 310018) 2(中国计量大学电子信息与通信研究所 浙江 杭州 310018) 3(浙江省电磁波信息技术与计量检测重点实验室 浙江 杭州 310018)

0 引 言

基于计算机辅助技术的医学影像分析为病灶识别与定位[1]、疾病诊断[2-3]等提供了可靠的信息基础。其中，医学图像配准是指对浮动图像(Moving image)进行空间变换，确定浮动图像与固定图像(Fixed image)的空间对应关系[4]，寻求最佳空间映射。

医学图像配准问题的核心在于如何确定图像间的相似性以及如何对图像进行空间映射[4]，在临床应用中，通常将同一个人在不同时间段所获得的图像进行配准或将多个人配准到标准图谱上。传统基于灰度值的配准方法[5]将配准问题抽象为最小化相似度测量以寻找最优参数的优化问题，然而，此种方法存在效果较差，无法取得最优解等局限性问题。近年来，由于神经网络具有能够从图像特征中学习复杂映射、自动优化配准参数、建立优化模型等优点，提供了用于处理配准问题的新思路，已被广泛用于医学图像配准中[6-8]。

基于编码-解码结构的卷积神经网络如U-net[9]能够自动对输入图像进行底层特征提取，其独特的网络结构在处理医学图像解剖结构相对固定、数据量较少的情况下具有较好的拟合性能。空间变换网络(Spatial Transform Networks, STN)[10]是一种结合卷积运算、仿射变换、空间采样的显式特征空间变换网络结构，其常用于特征提取网络之后，对配准场(Registration Field)进行空间变换得到配准图像(Registered Image)，以确定图像间的空间映射关系。基于空间变换网络的编码-解码结构是处理医学图像配准的常用网络框架，然而，此类网络具有没有充分考虑到各通道特征重要性不同、基于单路的特征空间变换对局部图像特征变换能力不足，以及特征与特征之间信息融合程度不高等问题，针对此问题，本文提出一种基于特征块局部空间变换的双路编码-解码网络用于配准问题，本文方法在图像特征块的基础上，进行了特征重要性加权和基于双通道的局部特征空间变换。

1 相关工作

1.1 医学图像配准

医学图像配准是后续图像融合，疾病诊断的基础，其目标是根据浮动图像与固定图像之间的相似性关系，确定两幅图像间的最佳映射。根据空间变换的类型，医学图像配准方法可分为刚体配准与非刚体配准[11]。刚体配准方法具有变换参数较少、运算效率高等优点，然而，大多数医学图像如心脏CT(Computed Tomography, CT)图像都存在非刚体变换，因此对于基于非刚体变换的配准方法研究较多。

传统基于灰度值的配准方法采用互相关系数(Cross Correlation)[12]等相似性测量指标，通过优化参数模型确定最终映射关系，此外，互信息(Mutual Information)[13]及其改进的条件互信息(conditional Mutual Information)[14]、最大梯度距离(Maximum Distance-Gradient, MDG)[15]、熵图像(Entropy Images)、流行学习(Manifold Learning)[16]等也常用于医学图像配准问题中。

神经网络本质上也是一种参数优化模型，其优点是能够通过最小化损失函数自动学习优化参数，建立图像间的配准模型。近年来，神经网络在处理配准问题上表现出较好的学习性能，如递归级联网络[17]通过网络级联的方式逐层对图像特征进行空间变换，在肝脏配准问题上表现出较好的泛化性能，递归配准网络[18]将递归神经网络引入图像配准问题中，将图像看作一系列序列模型。

医学图像配准常用的网络框架如图1所示，采用基础特征提取网络得到配准场，通过最小化损失函数得到浮动图像与固定图像间的空间映射参数。

图1 医学图像配准框架

1.2 基于空间变换网络的编码-解码配准模型

基于编码-解码结构的网络模型如U-net[9]、SegNet[19]等，利用跳跃连接(Skip Connection)结构融合同级之间不同分辨率的特征，利用高分辨率特征具有轮廓、梯度、颜色等信息，低分辨率特征具有抽象化特征信息，将两者的语义信息融合，使融合后的特征同时具有高层与底层特征的语义信息，能够对原始输入图像做出精细化特征表达。

图2 编码-解码网络结构

图2为编码-解码网络基本结构，其整体可分为三个组成部分：特征编码，特征解码，特征融合。一般来说，输入图像与输出图像具有相同或相近大小分辨率，输入图像经过卷积、激活函数、下采样等运算后得到编码后的图像特征，下采样的次数一般为3到4次。过多的下采样次数导致参数数量增加，计算速度减慢；过少的下采样次数导致提取到的特征不足以表征图像的深层细节。最大池化(max-pooling)是常见的下采样运算，可看作正则化的一种方式。卷积运算是实现图像特征提取的关键，其数学表达式如下：

(1)

式中:x为输入图像；k为卷积核；z为输出图像。基于编码-解码结构的卷积网络通常采用跳跃连接的方式进行特征融合，其数学表达如下：

(2)

式中:X表示编码层级的特征；Z表示解码层级的特征；O表示通道拼接后的融合特征。

图3 空间变换网络基本结构

图3为空间变换网络结构，空间变换网络是一种针对图像特征变换的通用网络模块，能对输入图像进行仿射变换，由于其采用卷积、全连接等运算，使得空间变换网络可以拟合复杂函数映射，输出变换参数，得到空间变换后的图像特征。空间变换网络整体由定位器、网格生成器、采样器三部分组成。

1) 定位器通过卷积等运算，输出输入特征与输出特征之间的变换参数：

(3)

2) 网格生成器根据仿射变换矩阵以及输出特征空间坐标，计算与之对应的输入特征空间坐标。

3) 采样器根据网格生成器得到的空间坐标点，利用双线性插值对输出特征进行填充：

(4)

2 基于特征块空间变换的配准模型

2.1 网络结构

基于空间变换网络的编码-解码网络没有考虑到不同通道的特征重要性不同，此外，空间变换网络针对整个输入特征，局部特征空间变换能力不足。针对此问题，本文提出了一种基于特征块的双通道空间变换网络(Feature Patch Based Dual-channel Spatial Transform Network, FPDC-STN)。结合编码-解码网络，本文的配准模型结构如图4所示。

图4 基于特征块的双通道空间变换网络配准模型

本文模型采用两个编码-解码结构的卷积神经网络分别对浮动图像和固定图像进行特征提取，同时将固定图像的特征编码部分与浮动图像的特征解码部分对齐，组成对称网络。为提高特征利用率，借鉴特征金字塔网络[20]思想，从固定图像和浮动图像的编码-解码网络中分别抽取了第2、第4、第6层特征，两者特征融合之前用二维卷积进行了特征空间转换。图4中，FW-Block为特征加权(Feature Weighting)模块，FPDC-STN为本文提出的特征块双通道变换模块，两者的网络结构如图5、图6所示。

特征加权模块对原始输入特征的通道加权可以分为4步：1) 对于输入特征，采用两个不同尺寸的分组卷积核得到两组不同的输出特征，并设置分组数K对两组输出特征进行分组，分别得到K个特征块；2) 利用式(5)对特征块进行全局平均池化(Global Average Pooling, GAP)，得到特征加权参数；3) 使用全连接和激活函数等对特征块进行通道加权，如式(6)所示；4) 分别对K个特征块进行合并，采用通道拼接的方式得到原始输入特征通道加权后的输出特征。

(5)

式中：H表示输入特征块的高；W表示输入特征块的宽；UC表示某个通道的特征块；ZC表示某个通道的输出特征块。

O=σ(g(Z,w))=σ(w2δ(w1Z))

(6)

式中：Z表示输出特征块；w1、w2表示全连接参数；δ表示非线性激活函数如ReLU等；σ表示Sigmoid函数。

借鉴空间变换网络结构的思想，本文采用了基于特征块空间变换的方式对原始输入特征进行空间变换。此外，针对原始空间变换网络感受野大小固定的局限性问题，本文还提出了一种基于双通道的空间变换网络，基于权值共享方式将特征块空间变换与双通道策略结合形成最终的基于特征块的双通道空间变换网络(FPDC-STN)。FPDC-STN对输入特征的空间变换可分为以下3步：1) 采用两个不同尺寸大小的卷积核(3×3，5×5)对原始输入特征进行分组卷积，然后对特征进行分块；2) 对每个特征块，分别进行空间定位-网格生成-采样操作，得到空间变换后的输出特征块；3) 对输出特征块进行合并，基于通道拼接的方式得到最终空间变换后的输出特征。

2.2 损失函数

本文所用的损失函数包括两部分，分别为衡量浮动图像与固定图像之间的相似性指标Lsin(φ)，配准场的正则化项Lsmooth(φ)。两者分别定义如下:

Lsim(F,R,φ)=-CC(F,R(φ))=

(7)

(8)

式中:CC(F,M(φ))表示局部相关系数(Correlation Coefficient)，衡量配准图像R(φ)与固定图像F之间的局部相关性；φ表示配准场(Registration Field);pi表示索引为i的像素。总损失函数定义为相似性指标Lsim(φ)与正则化项Lsmooth(φ)之和，如式(9)所示，λ为正则化参数。

L(F,M,φ)=Lsin(F,M,φ)+λLsmooth(φ)

(9)

2.3 整体算法流程

本文提出的基于局部空间变换网络的医学图像配准算法整体流程如图7所示。

图7 本文整体算法流程

Step1输入数据集划分参数train_rate、validate_rate、test_rate，将数据集划分为训练集、验证集、测试集。

Step2对训练数据进行数据增强，包括图像旋转、图像扭曲等。

Step3对待配准图像进行预配准、去噪等预处理操作。

Step4初始化网络参数、迭代次数Itr、正则化系数lamda等。

Step5更新网络参数，以及分别建立特征块双通道变换模块、特征加权模块、特征编码-解码网络。特征块双通道变换模块(其中split_num表示划分的特征块数量)主要伪代码表示如下(使用Pytorch框架编程)：

Apply 3×3 convolution and 5×5 convolution to input features;

Split input features into split_num;

# 空间定位

Location_module=nn.Squential(nn.Conv2d(1,8,3),

nn.MaxPool2d(2,2),

nn.Relu(True),

nn.Conv2d(8,10,5),

nn.MaxPool2d(2,2),

nn.ReLU(True))

#3*2 仿射变换回归器

fc_loc=nn.Sequential(

nn.Linear(10 * 3 * 3, 32),

nn.ReLU(True),

nn.Linear(32, 3 * 2)

)

# 组合特征块

Combine input features

特征加权模块表示如下；

# 划分特征块

Split input features into split_num;

# 平均池化

avg_pool = nn.AdaptiveAvgPool2d(1)

# 特征加权以及特征赋值

fc=nn.Sequential(

nn.Linear(channel, channel // reduc tion, bias=False),

nn.ReLU(inplace=True),

nn.Linear(channel // reduction, channel, bias=False),

nn.Sigmoid()

)

# 组合特征块

Combine input features

Step6网络前向传播获得预测值。

Step7用式(9)计算损失函数，并计算配准评价指标。

Step8使用梯度下降算法，如式(10)所示，网络反向传播。

(10)

Step9迭代次数或损失函数是否满足要求，若满足要求，则执行Step 10，若不满足要求，则执行Step 5。

Step10配准算法结束。

3 实 验

实验采用的硬件设备是显存大小为12 GB的Nvidia GeForce GT1080ti GPU以及Intel i7 8700 CPU，软件环境采用Python 3.6，机器学习处理包sklearn以及医学图像处理工具SimpleITK，配准实验基于深度学习框架Pytorch[27]实现。初始学习率设置为0.000 01，并使用初始动量为0.9的Adam优化器。迭代次数设置为2万次，并采用随机梯度下降(Stochastic Gradient Descent，SGD)更新网络参数，正则化系数λ设置为1.0。

3.1 数据集与实验参数

本文分别在脑部数据以及肝脏数据上进行了配准实验。肝脏数据集来自MSD[21]以及SLIVER[22]，脑部数据集来自ADNI[23]、ABIDE[24]、LBPA[25]。所有的数据重采样为132×132×132分辨率大小，并用医学图像处理工具FreeSurfer[26]进行了包括初步配准、脑部区域提取等标准预处理步骤。提取肝脏扫描数据、脑部扫描数据总计1 000个，其中800个用于训练配准模型，100个用于验证模型，100个用于测试模型。此外，为了增加实验数据，所有数据都进行一系列数据增强操作，如二维旋转、中心裁剪、数据扭曲等。实验样例图像如图8所示。

图8 实验样例图像

3.2 评价指标

本文采用Dice系数评价配准与固定图像之间的重合程度：

(11)

式中:Dice系数范围为[0,1]。当Dice系数为1时表示配准图像与固定图像完全匹配，当Dice系数为0时表示配准图像与固定图像完全不匹配。

3.3 脑图像配准实验

本文评估了不同实验参数对配准效果的影响，评价标准为验证集Dice系数、一段时间后评价指标是否变化、网络收敛时间。评估参数包括正则化系数λ、迭代次数、Adam优化器初始动量、学习率。结果如表1-表4所示。

表1 迭代次数对配准精度的影响

表2 正则化系数对配准精度的影响

表3 Adam初始动量对配准精度的影响

表4 学习率对配准精度的影响

根据表1可知，迭代次数为20 000次时，一段时间后验证集评价不再变化，故本文设置的迭代次数为20 000次。根据表2可知，正则化系数λ=1时，验证指标Dice系数最高，故本文设置的正则化系数为1。根据表3可知，Adam初始动量为0.90时，验证集Dice系数指标最高，故本文设置的Adam初始动量为0.90。根据表4可知，学习率为0.000 01时，配准网络已经收敛，并且验证集Dice系数指标最高，当学习率为0.000 1时，虽然配准网络已经收敛，但是其验证集Dice系数低于学习率为0.000 01时Dice系数，故本文设置的学习率为0.000 01。

图9为本文配准方法对于脑部图像的配准结果，此外，为了说明本文配准模型的有效性与准确性，与UtilzReg方法[28]以及FAIM方法[29]进行了对比。

图9 脑部图像配准结果

本文选取了两个切片的配准结果如图9第一行以及第二行所示，图9中第三行是第二行配准结果的局部放大图。从固定图像与浮动图像中可以看出，固定图像的外形轮廓与浮动图像有较大差别，两者内部组织结构也不相同，固定图像内部组织如脑白质、脑灰质等纹理较为简单，而浮动图像内部组织的纹理比较复杂，这对配准算法的空间变换性能提出了较高要求。

图9(c)为本文提出的配准方法的配准结果，本文方法可以针对浮动图像与固定图像之间的差异，对浮动图像的外部轮廓以及内部组织纹理进行准确的空间变换，这表明本文提出的基于特征块的双通道空间变换网络FPDC-STN以及特征通道加权能够捕捉图像细节特征，对图像局部做出准确配准。图9(d)为FAIM配准方法的配准结果，可以看出，FAIM方法能够根据固定图像与浮动图像之间的整体轮廓差异对浮动图像做出较为准确的空间变换。然而，相比本文配准方法，FAIM方法针对图像局部细节处理仍显不足，配准结果中局部细节与浮动图像差异不大，这说明FAIM方法不能准确地根据局部图像细节对浮动图像进行空间形变。图9(e)为UtilzReg方法的配准结果，整体来说，此方法与FAIM方法差距不大，在处理图像配准细节上较差。

表5为本文方法与FAIM方法以及UtilzReg方法的评价结果，本文进行了10次重复实验，对配准结果计算Dice指标的平均值以及标准差。可以看出，本文方法所取得的平均Dice指标相比FAIM方法提高2.7%，相比UtilzReg方法提高1.4%，在几种方法中取得最高的Dice系数。从配准稳定程度上看，本文的Dice指标的标准差最小，也是三种配准方法中最稳定的算法。

表5 各方法脑部配准评价

3.4 肝脏图像配准实验

图10为本文在肝脏图像上的配准结果，同时也与FAIM算法以及UtilzReg算法比较了肝脏图像配准的有效性与准确性，结果如表6所示。

图10 肝脏图像配准结果

表6 各方法肝脏配准评价

图10第一行与第二行为选取的部分切片图像配准结果，第三行为第二行配准结果的局部放大图。可以看出，本文配准方法能够较准确地处理图像细节，如图10(c)所示，针对固定图像中的局部亮点，本文配准方法能够对浮动图像进行局部空间变换，配准结果与固定图像差异不大。对于FAIM以及UtilzReg方法，两者能较好地对肝脏整体轮廓作出形变，但是在处理局部细节时仍然存在较多孤立像素点以及局部空间变换不足等问题。表6表明，与FAIM配准方法相比，本文方法平均Dice系数提高了8.5%，与UtilzReg方法相比，本文方法提高了3.6%。本文方法在保持最高平均Dice评教指标的同时，具有较好的稳定性。

4 结 语

本文在现有基于空间变换网络的编码-解码配准模型基础上，提出了针对特征融合的特征加权模块以及基于特征块的双通道空间变换网络FPDC-STN；针对空间变换网络无法有效解决图像局部空间变换、感受野有限，以及空间变换参数无法共享等问题，本文提出的基于特征块的双通道空间变换网络能够针对拆分的特征块进行准确局部变换；同时，利用双通道策略以及两种尺寸大小不同的分组卷积策略，共享空间变换网络，有利于减少网络参数量和提高网络的空间变换能力。此外，本文提出的特征加权模块以及特征块变换模块通过两个大小不同的卷积核，获取不同大小的感受野，本文提出的配准算法可以对自然图像和医学图像进行特征变换和特征加权，作为具体应用实例，本文在医学图像上验证了其配准有效性和准确性，进一步提升了图像的空间变换性能，提高了图像的配准精度。