改進(jìn)的數(shù)字重建影像快速生成算法

羅博博，贠照強(qiáng)，盧振泰，林仁回，張 娟，馮前進(jìn)

（南方醫(yī)科大學(xué) 生物醫(yī)學(xué)工程學(xué)院，廣東 廣州510515）

0 引 言

數(shù) 字 重 建 影 像 （digitally reconstructed radiographs，DRRs）是對CT圖像進(jìn)行模擬投影方式產(chǎn)生的虛擬X線圖像，廣泛應(yīng)用于計算機(jī)輔助外科手術(shù)、圖像引導(dǎo)介入治療和圖像引導(dǎo)放射治療等領(lǐng)域［1］。特別地，在骨科臨床，骨折手術(shù)治療涉及的解剖結(jié)構(gòu)相對復(fù)雜，常規(guī)手術(shù)風(fēng)險較大，應(yīng)用導(dǎo)航技術(shù)可提高手術(shù)的安全性和準(zhǔn)確性。DRRs的重建已經(jīng)成為臨床骨科手術(shù)導(dǎo)航中的一項十分關(guān)鍵的關(guān)鍵技術(shù)。

目前生成DRRs的算法主要包括：光線投影 （ray casting）［2］、搖晃拋雪球 （wobbled splatting）［3］、剪切變形分解（shear－warp factorization）［4］、跨 圖 形 學(xué) （the transgraph）［5］、自適應(yīng)蒙特卡洛體繪制 （adaptive monte carlo volume rendering）［6］、衰 減 場 （attenuation fields）［7］、逐 步 衰 減 場 （progressive attenuation fields）［8］等算法。DRRs的生成是一個耗時相當(dāng)長的過程，對于N×N×N的圖像數(shù)據(jù)，傳統(tǒng)的算法，如投影 （ray casting）算法、拋雪球 （SPLATTING）算法和剪切變形 （SHEAR－warp）等算法的計算復(fù)雜度為O（N3）［9］，為滿足實時性要求，Daniel B.Russakoff等提出了基于光場的數(shù)字重建影像方法，將光場理論引入DRRs的形成過程，大幅度提高了DRRs的生成速度，但該方法內(nèi)存消耗量大、建立的光場數(shù)據(jù)有大量冗余、預(yù)處理時間長；D.Ruijters等提出了基于圖形處理器 GPU （graphic processing unit）的DRRs生成方法［10］，大幅提升了DRRs圖像的生成效率。

而在臨床骨科實時圖像引導(dǎo)手術(shù)中，基于灰度的2D／3D醫(yī)學(xué)圖像配準(zhǔn)給臨床醫(yī)生提供一個到人體內(nèi)部虛擬的、非侵入式的窗口，使醫(yī)生能夠看到一個肉眼無法直接看到的解剖與手術(shù)器械的三維空間相對位置關(guān)系，使得手術(shù)精確、安全進(jìn)行。為了快速、準(zhǔn)確獲取病灶和手術(shù)器械的三維空間相對位置關(guān)系，常要生成數(shù)百張DRRs，并做到實時配準(zhǔn)。因此快速DRRs生成方法的研究成為骨科手術(shù)導(dǎo)航系統(tǒng)的關(guān)鍵問題。

