焦楓媛，桂志國，劉 祎，韓 意

(1. 中北大學(xué) 生物醫(yī)學(xué)成像與影像大數(shù)據(jù)山西省重點(diǎn)實(shí)驗(yàn)室，山西 太原 030051；2. 中北大學(xué) 信息與通信工程學(xué)院，山西 太原 030051；3. 中國人民解放軍63961部隊(duì)，北京 100012)

X射線計(jì)算機(jī)斷層成像(X-ray Computed Tomography, CT)技術(shù)在醫(yī)學(xué)領(lǐng)域獲得了廣泛應(yīng)用，在CT掃描過程中不可避免地會(huì)產(chǎn)生大量輻射劑量，對(duì)人體造成傷害[1]. 于是，低劑量CT技術(shù)應(yīng)運(yùn)而生，如何改善由于劑量降低帶來的圖像質(zhì)量下降問題也隨之成為當(dāng)前研究的一個(gè)熱點(diǎn)，目前的處理方法主要有3類： 投影域降噪[2,3]、統(tǒng)計(jì)迭代重建算法[4,5]及后處理算法[6,7]. 其中，后處理算法因其不依賴投影數(shù)據(jù)且易于操作和推廣等優(yōu)勢(shì)吸引了國內(nèi)外學(xué)者的關(guān)注.

國內(nèi)外眾多學(xué)者將經(jīng)典的圖像處理方法改進(jìn)后應(yīng)用于低劑量CT圖像處理，取得了較好的成果. Chen等[8]在肩部和腹部低劑量CT處理中，應(yīng)用了大尺度鄰域的加權(quán)平均算法，有效抑制了噪聲和偽影； Kang等[9]在低劑量CT冠脈數(shù)據(jù)處理中應(yīng)用了自適應(yīng)的三維塊匹配算法，既有效抑制了噪聲，又改善了左心室功能的評(píng)估； Chen等[10]提出的三維塊匹配算法對(duì)于提高圖像的對(duì)比噪聲比效果顯著； Wu等[11]提出的級(jí)聯(lián)卷積神經(jīng)網(wǎng)絡(luò)算法獲得了更高質(zhì)量的低劑量CT處理圖； Wolterink等[12]將生成式對(duì)抗網(wǎng)絡(luò)應(yīng)用于低劑量CT圖像去噪； Shan等[13]在低劑量CT圖像處理中將二維訓(xùn)練網(wǎng)絡(luò)的遷移學(xué)習(xí)與基于合同路徑的三維卷積編碼-解碼網(wǎng)絡(luò)相結(jié)合，有效保持了圖像的精細(xì)結(jié)構(gòu).

基于偏微分方程(Partial Differential Equation, PDE)的方法是圖像處理中一類重要方法. 1990年，Perona等[14]提出了著名的各項(xiàng)異性擴(kuò)散模型，即Perona-Malik (PM)模型，比較以往的各項(xiàng)同性擴(kuò)散模型，PM模型能在濾除噪聲的同時(shí)，更好地保持圖像邊緣. Rudin等[15]于1992年提出了全變分(Total Variatuon, TV)模型，在保持邊緣方面有更好的表現(xiàn). 為降低二階PDE所產(chǎn)生的塊狀效應(yīng)，You等[16]提出了一類四階PDE，這類算法處理結(jié)果圖中常常含有斑點(diǎn). 為克服以上PDE模型帶來的問題，Bai等[17]提出了分?jǐn)?shù)階PM(Fractional-order PM, FPM)模型； Zhang等[18]提出了分?jǐn)?shù)階TV(Fractional-order TV, FTV)模型； Wang等[19]將FPM與FTV模型加權(quán)組合并加以改進(jìn)提出了分?jǐn)?shù)階PMTV模型(Fractional-order PMTV, FPMTV)，這類分?jǐn)?shù)階PDE模型既能克服階梯效應(yīng)、斑點(diǎn)效應(yīng)，又能在一定程度上保持圖像的邊緣和細(xì)節(jié). Hossein等[20]將殘差局部冪應(yīng)用于圖像處理，通過殘差圖像的利用，有效提高了PDE模型降噪過程中細(xì)節(jié)和紋理的保持.

