基于殘差圖分解與平滑塊排列的低劑量CT投影降噪算法

王燕玲， 薛 敏， 邵燕靈， 劉 祎， 張 權， 陳 燕， 桂志國

(1. 中北大學 生物醫(yī)學成像與影像大數(shù)據(jù)山西省重點實驗室， 山西 太原 030051； 2. 中北大學 信息與通信工程學院， 山西 太原 030051； 3. 中北大學 理學院， 山西 太原 030051)

0 引 言

X射線計算機斷層成像(X-ray Computed Tomography, CT)由于其成像分辨率高的優(yōu)點已成為放射診斷領域內(nèi)不可缺少的一部分. 然而CT掃描中大量的輻射劑量會對人體造成難以避免的傷害[1]， 所以低劑量CT成像技術受到了越來越多的關注. 然而降低放射劑量往往會造成圖像質量的明顯退化， 因此如何有效抑制低劑量CT圖像中的偽影和噪聲成為目前研究中的一個熱點問題. 目前的處理方法主要集中在3個方面： ① 直接對濾波反投影(Filtered-Back Projection, FBP)算法重建的低劑量CT圖像進行后處理[2-3]， 降低圖像噪聲和偽影； ② 統(tǒng)計迭代重建算法[4-5]； ③ 先對低劑量CT投影數(shù)據(jù)濾波， 然后用FBP算法進行重建[6-9]. 第3類算法由于其具有FBP算法重建速度快、 應用技術成熟等優(yōu)勢而受到廣泛關注. Wang等[6]提出了懲罰重加權最小二乘(Penalized Reweighted Least Squares, PRWLS)算法， 并將其應用于低劑量CT圖像重建中， 較好地保持了圖像的邊緣. 劉祎等[7]提出了一種基于局部模糊熵的各項異性擴散(Local Fuzzy Entropy Anisotropic Diffusion Filter, LFEADF)投影降噪算法， 得到了高信噪比的重建圖像. 張權等[8]將結合非局部先驗模型的Bayesian降噪算法應用于低劑量CT圖像投影域降噪， 有效地保持了圖像的邊緣.

平滑塊排列(Smooth Patch Ordering, SPO)[10]算法是基于局部塊處理方法的新發(fā)展， 被應用于圖像去噪和修復中以提高圖像質量. SPO算法通過選取最短路徑的方式重排圖像塊， 然后應用一維濾波進一步處理圖像， 具有良好的去噪能力. 然而， 在濾除噪聲的同時， 不可避免地會丟失圖像中有用的結構成分， 而這部分信息往往在殘差圖像中可以被發(fā)現(xiàn).

如何更好地利用殘差圖像得到了越來越多學者的關注. 文獻[11]通過三維協(xié)同濾波和剩余塊的結構提取， 更好地保持了圖像結構特征. 文獻[12] 通過殘差圖像后處理的方式進一步提高了K-奇異值分解(K-Singular Value Decomposition, K-SVD)字典學習方法的去噪性能， 并發(fā)展了一種通用的遞歸算法來增強算法的圖像去噪能力[13]. 孫偉峰等[14]提出了一種多級殘差圖像濾波的非局部均值圖像去噪新方法， 有效提升了傳統(tǒng)非局部均值算法的性能.

低劑量CT圖像由于受到噪聲和偽影的影響， 圖像質量嚴重退化， 本文將具有強去噪性能的SPO算法應用于低劑量CT投影降噪中. 另外， 為進一步提升SPO算法性能， 充分利用殘差圖中的有用信息， 采用改進的形態(tài)分量分析(Morphological Component Analysis, MCA)[15]方法將SPO算法濾波后的殘差圖進行分解. 通過將殘差圖中提取的結構信息反饋到濾波投影， 使得獲得的低劑量CT重建圖像噪聲得到了有效抑制， 圖像細節(jié)、 邊緣以及結構等特征也得到了很好的保持. 基于不同體模的實驗結果驗證了本文所提基于殘差圖分解與平滑塊排列的投影恢復算法的有效性.

1 預備知識

1.1 噪聲模型

低劑量CT投影數(shù)據(jù)中含有明顯的噪聲， 通常使用的模型有高斯分布模型和泊松分布模型， 本文采用文獻[16]的研究結果， 低劑量CT投影數(shù)據(jù)經(jīng)系統(tǒng)校準和對數(shù)變換后近似服從非平穩(wěn)高斯分布， 均值和方差之間滿足如下非線性關系

(1)

1.2 平滑塊排列算法

文獻[10]提出的基于平滑塊排列的SPO算法中， 令A是一個N1×N2的圖像， 其中N1N2=N， 令B是A經(jīng)過損壞后獲得的圖像， 損壞后的圖像滿足如下表達式

b=Ma+v,(2)

