胡紫琪,謝 凱*,文 暢,李美然,賀建飚
(1.長江大學 電子信息學院,湖北 荊州 434023;2.長江大學 電工電子國家級實驗教學示范中心,湖北 荊州 434023;3.長江大學 西部研究院,新疆 克拉瑪依 834099;4.長江大學 計算機科學學院,湖北 荊州 434023;5.中南大學 計算機學院,長沙 410083)
計算機斷層掃描(Computed Tomography,CT)是最實用的成像方式之一,常用于臨床醫學,有助于疾病的診斷與治療;但CT 檢查產生的輻射遠高于普通的X 光檢查,通常情況下通過降低管電流、縮短曝光時間減少輻射對人體的危害,而低輻射劑量會導致CT 圖像中存在噪聲、邊緣不清晰、細節模糊等問題,不利于醫生診判病情[1]。
常用的低劑量CT(Low-Dose CT,LDCT)去噪方法有自適應非局部均值(No-Local Means,NLM)算法[2]。Feruglio等[3]使用三維塊匹配(Block Matching 3D,BM3D)去噪,能明顯改善圖像質量;但CT 圖像中的噪聲分布不均勻,處理完的圖像存在過度平滑和殘差等問題。傳統的醫學去噪算法在速度上有明顯優勢,且不需要大量的數據集,但去噪效果還有待提高。
近些年來,深度神經網絡技術快速發展,為醫學圖像處理提供了新思路和新方法[4-6]。Jain 等[7]用卷積神經網絡(Convolutional Neural Network,CNN)處理圖像去噪問題;Zhang 等[8]提出去噪網絡DNCNN(DeNoising CNN);Chen等[9]結合卷積自編碼器-解碼器與殘差學習提出RED-CNN(Residual Encoder-Decoder CNN),從低劑量估計正常劑量的CT 圖像。上述網絡結構顯著地提升了去噪效果,但以均方差或絕對值誤差作為損失函數,去噪后的圖像整體較為平滑,存在細節丟失、邊緣模糊的問題。生成對抗網絡(Generative Adversarial Network,GAN)[10]是基于博弈思想的網絡框架,由生成器和判別器組成,在對抗損失函數的迭代優化下,能顯著提升去噪效果;但原始GAN 訓練較為困難,存在模式崩塌、梯度消失等問題。研究學者結合GAN 與感知損失[11],提出WGAN(Wasserstein GAN)[12]解決GAN 訓練困難的問題,并采用WGAN-GP(WGAN-Gradient Penalty)[13]使網絡快速收斂。Yang 等[14]提出的WGAN-VGG(WGAN-Visual Geometry Group)也取得了良好的去噪效果,但該網絡只考慮了CT 圖像特征一致性,未考慮結構一致性。
醫生在使用CT 圖像診斷的過程中,需要觀察病變組織、正常組織及其他區域,對圖像的亮度及對比度有著較高的要求,因此在解決CT 圖像的噪聲問題后,還要對去噪后CT 圖像進行增強,包括邊緣增強和動態灰度增強。為解決以上問題,本文提出一種LDCT 圖像增強算法LDCT-IEA。該算法將WGAN-GP 與感知損失、結構損失相結合以提升去噪效果;對去噪后的CT 圖像分別進行動態灰度范圍增強和邊緣增強;最后通過非下采樣輪廓波變換Non-Subsampled Contourlet Transform,NSCT)分解增強后圖像,得到高低頻子圖,配對后的子圖用CNN 進行自適應融合,反變換后就得到最終增強結果。
本文算法LDCT-IEA 流程如圖1,包含圖像去噪、圖像增強和圖像融合三個部分:1)通過修改傳統GAN 的生成器、判別器及損失函數,對LDCT 圖像進行去噪;2)對去噪后的圖像進行動態灰度增強以及邊緣增強;3)使用NSCT 和CNN 融合增強后的圖像。

圖1 本文算法LDCT-IEA流程Fig.1 Flowchart of proposed algorithm LDCT-IEA
1.1.1 生成器
生成器是結構對稱的漏斗形全卷積網絡,包含編碼部分與解碼部分。編碼部分包含5 層Conv2D-LeakyReLU 組合,分別對應卷積層、激活函數層。卷積層的卷積核大小均為5×5,通道數均為128,固定步長為2,LeakyReLU 的α參數為0.2,卷積層不使用填充操作。解碼部分與編碼部分相對應,包含5 層Resize-Conv2D-LeakyReLU 組合(最后一層激活函數為Tanh),本文使用Resize+Conv2D 的操作代替反卷積操作,避免棋盤格效應[15]。該部分卷積層的卷積核大小、通道數與編碼部分一致,固定步長為1,采用填充操作,帶泄露修正線性單元(Leaky Rectified Linear Unit,LeakyReLU)的α參數為0.2。在編碼與解碼的過程中使用跳層連接,將編碼部分和解碼部分中具有相同分辨率的特征圖進行融合,充分利用圖像上下文信息和位置信息,網絡結構如圖2 所示。

