基于多图谱配准的海马自动分割

胡 昊，张明慧，卢振泰，冯前进，陈武凡

（南方医科大学 生物医学工程学院，广东 广州510515）

0 引 言

从临床研究可知，海马形态体积和功能的改变与许多神经系统疾病有关，如颞叶癫痫和阿尔兹海默病［1］。目前，对于阿尔兹海默病，只能通过早期诊断进而采取有效措施延缓疾病的进程。对海马体积的测量为辨别这种疾病的重要方法［2］。海马体积测量建立在海马分割的基础上，所以准确地从脑部图像中将海马分割出来显得尤为重要。

人工分割海马效率低，主观性强，半自动的海马分割方法可以减少人工分割的时间与负担，但在大样本的研究中，仍不能成为理想的分割方式。自动的图像分割方法［3，4］方为图像分割研究的发展趋势。海马形状不规则，边缘模糊，难以与周围其余组织区分开来，其分割是图像分割研究领域的难点。借助医学图谱提供的先验信息有助于更好地完成分割任务。配准技术有效地将医学先验知识融入分割过程，形成一种新的自动分割技术，称为基于配准的分割［5］。本文提出的基于多图谱配准的自动分割算法，在配准时不但考虑图像对应点的灰度信息，亦融入空间位置信息，实现灰度与空间的同步标准化，提高配准算法的有效性。配准结果的精确度得到保证，融合后的分割结果更接近真实值。

1 基于多图谱配准的海马分割方法

基于多图谱配准的分割方法可以分为基于单图谱、基于平均形状图谱和基于多图谱3种。其中，基于多图谱配准的分割方法准确度更高，适用于对海马进行分割［6］。本文介绍的是一种基于多图谱配准的自动分割算法。该类算法将预标记的图像向目标图像配准，建立对应关系并得到相应空间变换后，利用这个空间变换，将与预标记图像对应的分割结果即图谱，映射到目标图像，再选择合适的融合算法，生成目标图像的分割结果。算法流程如图1所示。

图1 算法流程

图像配准与图谱融合为基于多图谱配准分割算法的两个关键点。本文选择一种新的基于最小化残差复杂度（modified minimization of residual complexity）的配准方法［7］，并使用STAPLE（simultaneous truth and performance level estimation）算法对图谱进行融合［8］，得到最终分割结果。

1.1 最小化残差复杂度的图像配准算法

对磁共振图像而言，由于磁场的不均匀性，造成图像之间存在亮度不一致等问题，常用的配准方法无法达到较好的配准效果。最小化残差复杂度算法（minimization of residual complexity）适合处理空间灰度不均匀的问题。假设I与J 是待配准的两幅图像，则有

式中：I——参考图像，J——浮动图像，T——对浮动图像J 进行的形变，S——灰度校正场，η——零均值的高斯噪声。经过推导和演化，得能量函数如下

式中：N——体素个数，α——可调节的参数，qnT——对灰度校正场的约束演化而来的矩阵的特征向量。详细推导过程请参考文献［7］。

该方法对亮度不一致、明暗不均匀的图像配准有很好的效果。然而，由于磁共振扫描时间较长，头部运动增加，扫描结果易出现噪声和失真。此外，磁共振磁场本身存在不均匀性，不同仪器、不同时间、不同病人的扫描，相互之间存在较多偶然因素。噪声及这些偶然因素造成的突变值被融入残差r中，从而影响式（1）的能量函数，导致配准精度下降，鲁棒性较差。对于基于多图谱配准的分割方法，只有配准精度更高，最终的分割结果才能更精确。因此，本文对该方法进行改进，构造一个新的配准算法。

1.2 改进后的最小化残差复杂度的图像配准算法

式（1）中的r为参考图像与浮动图像的灰度值差，这种完全依赖图像灰度信息的配准方法没有最大化地利用图像提供的所有信息。改进后的新算法在考虑图像对应点的灰度信息的同时，融入空间位置信息。对不同像素点，根据其方向与空间信息的不同，赋予不同权重，进行空间标准化。权重参数计算步骤如图2所示。

图2 权重参数计算流程

对二维图像，邻域选择为8 邻域系统，如图3 所示。三维图像则为26邻域。以三维图像为例描述权重参数求解方法如下。

图3 8邻域系统

