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

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

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

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

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

1 理論與方法

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

CS理論指出，如果信號是稀疏可壓縮的，那么就可以使用遠(yuǎn)低于Nyquist頻率的抽樣頻率對原始信號進行抽樣，并且能精確恢復(fù)出該信號?；诖?，稀疏投影角度下的CT迭代圖像重建，可以看成是求解如下不適定線性方程組問題［13］。

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

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

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

1.2.1 投影驅(qū)動模型

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

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

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

1.2.2 正／反投影矩陣

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

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

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

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

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

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

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

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

1.3.2 最優(yōu)化計算

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

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

2 仿真實驗與結(jié)果分析

為研究CSVD算法的應(yīng)用性能，在不同極稀疏投影數(shù)據(jù)下進行了2組數(shù)值仿真實驗：圓周掃描扇束投影的二維CT圖像重建、圓周掃描錐束投影的三維CT圖像重建。上述CT掃描的幾何布局如圖4所示，另外，綜合考慮參數(shù)設(shè)置對成像質(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ù)下的二維扇束圖像重建

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

為進一步檢驗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)對3種算法（FBP算法、ART-TV 算法和CSVD算法）的重建質(zhì)量進行客觀評價，結(jié)果如圖6所示。實驗結(jié)果表明：CSVD算法在各項指標(biāo)上明顯優(yōu)于其他2種算法，說明CSVD算法在重建圖像精確度、噪聲抑制和圖像恢復(fù)程度等方面有著更好的性能。

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

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

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

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

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

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

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

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

3 結(jié) 論

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

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

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

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

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