降采樣是對過采樣或正常采樣的數(shù)據(jù)進(jìn)行二次抽取或者降低頻率二次采樣數(shù)據(jù)的過程，這樣可以減少計算量，節(jié)約計算時間，是實時處理中的一種常用的方法。通用并行計算架構(gòu) （compute unified device architecture）是一種將GPU作為數(shù)據(jù)并行計算設(shè)備的軟硬件體系，可以對數(shù)據(jù)進(jìn)行高度的并行計算，在處理圖形數(shù)據(jù)和復(fù)雜算法方面擁有比CPU更高的效率。為保證臨床骨科手術(shù)安全、順利進(jìn)行，采用只含骨性結(jié)構(gòu)的DRRs實現(xiàn)術(shù)中實時圖像引導(dǎo)，本文結(jié)合多種加速技術(shù)，提出一種改進(jìn)的算法——基于光線投影算法的DRRs生成原理，采用統(tǒng)一并行計算架構(gòu)CUDA加速，再引入大步長與只采樣骨性結(jié)構(gòu)相結(jié)合的降采樣，快速生成只含骨性結(jié)構(gòu)的DRRs，并采用B.SrinivasaReddy和 B.N.Chatterji提 出 的 基 于 Fourier－Meilin變換和相位相關(guān)相似性測度的醫(yī)學(xué)圖像配準(zhǔn)方法將通過改進(jìn)算法得到的DRRs和臨床X線圖像進(jìn)行醫(yī)學(xué)圖像配準(zhǔn)驗證。實驗結(jié)果表明，只采樣骨性結(jié)構(gòu)同時不超過3倍 （甚至5倍）原步長采樣時，改進(jìn)的算法能快速生成只含骨性結(jié)構(gòu)的DRRs，也能使得獲取的DRRs和臨床X線精確配準(zhǔn)。但采樣量需根據(jù)Volume數(shù)據(jù)大小而控制，Volume數(shù)據(jù)量越大，采樣范圍更大，改進(jìn)算法的提速優(yōu)勢更為明顯；如果超出理想采樣范圍，該改進(jìn)算法獲取的DRRs與臨床X線圖像配準(zhǔn)效果不理想，不宜應(yīng)用于臨床。

1 改進(jìn)的算法原理及方法

DRRs的生成是體繪制算法的一種典型應(yīng)用，它是模擬X線穿透CT體元，經(jīng)過衰減和吸收后投影到成像平面的過程，本文使用透視投影，如圖1所示。

DRRs生成過程中，如果入射光源的強(qiáng)度為I0，任何一光線穿過 “人體”后到達(dá)投影平面的強(qiáng)度為I，連續(xù)型積分表達(dá)式［11］為

式中：μ（x）——對應(yīng)組織的衰減系數(shù)，0——光線經(jīng)過體數(shù)據(jù)的穿入點(diǎn)，D——光線經(jīng)過體數(shù)據(jù)的穿出點(diǎn)。因為人體組織是非均勻的，不同的組織對應(yīng)不同的衰減系數(shù)，離散型表達(dá)式［11］為

下標(biāo)i為沿著射線路徑穿過的第i種組織的索引，xi是穿過第i種組織的有效長度，μi是對應(yīng)第i種組織的衰減系數(shù)，衰減系數(shù)μi與CT值的轉(zhuǎn)換公式參考如下［12］

圖1 光線投影算法

式中：μw——水的衰減系數(shù)，HUi——第i種組織的CT值，范圍為 ［－1000，3096］，單位為Hounsfield Unit。

光線投影算法生成DRRs的具體過程描述如下：光線穿過包含CT體數(shù)據(jù)的長方體，沿著光線的行進(jìn)方向以一定大小的步長采樣CT值，把得到的CT值轉(zhuǎn)換成相應(yīng)的衰減系數(shù)μi，累加此方向上的衰減系數(shù)，由此得到所有射線到達(dá)投影板上的強(qiáng)度I，最后形成DRRs。

