基于組合相似性測度的2d-3d醫(yī)學(xué)圖像并行配準(zhǔn)方法
【專利摘要】本發(fā)明公開了一種基于組合相似性測度的2D-3D醫(yī)學(xué)圖像并行配準(zhǔn)方法。該方法首先使用CUDA并行計算模型完成DRR圖像的快速生成過程，并組合差值絕對值和SAD與模式強度PI作為新的相似性測度在GPU上進(jìn)行并行計算，最后將組合相似性測度值傳遞到CPU上采用基于細(xì)菌趨化行為的果蠅優(yōu)化算法進(jìn)行優(yōu)化來尋找最優(yōu)配準(zhǔn)參數(shù)。通過實驗對本方法性能進(jìn)行驗證表明：由于本發(fā)明方法在GPU中實現(xiàn)DRR快速生成及混合相似性測度的計算，有效地提高了本發(fā)明方法的執(zhí)行速度，同時與單一相似性測度相比，本發(fā)明采用混合相似性測度提高了配準(zhǔn)結(jié)果的精確性。
【專利說明】基于組合相似性測度的2D-3D醫(yī)學(xué)圖像并行配準(zhǔn)方法

【技術(shù)領(lǐng)域】
[0001]本發(fā)明涉及醫(yī)學(xué)圖像快速配準(zhǔn)，尤其涉及的是一種基于組合相似性測度的2D-3D醫(yī)學(xué)圖像并行配準(zhǔn)方法。