由于低劑量CT圖像中存在大量組織結(jié)構(gòu)和細(xì)節(jié)，在去除噪聲的同時(shí)，最大限度地保持圖像細(xì)節(jié)對(duì)于醫(yī)療診斷具有重要意義，受文獻(xiàn)[20]啟發(fā)，本文將局部殘差冪作為紋理檢測(cè)算子，引入FPM模型的擴(kuò)散系數(shù)中，使其與梯度一起自適應(yīng)調(diào)節(jié)模型擴(kuò)散系數(shù)，使得獲得的低劑量CT處理結(jié)果圖在去除噪聲的同時(shí)更好地保持了圖像紋理和細(xì)節(jié)信息. 通過頭部體模和實(shí)際臨床數(shù)據(jù)實(shí)驗(yàn)驗(yàn)證了本文算法的有效性.

1 分?jǐn)?shù)階微積分

Grümwald-Letniklv (G-L)和Riemann- Liouville (R-L) 定義是分?jǐn)?shù)階微積分各類定義中應(yīng)用最廣泛的，本文將G-L定義應(yīng)用于低劑量CT圖像處理中，信號(hào)f(x)的α階微分的G-L定義如下[18]

(1)

2 基于殘差局部冪的分?jǐn)?shù)階各項(xiàng)異性擴(kuò)散算法

2.1 整數(shù)階各項(xiàng)異性擴(kuò)散算法

整數(shù)階各項(xiàng)異性擴(kuò)散算法，即PM模型的表達(dá)式為

(2)

式中，擴(kuò)散系數(shù)函數(shù)c可表示為

c(|?f|)=1/[1+|?f|2/k2] ,

(3)

式中：k為梯度閾值. PM模型雖然取得了很好的去噪效果，但是降噪后的圖像經(jīng)常會(huì)出現(xiàn)塊狀效應(yīng).

2.2 基于殘差局部冪的分?jǐn)?shù)階各項(xiàng)異性擴(kuò)散算法

為了克服階梯效應(yīng)，Bai等[17]在2007年提出了分?jǐn)?shù)階各項(xiàng)異性擴(kuò)散方程，即FPM模型，其能量函數(shù)可表示為

(4)

式中：Ω為圖像支撐，且p(|Dαf|)≥0是一個(gè)增函數(shù)，其與擴(kuò)散系數(shù)關(guān)系可表示為

(5)

為了解決上述問題，更好地輔助臨床診斷，本文將紋理檢測(cè)算子添加到分?jǐn)?shù)階各項(xiàng)異性擴(kuò)散方程的擴(kuò)散系數(shù)中. 由文獻(xiàn)[20]可知，殘差局部冪可作為一種紋理檢測(cè)算子，首先定義殘差圖像

fR=f0-f=fs+fn,

(6)

式中：fs表示紋理和細(xì)節(jié)成分；fn表示噪聲成分，由于兩者具有不相關(guān)性，因此

PR=Ps+Pn,

(7)

式中：PR、Ps和Pn分別為fR、fs和fn的局部冪，將PR做如下歸一化處理后可得

(8)

由此，分?jǐn)?shù)階各項(xiàng)異性擴(kuò)散模型中擴(kuò)散系數(shù)可寫為

(9)

式中：k1為常數(shù).

根據(jù)文獻(xiàn)[17]中分?jǐn)?shù)階微分的離散化方法，將(4)式中的分?jǐn)?shù)階微分替換為卷積積分后可表示為

(10)

式中：

根據(jù)梯度下降法，可以得到

(11)

將(vα*f)x和(vα*f)y離散化獲得

(12)

記

式中：ε為很小的整數(shù). 基于局部殘差冪的分?jǐn)?shù)階各項(xiàng)異性擴(kuò)散方程的離散化形式可表示為

