基于離散剪切波正則化的低劑量CT圖像統計重建算法*

張海燕,張立毅,2，孫云山,2

(1.天津大學電子信息工程學院,天津 300072;2.天津商業大學信息工程學院,天津 300134)

1 引言

由于CT技術具有快速、準確、無創傷、無痛苦等特點，被越來越多地應用于臨床診斷。但是，X射線在透射過程中會將部分能量轉移到人體，引起身體損傷甚至致癌[1，2]。近年來，由于CT的廣泛使用，輻射風險也越來越受到關注，低劑量CT問題逐漸成為研究熱點。

當前臨床上為獲得高質量的CT重建圖像，所用CT設備X射線管電流強度都較高，一般臨床診斷X射線管電流強度一般超過200 mA[3]，一次CT掃描通常需要上千視角[4]。常見的降低CT劑量的方法主要有兩種，一是降低X射線管電流強度來降低每個視角下的曝光劑量[5-10]。降低管電流強度能降低單次掃描的輻射劑量，但會使投影數據信噪比下降，噪聲強度呈指數倍增長。二是在不改變單個視角曝光劑量的前提下，減少掃描視角的數量[11 - 16]。掃描視角的減少可以降低一次掃描的總輻射劑量，但會使投影數據的數目遠遠小于待建CT圖像的像素數目，使圖像重建問題變成一個欠定問題。

對于第一種降低劑量的策略，通常對投影數據進行降噪預處理[5 - 8]或使用統計迭代方法SIR(Statistical Iterative Reconstruction)[9,10]。對于第二種降低劑量的策略，可以通過引入待建圖像的先驗知識來作為正則化項，使高度病態的不完備投影數據獲得穩定而準確的重建。其中最常見的是全變分TV(Total Variation)正則化迭代算法[11 - 13]。但是，當投影數據含噪聲時，利用上述算法得到的重建圖像中出現了較為嚴重的噪聲干擾和條形偽影，且部分細節被噪聲掩蓋，不能得到滿意的重建圖像。并且，臨床診斷CT圖像重建時使用TV最小化這一先驗仍然存在一些問題：首先，TV最小化表達圖像特征的能力有限，對其進行稀疏性約束時容易導致在消除噪聲的情況下損失信號；其次，TV最小化約束建立在圖像分片光滑的基礎上，實際中的CT圖像并不能精確滿足這一條件。因此，基于TV最小化約束的不完備數據CT重建中，重建結果常常存在圖像邊緣不清晰、表達細小結構能力差、高噪聲下產生塊狀偽影等問題[10,16]。因此，尋找更合適的圖像稀疏變換和稀疏表達方式，是進一步提高不完備投影數據重建質量的關鍵。

圖像的稀疏度對圖像重建質量有著重要的影響[17]，CT圖像中常常包含了大量的曲線和局部信息，所以需要尋找更易于有效地捕捉局部圖像特征的變換域。剪切波(Shearlet)是由Guo等[18，19]采用具有合成膨脹的仿射系統構造的一種多維函數，繼承了曲波變換和輪廓波變換平移不變性和方向選擇性的優點，通過對基函數縮放、剪切和平移等仿射變換，生成具有不同特征的剪切波基函數，是多尺度幾何分析的最新發展，可以對圖像的細節信息給出最優的表示性能。目前在圖像修復、壓縮、融合等領域取得了很好的效果[20 - 24]。對于二維信號，剪切波變換不僅可以檢測到所有的奇異點，而且還可以自適應跟蹤奇異曲線的方向，且隨著尺度參數變化，可精確描述函數的奇異性特征，從而獲得更好的圖像稀疏表示能力。因此，基于離散剪切波變換的稀疏表示約束比基于離散梯度變換的TV最小化約束具有更好的欠完備投影CT圖像重建性能。