上述DRRs生成過程中，每條穿透CT體元的模擬X線互不交叉，且得到投影板上像素值的計算形式相同，此特點(diǎn)正好與統(tǒng)一計算設(shè)備架構(gòu)CUDA的并行計算體系相吻合，因此，借助CUDA技術(shù)進(jìn)行高度并行計算，實現(xiàn)多條光線的同時繪制，大幅度加快了DRRs的形成。同時，穿透CT體元的任何模擬X線路徑都可表示為u→＋Δt v→，u→是模擬X線的起點(diǎn)，v→是模擬X線的方向，采樣開始于模擬X線進(jìn)入CT體元的起點(diǎn)tin，結(jié)束于光線穿出CT體元的終點(diǎn)tout。而降采樣是實時處理中一種常用方法。在體繪制中，它主要包括增大采樣步長Δt、增大起始點(diǎn)tin或減小終點(diǎn)tout、根據(jù)特征抽取特定信息和減少穿透CT體元的模擬X線的數(shù)量等方式。DRRs生成時，如果采用增大采樣步長Δt和只采樣骨頭組織相結(jié)合的降采樣，將大大減少DRRs過程的計算量，快速生成只含骨性結(jié)構(gòu)的DRRs。通過對頸部、胸部、腰部、腿部等多幅骨科臨床CT數(shù)據(jù)的骨窗和灰度直方圖分析，得知人體軟組織的灰度值小于860，骨組織的灰度值大于860。如果將整幅圖像的灰度值歸一化為0－1，采用灰度閾值采樣的方法，即采樣過程中只采樣灰度值大于0.2670的像素，同時增大采樣步長 （3倍甚至5倍原步長），結(jié)合CUDA硬件加速，快速獲取只含骨性結(jié)構(gòu)的DRRs。本文針對臨床骨科手術(shù)導(dǎo)航系統(tǒng)，該臨床應(yīng)用需獲取只含骨性結(jié)構(gòu)的DRRs來實現(xiàn)術(shù)中實時圖像引導(dǎo)。多次實驗證明，采用該改進(jìn)的算法——基于光線投影算法的DRRs生成原理，采用硬件CUDA加速，同時引入大步長與只采樣骨性結(jié)構(gòu)相結(jié)合的降采樣生成DRRs，可以大大減少重建DRRs的計算量，大幅度節(jié)約了生成DRRs的時間［13］，滿足了臨床骨科手術(shù)導(dǎo)航系統(tǒng)的實時性要求。

2 基于Fourier－Meilin變換和相位相關(guān)的圖像配準(zhǔn)

B.SrinivasaReddy和 B.N.Chatterji提出的 基于 Fourier－Meilin變換和相位相關(guān)相似性測度的醫(yī)學(xué)圖像配準(zhǔn)方法［14］具有運(yùn)算量小、抗干擾能力強(qiáng)、配準(zhǔn)精確度高和魯棒性強(qiáng)的等優(yōu)點(diǎn)，為驗證DRRs圖像質(zhì)量，本文采用此方法對獲取的DRRs和臨床的X線圖像進(jìn)行配準(zhǔn)。它將時域上的醫(yī)學(xué)圖像灰度經(jīng)過傅里葉變換到頻域上再求頻譜相位相關(guān)性——互功率譜特性，然后對目標(biāo)互功率譜進(jìn)行反傅里葉變換找到相應(yīng)的峰值位置，從而確定配準(zhǔn)參數(shù)。具體過程如下：

假設(shè)只存在平移變換時，f1（x，y）和f2（x，y）為兩個圖像信號，它們滿足關(guān)系 （即f2（x，y）是由f1（x，y）經(jīng)過平移得到）

根據(jù)傅里葉變換的性質(zhì)可得

F1（u，v）和F2（u，v）分別為f1（x，y）和f2（x，y）經(jīng)過傅里葉變換的結(jié)果，它們的互功率譜為

式中：F＊1（u，v）——F1（u，v）的復(fù)共軛，的傅里葉反變換為一個二維脈沖函數(shù)δ（x－x0，y－y0），相位相關(guān)法就是求取式（6）的反傅里葉變換，然后找到峰值位置最終確定配準(zhǔn)參數(shù)。

當(dāng)平移、旋轉(zhuǎn)、縮放參數(shù)同時存在時，設(shè)第一幅圖像為f1（x，y），第二幅圖像f2（x，y）是由第一幅圖像f1（x，y）平移、旋轉(zhuǎn)、縮放得到，可表示為

式中：σ——縮放尺度，θ0——旋轉(zhuǎn)角度，（x0，y0）——平移參量；配準(zhǔn)就是在參考圖像f1（x，y）中心處取一個圖像直角坐標(biāo)尺寸為N×N的小區(qū)域I1（x，y），配準(zhǔn)圖像f2（x，y）任何可能的位置截取一個與I1（x，y）大小相同的小區(qū)域I2（x，y），其次在對數(shù)極坐標(biāo)表示下，對I1（x，y）和I2（x，y）分別進(jìn)行傅里葉變換后求兩者的互功率譜，最后對互功率譜進(jìn)行傅里葉反變換，如果產(chǎn)生一個二維脈沖信號，則找出峰值位置估計出配準(zhǔn)參數(shù)；否則配準(zhǔn)失敗。

3 結(jié)果與分析

3.1 數(shù)字重建影像與采樣時間結(jié)果

