楊一鳴，劉 祎，桂志國(guó)

(中北大學(xué) 生物醫(yī)學(xué)成像與影像大數(shù)據(jù)重點(diǎn)實(shí)驗(yàn)室，山西 太原 030051)

0 引 言

近年來(lái)，計(jì)算機(jī)斷層掃描(Computed Tomography, CT)在肺癌檢查中廣泛應(yīng)用，主要原因是CT在檢測(cè)肺部小結(jié)節(jié)和早期肺癌中比胸部X-線(xiàn)檢查更準(zhǔn)確[1]. 數(shù)據(jù)顯示，一次常規(guī)CT檢查，受檢者將受到劑量約為1.5～20 mSv的X射線(xiàn)輻射，且CT掃描引發(fā)的X射線(xiàn)輻射在所有放射檢查產(chǎn)生的X射線(xiàn)輻射中比重較大，由于CT檢查可能導(dǎo)致受檢者余生患癌的風(fēng)險(xiǎn)加大，醫(yī)學(xué)檢測(cè)中已將重復(fù)CT掃描檢查列為醫(yī)源性輻射的主要輻射源[2]. 隨著診斷技術(shù)的不斷改進(jìn)，限制輻射劑量的低劑量CT(low-dose computed tomography, LDCT)已被用于肺癌高危人群的常規(guī)篩查，但低劑量CT又會(huì)造成圖像質(zhì)量下降. 如何在掃描劑量盡可能低的條件下，獲取與常規(guī)劑量掃描質(zhì)量相同或者更利于醫(yī)療診斷的CT圖像，自然就成為相關(guān)領(lǐng)域亟待研究的課題.

目前提高低劑量CT圖像質(zhì)量的方法可以分為三類(lèi)：投影數(shù)據(jù)校正法，迭代重建法和后處理方法. 低劑量CT圖像的后處理方法，是對(duì)已經(jīng)獲取的重建圖像進(jìn)行修正處理，這類(lèi)方法不依賴(lài)于投影數(shù)據(jù)，即直接處理濾波反投影算法重建的CT圖像，消除重建圖像中的條形偽影和噪聲. 后處理方法的難點(diǎn)是斑點(diǎn)噪聲和條形偽影分布于整個(gè)CT圖像上，不滿(mǎn)足特定的噪聲分布模型. 這種方法一般以標(biāo)準(zhǔn)劑量CT(Standard Dose Computed Tomography, SDCT)圖像為參考標(biāo)準(zhǔn)，一方面要保證處理后的結(jié)果不能丟失原來(lái)低劑量CT圖像中較小的組織結(jié)構(gòu)，另一方面不能在處理的過(guò)程中引入新的噪聲和偽影. 由于低劑量CT圖像后處理效果的好壞直接影響醫(yī)生診斷的準(zhǔn)確性，國(guó)內(nèi)外很多學(xué)者致力于對(duì)低劑量CT圖像偽影和噪聲抑制方面的研究. 最初，Rudin等[3]提出了總變分(Total Variation, TV)模型，TV去噪模型能夠較好地保護(hù)圖像邊緣，但是該模型僅沿圖像邊緣的切線(xiàn)方向擴(kuò)散，導(dǎo)致處理后的圖像產(chǎn)生階梯效應(yīng)，后來(lái)，Zhu Yining等[4]提出一種改進(jìn)的TV模型用于低劑量CT降噪，取得了不錯(cuò)的降噪效果； Zamyatin等[5]提出了一種自適應(yīng)多尺度總變分低劑量CT圖像濾波方法，在對(duì)低劑量CT圖像偽影抑制的同時(shí)，較好地保護(hù)了腹部CT圖像中組織結(jié)構(gòu)的邊緣； Liu Yi等[6]針對(duì)稀疏角度投影提出了基于中值先驗(yàn)約束的TV重建算法. 2005年，Buades等[7]首次提出了非局部均值(Non-Local Means, NLM)算法，該算法利用圖像具有重復(fù)結(jié)構(gòu)性質(zhì)，在全局范圍內(nèi)搜索與像素點(diǎn)所在圖像塊相應(yīng)的相似塊，再對(duì)其進(jìn)行加權(quán)平均來(lái)實(shí)現(xiàn)去噪. Chen Yang等[8]提出了利用結(jié)構(gòu)相似信息的大尺度非局部均值濾波來(lái)處理腹部低劑量CT圖像； Ha等[9]提出了通過(guò)小圖像塊數(shù)據(jù)庫(kù)與局部相似性匹配相結(jié)合的方法，對(duì)低劑量CT圖像后處理； 此外，Chen Yang等[10]提出了基于快速字典學(xué)習(xí)方法的腹部腫瘤低劑量CT圖像偽影抑制算法.