圖2 生成器結構Fig.2 Architecture of generator
1.1.2 判別器
判別器通過CNN 提取圖像特征,以判斷判別器輸入是否服從真實分布。因此設計了一個包含5 層Conv2DLayerNormalization-LeakyReLU 組合的卷積網絡,卷積核大小固定為5×5,首層卷積通道數為32,之后的卷積通道數均為前一層的2 倍,固定步長為2,不采用填充操作,LeakyReLU的α參數為0.2。特征圖經過卷積層處理后,通過全局平均池化層和全連接層輸出結果(全連接層的神經元節點數為1),結構如圖3 所示。

圖3 判別器結構Fig.3 Architecture of discriminator
1.1.3 損失函數
GAN 通過優化生成分布pg與真實分布pr之間的JS(Jensen-Shannon)散度來訓練網絡,使pg逼近pr;但pg與pr分布都是高維空間的低維流形,不存在重疊部分或重疊部分可以忽略,導致JS 散度為常數,在梯度求導的過程中梯度為0,產生梯度消失問題。為了從根本上解決GAN 訓練不穩定的問題,使用EM(Earth Mover)距離作為兩個分布的度量指標。在參數化的神經網中中,EM 表達式如下:

其中:Dθ為參數化的CNN。為使Dθ滿足1 階Lipschitz 約束條件,使用梯度懲罰項,定義如下:

其中x′為生成樣本xr與真實樣本xg的線性插值,定義如下:

結合EM 距離后本文判別網絡的損失函數如下:

其中:Dθ表示判別網絡;pg表示生成的常規劑量CT(Normal-Dose CT,NDCT)圖像分布,pr表示真實的NDCT 圖像分布;x′表示生成的NDCT 圖像與真實NDCT 圖像的線性插值;α表示權重系數,設為10。
本文生成網絡的損失函數為:

WGAN-GP 使得生成分布逼近真實分布,從而達到去噪效果;但圖像之間存在特征的差異,為保證圖像之間的內在相似性,本文使用感知損失優化生成器。在生成NDCT 圖像后,將其與對應的真實標簽同時送入已經訓練好的VGG16網絡中,分別在該網絡的淺層、中層及深層輸出特征差值之和,表達式如下:

其中:xr表示真實的NDCT 圖像,xg表示生成的NDCT 圖像;Wn、Hn、Cn分別表示某層卷積的特征圖的寬、高及通道數;F(·)表示特征提取。結合WGAN-GP對抗損失,目標函數為:

由于CNN 不是恒等映射,所以感知損失不能保證圖像間結構的一致性,即在F(xg)-F(xr)=0 的條件下,也不能保證xg=。為保證圖像結構的一致性,本文使用基于像素點的損失函數絕對值誤差來優化生成器:

這樣即滿足特征的相似性也保證圖像結構的一致性,能更好地保留圖像細節信息。總的目標函數如下:

其中:λ和β為歸一化因子,保證感知損失和絕對值誤差損失與對抗損失量綱的一致性。
為提升去噪后CT 圖像的視覺效果,對去噪后的圖像進行一定程度的增強,包括提升動態灰度范圍、增強邊緣輪廓。
1.2.1 動態灰度增強
Retinex[17]是對人類視覺系統進行建模的圖像增強算法,在此基礎上衍生出了許多有效的圖像增強算法,其中常用的SSR(Single Scale Retinex)、MSR(Multi-Scale Retinex)算法均存在偏色問題,為提升圖像的動態灰度范圍,本文使用基于MSR 的MSRCR(Multi-Scale Retinex with Color Restoration)[18]算法。MSR 算法的定義如下:

其中:R(x,y)為增強后CT 圖像;I(x,y)為待增強CT 圖像。通常情況下,環境函數F(x,y)為高斯卷積函數,定義如下:

其中c、λ使得以下條件成立:

在此MSR 基礎上,MSRCR 算法加入調整因子C,以此來解決降低色彩失真問題,定義如下:

其中:Ii(x,y)表示圖像的某個通道;α為受控的非線性強度;β表示增益常數。
1.2.2 邊緣增強
圖像邊緣是灰度值變化較大的位置,描述了圖像局部不連續的特征,也是區分各種組織、骨頭及鈣化斑塊的邊界。圖像邊緣增強使得模糊的邊界更加清晰,有利于區分各種組織部位。常用的一階微分算子存在邊緣定位效果較差、易丟失關鍵信息等不足;二階微分算子相較于一階可以提取更多邊緣信息,但是對噪聲較為敏感。因此本文使用非銳化掩膜法增強圖像邊緣。算法主要分為三個部分:平滑圖像、生成銳化模板以及銳化模板與原圖相加,公式如下:

其 中:f(x,y) 表示輸入圖像表示平滑后圖像;gmask(x,y)表示模板;λ表示控制增強效果的縮放因子。使用高斯濾波得到平滑結果,計算方法如下:

其中:s、t表示高斯模板坐標,模板中心為原點;G(s,t)表示二維高斯函數,表達式如式(18)。

其中σ表示高斯函數標準差。
為了使CT 圖像同時包含亮度增強和邊緣輪廓增強兩種增強效果,且避免過度增強和信息冗余的問題,本文將增強后的CT 圖像進行融合。常用的融合方法可分為空域法和變換域法。基于像素點的空域融合方法細節輪廓丟失嚴重,融合后圖像的信噪比和對比度降低;相較于空域法,變換域法是一種更為細致的融合方法,能在圖像的變換域上融合不同的高低頻信息。常用的融合算法是多尺度變換結合人工設計的融合規則,但是CT 圖像的輪廓、病變區等一些目標區域常常是不規則的形狀,而人工設計的融合規則一般基于頻域系數或者一些局部窗口,不是根據不規則目標區域設計的,局限性較大。因此本文采用多波段自適應融合的思想,使用CNN 的自適應特性,結合分解特性較好的NSCT,對增強后的CT 圖像在頻域上進行自適應融合。
1.3.1 非下采樣輪廓波變換
NSCT[19]是輪廓波變換的改進,利用非下采樣金字塔濾波器和非下采樣方向濾波器替換了原有濾波器,使該變換具有平移不變性,解決了輪廓波因下采樣操作導致的偽吉布斯現象。第一部分的非下采用金字塔分解使圖像具有多分辨率性,并輸出一個低頻子圖和一個高頻子圖;第二部分的非下采樣方向濾波器是實現圖像多方向性的關鍵,由雙通道扇形的非下采用濾波器組成,將高頻子圖分解為多個多方向子圖,每個子圖表示圖像在某個方向上的細節信息。若對原圖像進行m級分解,可得到1 個低頻子圖及個高頻子圖,其中ki表示在尺度下i下的多方向分解級數。分解流程如圖4 所示。

圖4 非下采樣輪廓波變換Fig.4 Non-subsampled contourlet transform
1.3.2 CNN自適應融合
為解決人工設計融合規則存在的問題,本文使用特征提取能力較強的CNN,設計了端到端的自適應融合方法,將分解后的高低頻系數子圖配對送入網絡得到融合結果。該網絡由兩個相同的CNN 構成,輸入大小均為512×512,并對兩個子圖進行歸一化操作,結果范圍為[0,1]。該網絡結構如圖5 所示,由4 個Conv-BN-ReLU 組合構成,其中卷積層的卷積核大小為5×5,步長為1,采用Padding 操作。

圖5 圖像融合網絡結構Fig.5 Architecture of image fusion network
經過特征提取后,兩個網絡分別給兩個子圖不同的權重分配,進而進行特征融合,再將融合結果進行反歸一化,通過非下采樣輪廓波反變換得到最終結果。公式表達如下:

其中:NSCTL(x)、NSCTH(x)表示非下采樣輪廓波變換后的低頻子圖和高頻子圖;fL、fH表示CNN 的映射函數;⊕表示特征融合;PL、PH、Pfinal分別表示低頻、高頻及最后的融合結果,NSCT()-1表示反變換。
本文實驗平臺如下:操作系統為Windows 10,顯卡為NVIDIA RTX2060,處理器為Intel Core i7-10750H;使用TensorFlow2.0 深度學習框架以及Python 3.7.4 搭建網絡模型。
實驗數據來自2016 年Low Dose CT Grand Challenge 比賽的公開數據集[20],為真實臨床數據,包括患者的腹部掃描CT以及對應的模擬1/4 劑量下的LDCT。本文隨機選取10 個匿名患者的CT 數據共2 780 對,CT 圖像大小為均為512×512,切片厚度為3 mm。訓練集、測試集劃分比例為8∶2,訓練集包含2 224 對數據,測試集包含556 對數據。在訓練階段,考慮到計算量問題,對訓練集的每對圖像進行5 次隨機裁剪,裁剪后的圖像尺寸均為64×64,當裁剪的圖像均值小于0.1時,對其進行重新裁剪,保證訓練數據的可靠性,最終得到13 900 對訓練數據。在測試階段,不使用裁剪操作,輸入圖像大小為512×512。
2.2.1 圖像去噪
為了評估本文去噪算法,選取一張典型的低劑量腹部CT 圖像,與常用的去噪算法BM3D、RED-CNN、WGAN-VGG作對比,并進行損失函數消融實驗。根據實驗結果調整優化后,實驗參數設置如下:訓練次數為500,訓練批次為64,生成器、判別器的學習率分別設置為0.000 1、0.000 4;優化器為Adam,其超參數設定β1=0.5、β2=0.9;歸一化參數β和λ分別為10 以及0.1。圖6 為不同算法的去噪結果對比。

圖6 各算法的去噪結果對比Fig.6 Comparison of denoising results of different algorithms
如圖6(a)、(b)所示,相較于NDCT 圖像,LDCT 圖像因受到量子噪聲的干擾,結構信息和細節信息損失較為嚴重,影響醫生的臨床診斷。由圖6(c)~(e)可知,BM3D、RED-CNN和WGAN-VGG 方法均在不同程度上抑制了噪聲,基于深度學習方法的RED-CNN 和WGAN-VGG 的去噪效果明顯優于BM3D,不過RED-CNN 的結果過于平滑,感興趣邊緣被模糊。
由圖6(f)可知,本文算法不僅有效去除了大部分噪聲及偽影,且保留了較多的細節信息。由于結合了結構損失函數,本文算法相較于WGAN-VGG,在感興趣區域保留了更多的結構信息。為了進一步說明本文去噪算法的有效性,對比了圖6 中所有去噪結果的第200 行輪廓線,如圖7 所示,結果直觀地展現了本文算法的去噪結果與NDCT 圖像的擬合程度最高,體現了本文算法的結構一致性。

圖7 圖6中所有去噪結果的第200行輪廓線對比Fig.7 Comparison of contour line of the 200th row among all denoising results in Fig.6
另外,本文還定量分析了不同算法的去噪性能,使用三種常用的去噪評價指標:峰值信噪比(Peak Signal-to-Noise Ratio,PSNR)、均方根誤差(Root Mean Square Error,RMSE)和結構相似度(Structural Similarity Index Measure,SSIM)作為評價標準,結果如圖8 所示。由圖8 可知:本文算法的PSNR取得了最大值,去噪效果優于對比算法;SSIM 也取得了最大值,RMSE 則取得了最小值,表明去噪后的結果與NDCT 圖像接近,體現了結構上的一致性。

圖8 不同算法的RMSE、PSNR、SSIM對比Fig.8 Comparison of RMSE,PSNR,SSIM of different algorithms
驗證損失函數去噪優勢的損失函數消融實驗結果如圖9 及表1 所示。

圖9 損失函數消融實驗結果對比Fig.9 Comparison of experimental results of loss function ablation

表1 消融實驗結果對比Tab.1 Comparison of ablation experimental results
在損失函數消融實驗上,LossGAN+Lossp與LossGAN+LossMAE兩種方法均在一定程度上去除了噪聲,但是前者更加注重圖像特征之間的相近,去噪后的CT 圖像在結構上有部分損失,后者更加關注圖像之間的結構一致性,圖像存在平滑現象。而本文方法的損失函數既考慮了CT 圖像間特征一致性又考慮了圖像間的結構一致性,在去噪效果及細節信息保留上均要優于上述兩種對比方法。表1 是消融實驗的評價指標對比,本文方法在三種指標上均取得最優。
2.2.2 圖像增強
動態灰度范圍增強實驗結果如圖10、11 所示,其中本文使用MSRCR 算法。可以看出,測試方法均增大了去噪后CT圖像的動態灰度范圍,亮度都有所增強。
圖10(a)為去噪后的CT 圖像(即圖6(f)),如圖10(b)~(f)及圖11(b)~(f)所示,基于直方圖均衡的方法擴大后的灰度范圍相較于其他方法要窄,整體亮度較低,且部分區域模糊;SSR 的結果整體較亮,導致各組織部位區分不明顯;MSR灰度分布不均勻導致背景亮度過高,存在色彩失衡現象;而MSRCR 灰度分布均勻,前景、背景分隔明顯,增強亮度的同時結構清晰。