本文按上述方法，在一臺CPU主頻為3.20GHZ、內(nèi)存為2G、顯卡為Nvidia Quadro 600、顯存為1G的PC機(jī)上處理了A、B、C三組16位CT數(shù)據(jù)，分辨率分別為250×512×512、831×512×512、1440×512×512，以下分別為A、B、C數(shù)據(jù)的DRRs效果圖和采樣時間的結(jié)果（λ＝0.01為原步長，B為只采樣骨性結(jié)構(gòu)）。

圖2（a）為只借助CUDA加速進(jìn)行原步長采樣獲取的A數(shù)據(jù)DRR結(jié)果，圖2（b）至圖2（f）為基于改進(jìn)的算法獲取的A數(shù)據(jù)DRRs結(jié)果。

圖2 不同采樣方式獲取的A數(shù)據(jù)DRRs

表1～表3為DRRs生成過程中不同采樣方式分別獲取一張A、B、C數(shù)據(jù)DRR的時間對照表。

表1 不同采樣方式獲取一張A數(shù)據(jù)DRR的時間對照

圖3（a）為只借助CUDA加速進(jìn)行原步長采樣獲取的B數(shù)據(jù)DRR結(jié)果，圖3（b）至圖3（f）為基于改進(jìn)的算法獲取的B數(shù)據(jù)DRRs結(jié)果。

圖3 不同采樣方式獲取的B數(shù)據(jù)DRRs

表2為B數(shù)據(jù)DRRs生成過程中不同采樣方式獲取一張DRR的時間對照表。

表2 不同采樣方式獲取一張B數(shù)據(jù)DRR的時間對照

圖4（a）為只借助CUDA加速進(jìn)行原步長采樣獲取的C數(shù)據(jù)DRR結(jié)果，圖4（b）至圖4（j）為基于改進(jìn)的算法獲取的C數(shù)據(jù)DRRs結(jié)果。

表3為C數(shù)據(jù)DRRs生成過程中不同采樣方式獲取一張DRR的時間對照表。

圖4 不同采樣方式獲取的C數(shù)據(jù)DRRs

表1結(jié)果顯示，3λ＋B采樣方式獲取A數(shù)據(jù)DRRs相比λ采樣方式獲取A數(shù)據(jù)DRRs的時間減少了53.60ms；表2結(jié)果顯示，3λ＋B采樣方式獲取B數(shù)據(jù)DRRs相比λ采樣方式獲取B數(shù)據(jù)DRRs的時間減少了67.90ms；表3結(jié)果顯示，5λ＋B采樣方式獲取C數(shù)據(jù)DRRs相比λ采樣方式獲取C數(shù)據(jù)DRRs的時間減少了100.60ms。由此得知，通過增大采樣步長與只采樣骨性結(jié)構(gòu)相結(jié)合的降采樣可以更大幅度提高DRRs的生成速度，且數(shù)據(jù)量越大，提速優(yōu)勢更明顯。

表3 不同采樣方式獲取一張C數(shù)據(jù)DRR的時間對照

3.2 不同采樣方式得到的DRRs與臨床X線圖像的配準(zhǔn)

改進(jìn)的算法可以大幅度提高DRRs的生成速度，同時保證了DRRs與臨床X線圖像的精確配準(zhǔn)。為驗證獲取的DRRs圖 像 質(zhì) 量，本 文 采 用 B.SrinivasaReddy和B.N.Chatterji提出的基于Fourier－Meilin變換和相位相關(guān)相似性測度的醫(yī)學(xué)圖像配準(zhǔn)方法對DRRs進(jìn)行配準(zhǔn)驗證，此處將A數(shù)據(jù)DRRs和A數(shù)據(jù)對應(yīng)的臨床X線、B數(shù)據(jù)DRRs和B數(shù)據(jù)對應(yīng)的臨床X線、C數(shù)據(jù)DRRs和C數(shù)據(jù)對應(yīng)的臨床X線配準(zhǔn) （用原步長采樣方式得到的模擬X線圖像代替臨床X線圖像，并設(shè)為參考圖像，降采樣得到的圖像為浮動圖像）。表4（a）顯示了A數(shù)據(jù)的DRRs和對應(yīng)臨床X線圖像配準(zhǔn)的結(jié)果，表4（b）顯示了B數(shù)據(jù)DRRs和對應(yīng)臨床X線圖像配準(zhǔn)的結(jié)果，表4（c）顯示了C數(shù)據(jù)的DRRs和對應(yīng)臨床X線圖像配準(zhǔn)的結(jié)果；圖5（a）、5圖（b）、圖5（c）分別顯示了采樣方式與配準(zhǔn)旋轉(zhuǎn)角度平均偏差、X方向平均偏差、Y方向平均偏差的關(guān)系。

