基于折射角的稀疏投影角度相襯CT圖像重建

司 凱,孫豐榮,宋尚玲,秦 峰

(1.山東大學信息科學與工程學院,濟南250100;2.山東大學第二醫(yī)院設(shè)備部,濟南250033)

基于折射角的稀疏投影角度相襯CT圖像重建

司 凱1,孫豐榮1,宋尚玲2,秦 峰1

(1.山東大學信息科學與工程學院,濟南250100;2.山東大學第二醫(yī)院設(shè)備部,濟南250033)

已有基于X射線吸收襯度機制的計算機斷層成像(CT)技術(shù)很難對由輕元素構(gòu)成的弱吸收物質(zhì)進行高質(zhì)量成像。X射線相位襯度CT成像是對弱吸收物質(zhì)具有超高分辨率的一種CT技術(shù),但該技術(shù)成像時間長、所需X射線輻射劑量大,不利于臨床推廣應(yīng)用,因此,研究稀疏投影角度條件下的X射線相位襯度CT圖像重建問題,基于壓縮感知圖像重建理論,使用折射角信息減少X射線輻射劑量,提出一種X射線相位襯度CT圖像重建算法。實驗結(jié)果表明,與濾波反投影算法相比,該算法在稀疏投影角度下可以得到較高質(zhì)量的重建圖像,在實際數(shù)據(jù)實驗中能獲得較高的峰值信噪比和數(shù)值準確性。

X射線相襯;計算機斷層成像;圖像重建;壓縮感知;折射角

1 概述

在1895年發(fā)現(xiàn)了X射線并拍攝了第一幅X光片后,隨著計算機技術(shù)的日漸成熟,X射線計算機斷層成像(Computerized Tomography,CT)技術(shù)應(yīng)運而生,并在醫(yī)學診斷、材料學、工業(yè)無損檢測和安保等領(lǐng)域得到了廣泛的應(yīng)用。傳統(tǒng)的X射線成像技術(shù)是基于物體各個部位對于硬X射線吸收量的不同來重建圖像,但是由輕元素構(gòu)成的弱吸收物質(zhì)的吸收量很小,成像效果不理想,因而限制了它的進一步應(yīng)

用。然而當X射線穿過物體時,它的強度和相移共同發(fā)生變化,在這個變化中本文使用δ代表X射線穿過物體發(fā)生的相移量,β則代表物體對X射線的吸收量。

研究發(fā)現(xiàn),對于人體肺部組織等弱吸收物質(zhì),X射線穿過時相移量δ的變化等級可達吸收量β的千倍以上［1］,使用硬X射線穿過物體的相移量來成像可以達到高密度分辨率,因此基于相移量的成像技術(shù)無疑更加適應(yīng)弱吸收物質(zhì)。在此背景下,文獻［2-3］分別于1995年和1996年在《Nature》上提出了硬X射線相位襯度(簡稱相襯)成像技術(shù),開創(chuàng)了基于相移量δ成像技術(shù)的研究領(lǐng)域。

近年來,X射線相襯成像技術(shù)發(fā)展出4種主要的方法:干涉儀成像法［4］、衍射增強成像(Diffraction Enhanced Imaging,DEI)［5］法、類同軸成像法［6］和光柵干涉儀成像法［7］。本文方法使用折射角信息并針對稀疏投影角度進行圖像重建,算法適應(yīng)于衍射增強成像法和光柵干涉儀成像法,可以改善傳統(tǒng)相襯成像技術(shù)的采集時間長、輻射劑量大等問題。

目前,關(guān)于X射線相襯CT圖像重建技術(shù)的研究主要有2個方向:第1個方向是重建出相位量的梯度▽δ［8-9］;第2個方向是重建出相位量δ［10-11］。文獻［12］分析將壓縮感知(Compressed Sensing,CS)理論與重建相位量梯度結(jié)合的可行性,提出利用經(jīng)典ART迭代方法求解該問題的算法,并利用實驗驗證了其算法的可靠性。文獻［13］在求解相位量梯度的問題中,也得到了較好的實驗結(jié)果。以上述研究為基礎(chǔ),本文使用折射角信息,基于恢復δ的重建策略,在CS圖像重建理論的框架下,采用基于距離驅(qū)動(Distance Driven,DD)的正/反投影運算,并引入有序子集(Ordered Subset,OS)等思想,提出一種基于CS的相襯CT迭代圖像重建算法。