字典學(xué)習(xí)和稀疏表示理論在圖像處理領(lǐng)域的應(yīng)用十分廣泛，用于圖像的去噪、 修復(fù)、 超分辨率重建[11]等方面. Aharon等[12]最早提出了K-奇異值分解(K-Singular Value Decomposition, K-SVD)算法，該算法利用逐列更新字典的方式，多次迭代建立了一個(gè)全局最優(yōu)的過(guò)完備字典，通過(guò)字典的冗余性對(duì)圖像進(jìn)行稀疏表示； Chen Yang等[13]提出了一種基于區(qū)別性字典的條形偽影抑制方法，該方法首先使用小波變換分解低劑量CT圖像，然后對(duì)水平、 垂直和對(duì)角方向的高頻圖像分別進(jìn)行對(duì)應(yīng)的區(qū)別性字典表達(dá)并去除偽影，最后再用傳統(tǒng)的字典學(xué)習(xí)方法去除殘留的噪聲，取得了很好的效果. 受這種方法啟發(fā)，本文提出了一種基于字典學(xué)習(xí)的偽影抑制算法，首先利用平穩(wěn)小波變換(Stationary Wavelet Transform, SWT)對(duì)低劑量CT圖像進(jìn)行單層分解，并對(duì)三個(gè)方向的高頻圖像利用在線(xiàn)字典學(xué)習(xí)的方法訓(xùn)練字典，然后利用等效視數(shù)(Equivalent Number of Looks, ENL)對(duì)字典進(jìn)行分區(qū)得到偽影字典和特征字典，并只對(duì)特征原子進(jìn)行稀疏編碼，經(jīng)過(guò)小波逆變換(Inverse Stationary Wavelet Transform, ISWT)后，再采用雙邊濾波器(Bilateral Filter, BF)對(duì)處理后的CT圖像進(jìn)行分解并訓(xùn)練高頻字典，最后通過(guò)等效視數(shù)摒棄偽影字典來(lái)去除高頻圖像殘留的偽影和噪聲，從而達(dá)到抑制條形偽影、 改善低劑量CT圖像質(zhì)量的目的. 將本文算法應(yīng)用于實(shí)際的低劑量CT圖像，與標(biāo)準(zhǔn)劑量CT圖像相比較，并與TV降噪算法、 K-SVD算法和三維塊匹配濾波(Block-Matching and 3-D Filtering, BM3D)算法[14]對(duì)比，本文算法的處理結(jié)果優(yōu)于其他算法.

1 字典學(xué)習(xí)理論

1.1 稀疏表示

(1)

(2)

式中：‖α‖0為向量α的l0范數(shù)，表示α中的非零元素個(gè)數(shù). 式(2)中允許y與Dα之間存在一定的誤差ε，但是式(2)中的l0范數(shù)最小化問(wèn)題在求解時(shí)是一個(gè)NP難題. 對(duì)于這種問(wèn)題的求解方法，可以利用正交匹配追蹤算法(Orthogonal Matching Pursuit，OMP)求近似解.

1.2 正交匹配追蹤算法

OMP算法[15]是在過(guò)完備字典中選擇較少的列向量，形成一個(gè)對(duì)觀(guān)測(cè)值的最優(yōu)近似表示. 基于稀疏度限制L的稀疏表示問(wèn)題可以描述為

(3)

1.3 在線(xiàn)字典學(xué)習(xí)

1) 圖像塊處理：假設(shè)Pi是從一幅圖像I中得到的第i個(gè)長(zhǎng)度為p的圖像塊，將Pi表示成向量形式. 因此，訓(xùn)練集中第i個(gè)圖像塊可以分別表示為Ii=PiI.

2) 初始化訓(xùn)練參數(shù)：設(shè)定字典在訓(xùn)練階段的正則化參數(shù)為λ，圖像塊集為Pn={P1,P2,…,Pn}，以Pn為輸入信號(hào).

3) 計(jì)算稀疏系數(shù)：由步驟(2)可知，迭代開(kāi)始時(shí)D=D0，用最小角回歸迭代方法求解，可使式(4) 的計(jì)算值最小.

(4)

式中：t表示訓(xùn)練字典的迭代次數(shù)，t={1,2,…,T}，根據(jù)輸入的圖像塊向量Ii更新稀疏系數(shù)θt.

Dt

(5)

為了保證字典中的基向量值在一定的范圍內(nèi)，對(duì)求得的D加一個(gè)約束條件D∈c，歸一化處理結(jié)果為

