極稀疏投影數(shù)據(jù)的CT圖像重建

武麗君，孫豐榮，楊江飛，于倩蕾，賀芳芳

（1.山東大學(xué) 信息科學(xué)與工程學(xué)院，青島266237； 2.山東大學(xué) 微電子學(xué)院，濟(jì)南250101；3.山東大學(xué)附屬山東省立醫(yī)院 醫(yī)學(xué)影像科，濟(jì)南250014）

傳統(tǒng)計(jì)算機(jī)斷層成像（CT）機(jī)在螺旋掃描情況下旋轉(zhuǎn)一周一般需要投影1 000～2 000次［1-2］。已有研究顯示，CT檢查的X-射線輻射可增加病人患癌癥的風(fēng)險(xiǎn)［3-4］，過(guò)量的X-射線輻射還會(huì)對(duì)人體產(chǎn)生不可逆的輻射損害如染色體變異等［5-6］。臨床上，可以通過(guò)減少投影角度數(shù)來(lái)降低X-射線對(duì)人體的傷害。減少投影角度數(shù)也能夠縮短掃描時(shí)間從而改善動(dòng)態(tài)器官的成像效果。因此，研究如何在稀疏投影數(shù)據(jù)情況下進(jìn)行CT圖像重建具有重要的臨床意義。

針對(duì)不完全投影數(shù)據(jù)的CT圖像重建算法可以分為2類：第1類只適用于特定的掃描結(jié)構(gòu)，首先將不完整的投影數(shù)據(jù)進(jìn)行插值補(bǔ)充完整，然后再采用濾波反投影（Filtered Back-Projection，F(xiàn)BP）算法進(jìn)行圖像重建；第2類是迭代類的CT圖像重建算法，如代數(shù)重建技術(shù)（Algebraic Reconstruction Technique，ART），但傳統(tǒng)迭代類算法的計(jì)算時(shí)間要求高，在早期CT中沒(méi)有得到應(yīng)用。近年來(lái)，隨著大規(guī)模并行計(jì)算技術(shù)的發(fā)展以及計(jì)算機(jī)硬件成本的降低，迭代類算法已經(jīng)成為相關(guān)領(lǐng)域研究人員和CT機(jī)生產(chǎn)廠商高度關(guān)注的熱點(diǎn)。這其中較典型的是Sidky等在2006年提出的TV（Total Variation）算法［7-8］，該算法利用梯度圖像稀疏性先驗(yàn)知識(shí)提高了稀疏投影條件下的重建圖像質(zhì)量。TV算法具有很高的實(shí)用性，但是TV算法的圖像平滑效應(yīng)會(huì)降低圖像的空間分辨率和對(duì)比度。TV算法是一種啟發(fā)式的最優(yōu)化計(jì)算方法，壓縮感知（Compressed Sensing，CS）理論［9］在一定程度上為TV算法良好的應(yīng)用性能提供理論支撐。隨后又出現(xiàn)諸多基于CS理論的CT迭代圖像重建算法，較知名的是美國(guó)威斯康辛大學(xué)Chen博士及其研究組利用動(dòng)態(tài)CT圖像序列相關(guān)性先驗(yàn)知識(shí)提出的先驗(yàn)圖像約束壓縮感知（Prior Image Constrained Compressed Sensing，PICCS）算法 。PICCS算法要求能夠從冗余或完整的投影數(shù)據(jù)集合中選取均勻采樣的投影數(shù)據(jù)并利用解析的FBP算法獲得所謂的先驗(yàn)圖像，這極大地限制了該算法的應(yīng)用范圍；同時(shí)，當(dāng)先驗(yàn)圖像和待重建圖像不完全對(duì)應(yīng)時(shí)，由該算法得到的重建圖像會(huì)出現(xiàn)偽影。隨后，Debatin等［11-12］將各向異性全變差作為目標(biāo)函數(shù)進(jìn)行CT迭代圖像重建，提出了ATV（Anisotropic Total Variation minimization）算 法 和GATV（Generalized Anisotropic Total Variation minimization）算法，在稀疏投影數(shù)據(jù)的情況下，ATV算法和GATV算法較經(jīng)典的TV算法得到的重建圖像邊界更清晰、細(xì)節(jié)更豐富。

