巖礦薄片偏光序列圖像的礦物顆粒分割方法①

吳耀坤,滕奇志,何海波,彭志偉

1(四川大學 電子信息學院 圖像信息研究所,成都 610065)

2(成都西圖科技有限公司,成都 610024)

圖像處理技術在巖礦薄片鑒定分析中的應用已經成為地質行業的常用手段.礦物顆粒分割是薄片鑒定和礦物成分識別的基礎,其效果直接影響后續鑒定分析的準確性.利用巖礦薄片在固定視域下的正交偏光序列圖像進行礦物顆粒的提取是一種新的研究方法,已取得一定成果.在實際生產應用中主要采用基于灰度閾值的分割算法[1]、基于統計區域合并(SRM)的算法[2]以及基于邊緣流[3]的算法.

超像素分割算法[4]的最大優勢在于邊緣精準定位以及算法復雜度低.本文使用基于熵率的超像素算法[5,6]對偏光序列融合圖像進行初步分割,通過快速區域合并算法[7,8]合并具有相似性的過分割區域,同時,基于礦物顆粒在正交偏光序列圖像中的變化規律,對結果進行二次合并,改善超像素算法過分割嚴重的問題,實現了對巖礦薄片圖像中顆粒的分割提取.本文分割算法的流程圖如圖1所示.

圖1 本文算法流程圖

1 基于偏光序列圖像的超像素分割

超像素[9]是指具備相似顏色、灰度、紋理等特征的相鄰像素組成的圖像塊,超像素分割實質上就是依據各種特征將相似或相同的像素聚集到一起的聚類問題.超像素分割算法在國內外已廣泛應用,主要分為基于圖論和梯度下降的兩類方法,本文采用基于圖論的熵率超像素分割算法對巖礦薄片偏光序列圖進行礦物顆粒分割.

1.1 熵率超像素分割

熵率法超像素[5]提出了將圖上隨機游走熵率和平衡項相結合的目標函數,并在聚類過程中不斷優化目標函數實現對圖像的分割.該方法分割出的超像素塊相對規則、均勻,其主要思想分為以下四個方面.

(1)圖的構造:采取映射的方式,將圖片映射至無向圖G=(V,E),其中V表示像素點集,E為邊集,并將像素間的相似性定義為邊的權值 ω .目標是尋找到集合A?E,使生成的無向圖G=(V,A)中包括K個連通的子圖.

(2)隨機游走熵率:熵率用來度量隨機過程的不確定性,可以很好的表征信息的冗余程度、內部結構的緊湊性,隨機游走模型的加入能夠獲取緊湊、均勻的超像素區域塊.無向圖G=(V,A)上的隨機游走熵率計算公式如下:

(3)平衡條件:在目標函數中添加平衡條件,可以約束限制熵率的隨機性.令A為 已選擇邊集,劃分出A中超像素塊集為為G=(V,A)中連通子集的數目,ZA為子圖聚類的概率分布:

由此,定義平衡函數為:

其中,ZA的熵H(ZA)可以使聚類結果有相似的尺寸,同時NA約束集群的數目.圖2中用不同粗細的線條連接代表不同的集群,可見圖2(a)中的平衡項大,集群分布更均勻.

圖2 平衡條件對集群的作用

(4)目標函數:組合隨機游走熵率和平衡條件,得到聚類最優化的目標函數如式(5)所示,式中集合滿足A?E,且 λ為正實數.聚類過程中通過最大化目標函數得到圖像的最佳分割結果.

1.2 基于偏光序列圖像的熵率超像素分割

巖礦薄片在單偏光下獲取的圖像中,礦物顆粒之間粘連比較嚴重、缺少清晰的邊界,當顆粒顏色相差不大時,很難辨認出具體是幾個顆粒,因此,給顆粒分割帶來極大的困難.在正交偏光鏡下,礦物顆粒在視域下會表現出一系列特有的光學特性[10],如消光特性、干涉色等.正交偏光下的消光特性[10]在非均質礦物中普遍存在,當正交偏光鏡旋轉一周,礦物共會出現四次消光現象,并且,不同礦物顆粒在不同偏光角度下消光狀態也可能不一樣,加上礦物顆粒在正交偏光下的干涉色,由此顯現出明顯的邊緣.本文所需的序列圖像由自主研制的偏光圖像采集與拼接系統獲取,先在單偏光下采集視域的一張圖像,然后保持顯微鏡載物臺固定,每旋轉 10°正交偏光鏡采集一張圖像,共旋轉 120°,共采集13張固定視域下的正交偏光圖像,圖3為部分序列圖像.

圖3 偏光序列圖像示意圖

考慮到礦物顆粒本身的復雜性,以及單個正交偏光角度下礦物顆粒的不完整性給分割帶來的影響,首先選取一個消光周期(90°)內的正交偏光序列圖像,即拍攝的前十張正交偏光序列圖像,依據式(6)進行圖像差分[11]融合,形成一張擁有視域中全部顆粒的正交偏光圖像,如圖4(a)所示.