本文利用CT圖像信號在剪切波變換中的稀疏特性，提出基于離散剪切波變換的CT圖像統計迭代重建算法，將離散剪切波變換作為正則化項，實現低管電流掃描、少視角掃描以及兩種情況同時存在的情況下的低劑量CT圖像重建。

2 基本原理

2.1 剪切波變換原理

剪切波是合成小波的一個特例，剪切波變換通過對一個函數ψ進行伸縮、平移和旋轉生成一組類似小波的多尺度幾何分析工具，具有簡單的數學結構。對于二維圖像，剪切波可以檢測到所有的奇異點，而且可以自適應跟蹤奇異曲線的方向，具有很強的方向敏感性。因此，函數ψ在剪切波系統中起著至關重要的作用。

假設函數ψ進行伸縮、平移和旋轉生成一組基用矩陣ψ表示，那么圖像f由剪切波進行稀疏表示可描述為如下優化問題:

s.t.f=ψα

(1)

其中，α為稀疏系數向量，其中只有少數非零元素，‖‖0表示l0范數。

2.2 CT圖像重建原理

CT不能直接觀察到人體內部結構，只能得到探測器記錄下的不同視角的X射線通過人體后的投影數據，然后利用重建算法得到掃描斷層的線性衰減系數圖。重建算法具體可分為解析法和迭代法兩類。解析重建算法在投影數據不滿足數據完備性條件時，不能精確重建斷層圖像，所以一次CT掃描通常需要上千個掃描視角[4]。迭代重建運算時間高于解析重建算法，但對投影數據的完備性要求較低，在數據缺失和信噪比較低情況下有較好的表現，因此是低劑量CT圖像重建通常使用的方法。

Af=p

(2)

那么，CT圖像重建的過程可以描述為已知p和A，求解f的一個逆問題。由于噪聲的存在，一般轉化為最小化如下目標函數:

(3)

其中，‖‖2表示向量的l2-范數。

通過迭代最小化式(3)實現CT圖像重建。

3 基于剪切波變換的低劑量CT圖像重建

剪切波變換具有很強的二維圖像稀疏表示及去噪性能，本文提出一種基于剪切波變換的低劑量CT圖像重建算法。即在求解式(3)時加入待建圖像的剪切波稀疏表示作為正則化項，以縮小解空間，提高重建質量，變為求解如下問題:

s.t.f=ψαandAf=p

(4)

式(4)是一個約束最優化問題，它所對應的增廣拉格朗日方程為:

(5)

其中，μ和η為增廣拉格朗日系數，取固定值；向量λ和β為增廣拉格朗日系數向量，分別隨稀疏表示誤差以及重建圖像的投影誤差的減小而增大。

當X射線管電流強度下降，投影數據信噪比降低時，考慮投影數據的統計特性可以提高算法的抗噪性能。根據低管電流條件下真實投影數據的實驗分析結果，文獻[4-8]提出經系統校準及對數變換后的低劑量CT投影數據可近似認為是理想投影數據經加性噪聲n污染的結果，其中n近似服從空間非平穩高斯分布，其各個投影數據上的噪聲均值為0，方差與各投影數據自身的統計均值呈現非線性解析關系，其解析式可用如下公式描述:

(6)

Figure 1 Value of si圖1 si的值

由式(6)可知，噪聲強度隨著投影數值的增長呈指數增長，為了降低噪聲對重建效果的影響，在式(5)的目標函數中加入統計加權:

wi=exp(-pi)

(7)

其中，pi為第i個探測器(Detector)上獲得的投影數據。那么式(5)的目標函數變為:

(8)

其中，拉格朗日系數μ和η取固定值，而向量λ和β分別隨稀疏表示誤差以及重建圖像的投影誤差的減小而增大，βi為β中的第i個元素，1≤i≤M。λ和β分別如式(9)和式(10)所示:

λt=λt-1-μ(ft-1-ψαt-1)

(9)

(10)

其中，t≥1表示迭代次數。