綜合已有文獻(xiàn)并與CT工程技術(shù)人員進(jìn)行交流，將［0，2π）范圍內(nèi)扇束／錐束掃描不超過(guò)20個(gè)投影角度稱為極稀疏投影，并著力于從極稀疏投影數(shù)據(jù)足夠精確地重建斷層圖像，從而能夠在顯著降低CT檢查X-射線輻射劑量的前提下，提供充分適宜影像學(xué)臨床診斷需求的重建圖像。為此，提出了投影驅(qū)動(dòng)的系統(tǒng)模型，并設(shè)計(jì)了一種將CT迭代圖像重建與CS理論相結(jié)合的重建算法；為檢驗(yàn)重建算法性能，分別對(duì)圓周掃描扇束和錐束投影得到的極稀疏投影數(shù)據(jù)進(jìn)行了圖像重建仿真實(shí)驗(yàn)；仿真實(shí)驗(yàn)結(jié)果表明，在極稀疏投影數(shù)據(jù)的條件下，算法數(shù)值精度高，計(jì)算復(fù)雜度低，內(nèi)存開(kāi)銷少，有很強(qiáng)的工程實(shí)用性。

1 理論與方法

1.1 基于壓縮感知的CT圖像重建

CS理論指出，如果信號(hào)是稀疏可壓縮的，那么就可以使用遠(yuǎn)低于Nyquist頻率的抽樣頻率對(duì)原始信號(hào)進(jìn)行抽樣，并且能精確恢復(fù)出該信號(hào)?；诖耍∈柰队敖嵌认碌腃T迭代圖像重建，可以看成是求解如下不適定線性方程組問(wèn)題［13］。

式中：Ψf為待重建圖像f的稀疏化表示。

1.2 投影驅(qū)動(dòng)的系統(tǒng)模型

在CT迭代圖像重建中，為了確定系統(tǒng)矩陣W 元素，必須建立合理的系統(tǒng)模型（即正／反投影運(yùn)算模型）。系統(tǒng)模型對(duì)CT迭代圖像重建的數(shù)值精度和重建圖像質(zhì)量都有著重要影響。目前，主流的系統(tǒng)模型主要分為3類：像素驅(qū)動(dòng)模型、射線驅(qū)動(dòng)模型和距離驅(qū)動(dòng)模型。像素驅(qū)動(dòng)模型和射線驅(qū)動(dòng)模型分別容易在投影域和圖像域引入高頻偽影；而距離驅(qū)動(dòng)模型則充分利用存儲(chǔ)空間的順序訪問(wèn)模式，算法復(fù)雜度相對(duì)較低，能夠避免圖像域偽影和投影域偽影［14-15］。結(jié)合像素驅(qū)動(dòng)和射線驅(qū)動(dòng)模型的優(yōu)點(diǎn)，并基于距離驅(qū)動(dòng)模型的基本思想，提出了如下投影驅(qū)動(dòng)的系統(tǒng)模型（即投影驅(qū)動(dòng)的正／反投影運(yùn)算模型）。

1.2.1 投影驅(qū)動(dòng)模型

圖1所示的二維扇束投影幾何描述了投影驅(qū)動(dòng)的正投影計(jì)算過(guò)程。在投影角度一定的情況下，將光源分別與某行像素邊界的中點(diǎn)與檢測(cè)器單元的邊界相連接，得到它們各自在公共軸（見(jiàn)圖1中的X軸）上的交點(diǎn)；然后在X軸上計(jì)算像素交點(diǎn)與相鄰檢測(cè)器交點(diǎn)的重疊長(zhǎng)度，并用歸一化的重疊長(zhǎng)度表示像素對(duì)檢測(cè)器正投影的貢獻(xiàn)；依次處理每一行像素。投影驅(qū)動(dòng)的反投影計(jì)算也是如此。至此，投影驅(qū)動(dòng)的正／反投影計(jì)算策略可概括為：一次迭代處理一個(gè)投影角度下的所有像素。

圖1也詳細(xì)地展示了檢測(cè)器邊界di、像素邊界pi、檢測(cè)器上投影值dij、像素值pij之間的關(guān)系。需要注意，投影驅(qū)動(dòng)的正／反投影中歸一化加權(quán)系數(shù)計(jì)算中的分母不同。正投影計(jì)算中，為了更精確地模擬線積分，將加權(quán)系數(shù)除以投影角度β的余弦值，投影驅(qū)動(dòng)的正投影運(yùn)算表示為