梯度图可以在一定程度上反映原图像灰度、位置的变化，如图4（b）所示，在各器官组织边缘部分，梯度值较大，而灰度均匀的同一组织内，梯度值则很小。由于磁共振扫描的特性，在两幅图像中，可能出现同样的组织有不同的灰度值，而相同的灰度值代表不同组织的情况。对梯度图而言，即使相同组织灰度值不同，梯度信息同样可以准确地反映该组织的位置结构信息。因此，权重参数ω的提出，为r融入了图像的空间位置信息，完成配准图像灰度与空间的同步标准化。这种方法可以有效降低噪声或其它偶然因素造成的图像灰度值突变的影响，具有更高的鲁棒性与配准精度。

图4 脑部MR 图像及对应梯度

两种配准方法均采用自由形变模型（free form transformation）对图像进行形变。自由形变模型的方法因为其逼真的生理结构模拟特性，被广泛地用于图像配准中［9］。此外，为避免陷入局部极值，均采用为由粗到细的多分辨率网格技术。将多分辨率网格与自由形变模型相结合，可以提高计算效率，保证配准结果的准确度。对形变场迭代求解时则选择梯度下降法为优化算法［10］。该方法工作量小，需要存储的变量少，对初始点的要求较低。

1.3 STAPLE融合算法

本文采用STAPLE算法实现图谱融合。STAPLE 算法的目标为利用函数f（D，Tr｜p，q）的似然函数估计性能水平参数（p，q）。f（D，Tr｜p，q）的最大化就实现了参数（p，q）的最优值。公式表示为

式中：D——N×R 的矩阵，表示图谱的每个像素在每个分割中的结果，R——图谱数量，N——每个图谱的体素个数，Tr——目标图像的真实分割结果。对每个像素用二进制0或1表示，其中1指示分割对象，0指示非分割对象。f（D，Tr｜p，q）表示概率群分布函数。其中pj表示敏感度，qj表示特异性

式中：pj和qj的取值范围为［0，1］。假设真实的分割结果、pj和qj、图谱的分割结果是条件独立的，即不同分割结果之间相互独立，没有影响，那么分割质量就可以用pj和qj来评价。

对式（3）可采用最大期望算法估计其值，详细步骤可参照文献［11，12］。

2 配准及分割结果的评价方法

选择合适的评价方法，可以帮助定量、科学、有效地对实验结果进行分析，并对不同实验方法进行深入比较，通过正确整合和严格评价数据来提取有用数据支持结论。

2.1 配准结果的评价

本文选择互信息作为配准结果的衡量标准。互信息是一种度量两个系统的统计相关性的量，其计算复杂度低，鲁棒性好。计算公式如下式

式中：PS、PA’——参考图像S和形变后的浮动图像A’的边缘概率密度分布函数，PA’S——联合概率密度分布函数。计算的互信息值越大，表示两幅图像的空间位置越一致，配准效果越好［13］。

假设用NRC方法配准后得到的图像计算的互信息值为mrcMI，用RC方法配准后得到的图像计算的互信息值为rcMI，由下式计算比较两种方法得到的互信息值

2.2 分割结果的评价

相似性测度Dice为分割结果的常用评价方法。相似性测度是对目标图像的真实分割结果和融合得到的分割结果的重叠率进行评价，计算公式如下

式中：Sseg——目标图像的真实分割结果，Sfus——融合得到的分割结果，V——体积大小。相似性测度Dice值越大，说明目标图像的真实分割结果和融合得到的分割结果的重叠率越高，结果越精确［14］。

3 实验结果及分析

磁共振图像具有非损伤性的优点，对软组织有极好的分辨力，被广泛用于临床医学研究。本文选择的18 个T1加权的磁共振脑部图像均来自美国的Henry Ford医院。18个图像分别由两台不同场强的磁共振扫描仪获得，部分病人患有颞叶癫痫症，出现海马体萎缩的现象，其余海马体的形态功能正常。

为增加实验结果的可信度，避免出现偶然性，随机抽取10个图像分别作为待分割图像，编号为01到10号，剩余8个图像编号为11 至18 号。用Matlab7.12.0在PC 机上进行了实验，并分别用最小化残差复杂度和改进后的最小化残差复杂度两种不同的配准方法进行对比。两种不同配准方法最终分别得到10个融合结果。

3.1 图像预处理