式中：M為逐點損壞圖像的線性算子；b和a分別為圖像B和A的列向量形式堆疊表示；v為加性高斯白噪聲.

(3)

1.3 形態(tài)分量分析方法原理

文獻[15]提出的MCA方法是一種新的信號分解方法. 假設圖像F是K個不同形態(tài)分量的線性組合， 其表達式為

(4)

式中：Fi為第i個形態(tài)分量. 假設當選取的基或過完備字典適當時， 每個形態(tài)分量都能在其下稀疏表示， 且每個形態(tài)分量的基或者字典又是互不相干的. 于是， 圖像分解問題可以轉化為能量最小化問題[17]， 其表達式為

MCA方法通常采用迭代方法來求解式(5)， 其中字典的選擇起著至關重要的作用.

1.4 殘差圖

殘差圖就是含噪圖與去噪圖之間的差值. 在理想情況下， 通過去噪算法濾波得到的殘差圖應該近似白噪聲. 而在實際應用中， 即使使用最先進的去噪算法， 其殘差圖也會包含部分細節(jié)信息[18]. 所以有效提取殘差圖中的有用信息， 對于增強去噪效果具有重要意義. 本文就是通過從殘差圖中提取結構信息而進一步增強SPO算法保持圖像結構的能力.

2 基于殘差圖分解與平滑塊排列的投影恢復方法

2.1 算法原理

在傳統(tǒng)的MCA方法中， 對于圖像的不同形態(tài)分量使用固定字典對其稀疏表示， 而固定字典一般都是根據(jù)經(jīng)驗來選擇的. 選擇合適的字典對于MCA分解起著十分關鍵的作用， 選擇不當?shù)淖值淇赡軐е滤惴ㄊ? 本文采用改進的MCA算法分解殘差圖， 具體方法描述如下：

第一， 從殘差圖Y中逐一像素提取圖像塊Ym，m=1,2,…,M，M為提取出來的圖像塊總數(shù).

第二， 由于在線字典學習算法通過塊坐標下降方法來更新字典， 可以快速得到與源圖像具有高度匹配性的字典. 因此為了更有效地提取低劑量CT投影殘差圖特征， 本文采用在線字典學習算法[19]從殘差圖中訓練得到字典D.

第三， 應用類似全變分的方法[20]測量字典中的每一個原子， 將字典D分為結構子字典Ds和噪聲子字典Dn.

對于字典D中任一個原子a， 其全變分測量表達式如下

式中：n為每個圖像塊的大小. 式(6)計算所得的值最大為1， 這些值可以反映每個原子的平滑度. 通過確定一個閾值L來區(qū)分噪聲原子和結構原子， 小于等于L的原子為結構原子， 大于L的原子為噪聲原子.

第四， 采用正交匹配(Orthogonal Matching Pursuit, OMP)算法[21]對每一個圖像塊Ym基于字典D進行稀疏編碼， 由此得到的稀疏表示系數(shù)記為α， 根據(jù)上述結構原子和噪聲原子的劃分可以將稀疏表示系數(shù)分別對應表示為αs和αn.

第五， 由結構字典和噪聲字典及其對應的稀疏表示系數(shù)分別重建結構部分和噪聲部分投影.

2.2 算法描述

3) 從殘差圖Y中逐一像素提取圖像塊Ym， 將在線字典學習方法應用于殘差圖訓練， 得到字典D， 使得這些字典能夠稀疏表示所有圖像塊；

4) 確定一個閾值L， 用類似全變分的方法將字典D分成兩部分： 結構子字典Ds和噪聲子字典Dn；

5) 采用OMP算法對每一個圖像塊Ym基于字典D進行稀疏編碼， 得到分別對應于結構子字典Ds和噪聲子字典Dn的稀疏表示系數(shù)αs和αn；

6) 由字典Ds和系數(shù)αs重建結構部分Ys， 由字典Dn和系數(shù)αn重建噪聲部分Yn；

3 實驗結果與討論

為了驗證算法的有效性， 本文采用改進的Shepp-Logan頭部體模和數(shù)字骨盆體模進行仿真實驗， 其大小均為256 mm×256 mm. 按照式(1)在理想投影數(shù)據(jù)中加入噪聲獲得的投影數(shù)據(jù)作為仿真的低劑量投影數(shù)據(jù)， 其中T=22 000， 頭部體模中fi=200， 骨盆體模中fi=100. 本文算法與PRWLS算法[6]、 LFEADF算法[7]和SPO算法[10]進行了比較. 仿真實驗采用FBP算法對各種方法濾波后的投影重建獲得低劑量CT重建圖像. 本文算法采用Matlab2012b編寫， 實驗硬件環(huán)境為32位Windows 7操作系統(tǒng)、 英特爾奔騰雙核G620 @ 2.60 GHz處理器、 4 GB內(nèi)存的計算機. 比較算法以及本文算法SPO濾波部分的參數(shù)根據(jù)相應參考文獻選擇， 本文算法頭部體模中閾值L=0.2， 骨盆體模中L=0.35.