2 CT圖像重建中的壓縮感知理論

在稀疏投影角度下的CT圖像重建屬于求解欠定線性系統(tǒng)的問題,即:

其中,A是m×n的系統(tǒng)矩陣,且m＜n;矩陣v∈Rn×1為待重建圖像;m維列向量u為投影數(shù)據(jù)。結(jié)合CS理論,如果待重建圖像v是S稀疏的,且測量矩陣A(CT圖像重建中的系統(tǒng)矩陣)滿足嚴格的等距同構(gòu)特性(Restricted Isometry Property,RIP)［14］,即對于1≤S≤n,可以定義一個等距同構(gòu)常數(shù)γS,滿足:

可以將式(1)轉(zhuǎn)化成下述有約束的l1范數(shù)最優(yōu)化問題［15-16］:

其中,Ψv為圖像v的稀疏化表示;Ψ為稀疏基矩陣。據(jù)此,當系統(tǒng)矩陣A滿足RIP時可以通過求解式(3)得到欠定線性系統(tǒng)式(1)理論上的精確結(jié)果。在實際應(yīng)用中,系統(tǒng)矩陣A很難滿足RIP性質(zhì),并且存在信號只能達到近似稀疏或投影數(shù)據(jù)含有噪聲等情況。雖然如此,鑒于求解有約束的l1范數(shù)最優(yōu)化問題的魯棒性,此領(lǐng)域的學者普遍將系統(tǒng)矩陣A近似的視作CS中的測量矩陣,并依然可以得到很好的重建效果［13,17］。

3 本文算法介紹

在本文算法中,為滿足式(3)中對目標函數(shù)的稀疏性要求,本文將x的梯度做為目標函數(shù)。在迭代重建過程中,本文采用了距離驅(qū)動的正/反投影運算,引入了自適應(yīng)性權(quán)重和OS思想并自適應(yīng)地選擇OS數(shù)量,使得算法的整體性能大大提高。在式(3)的最優(yōu)化環(huán)節(jié),可以采用任一種現(xiàn)有的尋優(yōu)算法,本文中使用模擬退火(Simulated Annealing, SA)算法。

3.1 折射角預(yù)處理

本文采用重建δ項的策略,如圖1所示,其中,l⊥代表射線傳播的垂直方向;p-q為旋轉(zhuǎn)坐標,將折射角定義為X射線穿過物體后發(fā)生偏折的角度θ,它與相位項的二維空間分布δ(x,y)的關(guān)系如下［2］:

圖1 折射角光路圖

在圖1的光路圖中,相位項δ滿足以下等式［18］:

(1)將采集到的折射角信息θ(p,β)與符號函數(shù)sgn(p)做卷積的預(yù)處理操作,得到投影數(shù)(即式(1)中的y)。

(2)采用本文重建算法對式(3)中有約束的l1范數(shù)最優(yōu)化問題進行求解。

3.2 啟發(fā)式求解有約束l1范數(shù)最優(yōu)化問題的方法

求解過程簡要偽代碼如下:

可以看出,迭代重建環(huán)節(jié)分為2個部分,粗略重建和最優(yōu)化,分別對應(yīng)本節(jié)啟發(fā)式策略中步驟(1)、步驟(2)的內(nèi)容。其中,粗略重建部分內(nèi)部包含IterationⅡ和IterationⅢ2種迭代。在最優(yōu)化部分,首先利用粗略重建后的圖像跟初始圖像計算偏差距離,然后使用SA算法得到最優(yōu)解。

3.2.1 有序子集思想和自適應(yīng)性權(quán)重

在IterationⅢ中,以往迭代類算法使用一個依據(jù)經(jīng)驗選擇的松弛因子來促進算法收斂,但是松弛因子的選定比較繁瑣,自適應(yīng)性差。本文利用迭代次數(shù)和當前重建圖像像素值提出一種自適應(yīng)性的權(quán)重w(3.2節(jié)求解過程簡要偽代碼第20行)替代松弛因子,可以增強本文算法對于不同投影數(shù)據(jù)的適應(yīng)能力。對于第m次內(nèi)部外層迭代:

