一种稳健的非刚性医学影像配准方法

李 琦,姬红兵,臧 博,刘 靳

(西安电子科技大学电子工程学院,陕西西安 710071)

一种稳健的非刚性医学影像配准方法

李 琦,姬红兵,臧 博,刘 靳

(西安电子科技大学电子工程学院,陕西西安 710071)

提出了一种基于局部线性嵌入(LLE)和改进的L-BFGS优化算法的非刚性配准新方法.该方法首先计算影像不同方向上的定序特征,用于补充传统互信息测度中缺失的空间结构信息;然后,运用LLE方法及其逆映射对高维定序特征进行降维和融合;进而结合影像灰度信息构造了一种基于混合熵的配准测度,有效保证了配准测度函数的光滑性和收敛性;最后,采用改进的L-BFGS优化方法搜索最优配准参数.多组仿真数据的测试结果表明,在噪声情况下,所提方法具有精度高、鲁棒性强的特点,优于现有几种方法.

非刚性配准;局部线性嵌入;定序特征

医学影像配准是医学影像处理领域的一个重要的和基本的研究课题,是医学影像融合、医学影像重建、影像与标准图谱的匹配等研究的基础[1-2].医学影像配准就是寻求一种最优空间变换,使两幅医学影像的对应点在给定相似性测度下达到空间位置和解剖位置的一致.医学影像的配准方法可以分为刚性配准和非刚性配准两大类.经过多年的发展,刚性配准算法已经成熟,可以达到很高的配准精度.然而,医学影像配准中刚性体的配准只是很小的一部分,许多重要的临床应用需要非刚性变换来描述影像之间的空间关系.

另外,因为噪声环境的复杂性,含噪医学影像的配准也受到很大关注.在医学影像的成像过程中,通常将成像必需的但对人体有害的强磁场、示踪剂和放射线等控制在一定范围内,加上成像模式本身的一些物理限制,往往导致影像模糊,并伴有噪声甚至伪影的出现.因此,在医学影像配准方法研究中,分析相似性测度对噪声的稳健性也是研究配准性能的一个重要方面.

当前,非刚性配准方法主要分为基于特征的方法和基于灰度的方法[1-4].基于特征的配准方法一般算法效率较高,但算法精度通常取决于特征提取的精度[3].基于灰度的配准度量影像灰度间的信息分享程度,具有较高的精度和鲁棒性,但计算量较大[4].在此类方法中,互信息由于不需要获取待配准影像的先验知识、无需预处理及人工干预,适用性强,得到了高度关注和广泛应用.然而,由于Shannon熵假设相邻像素灰度值不相关,互信息忽略了像素灰度分布的空间位置关系,导致了相似性测度存在过多局部极值、对噪声敏感等问题,最终影响了配准精度.当前的研究热点是建立将影像空间结构特征与互信息相结合的混合测度[5-9].表示空间结构信息的特征描述符很多,例如影像梯度、尺度不变特征变换(SIFT)、灰度共现矩阵(GLCM)以及定序特征等.Pluim等[6]将一个梯度相似性系数与归一化互信息相乘作为最终的配准测度,但影像梯度对噪声十分敏感,无法有效增强配准相似性测度的抗噪能力.SIFT特征对影像旋转、尺度缩放、亮度变化等保持不变性,但在其算法的关键步骤,利用关键点邻域像素的梯度方向分布特性,为每个关键点赋予方向;同时, SIFT丢弃低对比度的关键点以致对边缘模糊的目标无法准确提取特征点.因此,SIFT不适用于低分辨率的医学影像的配准.Rueckert等[7]提出了二阶互信息的概念,将相邻像素对的GLCM作为概率密度函数,扩展了边缘熵及联合熵的定义.然而,由于GLCM是像素距离和角度的矩阵函数,因此完整的GLCM计算,其参数的选取范围很广,计算量很大.文献[9]将影像定序特征引入到医学影像配准算法中,作为一种新的图像特征,定序特征包含了丰富的纹理、边缘等空间结构信息,具有很强的抗噪能力和通用性,同时不易受到影像亮度、对比度等因素的影响[10-11],能够有效补充互信息测度的不足.

