一種面向醫學圖像非剛性配準的多維特征度量方法

陸雪松　涂圣賢　張素

一種面向醫學圖像非剛性配準的多維特征度量方法

陸雪松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