基于分數(shù)階各向異性擴散和小波的MLEM低劑量CT重建算法*

焦楓媛， 韓 燮， 常 皓

(1. 中北大學 計算機與控制工程學院， 山西 太原 030051； 2. 中北大學 信息商務(wù)學院， 山西 晉中 030600；3. 中國人民解放軍63961部隊， 北京 100012)

基于分數(shù)階各向異性擴散和小波的MLEM低劑量CT重建算法*

焦楓媛1，2， 韓 燮1， 常 皓3

(1. 中北大學 計算機與控制工程學院， 山西 太原 030051； 2. 中北大學 信息商務(wù)學院， 山西 晉中 030600；3. 中國人民解放軍63961部隊， 北京 100012)

針對低劑量CT重建圖像受噪聲污染嚴重的問題， 提出了一種基于分數(shù)階各項異性擴散和小波的MLEM低劑量CT重建算法. 首先， 采用最大似然期望最大化(Maximum Likelihood Expectation Maximization， MLEM)算法重建低劑量投影數(shù)據(jù). 然后， 將小波變換應用于圖像， 使得圖像的低頻系數(shù)部分集中了主要信息， 而高頻系數(shù)部分集中了邊緣和噪聲. 最后， 在低頻系數(shù)部分進行基于差分的分數(shù)階各項異性擴散， 在高頻系數(shù)部分進行軟閾值處理. 采用Shepp-Logan和Hot-Cold模型進行低劑量CT圖像重建仿真， 實驗結(jié)果表明， 本文算法在降低噪聲和保持圖像細節(jié)方面都優(yōu)于其他算法.

分數(shù)階微分； 各項異性擴散； 小波變換； 低劑量CT； 圖像重建

0 引 言

計算機斷層掃描(Computed Tomography， CT)技術(shù)在疾病預防、 臨床診斷中有著巨大的作用. X射線放射劑量的增高可以提高重建圖像的質(zhì)量， 但這種高劑量帶來的高輻射又會對人體健康造成極大的危害； 而降低劑量， 又必然會使得重建圖像噪聲增加， 進而影響醫(yī)療診斷的準確度. 因此， 如何在降低輻射劑量的同時得到高質(zhì)量的CT重建圖像成為當前的一個熱門問題， 對該問題的研究也越來越受到人們的重視[1].

針對低劑量CT重建圖像的退化問題， 研究學者已經(jīng)提出了很多策略. 一是投影域降噪， 即先對圖像的投影數(shù)據(jù)進行降噪， 再對降噪后的數(shù)據(jù)進行重建. 如Wang等[2]通過在小波域中采用懲罰加權(quán)最小二乘濾波算法大大提高了重建圖像的質(zhì)量； Zhang等[3]提出的基于各項異性擴散的加權(quán)先驗方法實現(xiàn)了對貝葉斯統(tǒng)計的正弦圖的有效降噪， 而且邊緣和細節(jié)信息得到了很好的保持； 李凱旋等[4]提出的基于投影數(shù)據(jù)恢復導引的雙邊濾波權(quán)值優(yōu)化方法進一步保持了圖像細節(jié)信息. 二是圖像后處理方法， 即在重建后對得到的圖像進行降噪. 后處理的方法可以根據(jù)圖像的不同特征而采取不同的方法， 由于其方便實用而成為近年來研究的一個熱點. 許多有效的濾波算法被應用于低劑量CT圖像處理中， 并取得了較好的效果. 如Chen等[5-6]提出了一種基于字典學習的低劑量CT圖像偽影去除算法和降低腦部CT輻射劑量的字典學習算法， 有效改善了低劑量CT圖像的質(zhì)量； Rust等[7]利用非線性高斯濾波器鏈對重建圖像進行平滑處理； Paul等[8]提出的基于特征的低劑量胸部CT圖像噪聲降低方法， 可以對不同的組織結(jié)構(gòu)采取不同的降噪方法， 從而在減少偽影的同時很好地保持圖像的細節(jié)信息； Zhang等[9-10]首次將分數(shù)階微積分應用于醫(yī)學圖像處理， 并提出了兩個用于降低CT圖像金屬偽影的分數(shù)階方程， 使得圖像細節(jié)信息得到了進一步保持.

