原培新， 陳鼎夫

(東北大學(xué) 機(jī)械工程與自動化學(xué)院， 遼寧 沈陽 110819)

X射線檢測在國家公共安防領(lǐng)域扮演著特別重要的角色.它通過透射的方法掃描旅客的行李物品，并生成透照圖像使安檢人員可以在不打開旅客個(gè)人物品時(shí)查看是否攜帶違禁物.

現(xiàn)今安檢最常用的是雙能X射線檢測系統(tǒng).它通過產(chǎn)生高能與低能兩種不同能量的射線對同一物品進(jìn)行掃描，生成高、低能兩幅圖像，并各自包含部分微小細(xì)節(jié)，因此須同時(shí)對它們進(jìn)行處理[1].

X射線安檢圖像噪聲大、細(xì)節(jié)豐富，但最主要的問題在于其掃描出來的數(shù)據(jù)為16位數(shù)據(jù)，具有較高的動態(tài)范圍，但高動態(tài)范圍顯示器的造價(jià)成本昂貴，應(yīng)用于安檢設(shè)備性價(jià)比較低.若要將高動態(tài)范圍的圖像直接顯示在普通顯示器上，則需要進(jìn)行動態(tài)范圍的壓縮，降低檢測數(shù)據(jù)的位數(shù)，這勢必會導(dǎo)致數(shù)據(jù)信息的丟失，使顯示出來的圖像丟失部分細(xì)節(jié).同時(shí)，由于安檢圖像的對比度較低，圖像在黑暗區(qū)域、高亮區(qū)域以及物體重合區(qū)域顯示不清，增加了安檢人員的工作難度.因此，在壓縮X射線圖像動態(tài)范圍的同時(shí)，如何增強(qiáng)圖像細(xì)節(jié)信息、提高對比度是目前所需解決的重點(diǎn).本文提出了一種基于多尺度局部保邊(local edge-preserving,LEP)濾波的自適應(yīng)色調(diào)映射算法以解決以上問題.

1 雙能X射線的色調(diào)映射算法

現(xiàn)今對X射線圖像的動態(tài)范圍壓縮以及圖像增強(qiáng)的研究主要集中于工業(yè)探測領(lǐng)域以及醫(yī)療醫(yī)學(xué)領(lǐng)域[2]，而對X射線安檢圖像的處理研究較少.

常見的色調(diào)映射算法有線性壓縮法、對數(shù)函數(shù)壓縮法、梯度壓縮法[3]、基于雙邊濾波法等.線性壓縮法是將圖像每個(gè)像素點(diǎn)的灰度值線性縮小至原來的1/256，但這樣做會造成對比度降低，灰度過于集中；對數(shù)函數(shù)壓縮法能夠保持圖像的整體明亮效果，但會丟失大量細(xì)節(jié)信息；梯度壓縮法在處理完圖像之后，還需要進(jìn)行泊松方程的求解，較為費(fèi)時(shí)，不適用于安檢這種比較注重實(shí)時(shí)性的領(lǐng)域.

基于雙邊濾波[4]的算法以及基于指導(dǎo)濾波、基于加權(quán)最小二乘濾波[5]的算法都與本文算法有相似之處，它們都是通過濾波將原圖像分解為基礎(chǔ)層與細(xì)節(jié)層，不足之處在于雙邊濾波算法處理后會引起邊緣過度銳化，造成梯度反轉(zhuǎn)且會出現(xiàn)“階梯效應(yīng)”；指導(dǎo)濾波算法不能很好地保持邊緣，而加權(quán)最小二乘濾波算法的平滑效果較好，但對細(xì)節(jié)的提取較弱.

本文綜合上述問題，采用局部邊緣保持濾波，同時(shí)對高、低能圖像進(jìn)行多尺度濾波，分別得到1個(gè)基礎(chǔ)層以及3個(gè)細(xì)節(jié)層；對細(xì)節(jié)層進(jìn)行自適應(yīng)的增強(qiáng)，對基礎(chǔ)層進(jìn)行自適應(yīng)對數(shù)壓縮，并將處理后的各層圖像合并，進(jìn)行直方圖調(diào)整.最后將高、低能圖像融合得到細(xì)節(jié)清晰、對比度高、噪聲小的低動態(tài)范圍雙能X射線安檢圖像，具體算法流程如圖1所示.

