基于改進各項異性擴散和小波變換的MLEM低劑量CT重建算法

董嬋嬋，張 權，郝慧艷，劉 祎

（中北大學 電子測試技術國家重點實驗室，山西 太原030051）

0 引 言

醫(yī)學計算機斷層成像術（Computerized Tomography，CT）技術是通過測量物體不同角度下的射線投影，從而獲得物體橫截面信息的成像技術，被廣泛應用于醫(yī)療診斷、無損檢測等領域［1］.X射線的放射劑量越高，重建圖像的質量越好.但在進行CT檢查時，高劑量的輻射會對人體的健康造成危害；輻射劑量降低時，獲得的投影數據會伴隨著比較高的電子噪聲，重建出的圖像質量往往不能滿足圖像處理和醫(yī)療診斷的需求［2］.因此，在低劑量下重建出高分辨率和低噪聲的CT圖像的研究越來越受到關注.

在低劑量CT圖像重建過程中，由于探測器接受到的光子數較少，投影數據受噪聲污染嚴重，導致重建后圖像的質量較差［3］.研究至今，在對低劑量CT重建圖像進行降噪時，主要方法有兩類：一是先對圖像的投影數據進行投影域降噪，然后再用降噪后的投影數據進行重建；二是先重建，再對重建后得到的圖像進行圖像域降噪處理.到目前為止，針對投影域數據進行降噪的研究取得了很多成果.如Wang Jing等利用懲罰加權最小二乘法方法［4］，分別在圖像域、投影域以及在投影數據的K-L（Karhunen-Loeve）域中進行了懲罰加權最小二乘法（Penalized Weighted Least-Squares，PWLS），取得了不錯的效果；2008年，Lu Hongbing等人［5］在小波域中進行PWLS濾波算法，通過小波變換來對投影數據進行處理，大大提高了重建圖像的質量；李凱旋［6］等提出一種基于投影數據恢復導引的雙邊濾波權值優(yōu)化方法，在對噪聲去除的同時有效的保護了圖像的細節(jié)；Zhang Quan［7］等提出一種新穎的基于各項異性擴散加權先驗，通過對貝葉斯統計的正弦圖進行平滑，從而保持了圖像的邊緣.但對投影域降噪處理的方法計算復雜且計算量大，同時重建算法需要對系統的噪聲、探測器與X線源間校準等過程進行準確重建，增加了優(yōu)化難度，而優(yōu)化重建算法難度較高不利于在高速CT系統中應用［8］.而在圖像域進行降噪處理的方法比較方便直接，可以根據圖像不同的特性提出不同的降噪方法.因此，近年來在圖像域進行降噪處理成為了研究熱點，并取得了較好的效果.如Rust等利用非線性高斯濾波器鏈對重建圖像進行平滑處理［9］；Chen Yang［10］等通過使用一種新穎的非局部自適應加權非局部先驗統計重建方法，改善了低劑量CT圖像的質量.

針對以上研究分析，本文提出了一種基于小波收縮和差分曲率各項異性擴散的最大似然期望最大化（Maximum Likelihood Expectation Maximization，MLEM）的低劑量CT重建算法，先對投影數據進行重建，再將重建得到的圖像進行小波分解，對高頻系數進行閾值降噪處理，對小波低頻系數進行基于差分曲率的各項異性擴散.

1 噪聲模型

低劑量CT圖像的噪聲表現在圖像上是一些噪聲脈沖和條狀偽影，從而使圖像的信噪比下降嚴重，以及圖像中的一些細節(jié)被污染而無法顯現.低劑量投影數據的噪聲特點是低劑量CT圖像重建解決以上問題的關鍵.長期以來，投影數據的噪聲都被認為符合泊松分布，這是基于光子數目符合泊松分布而確定的.而Li［11］等人經過多次研究分析認為低劑量的投影數據通過系統校正以及對數變換之后，其均值與方差之間是非線性關系，近似服從非平穩(wěn)高斯分布.其滿足如下公式

式中：i=1，2，…，M，表示的是探測器的信道，M為信道的總數；λi和σ2i分別表示探測器信道處獲得的投影數據的均值和方差；fi表示信道探測器的參數；T是系統參數，它是一個用于描述系統校準過程的尺度系數.不同的CT采集系統，fi與T這兩個參數也不相同，但是對于一個給定的CT采集系統，fi與T是已知的.

2 基于小波變換和各項異性擴散的MLEM低劑量CT重建算法

2.1 MLEM重建方法