基于以上分析和研究， 本文提出了一種基于分數(shù)階各項異擴散和小波的MLEM 低劑量CT重建算法. 該算法首先對低劑量投影數(shù)據(jù)進行最大似然期望最大化算法重建， 然后對重建圖像進行小波分解， 在小波低頻系數(shù)部分， 采用基于方差的分數(shù)階各項異性擴散， 在高頻系數(shù)部分， 運用閾值降噪進行處理. 采用Shepp-Logan和Hot-Cold模型進行仿真實驗， 通過與其他算法進行比較， 本文算法重建圖像不但取得了最高的信噪比， 而且在降低噪聲的同時在視覺效果上也有效地保持了圖像的細節(jié)和邊緣信息.

1 分數(shù)階微積分

分數(shù)階微積分有空域和頻域兩大類定義， 其中應用最廣泛的是Grümwald-Letniklv (G-L)定義和 Riemann- Liouville (R-L) 定義. 本文在圖像處理中應用到的是G-L定義， 信號f(x)的α階微分的G-L定義如下[11]

2 基于分數(shù)階各項異性擴散和小波的MLEM低劑量CT重建算法

2.1 最大似然期望最大化重建算法

最大似然期望最大化(Maximum Likelihood Expectation Maximization, MLEM)算法通常會在重建過程中考慮到觀測數(shù)據(jù)的統(tǒng)計特性和系統(tǒng)的物理模型， 因此， 與濾波反投影(Filtexed Back Projection, FBP)重建算法比較， 該算法重建效果明顯更好， 其重建公式如下[12]

2.2 基于分數(shù)階各項異性擴散和小波的降噪方法

2.2.1 整數(shù)階各項異性擴散算法

Perona等[13]在1990年提出了各項異性擴散方程(PM模型)， 其表達式為

式中：c為擴散系數(shù)函數(shù)， 其表達式為

c(|f|)=1/[1+|

式中：k為梯度閾值. PM模型能夠在有效去除噪聲的同時保持圖像邊緣信息， 但去噪后的圖像往往會出現(xiàn)塊狀效應.

為了克服塊狀效應， You等[14]在2000年提出了一類四階偏微分方程， 其在一定程度上抑制了塊狀效應， 但是又產(chǎn)生了斑點效應.

2.2.2 基于方差的分數(shù)階各項異性擴散算法

為了解決以上問題， Bai等[15]在2007年提出了分數(shù)階各項異性擴散方程， 其能量函數(shù)表示為

當α=1時， 式(5) 近似于 PM模型， 當α=2時， 式(5)近似于四階擴散方程， 當1<α<2時， 式(5)可以看作二階和四階各項異性擴散的一個插值. 分數(shù)階各項異性擴散算法在降低噪聲的同時有效抑制了階梯效應和斑點效應.

為了更好地保持圖像細節(jié)和邊緣信息， 本文將可以較好區(qū)分細節(jié)和平坦區(qū)域的局部方差加入到分數(shù)階各項異性擴散方程的擴散系數(shù)中. 由文獻[16]可知歸一化的方差可表示為

其中：

此時， 分數(shù)階各項異性擴散模型中擴散系數(shù)為

式中：k1為常數(shù).

應用文獻[17]中分數(shù)階微分的離散化方法， 將分數(shù)階微分看作卷積積分， 于是式(5)變?yōu)?/p>

其中：

利用梯度下降法， 可以得到

首先對(vα*f)x和(vα*f)y進行離散化

記：

這里為了保證分母不為零，ε取很小的整數(shù). 基于方差的分數(shù)階各項異性擴散方程的離散化形式為

2.2.3 小波閾值降噪

噪聲和表示圖像特征的細節(jié)信息都屬于高頻成份， 因此在降低噪聲的圖像處理過程中往往會帶來細節(jié)信息的丟失. 針對這一問題， 人們在圖像去噪中廣泛采用小波變換， 因為其具有時頻局部化和多分辨率特性. 而在小波域中， 表示噪聲的系數(shù)會隨著尺度的增加而減小， 表示信號的系數(shù)會隨著尺度的增加而增加[18]. 所以， 小波分解之后， 絕對值比較大的高頻小波系數(shù)需要保留， 因為其主要由原始圖像信號提供； 絕對值比較小的高頻小波系數(shù)需要濾除， 因為其主要由噪聲提供. 在小波閾值降噪過程中， 通過對小波系數(shù)進行閾值處理可以使重建圖像中有用的信息得到很好的保持， 因此， 合適的閾值選擇至關(guān)重要.

2.3 重建算法