1.1 圖像預(yù)處理

首先分別將安檢機(jī)檢測到的高、低能數(shù)據(jù)代入至算法中，由于檢測數(shù)據(jù)為16位數(shù)據(jù)，數(shù)據(jù)范圍為[0,65 535]，因此首先對數(shù)據(jù)進(jìn)行歸一化處理，處理公式為

I=Iin/maxIin.

(1)

式中：Iin代表輸入圖像的灰度值；maxIin代表輸入圖像灰度值的最大值；I代表歸一化處理后的灰度值.

圖1 算法流程圖

1.2 局部保邊(LEP)濾波

局部保邊(LEP)濾波[6]與加權(quán)最小二乘濾波的想法一致，都希望濾波后的圖像在具有較大梯度區(qū)域以外的區(qū)域能夠盡可能平滑，同時(shí)使處理后的圖像與原圖像接近一致.但局部邊緣保持濾波對顯著邊緣的定義不同，它認(rèn)為顯著邊緣不僅僅是全局大梯度區(qū)域，在局部具有相對較大梯度的部分也被歸為圖像的顯著邊緣，因此在LEP濾波中，具有局部較大梯度的區(qū)域也被分解至基礎(chǔ)層.

局部邊緣保持濾波可以簡單理解為求式(2)的最小值：

(2)

式中：Bi,j表示輸出圖像的(i,j)像素點(diǎn)的灰度；Ii，j為輸入圖像(i,j)像素點(diǎn)的灰度；w為濾波窗口；Ii，j表示(i,j)像素點(diǎn)處I的梯度；Bi,j表示(i,j)像素點(diǎn)處B的梯度；α為可調(diào)節(jié)參數(shù)，β用于調(diào)節(jié)梯度敏感度.對輸入圖像進(jìn)行局部濾波以保留局部邊緣.

受指導(dǎo)濾波[7]的啟發(fā)，LEP濾波認(rèn)為，在濾波過程中輸入與輸出圖像呈線性關(guān)系，因此有

Bi，j=awIi,j+bw,(i,j)∈w.

(3)

聯(lián)立式(2)與式(3)并求偏導(dǎo)可得

(4)

在每一個(gè)濾波窗口w中都含有N個(gè)像素，而每個(gè)像素又會被N個(gè)窗口所包含.每一個(gè)窗口都會有一組aw，bw，因此對輸出Bi來說，也會有N種結(jié)果，對這N種結(jié)果取均值即可得到輸出.

(5)

1.3 多尺度分解

圖像多尺度分解是將圖像經(jīng)濾波器濾波后分解為一個(gè)基礎(chǔ)層以及多個(gè)細(xì)節(jié)層的過程.

常用的多尺度分解方法分為兩種:第一種為對原始輸入圖像進(jìn)行多次濾波，每次選擇不同大小的濾波窗口，之后對每次濾波后的細(xì)節(jié)層進(jìn)行處理并最終將它們?nèi)诤?；第二種方法是對原始圖像進(jìn)行迭代地濾波分解，即每次都是對上一次濾波后輸出的基礎(chǔ)層進(jìn)行再次濾波，并依次增大濾波窗口.兩種方法都能夠得到一系列逐步平滑的圖像，但針對安檢圖像，第一種方法的效果并不明顯，且易出現(xiàn)拼接痕跡，因此本文選擇迭代分解法.

對原始輸入圖像進(jìn)行一次LEP濾波后可得到一個(gè)輸出圖像，即為基礎(chǔ)層，將輸入圖像與輸出圖像作差可得一個(gè)細(xì)節(jié)層.基礎(chǔ)層保留著原始圖像的局部均值以及局部突出邊緣，而細(xì)節(jié)層包含著在零梯度周圍震蕩的信號.迭代地使用濾波器對圖像進(jìn)行濾波分解并依次增大濾波半徑可得一個(gè)基礎(chǔ)層以及多個(gè)細(xì)節(jié)層，實(shí)驗(yàn)表明：進(jìn)行3次迭代濾波后可近乎完全提取圖像細(xì)節(jié)，且運(yùn)行時(shí)間較短，因此將濾波尺度定為3，具體表達(dá)式如下：