[Dx(i+k,j)+Dy(i,j+k)].

(14)

2.3 算法描述

1) 通過FPM模型將噪聲和相關(guān)紋理進(jìn)行分離. 將低劑量CT圖像fs與FPM處理結(jié)果圖f0的差作為殘差圖，當(dāng)所得殘差的方差等于低劑量CT圖像噪聲方差的β倍時(shí)，F(xiàn)PM算法停止迭代.

2) 計(jì)算殘差局部冪. 將步驟1)停止迭代后獲得的殘差圖通過高斯濾波濾除部分噪聲后，再應(yīng)用式(8)獲得歸一化的殘差局部冪.

3) 將基于殘差局部冪的FPM模型應(yīng)用于低劑量CT圖像. 將步驟2)獲得的殘差局部冪加入分?jǐn)?shù)階各項(xiàng)異性擴(kuò)散模型擴(kuò)散系數(shù)中，根據(jù)式(9)～式(14)進(jìn)行計(jì)算，獲得最終處理結(jié)果圖.

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

為了驗(yàn)證所提算法的有效性，本文對(duì)Shepp-Logan頭部體模和實(shí)際臨床數(shù)據(jù)進(jìn)行仿真實(shí)驗(yàn). 頭部體模大小為256 mm×256 mm，實(shí)際臨床數(shù)據(jù)為DICOM醫(yī)學(xué)圖像庫中大小為512 mm×512 mm 的腹部高劑量CT與低劑量CT圖. 本文與FPM算法[17]、改進(jìn)的非局部均值(improved Non-Local Means, INLM)算法[21]和FPMTV算法[19]進(jìn)行了比較. 實(shí)驗(yàn)環(huán)境如下： 操作系統(tǒng)是64位Windows 10，處理器是Intel Core(TM)i7-8565U CPU @ 1.80 GHz 1.99 GHz，內(nèi)存為32 GB. 比較算法的參數(shù)根據(jù)相應(yīng)參考文獻(xiàn)選擇，本文算法參數(shù)設(shè)置為β=1.4，α=1.03，Δt=0.25，迭代次數(shù)在頭部體模實(shí)驗(yàn)中設(shè)為2 200次，在腹部數(shù)據(jù)中設(shè)為55次.

圖 1 頭部體模各算法處理結(jié)果圖及局部 感興趣區(qū)域放大圖對(duì)比Fig.1 Comparison of the denoised images and their zoom-in ROIs processed by a variety of algorithms on Shepp-Logan phantom

圖 1 給出了各種算法應(yīng)用于頭部體模所獲得的處理結(jié)果圖及局部感興趣區(qū)域放大圖對(duì)比. 從圖 1 可以看出，F(xiàn)PM和FPMTV算法處理結(jié)果圖不同程度地模糊了圖像的邊緣和細(xì)節(jié)，INLM算法處理結(jié)果圖在噪聲較多的區(qū)域仍然有明顯的偽影留存，本文算法在濾除噪聲的同時(shí)，很好地保持了圖像的邊緣和細(xì)節(jié).

圖 2 給出了各種算法處理腹部低劑量CT圖的結(jié)果對(duì)比，在實(shí)際臨床數(shù)據(jù)中，高、低劑量圖中均含有噪聲點(diǎn)，但低劑量圖中含有更多斑點(diǎn)噪聲，不利于腹部組織結(jié)構(gòu)的顯示. 通過比較各算法處理結(jié)果圖可以看出，本文算法在去除噪聲點(diǎn)的同時(shí)更好地保持了基本的組織結(jié)構(gòu).

圖 2 腹部低劑量CT圖各算法處理結(jié)果對(duì)比Fig.2 Comparison of the denoised images processed bya variety of algorithms onabdominal low-dose CT image

