一種自適應(yīng)低劑量CT圖像質(zhì)量改善算法

蔣慧琴，徐玉風(fēng)，馬 嶺，楊曉鵬，Toshiya Nakaguchi

(1.鄭州大學(xué) 信息工程學(xué)院，河南 鄭州 450001； 2.鄭州大學(xué)第一附屬醫(yī)院 醫(yī)學(xué)裝備部，河南 鄭州 450052； 3.日本千葉大學(xué) 先端醫(yī)工學(xué)研究中心, Chiba)

0 引言

為了解決CT檢查導(dǎo)致的輻射劑量過高問題，低劑量CT掃描逐漸成為臨床應(yīng)用熱點(diǎn).然而，低劑量CT掃描會(huì)導(dǎo)致圖像質(zhì)量退化，特別是LDCT圖像中的噪聲會(huì)嚴(yán)重影響診斷精度.因此，去除噪聲是改善低劑量CT圖像質(zhì)量的重要任務(wù).低劑量CT掃描所引起的噪聲主要是X射線量子噪聲[1].針對(duì)圖像中的細(xì)節(jié)和噪聲同屬高頻成分難以分離的經(jīng)典難題，小波變換等多尺度分析方法呈現(xiàn)出一定的優(yōu)勢(shì)[2].近幾年，許多學(xué)者已提出大量的小波域去除量子噪聲方法[3-4].主要包括兩種思路：其一，假設(shè)其服從Gaussian分布，利用閾值法去噪，如小波硬閾值、軟閾值和貝葉斯閾值法等[3].其二，假設(shè)其近似服從Poisson分布.首先利用方差穩(wěn)定變換將其轉(zhuǎn)化為近似服從Gaussian分布的噪聲，再利用去除Gaussian噪聲的常用方法進(jìn)行處理[4-7].隨著多尺度分析方法的發(fā)展，一些學(xué)者提出了基于剪切波變換的圖像去噪算法.例如，Dan等[8]提出的基于剪切波變換的貝葉斯最大后驗(yàn)估計(jì)方法，能有效去除因小波去噪產(chǎn)生的偽吉布斯效應(yīng)；胡海智等[9]提出的全變差正則化剪切波域收縮方法，能更有效地去除Gaussian噪聲并保留圖像中的細(xì)節(jié)信息等等.為了去除低劑量CT圖像中的量子噪聲，筆者假設(shè)量子噪聲近似服從Poisson分布，并利用Anscombe變換，將原圖像中的噪聲轉(zhuǎn)化為一個(gè)具有近似常數(shù)方差的Gaussian噪聲.然而，由于傳統(tǒng)Anscombe變換對(duì)低信噪比的圖像其噪聲方差無法趨于穩(wěn)定[10]，致使去噪效果欠佳.為了解決這一問題，文獻(xiàn)[11]通過最優(yōu)化Anscombe逆變換以改進(jìn)去噪效果；Zhang等[12]通過改進(jìn)Anscombe變換，使當(dāng)Poisson噪聲參數(shù)λ>0.1時(shí)，其噪聲方差可以趨于穩(wěn)定.但是，當(dāng)λ≤0.1時(shí)，該方法仍不能使噪聲方差穩(wěn)定到一個(gè)常數(shù)值.基于前期研究[7]可知，當(dāng)Poisson噪聲參數(shù)λ≤0.1時(shí)，如在應(yīng)用Anscombe變換的基礎(chǔ)上，再利用閾值法去除量子噪聲可以取得更好的效果，但其去噪效果依賴于噪聲方差值的估計(jì)精度.

鑒于此，筆者針對(duì)低信噪比的剪切波高頻系數(shù)子帶，提出基于剩余自相關(guān)功率的噪聲方差估計(jì)方法，并結(jié)合貝葉斯最大后驗(yàn)估計(jì)提取出非噪聲高頻系數(shù)，從而有效地提高低信噪比圖像的去噪效果.

1 剪切波變換

圖1 剪切波分解算法框架Fig.1 The block diagram of shearlet decomposition

圖2 剪切波分解的一層圖像Fig.2 The shearlet coefficients at one scale