Bl=Ll(Bl-1),

(6)

Dl=Bl-Bl-1，

(7)

Iout=B3+D1+D2+D3.

(8)

式中Bl,Dl,Iout分別代表基礎(chǔ)層、細(xì)節(jié)層和輸出圖像，B0為原始輸入圖像；Ll代表濾波函數(shù)；l代表濾波次數(shù)，l=1～3.

1.4 細(xì)節(jié)層自適應(yīng)增強(qiáng)

在圖像多尺度分解以后，需要對不同尺度下得到的細(xì)節(jié)層進(jìn)行增強(qiáng)，提升對比度，從而達(dá)到增強(qiáng)圖像細(xì)節(jié)信息的目的.

由于細(xì)節(jié)層包含許多在零梯度周圍震蕩的信息，因此希望找到一個(gè)函數(shù)，能夠壓縮遠(yuǎn)離0的較大值，同時(shí)增強(qiáng)較小值.此外，這樣的增強(qiáng)函數(shù)還應(yīng)使圖像中每一點(diǎn)的偏離大致相等，且為了避免梯度反轉(zhuǎn)，增強(qiáng)函數(shù)還應(yīng)為凸函數(shù)，并關(guān)于原點(diǎn)對稱.

滿足上述條件的函數(shù)有Sigmoid函數(shù)[8]，如式(9)所示：

(9)

其形狀大致如圖2所示.

圖2 Sigmoid函數(shù)

一般的Sigmoid函數(shù)都能夠增強(qiáng)細(xì)節(jié)層，但是一些殘留在細(xì)節(jié)層里的強(qiáng)邊緣也會隨著細(xì)節(jié)的增強(qiáng)而增強(qiáng)，當(dāng)它們與基礎(chǔ)層相加后就會引起較強(qiáng)的光暈現(xiàn)象，因此本文在Sigmoid函數(shù)里加入了梯度項(xiàng)，用來控制較大梯度的強(qiáng)邊緣帶來的影響，改進(jìn)函數(shù)如下：

(10)

本文算法應(yīng)用式(10)對圖像的細(xì)節(jié)層進(jìn)行增強(qiáng)，同時(shí)能夠盡可能避免光暈的出現(xiàn).

1.5 基礎(chǔ)層對數(shù)自適應(yīng)壓縮

本文算法對3次迭代后得到的基礎(chǔ)層采用自適應(yīng)對數(shù)函數(shù)進(jìn)行動態(tài)范圍的調(diào)整，原因在于對數(shù)函數(shù)符合人眼對物體亮度的感知，同時(shí)，物質(zhì)對X射線能量的衰減呈指數(shù)形式，因此在對數(shù)域上操作能夠得到更好的表現(xiàn).

Stockham于1972年首次提出應(yīng)用對數(shù)函數(shù)對圖像的亮度進(jìn)行調(diào)整[9]，但其算法由于設(shè)置了一個(gè)固定的對數(shù)基底，無法同時(shí)對黑暗及高亮區(qū)域進(jìn)行有效調(diào)整.Drago等提出了一種自適應(yīng)對數(shù)色調(diào)映射算法[10]，它可以根據(jù)輸入圖像各點(diǎn)像素亮度值的不同，自適應(yīng)地調(diào)整對數(shù)基底，具體變換公式如下：

(11)

式中：Lw為輸入圖像各點(diǎn)亮度值；Lw,max為最大亮度值；La為一比例因子，用于調(diào)整輸出圖像的亮度，這里取La=100 cd/m2；偏置函數(shù)b(x)用于自適應(yīng)調(diào)整對數(shù)基底：

(12)

式中：b為一關(guān)鍵系數(shù)，用于壓縮高亮區(qū)域值并增強(qiáng)黑暗區(qū)域細(xì)節(jié)的可見度，實(shí)驗(yàn)表明b的取值在0.7～0.9之間時(shí)，大部分圖像都能夠取得較好效果，這里取b=0.8.

1.6 直方圖均衡化

