基于混合蛙跳算法的SPECT－B超甲狀腺圖像配準(zhǔn)

鄭偉，孟繁婧，田華，郝冬梅，吳頌紅

（1.河北大學(xué) 電子信息工程學(xué)院，河北 保定 071002；2.河北大學(xué)附屬醫(yī)院 功能檢查科，河北 保定 071002）

世界權(quán)威組織統(tǒng)計(jì)顯示，甲狀腺癌在2010年已躍居女性惡性腫瘤的第6位.目前，針對(duì)甲狀腺腫瘤的影像檢查分為功能成像（SPECT，PET 等）和解剖結(jié)構(gòu)成像（B超，CT，MRI等）2類，單一檢測(cè)方式的診斷準(zhǔn)確率較低.圖像融合技術(shù)能夠?qū)碜圆煌Ｊ降膱D像中的信息結(jié)合起來，為外科手術(shù)的規(guī)劃和放射治療計(jì)劃的設(shè)計(jì)提供依據(jù).美國(guó)通用電氣公司生產(chǎn)的SPECT／CT 融合設(shè)備，將SPECT 的功能圖像CT 圖像的準(zhǔn)確定位相結(jié)合，確定腫瘤的位置和大小.但CT 放射性會(huì)削弱人體抗菌能力，而且對(duì)于小于2cm 的甲狀腺腫瘤，B超檢查檢出率可高達(dá)67%，檢出率明顯高于CT，而且B超無放射性損傷.根據(jù)臨床實(shí)際需求，提出了將甲狀腺B超圖像和SPECT 圖像進(jìn)行融合處理，克服單模醫(yī)學(xué)圖像信息單一、表征局限的缺點(diǎn)，提供涵蓋功能信息和結(jié)構(gòu)信息的圖像，提高對(duì)甲狀腺腫瘤良、惡性判斷的準(zhǔn)確率.

圖像配準(zhǔn)是融合的基礎(chǔ)，醫(yī)學(xué)圖像配準(zhǔn)就是指將來自不同模式的醫(yī)學(xué)圖像，通過空間變換使一幅醫(yī)學(xué)圖像與另一幅醫(yī)學(xué)圖像上的對(duì)應(yīng)點(diǎn)或?qū)?yīng)輪廓達(dá)到空間上的一致，或至少是所有具有診斷意義的點(diǎn)及手術(shù)感興趣的點(diǎn)都達(dá)到匹配.因研究對(duì)象，研究目的不同，配準(zhǔn)的方法也多種多樣，常用的圖像配準(zhǔn)方法有基于灰度的方法、變換域法和基于特征的方法.現(xiàn)有的SPECT／CT 融合設(shè)備為同機(jī)配準(zhǔn)，采用的是衰減校正質(zhì)量控制配準(zhǔn)方法［1］.文獻(xiàn)［2］提出基于人工免疫系統(tǒng)模型算法的點(diǎn)對(duì)點(diǎn)的配準(zhǔn)算法，證明其性能的優(yōu)越性和高精度性；文獻(xiàn)［3］對(duì)人腦的CT 和MR 采用一種新的非均勻采樣的方法提高了互信息的精確度；文獻(xiàn)［4］則對(duì)人腦的MR 和CT 圖像提出了一種改進(jìn)式同步擾動(dòng)隨機(jī)算法對(duì)配準(zhǔn)的過程進(jìn)行了優(yōu)化；文獻(xiàn)［5］采用大變形微分同胚尺度映射算法進(jìn)行配準(zhǔn)，高斯－牛頓策略則實(shí)現(xiàn)了更快的收斂速度；文獻(xiàn)［6］提出了基于穩(wěn)健點(diǎn)集的隨機(jī)全局優(yōu)化配準(zhǔn)；文獻(xiàn)［7］提出了一種改進(jìn)的人工魚群算法和powell算法結(jié)合完成了在低分辨率下的圖像配準(zhǔn)；文獻(xiàn)［8］首先把點(diǎn)集匹配問題轉(zhuǎn)化為解空間為仿射參數(shù)空間下的目標(biāo)函數(shù)優(yōu)化問題，然后運(yùn)用粒子群算法對(duì)相應(yīng)的變換參數(shù)進(jìn)行搜索，獲得問題最優(yōu)解；文獻(xiàn)［9］提出基于邊緣保護(hù)多尺度空間配準(zhǔn)以及自動(dòng)獲取非線性擴(kuò)散模型中平滑參數(shù)入的方法來提高配準(zhǔn)的精度和速度，避免局部極值.B超圖像和SPECT 圖像成像原理不同，同一器官表現(xiàn)出來的圖像特征差異較大，基于灰度的方法和變換域法不適用，而基于特征的方法適合用于有明顯特征的圖像，B超圖像和SPECT 圖像輪廓特征明顯，適合用該方法.