其中,t(r,c,d)為距離驅(qū)動的正投影運算中第(r,c)個像素對第d個檢測器的貢獻值。

3.2.2 距離驅(qū)動的正/反投影運算

在IterationⅢ中,現(xiàn)有的正/反投影運算主要分為3類,即基于像素驅(qū)動、基于射線驅(qū)動和基于距離驅(qū)動的方法。基于像素驅(qū)動的方法硬件實現(xiàn)簡單,但是存在過高的運算復雜度以及在正投影運算中容易引入高頻噪聲等弊端。基于射線驅(qū)動的方法在反投影運算中會產(chǎn)生摩爾紋偽影,其運算方式會導致不連續(xù)的內(nèi)存訪問機制從而影響處理速度。而文獻［20］提出的基于距離驅(qū)動的方法具有低算術(shù)復雜度、高度連續(xù)的內(nèi)存訪問機制和正/反投影運算中不易產(chǎn)生偽影等優(yōu)點［21］。鑒于這3類運算方法的優(yōu)缺點,本文采用基于距離驅(qū)動的正/反投影運算完成本文方法中的迭代環(huán)節(jié),以期獲得更高的算法通用性和高質(zhì)量重建結(jié)果。

3.3 算法實現(xiàn)過程

基于上述分析以及旨在重建δ項的初衷,本文將折射角信息預(yù)處理后做為投影數(shù)據(jù),使用OS思想

從時間向度看，任何國家和地區(qū)的社會管理模式與思想都不是一成不變的，都有一個生成與轉(zhuǎn)換的過程。任何社會管理思想都深深依存于它所產(chǎn)生的社會環(huán)境，目的也是為了解決其當下社會面臨的秩序、發(fā)展、社會經(jīng)濟文化等問題。從追求政治和行政獨立，到以績效為核心，再到更為重視社會公正，發(fā)展任務(wù)的轉(zhuǎn)換決定了不同時期西方社會管理的特點。

結(jié)合距離驅(qū)動正/反投影運算的CT迭代圖像重建技術(shù)得到粗略重建結(jié)果,然后通過SA算法得出最優(yōu)解。綜上,算法可分為兩部分,折射角預(yù)處理、迭代重建,其中,迭代重建分為粗略迭代重建和最優(yōu)化計算。算法描述如下:

(1)折射角預(yù)處理

將折射角信息利用3.1節(jié)中敘述的方法與符號函數(shù)進行卷積操作得到真實投影數(shù)據(jù)矩陣P,其角度數(shù)為N。

(2)迭代重建

假定待重建圖像為矩陣V,算法總的迭代次數(shù)i=1,2,…,I。對于每一次總的迭代,由粗略迭代重建和最優(yōu)化計算兩部分組成:

1)粗略迭代重建

此過程內(nèi)部包含2種迭代:一種是粗略迭代重建部分自身的迭代,表示為IterationⅡ(迭代次數(shù)m=1, 2,…,M);一種是每一次IterationⅡ下所有投影角度的迭代,表示為IterationⅢ(迭代次數(shù)n=1,2,…,N)。

每一次IterationⅢ,即在每一個投影角度下,根據(jù)3.2.1節(jié)中方法自適應(yīng)地劃分K個有序子集。對于第n次IterationⅢ,每個有序子集(k=1,2,…,K)都進行正向投影、作差并加權(quán)、反向投影以及圖像修正4個過程,其中,算法初始值V(i=1,m=1,n=1,k=1)=0,具體步驟如下:

②作差并加權(quán):將相應(yīng)的真實投影數(shù)據(jù)P(i,k)與w×PDD_OS作差,得到差矩陣PDiff。

③反向投影:使用距離驅(qū)動的反向投影運算對PDiff進行反向投影得到偏差修正圖像矩陣VDiff。

④圖像修正:將本次IterationⅢ的初始圖像與偏差修正圖像作和,得到第n次IterationⅢ中第k個有序子集的重建圖像,即V(i,m,n+1,k)=VDiff+V(i,m,n,k)。

2)最優(yōu)化計算

最優(yōu)化計算環(huán)節(jié)需要根據(jù)上述粗略迭代重建過程得出計算的范圍,所以本環(huán)節(jié)可分為求解偏差距離和尋優(yōu)兩部分。具體如下:

①求解偏差距離:

令i=i+1,并將V(i)作為本次總體迭代的初始值,重復粗略迭代重建和最優(yōu)化計算過程,直到i=I,程序終止,得到最終重建結(jié)果V=V(i=I)。

4 實驗結(jié)果與分析

對于X射線相襯設(shè)備減小X射線輻射劑量最直接的方式就是稀疏投影角度掃描。本文分別使用仿真和真實數(shù)據(jù)的實驗結(jié)果來驗證稀疏投影角度下本文方法的性能,并將其與濾波反投影(Filtered Back-projection,FBP)算法結(jié)果進行比較。欠采樣掃描是一種不滿足奈奎斯特抽樣定理的投影數(shù)據(jù)采集方式,本文對于X射線源從0°～180°旋轉(zhuǎn),主要使用3種欠采樣投影數(shù)據(jù)采集方式:(1)180個角度模式,采集步長為1°;(2)90個角度模式,采集步長為2°; (3)36個角度模式,采集步長為5°。

本文使用以下3種通用圖像質(zhì)量評判標準來衡量本文方法的性能:

(1)歸一化均方根誤差(Normalized Root Mean Square Error,NRMSE):

其中,Xu,v和Yu,v分別為重建圖像和參考圖像在(u,v)時的像素值,圖像的大小為(U,V)。NRMSE的值越小,說明重建圖像有越高的數(shù)值精確性。

(2)峰值信噪比(Peak Signal to Noise Ratio, PSNR):

其中,Xmax為重建圖像中像素值的最大值。PNSR值越大,說明重建圖像失真越小。

(3)統(tǒng)一質(zhì)量指標(Universal Quality Index, UQI),用于衡量2個圖像接近程度的指標,其值越高,重建圖像越接近參考圖像,最大值為1:

其中,μx,σx和μy,σy分別代表重建圖像和參考圖像的像素值均值和標準差;Cov(X,Y)代表兩圖像間的協(xié)方差。

4.1 仿真實驗分析

仿真實驗使用經(jīng)典Shepp-Logan人體頭部模型來模擬和重建,并修改其模型參數(shù)為相應(yīng)的相位參數(shù)來仿真純相位頭部模型［22］。

圖2分別使用本文算法和FBP算法重建圖像,其中,從左向右分別為180個、90個和36個投影角度下的重建圖像;重建圖像大小均為400×400像素;實驗中迭代次數(shù)I,M分別為6和4。從圖2可以看出,與FBP算法重建結(jié)果相比較,隨著投影角度的減少,本文算法重建的圖像質(zhì)量下降不明顯,依然可以保持較高的圖像質(zhì)量,且不產(chǎn)生顯著的偽影。

圖2 仿真實驗的重建結(jié)果

將依次分析本文算法和FBP算法的定量指標并做相應(yīng)對比。分別將相應(yīng)算法在360個角度下的重建結(jié)果(滿足數(shù)據(jù)完備性條件)作為相應(yīng)算法的參考圖像Y,結(jié)果如圖3所示。

圖3 仿真實驗結(jié)果分析

由圖3可以看出,隨著投影角度的減少,本文算法在重建圖像的信噪比、數(shù)值精確性以及與參考圖像的相似度等方面比FBP算法有著明顯的優(yōu)勢。在此基礎(chǔ)之上,在稀疏投影角度情況下,本文算法可以較高精度地恢復出原始圖像。

4.2 實際數(shù)據(jù)實驗分析

實際實驗投影數(shù)據(jù)來自日本高能加速器研究組織(High Energy Accelerator Research Organization)。同步輻射裝置中X射線的能量是35 kV,檢測器有450個。實驗樣品是一個有機玻璃體中包含13個圓柱狀的通道,通道中分別填充水、甘油、乙醇。裝置從上到下分別掃描500層,每層以0.5°的間隔掃描360次,共180°,每一層共采集360組數(shù)據(jù)。

本文中仿真實驗使用本文算法和FBP算法分別