在X射線安檢圖像產(chǎn)生以及處理過程中，必然會產(chǎn)生噪聲，因此應(yīng)該進(jìn)行直方圖均衡化[11]調(diào)

整.本文算法對經(jīng)過上述處理后的圖像的灰度級前后各0.5%進(jìn)行調(diào)整并對調(diào)整后圖像進(jìn)行歸一化操作，以減少噪聲對圖像的影響，同時(shí)提升主要像素的對比度.以工具箱為例，處理前后的效果如圖3所示；直方圖的對比如圖4所示.

圖3 直方圖均衡處理前后效果對比

圖4 圖像灰度直方圖調(diào)整前后對比

1.7 雙能X射線圖像融合

行李物品在經(jīng)過雙能安檢設(shè)備后，會在高能探測器以及低能探測器上分別出現(xiàn)一幅包含有物品信息的掃描圖像.有些內(nèi)部結(jié)構(gòu)較厚的地方可能在高能圖像中顯示較為清楚，而有些細(xì)微的結(jié)構(gòu)則可能在低能圖像中才能顯示清楚，因此分別對高、低能安檢圖像進(jìn)行上述處理以后，還需進(jìn)行圖像融合，將高、低能圖像所包含的有用信息結(jié)合起來，生成一幅便于安檢人員觀察的圖像.

常用的X射線融合方法有灰度最大值法、基于小波融合法[12]、加權(quán)平均法等，本文算法將采用主成分分析法對雙能圖像進(jìn)行按權(quán)重融合，具體步驟如下：

1)分別將待融合的高、低能圖像寫成列向量的形式，并組成為一個(gè)n行、2列的矩陣，其中n代表待融合圖像的總像素?cái)?shù).對上述矩陣求協(xié)方差矩陣，形式如下：

(13)

2)求取矩陣Z的特征值與特征向量.

3)設(shè)上述協(xié)方差矩陣中較大的特征值所對應(yīng)的特征向量為(a0,a1)T，則第一幅圖像的融合權(quán)重為a0/(a0+a1)，第二幅圖像的融合權(quán)重為a1/(a0+a1) .

采用上述方法對圖像進(jìn)行融合，耗時(shí)較短，符合安檢圖像處理實(shí)時(shí)性的要求，且沒有較明顯的拼接痕跡.

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

為驗(yàn)證本文算法的可行性，搭建實(shí)驗(yàn)平臺，采用Intel(R) Core(TM) i7-7700HQ CPU 2.80 GHz處理器，16 GB內(nèi)存，Windows 10及64位操作系統(tǒng)，編譯環(huán)境為MATLAB2018b.

實(shí)驗(yàn)設(shè)備采用沈陽地泰設(shè)備檢測有限公司的DEX6550雙視角X射線檢查設(shè)備，該設(shè)備為16位數(shù)據(jù)采集系統(tǒng)，實(shí)驗(yàn)數(shù)據(jù)均來自沈陽地鐵1，2號線的真實(shí)數(shù)據(jù)，原始圖像的動態(tài)范圍為[0,65 535].

針對本文算法所提出的一些可人為修改參數(shù)，取值如下：α=0.1，β=1，濾波次數(shù)l=3；濾波窗口w的半徑依次選為2，12，22；b=0.8.

2.1 主觀評價(jià)

對圖像最直接、最可靠的評價(jià)方法便是主觀評價(jià).對安檢圖像來講，主觀評價(jià)一幅圖像的好壞，主要是看安檢人員能否看清包裹內(nèi)的全部物品；當(dāng)物品之間有重疊或內(nèi)部結(jié)構(gòu)較厚時(shí)，是否能將它們的輪廓全部分清.本文選取了在安檢中較為常見的工具箱和手提包，處理前后的對比如圖5所示.

圖5 本文算法處理前后對比

由圖5a、圖5b可以很清晰地看到，在處理后，圖像整體的對比度效果提升，細(xì)節(jié)得到增強(qiáng)，左上方電路板的紋路顯示更加清晰，下方線圈的微小細(xì)節(jié)也得到顯示；由圖5c、圖5d可以看出，整個(gè)箱內(nèi)的重合物體在處理后，各自的邊緣都能得到較為清晰的顯示.