圖1 二維扇束投影示意圖與細(xì)節(jié)圖Fig.1 Schematic diagram and detailed diagram of two-dimensional fan-beam projection

1.2.2 正／反投影矩陣

在CT迭代圖像重建中，正／反投影矩陣將圖像域與投影域聯(lián)系起來(lái)。而且，上述正／反投影運(yùn)算模型（即系統(tǒng)模型）決定了正／反投影矩陣各元素的值。正投影運(yùn)算g＝Wf將圖像從圖像域映射到投影域，W 為正投影矩陣（即前文所述的系統(tǒng)矩陣）；f和g的含義前文已提及，分別為待重建圖像和投影數(shù)據(jù)列矢量。式（5）所示的反投影運(yùn)算將投影數(shù)據(jù)從投影域映射到圖像域，M 為反投影矩陣。

正／反投影矩陣W 和M 的結(jié)構(gòu)一致，但矩陣元素不同。圖2為正／反投影矩陣元素示意圖。結(jié)合圖2，進(jìn)一步闡述正／反投影矩陣各元素的含義。

正投影矩陣W 的結(jié)構(gòu)如圖2所示。圖中：ViewNum為投影角度數(shù)；R為像素行數(shù)；C為像素列數(shù)。Wθ為正投影矩陣W 在第θ個(gè)投影角度的子矩陣；Wθ（r）為在第θ個(gè)投影角度下第x行像素對(duì)所有檢測(cè)器單元的相對(duì)貢獻(xiàn)值，矩陣元素為式（3）中的加權(quán)系數(shù)。式（6）和式（7）清楚地揭示了以上3個(gè)矩陣的關(guān)系。 式（5）表明反投影矩陣M 的轉(zhuǎn)置MT直接參與反投影運(yùn)算，矩陣MT的結(jié)構(gòu)如圖2所示。與正投影類似，MTθ為矩陣MT在第θ個(gè)投影角度的子矩陣，MTθ（r）為在第θ個(gè)投影角度下，全部檢測(cè)器單元對(duì)第r行像素的相對(duì)貢獻(xiàn)值，矩陣元素為式（4）中的加權(quán)系數(shù)。以上3個(gè)矩陣的關(guān)系還可以由式（8）和式（9）表示。

1.3 壓縮感知的投影驅(qū)動(dòng)CT圖像重建

在CS理論框架下，綜合傳統(tǒng)CT迭代圖像重建技術(shù)的優(yōu)勢(shì)，并利用上述投影驅(qū)動(dòng)系統(tǒng)模型，設(shè)計(jì)了一種CT迭代圖像重建算法，稱作CS的投影驅(qū)動(dòng)CT圖像重建（Compressed Sensing View Driven CT image reconstruction，CSVD）算法，算法由基于投影驅(qū)動(dòng)模型的粗略圖像重建和最優(yōu)化計(jì)算2個(gè)環(huán)節(jié)組成。

圖2 正／反投影矩陣元素示意圖Fig.2 Schematic diagram of matrix elements of forward／backward projection

1.3.1 基于投影驅(qū)動(dòng)模型的粗略圖像重建

基于投影驅(qū)動(dòng)系統(tǒng)模型，迭代求解不適定的線性方程組g＝Wf，以獲得眾多滿足約束項(xiàng)的解。循環(huán)2還可以用圖3表示，每次循環(huán)依次處理每個(gè)投影角度下的投影數(shù)據(jù)；每個(gè)角度的投影數(shù)據(jù)都進(jìn)行正投影、作差、反投影和校正操作。圖中：Wθ和Mθ分別為上文提及的正／反投影矩陣的子矩陣；f（RR）［n，m，θ］為第n次總體迭代的第m次粗略重建子迭代中在第θ個(gè)投影角度下的圖像向量；Pθ為第θ個(gè)投影角度的實(shí)際投影數(shù)據(jù)；φθ為校正因子。

1.3.2 最優(yōu)化計(jì)算

圖3 循環(huán)2 Fig.3 Loop 2