表4 A、B、C數(shù)據(jù)DRRs分別與對應(yīng)臨床X線圖像的配準(zhǔn)結(jié)果

（c）C數(shù)據(jù)的DRRs和臨床X線圖像配準(zhǔn)結(jié)果

圖5 不同采樣方式與配準(zhǔn)偏差關(guān)系

表4（a）和圖5（a）～圖5（c）結(jié)果顯示，λ＋B、2λ＋B、3λ＋B方式獲取的A數(shù)據(jù)DRRs分別與A數(shù)據(jù)對應(yīng)的臨床X線圖像配準(zhǔn)的誤差和λ方式獲取的A數(shù)據(jù)DRRs與A數(shù)據(jù)對應(yīng)的臨床X線圖像配準(zhǔn)的誤差很相近（如果計算出的變換參數(shù)和真實值相差不超過1個像素或1度，即認(rèn)為配準(zhǔn)成功），以上誤差均在臨床可接受的范圍內(nèi)，分析可知增大采樣步長和只采樣骨頭相結(jié)合的降采樣不影響獲取的A數(shù)據(jù)DRRs和A數(shù)據(jù)對應(yīng)的臨床X線圖像的精確配準(zhǔn)，但實驗結(jié)果證明，如果繼續(xù)減少采樣量 （小于3λ＋B的采樣量），獲取的A數(shù)據(jù)DRRs和A數(shù)據(jù)對應(yīng)的臨床X線圖像的配準(zhǔn)將產(chǎn)生很大的偏差，很難保證實時圖像引導(dǎo)手術(shù)治療的準(zhǔn)確性，不宜采用；按上所述，表4（b）和圖5（a）～圖5（c）結(jié)果表明，增大采樣步長和只采樣骨頭相結(jié)合的降采樣不影響獲取的B數(shù)據(jù)DRRs和B數(shù)據(jù)對應(yīng)的臨床X線圖像的精確配準(zhǔn)，但實驗結(jié)果證明，如果繼續(xù)減少采樣量 （小于3λ＋B的采樣量），獲取的B數(shù)據(jù)DRRs和B數(shù)據(jù)對應(yīng)的臨床X線圖像的配準(zhǔn)將產(chǎn)生很大的偏差，很難保證實時圖像引導(dǎo)手術(shù)治療的準(zhǔn)確性，不宜采用；表4（c）和圖5（a）～圖5（c）結(jié)果表明，增大采樣步長和只采樣骨頭相結(jié)合的降采樣不影響獲取的C數(shù)據(jù)DRRs和C數(shù)據(jù)對應(yīng)的臨床X線圖像的精確配準(zhǔn)，但實驗結(jié)果證明，如果繼續(xù)減少采樣量 （小于5λ＋B的采樣量），獲取的C數(shù)據(jù)DRRs和C數(shù)據(jù)對應(yīng)的臨床X線圖像的配準(zhǔn)將產(chǎn)生很大的偏差，很難保證實時圖像引導(dǎo)手術(shù)治療的準(zhǔn)確性，不宜采用。

由上可知，當(dāng)采樣量控制在理想范圍內(nèi) （一般不低于3或5λ＋B的采樣量）時，改進(jìn)的算法可以更大幅度加快DRRs的形成，也能使得獲取的DRRs和臨床X線精確配準(zhǔn)；如果超出理想采樣范圍，獲取的DRRs不能保證和臨床X線的精確配準(zhǔn)。于此，采樣量的選擇需根據(jù)Volume數(shù)據(jù)量的大小而控制，需做到既能加快DRRs的形成，又不影響DRRs和臨床X線的精確配準(zhǔn)；并且Volume數(shù)據(jù)量越大，采樣范圍更大，改進(jìn)的算法的優(yōu)勢會更為明顯。