圖 1 給出了各種算法應用于頭部體模而獲得的濾波投影及重建圖像對比. 從圖1可以看出， PRWLS算法和LFEADF算法的濾波投影因噪聲影響顏色比較暗， 其重建圖明顯含有較多噪聲. 而SPO算法和本文算法的濾波投影顏色亮很多， 且其重建圖中噪聲大量減少. 進一步比較圖1(d1)和圖1(e1)可以看出， 本文算法在保持結構方面優(yōu)于SPO算法.

圖 1 頭部體模各種算法濾波投影及重建圖對比Fig.1 Comparison of the projections and the reconstructed images processed by a variety of algorithms on Shepp-Logan phantom

圖 2 給出了頭部體模低劑量噪聲投影與不同算法濾波投影的差值圖. 差值圖中含有的噪聲越多， 邊緣和細節(jié)信息越少， 表明算法效果越好. 圖2(a)中可以看到明顯的邊緣特征， 圖2(b)中邊緣特征雖然較圖2(a)中減少一些， 但是仍然能看到， 說明PRWLS算法和LFEADF算法在濾除噪聲的同時， 部分邊緣和結構也被濾除. 圖2(c)中隱約能看到部分邊緣信息， 圖2(d)中幾乎看不到邊緣殘留， 說明SPO算法和本文算法在保持邊緣和結構方面明顯優(yōu)于前兩個算法， 而本文算法表現(xiàn)又比SPO算法更好.

圖 2 頭部體模低劑量噪聲投影與不同算法濾波投影的差值圖對比Fig.2 Comparison of the image subtractions between low-dose noisy projection and the projection filtered by a variety of algorithms on Shepp-Logan phantom

圖 3 給出了各種算法應用于骨盆體模而獲得的重建圖及其局部感興趣區(qū)域(Regions of Interest, ROI)放大圖， 圖3(a)用矩形框標注了1個ROI. 從圖3可以看出， 本文算法得到的重建圖像去除了更多的噪聲， 更好地保持了圖像的細節(jié)、 邊緣以及結構等特征， 且對比其他算法， 本文算法重建效果更好.

圖 3 骨盆體模各種算法重建圖及局部感興趣區(qū)域放大圖對比Fig.3 Comparison of the reconstructed images and their zoom-in ROIs processed by a variety of algorithms on pelvis phantom

為了定量比較本文算法相較其他算法的優(yōu)越性， 采用均方根誤差(Root Mean Squared Error， RMSE)、 峰值信噪比(Peak Signal to Noise Ratio， PSNR)和結構相似度(Structural Similarity， SSIM)對各重建算法進行了定量描述， 其定義分別如下

(7)

式中：

σuoriginalu=Cov{uoriginal,u}=

在式(7)和式(8)中，un和uoriginaln分別為重建圖與原始圖的灰度值，N為重建圖的像素個數(shù). 式(9)中，c1和c2為根據(jù)文獻[22]設置的常數(shù).

表 1 為PPRWLS算法、 LFEADF算法、 SPO算法以及本文算法分別應用于頭部和骨盆體模而獲得的重建圖像相對于原始圖像的PSNR、 RMSE和SSIM值的比較. 從表 1 可以看出， 本文算法的PSNR值最高， RMSE值最低， 說明本文算法重建圖中含有的噪聲最少， 且最接近原始圖； 而本文算法的SSIM值最高， 說明本文算法重建圖的結構相似度最接近原始圖. 分析實驗中的各圖和表可知， 在主觀的視覺效果和客觀的定量評價方面， 本文算法都優(yōu)于其他算法， 既能降低低劑量CT投影的噪聲， 又能很好地保持圖像結構， 有效地改善了低劑量CT圖像的質量.

表 1 頭部和骨盆體模各種算法的客觀評價

4 結 論

本文提出了一種基于殘差圖分解與平滑塊排列的低劑量CT投影降噪算法. 該方法將SPO算法應用于低劑量CT圖像投影濾波， 并將濾波后所得的殘差圖通過改進的MCA算法分解為結構部分和噪聲部分. 為充分提取殘差圖中的有用信息， 將結構部分反饋于SPO濾波投影， 并對得到的新投影通過FBP算法獲得重建的低劑量CT圖像. 基于不同體模的客觀和主觀的實驗結果都表明， 本文算法無論在噪聲的抑制方面， 還是圖像細節(jié)、 邊緣以及結構等特征的保持方面都明顯優(yōu)于其他算法， 有效地增強了SPO算法性能， 使其更適應用于醫(yī)學圖像處理.