式(8)存在兩個變量f和α，本文使用交替最小化的方法來求解這兩個變量。具體過程如下：

步驟1初始化。

初始化CT圖像為投影數據經濾波反投影FBP(Filtered Back Projection)算法重建的圖像f0，μ和η分別取固定值25和210。

步驟2對ft-1(t≥1)進行稀疏表示。假設，這時求剪切波系數的目標函數為:

(11)

由于求解l0范數問題是NP-hard問題，計算量太大，用l1范數代替[26]，目標函數變為:

(12)

步驟3固定剪切波系數αt，更新重建圖像ft。

這時的優化目標為:

(13)

步驟4驗證是否滿足終止條件。

4 實驗分析

通過仿真實驗來定量分析本文提出的重建算法的性能。實驗對人體肺部進行掃描，并在相同條件下分別使用FBP、聯立迭代重建技術SIRT(Simultaneous Iterative Reconstruction Technique)、TV正則化以及本文提出算法進行CT圖像重建，以對比重建效果。

4.1 少視角掃描實驗

為了驗證少視角投影時的重建效果，本實驗使用一幅充足劑量CT圖像作為數字模體進行計算機仿真。模體圖像共512×512像素，如圖2所示。

Figure 2 Numerical phantom圖2 仿真模體

仿真使用扇形投影，888個等角探測器單元，X線源掃描軌跡半徑為400 mm，探測器陣列與X-射線源呈同心圓配置。臨床上一次CT掃描通常需要上千視角[4]，為了驗證本文算法在少視角掃描時的有效性，本實驗分別在0°～360°內均勻選取了180和90個視角，并對每組投影數據分別用FBP、SIRT、TV正則化以及本文提出算法進行了圖像重建，重建圖像均為512×512像素。如圖3所示，分別給出了用FBP、SIRT、TV正則化以及本文提出的算法分別對180、90個視角的投影數據進行重建得到的CT圖像。

Figure 3 Comparison of images reconstructed from limited-angle projection data圖3 少視角投影重建效果對比

從重建效果圖可以看出，隨著掃描視角減少，FBP算法偽影最嚴重，效果最差；SIRT算法比FBP算法重建效果稍好，但是偽影仍然較為嚴重，圖像模糊；TV正則化算法雖然沒有明顯的條紋狀偽影，但是圖像過于平滑，細節丟失較多；而本文算法無論在清晰度、對比度還是細節保留度上，效果都是最好的。當掃描視角降到90個時，雖然由于掃描視角過少而在邊緣出現了一些偽影，但左右肺部細節結構仍然比其它算法清晰。為了更清晰地看到細節保留效果，本文對圖3中各子圖分別取右肺的一部分進行放大對比，如圖4所示。

Figure 4 Detail comparison of right lung in images reconstructed from limited-angle projection data圖4 少視角投影重建圖像右肺部細節對比

從圖4可以看出，本文所提出算法對于細小的血管和氣管的重建最為清晰，細節最為明顯，重建效果最好。FBP和SIRT算法偽影嚴重，TV正則化算法雖然沒有明顯偽影，但是過于平滑，細節丟失過多。

為了對所提算法進行定量的評價，本文對圖3中的重建圖像與理想模體圖像從均方根誤差RMSE(Root Mean Square Error)和圖像質量評價指標SSIM(Structural SIMilarity)[27]兩個指標來判斷重建圖像質量，分別見表1和表2。SSIM反映了兩幅圖像之間的結構相似性，SSIM值越靠近1，圖像相似度越高。RMSE定義為:

(14)

Table 1 RMSE comparison among the FBP,SART,TVregularization and the proposed algorithm from limited-angleprojection in human lung simulation experiments

RMSE越小，說明重建圖像與理想模體的差距越小，由表1可知，本文所提出的算法的RMSE要遠遠低于另外三種算法。雖然隨著投影視角的減少，RMSE有所增加，但性能仍然好于其他三種算法。