經(jīng)過小波變換之后， 低頻部分集中了圖像的主要信息， 而邊緣和噪聲信息屬于高頻分量. 由于噪聲對分數(shù)階各項異性擴散影響較大， 所以在含有噪聲較少的低頻分量中采用分數(shù)階各項異性擴散來降噪， 就可以避免噪聲對其造成的影響. 另外， 分數(shù)階各項異性擴散需要較大的計算量， 而小波變換收斂速度很快. 因此， 本文結(jié)合了二者的優(yōu)點， 在降噪過程中運用了混合降噪方法， 即在低頻分量部分采用基于方差的分數(shù)階各項異性擴散， 在高頻分量部分采用小波變換軟閾值去噪. 本文重建算法的具體步驟如下：

1) 在理想投影數(shù)據(jù)中加入噪聲變?yōu)榈蛣┝客队皵?shù)據(jù)， 再根據(jù)式(2)對其進行最大似然期望最大化算法重建.

2)在每次迭代中， 將多尺度小波變換應用于MLEM重建圖像中， 由此生成的低頻分量為CAi， 高頻分量記為CHi,CVi,CDi， 其中i=1,2,3,…,n， 為分解尺度.

3) 選擇合適的閾值， 對高頻系數(shù)進行小波閾值降噪.

4) 將基于方差的分數(shù)階各項異性擴散降噪算法應用于低頻系數(shù)中. 由式(12)～式(14)， 可得迭代公式為

5) 進行小波反變換得到去噪后的圖像.

6) 不斷重復以上步驟， 最終獲得理想的重建圖像.

3 實驗結(jié)果及分析

3.1 重建圖像比較

為了直觀地證明本文算法的有效性， 仿真實驗中采用了2個模型數(shù)據(jù)： 模型1為Shepp-Logan頭部剖面圖模型， 模型2為模擬人體骨骼組織和軟組織的Hot-Cold 模型， 模型大小均為128×128， 系統(tǒng)矩陣大小為 (128×128) ×(128×128). Li等[19]的實驗結(jié)果表明： 低劑量投影數(shù)據(jù)的噪聲是期望和方差成非線性關(guān)系的高斯噪聲， 其表達式如下

如圖 1 所示， 本文除了將傳統(tǒng)的MLEM、 BI- MLEM(Block Iterative-Maximum Likelihood Expectation Maximization)[20]重建算法與本文算法進行比較， 還將經(jīng)典的PM、 TV算法應用到MLEM重建算法中， 分別形成了PM-MLEM和TV-MLEM算法， 并與本文算法進行了比較.

圖 1 Shepp-Logan模型比較實驗Fig.1 The comparative experiments on Shepp-Logan model

通過觀察圖 1 中Shepp-Logan頭部模型原圖及其重建圖像可以發(fā)現(xiàn)， MLEM和BI- MLEM重建圖像含有大量噪聲， 視覺效果非常差， 所以這兩種算法不適合應用于低劑量CT圖像重建. PM-MLEM和TV-MLEM重建圖像中包含更少的噪聲， 但是出現(xiàn)了明顯的塊狀效應. 從視覺效果上看， 本文算法重建圖像中包含的噪聲和條形偽影最少， 基本看不到塊狀效應， 圖像細節(jié)和紋理保持也比較好. 因此， 本文算法在低劑量CT圖像重建中的表現(xiàn)明顯優(yōu)于其它4種算法.

為了進一步驗證本文算法的可行性， 采用Hot-Cold 模型進行實驗. 從圖 2 中Hot-Cold 模型原圖及其重建圖像可以看到， MLEM重建圖像含有最多的噪聲， 重建效果很差. BI-MLEM重建圖像仍然含有大量噪聲， 而PM-MLEM和TV- MLEM重建圖像中噪聲大幅下降， 但是能看到明顯的條形偽影和塊狀效應. 本文算法重建圖像在噪聲降低和細節(jié)、 紋理保持方面都是最好的， 視覺效果最接近原圖. 比較其他算法， 本文算法在降低噪聲的同時保持了更多的細節(jié)信息， 因此更適合應用于低劑量CT圖像處理.

圖 2 Hot-Cold模型比較實驗Fig.2 The comparative experiments on Hot-Cold model

3.2 重建精度分析

從上述視覺效果分析可知， 本文算法對低劑量CT圖像有很好的處理效果. 為了進一步客觀上定量地說明本文算法在醫(yī)學圖像處理中的可行性， 采用歸一化均方誤差、 均方絕對誤差、 信噪比對重建圖像質(zhì)量進行定量描述. 其定義表述如下：