1 SPECT 和B超圖像特性分析

從原始的B超圖像（如圖1a所示）中可以看到甲狀腺的兩葉、腫瘤、頰部、氣管，其周圍的神經(jīng)、血管和甲狀軟骨等.SPECT 對(duì)于一些涼結(jié)節(jié)和冷結(jié)節(jié)不顯像，亦無法顯示氣管或病灶的解剖細(xì)節(jié)，圖像中只能觀察到甲狀腺區(qū)和其上方的放射線減淡區(qū)，即為涼結(jié)節(jié)（如圖1b所示）.醫(yī)學(xué)圖像配準(zhǔn)中所需的特征點(diǎn)通常為具有一定意義的點(diǎn)，而非純粹幾何意義上的點(diǎn)，河北大學(xué)附屬醫(yī)院專家在采集圖像過程中，通過對(duì)腫瘤的觸摸，標(biāo)記出涼結(jié)節(jié)的位置（如圖1c所示）.

圖1 原始的甲狀腺圖像Fig.1 Original images of thyroid

對(duì)比2種模式的圖像，2幅圖像都能顯示的解剖部位為甲狀腺和腫瘤，是具有臨床診斷意義的區(qū)域，可采用基于特征的配準(zhǔn)方法，將甲狀腺的蝴蝶形輪廓特征和腫瘤的圓形輪廓特征結(jié)合起來作為配準(zhǔn)的依據(jù).

2 SPECT 和B超圖像的配準(zhǔn)

采用基于特征的B超圖像和SPECT 圖像的配準(zhǔn)主要包括特征輪廓的提取，相似度測(cè)度，仿射變換模型的確定以及參數(shù)的優(yōu)化.

B超圖像灰度對(duì)比度低、亮度分布不均勻，本文采用基于活動(dòng)輪廓模型的圖割算法［10］得到了所需區(qū)域的分割圖（如圖2a所示）.SPECT 圖像可顯示的甲狀腺和腫瘤的邊界較為模糊，沒有一個(gè)可以視覺可觀的邊界，于是以腫瘤上標(biāo)記點(diǎn)為界限（如圖2b所示），分割出實(shí)驗(yàn)所需的甲狀腺和腫瘤的整體輪廓（如圖2c所示）.由圖2a可知，B超分割的輪廓線不圓滑，圖2c所示的SPECT 分割輪廓線較圓滑，且輪廓線不能達(dá)到完全對(duì)齊，于是不能采用基于輪廓線的配準(zhǔn)，而由于兩圖像輪廓的區(qū)域內(nèi)像素差異大，亦會(huì)造成誤配準(zhǔn)，于是我們將分割出的輪廓填充，以甲狀腺和腫瘤的輪廓區(qū)域的二值圖像分別作為參考圖像和待配準(zhǔn)圖像.填充后的結(jié)果如圖2d和圖2e所示.

圖2 分割輪廓和區(qū)域填充Fig.2 Contour image segmented and Contour region image filled

相似性測(cè)度是用來度量圖像間相似性的一種準(zhǔn)則，本文以甲狀腺和腫瘤輪廓區(qū)域二值圖像的平均互信息作為配準(zhǔn)測(cè)度.若A 代表參考圖像，B 代表待配準(zhǔn)圖像，它們之間的平均互信息I（A；B）為

其中，H（A）和H（B）分別是圖像A 和圖像B 的熵，H（AB）是它們的聯(lián)合熵，都可由A 和B 的聯(lián)合直方圖hAB求出，見式（2）

式中，A（i，j）和B（i，j）分別表示2幅圖像相同位置的坐標(biāo)點(diǎn)，n［A（i，j），B（i，j）］表示同一灰度級(jí)出現(xiàn)的次數(shù).對(duì)上式兩邊除以全部灰度值出現(xiàn)的次數(shù)和n，即可得到歸一化的聯(lián)合直方圖函數(shù)p（A，B）為

則B超圖像和SPECT 圖像輪廓區(qū)域二值圖像的互信息為