【背景技術(shù)】
[0002]圖像引導(dǎo)放療(Image-GuidedRad1therapy, IGRT)是繼三維適形放療(3DCRT)和調(diào)強放療(IMRT)之后的新的放療技術(shù)。在圖像引導(dǎo)放射治療中，為了確定分次治療的實際位置，物理師利用電子射野影像裝置等獲取X線射野影像。X線圖像反映了治療期間的實際位置，而治療前CT體數(shù)據(jù)代表治療的理想位置，通過X線射野圖像和CT體數(shù)據(jù)進(jìn)行2D-3D配準(zhǔn)，計算出治療床的偏移量，最后根據(jù)偏移量調(diào)整治療床位置，從而保證放療過程中擺位的精確性。
[0003]在2D-3D圖像配準(zhǔn)中，需要對三維CT體數(shù)據(jù)進(jìn)行重建生成DRR，由于DRR圖像的生成過程非常耗時，且在2D-3D圖像配準(zhǔn)算法中需要從不同角度生成多張DRR圖像，故DRR的生成過程直接影響了整個2D-3D配準(zhǔn)算法的效率，另外相似性測度的計算對整個配準(zhǔn)過程的運行效率和結(jié)果精確性也有重要的影響，因此如何提高2D-3D配準(zhǔn)算法的運行效率和精確性來滿足IGRT需要成為了很多學(xué)者研究的熱點。GraemeP.Penney等人對六種應(yīng)用到2D-3D醫(yī)學(xué)圖像配準(zhǔn)中的基于灰度的相似性測度進(jìn)行了比較，即使當(dāng)兩幅待配準(zhǔn)圖像中軟組織出現(xiàn)差異的情況下，模式強度和梯度差分兩種相似性測度仍舊能取得精確的配準(zhǔn)結(jié)果，并具有很好的魯棒性；W.Birkfellner等人提出了一種基于splat rendering的DRR圖像生成方法并應(yīng)用到2D-3D醫(yī)學(xué)圖像配準(zhǔn)中，在PC機(jī)上針對大小為30MB左右的體數(shù)據(jù)，繪制能用于2D-3D配準(zhǔn)的輕度模糊的DRR圖像的時間為10ms ；Daniel B.Russakoff等人提出了一種基于衰減域的DRR快速生成方法并將其應(yīng)用到2D-3D配準(zhǔn)中，其2D-3D配準(zhǔn)方法比使用基于光線投射生成DRR圖像的配準(zhǔn)算法在執(zhí)行效率和精確性方面都有一定的優(yōu)勢；OsamaDorgham等人使用八叉樹結(jié)構(gòu)來壓縮體數(shù)據(jù)，減少了投射光線和體數(shù)據(jù)內(nèi)部空間的交點的數(shù)量，從而提高了光線投射生成DRR圖像的速度，并將其應(yīng)用到2D-3D配準(zhǔn)中也能達(dá)到有效的配準(zhǔn)精度，能夠滿足校正病人擺位誤差的精度要求；以上配準(zhǔn)算法在采用不同方法改進(jìn)了 DRR圖像的生成過程，提高了配準(zhǔn)算法的效率，但不能滿足圖像引導(dǎo)放療中的實時性要求；kubiasA等人實現(xiàn)了一種基于GPU的2D-3D醫(yī)學(xué)圖像配準(zhǔn)算法，該算法將DRR生成和相似性測度計算都在GPU上完成，與傳統(tǒng)的基于CPU的配準(zhǔn)算法相比大大地提高了配準(zhǔn)速度，但由于受流水線的限制，GPU的加速能力仍舊沒有得到充分的發(fā)揮；0samaM.Dorgham等人在支持CUDA的GPU實現(xiàn)了 DRR圖像的并行加速生成，同時對光線投射過程進(jìn)行壓縮采樣等處理來進(jìn)一步提高DRR生成速度，最后將該DRR加速生成方法應(yīng)用到了 2D-3D醫(yī)學(xué)圖像配準(zhǔn)中，并且取得了臨床放療可接受的精度，但是該方法并沒有將相似性測度計算過程也搬到GPU上運行。
[0004]醫(yī)學(xué)圖像配準(zhǔn)是將一幅圖像的像素空間位置進(jìn)行變換使其與另一幅圖像像素的空間位置對齊的過程，是醫(yī)學(xué)圖像融合的前提和基礎(chǔ)，在病情診斷、圖像引導(dǎo)放療及手術(shù)實施中具有重要的應(yīng)用價值。在圖像引導(dǎo)放射治療應(yīng)用中2D-3D醫(yī)學(xué)圖像配準(zhǔn)是將放射治療前制定放療計劃用的CT體數(shù)據(jù)集與放射治療中代表病人當(dāng)前治療擺位的X線圖像進(jìn)行配準(zhǔn)，得到相應(yīng)的配準(zhǔn)參數(shù)來糾正放療實施過程中的擺位誤差。整個配準(zhǔn)過程包括三個主要的步驟:DRR圖像的生成、DRR圖像與X線圖像之間的相似度計算和相似性測度函數(shù)尋優(yōu)。
[0005]在2D-3D圖像配準(zhǔn)過程中，將實際拍攝的X線圖像作為參考圖像，計算機(jī)生成的DRR圖像作為浮動圖像，通過改變對CT體數(shù)據(jù)的虛擬透視位置及角度，獲取相應(yīng)的DRR圖像，并逐個與真實X線圖像進(jìn)行相似性測度的比較直到出現(xiàn)最優(yōu)解，此時CT數(shù)據(jù)的六自由度空間參數(shù)即為所求。2D-3D醫(yī)學(xué)圖像配準(zhǔn)的數(shù)學(xué)表示如下:
[0006]Γ* = argmax^lX^ (VolumeData^(I)

T
[0007]式中S是相似性測度函數(shù)，T是包含了六個自由度的空間變換矩陣，X是放療實施過程中代表病人當(dāng)前位置的X線圖像，VolumeData為三維CT體數(shù)據(jù)，T(VolumeData)即為對三維CT體數(shù)據(jù)進(jìn)行矩陣T所示的剛性變換后得到的DRR圖像，整個配準(zhǔn)過程就是尋求最佳空間變換T使得相似性測度函數(shù)S達(dá)到最大的過程。2D-3D配準(zhǔn)是剛體坐標(biāo)變換，研究對象的位置和方向變換采用六參數(shù)矢量來描述M7；,，即六自由度。其中，
Tx, Ty, Tz分別表示沿著X、Y、Z軸方向的平移量，Rx, Ry, Rz分別表示圍繞Χ、Υ、Ζ軸的旋轉(zhuǎn)量。


【發(fā)明內(nèi)容】

[0008]本發(fā)明針對圖像引導(dǎo)放療對2D-3D配準(zhǔn)算法的實時性和高精度的要求，且由于DRR生成過程的運行速度嚴(yán)重制約了 2D-3D醫(yī)學(xué)圖像配準(zhǔn)算法的效率，首先使用CUDA并行計算模型設(shè)計實現(xiàn)了 2D-3D醫(yī)學(xué)圖像配準(zhǔn)算法的DRR并行快速生成過程，并在此基礎(chǔ)上，由于SAD計算相對簡單且模式強度PI能有效提高配準(zhǔn)精度，本發(fā)明提出了一種基于SAD和PI的混合相似性測度函數(shù)并將該混合相似性測度函數(shù)也在CUDA上并行實現(xiàn)，最后在CPU上使用全局尋優(yōu)能力強且計算方便的細(xì)菌趨化果蠅優(yōu)化算法來進(jìn)行整個配準(zhǔn)過程的全局尋優(yōu)，從而在保證算法運行效率的前提下提高了配準(zhǔn)結(jié)果的精確性。
[0009]本發(fā)明的技術(shù)方案如下:
[0010]一種基于組合相似性測度的2D-3D醫(yī)學(xué)圖像并行配準(zhǔn)方法，采用基于SAD和PI相結(jié)合的相似性測度Mmw作為2D-3D醫(yī)學(xué)圖像配準(zhǔn)的目標(biāo)函數(shù)；其表達(dá)如式(I)所示:
[0011]
Mttew=~{-SAD + Pl)(I)
[0012]在式(I)中，平方差和SAD的計算方法如式(2)所示:
[0013]
I N-X M-1
SAD = ――Σ ΣI沖’ j.) -Ft {U7)1(2)
iV XjW 2=0 J=Q
[0014]在式⑵中，F(xiàn)T(i，j)是進(jìn)行空間變換后的浮動圖像，N表示圖像的寬度，M表示圖像的高度。在式(I)中，模式強度PI的計算方法如式(3)所示:
[0015]^=ΣΣ , h σ: , ,、、2(3)
ι卜σ +(/'/ (,，./)-々/ η))
[0016]其中:d2= (1-v)2+ (j-w)2 (4)
[0017]式⑶中，常數(shù)σ用來消除噪聲的干擾，通常取σ = 10，半徑r通常取3、4或5。式⑶中Id為將X線圖像Ix與DRR圖像Si.相減得到差值圖像，如式(5)所示:
[0018]Id = Ir-Slf (5)
[0019]所述的基于組合相似性測度的2D-3D醫(yī)學(xué)圖像并行配準(zhǔn)方法，包括以下步驟:
[0020]Stepl:獲取X線圖像作為配準(zhǔn)過程的參考圖像，獲取醫(yī)學(xué)圖像CT序列作為配準(zhǔn)過程的浮動圖像；
[0021]Step2:對X線圖像進(jìn)行去噪預(yù)處理，將醫(yī)學(xué)圖像CT序列圖像構(gòu)造生成三維體數(shù)據(jù)；
[0022]Step3:將處理后的X線圖像數(shù)據(jù)和三維體數(shù)據(jù)裝載到GPU中；
[0023]Step4:確定初始變換參數(shù)P。= (Tx, Ty, Tz, Rx, Ry, Rz)和CT值與衰減系數(shù)對應(yīng)關(guān)系的查找表，以及SAD、SFD參數(shù)；
[0024]Step5:根據(jù)當(dāng)前變換參數(shù)P在GPU中由多線程并行執(zhí)行內(nèi)核過程完成對三維體數(shù)據(jù)進(jìn)行光線投射產(chǎn)生DRR圖像；
[0025]Step6:在GPU中根據(jù)式(5)并行計算DRR圖像與X線圖像的混合相似性測度，并將計算得到的混合相似性測度值傳遞給CPU ；
[0026]Step7:在CPU中利用基于細(xì)菌趨化行為的果蠅優(yōu)化算法對相似性測度函數(shù)進(jìn)行尋優(yōu)；
[0027]StepS:如果當(dāng)前相似性測度函數(shù)已達(dá)到最優(yōu)，則輸出變換參數(shù)P作為最終的配準(zhǔn)結(jié)果；否則根據(jù)優(yōu)化策略來更新當(dāng)前的變換參數(shù)P，轉(zhuǎn)至Step5執(zhí)行。
[0028]所述St印7中的果蠅優(yōu)化算法優(yōu)化流程如下:
[0029]I)初始化果蠅種群規(guī)模SizeofPop、最大迭代次數(shù)MaxGen、初始迭代次數(shù)j = O及所有果蠅個體的位置Pi, Q = (τχ, Ty, Tz, Rx, Ry, Rz)等參數(shù)；
[0030]2)由于果蠅個體主要利用嗅覺來搜索食物，針對每個果蠅個體隨機(jī)賦予搜索食物的方向和距離，并采用式(7)在果蠅群體位置的基礎(chǔ)上添加隨機(jī)距離P來生成果蠅群體，其中Pu表示第i個果蠅在第j次迭代中的位置參數(shù)；
[0031]Pi； j = Pi, j+RandomDist (7)
[0032]3)當(dāng)前迭代次數(shù)j是否小于最大迭代次數(shù)MaxGen，如果是則執(zhí)行4)繼續(xù)迭代，如果否則終止優(yōu)化過程，轉(zhuǎn)至10)；
[0033]4)將每個果蠅個體的位置參數(shù)Pi,」=(Tx, Ty, Tz, Rx, Ry, Rz)傳遞到GPU中生成DRR圖像，生成后并在GPU中根據(jù)式(10)計算X線圖像和DRR圖像的混合相似性測度函數(shù)值Mi, j,將計算得到的相似性測度函數(shù)值Μ。傳回給CPU，其中Μ。表示第i個果蠅個體在第j次迭代中的相似性測度值；
[0034]5)根據(jù)式⑶對所有果蠅個體位置的對應(yīng)的相似性測度值Miij求出其最小值和最大值，從而針對所有果蠅個體分別找出當(dāng)前迭代中的混合相似性測度最優(yōu)的果蠅和混合相似性測度最差的果蠅；
[beslMeasure besifndex] = mini Μ..)(8)
[0035]]/ 、
\\vorsiMeasure worsiIndex\ = max( , I
[0036]6)根據(jù)式(9)記錄并保留當(dāng)前最佳相似性測度值bestMeasure與其Pb坐標(biāo),當(dāng)前最差混合相似性測度值worstMeasure與其對應(yīng)的Pw坐標(biāo)值；
MeasureBesi = bestMeasure
P=P
廠Ibbestlndex.1ζυλ、
[0037]<(9)
MeasureWorsi = worstMeasure

