基于小波和四階各向異性擴(kuò)散的MLEM低劑量CT重建算法

張　芳　桂志國(guó),2*　張　權(quán)　董嬋嬋

1(中北大學(xué)信息與通信工程學(xué)院電子測(cè)試技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室　山西 太原 030051)2(中北大學(xué)信息與通信工程學(xué)院儀器科學(xué)與動(dòng)態(tài)測(cè)試教育部重點(diǎn)實(shí)驗(yàn)室　山西 太原 030051)

?

基于小波和四階各向異性擴(kuò)散的MLEM低劑量CT重建算法

張芳1桂志國(guó)1,2*張權(quán)1董嬋嬋1

1(中北大學(xué)信息與通信工程學(xué)院電子測(cè)試技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室山西 太原 030051)2(中北大學(xué)信息與通信工程學(xué)院儀器科學(xué)與動(dòng)態(tài)測(cè)試教育部重點(diǎn)實(shí)驗(yàn)室山西 太原 030051)

在低劑量計(jì)算機(jī)斷層掃描CT(computedtomography)重建算法中，傳統(tǒng)的最大似然期望最大MLEM(MaximumLikelihoodExpectationMaximization)算法隨著迭代次數(shù)的增加會(huì)出現(xiàn)棋盤(pán)效應(yīng)而不能有效地抑制噪聲。針對(duì)上述問(wèn)題提出一種基于小波收縮和四階各向異性擴(kuò)散相結(jié)合的MLEM低劑量CT重建算法。該算法結(jié)合小波收縮和各向異性擴(kuò)散的優(yōu)點(diǎn)，在每次迭代中，對(duì)MLEM重建算法處理后的圖像進(jìn)行離散平穩(wěn)小波分解，在小波域的高頻部分進(jìn)行小波收縮，低頻部分使用降噪效果優(yōu)質(zhì)的四階各向異性擴(kuò)散進(jìn)行消噪，最后殘留的脈沖噪聲點(diǎn)通過(guò)中值濾波器進(jìn)行處理，從而進(jìn)一步優(yōu)化圖像。仿真實(shí)驗(yàn)結(jié)果表明，該算法可以有效地去除低劑量CT圖像的噪聲，且在保持圖像邊緣和細(xì)節(jié)信息方面有很好的表現(xiàn)，從而獲得高抗噪性能的圖像。

低劑量計(jì)算機(jī)斷層掃描四階各向異性擴(kuò)散圖像重建中值濾波

0　引　言

X射線計(jì)算機(jī)斷層CT成像以其成像分辨率高的優(yōu)點(diǎn)，在醫(yī)學(xué)臨床診斷中已獲得廣泛應(yīng)用。然而，因其較高的輻射劑量會(huì)對(duì)人體造成一定的傷害，所以在掃描時(shí)應(yīng)盡可能地降低射線劑量，即采用低劑量CT的成像技術(shù)。然而，低劑量CT重建會(huì)使圖像質(zhì)量發(fā)生明顯的退化，為了解決上述問(wèn)題,目前國(guó)內(nèi)外學(xué)者嘗試各種各樣的處理方法來(lái)提高低劑量CT重建圖像的質(zhì)量，使其達(dá)到滿意的效果。

