張 冉,王 雷,夏 威,高 欣
(1.中國科學院長春光學精密機械與物理研究所,吉林長春130033;2.中國科學院蘇州生物醫(yī)學工程技術(shù)研究所醫(yī)學影像室,江蘇蘇州215163;3.中國科學院大學,北京100049)
2D/3D圖像配準中的相似性測度和優(yōu)化算法
張 冉1,2,3,王 雷1,2,3,夏 威1,2,3,高 欣2
(1.中國科學院長春光學精密機械與物理研究所,吉林長春130033;2.中國科學院蘇州生物醫(yī)學工程技術(shù)研究所醫(yī)學影像室,江蘇蘇州215163;3.中國科學院大學,北京100049)
在手術(shù)引導治療中,2D/3D圖像配準能輔助醫(yī)生準確定位病人病灶,而準確的配準涉及相似性測度和優(yōu)化算法等眾多方面。為了研究相似性測度和優(yōu)化算法對2D/3D圖像剛性配準的影響,本文結(jié)合6種相似性測度和4種優(yōu)化方法在配準“金標準”數(shù)據(jù)上進行了2D/3D圖像配準實驗,并從配準成功率、平均迭代次數(shù)和平均配準時間三個方面對配準結(jié)果進行了對比研究。實驗結(jié)果表明,以模式強度為相似性測度,用Powell方法進行優(yōu)化搜索是最佳配準組合。并且,在不改變相似性測度條件下,Powell方法是所用優(yōu)化方法中配準效果最好的優(yōu)化方法。
2D/3D圖像配準;剛性配準;相似性測度;優(yōu)化方法;配準評估
在手術(shù)引導治療中,2D/3D圖像配準在輔助醫(yī)生準確定位病人病灶上是一種重要技術(shù)[1-3],被廣泛的應(yīng)用在圖像引導的微創(chuàng)手術(shù)、放射治療的計劃制定、術(shù)后治療的效果檢驗等介入手術(shù)方面。2D/3D圖像配準的引入可降低上述醫(yī)療技術(shù)的實現(xiàn)難度,減少介入性治療的創(chuàng)傷并提高手術(shù)精度。為了實現(xiàn)準確配準,科研工作者提出了大量的配準算法,基本可以分為基于特征[4-5]、基于灰度[4-7]和基于梯度[3]三類方法。其中,基于灰度的配準算法直接利用圖像的灰度值進行配準,無需進行分割處理和人工介入,可實現(xiàn)全自動配準且精度較高,因而成為2D/3D醫(yī)學圖像剛性配準的主要研究方向。該算法基于數(shù)字重建影像(Digitally Reconstructed Radiographs,DRR)技術(shù),首先采用Ray casting算法[4]將變換后的3D圖像投影為2D圖像即DRR圖像,選取適當?shù)南嗨菩詼y度作為目標函數(shù)來衡量DRR圖像與X-ray圖像之間的相似程度,然后對目標函數(shù)進行優(yōu)化,當函數(shù)達到最優(yōu)時即代表獲得最佳配準。
近年來,在基于灰度的配準算法中,科研人員對相似性測度和優(yōu)化方法做了一定的對比研究[4-5,8],以評定相似性測度和優(yōu)化算法對配準精度和魯棒性的影響。Penney等人[4]對歸一化相關(guān)(Normalized Cross Correlation,NCC)、互信息(Mutual Information,MI)、差值圖像的熵(Entropy of the Difference Image,EDI)、梯度相關(guān)(Gradient Correlation,GC)、梯度差分(Gradient Difference,GD)和模式強度(Pattern Intensity,PI)6種相似性測度進行了對比分析;Maes等人[5]研究了Powell方法、單純形法(Simplex,SMP)、梯度下降法(Steepest Gradient Descent,STD)、共軛梯度法(Conjugate Gradient,CJG)、擬牛頓法(Quasi Newton,QSN)和列文伯格-馬夸爾特法(Levenberg Marquardt,LVM)6種優(yōu)化算法對基于互信息的2D/3D圖像配準的影響,然而在2D/3D圖像配準中,僅僅考慮相似性測度或者優(yōu)化方法對配準結(jié)果的影響是不夠的,因此,本文結(jié)合6種相似性測度和4種優(yōu)化方法在公開發(fā)表的“金標準”數(shù)據(jù)[2]上進行2D/3D圖像剛性配準實驗,并利用Kraats等人[6]提出的方法從配準成功率、平均迭代次數(shù)和平均配準時間三個方面對配準結(jié)果進行評估,分析不同的相似性測度和優(yōu)化方法對圖像配準的影響,并給出通用性最好的組合,為2D/3D圖像剛性配準在臨床中的應(yīng)用提供理論基礎(chǔ)。
在2D/3D圖像剛性配準中,相似性測度是用來度量兩幅圖像配準程度的一個重要指標,可以衡量經(jīng)過空間變換后浮動圖像(DRR)與參考圖像(X-ray)在空間上的一致性程度。因此,相似性測度的選擇是否合理,直接影響到圖像配準的準確性和穩(wěn)定性。本文選取比較常用的6種相似性測度進行對比研究,它們是歸一化相關(guān)、互信息、歸一化互信息、梯度相關(guān)、梯度差分以及模式強度[4,7-10]。為了便于統(tǒng)一標識,定義IX( i,j)和IDRR( i,j)分別為X-ray圖像和DRR圖像中像素(i,j)處的灰度值。
2.1 歸一化互相關(guān)(NCC)
在2D/3D圖像剛性配準中,歸一化互相關(guān)測度是通過求兩幅圖像的相關(guān)系數(shù)來表征二者間的相關(guān)程度,用公式可表示為:
其中,IX,IDRR分別表示X-ray圖像和DRR圖像像素灰度的平均值。兩幅圖像越相似,該測度的值越接近1。
2.2 互信息(MI)
互信息是信息論里一種有用的信息度量,是指兩個事件集合之間的相關(guān)性,一般用熵來表示。圖像X的熵和圖像X,Y的聯(lián)合熵分別定義為:
其中,聯(lián)合概率分布p( x,y)用歸一化的聯(lián)合直方圖表示,p( x)、p( y)是p( x,y)的邊緣分布概率。
因此,兩幅圖像之間的互信息可定義為:
其中,()H X,()H Y分別是圖像X,Y的熵;H X,()Y是它們的聯(lián)合熵。在2D/3D圖像剛性配準中,互信息的值越大,兩幅圖像的相似性程度越高;當互信息達到最大值時,表示兩幅圖像達到最佳配準。
2.3 歸一化互信息(NMI)
由于互信息測度對配準圖像間的重疊程度較為敏感,研究人員提出了歸一化互信息的概念,其中Sthdholme等人[10]給出了其中一種形式,可表示為:
2.4 梯度相關(guān)(GC)
基于梯度的相似性測度首先要計算出兩幅圖像的梯度圖像(本文采用Sobel算子),然后根據(jù)公式(1)分別計算出梯度圖像在水平和豎直兩個方向的歸一化互相關(guān),它們的和就是兩幅圖像的梯度相關(guān),可表示為:
2.5 梯度差分(GD)
這種測度也是基于梯度圖像,在得到各向梯度圖像的基礎(chǔ)上,比較兩幅梯度圖像間的差異,實際上是比較兩幅圖像的邊緣方向信息。其計算公式如下:
其中,Aν和Ah是X-ray圖像分別在豎直和水平方向上的方差;s是衰減因子。IdiffV和IdiffH分別表示待配準圖像在豎直方向和水平方向上的梯度差值圖像。梯度差分測度通過使用1/(1+x2)的形式增強了對細線結(jié)構(gòu)的魯棒性[9],通過調(diào)整s,Aν和Ah的值可以使該測度更加平滑。
2.6 模式強度(PI)
該測度通過測量差值圖像(兩幅圖像的灰度差)中所存在的模式是否降到最低來判斷配準是否成功,其數(shù)學表達式為:
其中,σ,r是常量;Idiff( v,w)表示差值圖像中(v,w)處的灰度值。模式強度認為,當一個像素與其臨近的像素的值的差異顯著時,該像素便屬于一個模式,圖像配準的過程就是要盡量消除這種差異[7,9]。圖像達到最佳配準時,Idiff中待配準的模式會消失,其模式強度會最小,該測度的值最大。
除了相似性測度以外,優(yōu)化算法也會影響2D/ 3D圖像配準的精度。醫(yī)學圖像剛性配準在本質(zhì)上是一個多參數(shù)的最優(yōu)化問題,即尋找相似性測度達到最大時所對應(yīng)的6個空間變換參數(shù)的值,可由下式表示:
式中,C是由相似性測度構(gòu)成的目標函數(shù);IX和I3D分別是2D X-ray圖像和3D體數(shù)據(jù);μ是作用于體數(shù)據(jù)上的變換參數(shù),μ′是圖像配準時算法搜索到的最優(yōu)參數(shù)。
最優(yōu)化是一個迭代尋優(yōu)的過程,可以采用如下迭代格式[11]表示:
其中,dk表示第k次迭代時的搜索方向;αk是迭代步長。不同的優(yōu)化算法對應(yīng)著不同的搜索方向和搜索步長。文中選取4種優(yōu)化方法用于搜索圖像配準的最優(yōu)空間變換參數(shù),優(yōu)化方法具體如下[5,8,11-13]。
3.1 梯度下降法(GRadient Descent,GRD)
梯度下降法又稱最速下降法,其搜索方向沿著負梯度方向,目標函數(shù)可以最快的達到極小值,迭代格式可表示為:
3.2 非線性共軛梯度(Nonlinear Conjugate Gradient,NCG)
對于非線性共軛梯度,公式(13)中搜索方向dk可以表示為:
其中,βk公式有多種,Y H DAI,Y YUAN[12]給出βDYk的公式,Hestenes和Stiefel[8,12]給出βHSk的公式。本文采用的βk公式是二者的混合,即βk=分別由下面兩個公式求得:
3.3 Powell算法
Powell算法[11]是直接搜索法中一種比較有效的方法,它不需要目標函數(shù)的導數(shù)信息,只要求目標函數(shù)連續(xù)即可。由于原始Powell算法中無法保證搜索方向是線性無關(guān)的,因此本文中采用的是改進后的Powell算法[11]。
3.4 擬牛頓法(Quasi Newton,QN)
擬牛頓法是牛頓法的直接推廣,基本思想是構(gòu)造海森矩陣的近似矩陣Bk或逆海森矩陣的近似矩陣Hk,通過在試探點附近的二次逼近引導擬牛頓方程來確定搜索方向。擬牛頓方程主要有兩種形式:Hk+1yk=sk和Bk+1sk=y(tǒng)k,本文采用后者,其公式為:
其中,sk=xk+1-xk,yk=gk+1-gk,gk和gk+1分別表示x在第k和第(k+1)次的導數(shù)。
本文實驗采用MATLAB軟件平臺編程實現(xiàn)配準過程,同時結(jié)合CUDA技術(shù),電腦配置為Intel Xeon E5645,2.4 GHz處理器,內(nèi)存為48 GB,NVIDIA Tesla C2050顯卡,顯存為4 GB。
4.1 實驗數(shù)據(jù)
在配準實驗中,我們采用公開發(fā)表的2D/3D配準“金標準”數(shù)據(jù)[2],該數(shù)據(jù)是對豬的顱骨進行成像得到的,包含3D的CT體數(shù)據(jù)和2D的X-ray圖像。高斯濾波及各向同性采樣到1mm后,CT數(shù)據(jù)的大小為326×326×330,X-ray圖像尺度為410×410,如圖1(a)所示。由于配準數(shù)據(jù)同時包含大量的軟組織和硬組織,為提高配準精度,避免軟組織對剛性配準的影響,同時縮短配準時間,實驗選取包含軟組織較少的圓形區(qū)域作為目標區(qū)域,其直徑為200個像素。圖1(b)是去除標記點的X-ray圖像的興趣區(qū)。
圖1 豬顱骨X-ray圖像
4.2 配準評估
根據(jù)Kraats等人提出的2D/3D配準評估準則[6],本文從配準成功率(Success Rate,SR)、平均迭代次數(shù)以及平均配準時間三個方面對配準結(jié)果進行評估。配準成功率可以用平均目標配準誤差(mean target registration error,mTRE)來計算,mTRE公式如下:
其中,T和Tgold分別表示配準后計算得到的變換矩陣T(由配準算法決定)和金標準變換矩陣Tgold,pn為在CT體數(shù)據(jù)中隨機選取的目標點集。對“金標準”變換參數(shù)進行隨機擾動后得到初始變換參數(shù),使得初始mTRE分布在0~10 mm之間,并在每1 mm的間隔內(nèi)再隨機選取10個初始位置,這樣就產(chǎn)生了總共100次初始變換。
4.3 參數(shù)設(shè)置
本文中所有參數(shù)均是參照以往的文獻設(shè)置,對于梯度差分(GD),公式(7)中參數(shù)Aν和Ah均設(shè)為10,公式(8)和公式(9)中參數(shù)s均設(shè)為0.2。模式強度(PI)中參數(shù)σ設(shè)為10,r設(shè)為3,s設(shè)為0.2。
4.4 實驗結(jié)果
文中對6種相似性測度和4種優(yōu)化算法的不同組合進行配準對比,結(jié)果如表1和表2所示。由于“金標準”數(shù)據(jù)含有大量的軟組織并且在使用Ray Casting時沒進行優(yōu)化處理來生成DRR圖像,因此本文選取的配準成功閾值為8 mm,即最終目標配準誤差(mTRE)小于8 mm視為配準成功。
表1 所有相似性測度和優(yōu)化算法所進行配準組合的成功配準率
由表1可以看出,配準成功率范圍為0.61~0.81,最高成功率(0.81)的配準組合是PI和Powell組合,而最低成功率(0.61)的配準組合是MI和NCG組合。從結(jié)果中還可以看出,相同的相似性測度,Powell是成功率最高的優(yōu)化方法,而NCG的成功率最低。
除了配準成功率,本文還從平均配準時間和平均迭代次數(shù)上對2D/3D剛性配準進行了評估,其中,平均配準時間是指一次完整配準所需的平均時間,平均迭代次數(shù)是指一次完整配準所需的平均迭代次數(shù)。
表2 所有配準組合的平均配準時間(t)和平均迭代次數(shù)(N)
從表2可以看出,后三種相似性測度(GC、GD以及PI)無論用什么優(yōu)化策略進行配準所需時間都大于前面三種(NCC、MI以及NMI)。在所有的優(yōu)化方法中,Powell所需的平均配準時間和平均迭代次數(shù)都是最少的。因此,在2D/3D圖像剛性配準中,Powell成了最好的優(yōu)化方法。
本文對6種相似性測度和4種優(yōu)化方法的不同組合進行了2D/3D剛性配準對比研究,并從配準成功率、平均配準時間和平均迭代次數(shù)上進行了配準評估。結(jié)果表明,除了相似性測度會對2D/3D剛性配準產(chǎn)生影響外,優(yōu)化方法的選擇同樣會影響配準效果。就本文而言,以PI(模式強度)作為相似性測度,用Powell進行優(yōu)化搜索是最好的配準組合。此外,結(jié)果還表明Powell方法是2D/3D圖像剛性配準最好的優(yōu)化方法。
[1] PMarkelj,D Toma?evi?,B Likar,etal.A review of3D-2D registration methods for image-guided interventions[J]. Medical Image Analysis,2012,16(3):642-661.
[2] PMarkelj,F(xiàn) Pernu?,SA Pawiro,et al.Validation for 2D/3D registration I:a new gold standard data set[J].Medical physics,2011,38(3):1481-1490.
[3] Barbara R?per,Nassir Navab,Wolfgang Wein.2D/3D registration based on volume gradients[C]//Medical Imaging International Society for Optics and Photonics,2005:144-150.
[4] Jürgen Weese,Graeme P Penney,John A Little,et al.A comparison of similaritymeasures for use in 2D3D medical image registration[J].IEEE Transactions on Medical Imaging,1998,17(4):586-595.
[5] FMaes,Dirk Vandermeulen,Paul Suetens.Comparative evaluation ofmultiresolution optimization strategies formultimodality image registration bymaximization ofmutual information[J].Medical Image Analysis,1999,3(4):373-386.
[6] Graeme P Penney,Everine B van de Kraats,Dejan Toma?evi?.Standardized evaluationmethodology for 2D-3D registration[J].IEEE Transactions on Medical Imaging,2005,24(9):1177-1189.
[7] Liang Wei.2D-3D registration ofmedical image[D]. Nanjing:Southeast University,2004.(in Chinese)梁瑋.2D-3D醫(yī)學圖像配準研究[D].南京:東南大學,2004.
[8] IM Jvan der Bom,SKlein,M Staring.Evaluation of optimizationmethods for intensity-based 2D-3D registration in x-ray guided interventions[C]//SPIEMedical Imaging International Society for Optics and Photonics,2001,7962(23):1-15.
[9] Wang Kai.The study ofmultimodality medical image regitration[D].Changsha:Central South University,2008.(in Chinese)王凱.多模態(tài)醫(yī)學圖像配準研究[D].長沙:中南大學,2008.
[10]D L G Hill,C Studholme,D JHawkes.An overlap invariant entropymeasure of 3D medical image alignment[J]. Pattern Recognition,1999.
[11]Sun Jingming,Ling Yingchun.Optimal design ofmachine[M].Beijing:China Machine Press,2006.(in Chinese)孫靖明,梁迎春.機械優(yōu)化設(shè)計[M].北京:機械工業(yè)出版社,2006.
[12]Y H DAI,Y YUAN.An efficienthybrid conjugate gradient method for unconstrained optimization[J].Annals of Operations Research,2001,103(1):33-47.
[13]Ma Changfeng.Optimizationmethods and the design of its matlab program[M].Beijing:Science Press,2010.(in Chinese)馬昌鳳.最優(yōu)化方法及其Matlab程序設(shè)計[M].北京:科學出版社,2010.
Comparison of sim ilarity measurement and optim ization methods in 2D/3D image registration
ZHANG Ran1,2,3,WANG Lei1,2,3,XIAWei1,2,3,GAO Xin2
(1.Changchun Institute of Optics,F(xiàn)ine Mechanics and Physics,Chinese Academy of Sciences,Changchun 130033,China;2.Medical Imaging Department,Suzhou Institute of Biomedical Engineering and Technology,Chinese Academy of Sciences,Suzhou 215163,China;3.University of the Chinese Academy of Sciences,Beijing 100049,China)
In surgical guide treatment,2D/3Dmedical image registration can provide the precise position of patient for surgeon.Accurate registration involvesmany aspects,such as similaritymeasurements and optimizationmethods.In order to investigate the influence of similarity measurements and optimization methods on 2D/3D image registration,a comparison of six similaritymeasurements in combination with four optimizationmethods is performed using the public and available porcine skull phantom datasets from Medical University Vienna.Comparison is performed for the registration results based on success rate,the number of iterations and execution time.The results show that themost accuracy registration is obtained by pattern intensity combined with Powell.Furthermore,the best 2D/3D registration results are obtained by Powell search strategy with fixed similaritymeasurement.
2D/3D image registration;rigid registration;similaritymeasures;optimization;registration evaluation
TP391
A
10.3969/j.issn.1001-5078.2014.01.022
1001-5078(2014)01-0098-05
國家自然科學基金項目(No.81000651);蘇州市科技計劃項目(No.SH201210);江蘇省科技計劃項目(No.BL2012049)資助。
張 冉(1988-),女,碩士研究生,研究方向為醫(yī)學圖像處理。
2013-05-20