參考文獻：

[1] Brenner D J, Hall E J. Computed tomography——an increasing source of radiation exposure[J]. New England Journal of Medicine, 2007, 357(22): 2277-2284.

[2] Chen Y, Shi L Y, Feng Q J, et al. Artifact suppressed dictionary learning for low-dose CT image processing[J]. IEEE Transactions on Medical Imaging, 2014, 33(12): 2271-2292.

[3] Wang Y L, Shao Y L, Gui Z G, et al. A novel fractional-order differentiation model for low-dose CT image processing[J]. IEEE Access, 2016, 4: 8487-8499.

[4] Chen Y, Ma J H, Feng Q J, et al. Nonlocal prior bayesian tomographic reconstruction[J]. Journal of Mathematical Imaging and Vision, 2008, 30(2): 133-146.

[5] Liu J, Hu Y N, Yang J, Chen Y, et al. 3D feature constrained reconstruction for low dose CT imaging[J]. IEEE Transactions on Circuits and Systems for Video Technology, 2016, 99: 1-14.

[6] Wang J, Li T F, Lu H B, et al. Penalized weighted least-squares approach to sinogram noise reduction and image reconstruction for low-dose X-ray computed tomography[J]. IEEE Transactions on Medical Imaging, 2006, 25(10): 1272-1283.

[7] 劉祎， 張權， 桂志國. 基于模糊熵的低劑量CT投影降噪算法研究[J]. 電子與信息學報， 2013， 35(6): 1421-1427.

Liu Yi, Zhang Quan, Gui Zhiguo. Noise reduction for low-dose CT sinogram based on fuzzy entropy[J]. Journal of Electronics & Information Technology, 2013, 35(6): 1421-1427. (in Chinese)

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

Zhang Quan, Luo Limin, Gui Zhiguo. Improved nonlocal prior based Bayesian sonogram smoothing algorithm for low-dose CT[J]. Journal of Southeast University, 2014, 44(3): 499-503. (in Chinese)

[9] Liu J, Ma J H, Zhang Y, et al. Discriminative feature representation to improve projection data inconsistency for low dose CT imaging[J]. IEEE Transactions on Medical Imaging, 2017, 36(12): 2499-2509.

[10] Ram I, Elad M, Cohen I. Image processing using smooth ordering of its patches[J]. IEEE Transaction on Image Processing, 2013, 22(7): 2764-2774.

[11] Pyo Y, Park R H, Chang S. Noise reduction in high-iso images using 3-D collaborative filtering and structure extraction from residual blocks[J]. IEEE Transactions on Consumer Electronics, 2011, 57(2): 687-695.

[12] Romano Y, Elad M. Improving K-SVD denoising by post-processing its method noise[C]. IEEE International Conference on Imaging Processing, Melbourne, IEEE, 2013: 435-439.

[13] Romano Y, Elad M. Boosting of image denoising algorithms[J]. SIAM Journal on Imaging Sciences, 2015, 8(2): 1187-1219.

[14] 孫偉峰， 戴永壽. 采用多級殘差濾波的非局部均值圖像去噪方法[J]. 電子與信息學報， 2016， 38(8)： 1999-2006.

Sun Weifeng, Dai Yongshou. Non-local means image denoising with multi-stage residual filtering[J]. Journal of Electronics & Information Technology, 2016, 38(8): 1999-2006. (in Chinese)

[15] Starck J K, Elad M, Donoho D L. Redundant muitiscale transforms and their application for morphological component analysis[J]. Advance in Imaging and Electron physics, 2004, 88(4): 1-64.

[16] Lu H B, Hsiao I T, Li X, et al. Noise properties of low-dose CT projections and noise treatment by scale transformations[C]. Nuclear Science Symposium Conference Record, San Diego, IEEE, 2001: 1662-1666.

[17] Cui X Y, Gui Z G, Zhang Q, et al. Learning-based artifact removal via image decomposition for low-dose CT image processing[J]. IEEE Transactions on Nuclear Science, 2016, 63(3): 1860-1873.

[18] Kumar B K S. Image denoising based on non-local means filter and its method noise thresholding[J]. Signal, Image and Video Processing, 2012, 7(6): 1211-1227.

[19] Mairal J, Bach F, Ponce J, et al. Online learning for matrix factorization and sparse coding[J]. Journal of Machine Learning Research, 2010, 11(1): 19-60.

[20] Elad M. Sparse and redundant representations: from theory to applications in signal and image processing[M]. New York: Springer Publishing Company Incorporated, 2010.

[21] Pati Y C, Rezaiifar R, Krishnaprasad P S. Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[C]. Proceedings of the 27th Annual Asilomar Conference on Signals, Systems, and Computers, Pacific Grove, IEEE, 1993: 40-44.

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