針對(duì)低劑量CT重建圖像的退化問(wèn)題，目前主要方法可歸于以下兩種方式：1) 在重建后對(duì)得到的圖像進(jìn)行圖像域降噪。2) 為先對(duì)圖像的投影數(shù)據(jù)進(jìn)行投影域降噪，然后再用降噪后的投影數(shù)據(jù)進(jìn)行重建。近年來(lái)，針對(duì)前者方法學(xué)者們提出很多方法，且取得了令人滿意的效果。如ChenYang[1]等通過(guò)使用一種新穎的非局部自適應(yīng)加權(quán)非局部先驗(yàn)統(tǒng)計(jì)重建方法改善了低劑量CT圖像的質(zhì)量；Rust[2]等通過(guò)使用非線性高斯濾波器鏈對(duì)低劑量CT圖像進(jìn)行濾波，從而得到保護(hù)了邊緣信息的降噪處理；LuiD[3]提出一種新穎的噪聲補(bǔ)償CT重建方法，從而提高了重建圖像的信噪比。對(duì)于后者方法，即先對(duì)圖像的投影數(shù)據(jù)進(jìn)行投影域降噪，然后再用降噪后的投影數(shù)據(jù)進(jìn)行重建，也取得了很大的成果。如由于傳統(tǒng)的濾波反投影算法FBP(FilteredBackProjection)的算法簡(jiǎn)單，且速度較快，故該類算法一般是利用低劑量投影數(shù)據(jù)的特點(diǎn)設(shè)計(jì)濾波算法，濾除噪聲后再利用FBP進(jìn)行重建，從而獲得良好的CT圖像質(zhì)量。GuiZhiguo[4]等使用模糊濾波器對(duì)低劑量CT的投影數(shù)據(jù)進(jìn)行降噪處理。ZhangQuan[5]等對(duì)基于各向異性擴(kuò)散加權(quán)先驗(yàn)的貝葉斯正弦圖進(jìn)行平滑。劉祎[6]等通過(guò)把局部模糊熵的自適應(yīng)算法加入到投影數(shù)據(jù)中，從而得到高信噪比的重建圖像。馬建華[7]等通過(guò)對(duì)投影數(shù)據(jù)進(jìn)行統(tǒng)計(jì)建模，然后采用貝葉斯最大后驗(yàn)估計(jì)的方法，將其非局部的先驗(yàn)信息加入到投影數(shù)據(jù)中，從而實(shí)現(xiàn)投影數(shù)據(jù)降噪的效果。WangJing[8]等利用懲罰加權(quán)最小二乘法方法，分別在圖像域、投影域以及在投影數(shù)據(jù)的K-L(Karhunen-Loeve)域中進(jìn)行了懲罰加權(quán)最小二乘法PWLS(PenalizedWeightedLeast-Squares)，取得了不錯(cuò)的效果。

雖然投影域降噪的重建方法取得了較大的發(fā)展，但由于在投影域進(jìn)行降噪后，即使投影域中很小的點(diǎn)噪聲映射到圖像域，也會(huì)形成條形偽影。因此本文提出一種基于四階各向異性擴(kuò)散和小波收縮的的低劑量CT重建算法。該算法重建出來(lái)的圖像不僅能夠有效地保持圖像邊緣的細(xì)節(jié)信息，而且不存在各向異性擴(kuò)散存在的明顯階梯效應(yīng)。

1　本文重建算法

1.1MLEM重建方法

最大似然期望最大化MLEM算法由于在重建過(guò)程中既考慮了系統(tǒng)的物理模型，又考慮了觀測(cè)數(shù)據(jù)的統(tǒng)計(jì)特性，因而被廣泛的使用，其重建公式如下式：

(1)

1.2重建過(guò)程中的降噪方法

1.2.1四階偏微分方程去噪算法

基于偏微分方程PDE(PartialDifferentialEquation)的各向異性擴(kuò)散的去噪方法，可以滿足整幅圖像中需要去噪的不同強(qiáng)度的需求，因此是一種自適應(yīng)的去噪技術(shù)。該方法的原理是在圖像的平滑區(qū)域增加平滑的強(qiáng)度，而在圖像的邊緣位置則適當(dāng)?shù)販p弱平滑強(qiáng)度從而保持圖像的邊緣區(qū)域，從而可以更好地在濾除噪聲的同時(shí)抑制圖像邊緣的過(guò)平滑。由于二階PDE降噪會(huì)出現(xiàn)“階梯”效應(yīng)，因此后人提出了四階PDE模型[9]，其表達(dá)式為：

(2)

其中，fηη和fξξ為圖像梯度方向和切線方向的二階導(dǎo)數(shù)，表達(dá)式如下：

(3)

(4)

式(2)的離散化形式為：

(5)

四階PDE既可以克服二階PDE產(chǎn)生的“階梯”效應(yīng)，又可以根據(jù)梯度和切線方向進(jìn)行不同程度的擴(kuò)散，進(jìn)而有效地保持邊緣。

1.2.2小波變換

目前，圖像降噪的方法有很多，小波變換作為一種新的多分辨分析方法，在圖像去噪領(lǐng)域被廣泛應(yīng)用。小波變換能將信號(hào)的能量集中到少數(shù)小波系數(shù)上，而白噪聲在任何正交基上的變換仍然是白噪聲，并且有著相同的幅度。文獻(xiàn)[10]分析小波系數(shù)在層間存在較強(qiáng)的持續(xù)性，小波域中信號(hào)系數(shù)隨著尺度增加而增加，而噪聲系數(shù)隨著尺度的增加而減小。由于噪聲和圖像邊緣主要集中在圖像的高頻部分，因此可以采用小波收縮法達(dá)到圖像去噪的目的。