Table 2 SSIM comparison among the FBP,SIRT,TVregularization and the proposed algorithm from limited-angleprojection in human lung simulation experiments

由表2可以看到，本文所提出的算法得到的重建圖像與理想模體圖像的結構相似度最高。

由上述實驗結果可知，傳統的解析重建算法或迭代重建算法在投影數據不滿足數據完備性條件時，不能精確重建斷層圖像。本算法在掃描視角降到90個時，仍然能重構出與理想模體結構相似度較高的圖像，重建質量比傳統的迭代重建算法以及目前受到廣泛關注的TV最小化迭代算法都要好。本文提出的算法在大大降低了輻射劑量的同時保證了重建圖像質量。

4.2 低管電流掃描實驗

為了驗證低X-射線管電流強度投影時的重建效果，本實驗在投影數據中加入了文獻[5-8]提出的在10 mAs的管電流強度掃描得到的噪聲模型。

本實驗仍使用圖2所示的數字模體進行計算機仿真，模體圖像共512×512像素。仿真使用扇形投影，共888個等角探測器單元，旋轉中心與探測器的距離是400 mm，探測器陣列與X-射線源呈同心圓配置。本實驗在0°～360°內均勻選取了800個視角的投影數據，加入均值為0，方差如式(6)所示的高斯噪聲，并對每組投影數據分別用FBP、SIRT、TV正則化以及本文提出算法進行了圖像重建，重建圖像均為512×512像素。如圖5所示。

Figure 5 Effect of images reconstructed from low X-ray current intensity scanned projection data圖5 低X-射線強度管電流投影重建效果

由于噪聲方差與投影強度呈指數關系，圖5所示的FBP、SIRT以及TV正則化算法得到的重建結果均受噪聲影響較大，重建圖像模糊。而本文所提出的算法由于考慮了投影數據的統計特性，在數據保真項中加入了統計加權wi，并且正則化系數λ和β在重建過程中是自適應變化的，降低了噪聲對重建結果的影響，并且剪切波變換本身具有很強的去噪性能，因此本文所提出的算法重建出的圖像質量遠遠高于其他重建算法。在10 mAs電流強度800個掃描視角時重建圖像仍然清晰，和傳統方法所需要的200 mA的電流強度及上千掃描視角相比[3，4]，在保證重建質量的同時，輻射強度將為不足傳統方法的10%。

為了更清晰地看到細節保留效果，本文對圖5中各子圖分別取右肺的一部分進行放大對比，如圖6所示。

Figure 6 Detail comparison of right lung in images reconstructed from low X-ray current intensity scanned projection data圖6 低X-射線管電流投影重建圖像右肺部細節對比

從圖6可以看出，本文所提出算法對于細小的血管和氣管的重建最為清晰，細節最為明顯，重建效果最好。FBP和SIRT算法偽影嚴重、圖像模糊，幾乎無法顯示任何細節。TV算法雖然偽影較輕，但是過于平滑，細節丟失嚴重。本文所提出的算法在10 mAs電流強度800個掃描視角時，重建圖像細節清晰，具有較強抗噪性能。

本實驗所得的重建圖像和理想模體圖像的均方根誤差和結構相似度見表3。

Table 3 RMSE and SSIM comparison among the FBP,SART,TV regularization and the proposed algorithms from10 mA X-ray current intensity scanned projection data

由表3可知，本文所提出的算法在X-射線管電流強度降低、投影數據信噪比下降時重建的圖像的RMSE要遠遠低于另外三種算法，而結構相似度要遠遠高于其它三種算法。

5 結束語