MLEM算法由于在重建過程中既考慮了系統的物理模型，又考慮了觀測數據的統計特性，重建出的圖像要優(yōu)于FBP重建出的圖像，其重建公式為［12］

2.2 重建過程中的降噪方法

2.2.1 各項異性擴散算法

傳統P-M模型的表達式為［13］

式中：ct為擴散系數函數，其表達式為

式中：K為常數.傳統P-M模型僅利用梯度來控制模型的擴散速度，雖然可以有效地去除噪聲，但圖像中有些平坦區(qū)域的梯度與細節(jié)處的梯度相差無幾，因此只用梯度來檢測邊緣是不夠的.2.2.2 基于差分曲率的各項異性擴散

文獻［14］提出的差分曲率可以較好地區(qū)分細節(jié)和平坦區(qū)域、獨立噪聲，其表達式為

式中：fηη和fξξ分別表示圖像f在梯度方向和水平方向的二階導數，其表達式為

由式（5）可知：在邊緣處，|fηη|較大，|fξξ|較小，其差分曲率D值較大；在平坦區(qū)域，|fηη|和|fξξ|都較小，所以D值較??；在獨立噪聲點處，|fηη|和|fξξ|都較大，且?guī)缀跸嗟?，所以D較小.因此根據差分曲率D值大小，可以很好地區(qū)分邊緣和平坦區(qū)域、獨立噪聲點.因此，將差分曲率引入后的擴散系數函數為

其中，dt，N為歸一化差分曲率

此時，各項異性擴散方程為

2.2.3 小波閾值降噪

小波變換是用小波系數的形式來表示信號，描述不同尺度的信號變化情況.由于小波變換具有時頻局部化和多分辨率的特性，因此圖像去噪在小波域中的研究成為熱點.小波變換能將信號的能量集中到少數小波系數上［15-16］，而白噪聲在任何正交基上的變換仍然是白噪聲，并且有著相同的幅度.相對而言，信號的小波系數值必然大于那些能量分散且幅值較小的噪聲的小波系數值.文獻［17］中分析了小波系數在層間存在較強的持續(xù)性，小波域中信號系數隨著尺度增加而增加，而噪聲系數隨著尺度的增加而減小.由于噪聲主要集中在圖像的高頻部分，因此選擇一個合適的閾值，對小波系數進行閾值處理，從而達到去除噪聲而保留有用信號的目的.

2.3 算法描述

一方面，圖像經過小波變換后，邊緣和噪聲是高頻分量，圖像的主要信息都集中在低頻部分.由于各項異性擴散對噪聲較敏感，而低頻分量中的噪聲較少，因此在低頻部分進行各項異性擴散可以減小噪聲對其的影響.另一方面，由參考文獻［17］可知，小波收縮方法收斂速度快，而各項異性擴散需要較多次迭代，計算量較大.由于經過多尺度小波分解后，低頻部分的大小遠小于原圖像.因此，本文在降噪部分采用混合去噪算法，即在高頻分量部分采用小波收縮軟閾值去噪，在低頻采用基于差分曲率的各項異性擴散.算法的具體步驟如下：

1）低劑量投影數據是按式（1）的方法加在理想投影數據上，然后按式（2）進行MLEM算法重建.

2）在每次重建迭代中，首先對上步重建后的圖像信號進行多尺度小波變換，生成相應的低頻分量CAi，高頻分量CHi，CVi，CDi，i=1，2，3，…，n為分解尺度.

3）對高頻系數進行軟閾值處理，去除噪聲；對低頻系數用基于差分曲率各項異性擴散算法進行圖像降噪處理.

4）在小波域進行圖像降噪處理后，進行小波反變換得到去噪后的圖像.

5）對各項異性擴散處理后的脈沖噪聲進行中值濾波處理：由文獻［18］可知，低劑量CT重建圖像的噪聲還表現為一些脈沖噪聲，各項異性擴散降噪技術對重建的圖像進行降噪后，僅可以平滑圖像的小梯度區(qū)域，而相對于周圍區(qū)域的大梯度區(qū)域則保持不變，這些大梯度可能是邊緣，也可能是圖像的峰值噪聲.而中值濾波器只會對大噪聲峰值產生的大梯度起作用，邊緣產生的大梯度將不會受到影響.因此在低劑量CT重建時，低噪聲可以由各項異性擴散平滑，而大噪聲等脈沖噪聲則由中值濾波器所消除.中值濾波公式為fn+1i，j=Median（fn+1i，j，w），其中w是中值算子的窗口.

6）重復步驟1）～5）直至得到最終的重建圖像.