f(k)=s(k)+n(k)k=0,1,2…,N-1

(6)

對(duì)f(k)做離散小波變換，可得：

wf(j,k)=ws(j,k)+wn(j,k)

(7)

其中，wf(j,k)、ws(j,k)和wn(j,k)分別表示包含噪聲信號(hào)、原始信號(hào)和噪聲信號(hào)在第j層上的小波系數(shù)。由式(7)可以看出，觀測(cè)信號(hào)的小波系數(shù)wf(j,k)，由原始信號(hào)小波系數(shù)ws(j,k)和噪聲信號(hào)小波系數(shù)wn(j,k)兩部分組成。對(duì)包含噪聲信號(hào)的小波系數(shù)，若它的值大于指定的閾值，則此系數(shù)為信號(hào)分量，予以保留；若其值小于閾值，則認(rèn)為此系數(shù)為噪聲分量，濾除這樣的系數(shù)就可達(dá)到降噪效果。

上述方法為小波閾值去噪，又稱小波收縮[11]，它對(duì)小波系數(shù)進(jìn)行統(tǒng)一處理，在最小均方誤差下達(dá)到近似最優(yōu)。和正交小波變換相比，平穩(wěn)小波變換具有平移不變性而更適合用于圖像的處理[12，13]。因此本文選用平穩(wěn)小波變換對(duì)圖像進(jìn)行處理。其中2j為相應(yīng)的尺度。

1.3本文提出的重建算法

由于小波變換具有良好的時(shí)頻局部特性，小波變換之后，低頻部分所含的噪聲比較少。高頻部分主要包含了圖像的噪聲和邊緣，且噪聲對(duì)應(yīng)的小波系數(shù)幅值在噪聲處比較小，在邊緣處比較大。四階偏微分對(duì)噪聲的敏感性高，噪聲越少，降噪效果越好。小波閾值收縮對(duì)處理此處高頻噪聲有很好的效果，故把四階偏微分降噪用在低頻處理，小波閾值用在高頻處理是合適的。根據(jù)上述分析，本文的流程圖1所示。

圖1　本文算法流程圖

本文的具體重建算法如下：

1)MLEM重建算法：

(8)

2) 小波收縮和四階各向異性降噪處理：

在每次重建迭代中，首先對(duì)上步重建后的圖像信號(hào)進(jìn)行多尺度小波變換，生成相應(yīng)的低頻分量CAi，高頻分量CHi,CVi,CDi，i=1,2,3,…,n為分解尺度。然后，對(duì)高頻系數(shù)進(jìn)行軟閾值處理，去除噪聲；對(duì)低頻系數(shù)用基于四階各向異性擴(kuò)散算法進(jìn)行圖像降噪處理。所需公式如下：

(9)

(10)

其中，式(9)為軟閾值收縮函數(shù)，ω表示含噪的小波系數(shù)，δ(ω)表示去噪后的小波系數(shù)，T為閾值，經(jīng)過(guò)與其他幾種常用的閾值做對(duì)比試驗(yàn)，證明本文選取的長(zhǎng)度對(duì)數(shù)閾值效果最優(yōu)且能達(dá)到令人滿意的效果。故本文使用長(zhǎng)度對(duì)數(shù)閾值，即對(duì)小波分解的高頻小波系數(shù)按式(9)進(jìn)行處理，對(duì)低頻部分按式(10)進(jìn)行四階各向異性擴(kuò)散處理，最后進(jìn)行平穩(wěn)小波逆變換重構(gòu)，得到優(yōu)化后的離散圖像。

3) 進(jìn)一步中值濾波處理：

由文獻(xiàn)[14]可知，低劑量CT重建圖像的噪聲還表現(xiàn)為一些脈沖噪聲。各向異性擴(kuò)散降噪技術(shù)對(duì)重建的圖像進(jìn)行降噪后，僅可以平滑圖像的小梯度區(qū)域，而相對(duì)于周圍區(qū)域的大梯度區(qū)域則保持不變，這些大梯度可能是邊緣，也可能是圖像的峰值噪聲。而中值濾波器只會(huì)對(duì)大噪聲峰值產(chǎn)生的大梯度起作用，邊緣產(chǎn)生的大梯度將不會(huì)受到影響。因此在低劑量CT重建時(shí)，低噪聲可以由各向異性擴(kuò)散平滑，而大噪聲等脈沖噪聲則由中值濾波器所消除。