医学影像多方向的定序特征是高维的,包含了很多不重要的冗余信息,文献[9]中采用主成分分析(PCA)方法进行降维处理.但是由于高维定序特征本质上是非线性的,PCA方法很难发掘出其中非线性分布的几何结构和相关性.局部线性嵌入(LLE)是一种无监督的非线性学习方法,以保存数据局部邻域间相互关系的方式把高维数据映射到一个低维全局坐标系下,能够揭示数据的全局非线性结构[12-13].因此,笔者采用LLE及其逆映射对多方向的定序特征进行降维和融合,并在此基础上结合影像灰度信息构造了一种基于混合熵的配准测度[9],利用定序特征对噪声的强鲁棒性有效保证了配准测度函数的光滑性和收敛性,大大提高了配准精度.

1 空间结构特征

使用一组和差定序滤波器[11](相位分别为θ=0,π/12,π/6,π/4,π/3,5π/12)对一幅T2-MR影像进行6个方向滤波,如图1所示.图1(a)为待滤波影像,图1(b)～(g)分别为6个方向的滤波结果.可以看到,一幅医学影像多个方向上的定序特征所包含的信息不尽相同,但也存在极强的相关性.为了最小化不同方向的定序特征中所包含的冗余信息,同时降低配准算法的计算代价,文献[9]采用PCA方法进行降维处理.但是由于高维定序特征是影像多方向的滤波结果,因此本质上是非线性的,而PCA方法作为一种线性的数据降维方法并不能准确揭示出非线性数据的内在几何结构和相关性.若采用PCA方法处理非线性数据,很有可能会丢失一些重要的信息,从而影响最终的配准精度.因此,文中采用更为适合的LLE方法将高维观测空间中的多方向定序特征映射为低维嵌入空间中的对应点集.接下来,如何通过这些低维嵌入空间中的对应点集获得高维观测空间中医学影像的空间结构特征是一个关键问题.众所周知,样本的聚类中心通常最具代表性,因此与低维嵌入空间中的中心点对应的高维观测空间中的影像也应该最具代表性,即能够描述医学影像整体空间结构特征的多方向定序特征的融合.文献[13]在LLE方法的基础上提出一种从低维嵌入空间向高维空间映射的方法,并在多姿态人脸图像的重构实验中得到有效验证.通过这一逆向映射方法,与低维嵌入空间中的中心点对应的高维观测空间中的影像得以重构生成.具体步骤如下:

(1)降维.设对一幅大小为M×N的医学影像I(x,y)进行n个方向的和差定序滤波,得到的n个方向的影像定序特征Y(x,y,θi)(i=1,2,…,n),大小也为M×N.将这些影像定序特征作为n个样本,分别按列连接成M×N维的n个影像向量Yi∈RM×N(i=1,2,…,n).参数k为最近邻数,d为低维嵌入空间的维数.

图1 影像定序特征的提取及融合

定义代价函数为

固定权值W,由式(3)计算d维嵌入矢量{Xi}i=1,2,…,n∈Rd,即

(2)计算低维嵌入空间中Xi(i=1,2,…,n)的中心为

(3)融合.用其k邻域的加权和表示Xcenter,优化代价函数

得到权值矩阵

在LLE方法中,低维空间中的局部邻域关系在高维空间中保持不变.因此,与低维嵌入空间中的Xcenter对应的高维空间中的影像YF,可由式(6)得到的低维空间中的权值W重构,即

将YF按列还原为M×N的影像,该影像即为I(x,y)融合后的定序特征,如图1(i)所示.为了便于比较,文献[9]通过PCA方法得到的定序特征在图1(h)给出.可以明显看出,通过LLE方法得到的定序特征比PCA方法得到的定序特征包含更清晰、更丰富的空间信息.采用均值、标准差等图像融合的指标客观评价两种方法,比较结果如表1所示.除了熵的值几乎相等,在其他指标上均是LLE方法明显优于PCA方法的.因此, LLE方法比PCA方法更能揭示影像多方向定序特征的非线性分布,由文中方法得到的定序特征中,冗余信息被最小化,互补信息得以保留并加强,同时不易受到噪声的影响,鲁棒性强.

将待配准的两幅医学影像中分辨率较高的一幅影像设为参考影像R,其表示影像整体空间结构的定序特征设为X,另外一幅影像设为浮动影像F.这里采用一种基于混合熵的配准相似性测度[9],定义为

其中,H(F)为F的边缘熵;H(R,X)为R和X的二阶联合熵;H(R,F,X)表示R、F和X的三阶联合熵.该测度将R和X的联合(R,X)作为一个整体与F进行配准,配准过程视为使得(R,X)和F之间互信息最大的参数求解过程.在基于混合熵的配准测度中,影像定序特征的引入不仅增加了互信息中缺失的空间信息,使得配准测度能够同时衡量待配准影像对应灰度信息和空间信息的一致性,而且能够有效抑制噪声的影响,使配准测度函数更加平滑,有效避免了在配准搜索过程中陷入局部极值,提高了配准的精度.