3 實驗結果與討論

為了驗證算法的有效性，本文采用模擬人體骨骼組織和軟組織的Hot-Cold模型作為實驗模型進行低劑量CT圖像重建的仿真，本文與BIMART（Block Iterative-Multiply algebraic reconstruction technique），BI-MLEM（Block Iterative-Maximum Likelihood Expe-ctation Maximization）［19］進行了比較，且為了進一步說明基于小波收縮和差分曲率的各項異性擴散算法在MLEM的低劑量CT重建算法的有效性，本文將傳統P-M、方差應用到MLEM算法中，并和本文方法進行了比較.體模的大小為128 mm×128 mm，按照式（1）的方法向理想投影數據加入期望與方差的非線性關系的高斯噪聲來仿真低劑量的投影數據，其中fi=200，T=1 200；實現本文算法的計算機操作系統為32位Microsoft Windows XP Professional 2002，處理器為英特爾奔騰雙核E5300＠2.60 GHz，內存為2 G.計算結果如圖1所示.

圖1 各種算法重建圖像的對比 Fig.1 Comparison of reconstructed image processed by various algorithms

由圖1（b）可以明顯看出，MLEM重建算法的重建結果圖中含有較多噪聲，即該算法不能有效地解決低劑量重建圖像受噪聲污染嚴重的問題.圖1（c）和圖1（d）與圖1（b）相比可以看出：算法BI-MART和BI-EMML比MLEM算法重建后的圖像包含更少的噪聲，處理效果較好.將圖1（g）的結果與圖1（b）、（c）、（d）、（e）和（f）的處理結果相比可以看出：本文算法包含的噪聲更少，圖像更清晰，重建結果明顯優(yōu)于其它4種算法.綜上所述，本文算法可以有效地解決低劑量重建圖像的噪聲問題，且可以在光滑去噪的同時較好地保持圖像的紋理和邊緣信息.

圖2 肩部仿真模型重建 Fig.2 Comparison diagram of algorithm results

圖2為肩部仿真模型的重建結果.將圖2（e）與圖2（a）、（b）、（c）和（d）相比，可以清晰地看出：（e）中包含的噪聲較少，圖像較平滑、清晰，圖像質量較高，即本文算法處理后的圖像明顯優(yōu)于其它算法處理后的圖像.因此，本文算法可以有效地解決低劑量重建圖像的噪聲問題，且可以在光滑去噪的同時較好地保持圖像的紋理和邊緣信息.

從上述分析可得，本文提出的算法在去噪能力與保持邊緣和細節(jié)能力方面，都明顯優(yōu)于其他重建方法.為了更好地證明算法的有效性，本文除了采用主觀觀察方法之外，對重建的結果也進行了客觀、定量的分析.采用歸一化均方誤差、均方絕對誤差、歸一化均方距離、信噪比等評價參數對其進行評價.其定義分別為

1）歸一化均方誤差（Normalized Mean Squared Error，NMSE）

2）均方絕對誤差（Mean Absolute Error，MAE）

3）歸一化均方距離（Normalized Mean Square Distance，NMSD）

4）信噪比（Signal to Noise Ratio，SNR）

在公式（11）～（13）中，Fi和qi分別表示重建圖像與原始圖像的第i個像素的灰度值；Mi和mi分別表示重建圖像與原始圖像的均值；M和N分別為圖像的行數和列數.

本文提出的新算法與各種算法比較的客觀評價結果如表1所示.

表1 各種算法的客觀評價 Tab.1 The objective evaluation of various algorithms

由表1可以明顯看出，本文算法的歸一化均方誤差、均方絕對誤差、歸一化均方距離比其他幾種方法都要小，說明了重建圖像與原始圖像的偏差和誤差都較小，相似度較高.信噪比比其他的幾種算法都大，說明重建圖像中含有的噪聲最少.由以上分析可知，本文算法得到的重建圖像更接近原始理想圖像.因此結合圖2和表1可知，無論在視覺方面還是在定量評價方面，本算法在CT重建中是可行且有效的.

為了更好地比較4種算法的降噪效果，本文給出了肩部仿真模型的以上重建圖像與原始理想圖像的側面輪廓線的比較圖，選取第65行的灰度值進行比較，如圖3所示.

圖3 各種算法重建圖像的截面圖對比 Fig.3 Comparison of profiles of reconstructed image processed by various algorithms