fi,jn + 1=Median(fi,jn + 1,w)

(11)

所取窗口太大，不僅會(huì)增加算法的運(yùn)算量，而且會(huì)使圖像過(guò)平滑，所取窗口太小，會(huì)達(dá)不到降噪的效果，經(jīng)過(guò)反復(fù)試驗(yàn)，這里取w為3×3的窗口。

4) 重復(fù)以上步驟，反復(fù)調(diào)試，選取視覺(jué)效果最優(yōu)的結(jié)果作為最終的重建圖像。本文提出算法的部分代碼如下：

fork=1:150

F=(F·*(ggg*(p./(G*F))));

g=reshape(F,128,128);

%對(duì)MLEM重建后的圖像進(jìn)行平穩(wěn)離散小波變換

[swa,swh,swv,swd] =swt2(x,3,′db1′);

[thr] =ddencmp(′den′,′wv′,x);

sorh= ′s′;

%對(duì)高頻成分進(jìn)行軟閾值處理

dswh=wthresh(swh,sorh,thr);

dswv=wthresh(swv,sorh,thr);

dswd=wthresh(swd,sorh,thr);

h=0.8;

%對(duì)低頻成分進(jìn)行四階各向異性降噪處理

fork1=1:3

g=swa(:,:,k1);

k=4;

fori=2:M-1

forj=2:N-1

gx=(g(i+1,j)-g(i-1,j))/2;

gy=(g(i,j+1)-g(i,j-1))/2;

gxx=g(i+1,j)+g(i-1,j)-2*g(i,j);

gyy=g(i,j+1)+g(i,j-1)-2*g(i,j);

gxy=(g(i+1,j+1)+g(i-1,j-1)-g(i+1,j-1)-g(i-1,j+1))/4;

tidu=sqrt(gx*gx+gy*gy);

gc=(gx*gx*gxx+2*gy*gx*gxy+gy*gy*gyy)/(tidu*tidu+eps);

gs=(gy*gy*gxx-2*gy*gx*gxy+gx*gx*gyy)/(tidu*tidu+eps);

C=(k*k)/(k*k+(tidu*tidu));

l(i,j)=C*C*gc+C*gs;

end

end

fori=2:M-2

forj=2:N-2

t(i,j)=l(i+1,j)+l(i-1,j)+l(i,j+1)+l(i,j-1)-4*l(i,j);

g(i,j)=g(i,j)-dt*t(i,j);

end

end

swa(:,:,k1)=g;

end

xd=iswt2(swa,dswh,dswv,dswd,′db1′);

2　實(shí)驗(yàn)結(jié)果與分析

2.1重建圖像比較

本文算法及所有比較算法的實(shí)驗(yàn)仿真環(huán)境為：Windows7 旗艦版32位SP1(DirectX11)操作系統(tǒng)，英特爾Celeron(賽揚(yáng))E3300@2.50GHz雙核處理器，2GB內(nèi)存。編程工具使用MATLAB7.6.0(R2008a)。本文選取大小為128mm×128mm的Sheep-Logan頭部剖面圖模型作為實(shí)驗(yàn)對(duì)象，其灰度范圍為[0,255]，如圖2(a)所示。實(shí)驗(yàn)的無(wú)噪聲投影數(shù)據(jù)通過(guò)在180個(gè)角度中均勻選取128個(gè)投影角度，每個(gè)投影角度下有128條射線的方法來(lái)獲得，故本文選取大小為(128×128)×(128×128)的系統(tǒng)矩陣，實(shí)驗(yàn)中采用平行投影方式，在180個(gè)角度中均勻采取128個(gè)投影方向，每個(gè)方向分布128個(gè)探測(cè)器對(duì)。本文按照下式向理想投影數(shù)據(jù)加入期望和方差成非線性關(guān)系的高斯噪聲來(lái)仿真低劑量CT投影數(shù)據(jù)：

σi2=kiexp(λi/T)

(12)

