極稀疏投影數據的CT圖像重建

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

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

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

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

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

1 理論與方法

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

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

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

1.2 投影驅動的系統模型

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

1.2.1 投影驅動模型

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

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

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

1.2.2 正／反投影矩陣

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

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

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

1.3 壓縮感知的投影驅動CT圖像重建

在CS理論框架下，綜合傳統CT迭代圖像重建技術的優勢，并利用上述投影驅動系統模型，設計了一種CT迭代圖像重建算法，稱作CS的投影驅動CT圖像重建（Compressed Sensing View Driven CT image reconstruction，CSVD）算法，算法由基于投影驅動模型的粗略圖像重建和最優化計算2個環節組成。

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

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

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

1.3.2 最優化計算

圖3 循環2 Fig.3 Loop 2

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

2 仿真實驗與結果分析

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

圖4 圓周CT掃描系統結構示意圖Fig.4 Schematic diagram of circular CT scanning system

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

2.1 極稀疏投影數據下的二維扇束圖像重建

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

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

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

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

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

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

2.2 極稀疏投影數據下的三維錐束圖像重建

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

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

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

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

3 結 論

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

1）數值精度高。仿真實驗結果表明：極稀疏投影角度數下的重建圖像精確地再現了參考圖像的圖像結構及像素分布，另外CSVD算法的各項圖像質量衡量指標明顯優于ART-TV算法和FBP算法。

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

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

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