待配準(zhǔn)B超圖像是三維甲狀腺圖像的斷層圖像，由于甲狀腺內(nèi)部無相對(duì)運(yùn)動(dòng)，可以把甲狀腺近似的看做剛體運(yùn)動(dòng)，以SPECT 圖像中醫(yī)生標(biāo)出的腫瘤大小為基準(zhǔn)同比例縮放，并把圖像均裁至208＊208，使待配準(zhǔn)的圖像經(jīng)過預(yù)處理后具有相同的空間比例.設(shè)待配準(zhǔn)圖像未旋轉(zhuǎn)時(shí)中心坐標(biāo)為（a，b），旋轉(zhuǎn)后的中心坐標(biāo)為（c，d），若對(duì)圖像進(jìn)行水平平移Δx 個(gè)像素，垂直平移Δy 個(gè)像素，并以圖像中心為基點(diǎn)旋轉(zhuǎn)θ弧度，旋轉(zhuǎn)后新圖像左上角為原點(diǎn)，則SPECT 圖像和B超圖像的仿射變換模型為

B超圖像和SPECT 圖像的配準(zhǔn)問題就是尋求參數(shù)Δx，Δy，θ 的解，代入仿射變換模型T，使互信息MI（T）最大，即可得2種模式圖像的配準(zhǔn)結(jié)果.多參數(shù)優(yōu)化問題成為決定配準(zhǔn)精度的核心因素，優(yōu)化策略的選取尤為重要.

3 利用混合蛙跳算法優(yōu)化參數(shù)

混合蛙跳算法（shufled frog leaping algorithm，SFLA）是2003年由Eusuf和Lansey提出的一種基于群體智能的后啟發(fā)式計(jì)算技術(shù)［11］.它結(jié)合了基于模因進(jìn)化的模因演算法（MA，memeticalgorithm）和基于群體行為的粒子群算法（PSO，particle swarm optimization）2種群智能優(yōu)化算法的優(yōu)點(diǎn)，模擬青蛙群體尋找食物時(shí)，按族群分類進(jìn)行思想傳遞的過程，關(guān)鍵是全局搜索策略和局部深度搜索策略的完美結(jié)合使得最后的解能夠跳出局部極值點(diǎn)，向著全局最優(yōu)的方向進(jìn)行搜索，具有概念簡(jiǎn)單、調(diào)整的參數(shù)少、計(jì)算速度快、全局搜索尋優(yōu)能力強(qiáng)，易于實(shí)現(xiàn)等特點(diǎn).

首先，隨機(jī)生成U 只青蛙組成初始群體，nDim 表示變量的個(gè)數(shù)即解空間的維數(shù)，因需要尋優(yōu)的參數(shù)為Δx，Δy，θ，所以nDim＝3.第i只青蛙可以表示成U（i）＝.然后，把代入放射變換模型對(duì)待配準(zhǔn)圖像進(jìn)行變換，再計(jì)算參考圖像和變換后待配準(zhǔn)圖像的最大互信息，即為每只青蛙的適應(yīng)度，用fUi表示，并將種群內(nèi)青蛙個(gè)體按適應(yīng)度降序排列，并記錄全局最佳解Ug，即為U1.

將整個(gè)青蛙群體分成m 個(gè)族群，根據(jù)式（6）的分組方式每個(gè)族群包含n 只青蛙，并滿足關(guān)系U＝m×n.

對(duì)每個(gè)子族群進(jìn)行局部深度搜索，即對(duì)k循環(huán)進(jìn)行更新操作，根據(jù)混合蛙跳規(guī)則［12］（如圖3所示）青蛙的位置進(jìn)行更新，Ds表示移動(dòng)的距離，U′w表示更新后的青蛙，更新策略為

其中，rand（）表示0和1之間的隨機(jī)數(shù)，同時(shí)更新的距離Ds必須在可行域內(nèi).其中，Ub表示子群中的最佳解，Uw代表子群中的最差解.更新后，如果得到的解U′w優(yōu)于原來的解Uw，則令U′w＝Uw.如果沒有改進(jìn)，則用全局最優(yōu)解Ug取代原來的解Uw重復(fù)執(zhí)行更新策略.如果仍然沒有改進(jìn)，則隨機(jī)產(chǎn)生一個(gè)新的解Urandom取代原來的解Uw.當(dāng)所有子群內(nèi)部更新完成后，對(duì)子群的青蛙重新混合并排序進(jìn)行分組和對(duì)子群的內(nèi)部搜索，如此反復(fù)直到收斂到最優(yōu)解或達(dá)到最大進(jìn)化代數(shù)為止.