本文針對低劑量CT投影數據欠完備或低信噪比的特點，根據投影數據所特有的空間非平穩噪聲特性，提出一種基于離散剪切波正則化的低劑量CT統計迭代重建算法。在數據保真項加入符合投影數據統計特性的加權系數，將待建圖像在剪切波域可以進行稀疏表示作為正則化項加入目標函數，并設置自適應的正則化參數以增加算法的魯棒性。本算法能夠取得比TV正則化方法更好的重建效果，更適于低劑量不完備投影數據的重建。對所提出重建方法做了充分的實驗驗證，數據表明本方法在輻射劑量降低到不足原來10%時仍能準確重建CT圖像。

[1] Luo Li-min, Hu Yi-ning, Chen Yang. Research status and prospect for low-dose CT imaging [J].Journal of Data Acquisition & Processing,2015,30(1):24-34.(in Chinese)

[2] Sodickson A, Baeyens P F, Andriole K P,et al.Recurrent CT,cumulative radiation exposure,and associated radiation-induced cancer risks from CT of adults[J].Radiology,2009,251(1):175.

[3] Hsieh J.Adaptive streak artifact reduction in computed tomography resulting from excessive x-ray photon noise[J].Medical Physics,1998,25(11):2139-2147.

[4] Hara A K,Paden R G,Silva A C,et al.Iterative reconstruction technique for reducing body radiation dose at CT:Feasibility study[J].American Journal of Roentgenology,2009,193(3):764-771.

[5] Zhang Yuan-ke, Zhang Jun-ying,Lu Hong-bing.Parameter estimation and statistical noise reduction for low-dose CT sinogram [J].Journal of Xidian University(Natural Science Edition),2011,38(3):99-106.(in Chinese)

[6] Wang J,Lu H,Wen J,et al.Multiscale penalized weighted least-squares sinogram restoration for low-dose X-ray computed tomography[J].IEEE Transactions on Biomedical Engineering,2008,55(3):1022-1031.

[7] Zhang H,Zhang L,Sun Y,et al.Projection domain denoising method based on dictionary learning for low-dose CT image reconstruction[J].Journal of X-ray Science and Technology,2015,23(5):567-578.

[8] Zhang Quan,Luo Li-min,Gui Zhi-guo.Improved nonlocal prior-based Bayesian sinogram smoothing algorithm for low-dose CT[J].Journal of Southeast University (Natural Science Edition),2014,44(3):499-503.(in Chinese)

[9] Wang Li-yan,Wei Zhi-hui.Linearized Bregman iterations for low-dose CT reconstruction [J].Journal of Electronics & Information Technology,2013,35(10):2418-2424.(in Chinese)

[10] Xu Q,Yu H,Mou X,et al.Low-dose X-ray CT reconstruction via dictionary learning[J].IEEE Transactions on Medical Imaging,2012,31(9):1682-1697.

[11] Candes E J,Romberg J,Tao T.Robust uncertainty principles:Exact signal reconstruction from highly incomplete frequency information [J].IEEE Transactions on Information Theory,2006,52(2):489-509.

[12] Qian Shan-shan,Huang Jing,Ma Jian-hua,et al.Nonmonotone total variation minimization based projection restoration for low-dose CT reconstruction [J].Acta Electronica Sinica,2011,39(7):1702-1707.(in Chinese)

[13] Mao Bao-lin, Chen Xiao-zhao,Xiao Da-yu,et al.Ordered subset image reconstruction studied by means of total variation minimization and fast first-order method in low dose computed tomograhpy [J].Acta Physica Sinica,2014,63(13):138701.(in Chinese)

[14] Yu H,Wang G.A soft-threshold filtering approach for reconstruction from a limited number of projections[J].Physics in Medicine and Biology,2010,55(13):3905-3916.

[15] Wang Lin-yuan, Liu Hong-kui,Li Lei,et al.Review of sparse optimization-based computed tomography image reconstruction from few-view projections [J].Acta Physica Sinica,2014,63(20):208702.(in Chinese)

[16] Rantala M,V?nsk? S,J?rvenp?? S,et al.Wavelet-based reconstruction for limited-angle X-ray tomography[J].IEEE Transactions on Medical Imaging,2006,25(2):210-217.

