一种面向医学图像非刚性配准的多维特征度量方法

陆雪松　涂圣贤　张素

一种面向医学图像非刚性配准的多维特征度量方法

陆雪松1涂圣贤2张素2

医学图像的非刚性配准对于临床的精确诊疗具有重要意义.待配准图像对中目标的大形变和灰度分布呈各向异性给非刚性配准带来困难.本文针对这个问题，提出基于多维特征的联合Renyi α-entropy度量结合全局和局部特征的非刚性配准算法.首先，采用最小距离树构造联合Renyi α-entropy，建立多维特征度量新方法.然后，演绎出新度量准则相对于形变模型参数的梯度解析表达式，采用随机梯度下降法进行参数寻优.最终，将图像的Canny特征和梯度方向特征融入新度量中，实现全局和局部特征相结合的非刚性配准.通过在36对宫颈磁共振（Magnetic resonance，MR）图像上的实验，该方法的配准精度相比较于传统互信息法和互相关系数法有明显提高.这也表明，这种度量新方法能克服因图像局部灰度分布不一致造成的影响，一定程度地减少误匹配，为临床的精确诊疗提供科学依据.

非刚性配准，联合Renyi α-entropy，最小距离树，局部特征，自由形变模型

引用格式陆雪松，涂圣贤，张素.一种面向医学图像非刚性配准的多维特征度量方法.自动化学报，2016，42（9）：1413-1420

近些年来，随着医疗成像技术的进步，医学图像的质量和数量迅猛提升［1］.作为基础环节的医学图像配准技术，已成为临床诊疗的坚实依靠.通常，待配准图像对中目标的“运动”大多由患者体位变化、人体的消化和排泄器官的填充程度、呼吸等情况造成［2-3］.例如，盆腔部的宫颈会因治疗和周围膀胱、直肠等器官的扩张、收缩等有强烈变化（如图1所示）.因此，对于这些“运动”必须采用非刚性变换模型进行描述.

目前，医学图像的非刚性配准大都包括变换模型、度量准则、图像插值和优化方法四个部分［4］.其中相似性度量起着非常关键的作用［5］，它的选取都是以图像中所包含的特征为依据.传统互信息（Mutual information，MI）［6］和互相关系数（Correlation coefficient，CC）［7］直接利用图像灰度信息，不需要对图像做分割、特征提取等预处理，可以完全实现自动化，因而得到了广泛的应用和发展.然而，一些因素严重阻碍了这些配准度量的进一步发展.例如，待配准图像对中的像素呈现出很高的各向异性，组织器官的形状和大小有巨大变化等.

图1　膀胱和宫颈放射治疗前后磁共振（Magnetic resonance，MR）图像对比Fig.1Contrast of magnetic resonance（MR）images before and after radiotherapy at the bladder and the cervix

对于此种问题，一些研究者在相似性度量中融合了局部特征.Mellor等［8］提取图像局部结构的相位构造相位互信息代价函数，对医学图像进行非刚性配准.Loeckx等［9］在信息理论中从互信息演绎到条件互信息，用于图像的非刚性配准.实验表明，这些将图像划分为小子块的方法有时会受到采样点稀疏问题的影响.Neemuchwala等［10］由Renyi α-entropy提出一种多特征互信息α-MI，用于乳腺超声图像的非刚性配准，能将提取的多种特征融入其中，是传统互信息在多维特征空间上的一种推广. Staring等［11］将图像的局部灰度和梯度特征融入到α-MI中，对宫颈MR图像进行非刚性配准，实验结果显示其配准效果明显好于仅靠灰度的传统互信息.

多特征互信息采用熵图方法来度量多维特征［12］，尤其是局部特征，对于局部区域的像素各向异性和大形变是有效的.一般地，两种熵图最小距离树（Minimum spanning tree，MST）和K—最近邻图（K-nearest neighbors graph，KNNG）在图像配准中应用较为常见.多特征互信息度量值的大小往往与距离图的总长度有关，而且大多采用梯度下降法进行参数寻优.但是熵图的构造较为耗时，进行非刚性配准的速度较慢.另外，针对不同的配准对象，多特征互信息需要采用不同的局部特征才能获得较好的配准结果.

本文针对宫颈MR图像的非刚性配准难题，采用两幅图像的联合Renyi α-entropy度量其多维局部特征，提高组织器官和病灶区域的配准精度.主要优势体现在：1）采用最小距离树构造联合Renyi α-entropy，实现快速的自动配准；2）演绎出这种新的度量准则相对于形变模型参数的梯度解析表达式，以便利于采用梯度下降法解决配准问题；3）将图像的Canny特征和梯度方向特征融入度量准则中，对于因像素的各向异性、组织大形变等引起的配准难题进行了有益的探索.