在選定了混合蛙跳的優(yōu)化方法后，有5個(gè)參數(shù)需要調(diào)整和確定［13］：子種群的個(gè)數(shù)、每個(gè)子群蛙的個(gè)數(shù)、每個(gè)子種群中的進(jìn)化迭代次數(shù)、允許每只青蛙移動(dòng)的最大距離、允許整個(gè)種群進(jìn)化的代數(shù)次數(shù).根據(jù)多次實(shí)驗(yàn)確定的參數(shù)可知，子群的個(gè)數(shù)和每個(gè)子群蛙的個(gè)數(shù)要選擇合適的大小，子群蛙的個(gè)數(shù)太小就會(huì)喪失局部搜索的優(yōu)點(diǎn)，還要保證初始種群的容量，容量越大，則能夠找到全局最優(yōu)的概率越大，經(jīng)實(shí)驗(yàn)測(cè)試，選用4＊5青蛙群體.子群的進(jìn)化迭代次數(shù)太小不益于子種群信息交流，太大則容易陷入局部最優(yōu)，經(jīng)實(shí)驗(yàn)測(cè)得迭代次數(shù)5至8最恰當(dāng).允許青蛙改變的距離Dmax控制的是算法進(jìn)行全局搜索的能力，所以Dmax太小，則會(huì)降低全局搜索的能力，使算法容易陷入局部最優(yōu)值，如果Dmax過大，則容易使算法錯(cuò)過最優(yōu)解，經(jīng)實(shí)驗(yàn)測(cè)試，20至34最適合.一般來說，當(dāng)循環(huán)進(jìn)化到一定次數(shù)后，代表最好解的青蛙的適應(yīng)值就不再改變了，算法即可以此為條件停止了.圖4橫坐標(biāo)代表進(jìn)化迭代的次數(shù)，縱坐標(biāo)代表全局最好青蛙的適應(yīng)值，可以看出，在45代的時(shí)候適應(yīng)值就不再改變了，并基本保持穩(wěn)定，設(shè)定迭代次數(shù)為45即可.

圖3 混合蛙跳算法迭代Fig.3 Digaram of shuffled frog leaping rule

圖4 混合蛙跳規(guī)則示意Fig.4 Iterative map of shuffled frog leaping algorithm

4 實(shí)驗(yàn)結(jié)果與分析

實(shí)驗(yàn)環(huán)境為Matlab7.1，Dell，CPU2.53GB，內(nèi)存2GB.采用的所有數(shù)據(jù)來源于河北大學(xué)附屬醫(yī)院，SPECT 圖像采自GE Infina Hawkeye 4SPECT－CT 單光子發(fā)射斷層儀，B 超圖像采自Voluson E8 三維彩色超聲診斷儀，這2幅圖像均按每4mm 一斷層，同一時(shí)期取至同一病人甲狀腺的同一層面.

為了驗(yàn)證本方法的有效性，先以2幅SPECT 圖像為例，一幅為參考圖像，另一幅圖像是人為改變水平平移量為10個(gè)像素，垂直平移量為20個(gè)像素，旋轉(zhuǎn)角度為10個(gè)角度的待配準(zhǔn)圖像.利用本方法配準(zhǔn)前后的數(shù)據(jù)如表1所示.

表1 2幅SPECT圖像實(shí)驗(yàn)結(jié)果Tab.1 Results of two SPECT images

由表1數(shù)據(jù)可以看到，基于輪廓特征點(diǎn)最大互信息的圖像配準(zhǔn)方法誤差可以保證在3個(gè)像素內(nèi)，可以達(dá)到很好的配準(zhǔn)效果，配準(zhǔn)精度較高，且有較好的穩(wěn)定性.

在調(diào)至相同迭代次數(shù)的條件下，再與其他算法的配準(zhǔn)效果加以比較，結(jié)果如圖5所示.第1行分別表示配準(zhǔn)前，粒子群算法（PSO），蟻群算法（Ant Colony），本文算法配準(zhǔn)后的輪廓疊加圖，第2行為對(duì)應(yīng)的最后融合結(jié)果.

比較3幅結(jié)果圖，可以明顯看出，粒子群算法和蟻群算法并未達(dá)到甲狀腺和腫瘤位置的對(duì)齊，而本文算法具有較好的配準(zhǔn)效果，混合蛙跳算法可以基本達(dá)到甲狀腺和腫瘤的對(duì)齊，能夠滿足臨床基本診斷需要.

對(duì)這3種算法的配準(zhǔn)時(shí)間和配準(zhǔn)精度進(jìn)行了比較，結(jié)果如表2所示.

圖5 各種算法的配準(zhǔn)比較Fig.5 Registration contrast of different algorithms

表2 3種配準(zhǔn)算法結(jié)果比較Tab.2 Results contrast of three registration algorithms