從圖3中可以清晰地看出，本文算法重建出的圖像與理想圖像的吻合度最高，更接近理想圖像，即與其他幾種算法相比，該重建圖的噪聲波動最小，去除噪聲以及保留圖像的細節(jié)信息的效果最好.

4 結 論

本文提出了一種基于小波變換和差分曲率各項異性擴散的MLEM的低劑量CT重建算法.該算法是采用基本的MLEM算法對低劑量投影數據進行重建，然后對重建后的圖像在小波域進行降噪處理.通過與其他重建方法以及不同降噪方法的仿真實驗相比較，無論是從主觀的視覺效果還是從客觀的質量評價，均表明本文提出的算法能夠有效地去除噪聲并保持圖像的細節(jié)信息.算法可以應用于醫(yī)學CT領域中，能夠在降低輻射劑量的同時，有效地重建出符合要求的圖像.因此，算法在醫(yī)學成像領域具有較好的應用前景.

［1］唐瑜，崔學英，張權，等.基于偏微分方程的低劑量CT投影降噪算法［J］.計算機應用與軟件，2014，31（3）：179-183.Tang Yu，Cui Xueying，Zhang Quan，et al.Denoising algorithm for low-dose CT projection based on partial differential equation［J］.Computer Applications and Software，2014，31（3）：179-183.（in Chinese）

［2］高志凌，劉祎，桂志國.基于模糊數學的低劑量CT投影域降噪算法［J］.測試技術學報，2011，25（6）：477-482.Gao Zhiling，Liu Yi，Gui Zhiguo.Noise reduction based on fuzzy theory for low-does CT sinograms［J］.Journal of Test and Measurement Technology，2011，25（6）：477-482.（in Chinese）

［3］劉祎，張權，桂志國.基于模糊熵的低劑量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 ffInformation Technology，2013，35（6）：1421-1427.（in Chinese）

［4］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.

［5］Wang Jing，Lu Hongbing.Multiscale penalized weighted least-squares sonogram restoration for low-dose Xray computed tomography［J］.IEEE Transactions on Biomedical Engineering，2008，55（3）：1022-1032.

［6］李凱旋，黃靜，馬建華，等.低劑量CT重建中的雙邊濾波權值優(yōu)化新方法［J］.電路與系統學報，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）

［7］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.

［8］高志凌.基于投影域降噪的低劑量CT重建算法研究［D］.太原：中北大學，2012.

［9］Rust G F，Aurich V，Reiser M.Noise dose reduction and image improvements in screening virtual colonoscopy with tube currents of 20 m As with nonlinear Gaussian filter chains［C］.Medical Imaging 2002 Conference.New York：IEEE，2002：186-197.

［10］Chen Y，Gao D Z，Nie C，et al.Bayesian statistical reconstruction for low-dose X-ray computed tomography using an adaptive-weighting local nonprior［J］.Computerized Medical Imaging and Graphics，2009，33（7）：495-500.

［11］Li T，Li X，Wang J，et al.Nonlinear sinogram smoothing for low-dose X-ray CT［J］.IEEE Transactions on Nuclear Science，2004，51（5）：2505-2513.

［12］張權，付學敬，桂志國，等.基于小波和各向異性擴散的PET圖像MLEM重建算法［J］.計算機應用與軟件，2013，30（11）：50-53.Zhang Quan，Fu Xuejing，Gui Zhiguo，et al.MLEMreconstrution 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］Chen Q，Montesinos P，Sun Q，et al.Adaptive total variation denoising based on difference curvature［J］.Image and Vision Computing，2010，28（3）：298-306.

［15］Zhong J，Sun H.Wavelet-based multiscale anisotropic diffusion with adaptive statistical analysis for image restoration［J］.IEEE Transactions on Circuits and Systems，2008，55（9）：2716-2725.

［16］王知強.基于小波收縮與非線性擴散的去噪算法［J］.計算機工程，2011，37（7）：249-252.Wang Zhiqiang.Denoising algorithm based on wavelet shrinkage and nonlinear diffusion［J］.Computer Engineering，2011，37（7）：249-252.（in Chinese）

［17］Juan Liu，Moulin，P.Information-theoretic analysis of interscale and intrascale dependencies between image wavelet coefficients［J］.IEEE Transactions Image Processing，2001，10（11）：1647-1658.

［18］Ling J，Bovik A C.Smoothing low-SNR molecular images via anisotropic median-diffusion［J］.IEEE Transactions on Medical Imaging，2002，21（4）：377-384.

［19］Byrne C L.Block-iterative methods for image reconstruction from projections［J］.IEEE Transactions on image processing，1996，5（5）：792-794.