P - P

w worstlndex.j
[0038]7)根據(jù)式(10)來計算果蠅群體的平均相似性測度值Mavg，然后根據(jù)平均相似性測度值Mavg利用式(11)來計算相似性測度方差σ2，并根據(jù)相似性測度方差σ 2來判斷果蠅群體朝最優(yōu)個體的方向移動還是朝遠(yuǎn)離最差個體的方向移動；

ISLeOfPop
[0039]V/, =——— Σ V/(丨(”

SizeOfPop ,--1

1SneOfPon
_] σ2§ (H)2⑴)
[0041]8)根據(jù)混合相似性測度方差來判斷果蠅群體執(zhí)行吸引操作還是逃離操作，從而來豐富種群多樣性。若。2古0，則按式(12)更新所有果蠅個體位置Piij;否則按式(13)更新所有果蠅個體位置Pu;
[0042]Pij = PjPb (12)
[0043]Pij = Pij-Pw (13)
[0044]9)當(dāng)前迭代次數(shù)自增1，即j = j+1，轉(zhuǎn)至3)執(zhí)行；
[0045]10)輸出混合相似性測度函數(shù)的全局最優(yōu)值及其對應(yīng)的果蠅個體的位置，該位置參數(shù)即為2D-3D配準(zhǔn)算法所求的最佳配準(zhǔn)變換參數(shù)。
[0046]該算法首先使用CUDA并行計算模型完成DRR圖像的快速生成過程，并組合差值絕對值和SAD與模式強度PI作為新的相似性測度在GPU上進(jìn)行并行計算，最后將組合相似性測度值傳遞到CPU上采用基于細(xì)菌趨化行為的果蠅優(yōu)化算法進(jìn)行優(yōu)化來尋找最優(yōu)配準(zhǔn)參數(shù)。通過實驗對算法性能進(jìn)行驗證表明:由于本發(fā)明算法在GPU中實現(xiàn)DRR快速生成及混合相似性測度的計算，有效地提高了本發(fā)明算法的執(zhí)行速度，同時與單一相似性測度相比，本發(fā)明采用混合相似性測度提高了配準(zhǔn)結(jié)果的精確性。