[17] Guo De-quan, Yang Hong-yu,Liu Dong-quan,et al.Overview on sparse image denoising [J].Application Research of Computers,2012,29(2):406-413.(in Chinese)

[18] Guo K,Labate D.Optimally sparse multidimensional representation using shearlets[J].SIAM Journal on Mathematical Analysis,2007,39(1):298-318.

[19] Easley G,Labate D,Lim W Q.Sparse directional image representations using the discrete shearlet transform[J].Applied and Computational Harmonic Analysis,2008,25(1):25-46.

[20] Lim W Q.The discrete shearlet transform:A new directional transform and compactly supported shearlet frames[J].IEEE Transactions on Image Processing,2010,19(5):1166-1180.

[21] Wang Lei, Li Bin,Tian Lian-fang.Medical image fusion based on shift-invariant shearlet transformation [J].Journal of South China University of Technology(Natural Science Edition),2011,39(12):145-152.(in Chinese)

[22] Zhao J,Lü L,Sun H.Multi-threshold image denoising based on shearlet transform[J].Applied Mechanics and Materials,2010,29-32:2251-2255.

[23] Guo K,Labate D.Characterization and analysis of edges using the continuous shearlet transform[J].SIAM Journal on Imaging Sciences,2009,2(3):959-986.

[24] Deng Cheng-zhi,Liu Juan-juan,Wang Sheng-qian,et al.Feature retained image inpainting based on sparsity regularization [J].Optics and Precision Engineering,2013,21(7):1906-1913.(in Chinese)

[25] Kalender W A. Computed tomography:Fundamentals,system technology,image quality,applications [M].Hoboken, NJ:John Wiley & Sons Inc.,2011.

[26] Chen S S,Donoho D L,Saunders M A.Atomic decomposition by basis pursuit[J].SIAM Review,2001,43(1):129-159.

[27] Wang Z,Bovik A C,Sheikh H R,et al.Image quality assessment:From error visibility to structural similarity[J].IEEE Transactions on Image Processing,2004,13(4):600-612.

附中文參考文獻：

[1] 羅立民,胡軼寧,陳陽.低劑量 CT 成像的研究現狀與展望[J].數據采集與處理,2015,30(1):24-34.

[5] 張元科,張軍英,盧虹冰.低劑量 CT 圖像模型參數估計及統計去噪研究[J].西安電子科技大學學報(自然科學版),2011,38(3):99-106.

[8] 張權,羅立民,桂志國.基于改進非局部先驗的 Bayesian 低劑量 CT 投影平滑算法[J].東南大學學報(自然科學版),2014,44(3):499-503.

[9] 王麗艷,韋志輝.低劑量 CT 的線性 Bregman 迭代重建算法[J].電子與信息學報,2013,35(10):2418-2424.

[12] 錢姍姍, 黃靜,馬建華,等.基于投影數據非單調性全變分恢復的低劑量 CT 重建[J].電子學報,2011,39(7):1702-1707.

[13] 毛寶林,陳曉朝,孝大宇,等.基于全變分最小化和快速一階方法的低劑量 CT 有序子集圖像重建[J].物理學報,2014,63(13):138701.

[15] 王林元,劉宏奎,李磊,等.基于稀疏優化的計算機斷層成像圖像不完全角度重建綜述[J].物理學報,2014,63(20):208702.

[17] 郭德全,楊紅雨,劉東權,等.基于稀疏性的圖像去噪綜述[J].計算機應用研究,2012,29(2):406-413.

[21] 王雷,李彬,田聯房.基于平移不變剪切波變換的醫學圖像融合[J].華南理工大學學報(自然科學版),2011,39(12):145-152.

[24] 鄧承志,劉娟娟,汪勝前,等.保留結構特征的稀疏性正則化圖像修復[J].光學精密工程,2013,21(7):1906-1913.