圖 3 給出了腹部低劑量CT圖與不同算法處理結(jié)果圖的差值圖. 從圖3(a)和圖3(c)可以看出，F(xiàn)PM和FPMTV算法在噪聲被濾除的同時(shí)，部分邊緣和結(jié)構(gòu)也被濾除，從圖3(b)和圖3(d)可以看出，INLM和本文算法在保持邊緣和細(xì)節(jié)方面明顯優(yōu)于前兩個(gè)算法，進(jìn)行細(xì)微比較則可以發(fā)現(xiàn)本文算法殘差圖中幾乎看不到殘留的邊緣和細(xì)節(jié)信息.

為了對(duì)本文算法的有效性進(jìn)行定量驗(yàn)證，在Shepp-Logan頭部體模實(shí)驗(yàn)中，采用峰值信噪比(Peak Signal to Noise Ratio，PSNR)和結(jié)構(gòu)相似度(Structural Similarity，SSIM)對(duì)各算法進(jìn)行定量描述. 在實(shí)際臨床數(shù)據(jù)實(shí)驗(yàn)中，對(duì)于病人呼吸和運(yùn)動(dòng)等因素帶來的高、低劑量CT圖難以完全匹配問題，采用標(biāo)準(zhǔn)差(Standard Deviation, STD)進(jìn)行客觀評(píng)價(jià). 各評(píng)價(jià)指標(biāo)定義為

(7)

(8)

式中：

σuoriginalu=Cov{uoriginal,u}=

(9)

圖 3 腹部低劑量CT圖與不同算法處理結(jié)果圖的差值圖對(duì)比Fig.3 Comparison of the image subtractions between the abdominal low-dose CT image and the denoised images by a variety of algorithms on abdominal low-dose CT image

表 1 為FPM、INLM、FPMTV以及本文算法分別應(yīng)用于頭部體模而獲得的處理結(jié)果圖相對(duì)于原始圖像的PSNR和SSIM值的比較. 由表1很容易看出，本文算法的PSNR值最高，說明本文算法處理結(jié)果圖中含有最少的噪聲； 本文算法的SSIM值最高，說明本文算法處理結(jié)果圖與原圖最接近.

表 1 頭部體模各種算法的客觀評(píng)價(jià)Tab.1 The objective evaluation of a variety of algorithms on Shepp-Logan phantom

圖 4 為低、高劑量腹部CT圖像以及不同算法處理結(jié)果圖的4個(gè)ROIs的STD值柱狀圖，圖 4 左上角的腹部圖用方框標(biāo)注了4個(gè)ROIs. 通過柱狀圖可以直觀地看出，在4個(gè)ROIs中，本文提出的FPMLP算法對(duì)應(yīng)的STD值均最接近高劑量CT圖，這表明本文算法處理結(jié)果圖與高劑量CT圖最接近. 結(jié)合圖 1～圖 4 和表 1 可知，在視覺效果和定量評(píng)價(jià)方面，本文算法都優(yōu)于其他算法，既有效降低了低劑量CT圖像的噪聲和偽影，又較好地保持了CT圖像的紋理、細(xì)節(jié)等結(jié)構(gòu)特征.

圖 4 低、高劑量CT圖像以及各算法處理結(jié)果圖在腹部圖方框標(biāo)注的4個(gè)ROIs的STD值柱狀圖

4 結(jié) 論

本文提出了一種基于殘差局部冪的各項(xiàng)異性擴(kuò)散低劑量CT降噪算法. 該方法將殘差局部冪作為紋理檢測(cè)算子加入分?jǐn)?shù)階各項(xiàng)異性擴(kuò)散模型的擴(kuò)散系數(shù)中，充分利用了殘差圖像中的紋理、細(xì)節(jié)和邊緣信息，有效增強(qiáng)了分?jǐn)?shù)階各項(xiàng)異性擴(kuò)散算法的降噪性能. 基于Shepp-Logan和實(shí)際腹部醫(yī)療數(shù)據(jù)的仿真實(shí)驗(yàn)結(jié)果表明，本文算法在抑制噪聲和保持圖像紋理、細(xì)節(jié)和邊緣方面均優(yōu)于傳統(tǒng)的低劑量CT降噪算法.