表1 PCA方法和LLE方法的性能比较

2 配准方法

任何非刚性配准方法都可以描述为3个部分:联系浮动影像与参考影像的空间变换,测量浮动影像与参考影像相似程度的相似性测度,确定相似性测度中最优变换参数的优化方法.文中配准方法的相似性测度在上节中已给出相应的数学描述,以下对文中方法的空间变换与优化方法进行介绍.

2.1 空间变换模型

在医学影像非刚性配准算法中,为适合各种复杂的组织形变建立合理的空间变换模型是非常活跃的研究领域[1-2].一直以来,研究人员提出了各种方法.基于基函数的方法使用基函数的线性组合来描述变形域,然而它补偿解剖形变的能力非常有限,通常用于全局粗配准.基于薄板样条的方法是弹性形变中弹性薄片应力最小的一种插值模型,但是因为每一个控制点对变换都具有全局影响,所以很难模拟局部变形,不适合具有局部几何差异的图像配准.同薄板样条相比,B样条可以控制局部变形,改变控制点只影响它附近局部邻域的形状改变,其中基于B样条的自由形变模型(FFD)是目前使用较多的一种样条配准方法[14-15],其基本思想是:变形操作是作用于物体所嵌入的变形空间,而不是直接作用于物体;如果变形空间被改变了,则嵌入其中的物体自然随之改变.基于B样条的FFD模型是通过操纵分布在图像区域上的控制点网格来完成非线性变形的,以下是对该模型的数学描述.

对于二维情形,定义平面上的矩形域为

令Φ表示覆盖在矩形域Ω上由nx×ny个控制点φi,j所组成的网格,每个控制点之间的间距为δ,则该位移场定义为

其中,0≤u≤1.

B样条具有局部特性,即控制点位置的改变仅影响该控制点周围邻域网格内像素点的位移.因此,更新某个控制点的位置时,仅需更新该控制点邻域网格内的点;进行相似度计算时,也只需计算两幅影像对应控制点邻域的相似度值,因此算法的计算量大大减小.在文中配准方法中,首先对待配准影像进行刚性全局粗配准,再对发生局部形变的感兴趣区域(ROI)用基于B样条的FFD模型进行非线性空间变换.

2.2 改进的L-BFGS优化方法

在确定了空间变换类型和相似性测度函数之后,配准过程实质上就是一个多参数的优化求解过程,目的是寻找最佳空间变换参数,使得配准相似性测度达到最大值.所以优化方法的选择十分重要,快速有效的优化算法可以大大节省运行时间,为实时图像处理提供可能性.

L-BFGS优化方法是在BFGS算法的基础上发展起来的一种算法,非常适合求解高维度和大规模的无约束优化问题,它并不指定Hessian矩阵的形式和存储,直接计算搜索方向,大大降低了对内存的需求,提高了计算速度[16].在L-BFGS算法中,位置迭代公式为

其中,αk为步长因子,dk为搜索方向.

但在实际的非刚性医学影像配准过程中,L-BFGS优化方法存在以下缺陷:

(1)当待配准的两幅医学影像中存在较大非线性形变、局部形变不连续、灰度突变或者存在噪声等情况时,配准相似性测度会出现异常波动、起伏剧烈、局部极值增多、光滑性和收敛性恶化,从而导致优化算法在开始时,位置迭代过程中αk数值很大,最终很有可能会错过最优空间变换参数.

(2)在优化过程中,L-BFGS优化方法无法摆脱局部极值的吸引,因此有必要在保留当前最优位置的同时引入变异算子.但是,若将dk全部重新初始化,将会完全破坏当前搜索,使收敛速度大大减缓.

针对上述情况,将L-BFGS优化方法进行如下改进:

(1)在位置迭代更新前,对αk设定阈值进行限制,防止误配情况的发生.

(2)若当前位置xk变化极小时,按随机几率将dk中少部分维数重新初始化,以此增强全局搜索能力,使算法不易陷入局部极值,提高搜索精度.

最终,满足改进的L-BFGS优化方法迭代终止条件的空间变换参数输出为最优配准参数,完成配准.

3 实验结果与分析