(6)

5) 最后生成字典：達(dá)到迭代次數(shù)T后，即可獲得過(guò)完備字典DT.

2 本文算法

2.1 平穩(wěn)小波變換

在經(jīng)典的小波變換算法中，圖像經(jīng)過(guò)高通濾波和低通濾波后都要經(jīng)過(guò)下采樣，使分解后的各子帶變?yōu)樵瓐D像大小的一半，但這樣容易造成圖像邊緣的小波系數(shù)部分信息丟失，給重構(gòu)圖像帶來(lái)不穩(wěn)定性. 而平穩(wěn)小波變換是對(duì)經(jīng)典小波變換的改進(jìn)，其最大的特點(diǎn)是冗余性和平移不變性. 本文算法中首先采用平穩(wěn)小波變換對(duì)低劑量CT圖像進(jìn)行分解，通過(guò)平穩(wěn)小波變換將LDCT圖像分成低頻和高頻兩部分，其中高頻又包含了三個(gè)部分：水平高頻部分，垂直高頻部分和對(duì)角高頻部分. 本文算法中使用哈爾(Haar)小波基進(jìn)行分解，因其速度快且可以避免Gibbs振蕩現(xiàn)象.

2.2 基于等效視數(shù)的字典分區(qū)方法

在對(duì)質(zhì)量退化的低劑量CT圖像I進(jìn)行平穩(wěn)小波變換(SWT)之后，將其分解成低頻分量Ic a和三個(gè)高頻分量：水平分量Ic hd、 垂直分量Ic vd和對(duì)角分量Ic dd； 然后把高頻部分分別運(yùn)用在線(xiàn)字典學(xué)習(xí)方法訓(xùn)練出三個(gè)方向上的字典Dc hd，Dcvd和Dc dd； 圖 1 為構(gòu)成字典的特征原子及偽影原子，均由水平分量字典隨機(jī)抽取得到，相比特征原子圖像塊中的圖像信息，偽影原子圖像塊中的圖像信息變化更顯混亂.

圖 1 構(gòu)成字典的原子Fig.1 Atoms from dictionary

(7)

(8)

(9)

圖 2 三個(gè)方向高頻分量的字典分區(qū)示意圖Fig.2 Illustration of the dictionaries for the high frequency bands with different orientations

2.3 殘留偽影的抑制

圖 3 算法各步驟處理結(jié)果局部放大圖Fig.3 Local amplified images of processing results by each step of the algorithm

本文算法的流程示意圖見(jiàn)圖 4，算法的具體步驟如下：

圖 4 本文算法的流程示意圖Fig.4 Flow schematic of the algorithm

Step 1：對(duì)低劑量CT圖像I進(jìn)行平穩(wěn)小波變換，分解成低頻分量Ic a與三個(gè)高頻分量Ic hd,Icvd和Ic dd；

Step 2：設(shè)定字典訓(xùn)練階段用到的正則化參數(shù)λ和圖像塊大小q，輸入圖像的圖像塊集Pn，初始字典D0，迭代次數(shù)T，利用式(4)進(jìn)行字典訓(xùn)練，利用式(5)迭代更新字典；

Step 3：利用式(7)計(jì)算圖像塊的等效視數(shù)，將字典分區(qū)成偽影原子字典和特征原子字典；

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

為了驗(yàn)證本文提出算法的有效性，采用實(shí)際的臨床CT圖像進(jìn)行實(shí)驗(yàn)，CT圖像由一臺(tái)西門(mén)子SOMATOM Sensation 16排CT得到，圖像大小為512×512. 其中LDCT圖像由40 mAs電流掃描，SDCT圖像由160 mAs電流掃描，重建算法為FBP. 選取如圖5所示的兩張LDCT圖像進(jìn)行測(cè)試，偽影抑制的結(jié)果與SDCT圖像做對(duì)比，兩組圖都用腹窗顯示(窗寬350 HU，窗位50 HU). LDCT圖像選用B70f卷積核，SDCT圖像則選用臨床上腹部掃描和胸部掃描常用的B30f卷積核. 實(shí)現(xiàn)算法的工具為MATLAB R2014a，計(jì)算機(jī)的配置為CPU：Intel CoreTMi7-4771 QM CPU@3.50 GHz，內(nèi)存8 192 Mb RAM. 本文算法中涉及到7個(gè)參數(shù). 在線(xiàn)字典學(xué)習(xí)方法有4個(gè)參數(shù)，正規(guī)化參數(shù)λ取值為0.1，圖像塊的大小為16×16，字典訓(xùn)練塊數(shù)量即為：(512-15)×(512-15)=247 009，字典中原子的個(gè)數(shù)為1 500，迭代的次數(shù)為90. 圖像塊的大小取為16×16，一方面有利于字典學(xué)習(xí)時(shí)提取圖像的細(xì)節(jié)特征，另一方面有利于運(yùn)用等效視數(shù)的方法區(qū)分偽影原子和特征原子. 運(yùn)用分區(qū)的子字典稀疏表示時(shí)，所用OMP方法的稀疏限制水平參數(shù)L取為3. 雙邊濾波器中有3個(gè)參數(shù)，一個(gè)是濾波器的窗口的大小，通常設(shè)置為10×10，另外兩個(gè)參數(shù)是空間標(biāo)準(zhǔn)差σd和密度標(biāo)準(zhǔn)差σr，在估計(jì)了幾個(gè)可能的值后，發(fā)現(xiàn)σd=6，σr=0.2時(shí)能去除低劑量CT圖像中的大部分的條形偽影.