圖10 動態灰度增強結果Fig.10 Dynamic gray-scale enhancement results

圖11 灰度直方圖統計Fig.11 Gray histogram statistics
邊緣增強實驗對比了一些常用的邊緣增強算法,如:Robert、Sobel、Laplacian、LOG(Laplacian Of Gaussian),實驗結果如圖12、13 所示,圖12(a)為去噪后的CT 圖像(即圖6(f))。從圖12 可以看出,測試方法都實現了邊緣增強。結合相應圖13 的三維灰度熱力圖(顏色的深淺表示像素值的大小,顏色越深該位置的灰度值越大)可知:Robert、Sobel算子邊緣定位不準確,導致外圍輪廓增強較為明顯,其他結構細節部分效果較差,且存在邊緣輪廓模糊的現象;Laplacian、LOG 算子相較于一階算子邊緣定位更加準確,邊緣增強效果較為明顯,但放大了噪聲信號,且但存在過度增強現象;非銳化掩膜在灰度變化較小的邊界處增強效果明顯,整體的邊緣增強效果較好,圖像結構清晰。

圖12 邊緣增強結果Fig.12 Edge enhancement results

圖13 三維熱力圖Fig.13 Three-dimensional heat maps
2.2.3 圖像融合
圖像融合實驗使用非下采用輪廓波結合CNN 的方法與常用的融合方法,如像素點融合、頻域融合及CNN 融合作對比。為提升融合效果,實驗中對圖像進行兩級非下采樣輪廓波變換,分解后包含1 張低頻子圖及6 張高頻子圖。為了訓練本文的CNN,手動對NDCT 實驗數據集同時進行動態灰度增強及對比度增強,調整窗位(window level)及窗寬(window width),設置為45 和350,分別訓練高低頻融合網絡,使用MSE 作為優化損失函數,優化器、學習率、訓練批次、訓練次數分別為Adam、0.000 1、8、100,實驗結果如圖15所示。從圖15 可以看出,基于像素點及非下采樣輪廓波融合的結果整體噪聲較為明顯,PSNR 下降;而其他兩種方法的融合結果更為細致。為客觀分析融合效果,除PSNR、SSIM、RMSE 外,還使用了基于圖像特征的三種評價指標,包括:平均梯度(Average Gradient,AG)、空間頻率(Spatial Frequency,SF)和標準差(Standard Deviation,SD)。AG 用于衡量圖像的清晰度,SF 用于衡量圖像細節信息豐富度,SD 則反映了灰度的離散情況。分析圖15、16 可知,本文方法在PSNR、SSIM、RMSE上的結果分別為33.015 5 dB、0.918 5、5.99,在AG、SF、SD 上的結果分別為109.662 7、43.407 6、28.294 2,對比其他融合方法,上述參數均為最優值。

圖14 融合結果對比Fig.14 Comparison of fusion results

圖15 融合結果的SSIM、PSNR、RMSE對比Fig.15 Comparison of RMSE,PSNR,SSIM of fusion results

圖16 融合結果的AG、SF、SD對比Fig.16 Comparison of AG,SF,SD of fusion results
相較于RED-CNN 及WGAN-VGG,本文方法LDCT-IEA 的RMSE 分別降低了2.68%、3.63%,PSNR 分別提高了2.088%、2.358%,SSIM 分別提高了1.413%、1.773%。綜合實驗結果可知,本文方法不僅去除了LDCT圖像中的噪聲,且增強了CT圖像的動態灰度范圍及輪廓邊緣,使CT 圖像更加清晰,視覺效果更好,更加有利于醫生診斷病情。
本文提出了生成對抗網絡下的LDCT 圖像增強算法,不僅考慮了LDCT 圖像的噪聲影響,還考慮了亮度及對比度。通過綜合改進WGAN-GP 算法,有效地提升對LDCT 圖像的去噪效果,較好地保留了圖像細節信息;在圖像融合方面,為使去噪增強后的CT 圖像能夠有效實現信息互補,消除過度增強,將非下采樣輪廓波結合CNN 實現頻域上子圖的自適應融合,融合后不僅有效地減少了噪聲,且明顯地提升了動態灰度范圍、增強了圖像的邊緣信息。實驗結果也驗證了本文方法的有效性和可行性。
今后還可以對算法進行更深入的探索來適應未來更加多樣化的CT 數據集;將算法擴展到X 光片、核磁共振等醫學圖像,為醫療輔助提供可靠的算法。