對最具代表性的第0層和第300層投影數(shù)據(jù)進行重建,采用等距間隔(分別為1°,2°,5°)獲得180個、90個、36個角度進行重建。圖4為實際數(shù)據(jù)實驗重建圖像,其中,從左至右依次為180個、90個和36個角度下的重建圖像。從圖4可以看出,在投影角度減少的情況下,本文算法在偽影控制方面比FBP算法擁有更好的效果,且可以得到不輸FBP算法的圖像重建效果。在衡量實際數(shù)據(jù)實驗重建圖像質(zhì)量指標中,方便起見,本文只對第300層的結(jié)果進行分析并與FBP指標對比,結(jié)果如圖5所示。

圖4 實際數(shù)據(jù)實驗重建圖像

圖5 實際數(shù)據(jù)實驗結(jié)果分析

從圖5可以看出,在實際數(shù)據(jù)實驗中,較之FBP的分析結(jié)果,本文算法同樣表現(xiàn)出更高的信噪比以及數(shù)值準確性,此外,與參考圖像相似度方面也略勝一籌。證明本文算法在處理實際數(shù)據(jù)的圖像重建問題中也有著不遜于FBP算法的性能。

5 結(jié)束語

基于壓縮感知理論和恢復相位的策略,本文提出一種相襯CT迭代圖像重建算法。該算法將折射角信息經(jīng)預(yù)處理后,迭代重建出最終結(jié)果。在迭代重建部分,采用有序子集思想和自適應(yīng)性權(quán)重,以增加對噪聲的平滑作用和對不同投影數(shù)據(jù)的適應(yīng)性。實驗結(jié)果表明,對于欠采樣的仿真數(shù)據(jù)和實際數(shù)據(jù),本文算法有較好的圖像重建質(zhì)量和性能。由于目標函數(shù)帶來的重建圖像邊緣銳度欠佳,因此后期需要改進目標函數(shù)構(gòu)造來進一步優(yōu)化圖像質(zhì)量。

［1]陳建文,高鴻奕,李儒新,等.X射線相襯成像［J］.物理學進展,2005,25(2):175-194.

［2]Davis T J,Gao D,Gureyev T E,et al.Phase-contrast Imaging of Weakly Absorbing Materials Using Hard X-rays［J］.Nature,1995,373(6515):595-598.

［3]Wilkins S W,Gureyev T E,Gao D,et al.Phase-contrast Imaging Using Polychromatic Hard X-rays［J］.Nature, 1996,384(6607):335-338.

［4]Momose A.DemonstrationofPhase-contrastX-ray Computed Tomography Using an X-ray Interferometer［J］.Nuclear Instruments and Methods in Physics Research SectionA:Accelerators,Spectrometers,Detectorsand Associated Equipment,1995,352(3):622-628.

［5]Chapman D,Thomlinson W,Johnston R E,etal.Diffraction Enhanced X-ray Imaging［J］.Physics in Medicine and Biology,1997,42(11).

［6]Snigirev A,SnigirevaI,KohnV,etal.Onthe Possibilities of X-ray Phase Contrast Microimaging by CoherentHigh-energySynchrotronRadiation［J］.Review ofScientificInstruments,1995,66(12): 5486-5492.

［7]Pfeiffer F,Bech M,Bunk O,et al.Hard-X-ray Dark-field Imaging Using a Grating Interferometer［J］.Nature Materials,2008,7(2):134-137.

［8]Wang Zhentian,Zhang Li,Huang Zhifeng,et al.An ART Iterative Reconstruction Algorithm for Computed Tomography of Diffraction Enhanced Imaging［J］.Chinese Physics C,2009,33(11):975-980.

［9]張 凱,朱佩平,黃萬霞,等.代數(shù)迭代重建算法在折射襯度CT中的應(yīng)用［J］.物理學報,2008,57(6): 3410-3418.

［10]李 鏡,劉文杰,朱佩平,等.基于光柵相襯成像的扇束螺旋CT重建算法［J］.光學學報,2010,30(2): 421-427.

［11]Zhu Peiping,Zhang Kai,Wang Zhili,et al.Low-dose, Simple,and Fast Grating-based X-ray Phase-contrast Imaging［J］.Proceedings of the National Academy of Sciences,2010,107(31):13576-13581.

［12]李 鏡,孫 怡.基于L1范數(shù)的微分相位襯度CT稀疏角度重建算法［J］.光學學報,2012,32(3):70-76.