本文第1节介绍非刚性配准的整体框架；第2节分析现有特征描述子的一些特点，在此基础之上给出本文所用多维特征的详细描述；第3节直接阐述基于最小距离树的联合Renyi α-entropy度量方法；第4节说明该度量方法相对于形变模型的梯度解析表达式，以便于配准的参数寻优；第5节展示实验结果；第6节进行总结.

1　非刚性配准的整体框架

医学图像配准的目的是发现一个最佳变形场T：（x，y，z）→（x＇，y＇，z＇），能将浮动图像I（x，y，z）中的点映射到参考图像I（x＇，y＇，z＇）中的对应点.这里，采用一个结合的变形场

其中，全局变形场Tglobal是仿射模型，局部变形场Tlocal是基于B样条的自由形变（Free-form deformation，FFD）模型［13］.仿射配准的结果被看成采用FFD模型配准的初始参数.

图2显示了整体的配准流程.对于全局变形场的粗对齐，仿射模型可表达为

其中，12个形变参数可通过梯度下降法对传统互信息度量进行迭代寻优获取.而对于局部变形场的精配准，基于B样条的FFD模型可表达为

图2　非刚性配准的框架图Fig.2The framework of nonrigid registration

其中，Bl表示B样条函数的第l阶基函数.

2　图像的多维特征

对于待配准图像对中存在目标大形变和灰度分布各向异性等问题，一些局部结构描述子被提出并被验证其有效性.例如，SIFT（Scale invariant feature transform）［14］、SURF（Speeded up robust features）［15］和DAISY［16］等.但是它们中的大多数都是稀疏、高维描述子，在工程实践中很难被应用到新配准度量中.与α-MI类似，良好的多维特征对于Renyi α-entropy相似性度量也是至关重要的.对于宫颈MR图像的配准，包含图像梯度特征的Cartesian结构特征集往往被使用［11］.但是由于它的维度较高、运算量较大，在工程实际中需要采用主成分分析（Principal component analysis，PCA）等方法进行降维处理.

为减轻非刚性配准负荷，本文采用一种混合的特征矢量计算Renyi α-entropy.首先，记配准原灰度图像为I（x），那么其高斯平滑图像可表示为Bσ（x），σ=1，2.对于图像的边界特征，采用Canny边缘［17］检测器Eσ（x），尺度也取为σ=1，2.再者，类似于SIFT算子的梯度方向特征对配准有帮助.假设I（x）在尺度σ下像素点x的梯度矢量被定义为（v1，v2，v3），那么两个方向角θσ（x）和φσ（x）可以表示为

其中的一些特征实例见图3所示.

3　基于最小距离树的联合Renyiα-entropy度量准则

假设随机变量x的概率密度函数为f（x），那么其α阶的Renyi entropy可定义为［18］

当α→1，Renyi α-entropy就趋近到香农熵.通常，多变量分布的熵大都可以采用一个最小图的长度进行估计.基于构造图的代价较低，一些研究者选择K—最近邻图来获取密度函数.然而，另外一些学者却采用最小距离树直接计算熵值［19］.他们争辩这种方法不需要估计概率密度，而且需要计算的边的数量比K—最近邻图少.因此，本文选用最小距离树构建联合Renyi α-entropy度量准则.

假设在一个d维特征空间有n个独立的多变量点，由它们构成一个图.其中顶点集合为V=x1，x2，···，xn，边集合为E=（xi，xj）；1≤i＜j≤n.那么此图的最小距离树长度可表示为

图3　特征图像实例Fig.3Examples of the image features

其中，eij=‖xi-xj‖表示点i和j之间的欧氏距离，T是包含所有顶点的全部图的集合，权重γ∈（0，d）.于是，概率密度函数f（x）与最小距离树长度之间的关系有

其中，α=（d-γ）/d，c是常量.因此，Renyi αentropy可表示为［20］

那么，它们的联合Renyi α-entropy度量就可表示为

4　非刚性配准的优化策略

通常，基于梯度的参数寻优方法比非梯度方法速度快.本文采用随机梯度下降法［21］进行非刚性配准，因此需要获得Renyi α-entropy度量相对于

FFD模型的梯度解析表达式.为更清晰地表达，假定

那么，联合Renyi α-entropy的梯度为

接着，Lfm（µ）的梯度可写为

5　配准实验结果

