基于殘差局部冪的分數(shù)階各項異性擴散低劑量CT降噪算法

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

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

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

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

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

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

1 分數(shù)階微積分

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

(1)

2 基于殘差局部冪的分數(shù)階各項異性擴散算法

2.1 整數(shù)階各項異性擴散算法

整數(shù)階各項異性擴散算法，即PM模型的表達式為

(2)

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

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

(3)

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

2.2 基于殘差局部冪的分數(shù)階各項異性擴散算法

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

(4)

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

(5)

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

fR=f0-f=fs+fn,

(6)

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

PR=Ps+Pn,

(7)

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

(8)

由此，分數(shù)階各項異性擴散模型中擴散系數(shù)可寫為

(9)

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

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

(10)

式中：

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

(11)

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

(12)

記

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

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

(14)

2.3 算法描述

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

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

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

3 實驗結(jié)果及分析

為了驗證所提算法的有效性，本文對Shepp-Logan頭部體模和實際臨床數(shù)據(jù)進行仿真實驗. 頭部體模大小為256 mm×256 mm，實際臨床數(shù)據(jù)為DICOM醫(yī)學圖像庫中大小為512 mm×512 mm 的腹部高劑量CT與低劑量CT圖. 本文與FPM算法[17]、改進的非局部均值(improved Non-Local Means, INLM)算法[21]和FPMTV算法[19]進行了比較. 實驗環(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ù)相應參考文獻選擇，本文算法參數(shù)設置為β=1.4，α=1.03，Δt=0.25，迭代次數(shù)在頭部體模實驗中設為2 200次，在腹部數(shù)據(jù)中設為55次.

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

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

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

圖 2 腹部低劑量CT圖各算法處理結(jié)果對比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算法在噪聲被濾除的同時，部分邊緣和結(jié)構(gòu)也被濾除，從圖3(b)和圖3(d)可以看出，INLM和本文算法在保持邊緣和細節(jié)方面明顯優(yōu)于前兩個算法，進行細微比較則可以發(fā)現(xiàn)本文算法殘差圖中幾乎看不到殘留的邊緣和細節(jié)信息.

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

(7)

(8)

式中：

σuoriginalu=Cov{uoriginal,u}=

(9)

圖 3 腹部低劑量CT圖與不同算法處理結(jié)果圖的差值圖對比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以及本文算法分別應用于頭部體模而獲得的處理結(jié)果圖相對于原始圖像的PSNR和SSIM值的比較. 由表1很容易看出，本文算法的PSNR值最高，說明本文算法處理結(jié)果圖中含有最少的噪聲； 本文算法的SSIM值最高，說明本文算法處理結(jié)果圖與原圖最接近.

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

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

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

4 結(jié) 論

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