［13]秦 峰,孫豐榮,宋尚玲,等.基于壓縮感知的微分相襯CT迭代圖像重建［J］.計算機應(yīng)用,2013,33(6): 1732-1736.

［14]Candes E J,TerenceT.DecodingbyLinearProgramming［J］.IEEE Transactions on Information Theory, 2005,51(12):4203-4215.

［15]Candes E J,Romberg J K,Terence T.Stable Signal Recovery from Incomplete and Inaccurate Measurements［J］.Communications on Pure and Applied Mathematics,2006,59(8): 1207-1223.

［16]Candès E J,Romberg J,Terence T.Robust Uncertainty Principles:Exact Signal Reconstruction from Highly Incomplete Frequency Information［J］.IEEE Transactions on Information Theory,2006,52(2):489-509.

［17]Candès E J,Wakin M B.An Introduction to CompressiveSampling［J］.IEEESignalProcessing Magazine,2008,25(2):21-30.

［18]Sunaguchi N,Yuasa T,Huo Qingkai,et al.Convolution Reconstruction Algorithm for Refraction-contrast Computed Tomography Using a Laue-case Analyzer for Dark-field Imaging［J］.Optics letters,2011,36(3): 391-393.

［19]Hudson H M,Larkin R S.Accelerated Image Reconstruction Using Ordered Subsets of Projection Data［J］.IEEE Transactions on Medical Imaging,1994,13(4):601-609.

［20]De M B,Basu S.Distance-driven Projection and Backprojection［C］//Proceedings of IEEE Nuclear Science Symposium Conference.［S.l.］:IEEE Press,2002: 1477-1480.

［21]De M B,BasuS.Distance-drivenProjectionand Backprojection in Three Dimensions［J］.Physics in Medicine and Biology,2004,49(11):2463-2475.

［22]李濤濤,李 華,刁麓弘.復雜物體的X射線相襯成像三維仿真［J］.系統(tǒng)仿真學報,2009,21(1):32-35.

編輯 劉 冰

Sparse Projection Angle Phase Contrast CT Image Reconstruction Based on Refractive Angle

SI Kai1,SUN Fengrong1,SONG Shangling2,QIN Feng1
(1.School of Information Science and Engineering,Shandong University,Jinan 250100,China;
2.Department of Equipment,The Second Hospital of Shandong University,Jinan 250033,China)

It is not easy for conventional X-ray absorption contrast based Computed Tomography(CT)imaging technique to get a high quality image from the weak absorption materials comprising of light elements.X-ray phase contrast CT imaging is a CT imaging technique aiming at weak absorption materials and possessing ultrahigh resolution.However,X-ray phase contrast CT imaging takes too long time and requires too high X-ray radiation dose,which can set back the promotion of clinical applications.So it has significant meaning to study the image reconstruction for phase contrast CT using few projection views.This paper is under the theory of Compressed Sensing(CS),uses refractive angle information,aims to reduce the X-ray radiation dose,proposes a X-ray phase contrast CT image reconstruction algorithm.Experimental results show that compared with Filtered Back-projection(FBP)algorithm,this algorithm can use few projection views to obtain reconstructed images of higher quality,and can get higher Peak Signal to Noise Ratio(PSNR) and numerical accuracy in the actual data experiment.

X-ray phase contrast;Computed Tomography(CT);image reconstruction;Compressed Sensing(CS); refractive angle

司 凱,孫豐榮,宋尚玲,等.基于折射角的稀疏投影角度相襯CT圖像重建［J］.計算機工程,2015, 41(3):262-268.

英文引用格式:Si Kai,Sun Fengrong,Song Shangling,et al.Sparse Projection Angle Phase Contrast CT Image Reconstruction Based on Refractive Angle［J］.Computer Engineering,2015,41(3):262-268.

1000-3428(2015)03-0262-07

:A

:TP391.41

10.3969/j.issn.1000-3428.2015.03.049

國家自然科學基金資助項目(61071053);山東省自然科學基金資助項目(ZR2014FM006)。

司 凱(1989-),男,碩士研究生,主研方向:醫(yī)學圖像處理;孫豐榮,教授;宋尚玲,博士;秦 峰,碩士研究生。

2014-03-03

:2014-05-18E-mail:sunfr@sdu.edu.cn