我们将这种新算法实现在Elastix［22］中，它是一种基于ITK（Insight toolkit）［23］框架的开源配准软件包.19个病人的宫颈MR三维数据被用于非刚性配准实验，图像大小为512像素×512像素× 30像素，空间分辨率为0.625mm×0.625mm× 4.5mm.在为期五周的放射治疗中，每个病人每周扫描一次.在配准前，图像大小被裁剪为210像素×250像素×30像素.为便于算法比较，配准执行在第2周和第3周（18对）、第3周和第4周（18对）的扫描数据上，共有36对三维图像数据参加实验.

待配准图像中的临床靶区（Clinical target volume，CTV）、膀胱（Bladder）和直肠（Rectum）三个区域的人工分割结果被用于配准质量的评估.通过配准结果获得的形变场，将浮动图像的人工分割结果进行形变，可作为参考图像的自动分割结果.通过对比参考图像的人工分割与自动分割结果，就能评估配准算法的性能.通常，作为重叠率的测量，参考图像的人工分割与自动分割结果之间的Dice相似系数［24］被计算.DSC=0表示两者之间没有重叠，DSC=1表示两者完美的对齐.当两种配准算法进行比较时，双边的Wilcoxon检验［25］被执行在对应的DSC值上.一个p＜0.05的值表示在统计上有明显的差异.

在采用B样条的FFD模型进行局部非刚性配准中，为避免局部极值的产生，一个三级多分辨率策略被使用.这里，高斯平滑而不是降采样被用到浮动图像中.对于x和y轴方向，平滑参数σ=4.0，2.0，1.0；对于z轴方向，平滑参数σ=2.0，1.0，0.5.至于B样条的控制点，三个分辨率水平下的网格控制点间距分别为80mm、40mm、20mm.在随机梯度下降的参数优化中，我们选择A=50，τ=0.602，a =2000.在所有的配准实验中，600次迭代被设置，每次迭代的随机采样点数被设置为N=5000.另外，新度量方法的α被设置为0.99.

在36对配准实验中，首先进行全局的仿射粗对齐.然后，在此基础上分别采用互相关系数、互信息度量、基于图的广义互信息（Generalized mutual information，GMI）［11］和新配准算法进行局部的精配准.同时，后两种方法都采用本文介绍的多维特征.图4显示了所有配准结果的DSC值比较.从图4可以看出，仿射粗配的结果最差.互信息与互相关系数的配准结果相近，在CTV区域互信息比互相关系数有明显提高，其他区域没有明显变化.GMI和本文方法的配准结果最好，但在直肠区域后者比前者有明显提高.相比较于互信息，本文方法在三个解剖区域都有不同程度的明显提高.特别是膀胱区域，重叠率中值从0.75增长到0.81（p=1.4×10-6）.所有配准DSC的均值和方差如表1.图5是一个典型的配准结果示例，其中图5（c）显示了依据互信息配准结果的融合图像，图5（d）显示了依据本文方法配准结果的融合图像.我们不难发现参考图像图5（a）和浮动图像图5（b）中的膀胱区域有很大区别，它们之间有着巨大的形变和灰度不一致.采用传统互信息在该区域明显失配，而本文方法的配准结果尽管不是非常完美，但是效果更好.

6　结论

本文将两幅图像的联合Renyi α-entropy引入多维特征的度量，通过最小距离树进行模型构建，选择图像的Canny和梯度方向等局部特征抑制误匹配，采用随机梯度下降法实现快速自动的非刚性配准.最后的实验结果表明对于大形变和灰度分布呈各向异性等复杂情况，本文方法比传统互信息法和互相关系数法表现出更好的性能.尽管本文的实验主要集中在MR单模态图像上，但是我们相信这种新方法也适合多模态图像，例如CT和MR图像间的配准.

类似宫颈MR等医学图像的非刚性配准是一项非常具有挑战性的任务，更快的配准速度和更高的配准精度在临床应用中仍然是需要的.根据不同的配准对象，发展更鲁棒和更尖端的局部结构描述子是我们进一步工作的发展方向.同时，更快的最小距离树初始图构造方法也将使新算法在工程实践中更具价值.

表1　所有配准DSC的均值和方差Table 1The mean and variance of DSC about all registration results

图4　仿射粗配、互相关系数、互信息、基于图的广义互信息和本文方法的重叠率对比boxplot图（对于每一个解剖结构，最左边的列显示了仿射粗配的结果，第2列显示了互相关系数的结果，第3列显示了互信息的结果，第4列显示了基于图的广义互信息的结果，最右边的列显示了本文方法的结果.一个星号表示两种方法重叠率中值在统计上的明显差异.）Fig.4The boxplot of overlap scores using affine matching，CC，MI，GMI and our proposal（For each anatomical structure，the leftmost column shows the result for affine matching.The second column shows the result for CC.The third column shows the result for MI.The fourth column shows the result for GMI.The rightmost column shows the result for our proposal.A star indicates a statistical significant difference of the median overlaps of the two methods.）