其中，i=1,2，…，N表示探測(cè)器信道，N表示信道總數(shù)，λi表示第i個(gè)探測(cè)器獲得投影數(shù)據(jù)的均值，σi2表示第i個(gè)探測(cè)器獲得投影數(shù)據(jù)的方差，ki表示第i個(gè)探測(cè)器的參數(shù)，T表示用于描述掃描系統(tǒng)校準(zhǔn)過(guò)程的系統(tǒng)參數(shù)。對(duì)于給定的CT采集系統(tǒng)，ki與T是給定的。ki越大，T越小，低劑量CT重建圖像的噪聲越大；反之ki越小，T越大，低劑量CT重建圖像的噪聲越小。若噪聲太大，會(huì)使圖像失真，從而得不到滿意的效果；若噪聲太小，不能達(dá)到和實(shí)際數(shù)據(jù)的比擬，失去了研究的實(shí)際意義。為了更好地模擬實(shí)際低劑量CT的臨床數(shù)據(jù)，經(jīng)過(guò)大量的實(shí)驗(yàn)經(jīng)驗(yàn)與對(duì)比，本文噪聲的參數(shù)選?。簁i=200，T=12 000。為了驗(yàn)證本文算法的有效性，將本文算法與傳統(tǒng)的MLEM、BI-PLS(blockiteration-PenalizedLeastSquares)[15]、OS-PML-OSL(OrderedSubsets-PenalizedMaximumLikelihood-OneStepLate)和在每次重建迭代中使用基于四階的各向異性擴(kuò)散降噪的算法進(jìn)行比較。各種算法的對(duì)比結(jié)果如圖2所示。

圖2　各種算法的對(duì)比結(jié)果

圖2中，(a)為原始圖像，(b)采用傳統(tǒng)的MLEM重建算法，(c) 采用基于塊迭代的BI-PLS重建算法，(d)采用基于有序子集的OS-PML-OSL重建算法。為了把本文算法和降噪效果很好的各向異性擴(kuò)散算法進(jìn)行比較，(e)在每次迭代中使用基于四階的各向異性擴(kuò)散降噪算法對(duì)重建圖像進(jìn)行處理，(f)為本文提出的算法。各種算法中涉及到的各種參數(shù)以及迭代次數(shù)，均為反復(fù)實(shí)驗(yàn)后得到的最優(yōu)值，從結(jié)果圖中可得，傳統(tǒng)的MLEM重建圖像的質(zhì)量最差；BI-PLS和OS-PML-OSL算法的結(jié)果圖可以達(dá)到降噪的效果，但同時(shí)也模糊了圖像的邊緣和細(xì)節(jié)；基于四階的各向異性擴(kuò)散降噪算法對(duì)重建圖像的噪聲進(jìn)行了一定程度的抑制且可獲得比較清晰的圖像，但是圖像中存在一些比較明顯的塊狀陰影；本文算法在消除噪聲的同時(shí)很好地保持了圖像的邊緣和細(xì)節(jié)信息，使圖像達(dá)到比較優(yōu)質(zhì)的效果，從視覺(jué)上分析，重建效果和其他幾種比較算法相比達(dá)到最優(yōu)，初步說(shuō)明了本文算法的有效性。

2.2重建精度比較

從2.1節(jié)分析可知，本文提出的算法對(duì)低劑量CT的重建圖像有很好的降噪以及保持邊緣和細(xì)節(jié)信息的能力，為了進(jìn)一步說(shuō)明本文算法的有效性，本文采用歸一化均方距離、均方絕對(duì)誤差、信噪比對(duì)重建圖像質(zhì)量進(jìn)行定量描述。其定義如下：

1) 歸一化均方距離NMSD(NormalizedMeanSquareDistance):

(13)

2) 均方絕對(duì)誤差MAE(MeanAbsoluteError):

(14)

3) 信噪比SNR(SignaltoNoiseRatio):

(15)

其中，J表示圖像像素點(diǎn)的總數(shù)，F(xiàn)i和fi分別表示重建圖像與原始圖像的第i個(gè)像素的灰度值，Mi和mi分別表示重建圖像與原始圖像的平均值。這些指標(biāo)分別從不同的方面評(píng)價(jià)重建圖像與原始圖像的接近程度以及重建圖像的質(zhì)量，表1為本文算法與其他幾種比較算法的客觀評(píng)價(jià)結(jié)果。

表1　各種算法的客觀評(píng)價(jià)