圖 5 測(cè)試圖像Fig.5 Test images

圖 6(c)～(f)給出了采用TV算法，K-SVD算法，BM3D算法處理結(jié)果以及本文算法對(duì)測(cè)試圖a的結(jié)果對(duì)比，圖6(a)和(b)分別為腹部LDCT圖像和腹部SDCT圖像. 圖 7 給出了圖 6 中各圖對(duì)應(yīng)的局部放大結(jié)果. 對(duì)比圖7中(a1)和(b1)可以看出，相比于SDCT圖像，LDCT圖像中的偽影比較嚴(yán)重，解剖結(jié)構(gòu)的細(xì)節(jié)分辨率低； 對(duì)比圖 7 中(c1)和(b1)可以看出，TV算法雖然有效地去除了條形偽影，但是由于圖像過(guò)度平滑造成細(xì)節(jié)信息損失； 對(duì)比圖 7 中(d1)和(b1)可以看出，K-SVD算法去除了大量的噪聲和偽影，但是圖像有部分邊緣變得模糊； 觀(guān)察圖 7 中(e1)可以看出，BM3D算法有效抑制了圖像中的噪聲和偽影，但是對(duì)比圖 7 中(b1) 的部分細(xì)節(jié)發(fā)現(xiàn)BM3D算法處理后的圖像過(guò)度平滑造成了邊緣和細(xì)節(jié)信息的丟失. 從對(duì)比結(jié)果來(lái)看，本文的算法能有效地抑制偽影，在圖像的邊緣和細(xì)節(jié)保持方面優(yōu)于其他算法，結(jié)果圖與SDCT圖像最為接近.

圖 6 測(cè)試圖a的各種算法處理后結(jié)果圖Fig.6 Various algorithms processed result images of test image a

圖 7 中的(a2)～(f2) 給出了圖 6 中各圖對(duì)應(yīng)的肝臟腫瘤部分局部放大結(jié)果，圖 7 中(a2)中箭頭所指為肝臟腫瘤部分. 對(duì)比圖 7 中(a2)和(b2)可以看出，相比SDCT圖像中的腫瘤部分，LDCT圖像難以找到準(zhǔn)確的腫瘤部分，噪聲影響比較嚴(yán)重； 對(duì)比圖 7 中(c2), (d2), (e2)和(f2)可以看出，圖 7 中(f2)與SDCT圖像的腫瘤部分最接近，效果最好.

圖 7 測(cè)試圖a結(jié)果圖對(duì)應(yīng)的局部放大結(jié)果Fig.7 Partial magnification of result images of the test image a

圖 8 給出了測(cè)試圖b的各種算法處理結(jié)果圖，圖8(a)～(f)分別是LDCT圖像、 SDCT圖像、 TV算法處理結(jié)果、 K-SVD算法處理結(jié)果、 BM3D算法處理結(jié)果和本文算法處理結(jié)果圖.

圖 9 給出了圖 8 中各圖對(duì)應(yīng)的局部放大結(jié)果. 對(duì)比圖9(c)和(b)，TV算法會(huì)由于圖像過(guò)度平滑造成細(xì)節(jié)信息損失； 對(duì)比圖9(d)和(b)可以看出K-SVD算法去除了大量的噪聲和偽影，但是圖像還是存在部分邊緣模糊的情況； 由圖9(e)可以看出，BM3D算法有效抑制了圖像中的噪聲和偽影，但是與圖9(d)對(duì)比會(huì)發(fā)現(xiàn)BM3D算法處理后的圖像細(xì)節(jié)部分比較好. 而本文的算法在抑制噪聲和偽影的同時(shí)有效地保持了組織結(jié)構(gòu)分辨度，展現(xiàn)出了真實(shí)SDCT圖像的紋理部分. 處理后的測(cè)試圖b 中組織部分保留的細(xì)節(jié)紋理較多且清晰，因此從視覺(jué)效果來(lái)看，本文算法的處理結(jié)果相比其他算法最好.