18个图像分别由两台不同场强的磁共振扫描仪获得，图像的分辨率存在较大差异，配准前需进行预处理。处理方法为以图01为基准，将其余17个图像归一化到基准图像的空间上，使它们有一致的原点与尺寸。此外，图像的对比度也存在一定差别，需对图像的灰度值进行调节，使得18个图像的灰度值和对比度处于同一水平。

本算法是对目标图像的海马进行分割，配准部位不必涵盖整个脑部。因此，为提高配准效率，降低算法所需时间，以海马为中心进行ROI提取。处理方法为减小图像的体积，提取一个以海马为中心的长方体。提取后的图像大小均为64×128×64。

图5为01 号图像分别在预处理前及预处理后，横断面、冠状面和矢状面的二维图像，其中第一行为预处理前，第二行为预处理后。箭头指示区域为图像对应的海马。

图5 预处理前及预处理后的01图像

经过预处理后图像以海马为中心，去除与分割目标无关的大部分组织，能更好地对海马区域进行配准，大大降低算法所需时间。

3.2 基于最小化残差复杂度算法的配准

为进行对照实验，采用了最小化残差复杂度和改进后的最小化残差复杂度两种方法进行配准。其中，两种配准方法的参数选择如下：网格大小为8，参数α为：0.01，参数λ为：0.001，其中λ为式（1）中约束项qnT的参数。

图6为式（4）计算结果的条形图，横坐标表示不同目标图像的十组配准结果，每组分别计算得到17个不同值，纵坐标表示MIdif 值。

图6 MIdif 值条形图

由式（5）可知，MIdif 值若为正数，表示由式（4）计算的mrcMI值大于rcMI 值。从图6可观察到，MIdif 值多为正数，说明MRC方法配准得到的图像计算的互信息值比RC方法配准得到的图像计算的互信息值要大，配准效果更好。

配准完成后，每组分别得到17个形变场，利用该形变场将已分割的海马映射到目标图像上，即可得到17个形变后的海马。图8为07号图像作为待分割目标图像时，映射后得到的目标图像的17个分割结果。其中第一幅为07号图像对应的真实海马，后17幅为映射形变后的海马。图7只显示横断面。

图7 17号图像的真实分割图及映射后的各海马

观察图7可发现，17个初步分割结果在形状位置上都与真实值有一定差异，若选择其中一个图像进行单图谱配准分割，得到的结果必不能作为有效的海马分割图。因此，本文采用多图谱配准分割方法，对17个初步分割结果进行融合，生成一个与真实值更为接近、更精确的海马图。

3.3 基于STAPLE算法的图谱融合

本文采用的融合方法为STAPLE 算法。参数选择如下：收敛因子为0.00001，初始化类型为0，插值类型选择邻域插值法。其中初始化类型的选择与最大期望算法的初始估计值有关。对两种不同配准方法分别得到的十组结果融合后，最终分割的海马如图8所示。其中第一行与第四行为目标图像的真实分割结果，第二行与第五行为用RC方法进行配准融合后得到的结果，第三行与第六行为用MRC方法进行配准融合后得到的结果。第一列第一行至第三行为第一组分割结果对比，第二列第一行至第三行为第二组分割结果对比，以此类推。

由图8可观察到，MRC配准方法融合得到的海马，在形状、大小上都与真实分割结果更为接近。说明基于该配准方法将图谱映射到目标图像上得到的结果与RC 配准方法相比更精确，更接近目标图像。

通过图8可对分割结果进行主观定性分析，但定量分析的方法更为科学，在对实验结果的评价中更为重要。通过式（6）计算两种不同配准方法最后融合的结果与待分割的目标图像的真实分割数据的Dice值，计算结果见表1。

图8 目标图像的真实分割海马与两种不同配准方法融合结果

表1 两种不同配准方法融合结果Dice值

对比表1数据可以发现，用RC方法配准后融合的结果的Dice值有两组大于0.8，而用MRC方法配准后融合的结果的Dice值均在0.8以上。改进后的配准方法明显地提高了海马分割的精度。

为更直观地对表1数据进行对比分析，分别计算两种方法Dice值的均值并作折线，如图9所示。

图9 两种不同配准方法融合结果Dice值折线

观察图9可以看出，MRC配准方法融合得到的十组分割结果Dice均值远大于RC 配准方法融合得到的结果，实心圆点折线位于空心圆点折线之上，且更平缓，即分割结果更精确、更稳定。由此说明MRC 配准方法的精确度更高，本文选用的基于MRC的配准方法的多图谱配准分割算法是一种有效的海马分割算法。