1) 歸一化均方誤差(NMSE)

2) 均方絕對誤差(MAE)

).

(17)

3) 信噪比(SNR)

式中：Mi和mi分別為重建圖像與真實圖像的均值；Fi和qi分別為重建圖像與真實圖像的第i個像素的灰度值；M和N分別為圖像的行數(shù)和列數(shù). NMSE和MAE值越小， 重建圖像越接近于真實圖像； SNR的值越高， 重建圖像質(zhì)量越好.

本文算法與其他幾種算法應用于Shepp-Logan模型和Hot-Cold 模型的客觀評價結(jié)果如表 1 所示.

表 1 Shepp-Logan和Hot-Cold 模型各種算法的客觀評價

從表 1 可以明顯看出， 本文算法的信噪比SNR最大， 說明重建圖像中含有最少的噪聲， 而本文算法的歸一化均方誤差NMSE和均方絕對誤差MAE均為最小， 說明本文算法重建圖像與真實圖像最接近. 因此， 本文算法重建圖像質(zhì)量是最好的. 從圖1、 圖2和表1可以看出， 本文算法無論在視覺效果還是在定量分析方面均優(yōu)于其他算法.

為了更加直觀地比較本文算法與其他4種算法重建圖像與原始圖像之間的差距， 本文給出了Shepp-Logan模型的原始理想圖像與另外幾種重建算法結(jié)果圖像的側(cè)面輪廓線比較圖， 并比較了第65行的灰度值， 如圖 3 所示.

圖 3 各種算法重建圖像第65行截面圖對比Fig.3 Comparison of the profiles on line 65 of processed image by using various algorithms

從圖 3 可以看出， 與其他幾種算法相比， 本文算法重建出的圖像最接近理想圖像， 并且在降低噪聲和保持圖像細節(jié)方面效果都最好.

4 結(jié) 論

本文基于分數(shù)階各項異性擴散和小波閾值， 提出了一種基于分數(shù)階各項異性擴散和小波的MLEM低劑量CT重建算法. 由于分數(shù)階各項異性擴散對噪聲比較敏感， 所以本文算法在對低劑量投影數(shù)據(jù)進行MLEM重建之后， 對重建圖像進行小波分解， 然后在低頻分量部分采用基于方差的分數(shù)階各項異性擴散， 在高頻分量部分采用小波變換軟閾值去噪. 與此同時， 傳統(tǒng)的各項異性擴散容易引起塊狀效應， 本文采用分數(shù)階各項異性擴散， 有效去除了塊狀效應. 另外， 本文在擴散系數(shù)中加入方差項， 進一步保持了圖像細節(jié)和紋理. 仿真實驗結(jié)果表明， 本文算法在降低噪聲和保持圖像細節(jié)方面都優(yōu)于其他算法. 因此， 基于分數(shù)階各項異性擴散和小波的MLEM重建算法適合應用于醫(yī)學圖像處理.

[1] Bian J G, Yang K, Boone J M, et al. Investigation of iterative image reconstruction in low-dose breast CT[J]. Physics in Medicine and Biology, 2014, 59(11)： 2659-2685.

[2] Wang J, Lu H B. Multiscale penalized weighted least-squares sonogram restoration for low-dose X-Ray computed tomography[J]. IEEE Transactions on Biomedical Engineering, 2008, 55(3)： 1022-1032．

[3] Zhang Q, Gui Z G, Chen Y, et al. Bayesian sinogram smoothing with an anisotropic diffusion weighted prior for low-dose X-ray computed tomography[J]. Optik-International Journal for Light and Electron Optics, 2013, 124(17)： 2811-2816.

[4] 李凱旋, 黃靜, 馬建華, 等. 低劑量CT重建中的雙邊濾波權(quán)值優(yōu)化新方法[J]. 電路與系統(tǒng)學報， 2012， 17(2)： 94-99.

Li Kaixuan, Huang Jing, Ma Jianhua, et al. Bilateral filtering weights optimization in low-dose computed tomography image reconstruction[J]. Journal of Circuits and Systems, 2012, 17(2)： 94-99. (in Chinese)

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

[6] Chen Y, Shi L, Yang J, et al. Radiation dose reduction with dictionary learning based processing for head CT[J]. Australasian Physical & Engineering Sciences in Medicine, 2014, 37(3)： 483-493.