为了对比分析文中方法的有效性,选取相同算法框架的互信息测度[4]、互信息与影像梯度结合测度[6]、二阶互信息测度[8]和文中式(8)定义的测度(以下简称为测度1～测度4)对医学影像进行配准实验.实验环境为Pentium(R)Dual-Core CUP 2.50 GHz和RAM 2 GB的PC机,操作系统为Windows XP,实验用Matlab7.6编程实现.

3.1 测试不同优化方法的搜索性能

将同一病人脑部的不同切片分别作为浮动影像和参考影像,大小均为256×256,如图2所示.对测度2分别采用L-BFGS方法和文中改进的L-BFGS方法进行最优配准参数搜索,对比结果如图3所示.可以看出,由于测度2光滑性较差,L-BFGS优化方法的位置迭代更新步长过大,并且最终陷于局部极值,导致配准结果失真;而改进的L-BFGS优化方法,由于能及时判断是否已收敛于局部最优点,并通过变异操作迅速摆脱它的束缚,提高了搜索成功率.实验表明,该改进算法是非常有效的.

图2 配准数据1

3.2 测试不同配准测度的精度和鲁棒性

为验证4种配准测度对噪声的鲁棒性,在图2两幅影像中均加入均值为0、方差为0.01的高斯噪声,采用改进的L-BFGS优化方法进行最优配准参数搜索,配准结果如图4所示.

因受噪声影响,直观上并不能明显看出4种测度配准结果的优劣.为客观评价4种测度的配准精度,采用灰度差平方和(SSD)和相关系数(CC)作为评价标准[14],SSD数值越大,表明配准后影像与参考影像差值越大,配准精度低;CC值越大,表明配准后影像与参考影像相似度越高,配准精度高.4种配准测度的精度比较结果如表2所示.可以看出,在相同噪声条件下,测度2配准效果最差,测度1和测度3相差不大,文中测度始终保持较小的误差,配准精度高于其他测度.

表2 N(0,0.01)高斯噪声条件下4种配准测度的精度对比

图3 两种优化方法的实验结果比较

图4 N(0,0.01)高斯噪声条件下4种配准测度的实验结果比较

3.3 测试不同配准测度的运行时间

实验采用一对CT/MR影像,如图5所示,大小均为256×256.以CT影像作为参考影像,MR影像作为浮动影像.4种配准测度单次运行时间比较结果见表3.可以看到,测度1用时最少,其他改进测度的运算时间均有所增加,其中测度1、测度2定义中只包含二维联合概率分布,节约了较多的内存,运行时间较短;测度3定义中包含四维联合概率分布,运行时间最长;文中测度定义中虽然只包含三维联合概率分布,但LLE算法的正向与反向映射还是消耗了一些运算时间,故文中测度算法效率仅优于测度3,仍需进一步提高.

图5 配准数据2

表3 4种配准测度的运行时间对比

4 结束语

提出了一种基于LLE和改进的L-BFGS优化的非刚性医学影像配准方法,运用基于B样条的FFD模型进行空间变换;针对高维定序特征,采用LLE及其逆映射进行降维和融合处理;得到的医学影像整体定序特征与影像灰度构成特征空间,基于混合熵的配准测度通过衡量两幅影像对应灰度信息和空间信息的一致性,使配准测度函数更加平滑,减少了配准过程中可能陷入的局部极值;在最优配准参数的优化搜索过程中,改进的L-BFGS优化方法采用变异操作增强了跳出局部最优解的能力.实验表明,文中方法能有效抑制噪声影响,配准精度较高.后续作者将致力于有效降低算法执行时间,使其能够满足临床的需要.

[1]Sotiras A,Daratzikos C,Paragios N.Deformable Medical Image Registration:a Survey[J].IEEE Transactions on Medical Imaging,2013,32(7):1153-1190.

[2]Liao S,Chung A C S.Nonrigid Brain MR Image Registration Using Uniform Spherical Region Descriptor[J].IEEE Transactions on Image Processing,2012,21(1):157-169.

[3]Wang Shanhu,You Hongjian,Fu Kun.BFSIFT:a Novel Method to Find Feature Matches for SAR Image Registration [J].IEEE Geoscience and Remote Sensing Letters,2012,9(4):649-653.

[4]Pluim J P W,Maintz J B A,Viergever M A.Mutual-information-based Registration of Medical Images:a Survey[J]. IEEE Transactions on Medical Imaging,2003,22(8):986-1004.

[5]So R W K,Tang T W H,Chung A C S.Non-rigid Image Registration of Brain Magnetic Resonance Images Using Graph-cuts[J].Pattern Recognition,2011,44(10-11):2450-2467.