從表1可得出本文算法的信噪比最大，其他指標(biāo)值均比其他的幾種比較算法小。該結(jié)論說(shuō)明本文算法的重建圖像和原始圖像最為接近。因此在定量評(píng)價(jià)方面，同樣表明本算法在低劑量CT重建中的可行性。

圖3給出了本文所用的Sheep-Logan模型的原始理想圖像與在以上各種算法的重建圖像的側(cè)面輪廓線的比較圖。從圖中可以看出本文算法的重建圖像與原始圖像的吻合度最高，最接近于理想圖像，具有最小的噪聲波動(dòng)，可以在保持圖像邊緣和細(xì)節(jié)信息的同時(shí)較好地解決低劑量CT重建圖像的噪聲問(wèn)題。

圖3　各種算法第65行側(cè)面輪廓線的對(duì)比結(jié)果

3　結(jié)　語(yǔ)

本文針對(duì)傳統(tǒng)的MLEM重建算法收斂速度慢，不能很好地控制噪聲的問(wèn)題，提出了一種基于小波收縮和四階各向異性擴(kuò)散的MLEM的低劑量CT重建算法。其中小波收縮可以在去除低劑量圖像噪聲的同時(shí)很好地保持圖像的邊緣和細(xì)節(jié)信息；基于四階的各向異性擴(kuò)散可以有效地判斷圖像的背景和邊緣，然后進(jìn)行不同程度的擴(kuò)散；在低劑量CT重建中，小噪聲可以由上述方法處理，而大噪聲等脈沖噪聲則可以利用中值濾波器去除。本文算法結(jié)合了以上各種算法的優(yōu)勢(shì)，既有效地保持了圖像的細(xì)節(jié)和邊緣信息，又很好地解決了低劑量圖像的噪聲問(wèn)題。實(shí)驗(yàn)結(jié)果表明，該算法在主觀效果和客觀效果來(lái)看，均說(shuō)明本文算法是可行的。在臨床中該算法可以提高低劑量CT重建圖像的質(zhì)量，使其達(dá)到滿意的效果，使醫(yī)生在患者接受較少的輻射劑量下便可進(jìn)行良好的診斷和治療。本算法不僅可以應(yīng)用于醫(yī)學(xué)CT的醫(yī)學(xué)成像領(lǐng)域，也可用于工業(yè)CT領(lǐng)域，同時(shí)在醫(yī)學(xué)PET和工業(yè)無(wú)損檢測(cè)中也有良好的應(yīng)用前景。

[1]ChenY,GaoDZ,NieC,etal.Bayesianstatisticalreconstructionforlow-doseX-raycomputedtomographyusinganadaptive-weightinglocalnonprior[J].ComputerizedMedicalImagingandGraphics,2009,33(7): 495-500.

[2]RustGF,AurichV,ReiserM.Noisedosereductionandimageimprovementsinscreeningvirtualcolonoscopywithtubecurrentsof20mAswithnonlinearGaussianfilterchains[C]//MedicalImaging2002Conference,NewYork:IEEE,2002:186-197.

[3]LuiD,CameronA,ModhafarA,etal.Low-dosecomputedtomographyviaspatiallyadaptiveMonte-Carloreconstruction[J].ComputerizedMedicalImagingandGraphics,2013,37(7-8): 438-449.

[4]GuiZG,LiuY.Noisereductionforlow-doseX-raycomputedtomographywithfuzzyfilter[J].Optik-InternationalJournalforLightandElectronOptics,2012,123(13):1207-1211.

[5]ZhangQ,GuiZG,ChenY,etal.Bayesiansinogramsmoothingwithananisotropicdiffusionweightedpriorforlow-doseX-raycomputedtomography[J].Optik-InternationalJournalforLightandElectronOptics,2013,124(17):2811-2816.

[6] 劉祎，張權(quán)，桂志國(guó).基于模糊熵的低劑量CT投影降噪算法研究[J].電子與信息學(xué)報(bào),2013,35(6):1421-1427.

[7] 馬建華，黃靜，陳陽(yáng)，等.基于廣義Gibbs先驗(yàn)的低劑量X-CT優(yōu)質(zhì)重建研究[J].計(jì)算機(jī)工程與應(yīng)用,2008,44(16):4-6,93.