致谢

感谢荷兰Leiden大学Marius Staring博士提供数据支持.

图5　配准结果示例，参考图像与被形变的浮动图像采用Checkerboard模式融合Fig.5Demonstration of the registration result，and the reference image is combined with the deformed moving image，using a Checkerboard pattern

References

1 Savage N.Technology：multiple exposure.Nature，2013，502（7473）：S90-S91

2 Sotiras A，Davatzikos C，Paragios N.Deformable medical image registration：a survey.IEEE Transactions on Medical Imaging，2013，32（7）：1153-1190

3 Schmidt-Richberg A，Werner R，Handels H，Ehrhardt J. Estimation of slipping organ motion by registration with direction-dependent regularization.Medical Image Analysis，2012，16（1）：150-159

4 Zhang Gui-Mei，Cao Hong-Yang，Chu Jun，Zeng Jie-Xian. Non-rigid image registration based on low-rank Nystr¨om approximation and spectral feature.Acta Automatica Sinica，2015，41（2）：429-438（张桂梅，曹红洋，储王君，曾接贤.基于Nystr¨om低阶近似和谱特征的图像非刚性配准.自动化学报，2015，41（2）：429-438）

5 Zhou Bin，Song Xiao，Huang Xiao-Yang，Wang Bo-Liang. A rigid registration method based on joint distribution histogram and its application on liver CT images.Journal of Xiamen University（Natural Science），2015，54（3）：397-403（周斌，宋晓，黄晓阳，王博亮.联合直方图的刚性配准方法及其在肝脏CT图像配准中的应用.厦门大学学报（自然科学版），2015，54（3）：397-403）

6 Xia Wei，Gao Xin，Wang Lei，Zhou Zhi-Yong，Zhang Ran. 3D deformable registration algorithm of CT lung images based on point set and mutual information.Journal of Jiangsu University（Natural Science Edition），2014，35（5）：558-563（夏威，高欣，王雷，周志勇，张冉.基于点集与互信息的肺部CT图像三维弹性配准算法.江苏大学学报（自然科学版），2014，35（5）：558-563）

7 Cachier P，Bardinet E，Dormont D，Pennec X，Ayache N. Iconic feature based nonrigid registration：the PASHA algorithm.Computer Vision and Image Understanding，2003，89（2-3）：272-298

8 Mellor M，Brady M.Phase mutual information as a similarity measure for registration.Medical Image Analysis，2005，9（4）：330-343

9 Loeckx D，Slagmolen P，Maes F，Vandermeulen D，Suetens P.Nonrigid image registration using conditional mutual information.IEEE Transactions on Medical Imaging，2010，29（1）：19-29

10 Neemuchwala H，Hero A，Carson P.Image matching using alpha-entropy measures and entropic graphs.Signal Processing，2005，85（2）：277-296

11 Staring M，van der Heide U A，Klein S，Viergever M A，Pluim J P W.Registration of cervical MRI using multifeature mutual information.IEEE Transactions on Medical Imaging，2009，28（9）：1412-1421

12 Rivaz H，Karimaghaloo Z，Collins D L.Self-similarity weighted mutual information：a new nonrigid image registration metric.Medical Image Analysis，2014，18（2）：343-358

13 Feng Zhao-Mei，Dang Jun，Cui Xiao-Yao，Jiao Yang.B-spline free-form deformation based 3-D non-rigid registration of medical images.Imaging Science and Photochemistry，2014，32（2）：200-208（冯兆美，党军，崔崤山尧，焦阳.基于B样条自由形变三维医学图像非刚性配准研究.影像科学与光化学，2014，32（2）：200-208）

14 Lowe D G.Distinctive image features from scale-invariant keypoints.International Journal of Computer Vision，2004，60（2）：91-110

15 Bay H，Tuytelaars T，van Gool L.SURF：speeded up robust features.In：Proceedings of the 9th European Conference on Computer Vision.Graz，Austria：Springer，2006.404-417

16 Tola E，Lepetit V，Fua P.DAISY：an efficient dense descriptor applied to wide-baseline stereo.IEEE Transactions on Pattern Analysis and Machine Intelligence，2010，32（5）：815-830

17 Canny J.A computational approach to edge detection. IEEE Transactions on Pattern Analysis and Machine Intelligence，1986，PAMI-8（6）：679-698