由表2數(shù)據(jù)也可證明雖然粒子群算法時(shí)間較快，但配準(zhǔn)精度不高，且實(shí)驗(yàn)結(jié)果具有隨機(jī)性；蟻群算法配準(zhǔn)精度比粒子群算法精度高，但優(yōu)化時(shí)間長(zhǎng).本文算法比蟻群算法有更高的配準(zhǔn)精度，比粒子群算法有更快的配準(zhǔn)速度.

5 結(jié)論

提出了一種基于混合蛙跳算法的甲狀腺SPECT－B超圖像配準(zhǔn)方法，通過多次實(shí)驗(yàn)證明，該方法對(duì)甲狀腺的SPECT 圖像和B超圖像配準(zhǔn)有較好的效果，不僅能夠有效地防止配準(zhǔn)結(jié)果陷入局部最優(yōu)當(dāng)中，還有很強(qiáng)的魯棒性，能夠?yàn)榧谞钕俚腟PECT 圖像和B超圖像的融合提供滿足要求的配準(zhǔn)圖像，對(duì)甲狀腺腫瘤的臨床診斷具有較高的參考價(jià)值.

［1］ SUH J W，KWON O K，SCHEINOST D，et al.CT－PET weighted image fusion for separately scanned whole body rat［J］.Med Phys，2012，39（1）：533－542.

［2］ DELIBASIS K K，ASVESTAS P A，MATSOPOULOS G K.Automatic point correspondence using an artificial immune system optimization technique for medical image registration［J］.Computerized Medical Imaging and Graphics，2011，35（1）：31－41.

［3］ FREIMAN M，WERMAN M，JOSKOWICZ L A.curvelet－based patient－specific prior for accurate multi－modal brain image rigid registration［J］.Medical Image Analysis，2011，15（1）：125－132.

［4］ KHADER M，HAMZA A B.An information－theoretic method for multimodality medical image registration［J］.Expert Systems with Applications，2012，39（5）：5548－5556.

［5］ ASHBURNER J，F(xiàn)RISTON K J.Diffeomorphic registration using geodesic shooting and Gauss－Newton optimization［J］.NeuroImage，2011，55（3）：954－967.

［6］ PAPAZOV C，BURSCHKA D.Stochastic global optimization for robust point set registration［J］.Computer Vision and Image Understanding，2011，115（12）：1598－1609.

［7］ 趙海峰，姚麗莎，羅斌，等.改進(jìn)的人工魚群算法和Powell法結(jié)合的醫(yī)學(xué)圖像配準(zhǔn)［J］.西安交通大學(xué)報(bào)，2011，45（4）：46－52.ZHAO Haifeng，YAO Lisha，LUO Bin，et al.Registration of multi－resolution medical images using a modified artificial fish－swarm algorithm combined with Powell's Method［J］.Journal of Xi'an Jiaotong University，2011，45（4）：46－52.

［8］ 譚志國(guó)，魯敏，任戈，等.匹配與姿態(tài)估計(jì)的粒子群優(yōu)化算法［J］.中國(guó)圖象圖形學(xué)報(bào)，2011，16（4）：640－646.TAN Zhiguo，LU Min，REN Ge，et al.Particle swarm optimization based pose and correspondence estimation［J］.Journal of Image and Graphics，2011，16（4）：640－646.

［9］ 李登旺，王洪君，尹勇，等.基于邊緣保護(hù)多尺度空間的醫(yī)學(xué)圖像配準(zhǔn)方法［J］.模式識(shí)別與人工智能，2011，24（1）：117－122.LI Dengwang，WANG Hongjun，YIN Yong，et al.Multiscale registration based on edge－preserved scale space for medical Images［J］.Pattern Recognition and Artificial Intelligence，2011，24（1）：117－122.

［10］ XU Ning，AHUJA N，BANSAL R.Object segmentation using graph cuts based active contours［J］.Computer Vision and Image Understanding，2007，107（3）：210－224.

［11］ EUSUFF M M，LANSEY K E.Optimization of water distribution network design using the shuffled frog leaping algorithm［J］.Journal of Water Resources Planning and Management，2003，129（3）：210－225.

［12］ 崔文華，劉曉冰，王偉，等.混合蛙跳算法研究綜述［J］.控制與決策，2012，27（4）：481－486，493.CUI Wenhua，LIU Xiaobing WANG Wei，et al.Survey on shuffled frog leaping algorithm［J］.Control and Decision，2012，27（4）：481－486，493.

［13］ ELBELTAGI E，HEGAZY T，GRIERSON D.Comparison among five evolutionary－based optimization algorithm［J］.Advanced Engineering Informaties，2001，19（1）：43－53.