2 改進(jìn)的噪聲方差估計(jì)算法

雖然Donoho提出的噪聲方差計(jì)算方法(MAD方法)被廣泛應(yīng)用,但因其是針對(duì)小波域估計(jì)高斯噪聲方差而設(shè)計(jì)，存在一定誤差.為此，剩余自相關(guān)功率(RAP)[10]與貝葉斯最大后驗(yàn)估計(jì)相結(jié)合的改進(jìn)噪聲方差估計(jì)精度的方法，其實(shí)現(xiàn)過程如下.

2.1 貝葉斯最大后驗(yàn)估計(jì)

設(shè)

g=f+w，

(1)

式中：f和g分別表示原圖像和加噪后的圖像；w表示均值為0的Gaussian白噪聲.

設(shè)

y=x+n，

(2)

基于貝葉斯最大后驗(yàn)估計(jì)計(jì)算該近似值的步驟如下.

(1)用Donoho提出的估計(jì)方法計(jì)算噪聲方差近似值:

(3)

其中，yi,j表示選取的高頻子帶系數(shù).

(2)計(jì)算提取非噪聲剪切波高頻子帶系數(shù)的閾值如下：

(4)

(5)

(6)

(3)用貝葉斯最大后驗(yàn)估計(jì)法提取高頻子帶中非噪聲剪切波系數(shù)近似值：

(7)

2.2 噪聲方差估計(jì)

對(duì)CT圖像加入均值為0的Gaussian白噪聲，用非下采樣剪切波對(duì)圖像進(jìn)行尺度為4的分解，每層選定一個(gè)方向.由于隨著分解層數(shù)的增加，高頻子帶的圖像信息逐漸減少，而噪聲信息逐漸增多，且第4層高頻子帶的信噪比最低，故筆者分別用MAD方法和剩余自相關(guān)功率方法計(jì)算含噪圖像第4層高頻細(xì)節(jié)子帶的噪聲方差，并對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行比較，其實(shí)現(xiàn)步驟如下.

(1)選取噪聲方差的候選值.利用式(3)計(jì)算所選高頻子帶的噪聲方差σm，并且以該值為中心,間隔為0.1, 選取30個(gè)噪聲方差候選值∑=[σ1,σ2,…,σm,…,σ30].

(2)計(jì)算所選高頻子帶所需的閾值.針對(duì)選定的每一個(gè)噪聲方差候選值σm，利用式(4)所示的貝葉斯最大后驗(yàn)估計(jì)方法計(jì)算閾值[8].

(8)

(9)

其中,N是自相關(guān)中點(diǎn)的數(shù)量.當(dāng)真正的噪聲方差為20時(shí)，基于不同的噪聲方差候選值計(jì)算出其對(duì)應(yīng)的Pm，其對(duì)應(yīng)關(guān)系如圖3所示.從圖3可以看出，其對(duì)應(yīng)的Pm在接近真正的噪聲方差時(shí)有一個(gè)突變，且之后Pm趨于穩(wěn)定.

(5)計(jì)算噪聲方差.為了捕捉Pm的這種特性，記

Dm=Pm+1-Pm.

(10)

不同的噪聲方差候選值與Dm的對(duì)應(yīng)關(guān)系如圖4所示，從圖4可以看出，真正的噪聲方差即當(dāng)Dm達(dá)到最大值時(shí)對(duì)應(yīng)的σm右側(cè)一個(gè)較小的值，這個(gè)值的估計(jì)公式為:

mmax=argmmaxDm.

(11)

(12)

σ(RAP)=σm*.

(13)

圖3 當(dāng)噪聲方差為20時(shí)的剩余自相關(guān)功率Fig.3 Residuals autocorrelation power (RAP) when the true variance is 20

圖4 對(duì)應(yīng)圖3不同的RAP由式(10)所得的DmFig.4 Difference of RAP, Dm in(10),for the RAPs in Fig.3