式中：dA［n］為第n次總體迭代中計(jì)算得到的平衡因子；第n次總體迭代的第k次最優(yōu)化子迭代中的圖像向量表示為f（GRAD）［n，k］；α為步長(zhǎng)因子。

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

為研究CSVD算法的應(yīng)用性能，在不同極稀疏投影數(shù)據(jù)下進(jìn)行了2組數(shù)值仿真實(shí)驗(yàn)：圓周掃描扇束投影的二維CT圖像重建、圓周掃描錐束投影的三維CT圖像重建。上述CT掃描的幾何布局如圖4所示，另外，綜合考慮參數(shù)設(shè)置對(duì)成像質(zhì)量等多方面的影響，掃描參數(shù)的設(shè)置以及仿真模型的信息如表1所示。需要注意，圓周軌跡的錐束極稀疏投影掃描方式顯然不滿足Tuy條件［16］，故此算法屬于近似的錐束CT圖像重建算法。

圖4 圓周CT掃描系統(tǒng)結(jié)構(gòu)示意圖Fig.4 Schematic diagram of circular CT scanning system

表1 CT掃描參數(shù)和仿真模型信息Tab1e 1 CT scanning parameters and simu1ation mode1information

2.1 極稀疏投影數(shù)據(jù)下的二維扇束圖像重建

仿真實(shí)驗(yàn)先通過(guò)對(duì)圖5（a）所示的Shepp-Logan模型掃描獲得投影數(shù)據(jù)，再使用CSVD算法對(duì)該模型在不同投影角度數(shù)下進(jìn)行二維CT圖像重建，重建結(jié)果如圖5（b）所示。為了更精確地比較重建圖像與原始圖像的像素值，圖5（c）給出兩者在水平方向上圖像的中心像素值分布。重建圖像和像素值分布曲線顯示，36和20個(gè)投影角度下的重建圖像與參考圖像幾乎相同、像素曲線幾乎重合；18個(gè)投影角度下的重建圖像部分邊緣模糊且部分區(qū)域有塊狀偽影，但是各個(gè)組織結(jié)構(gòu)都能清晰分辨，像素值分布曲線在像素值跳變部位有沖激，其他部分與參考圖像曲線幾乎重合。

圖5 Shepp-Logan模型和重建圖像及其中心行像素值比較Fig.5 Shepp-Logan model and reconstructed images，as well as comparison of pixel values of center row between them

為進(jìn)一步檢驗(yàn)CSVD算法在極稀疏投影數(shù)據(jù)下的重建性能，使用標(biāo)準(zhǔn)均方根誤差（Normalized Root Mean Squared Error，NRMSE）、峰值信噪比（Peak Signal to Noise Ratio，PSNR）、全局質(zhì)量指標(biāo)（Universal image Quality Index，UQI）［17］3種衡量標(biāo)準(zhǔn)對(duì)3種算法（FBP算法、ART-TV 算法和CSVD算法）的重建質(zhì)量進(jìn)行客觀評(píng)價(jià)，結(jié)果如圖6所示。實(shí)驗(yàn)結(jié)果表明：CSVD算法在各項(xiàng)指標(biāo)上明顯優(yōu)于其他2種算法，說(shuō)明CSVD算法在重建圖像精確度、噪聲抑制和圖像恢復(fù)程度等方面有著更好的性能。

最后給出以上3種算法在20個(gè)投影角度下的3種衡量指標(biāo)具體數(shù)值以及重建時(shí)間，如表2所示。與FBP算法相比，CSVD算法耗費(fèi)了較多的時(shí)間，但是CSVD算法在不完全投影數(shù)據(jù)下的重建性能明顯優(yōu)于FBP算法。與ART-TV算法相比，CSVD算法的重建效率提高了2倍多，同時(shí)從以上3種衡量指標(biāo)可以看出CSVD算法的重建質(zhì)量有明顯提升。

圖6 各個(gè)投影角度數(shù)下不同算法的重建圖像質(zhì)量對(duì)比Fig.6 Comparison of reconstructed image quality of different algorithms under various projection angles

表2 各個(gè)投影角度數(shù)下不同算法的重建結(jié)果Tab1e 2 Reconstruction resu1ts of different a1gorithms under various projection ang1es

2.2 極稀疏投影數(shù)據(jù)下的三維錐束圖像重建

