駱簡,滕奇志,何海波
(1.四川大學電子信息學院圖像信息研究所,成都610065;2.成都西圖科技有限公司,成都610065)
在石油地質勘探行業中,通常使用巖石樣本制成巖石薄片研究地質油氣分布情況。巖石薄片在長期保存中容易受到損壞,而顯微成像技術為巖石薄片的數字化信息保存提供了可能,對于后續巖石薄片的研究具有重要意義。由于顯微鏡視域的限制,小視域下無法觀察巖石薄片的整體顆粒分布情況,大視域下無法看清巖石薄片的細節紋理[1]。為了從整體上對巖石薄片進行分析研究,需要在多個連續的視域下采集圖像并進行拼接。在圖像采集過程中,由于自動聚焦存在誤判現象、或因載物臺的抖動以及巖石薄片的不平整,導致采集到的某個視域圖像模糊或局部區域發虛,最終造成拼接好的全薄片圖出現區塊狀模糊,不能較好地保留巖石薄片信息。為了較好地保留巖石薄片信息,有必要在進行巖石圖像拼接之前對模糊圖進行復原增強處理。
本文在USM(Unsharp Mask)算法[2]的基礎上進行改進,通過Canny 邊緣算子提取擁有較少噪聲的真實邊緣,在Lab 彩色空間[3]的L 通道對模糊圖片進行邊緣銳化處理,利用動態增益參數抑制噪聲的同時達到局部區域圖像增強[4]。
Lab 色彩空間是基于CIE 色彩測量系統的一種色彩模式,其本身基于人類視覺,是獨立于硬件設備的,它彌補了RGB 模式依賴于硬件設備特性的缺陷。自然界中任何一點色彩都可以在Lab 空間中表達出來,RGB 所能描述的色彩都能映射到Lab 空間,Lab 比RGB 更適合進行數字圖像處理。Lab 空間基于亮度通道L、紅綠通道a 和黃藍通道b 三個通道,各個通道的相關性較弱[3],有利于后續巖石圖像在L 通道的增強處理,保證了巖石薄片的色相不會受到影響。因此,我們將巖石圖像從RGB 色彩空間映射到Lab 色彩空間進行處理。而RGB 空間無法直接轉換成Lab 空間,需要先轉換成XYZ 空間再轉換成Lab 空間。
(1)RGB 空間轉換到XYZ 空間[5]:

則:

(2)XYZ 空間轉換到Lab 空間[6]:

(3)轉化Lab 的范圍:我們處理的是8 位圖,0≤L ≤100,-128 ≤a ≤127,-128 ≤b ≤127,式(7)將其轉換到0-255 范圍內處理:

Canny 邊緣檢測算子非常優秀,它能夠準確地捕獲圖像中盡可能多的邊緣并且盡可能地濾除圖像的噪聲產生的假邊緣,同時能夠準確定位在檢測所得邊緣的中心,并且只返回一個邊緣響應點。因此,Canny 邊緣檢測具有低差錯、易定位、單一響應等特點。Canny 邊緣檢測算法[2]由下列基本步驟組成:
(1)高斯平滑輸入圖像
輸入圖像為s(x,y),高斯函數G(x,y)由式(8)表示:

通過式(9)得到高斯平滑后的圖像B(x,y),其中(*)表示卷積運算:

(2)計算平滑圖像的梯度幅值和角度
梯度的幅度和角度能反映出每一點處的邊緣強度和方向,結合式(10)和式(11)得到幅值M(x,y)和邊緣(法向)方向α(x,y):

(3)對圖像的梯度幅值應用非最大抑制
M(x,y)是使用梯度產生的,其局部最大值周圍通常包含更寬的范圍,需要使用非最大抑制進行邊緣細化。非最大抑制方案詳細如下:
①對于通過3×3 領域中間點的邊緣,可以劃分四個方向:水平d1、垂直d2、+45 度d3、-45 度d4。如表1所示,如果邊緣法向方向的范圍是從67.5°到112.5°,或從247.5°到292.5°,則我們稱該邊緣方向為垂直邊緣方向;其他方向依此類推。
②尋找最接近α(x,y)的方向dk。
③如果M(x,y)的值同時大于沿dk方向的領域梯度幅度,則令gN( )x,y=M(x,y);否則,抑制非最大邊緣點的M(x,y),為非最大邊緣抑制后的圖像梯度幅值。
(4)閾值處理和連接邊緣。
Canny 算子采用雙閾值對gN進行閾值處理,高閾值TH和低閾值比例一般為2:1 或3:1。

通過式(12)獲取強邊緣,(13)獲取弱邊緣。閾值處理后,gNH中的高像素被認為有效邊緣像素并對其進行標記。閾值TH決定著gNH邊緣縫隙的大小。連續的邊緣形成步驟如下:
①在gNH中定位下一個未到達的邊緣像素p。
②在gNL將所有弱像素記為有效,用8 連通法連接到p。
③如果gNH中的所有非零像素已到達,則進行第4步,否則轉到第1 步。
④將gNL中其他未標記像素認為無效,置零。
最后,將gNL中所有非零像素附加到gNH中,用Canny 算子形成最終的輸出邊緣圖像。
傳統的USM 算法主要用于印刷和出版界[2],先將原圖低通濾波,然后確定原圖和低通濾波圖的差異,并利用尺度因子放大補償該差異以增加清晰度。USM 算法詳細描述如下:
(1)先將原圖s(x,y)卷積一個低通函數w(x,y)得到模糊圖blur(x,y):