為了驗(yàn)證該方法的準(zhǔn)確性，對(duì)常規(guī)劑量的CT圖像加入均值為0、方差不同的Gaussian白噪聲進(jìn)行一系列的仿真實(shí)驗(yàn)，實(shí)驗(yàn)結(jié)果如表1所示，表示同一幅CT圖像加入的Gaussian白噪聲方差分別為5、10、15、20、25時(shí)的實(shí)驗(yàn)結(jié)果.

表1 同一圖像Gaussian白噪聲方差不同時(shí)實(shí)驗(yàn)結(jié)果對(duì)比Tab.1 Comparison of results in the same image with different White Gaussian noise variances

由表1的實(shí)驗(yàn)結(jié)果可以看出，對(duì)不同的Gaussian白噪聲方差，本文算法的估計(jì)精度均高于文獻(xiàn)[3]中的MAD方法.

3 本文算法

筆者提出的自適應(yīng)低劑量CT圖像質(zhì)量改善算法，其算法框圖如圖5所示，主要包括5個(gè)步驟.

圖5 本文算法框圖
Fig.5 The diagram of the proposed algorithm

步驟1：Anscombe變換.首先采用Anscombe變換將低劑量CT圖像中的量子噪聲轉(zhuǎn)變?yōu)榫哂薪瞥?shù)方差Gaussian噪聲[10].變換公式為：

(14)

步驟2：剪切波分解.將Anscombe變換后的圖像進(jìn)行4層剪切波分解，每層分解為8個(gè)方向，得到1個(gè)低頻圖像和32個(gè)高頻圖像.

步驟3：提取非噪聲系數(shù).由于隨著分解層數(shù)的增加，高頻子帶的信噪比逐漸降低，故利用提出的噪聲方差估計(jì)方法對(duì)信噪比較低的第3、4層高頻子帶估計(jì)其噪聲方差.然后基于噪聲方差估值，結(jié)合貝葉斯最大后驗(yàn)估計(jì)方法提取非噪聲高頻系數(shù).

步驟4：重構(gòu)圖像.基于步驟3獲得的非噪聲高頻系數(shù)和低頻系數(shù)進(jìn)行剪切波逆變換，獲得去噪圖像；對(duì)去噪后的圖像實(shí)行Anscombe逆變換，得到重構(gòu)圖像.公式為：

(15)

4 仿真實(shí)驗(yàn)

筆者從定量評(píng)價(jià)和視覺效果來驗(yàn)證提出算法的有效性.采用PSNR(peak signal to noise ratio)值和MSSIM(mean structure similarity)值[11]為定量評(píng)價(jià)指標(biāo):

(16)

(17)

4.1 定量評(píng)價(jià)

對(duì)取自鄭州大學(xué)第一附屬醫(yī)院的一系列含有病變的、大小為512×512的常規(guī)劑量的胸部CT圖像加入?yún)?shù)λ≤0.1的Poisson噪聲，并進(jìn)行仿真實(shí)驗(yàn)，采用PSNR和MSSIM來定量評(píng)價(jià)去噪效果,實(shí)驗(yàn)結(jié)果如表2和圖6所示.

表2表示對(duì)10幅常規(guī)劑量的CT圖像添加強(qiáng)度參數(shù)λ=0.01的Poisson噪聲，利用本文方法與文獻(xiàn)[3]方法對(duì)噪聲圖像進(jìn)行去噪的實(shí)驗(yàn)結(jié)果.

表2 當(dāng)Poisson強(qiáng)度λ=0.01時(shí)的PSNR和MSSIM比較Tab.2 Comparison of PSNR and MSSIM when the Poisson noise intensity is 0.01 dB

圖6表示對(duì)一幅CT圖像添加強(qiáng)度參數(shù)λ為0.005、0.008、0.010、0.050、0.080和0.100的Poisson噪聲，本文方法與文獻(xiàn)[7]方法對(duì)噪聲圖像去噪的實(shí)驗(yàn)結(jié)果.圖6(a)表示隨著噪聲量的變化，兩種算法的重構(gòu)圖像與噪聲圖像的PSNR折線圖，圖6(b)表示隨著噪聲圖像的MSSIM變化，兩種算法的重構(gòu)圖像的MSSIM折線圖.