[7] Rust G F, Aurich V, Reiser M. Noise dose reduction and image improvements in screening virtual colonoscopy with tube currents of 20 mAs with nonlinear Gaussian filter chains[C]. Medical Imaging 2002 Conference, New York： IEEE, 2002： 186-197.

[8] Paul N S, Blobelj J, Prezelj E, et al. The reduction of image noise and steak artifact in the thoracic inlet during low dose and ultra-low dose thoracic CT[J]. Physics in Medicine and Biology, 2010, 55(5)： 1363-1380.

[9] Zhang Y, Pu Y F, Hu J R, et al. A new CT metal artifacts reduction algorithm based on fractional-order sinogram inpainting[J]. Journal of X-Ray Science and Technology, 2011, 19(3)： 373-384.

[10] Zhang Y, Pu Y F, Hu J R, et al. Efficient CT metal artifacts reduction based on fractional-order curvature diffusion[J]. Computational & Mathematical Methods in Medicine, 2011(12)： 173748.

[11] Pu Y F, Wang W X, Zhou J L, et al. Fractional differential approach to detecting texture features of digital image and its fractional differential filter implementation[J]. Science China Information Sciences, 2008, 51(9)： 1319-1339.

[12] 張權(quán), 付學敬, 桂志國, 等. 基于小波和各向異性擴散的PET圖像MLEM重建算法[J]. 計算機應用與軟件, 2013, 30(11)： 50-53.

Zhang Quan, Fu Xuejing, Gui Zhiguo, et al. MLEM reconstrution algorithm for PET image based on wavelet and anisotropic diffusion[J]. Computer Applications and Software, 2013, 30(11)： 50-53. (in Chinese)

[13] Perona P, Malik J. Scale space and edge detection using anisotropic diffusion[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1990, 12(7)： 629-639.

[14] You Y L, Kaveh M. Fourth-order partial differential equations for noise removal[J]. IEEE Transactions on Image Processing, 2000, 9(10)： 1723-1730.

[15] Bai J, Feng X C. Fractional-order anisotropic diffusion for image denoising[J]. IEEE Transactions on Image Processing, 2007, 16(10)： 2492-2502.

[16] Chao S M, Tsai D M. An improved anisotropic diffusion model for detail and edge-preserving smoothing Pattern Recogn[J]. Pattern Recognition Letters, 2010, 31(13)： 2012-2023.

[17] 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.

[18] Liu J, Moulin P. Information-theoretic analysis of interscale and intrascale dependencies between image wavelet coefficients[J]. IEEE Transactions on Image Processing, 2001, 10(11)： 1647-1658.

[19] Li T, Li X, Wang J, et al. Nonlinear sonogram smoothing for low-dose X-ray CT[J]. IEEE Transactions on Nuclear Science, 2004, 51(5)： 2505-2513.

[20] Byrne C L. Block-iterative methods for image reconstruction from projections[J]. IEEE Transactions on Image Processing, 1996, 5(5)： 792-794.

TheMLEMLow-DoseCTReconstructionAlgorithmBasedonFractional-OrderAnisotropicDiffusionandWavelet

JIAO Feng-yuan1,2, HAN Xie1, CHANG Hao3

(1. School of Computer Science and Control Engineering, North University of China, Taiyuan 030051, China；2. College of Information and Business, North University of China, Jinzhong 030600, China; 3. PLA Unit 63961, Beijing 100012, China)

Concerning the noise problem about reconstructed images of low-dose CT, a MLEM low-dose CT reconstruction algorithm based on fractional-order anisotropic diffusion and wavelet was proposed. Firstly, the MLEM algorithm was used to do reconstruction of the low-dose projection data. Secondly, by the wavelet transformation, the main information of image was concentrated in the low frequency coefficients, and the edge and noise was stored in the high frequency coefficients. Finally, the low frequency coefficients of the image were filtered by fractional-order anisotropic diffusion based on difference, and the high frequency coefficients were processed by the soft threshold function. Shepp-Logan and Hot-Cold model were used for low-dose CT reconstruction simulation. The experimental results demonstrate that compared with other methods, the proposed algorithm achieves superior performance in terms of both noise suppression and detail preservation.

fractional-order differentiation; anisotropic diffusion; wavelet transform; low-dose computed tomography(CT); image reconstruction

1673-3193(2017)03-0348-06

2016-08-16

山西省重點研發(fā)計劃資助項目(201603D121012)

焦楓媛(1992-)， 女， 碩士生， 主要從事圖像處理及計算軟件的研究.

TP391.41

A

10.3969/j.issn.1673-3193.2017.03.017