雙能X射線高動態范圍安檢圖像壓縮算法

原培新， 陳鼎夫

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

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

現今安檢最常用的是雙能X射線檢測系統.它通過產生高能與低能兩種不同能量的射線對同一物品進行掃描，生成高、低能兩幅圖像，并各自包含部分微小細節，因此須同時對它們進行處理[1].

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

1 雙能X射線的色調映射算法

現今對X射線圖像的動態范圍壓縮以及圖像增強的研究主要集中于工業探測領域以及醫療醫學領域[2]，而對X射線安檢圖像的處理研究較少.

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

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

本文綜合上述問題，采用局部邊緣保持濾波，同時對高、低能圖像進行多尺度濾波，分別得到1個基礎層以及3個細節層；對細節層進行自適應的增強，對基礎層進行自適應對數壓縮，并將處理后的各層圖像合并，進行直方圖調整.最后將高、低能圖像融合得到細節清晰、對比度高、噪聲小的低動態范圍雙能X射線安檢圖像，具體算法流程如圖1所示.

1.1 圖像預處理

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

I=Iin/maxIin.

(1)

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

圖1 算法流程圖

1.2 局部保邊(LEP)濾波

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

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

(2)

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

受指導濾波[7]的啟發，LEP濾波認為，在濾波過程中輸入與輸出圖像呈線性關系，因此有

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

(3)

聯立式(2)與式(3)并求偏導可得

(4)

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

(5)

1.3 多尺度分解

圖像多尺度分解是將圖像經濾波器濾波后分解為一個基礎層以及多個細節層的過程.

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

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

Bl=Ll(Bl-1),

(6)

Dl=Bl-Bl-1，

(7)

Iout=B3+D1+D2+D3.

(8)

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

1.4 細節層自適應增強

在圖像多尺度分解以后，需要對不同尺度下得到的細節層進行增強，提升對比度，從而達到增強圖像細節信息的目的.

由于細節層包含許多在零梯度周圍震蕩的信息，因此希望找到一個函數，能夠壓縮遠離0的較大值，同時增強較小值.此外，這樣的增強函數還應使圖像中每一點的偏離大致相等，且為了避免梯度反轉，增強函數還應為凸函數，并關于原點對稱.

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

(9)

其形狀大致如圖2所示.

圖2 Sigmoid函數

一般的Sigmoid函數都能夠增強細節層，但是一些殘留在細節層里的強邊緣也會隨著細節的增強而增強，當它們與基礎層相加后就會引起較強的光暈現象，因此本文在Sigmoid函數里加入了梯度項，用來控制較大梯度的強邊緣帶來的影響，改進函數如下：

(10)

本文算法應用式(10)對圖像的細節層進行增強，同時能夠盡可能避免光暈的出現.

1.5 基礎層對數自適應壓縮

本文算法對3次迭代后得到的基礎層采用自適應對數函數進行動態范圍的調整，原因在于對數函數符合人眼對物體亮度的感知，同時，物質對X射線能量的衰減呈指數形式，因此在對數域上操作能夠得到更好的表現.

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

(11)

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

(12)

式中：b為一關鍵系數，用于壓縮高亮區域值并增強黑暗區域細節的可見度，實驗表明b的取值在0.7～0.9之間時，大部分圖像都能夠取得較好效果，這里取b=0.8.

1.6 直方圖均衡化

在X射線安檢圖像產生以及處理過程中，必然會產生噪聲，因此應該進行直方圖均衡化[11]調

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

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

圖4 圖像灰度直方圖調整前后對比

1.7 雙能X射線圖像融合

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

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

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

(13)

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

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

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

2 實驗結果對比與分析

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

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

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

2.1 主觀評價

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

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

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

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

2.2 客觀評價

為了避免受到光照等周圍環境以及觀察者本身狀態、心理的影響，本文在對算法進行主觀評價的同時，還將對調整后的圖像進行客觀評價[13]，采用的評價標準有信息熵、對比度提升指數以及峰值信噪比.

圖6 不同算法對比

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

信息熵的計算公式為

(14)

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

對比度提升指數的計算公式為

(15)

(16)

(17)

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

峰值信噪比的計算公式為

(18)

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

本文利用上述3個客觀評價標準對處理后圖像進行評價，并對比本文算法與其他算法之間的差異，同時，為了規避實驗數據較少給結果帶來的偶然性，另外選取三幅圖片進行實驗，所得結果如表1所示.

表1 不同算法客觀評價

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

3 結 論

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