【專利附圖】

【附圖說明】
[0047]圖1是基于CUDA實現(xiàn)DRR生成算法的主要流程；
[0048]圖2是本發(fā)明2D-3D配準(zhǔn)算法的主要流程。

【具體實施方式】
[0049]以下結(jié)合具體實施例，對本發(fā)明進(jìn)行詳細(xì)說明。
[0050]實施例1
[0051]CUDA是一種不需要借助圖形學(xué)API就能使用類C語言進(jìn)行通用計算的開發(fā)環(huán)境和軟件體系結(jié)構(gòu)，它的出現(xiàn)有效地發(fā)揮了 GPU的并行計算性能，針對本質(zhì)具有并行特征的算法能夠有效地提高算法的執(zhí)行效率。在基于光線投射的DRR生成算法中，每束光線穿透體數(shù)據(jù)來模擬X線入射人體后的衰減過程相互獨立，具有良好的并行性，因此該過程可通過設(shè)計在GPU上執(zhí)行的內(nèi)核函數(shù)來實現(xiàn)。
[0052]由于生成DRR的輸入數(shù)據(jù)是醫(yī)學(xué)CT圖像序列，其分辨率通常為512*512，且進(jìn)行2D-3D配準(zhǔn)過程中需要使用的DRR的分辨率也為512*512，因此可從光源發(fā)出512*512條射線，每條射線對應(yīng)DRR圖像的一個像素，將一條射線穿透體數(shù)據(jù)進(jìn)行采樣并累加計算CT值的過程封裝在內(nèi)核函數(shù)中由一個線程完成，最終所有射線穿過體數(shù)據(jù)進(jìn)行光線衰減產(chǎn)生二維的DRR圖像。為了充分發(fā)揮GPU的并行計算能力，可將整個DRR圖像屏幕作為一個網(wǎng)格，并將屏幕分為若干塊區(qū)域，屏幕上每個像素點對應(yīng)一個線程的執(zhí)行結(jié)果，這樣可保證整個并行計算過程充分利用CUDA的三層體系結(jié)構(gòu)。因此根據(jù)上述分析確定線程分配方案:網(wǎng)格大小為32*32，線程塊大小為16*16，每條光線在X和Y方向上都有32*16 = 512條射線，每條光線路徑上的采樣點計算和CT值累加計算過程由一個線程獨立執(zhí)行內(nèi)核函數(shù)完成。在完成線程分配和數(shù)據(jù)存儲之后，根據(jù)每條光線的投射采樣和CT值衰減累加計算的過程來設(shè)計內(nèi)核函數(shù)，在CUDA中多線程來執(zhí)行內(nèi)核函數(shù)實現(xiàn)光線投射過程的并行計算?；贑UDA實現(xiàn)DRR生成算法的主要流程如圖1所示。
[0053]相似性測度函數(shù)是醫(yī)學(xué)圖像配準(zhǔn)過程中非常重要的一步，相似性測度的選擇直接影響到配準(zhǔn)結(jié)果的精確性，因此尋找一種良好的相似性測度函數(shù)對整個2D-3D配準(zhǔn)非常重要。在2D-3D醫(yī)學(xué)圖像配準(zhǔn)算法中，常用的相似性測度有差值平方和SSD、差值絕對值和SAD、互信息M1、歸一化互信息匪1、模式強度P1、梯度相關(guān)系數(shù)GC及梯度差分⑶等。由于SAD只需要兩幅圖像對應(yīng)像素求差，不需要進(jìn)行導(dǎo)數(shù)運算，計算十分簡單，在同模態(tài)圖像配準(zhǔn)中使用比較廣泛，且在2D-3D配準(zhǔn)中的兩幅待配準(zhǔn)圖像X線圖像和DRR圖像的成像原理相似，可以看做同模態(tài)圖像，故可采用SAD作為相似性測度函數(shù)；另外模式強度PI考慮了兩幅待配準(zhǔn)圖像的像素灰度的同時還考慮了圖像的空間信息，配準(zhǔn)精度高且計算相對方便，因此為了獲取更好的配準(zhǔn)精度，結(jié)合待配準(zhǔn)圖像的空間信息和像素強度信息，在本發(fā)明算法中采用了基于SAD和PI相結(jié)合的相似性測度Mmw作為本發(fā)明2D-3D醫(yī)學(xué)圖像配準(zhǔn)的目標(biāo)函數(shù)。其數(shù)學(xué)表達(dá)如式(2)所示:
[0054]
W =^(-SAD + PI)(2)
2
[0055]在式(2)中，平方差和SAD的計算方法如式(3)所示:
[0056]
I N-XM-X
[0057]在式(3)中，F(xiàn)T(i，j)是進(jìn)行空間變換后的浮動圖像，N表示圖像的寬度，M表示圖像的高度。在式(2)中，模式強度PI的計算方法如式(4)所示:
[0058]