[6]Pluim J P W,Maintz J B A,Viergever M A.Image Registration by Maximization of Combined Mutual Information and Gradient Information[J].IEEE Transactions on Medical Imaging,2000,19(8):809-814.

[7]Rueckert D,Clarkson M J,Hill D L G,et al.Non-rigid Registration Using Higher-order Mutual Information[C]// Proceedings of SPIE—the International Society for Optical Engineering.Bellingham:SPIE,2000:438-447.

[8]Chen B J,Li J L,Chen G.Study of Medical Image Registration Based on Second-Order Mutual Information.[C]// Proceedings of the 1st International Conference on Bioinformatics and Biomedical Engineering.Piscataway:IEEE,2007: 956-959.

[9]李琦,姬红兵,同鸣.一种多模态医学影像鲁棒配准方法[J].西安电子科技大学学报,2011,38(6):30-36.

Li Qi,Ji Hongbing,Tong Ming.Robust Registration of Multimodality Medical Images Based on the Principal Ordinal Feature and Hybrid Entropy[J].Journal of Xidian University,2011,38(6):30-36.

[10]Tan T N,Sun Z N.Ordinal Representations for Biometrics Recognition[C]//Proceedings of 15th European Signal Processing Conference.Poland:EURASIP,2007:35-39.

[11]Han Y F,Tan T N,Sun Z N,et al.Embedded Palmprint Recognition System on Mobile Devices.[C]//Lecture Notes in Computer Science:4642.Heidelberg:Springer Verlag,2007:1184-1193.

[12]Roweis S T,Saul L K.Nonlinear Dimensionality Reduction by Locally Linear Embedding[J].Science,2000,290 (5500):2323-2326.

[13]Zhang C S,Wang J,Zhao N Y,et al.Reconstruction and Analysis of Multi-pose Face Images Based on Nonlinear Dimensionality Reduction[J].Pattern Recognition,2004,37(2):325-336.

[14]Rueckert D,Sonoda L I,Hayes C,et al.Non-rigid Registration Using Free-form Deformations:Application to Breast MR Images[J].IEEE Transactions on Medical Imaging,1999,18(8):712-721.

[15]Díez Y,Oliver A,LladóX,et al.Revisiting Intensity-based Image Registration Applied to Mammography[J].IEEE Transactions on Information Technology in Biomedicine,2011,15(5):716-725.

[16]Morales J L,Nocedal J.Remark on“Algorithm 778:L-BFGS-B:Fortran Subroutines for Large-scale Bound Constrained Optimization"[J].ACM Transactions on Mathematical Software,2011,38(1):7.

(编辑:齐淑娟)

Non-rigid registration of medical images based on local linear embedding and improved L-BFGS optimization

LI Qi,JI Hongbing,ZANG Bo,LIU Jin
(School of Electronic Engineering,Xidian Univ.,Xi’an 710071,China)

Non-rigid registration of medical images has become a challenging task in medical image processing and applications.In this paper,we propose a local linear embedding(LLE)and improved LBFGS(limited-memory Broyden Fletcher Goldfarb Shanno)optimization based registration method.With abundant spatial information and good stability in noisy environment,the ordinal features are computed on different orientations to represent spatial information in medical images.For high dimensional ordinal features,the LLE algorithm is used for dimensionality reduction and the inverse mapping of LLE is used to fuse complementary information together.Then a hybrid entropy based similarity measure which integrates image intensity with ordinal feature is chosen as the registration function.Finally an improved L-BFGS algorithm is used to search for the optimal registration parameters.We evaluate the effectiveness of the proposed approach by applying it to the simulated brain image data.Experimental results show that the proposed registration algorithm is less sensitive to noise in images.Compared with some traditional methods,the proposed algorithm is of higher precision and better robustness.

non-rigid registration;local linear embedding;ordinal feature

TP391.4

A

1001-2400(2014)05-0054-07

2013-05-13< class="emphasis_bold">网络出版时间:

时间:2014-01-12

国家自然科学基金资助项目(61101246);中央高校基本科研业务费专项资金资助项目(JB140209,72125748)

李 琦(1979-),女,讲师,西安电子科技大学博士研究生,E-mail:qili@xidian.edu.cn.

http://www.cnki.net/kcms/doi/10.3969/j.issn.1001-2400.2014.05.010.html

10.3969/j.issn.1001-2400.2014.05.010