[8]WangJ,LiTF,LuHB,etal.PenalizedWeightedLeast-SquaresApproachtoSinogramNoiseReductionandImageReconstructionforLow-DoseX-RayComputedTomography[J].IEEETransactionsonmedicalimaging,2006,25(10):1272-1283.

[9]HajiaboliMR.AnAnisotropicFourth-OrderDiffusionFilterforImageNoiseRemoval[J].InternationalJournalofComputerVision,2011,92(2):177-191

[10]LiuJ,MoulinP.Information-theoreticanalysisofinterscaleandintrascaledependenciesbetweenimagewavelet-coefficients[J].IEEETransImageProcessing,2001,10(11):1647-1658.

[11]LeeK,VidakovicB.Semi-supervisedwaveletshrinkage[J].ComputationalStatistics&DataAnalysis,2012，56(6)：1681-1691.

[12]DengJM,YueHZ,ZhuoZZ,etal.AstationarywavelettransformbasedapproachtoregistrationofplanningCTandsetupconebeam-CTimagesinradiotherapy[J].JournalofMedicalSystems,2014,38(5)：1-9.

[13]JumahAA.Denoisingofanimageusingdiscretestationarywavelettransformandvariousthresholdingtechniques[J].JournalofSignalandInformationProcessing,2013,4(1)：33-41.

[14]LingJ,BovikAC.Smoothinglow-SNRmolecularimagesviaanisotropicmedian-diffusion[J].IEEETransactionsonMedicalImaging,2002,21(4):377-384.

[15]ByrneCL.Block-Iterativemethodsforimagereconstructionfromprojections[J].IEEETransactionsonimageprocessing,1996,5(5):792-794.

MLEMLOW-DOSECTRECONSTRUCTIONBASEDONWAVELETANDFOURTH-ORDERANISOTROPICDIFFUSION

ZhangFang1GuiZhiguo1,2*ZhangQuan1DongChanchan1

1(National Key Laboratory for Electronic Measurement Technology,School of Information and Communication Engineering,North University of China,Taiyuan 030051,Shanxi,China)2(Key Laboratory of Instrumentation Science and Dynamic Measurement,School of Information and Communication Engineering,North University of China,Taiyuan 030051,Shanxi,China)

Inlow-doseCT(computedtomography)reconstructionalgorithm,traditionalmaximumlikelihoodexpectationmaximisation(MLEM)algorithmwillappearchessboardeffectalongwiththeincreaseofthenumberofiterations,thuscannoteffectivelysuppressnoises.Forthisproblem,thispaperproposesalowdoseCTreconstructionalgorithm,whichisbasedonthecombinationofwaveletshrinkageandfourth-orderanisotropicdiffusion.Itcombinestheadvantagesofwaveletshrinkageandanisotropicdiffusion,ineachiteration,itconductsthediscretestationarywaveletdecompositionontheimageprocessedwithMLEMreconstructionalgorithm,inhighfrequencypartofthewaveletdomainitshrinksthewavelet,inlowfrequencypartitusesfourth-orderanisotropicdiffusion,whichhashighqualityeffectindenoising,toeliminatenoises,itprocessesthefinalresidualpulsenoisepointswiththemedianfilter,sothatfurtheroptimisestheimage.Simulationexperimentresultsshowthattheproposedalgorithmcaneffectivelyremovethenoiseinlow-doseCTimage,andhasgoodperformancesinkeepingboththeimageedgesanddetailedinformation,therebygainstheimagewithhighanti-noiseperformance.

Low-dosecomputedtomographyFourth-orderAnisotropicdiffusionImagereconstructionMedianfilter

2014-09-03。國(guó)家自然科學(xué)基金項(xiàng)目(61071192，61271357，61171178)；山西省國(guó)際合作項(xiàng)目(2013081035)；山西省研究生優(yōu)秀創(chuàng)新項(xiàng)目(2009011020-2，20123098)；中北大學(xué)第十屆研究生科技基金項(xiàng)目(20131035)；山西省高等學(xué)校優(yōu)秀青年學(xué)術(shù)帶頭人支持計(jì)劃資助項(xiàng)目；中北大學(xué)2013年?？茖W(xué)基金計(jì)劃。張芳，碩士生，主研領(lǐng)域：基于低劑量CT的圖像重建。桂志國(guó)，教授。張權(quán)，副教授。董嬋嬋，碩士生。

TP391

ADOI:10.3969/j.issn.1000-386x.2016.03.043