其中,g(x,y)為融合后圖像,fi(x,y)是序列圖中第i張正交偏光圖像,且g(x,y)的R、G、B分量gR(x,y)、gG(x,y)、gB(x,y)均已映射至 0 ～255.

為了減弱R、G、B分量間的相互影響,提高礦物顆粒分割的效果,本文將融合后的圖像g(x,y)轉化至Lab空間.對Lab空間下的融合圖像采用1.1節中的算法進行圖像分割,其分割結果如圖4(b)所示.

圖4 熵率超像素分割

2 礦物顆粒目標的區域合并

2.1 基于融合圖像與單偏光圖像的快速區域合并

由圖4(b)的結果可以看出,利用超像素分割算法可以完整的保留礦物顆粒的邊緣信息,但是使礦物顆粒出現過分割現象.一般而言,屬于同一個礦物顆粒的區域具有相似性,利用其相似性將同質區域合并.本文首先采用文獻[7]中提出的基于RAG(區域鄰接圖)和NNG(最近鄰接圖)數據結構的快速區域合并算法對同質區域進行合并,其中合并判定準則結合了目標區域融合圖像的差分直方圖和單偏光圖像的灰度直方圖信息.

(1)相似性合并準則的確定

① 融合圖像的差分直方圖:差分直方圖反映了正交偏光下礦物顆粒圖像在一定方向、一定距離上相鄰點灰度之間差的概率分布[12,13].在正交偏光的融合圖像中,礦物顆粒本身的形態、紋理比較復雜,較顏色直方圖而言,差分直方圖具有顆粒目標區域圖像的綜合紋理信息,因此,提取偏光序列融合圖像的差分直方圖作為相似度主要判定準則.

設一幅具有M(M=1,2,···,L)灰度級的圖像F,區域D內的兩個像素f1=f(x,y)、f2=f(x+dx,y+dy),(dx,dy)∈D,對于偏移(dx,dy)的像素差分按式(7)計算如下:

由此可得歸一化的差分直方圖如下:

式中,hd(i)為區域D中d(x,y)的分布,N為區域D的像素總數.在計算每一個區域的差分直方圖矩陣時,本文在每個像素的8鄰域中選取四個方向(dx,dy)=(d,0),(d,d),(0,d),(-d,d)進行合并計算.

統計學中,可以用Bhattacharyya系數[14]計算兩個統計樣本的重疊量,從而進行兩組樣本的相關性測量,其計算公式如式(9)所示.在進行直方圖的相似性計算時,該計算方法可以得到最好效果.

式中,p(x)、q(x)為兩組不同的概率分布,且0＜BC＜1 .

因此,根據式(9),定義融合圖像差分直方圖的相似度合并準則hC為:

其中,Pdp、Pdq分別為兩不同區域p、q的差分直方圖.

② 單偏光灰度直方圖:單偏光圖像中,同一礦物顆粒所在區域的灰度值差異變化不大,可將區域的灰度值作為次要判定準則.直方圖統計是提取區域灰度特征的常用方法,根據式(9)可得單偏光下灰度相似性合并準則hA為:

其中,Hp、Hq分別為兩不同區域p、q的灰度直方圖.

綜合①、②定義本文方法中使用的相似性合并準則如式(12)所示:

(2)快速區域合并算法

區域鄰接圖(RAG)的數據結構可以表示圖像中全部區域之間的鄰近關系,利用該數據結構并依據(1)中討論的相似性準則實現區域合并,同時,使用最近鄰接圖(NNG)實現基于RAG區域合并算法的加速優化.具體操作過程如下:

① 初始化數據結構:計算所有超像素區域下對應融合圖像的差分直方圖以及單偏光圖像的灰度直方圖.將全部的超像素區域作為單獨的節點,按照區域鄰接圖(RAG)的規則生成一個無向鄰接圖,并據此生成對應的最近鄰接圖(NNG);

② 搜索鄰接情況:采用基于NNG的搜索方法,對于每一個區域確定其所有的鄰接區域;

③ 合并:依據(1)中的相似性合并準則計算兩鄰接區域的相似性,若達到閾值,將兩鄰接區域合并成一個區域;

④ 更新數據結構:合并之后計算該區域相應的差分直方圖以及灰度直方圖數據,并更新RAG和NNG;

⑤ 重復②～④,直至沒有區域需要合并為止.

圖5為合并前后對比圖.可以看出,使用該方法基本可以實現大部分礦物顆粒的過分割合并,但由于礦物顆粒本身的復雜性,仍有部分礦物顆粒的過分割區域未實現合并.本文在此結果基礎上,利用礦物顆粒在正交偏光序列圖像中的變化特性對未實現合并的礦物顆粒再次合并,最終達到較為理想的結果.

圖5 初步合并效果

2.2 基于正交偏光序列圖像的二次區域合并