2
Ρ‘ΣΣ , (ι (.σ, “、ν⑷
廣'一 +(/,'(Μ)-7"、'、'1'))
[0059]其中:d2= (1-v)2+(j-w)2 (5)
[0060]式(4)中，常數(shù)σ用來消除噪聲的干擾，通常取σ = 10，半徑r通常取3、4或5。式⑷中Id為將X線圖像Ix與DRR圖像Si.相減得到差值圖像，如式(6)所示:
[0061]Id = Ir-Slf (6)
[0062]從式(4)可以看出:當(dāng)該模式趨近于零時，模式強度測度的值趨近于一個最大值；當(dāng)該模式增大時，模式強度測度趨于零。且模式強度采用函數(shù)I/(1+x2)的形式進(jìn)行構(gòu)造，根據(jù)該函數(shù)的漸進(jìn)性質(zhì)可知，當(dāng)兩幅圖像的灰度差別很大時對該測度影響很小，從而當(dāng)參考圖像和浮動圖像之間具有較大灰度差別時模式強度測度也能體現(xiàn)出良好的適應(yīng)性。
[0063]新的相似性測度Mmw在保持SAD計算相對簡單的前提下引入了模式強度PI，可以有效地提高配準(zhǔn)結(jié)果的精確性。本發(fā)明算法在對CT序列圖像進(jìn)行變換并生成DRR圖像后，根據(jù)X線圖像與DRR圖像在GPU上通過執(zhí)行內(nèi)核函數(shù)來計算混合相似性測度Mnew,將計算得到的混合相似性測度值傳到CPU上，供優(yōu)化算法完成配準(zhǔn)過程的全局優(yōu)化。
[0064]為了提高2D-3D醫(yī)學(xué)圖像配準(zhǔn)算法執(zhí)行效率的同時兼顧配準(zhǔn)結(jié)果的精度，本發(fā)明提出了一種基于CUDA并行計算模型和組合相似性測度的2D-3D醫(yī)學(xué)圖像快速配準(zhǔn)算法。本發(fā)明算法首先在支持CUDA并行計算模型的GPU上完成DRR圖像的生成過程；并鑒于SAD計算簡單的特點與模式強度PI能有效提高配準(zhǔn)精度的優(yōu)勢，將兩者進(jìn)行組合作為新的相似性測度，并在GPU上完成參考圖像X圖像和DRR圖像的組合相似性測度的計算；最后將混合相似性測度值傳遞到CPU上采用基于細(xì)菌趨化行為的果蠅優(yōu)化算法來尋找最佳配準(zhǔn)變換參數(shù)。
[0065]本發(fā)明2D-3D醫(yī)學(xué)圖像配準(zhǔn)算法的主要流程如下:
[0066]Stepl:獲取X線圖像作為配準(zhǔn)過程的參考圖像，獲取醫(yī)學(xué)圖像CT序列作為配準(zhǔn)過程的浮動圖像；
[0067]Step2:對X線圖像進(jìn)行去噪預(yù)處理，將醫(yī)學(xué)圖像CT序列圖像構(gòu)造生成三維體數(shù)據(jù)；
[0068]Step3:將處理后的X線圖像數(shù)據(jù)和三維體數(shù)據(jù)裝載到GPU中；
[0069]St印4:確定初始變換參數(shù)P。= (Tx, Ty, Tz, Rx, Ry, Rz)和CT值與衰減系數(shù)對應(yīng)關(guān)系的查找表，以及SAD、SFD參數(shù)；
[0070]Step5:根據(jù)當(dāng)前變換參數(shù)P在GPU中由多線程并行執(zhí)行內(nèi)核過程完成對三維體數(shù)據(jù)進(jìn)行光線投射產(chǎn)生DRR圖像；
[0071]Step6:在GPU中根據(jù)式(6)并行計算DRR圖像與X線圖像的混合相似性測度，并將計算得到的混合相似性測度值傳遞給CPU ；
[0072]Step7:在CPU中按照實施例2的優(yōu)化流程利用基于細(xì)菌趨化行為的果蠅優(yōu)化算法對相似性測度函數(shù)進(jìn)行尋優(yōu)；
[0073]StepS:如果當(dāng)前相似性測度函數(shù)已達(dá)到最優(yōu)，則輸出變換參數(shù)P作為最終的配準(zhǔn)結(jié)果；否則根據(jù)優(yōu)化策略來更新當(dāng)前的變換參數(shù)P，轉(zhuǎn)至Step5執(zhí)行。
[0074]本發(fā)明2D-3D配準(zhǔn)算法的主要流程如圖2所示。
[0075]實施例2
[0076]醫(yī)學(xué)圖像配準(zhǔn)本質(zhì)上是將相似性測度函數(shù)作為目標(biāo)函數(shù)尋求最優(yōu)解的過程。由于很多相似性測度函數(shù)都具有多個局部極值，因此在醫(yī)學(xué)圖像配準(zhǔn)過程中使用一種具有很強全局尋優(yōu)能力的智能優(yōu)化算法，對尋找全局最優(yōu)值和提高配準(zhǔn)結(jié)果精度具有重要的意義。
[0077]果蠅優(yōu)化算法(FruitFlyAlgorithm, F0A)是由中國臺灣學(xué)者潘文超教授提出的一種模擬果蠅覓食行為的全局優(yōu)化算法，該算法具有參數(shù)少、運行時間短、容易理解和實現(xiàn)等特點，被廣泛應(yīng)用于科學(xué)和工程很多領(lǐng)域。但是針對醫(yī)學(xué)圖像配準(zhǔn)中的具有多極值特征的相似性測度函數(shù)進(jìn)行優(yōu)化時，F(xiàn)OA與遺傳算法GA、粒子群算法PSO —樣容易陷入局部最優(yōu)，從而降低了醫(yī)學(xué)圖像配準(zhǔn)的精確性。韓俊英等將細(xì)菌趨化過程中的吸引與排斥行為引入到FOA中，提出了基于細(xì)菌趨化的果蠅優(yōu)化算法(BCFOA)，有效地解決了容易陷入局部最優(yōu)的問題。
[0078]本發(fā)明在2D-3D配準(zhǔn)算法的優(yōu)化過程中采用基于細(xì)菌趨化的果蠅優(yōu)化算法BCFOA來實現(xiàn)混合相似性測度函數(shù)尋優(yōu)，主要優(yōu)化流程如下:
[0079]I)初始化果蠅種群規(guī)模SizeofPop、最大迭代次數(shù)MaxGen、初始迭代次數(shù)j = O及所有果蠅個體的位置Pi, Q = (τχ, Ty, Tz, Rx, Ry, Rz)等參數(shù)；
[0080]2)由于果蠅個體主要利用嗅覺來搜索食物，針對每個果蠅個體隨機(jī)賦予搜索食物的方向和距離，并采用式(7)在果蠅群體位置的基礎(chǔ)上添加隨機(jī)距離來生成果蠅群體，其中Pu表示第i個果蠅在第j次迭代中的位置參數(shù)；
[0081]Pi,」=Pi； j+RandomDist (7)
[0082]3)當(dāng)前迭代次數(shù)j是否小于最大迭代次數(shù)MaxGen，如果是則執(zhí)行4)繼續(xù)迭代，如果否則終止優(yōu)化過程，轉(zhuǎn)至10)；
[0083]4)將每個果蠅個體的位置參數(shù)Pi,」=(Tx, Ty, Tz, Rx, Ry, Rz)傳遞到GPU中生成DRR圖像，生成后并在GPU中根據(jù)式(10)計算X線圖像和DRR圖像的混合相似性測度函數(shù)值Mi, j,將計算得到的相似性測度函數(shù)值Μ。傳回給CPU，其中Μ。表示第i個果蠅個體在第j次迭代中的相似性測度值；
[0084]5)根據(jù)式⑶對所有果蠅個體位置的對應(yīng)的相似性測度值Miij求出其最小值和最大值，從而針對所有果蠅個體分別找出當(dāng)前迭代中的混合相似性測度最優(yōu)的果蠅和混合相似性測度最差的果蠅；
[bestMeasure be.sl Index] = mi η I Μ..)(8)
[0085]/ 、
\y\-orsiMeasure worsilndex\ = max (Mj j j
[0086]6)根據(jù)式(9)記錄并保留當(dāng)前最佳相似性測度值bestMeasure與其Pb坐標(biāo),當(dāng)前最差混合相似性測度值worstMeasure與其對應(yīng)的Pw坐標(biāo)值；
MeasureBesi = bestMeasure
P=P
「Ib bestlndex、j
[0087]<(9)
MeasiireWorsi = worstMeasure
P=P