4 結(jié)束語

本文提出了一種改進(jìn)的數(shù)字重建影像生成算法——基于光線投影算法的DRRs生成原理，采用統(tǒng)一并行計算架構(gòu)CUDA加速，再引入大步長和只采樣骨性結(jié)構(gòu)相結(jié)合的降采樣，快速生成只含骨性結(jié)構(gòu)的DRRs，并將改進(jìn)算法獲取的DRRs和臨床X線圖像進(jìn)行了醫(yī)學(xué)圖像配準(zhǔn)驗證。實驗結(jié)果表明，當(dāng)采樣量控制在理想范圍內(nèi)時，本文采用改進(jìn)的算法獲取DRRs不但在提速上優(yōu)勢很明顯而且能保證獲取的DRRs與臨床X線圖像的精確配準(zhǔn)，從而實現(xiàn)實時圖像引導(dǎo)臨床骨科等手術(shù)，滿足了臨床骨科手術(shù)導(dǎo)航系統(tǒng)等臨床應(yīng)用的實時性和精確性的要求。

［1］Markelj P，TomaeviD，Likar B，et al.A review of 2D／3D registration methods for image－guided interventions ［J］.Medical Image Analysis，2012，16 （3）：642－661.

［2］Deutschmann H，Steininger P，Nairz O，et al. “Augmented reality”in conventional simulation by projection of 3－D structures into 2－D images ［J］.Strahlentherapie und Onkologie，2008，184 （2）：93－99.

［3］Jakob Spoerk，Helmar Bergmann，F(xiàn)elix Wanschitz，et al.Fast DRR splat rendering using common consumer graphics hardware ［J］.Med Phys，2007，34 （11）：4302－4308.

［4］Hadjira Bentoumi，Pascal Gautron，Kadi Bouatouch.GPU－based volume rendering for medical imagery ［J］.World Academy of Science，Engineering and Technology，2010，4 （1）：419－125.

［5］Gu Y，Wang C.Transgraph：Hierarchical exploration of transition relationships in time－varying volumetric data ［J］.IEEE Trans Vis Comput Graph，2011，17 （12）：2015－2024.

［6］Gao Y，Haapasalo M，Shen Y，et al.Development of virtual simulation platform for investigation of the radiographic features of periapical bone lesion ［J］.Journal of Endodontics，2010，36 （8）：1404－1409.

［7］Baumhauer M，F(xiàn)euerstein M，Meinzer H P，et al.Navigation in endoscopic soft tissue surgery：Perspectives and limitations［J］.Journal of Endourology，2008，22 （4）：751－766.

［8］Mori S，Kobayashi M，Kumagai M，et al.Development of a GPU－based multithreaded software application to calculate digitally reconstructed radiographs for radiotherapy ［J］.Radiological Physics and Technology，2009，2 （1）：40－45.

［9］Li X，Zhou L，Zhen X，et al.The generation of digitally reconstructed radiographs with six parameters ［C］／／4th International Conference on Bioinformatics and Biomedical Engineering.IEEE，2010：1－4.

［10］Fluck O，Vetter C，Wein W，et al.A survey of medical image registration on graphics hardware ［J］.Computer Methods and Programs in Biomedicine，2011，104 （3）：45－57.

［11］Ruijters D，ter Haar Romeny B M，Suetens P.GPU－accelerated digitally reconstructed radiographs ［C］／／Proceedings of the Sixth IASTED International Conference on Biomedical Engineering，2008：431－435.

［12］Jaffray D，Kupelian P，Djemil T，et al.Review of imageguided radiation therapy ［J］.Expert Review of Anticancer Therapy，2007，7 （1）：89－103.

［13］Osama M Dorgham，Stephen D Laycock，Mark H Fisher.GPU accelerated generation of digitally reconstructed radiographs for 2D／3Dimage registration ［J］.IEEE transactions on Biomedical Engineering，2012，59 （9）：2594－2603.

［14］Lin Y H，Chen C H.Template matching using the parametric template vector with translation，rotation and scale invariance［J］.Pattern Recognition，2008，41 （7）：2413－2421.