與二維扇束圖像重建實(shí)驗(yàn)類似，首先對(duì)自定義的三維Shepp-Logan模型掃描獲得投影數(shù)據(jù)，該模型透視圖如圖7（a）所示；其次，使用CSVD算法對(duì)該模型分別在20和18個(gè)投影角度數(shù)下進(jìn)行三維圖像重建，獲得中間6層（125層至130層）的橫斷面切片，選取其中的第125層和128層的重建結(jié)果作為展示，如圖7所示，其中圖7（b）、圖7（c）分別為20和18個(gè)投影角度下的重建圖像，同時(shí)還比較了第128層重建圖像與參考圖像在水平方向上圖像中心像素值，得到的像素分布曲線與圖5（c）類似，這里不再列出；最后，為了更客觀地評(píng)價(jià)重建圖像的質(zhì)量，使用NRMSE、PSNR、UQI三種衡量標(biāo)準(zhǔn)對(duì)第128層重建圖像質(zhì)量進(jìn)行分析如表3所示。

重建結(jié)果表明：在極稀疏投影數(shù)據(jù)的條件下，CSVD算法表現(xiàn)出良好的重建性能。在僅有18個(gè)投影角度的前提下，除重建圖像邊緣較模糊外，各個(gè)組織結(jié)構(gòu)仍能清晰分辨，像素值分布曲線重合度高，且在NRMSE、PSNR、UQI這3項(xiàng)指標(biāo)上與20個(gè)投影角度相當(dāng)。

圖7 Shepp-Logan模型透視圖和重建圖像（第125層和128層）Fig.7 Perspective of Shepp-Logan model and reconstructed images（125th and 128th layers）

表3 CSVD算法重建圖像質(zhì)量評(píng)價(jià)Tab1e 3 Qua1itv eva1uation of reconstructed images using CSVD a1gorithm

3 結(jié) 論

研究如何在極稀疏投影數(shù)據(jù)情況下進(jìn)行CT迭代圖像重建具有重要的臨床意義。理論分析與仿真實(shí)驗(yàn)結(jié)果都表明，CSVD算法能夠從極稀疏投影數(shù)據(jù)足夠精確地重建斷層圖像，從而將能夠在顯著降低CT檢查X-射線輻射劑量的前提下，提供充分適宜影像學(xué)臨床診斷需求的重建圖像。CSVD算法良好的圖像重建性能主要體現(xiàn)在：

1）數(shù)值精度高。仿真實(shí)驗(yàn)結(jié)果表明：極稀疏投影角度數(shù)下的重建圖像精確地再現(xiàn)了參考圖像的圖像結(jié)構(gòu)及像素分布，另外CSVD算法的各項(xiàng)圖像質(zhì)量衡量指標(biāo)明顯優(yōu)于ART-TV算法和FBP算法。

2）計(jì)算復(fù)雜度低。由前文的理論分析可知，CSVD算法一次迭代會(huì)處理一個(gè)投影角度，且在一個(gè)投影角度下，只需一次迭代就可以獲得一行像素對(duì)所有檢測(cè)器單元的貢獻(xiàn)（正投影過(guò)程），減少了大量不必要的遍歷運(yùn)算，使得計(jì)算復(fù)雜度大幅度降低。

3）內(nèi)存開(kāi)銷少。在醫(yī)學(xué)成像應(yīng)用中，系統(tǒng)矩陣的規(guī)模通常都很大，致使基于其他系統(tǒng)模型的CT迭代圖像重建一般都需要存儲(chǔ)系統(tǒng)矩陣，而基于投影驅(qū)動(dòng)系統(tǒng)模型的圖像重建則不需要對(duì)系統(tǒng)矩陣進(jìn)行存儲(chǔ)，大大減小了內(nèi)存的開(kāi)銷，CSVD算法有著很強(qiáng)的工程實(shí)用性。

總之，在極稀疏投影數(shù)據(jù)的條件下，CSVD算法數(shù)值精度高，計(jì)算復(fù)雜度低，內(nèi)存開(kāi)銷少，有很強(qiáng)的工程實(shí)用性。這給相關(guān)領(lǐng)域的科研工作者在極稀疏投影數(shù)據(jù)情況下進(jìn)行CT迭代圖像重建（從而可以降低CT檢查的X-射線輻射劑量）提供了一條切實(shí)的技術(shù)途徑。