2.2.1 顆粒目標圖層的初步提取

采用熵率超像素算法可以得到良好的顆粒邊界,其中碎屑巖顆粒的邊界也較為理想.由于最終需要得到顆粒目標圖形層進行粒度分析、成分識別,因此,利用礦物顆粒在序列圖中的亮度特性采用如下步驟初步獲取顆粒目標圖層:

① 對2.1節的邊緣結果進行目標反轉操作得到全部目標區域;

② 依據式(13)對①中的所有區域進行消光性篩選,將亮度差低于閾值T的區域置為背景色,獲得礦物顆粒區域;

其中,Imax、Imin分別為目標區域在正交偏光的前后序列圖中的最大亮度均值與最小亮度均值,T為亮度篩選閾值條件.經大量實驗,T取25篩選效果最佳;

③ 對②中的結果進行去噪、填孔、數字形態學、邊界平滑等操作.

根據上述步驟初步得到礦物顆粒目標的圖形層圖像如圖6所示.

圖6 亮度篩選結果

2.2.2 顆粒目標圖層的消光性合并

經過區域亮度篩選可基本獲取顆粒目標區域,但仍存在部分顆粒被分為若干個,因此,基于圖6的結果,利用礦物顆粒目標在序列圖中的消光變化特性繼續對仍存在的過分割顆粒目標區域進行相似性合并.

對于圖7(a)中的完整顆粒區域W以及其在圖7(b)中的過分割區域A、B,每個目標區域在正交偏光序列圖像下的亮度均值變化趨勢如圖8所示.

圖7 過分割顆粒示意圖

圖8 同一礦物顆粒的亮度變化趨勢

雖然礦物顆粒出現了過分割現象,由于其存在共同的消光特性,亮度均值在序列圖像中具有相同的變化趨勢.在統計學中,相關系數用來描述兩個序列之間的相關程度,設兩個亮度均值序列為a=(a1,a2,···,aN),b=(b1,b2,···,bN),則以式(14)計算序列的a、b的相關系數 βab:

其中,cov(a,b)為兩個序列的協方差,σa、σb分別為序列a、b的標準差.

因此,本文以兩區域在正交偏光序列圖下亮度均值序列的相關系數作為區域消光特性的合并準則,當兩區域的相關系數 β小于閾值 βT時,判定兩區域為同一顆粒目標,將兩目標區域進行合并.對初步提取的顆粒圖層中所有顆粒目標進行合并預測,直至無顆粒目標需要合并為止,得到最終的分割結果.

3 實驗結果分析

本文基于熵率超像素算法,針對礦物顆粒在正交偏光前后序列圖像中的變化特性完成礦物顆粒圖層的提取.本文實驗程序使用C++語言,在VS2015平臺下基于MFC框架編寫實現.

經過大量薄片圖像的測試,本文的分割算法能夠取得較好的分割效果.以圖3中礦物薄片的原始偏光序列圖像為例,將本文算法與文獻[2]、文獻[3]中的礦物顆粒分割算法的礦物顆粒提取結果進行對比分析,如圖9所示.其中,圖9(b)為本文算法得到的顆粒圖層,圖9(c)、圖9(d)分別為SRM和邊緣流算法得到的顆粒目標圖層.

圖9 本文算法與其他算法分割結果圖對比

結合該薄片的原始偏光序列圖像(圖3)以及原始融合圖像(圖9(a))可以看出,文獻[2]中利用SRM算法得到的分割結果容易出現大量的欠分割現象,使很多礦物顆粒目標粘連在一起,并且一些礦物顆粒的邊緣定位不夠準確;文獻[3]中利用邊緣流算法得到的分割結果大部分比較理想,由于邊緣流算法對邊界定位的敏感性,對于表面紋理特征較為復雜的堿性長石易出現過分割情況且不易融合,使提取區域不夠完整,如圖9(d)中標注出的顆粒目標.本文的分割算法可以較為完整地提取視域中的礦物顆粒,改善了堿性長石的過分割現象,同時,對于粒徑較大的碎屑巖顆粒也有了較為完整的提取結果,符合實際生產需求.

4 結論

本文針對固定視域下巖礦薄片的偏光序列圖像,提出一種礦物顆粒目標提取分割的方法:首先將巖礦薄片在正交偏光下的序列圖像進行差分融合后再基于熵率超像素算法分割礦物顆粒,然后提取分割區域在融合圖像下的差分直方圖以及單偏光圖像下的灰度直方圖作為相似性特征,對過分割區域進行快速區域合并,最后通過礦物顆粒在正交偏光下的消光特性對分割區域篩選后再次進行顆粒融合,從而得到最終的顆粒分割結果.

實驗結果表明,該算法綜合了視域下偏光序列圖的特征信息,能夠較為準確的提取礦物顆粒目標,并且,對碎屑巖顆粒的提取效果有一定的提升,為礦物顆粒的分割提取提供了一個新的方向.