圖6 不同Poisson噪聲強(qiáng)度時(shí)的PSNR和MSSIM對(duì)比圖Fig.6 Comparison of PSNR and MSSIM in image with different Poisson noise intensity

由表2可知，本文方法與文獻(xiàn)[3]方法相比，PSNR和MSSIM分別提高了52.2%和34.9%.從圖6可以看出，對(duì)參數(shù)λ≤0.1的Poisson噪聲，本文算法所得重構(gòu)圖像的PSNR和MSSIM均高于文獻(xiàn)[3]方法；本文方法在噪聲量較多的情況下，去噪效果尤為突出.

表2和圖6的實(shí)驗(yàn)結(jié)果表明，對(duì)參數(shù)λ≤0.1的Poisson噪聲，無論是同一噪聲強(qiáng)度下的不同圖像還是多種噪聲強(qiáng)度下的同幅圖像，與文獻(xiàn)[3]方法相比，本文算法所得重構(gòu)圖像的PSNR和MSSIM均有所提高.由此得出，本文算法具有較強(qiáng)的自適應(yīng)性和魯棒性.

4.2 視覺效果

本文視覺效果評(píng)價(jià)采用在日本千葉大學(xué)先端醫(yī)、工學(xué)研究中心，以體膜為掃描對(duì)象，進(jìn)行低劑量掃描所采集的降低不同劑量的多幅實(shí)際CT圖像.低劑量掃描的主要原理是利用放射線劑量與管電流的線性關(guān)系，當(dāng)保持其他參數(shù)不變時(shí)，通過降低管電流和管電壓來降低放射劑量.圖7(ai)(i=1,2)所表示圖像的掃描參數(shù)設(shè)置分別為管電壓120 kV，掃描層厚為0.5 mm，管電流分別為100 mA和150 mA；圖7(a3)表示圖像的取像參數(shù)設(shè)置分別為管電壓為120 kV，管電流為200 mA，掃描層厚為1 mm.

從圖7(ai)(i=1,2,3)可知，隨著管電流的減少，噪聲量越來越多.分別用本文方法與文獻(xiàn)法對(duì)這些圖像進(jìn)行仿真實(shí)驗(yàn)，實(shí)驗(yàn)結(jié)果分別由圖7(bi)和圖7(ci)(i=1,2,3)表示.圖7(bi)(i=1,2,3)為利用文獻(xiàn)[3]方法所獲得的去噪圖像，圖7 (ci)(i=1, 2,3)分別表示本文方法的重構(gòu)圖像.由圖7可以看出，本文算法不僅能有效去除噪聲，而且與文獻(xiàn)[3]方法相比能更好地保留圖像中的邊緣信息和紋理信息.

上述實(shí)驗(yàn)的定量分析和視覺效果表明，本文算法具有自適應(yīng)地去除參數(shù)λ≤0.1的Poisson噪聲的能力，而且該方法有效提升了重構(gòu)圖像的質(zhì)量.

圖7 低劑量CT圖像去噪效果對(duì)比Fig.7 The comparison of de-noising effects for low-dose CT

5 結(jié)論

筆者利用剪切波良好的稀疏特性，將Anscombe變換與改進(jìn)的自適應(yīng)噪聲方差估計(jì)算法相結(jié)合，提出了一種有效改善低劑量CT圖像質(zhì)量的方法.對(duì)含有λ≤0.1的Poisson噪聲的CT圖像進(jìn)行的定量評(píng)價(jià).結(jié)果表明，本文方法與文獻(xiàn)[3]方法相比，其PSNR平均提高52.2%，MSSIM提高34.9%.同時(shí)，利用本文方法和文獻(xiàn)[3]方法對(duì)多幅實(shí)際低劑量CT圖像進(jìn)行去噪的對(duì)比實(shí)驗(yàn)結(jié)果也表明，本文方法較文獻(xiàn)[3]方法可以更好地保留圖像細(xì)節(jié)信息.因此，本文算法具有自適應(yīng)能力強(qiáng)、重構(gòu)圖像質(zhì)量高的特點(diǎn)，具有一定的臨床應(yīng)用價(jià)值.