表1 邊緣方向的劃分

(2)將原圖與低通圖作差就得到非掩蔽模板m(x,y):

非掩蔽模板包括圖像的細節邊緣,噪聲以及過沖與下沖點等信息。
(3)將上述模板乘以尺度因子k 加到原圖上得到銳化圖d(x,y):

根據線性系統卷積理論[2],并充分利用模板的對稱性,式(15)可簡化為加權求和的過程,即卷積運算滿足下式:

尺度因子k 控制著圖像銳化的程度。當k=1 時,對原圖像的處理稱為非銳化掩蔽;當k >1時,對原圖像的處理稱為高提升濾波;選擇k <1 則不強調USM 模板的貢獻[7]。
傳統的USM 算法是通過原圖減去低通濾波圖得到需要增強的全局像素區域,增強系數為一固定值k,在增強圖像邊緣與細節的同時噪聲也得到同樣倍數的放大;當k 值較大時,增強后的圖像會出現負灰度值和超過255 灰度值的像素點,圖像局部將被過增強而出現暈輪現象[7-8]。
文獻[10]中先計算出3×3 像素領域的局部方差值v,然后對中心像素點進行分類。設置兩個正的閾值τ1、τ2和三個常數增益系數1、αdh、αdl,其中τ1<τ2,1 <αdl<αdh。如果v <τ1,將中心像素歸為像素平滑區域,將增益系數設為常數1;如果τ1≤v <τ2,將中心像素歸為中等比對度區域,將增益系數設為常數αdh;如果v ≥τ2,將中心像素歸為像素劇變區域,將增益系數設為常數αdl。此算法中,雖然考慮到局部區域的像素變化情況,但是并沒有消除大部分噪聲增強的影響;閾值和增益系數參數選擇對圖像的增強也有較大影響。此外,對高分辨率的彩色圖像進行處理,需要采用更大的模板和局部區域動態增益。
為了克服傳統USM 算法和文獻[10]算法的缺點,我們對USM 算法進行改進,采用Canny 邊緣算子,在抑制噪聲的同時選擇性地提取出真實有效的邊緣像素[9];此外,采用動態增益系數來控制高頻和細節部分增強的程度,滿足圖像局部區域動態增強[10]??紤]到圖像局部區域的差異性可以用局部標準差LSD 來表征,我們將全局的標準差GSD 與LSD 的比值作為增益系數,保留尺度因子K(本文取值0.24)控制對比度拉伸的程度。改進后的USM 圖像銳化的詳細步驟如下:
(1)選取大小為(2n+1)×(2n+1)的像素區域,通過式(19)計算局部區域的像素平均值m(i,j),通過區域中間像素f(i,j)的窗口滑動,最終獲取原圖的模糊圖。

(2)采用兩次Canny 邊緣算子,適當設置閾值TH和TL,分別提取出原圖和模糊圖的邊緣,分別記為E1(x,y),E2(x,y),從而獲取需要增強的像素區域E(x,y),即:

(3)計算(2n+1)×(2n+1)局部區域的標準差σ(i,j)和全局區域的標準差GSD,獲得動態增益系數矩陣DG。

(4)計算增強后的圖像d(x,y),將負灰度值像素置為零,超過255 灰度值像素置為255:

(5)增益系數矩陣DG 由全局區域標準差和局部區域標準差的比值計算得出,反映了局部區域像素與全局區域像素對比度的差異。在圖像的邊緣或者其他變化劇烈的像素區域,局部均方差比較大,DG 的值比較小才不會導致圖像過增強;在圖像灰度變化緩慢的區域,局部均方差就會很小,DG 值過大會在增強圖像的同時引起噪聲的放大,因此需要給DG(i,j)設置閾值。經過大量的實驗,閾值為11 時多數圖像增強的效果比較好。為了防止緩慢變化區域過度增強而引起噪聲 的 放 大,當DG(i,j)>11 時 對 其 取 值 為1,當DG(i,j)<=11 時對其取值不變。對于少數增強效果不好的圖像,可以調節DG 的閾值,直到對圖像增強效果滿意為止。
為了衡量圖像的處理效果,可以利用一些客觀分析函數作為評價指標。常見的客觀評價指標主要有均方誤差(MSE)、峰值信噪比(PSNR)、結構相似性(SSIM)函數[11-12]。
均方誤差是計算已處理圖和參考圖像素差值的均方值,然后通過MSE 值來衡量已處理圖與參考圖的接近程度,MSE 值越接近于0,表示已處理圖越接近于參考圖像。MSE 值可通過式(24)獲得:

信噪比函數PSNR 一般是用于信號最大值和背景噪聲之間的衡量,可以用來評價已處理圖像的效果,其值可以結合式(24)、式(25)獲得:

其中,d(i,j)為待評價圖像,s(i,j)為參考圖像,圖像的大小是M×N,n 代表處理圖像的位深,此處為8。PSNR 的單位為db,值越大代表增強處理后的圖像效果越好。
結構相似性SSIM 函數也是一種圖像質量評價標準,分別從亮度L(X,Y)、對比度C(X,Y)、結構性S(X,Y)三方面度量圖像處理后的質量。

其中,C1=6.5025,C2=58.5225,C3=29.2613,μY表示圖像X,Y 的均值,σX,σY表示圖像X,Y 的方差,σXY表示圖像X 和Y 的協方差:

SSIM 值域為[0,1],值越大,表示處理圖效果越好。SSIM 值通過式(32)獲得:

在實際應用中,將圖像分成N 塊,采用高斯加權計算每一塊的SSIM,最后取得平均值MSSIM(X,Y)。

為了檢驗本文改進算法的效果,采用Canon EOS 600D 相機連接顯微鏡拍攝了多組巖石薄片圖進行實驗;每一組圖都由一張清晰圖和3 張失焦模糊圖組成,其中清晰圖作為參考圖像,失焦模糊圖作為待處理圖像。然后采用傳統USM 算法、文獻[10]中改進的USM算法、以及本文USM 改進算法對采集的失焦模糊圖進行增強處理,最后結合參考圖像進行對比分析。本文分別在單偏光和正交偏光模式下處理分辨率為5184×2912 的巖石圖像,實驗結果如圖1、圖2。
圖1(a)為單偏光模式下采集的巖石薄片某視域清晰圖,圖1(b)為與圖1(a)相同視域采集的失焦模糊圖。圖1(b)進行傳統USM 算法增強處理的圖,可以看出,與模糊圖相比清晰度有明顯的提升,但與原圖相比噪聲也明顯得到放大。圖1(d)為對圖1(b)進行文獻[10]改進USM 算法增強處理圖,由于考慮到局部區域的變化情況,清晰度比傳統USM 算法增強圖有所提高,對比度也得到較好的拉伸,但還是存在明顯的噪聲;圖1(e)為對圖1(b)進行本文改進USM 算法增強處理的圖,清晰度比傳統USM 算法處理圖有較大提升,與文獻[10]算法處理圖對比清晰度變化不大,但卻很好地抑制了噪聲。

圖1 單偏光圖像處理結果

圖2 正交偏光圖像處理結果
圖2 (a)為正交偏光模式下采集的巖石薄片某視域清晰圖,圖2(b)為與圖2(a)相同視域采集的失焦模糊圖。圖2(c)為對圖2(b)進行傳統USM 算法增強處理的圖,與模糊圖相比清晰度有一定的提升,但與原圖相比噪聲也明顯得到放大,表現為許多單獨的噪點、小的巖石顆粒粘連和部分巖石顆粒邊緣出現暈輪。圖2(d)為對圖2(b)進行文獻[10]改進USM 算法增強處理圖,與失焦模糊圖和傳統USM 算法增強圖相比清晰度得到更大提高,對比度得到更大拉伸,但噪聲也得到明顯放大。圖2(e)為對圖2(b)進行本文改進USM 算法增強處理的圖,與模糊圖相比清晰度具有明顯的提升,與傳統USM 增強圖相比清晰度也有較大提升,暈輪現象得以消除,巖石顆粒粘連性減弱,巖石顆粒邊緣的暈輪被消除,對比度也略微拉伸;與文獻[10]算法處理圖對比清晰度變化不大,但噪聲卻明顯得到有效地抑制。
以采集的清晰圖為參考圖,分別計算出失焦模糊圖、本文算法增強圖、傳統USM 算法增強圖、文獻[10]算法增強圖相對于參考圖的MSE 值、PSNR 值和MSSIM 值。表2 和表3 為上述兩組實驗的客觀評價函數值的統計結果。

表2 單偏光圖處理結果客觀評價

表3 正交偏光圖處理結果客觀評價
從表中可以看出如下大小關系,MSE 值:模糊圖>傳統USM 算法處理圖>文獻[10]算法處理圖>本文算法;PSNR 值:模糊圖<傳統USM 算法處理圖<文獻[10]算法處理圖<本文算法;MSSIM 值:模糊圖<傳統USM算法處理圖<文獻[10]算法處理圖<本文算法。實驗結果表明,本文改進USM 算法在提升圖像清晰度的同時能有效地抑制噪聲增強,處理圖效果更接近于參考圖,這與主觀分析結果一致。
本文提出一個針對巖石薄片圖像的銳化算法,通過Canny 邊緣算子提取出需要增強的有效邊緣部分,采用基于局部標準差的局部區域動態增益函數DG,給DG 設置閾值控制圖像增強的范圍。實驗結果表明,本文算法對模糊巖石薄片有較好的增強效果,可以克服傳統USM 算法增強邊緣的同時放大噪聲的缺點,消除大部分的暈輪現象;克服文獻[10]算法不能有效抑制噪聲的缺點;可以較好地保留整個巖石薄片的紋理信息以及提高全薄片圖的質量。