圖 8 測(cè)試圖b的各種算法處理后結(jié)果圖 Fig.8 Various algorithms processed result images of test image b

圖 9 測(cè)試圖b結(jié)果圖對(duì)應(yīng)的局部放大結(jié)果Fig.9 Partial magnification of result images of the test image b

為了對(duì)算法結(jié)果進(jìn)行客觀(guān)、 定量的評(píng)價(jià)，本文采用較常用的峰值信噪比、 均方根誤差和結(jié)構(gòu)相似性來(lái)評(píng)價(jià)處理后CT圖像的質(zhì)量. 這些參數(shù)的定義如下：

1) 峰值信噪比(Peak Signal-to-Noise Ratio, PSNR)

式中：M×N表示圖像的大?。籿是SDCT圖像；u是偽影抑制后的圖像. PSNR是圖像信噪比變化情況的統(tǒng)計(jì)平均，其值越大，圖像失真越少，這里取Peak=255.

2) 均方根誤差(Root Mean Square Error, RMSE)

(12)

均方根誤差是用來(lái)衡量u與u0之間的偏差，其值越小，圖像u失真越小.

3) 結(jié)構(gòu)相似性(Structural Similarity Index Measurement, SSIM)

SSIM(v,u)=[L(v,u)]α×[C(v,u)]β×

[S(v,u)]γ，

(13)

(14)

(15)

(16)

表 1 和表 2 分別給出了不同算法對(duì)測(cè)試圖a和測(cè)試圖b處理后的質(zhì)量評(píng)估參數(shù)，表 1 和表 2 由圖 10 分別截取SDCT圖像的感興趣區(qū)域(Region of Interest, ROI)計(jì)算得到，每張測(cè)試圖像分別截取兩個(gè)感興趣區(qū)域.

表 1 測(cè)試圖a的質(zhì)量評(píng)估參數(shù)

表 2 測(cè)試圖b的質(zhì)量評(píng)估參數(shù)

由表 1 和表 2 分析可知：本文提出算法與其他算法相比，結(jié)構(gòu)相似性SSIM和峰值信噪比PSNR都較大，均方根誤差RMSE較小. 說(shuō)明本算法處理后的圖像質(zhì)量較好，能夠有效地抑制偽影，并保留圖像的邊緣等重要細(xì)節(jié)特征. 因此，在視覺(jué)效果和定量評(píng)價(jià)方面，所得結(jié)果都表明本文提出的算法對(duì)低劑量CT圖像的偽影抑制是有效的.

表 3 給出了不同算法對(duì)測(cè)試圖a和測(cè)試圖b處理后的運(yùn)行時(shí)間，可以看出，本文所提算法的計(jì)算速度優(yōu)于K-SVD算法，但是相比于BM3D算法，運(yùn)行時(shí)間還是相對(duì)較長(zhǎng)，對(duì)算法加速方面而言，可作為下一步研究的方向.

圖 10 選取SDCT圖像的感興趣區(qū)域Fig.10 Select the region of interest of the SDCT image

4 結(jié) 論

本文提出了一種基于字典學(xué)習(xí)與等效視數(shù)的偽影抑制算法，利用FBP算法獲得低劑量CT重建圖像，再對(duì)CT圖像進(jìn)行降噪處理. 本文算法在處理低劑量CT圖像的過(guò)程中，每一步驟得到的高頻分量都運(yùn)用字典學(xué)習(xí)方法訓(xùn)練得到字典，再運(yùn)用等效視數(shù)判別法對(duì)字典進(jìn)行分區(qū)得到偽影字典和特征字典，并只對(duì)特征原子進(jìn)行稀疏編碼，從而達(dá)到去除條形偽影的目的. 實(shí)驗(yàn)結(jié)果表明，本文算法可以在抑制條形偽影的同時(shí)，保留較多的細(xì)節(jié)信息，并且更加接近原始SDCT圖像. 從視覺(jué)效果和質(zhì)量評(píng)價(jià)可以看出，本文算法的偽影抑制效果優(yōu)于TV算法、 K-SVD算法和BM3D算法. 但是，本文所提算法中一些參數(shù)的選擇，仍需憑借經(jīng)驗(yàn)調(diào)整，這也是后續(xù)研究中需要解決的一個(gè)問(wèn)題.