18 Neemuchwala H F，Hero A O.Image registration in high dimensional feature space.In：Proceedings of the 2005 SPIE 5674，Computational Imaging III.San Jose，CA：SPIE，2005.99-113

19 Sabuncu M R，Ramadge P J.Using spanning graphs for efficient image registration.IEEE Transactions on Image Processing，2008，17（5）：788-797

20 Hero A O，Ma B，Michel O J J，Gorman J.Applications of entropic spanning graphs.IEEE Signal Processing Magazine，2002，19（5）：85-95

21 Klein S，Staring M，Pluim J P W.Evaluation of optimization methods for nonrigid medical image registration using mutual information and B-splines.IEEE Transactions on Image Processing，2007，16（12）：2879-2890

22 Klein S，Staring M，Murphy K，Viergever M A，Pluim J P W.Elastix：a toolbox for intensity-based medical image registration.IEEE Transactions on Medical Imaging，2010，29（1）：196-205

23 Ibanez L，Schroeder W，Ng L，Cates J.The ITK software guide version 2.4［Online］，available：http：//www.itk.org，September 16，2014

24 Dice L R.Measures of the amount of ecologic association between species.Ecology，1945，26（3）：297-302

25 Wilcoxon F.Individual comparisons by ranking methods. Biometrics Bulletin，1945，1（6）：80-83

陆雪松中南民族大学生物医学工程学院副教授.主要研究方向为医学影像配准、分割和辅助诊断.本文通信作者.

E-mail：xslu-scuec@hotmail.com

（LU Xue-SongAssociate professor at the College of Biomedical Engineering，South-Central University for Nationalities.His research interest covers medical image registration，segmentation，and assisted diagnosis.Corresponding author of this paper.）

涂圣贤上海交通大学生物医学工程学院东方学者特聘教授.主要研究方向为心血管成像与定量分析.

E-mail：sxtu@sjtu.edu.cn

（TU Sheng-XianProfessor of specialappointmentattheSchoolof Biomedical Engineering，Shanghai Jiao Tong University.His research interest covers cardiovascular imaging and quantitative analysis.）

张素上海交通大学生物医学工程学院副教授.主要研究方向为医学影像处理与分析.

E-mail：suzhang@sjtu.edu.cn

（ZHANG SuAssociate professor at the School of Biomedical Engineering，Shanghai Jiao Tong University.Her research interest covers medical image processing and analysis.）

A Metric Method Using Multidimensional Features for Nonrigid Registration of Medical Images

LU Xue-Song1TU Sheng-Xian2ZHANG Su2

Nonrigid registration of medical images has great significance for accurate diagnosis and therapy in clinic.It is difficult to register the images containing large deformation of object region and data anisotropy.According to this problem，an algorithm of nonrigid registration based on joint Renyi α-entropy is proposed in this paper，which combines global features with local features.Firstly，minimum spanning tree is employed for construction of joint Renyi α-entropy. A new metric is built on multidimensional features.And then，the analytical derivative of the new metric with respect to the parameters of deformation model is derived，in order to find the optima by a stochastic gradient descent method. Finally，Canny feature and gradient orientation feature of images are merged into the new metric，which implements nonrigid registration including global and local features.Experiments are performed on 36 cervical magnetic resonance（MR）image pairs.Compared to the traditional mutual information and correlation coefficient，the registration accuracy is improved significantly.It also manifests that the proposed method is able to overcome the adverse effects of local intensity inhomogeneity，and provides scientific evidence for accurate diagnosis and therapy in clinic，due to reducing mismatch in some degree.

Nonrigid registration，joint Renyi α-entropy，minimum spanning tree，local feature，free-form deformation model

Manuscript October 8，2015；accepted March 10，2016

10.16383/j.aas.2016.c150608

Lu Xue-Song，Tu Sheng-Xian，Zhang Su.A metric method using multidimensional features for nonrigid registration of medical images.Acta Automatica Sinica，2016，42（9）：1413-1420

2015-10-08录用日期2016-03-10

国家自然科学基金（61002046），国家民委科研项目（14ZNZ024）资助

Supported by National Natural Science Foundation of China（61002046）and Scientific Research Projects by the State Ethnic Affairs Commission of China（14ZNZ024）

本文责任编委黄庆明

Recommended by Associate Editor HUANG Qing-Ming

1.中南民族大学生物医学工程学院武汉4300742.上海交通大学生物医学工程学院上海200240

1.College of Biomedical Engineering，South-Central University for Nationalities，Wuhan 4300742.School of Biomedical Engineering，Shanghai Jiao Tong University，Shanghai 200240