同時(shí)，將本文算法與對數(shù)壓縮、基于Sigmoid函數(shù)壓縮以及基于LEP濾波壓縮進(jìn)行對比，如圖6所示，也可以清晰看出本文算法對圖像的改善.

2.2 客觀評價(jià)

為了避免受到光照等周圍環(huán)境以及觀察者本身狀態(tài)、心理的影響，本文在對算法進(jìn)行主觀評價(jià)的同時(shí)，還將對調(diào)整后的圖像進(jìn)行客觀評價(jià)[13]，采用的評價(jià)標(biāo)準(zhǔn)有信息熵、對比度提升指數(shù)以及峰值信噪比.

圖6 不同算法對比

信息熵[14]描述的是一幅圖像所含信息的豐富程度，其值越大，則表示圖像中所包含的信息越豐富；對比度提升指數(shù)[15]是針對射線圖像對比度普遍較低所提出的用來評價(jià)X射線在增強(qiáng)后對比度變化的指標(biāo)；峰值信噪比[16]用來表示圖像處理后的失真情況，其值越大，圖像失真越少.

信息熵的計(jì)算公式為

(14)

式中：Pi表示灰度級為i的概率；G表示圖像的總灰度級.

對比度提升指數(shù)的計(jì)算公式為

(15)

(16)

(17)

式中：Cori,Cpro分別表示處理前與處理后圖像；m(·)為取均值；fmax(i,j)，fmin(i,j)表示輸入圖像的像素點(diǎn)(i,j)在3×3窗口內(nèi)像素最大、最小值；gmax(i,j)，gmin(i,j)表示增強(qiáng)后像素的最大、最小值.

峰值信噪比的計(jì)算公式為

(18)

式中：M×Q表示原圖尺寸(像素?cái)?shù))；f(i,j)，g(i,j)分別代表處理前后圖像的灰度值.

本文利用上述3個(gè)客觀評價(jià)標(biāo)準(zhǔn)對處理后圖像進(jìn)行評價(jià)，并對比本文算法與其他算法之間的差異，同時(shí)，為了規(guī)避實(shí)驗(yàn)數(shù)據(jù)較少給結(jié)果帶來的偶然性，另外選取三幅圖片進(jìn)行實(shí)驗(yàn)，所得結(jié)果如表1所示.

表1 不同算法客觀評價(jià)

由表1可以看出，相較于其他色調(diào)映射算法，經(jīng)本文所提算法進(jìn)行動態(tài)范圍調(diào)整后的圖像，其信息熵、對比度提升指數(shù)、峰值信噪比都有顯著提高，同時(shí)本文算法的平均處理時(shí)間在1.2 s左右，符合安檢工作的實(shí)時(shí)性.這表明本文算法不僅能夠有效壓縮動態(tài)范圍，同時(shí)能夠增強(qiáng)圖像所蘊(yùn)含的細(xì)節(jié)信息、提高圖像對比度及整體明暗程度，得到比其他算法更好的低動態(tài)范圍雙能X射線安檢圖像.

3 結(jié) 論

針對高動態(tài)范圍雙能X射線圖像的特點(diǎn)，本文從實(shí)際出發(fā)，結(jié)合安檢工作中對輸出圖像的實(shí)際需求，應(yīng)用局部保邊濾波器對輸入圖像進(jìn)行迭代濾波，并對得到的細(xì)節(jié)層以及基礎(chǔ)層同時(shí)進(jìn)行自適應(yīng)處理，之后將它們相加.最后進(jìn)行直方圖均衡化操作并按照權(quán)重融合高、低能圖像，得到低動態(tài)范圍雙能X射線安檢圖像.實(shí)驗(yàn)結(jié)果表明，本文算法在壓縮動態(tài)范圍的同時(shí)，能夠改善安檢圖像的缺點(diǎn)，得到細(xì)節(jié)清晰、對比度高、更加便于安檢人員觀察的圖像，符合實(shí)際要求.后續(xù)可將本算法針對實(shí)時(shí)線掃描系統(tǒng)，并結(jié)合雙視角進(jìn)行改進(jìn)，得到實(shí)時(shí)無痕的清晰圖像.