4 结束语

本文提出一种基于多图谱配准的海马分割算法，该算法在配准时融入图像的空间与方向信息，增加配准算法的鲁棒性与有效性，从而实现更精确的海马分割。从实验结果可以观察到，该分割算法更稳定，对海马分割的精度更高。由此表明本文提出的基于NRC配准方法的多图谱配准分割算法是一种有效的海马分割算法。分析实验结果Dice值可知，分割结果仍有一定的提高空间。后续工作中，可增加图谱数量并进行择优筛选，或采用多种不同的融合算法，进行对比研究，找到最适合进行海马分割的基于多图谱配准的分割方法。

［1］Yushkevich Paul A，Wang Hongzhi，Pluta Jogn，et al.Nearly automatic segmentation of hippocampal subfields in vivo focal T2－weighted MRI ［J］.Neuroimage，2010，53 （4）：1208－1224.

［2］Morra Jonathan H，Tu Zhuowen，Apostolova Liana G，et al.Validation of a fully automated 3D hippocampal segmentation method using subjects with Alzheimer’s disease，mild cognitive impairment，and elderly controls ［J］.Neuroimage，2008，43（1）：59－68.

［3］Chupin Marie，Gérardin Emilie，Cuingnet Rémi，et al.Fully automatic hippocampus segmentation and classification in Alzheimer’s disease and mild cognitive impairment applied on data from ADNI［J］.Hippocampus，2009，19 （6）：579－587.

［4］Leemput Koen Van，Bakkour Akram，Benner Thomas，et al.Automated segmentation of hippocampal subfields from Ultrahigh resolution in vivo MRI ［J］.Hippocampus，2009，19（6）：549－557.

［5］Chen Wenyan，Li Shutao，Jia Fucang，et al.Segmentation of hippocampus based on ROI atlas registration ［R］.International Symposium on IT in Medicine and Education，2011：226－230.

［6］Shan Liang，Charles Cecil，Niethammer Marc.Automatic multi－atlas－based cartilage segmentation from knee MR images［R］.9th IEEE International Symposium on Biomedical Imaging，2012：1028－1031.

［7］Myronenko Andriy，Song Xubo.Image registration by minimization of residual complexity ［R］.IEEE Conference on Computer Vision and Pattern Recognition，2009：49－56.

［8］Andrew J A，Bennett A L.Non－local STAPLE：an intensitydriven multi－atlas rater model［J］.Med Image Comput Comput Assist Interv，2012，15 （Pt 3）：426－434.

［9］Rueckert Daniel，Aljabar Paul，Heckemann Rolf A，et al.Diffeomorphic registration using B－splines ［J］.Med Image Comput Comput Assist Interv，2006，9 （2）：702－709.

［10］LI Bin，OU Shanxing，TIAN Lianfang，et al.Thorax multimodal medical image registration based on adaptive free－form deformation and gradient descent［J］.Application Research of Computers，2009，26 （10）：3978－3982 （in Chinese）. ［李彬，欧陕兴，田联房，等.基于自适应自由变形法和梯度下降法的胸部多模医学图像配准 ［J］.计算机应用研究，2009，26 （10）：3978－3982.］

［11］Wikipedia.Expectation maximization，http：／／en.wikipedia.org／wiki／Expectation＿maximization［DB／OL］.2013－7－20（in Chinese）.［维基百科.最大期望算法，http：／／zh.wikipedia.org／wiki／最大期望算法［DB／OL］.2013－7－20.］

［12］CHEN Wenyan，LI Shutao.Segmentation of hippocampus based on ROI atlas registration ［D］.Changsha：Hunan University，2012 （in Chinese）.［陈雯艳，李树涛.基于ROI多图谱配准的海马磁共振图像分割 ［D］.长沙：湖南大学，2012.］

［13］SHI Yingqi，GU Lixu.Optimizied multistage medical image registration method based on mutual information ［J］.Computer Engineering，2006，32 （22）：187－190 （in Chinese）.［施颖琦，顾力栩.基于互信息多步骤优化的医学图像配准［J］.计算机工程，2006，32 （22）：187－190.］

［14］Han Xiao，Fischl Bruce.Atlas renormalization for improved brain MR image segmentation across scanner platforms ［J］.IEEE Transactions on Medical imaging，2007，26 （4）：479－486.