w worstlndex yj
[0088]7)根據(jù)式(10)來計算果蠅群體的平均相似性測度值Mavg，然后根據(jù)平均相似性測度值Mavg利用式(11)來計算相似性測度方差σ2，并根據(jù)相似性測度方差σ 2來判斷果蠅群體朝最優(yōu)個體的方向移動還是朝遠(yuǎn)離最差個體的方向移動；

【權(quán)利要求】
1.一種基于組合相似性測度的2D-3D醫(yī)學(xué)圖像并行配準(zhǔn)方法，其特征在于， 采用基于SAD和PI相結(jié)合的相似性測度Mnew作為2D-3D醫(yī)學(xué)圖像配準(zhǔn)的目標(biāo)函數(shù)；其表達(dá)如式(I)所示: Maey =-[-SAD PI)(I) 在式U」中，平方差和SAD的計算方法如式(2)所不:
I N-1M-1 SAD =
iv XiW ,.=ο J=O 在式(2)中，F(xiàn)T(i，j)是進(jìn)行空間變換后的浮動圖像，N表示圖像的寬度，M表示圖像的高度；在式(I)中，模式強度PI的計算方法如式(3)所示:P -yy-----⑴
r/-σ- + (Iil (/，/.) — /(/ (V’，U' )j 其中-A2 = (1-v)2+(j-w)2 (4) 式(3)中，常數(shù)σ用來消除噪聲的干擾，通常取σ = 10，半徑r通常取3、4或5;式(3)中Id為將X線圖像Ix與DRR圖像SIddk相減得到差值圖像，如式(5)所示:
Id = Ir-Slf (5) 所述的基于組合相似性測度的2D-3D醫(yī)學(xué)圖像并行配準(zhǔn)方法，包括以下步驟: Stepl:獲取X線圖像作為配準(zhǔn)過程的參考圖像，獲取醫(yī)學(xué)圖像CT序列作為配準(zhǔn)過程的浮動圖像； Step2:對X線圖像進(jìn)行去噪預(yù)處理，將醫(yī)學(xué)圖像CT序列圖像構(gòu)造生成三維體數(shù)據(jù)； Step3:將處理后的X線圖像數(shù)據(jù)和三維體數(shù)據(jù)裝載到GPU中； Step4:確定初始變換參數(shù)Ptl = (Tx, Ty, Tz, Rx, Ry, Rz)和CT值與衰減系數(shù)對應(yīng)關(guān)系的查找表，以及SAD、SFD參數(shù)； Step5:根據(jù)當(dāng)前變換參數(shù)P在GPU中由多線程并行執(zhí)行內(nèi)核過程完成對三維體數(shù)據(jù)進(jìn)行光線投射產(chǎn)生DRR圖像； Step6:在GPU中根據(jù)式(5)并行計算DRR圖像與X線圖像的混合相似性測度，并將計算得到的混合相似性測度值傳遞給CPU ； Step7:在CPU中利用基于細(xì)菌趨化行為的果蠅優(yōu)化算法對相似性測度函數(shù)進(jìn)行尋優(yōu)；StepS:如果當(dāng)前相似性測度函數(shù)已達(dá)到最優(yōu)，則輸出變換參數(shù)P作為最終的配準(zhǔn)結(jié)果；否則根據(jù)優(yōu)化策略來更新當(dāng)前的變換參數(shù)P，轉(zhuǎn)至Step5執(zhí)行。
2.根據(jù)權(quán)利要求1所述的方法，其特征在于，所述Step7中的果蠅優(yōu)化算法優(yōu)化流程如下: 1)初始化果蠅種群規(guī)模SizeofPop、最大迭代次數(shù)MaxGen、初始迭代次數(shù)j= O及所有果蠅個體的位置Pi, ο = (Tx, Ty, Tz, Rx, Ry, Rz)等參數(shù)； 2)由于果蠅個體主要利用嗅覺來搜索食物，針對每個果蠅個體隨機(jī)賦予搜索食物的方向和距離，并采用式(6)在果蠅群體位置的基礎(chǔ)上添加隨機(jī)距離來生成果蠅群體，其中Pu表示第i個果蠅在第j次迭代中的位置參數(shù)；
Pi； j = Pi, j+RandomDist (6) 3)當(dāng)前迭代次數(shù)j是否小于最大迭代次數(shù)MaxGen，如果是則執(zhí)行4)繼續(xù)迭代，如果否則終止優(yōu)化過程，轉(zhuǎn)至10)； 4)將每個果蠅個體的位置參數(shù)Pu= (Tx, Ty, Tz, Rx, Ry, Rz)傳遞到GPU中生成DRR圖像，生成后并在GPU中根據(jù)式(I)計算X線圖像和DRR圖像的混合相似性測度函數(shù)值Miij，將計算得到的相似性測度函數(shù)值Μ。傳回給CPU，其中Μ。表示第i個果蠅個體在第j次迭代中的相似性測度值； 5)根據(jù)式(7)對所有果蠅個體位置的對應(yīng)的相似性測度值Μ。求出其最小值和最大值，從而針對所有果蠅個體分別找出當(dāng)前迭代中的混合相似性測度最優(yōu)的果蠅和混合相似性測度最差的果蠅；
[hesiMeasure he.sifndcx] = m.1.n (M;.j)(7)
<
[worstMeasure worsilndex^ = max ,) 6)根據(jù)式⑶記錄并保留當(dāng)前最佳相似性測度值bestMeasure與其Pb坐標(biāo)，當(dāng)前最差混合相似性測度值worstMeasure與其對應(yīng)的Pw坐標(biāo)值；
MeasureBesl = best M easure
P=P
b bestlndex.j(8)
MeasureWors1- worstMeasure
p=p
w worstlndex, j 7)根據(jù)式(9)來計算果蠅群體的平均相似性測度值Mavg，然后根據(jù)平均相似性測度值Mavg利用式(10)來計算相似性測度方差σ2，并根據(jù)相似性測度方差σ 2來判斷果蠅群體朝最優(yōu)個體的方向移動還是朝遠(yuǎn)離最差個體的方向移動；
ISizeOfPop M =- ^ M-.(9)
噸 SizeOfPop il
ISizeOfPop0 σ2 =- V (Μ - )':(10)
SizeOfPop tT V '' g), 8)根據(jù)混合相似性測度方差來判斷果蠅群體執(zhí)行吸引操作還是逃離操作，從而來豐富種群多樣性；若σ 2古0，則按式(11)更新所有果蠅個體位置Pu ;否則按式(12)更新所有果蠅個體位置Pu;
Pi, j = PiiJ+Pb (H)
Pi, J = Pi, J-Pw (12) 9)當(dāng)前迭代次數(shù)自增1，即j= j+1，轉(zhuǎn)至3)執(zhí)行； 10)輸出混合相似性測度函數(shù)的全局最優(yōu)值及其對應(yīng)的果蠅個體的位置，該位置參數(shù)即為2D-3D配準(zhǔn)算法所求的最佳配準(zhǔn)變換參數(shù)。
【文檔編號】G06T7/00GK104134210SQ201410351843
【公開日】2014年11月5日 申請日期:2014年7月22日 優(yōu)先權(quán)日:2014年7月22日
【發(fā)明者】黨建武, 杜曉剛, 王陽萍, 楊景玉, 王松, 陳永, 王冰 申請人:蘭州交通大學(xué)