董嬋嬋，張 權(quán)，郝慧艷，劉 祎

（中北大學(xué) 電子測(cè)試技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，山西 太原030051）

0 引 言

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

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

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

1 噪聲模型

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

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

2 基于小波變換和各項(xiàng)異性擴(kuò)散的MLEM低劑量CT重建算法

2.1 MLEM重建方法

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

2.2 重建過程中的降噪方法

2.2.1 各項(xiàng)異性擴(kuò)散算法

傳統(tǒng)P-M模型的表達(dá)式為［13］

式中：ct為擴(kuò)散系數(shù)函數(shù)，其表達(dá)式為

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

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

式中：fηη和fξξ分別表示圖像f在梯度方向和水平方向的二階導(dǎo)數(shù)，其表達(dá)式為

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

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

此時(shí)，各項(xiàng)異性擴(kuò)散方程為

2.2.3 小波閾值降噪

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

2.3 算法描述

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

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

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

3）對(duì)高頻系數(shù)進(jìn)行軟閾值處理，去除噪聲；對(duì)低頻系數(shù)用基于差分曲率各項(xiàng)異性擴(kuò)散算法進(jìn)行圖像降噪處理.

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

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

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

3 實(shí)驗(yàn)結(jié)果與討論

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

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

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

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

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

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

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

2）均方絕對(duì)誤差（Mean Absolute Error，MAE）

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

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

在公式（11）～（13）中，F(xiàn)i和qi分別表示重建圖像與原始圖像的第i個(gè)像素的灰度值；Mi和mi分別表示重建圖像與原始圖像的均值；M和N分別為圖像的行數(shù)和列數(shù).

本文提出的新算法與各種算法比較的客觀評(píng)價(jià)結(jié)果如表1所示.

表1 各種算法的客觀評(píng)價(jià) Tab.1 The objective evaluation of various algorithms

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

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

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

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

4 結(jié) 論

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

［1］唐瑜，崔學(xué)英，張權(quán)，等.基于偏微分方程的低劑量CT投影降噪算法［J］.計(jì)算機(jī)應(yīng)用與軟件，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］高志凌，劉祎，桂志國(guó).基于模糊數(shù)學(xué)的低劑量CT投影域降噪算法［J］.測(cè)試技術(shù)學(xué)報(bào)，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］劉祎，張權(quán)，桂志國(guó).基于模糊熵的低劑量CT投影降噪算法研究［J］.電子與信息學(xué)報(bào)，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重建中的雙邊濾波權(quán)值優(yōu)化新方法［J］.電路與系統(tǒng)學(xué)報(bào)，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］.太原：中北大學(xué)，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］張權(quán)，付學(xué)敬，桂志國(guó)，等.基于小波和各向異性擴(kuò)散的PET圖像MLEM重建算法［J］.計(jì)算機(jī)應(yīng)用與軟件，2013，30（11）：50-53.Zhang Quan，F(xiàn)u 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］王知強(qiáng).基于小波收縮與非線性擴(kuò)散的去噪算法［J］.計(jì)算機(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.