岩石薄片正交偏光融合图像的颗粒分割方法*

钟 逸，熊淑华，滕奇志，吴耀坤，龚 剑

(1.四川大学 电子信息学院 图像信息研究所，四川 成都 610065； 2.成都西图科技有限公司，四川 成都 610065)

0 引言

在石油地质行业中，对岩石薄片中的矿物颗粒进行粒度分析具有重大意义。利用岩石薄片在正交偏光下的图像进行分析已取得一定成果[1-2]。目前主要采用的方法有基于统计区域合并(SRM)[1]的算法、基于边缘流[2]的算法以及基于灰度阈值的分割算法[3]等。由于实际应用中岩石薄片颗粒的存在形式是复杂多变的，目前所采用的图像分割算法仍不能对工业矿场中采集得到的岩石薄片颗粒进行完整的分割，因此，如何准确地分割出岩石薄片中的矿物颗粒是行业中待解决的一大难题。

超像素分割算法[4]是近年来图像分割领域最热门的方法之一，它可以获取图像的冗余信息[4]，处理图像结果准确、算法复杂度低。本文提出了一种新的分割方法，首先采用基于Fast-MBD变换的超像素算法[5-7]对正交偏光序列最大值融合图像进行初步分割，结合矿物颗粒在正交偏光融合图像中的亮度特性[8]进行筛选，再以快速区域合并算法[9-10]的思想对过分割的矿物颗粒进行合并，实现了对岩石薄片图像中矿物颗粒更为完整的分割提取。

1 岩石薄片正交偏光融合图像

在单偏光下，由于矿物晶体的颜色相似、边界模糊、粘连性等原因，导致利用图像处理技术难以提取颗粒边界。非均质矿物在正交偏光下呈现一系列消光特性[8]，如消光性、干涉色等。本文所用的岩石薄片正交偏光序列图像由笔者所在课题实验组自主研制的偏光图像采集与拼接系统采集获得，首先在单偏光镜下采集一张图像，然后保持摄像头和显微镜载物台固定，切换到正交偏光镜下，每转过15°采集一张图像，共旋转120°(一个消光周期为90°)，总共采集9张相同视域下的正交偏光图像，图1所示为部分角度正交偏光下的序列图像，从图中看到，矿物颗粒在某些角度正交偏光下呈现出明显的边缘。

图1 偏光序列图像示意图

为将序列图中的所有矿物颗粒同时呈现，依据式(1)进行序列图像最大值融合，形成一张拥有视域中全部颗粒的正交偏光最大值融合图，如图2所示。

(1)

其中，P(x,y)为融合后的结果图，fi(x,y)是序列图中的第i张正交偏光图像。

图2 正交偏光融合图像

2 基于正交偏光融合图像的超像素分割

超像素是指受一定相似性准则约束的分割区域[4]，超像素分割实质上就是依据各种特征将相似或相同的像素聚集到一起的聚类问题，其最大的优势在于对分割目标的边缘拟合精准。超像素分割算法在图像分割领域已广泛应用并有一些改进，主要分为两类：基于图论的方法和基于梯度上升的方法，本文采用基于梯度上升的Fast-MBD变换超像素分割算法[5]进行图像的初始分割。

2.1 基于Fast-MBD变换的超像素分割

Minimun Barrier Distance (MBD)变换对像素值波动具有很强的稳定性，因此可以有效地应用于原始像素，无需区域抽象。对于一幅二维单通道数字图像I和种子集S，图像I中的一条路径p={p(0),p(1),…,p(k)}是由一系列像素所构成的序列，而序列当中任意相邻的像素对在图像I中是相邻的。距离变换的目的在于计算距离图D，这样对于像素n，其所对应的距离图表示为：

(2)

其中，∏S,n是像素n到种子像素集S的所有路径集合，F(p)是路径成本函数，其定义取决于具体应用。

文献[5]中引入的MBD变换被证明在图像分割中是有利的，其对噪声和模糊有很强的鲁棒性。对于一幅三通道彩色图像I，其各通道分量表示为I1,I2,I3，定义路径成本函数如下：

(3)

其中d(p0,pn)表示路径的两个端点p0和pn之间的空间欧氏距离，α表示空间距离相对于沿路径Barrier Distance的权重参数[6]。相关研究结果表明，MBD在图像分割上对噪声和模糊有很强的鲁棒性。然而，MBD变换的计算效率使其不适于实际应用，因此，文献[5]中使用的Fast-MBD变换算法，其采用正向光栅扫描和反向光栅扫描搜索访问每个像素n，如图3所示，并运用像素n的半邻域对应的相邻像素m迭代更新像素n的最小路径成本D(n)，更新公式如下：

(4)

F(p(m)·[m,n])=max{U(m),I(n)}-min

{L(m),I(n)}

(5)

其中，p(m)表示像素m到种子点集合的路径。[m,n]表示从m到n的边,p(m)·[m,n]表示将边[m,n]附加到p(m)至像素点n所形成的路径，也就是像素n到种子点的路径。U(m)和L(m)分别表示路径p(m)中的最大像素灰度值和最小像素灰度值。因此，通过使用两个辅助映射U和L，可以有效地计算新的F(p(m)·[m,n])。当F(p(m)·[m,n])

(6)

图3 光栅扫描搜索示意

2.2 正交偏光融合图像的Fast-MBD变换的超像素分割

由于不同的矿物颗粒在正交偏光下呈现出复杂的光学特性，在颗粒表面会形成不同的纹理特征，同时，正交偏光下一些碎屑和填隙物在图像中类似“噪声”，这些“噪声”在颗粒分割时形成高频信息。为了减弱此类“噪声”高频信息的影响，对融合图像的滤波预处理是有必要的。双边滤波[11]是一种非线性滤波方法，同时将图像的空间信息和像素值的相似性纳入考虑范围，以达到保留边界信息并去噪的目的。

拉普拉斯算子[12]是一种二阶边缘检测方法，用来提取图像边缘，对于低纹理区域进行完全平滑，将边缘与原图相加，就可以达到锐化的效果。本文采用双边滤波结合拉普拉斯算子对图2的融合图像进行预处理，滤波结果如图4(a)所示。对预处理后的融合图像采用2.1节中的算法进行图像分割，其分割效果如图4(b)所示。

图4 滤波预处理及Fast-MBD超像素分割

3 矿物颗粒目标的过分割合并

从图4(b)的结果可以看出，利用超像素分割算法可以完整地保留矿物颗粒的边缘信息，但是矿物颗粒会出现过分割的现象，同时，矿物颗粒目标以外的背景区域也会有过分割现象。首先需要将过分割的矿物颗粒筛选出来，再对颗粒目标进行合并。由于同一颗粒颜色一致且均匀，因此，利用目标区域正交偏光融合图像和单偏光图像的颜色相似性合并同一颗粒目标的过分割区域是可行的。本文结合文献[9]中提出的基于RAG(区域邻接图)和NNG(最近邻接图)数据结构的快速区域合并算法的思想对颗粒进行合并。

3.1 矿物颗粒目标筛选

采用基于Fast-MBD的超像素算法得到了很好的颗粒边界，对于正交偏光下部分碎屑岩颗粒的边划分界也比较理想。利用矿物颗粒在融合图中呈现的亮度特性，在正交偏光融合图中采用以下步骤筛选获取过分割的颗粒目标图形层：

(1)对2.2节过分割后的边缘结果使用目标反转操作得到初步的全部目标区域。

(2)通过对正交偏光融合图各颜色通道进行直方图统计，计算各颜色分量的最大类间方差(Otsu)[13]值，依据式(7)对各分量Otsu值取不同的权重，最终得到融合图的Otsu：

(7)

(3)对步骤(1)中取得的所有区域进行颗粒目标筛选，计算每个目标的平均灰度值，以式(8)为基准将筛选的非目标区域置为背景色：

(8)

式中，Iobj为目标区域在正交偏光的平均灰度值；ξ为亮度筛选阈值系数，且0<ξ<1。

(4)对步骤(3)中得到的颗粒目标区域图形层进行去噪、填孔、边界平滑等一系列数学形态学[14-15]操作，得到过分割的矿物颗粒图层，如图5(a)所示。

3.2 矿物颗粒目标合并

3.2.1 相似性合并准则的确定

(1)融合图像的颜色直方图。颜色直方图包含了正交偏光融合图中矿物颗粒完整的像素数据信息，反映了该图像像素值的分布情况，呈现矿物颗粒各方面的特征。在此利用矿物颗粒区域的颜色直方图来作为图像相似度的判定准则。

对于直方图的相似度比较，本文采用卡方比较(Chi-Square)。卡方比较来源于统计学中的卡方检验，卡方检验就是统计样本的实际值与推断值之间的偏离程度，计算公式如式(9)所示：

(9)

其中，p(x)、q(x)为两组不同的概率分布，且0<χ2<1。

根据式(9)，定义融合图像颜色直方图的相似度合并准则卡方系数χC2为：

(10)

其中，H1、H2分别为两个相邻区域所对应在融合图中的颜色直方图，L为图像的灰度级。在计算之前已对两区域的颜色直方图H1、H2做了归一化处理，其归一化公式如式(11)所示：

(11)

(2)单偏光图的颜色直方图。在单偏光图像中，由于同一矿物颗粒所在区域的RGB值差异很小，因此，将区域的颜色直方图纳入相似度的判定准则，根据式(9)可得单偏光图下颜色相似性合并准则χA2为：

(12)

式中，H1、H2分别为两个相邻区域所对应在单偏光图中的颜色直方图，同样，H1、H2在计算前已根据公式(11)做了归一化处理。

综合融合图颜色直方图以及单偏光图颜色直方图信息，定义本文方法中使用的相似性合并准则如式(13)所示：

χ2=βχC2+(1-β)χA2

(13)

其中，β为权重系数，0<β<1。经大量实验，β取值为0.68时合并效果最佳。

3.2.2 目标区域合并算法

区域邻接图(RAG)数据结构可以表示颗粒目标筛选后所有目标之间的邻近关系，采用该数据结构的思路判断所有目标的邻接区域，并依据3.2.1节中讨论的直方图相似性准则实现颗粒目标的合并，同时，采用最近邻接图(NNG)优化颗粒目标合并算法的速度。具体操作过程如下：

(1)初始化数据结构：计算所有目标区域下对应融合图的颜色直方图以及单偏光图像的颜色直方图。将全部的目标区域作为单独的节点，结合区域邻接图(RAG)的思路生成目标区域邻接映射图，根据此图生成对应的最近邻接图(NNG)。

(2)区域合并：依据3.2.1节中的直方图相似性准则计算两邻接区域的相似性，若达到阈值，则将两邻接区域合并为一个区域。

(3)更新数据结构：合并之后按步骤(1)更新计算目标区域下对应融合图的颜色直方图以及单偏光图像的颜色直方图，并更新目标区域邻接映射图和NNG。

(4)重复步骤(2)～(3)，直至没有目标区域需要合并为止。

图5(b)、(c)为颗粒目标合并前后的局部对比图。由对比图看出，该方法可以实现过分割颗粒目标的合并，最终达到岩石薄片中矿物颗粒较为理想的分割提取。

图5 矿物颗粒目标筛选及合并效果

4 实验结果分析

本文基于Fast-MBD超像素算法，针对矿物颗粒在正交偏光融合图像中的亮度特性并结合单偏光图像完成矿物颗粒图层的提取。

经过对大量岩石薄片图像的测试，结果显示，本文的分割算法能够取得较好的矿物颗粒分割效果。以图1中岩石薄片的正交偏光序列图像为例，将本文方法与文献[1]、文献[2]中的矿物颗粒分割算法的效果进行对比分析，如图6所示。其中，图6(b)为本文算法得到的矿物颗粒目标图层，图6(c)、图6(d) 分别为文献[1]中使用的SRM算法和文献[2]中使用的边缘流算法得到的矿物颗粒目标图层。

图6 本文算法与其他算法分割结果图对比

根据该薄片的原始偏光序列图像(图1)以及其最大值融合图像(图6(a))可以看出，文献[1]中使用SRM算法得到的结果容易出现欠分割的现象，而且由于岩石薄片图像的复杂性，分割结果中一些矿物颗粒的边界贴合不太准确；文献[2]中采用边缘流算法得到的分割结果对于部分完整的石英颗粒而言比较理想，但由于长石、岩屑等颗粒表面纹理的复杂性，而边缘流算法对图像高频信息敏感，使得此类颗粒容易出现过分割现象，从而导致提取颗粒的不够完整。本文提出的算法可以较为完整地分割提取的矿物颗粒，优化了表面纹理复杂的长石、岩屑颗粒等矿物颗粒的过分割现象。

5 结论

本文利用同一视域下岩石薄片的单偏光图像和正交偏光序列图像，提出一种新的矿物颗粒分割方法：首先将岩石薄片在偏光采集系统中获得的序列图像进行最大值融合，对经双边滤波结合拉普拉斯算子锐化预处理后的融合图像再采用基于Fast-MBD的超像素算法进行矿物颗粒分割，然后针对矿物颗粒在正交偏光融合图下的亮度特性对颗粒目标区域筛选，最后结合目标区域在融合图下的颜色直方图以及单偏光图下的颜色直方图作为相似性特征，对筛选后的过分割矿物颗粒进行合并，得到最终矿物颗粒分割的结果。

实验结果表明，本文算法结合了同一视域下正交偏光融合图中矿物颗粒的亮度特性以及对应单偏光图的信息，能够较为精确地分割提取出岩石薄片中的矿物颗粒目标，得到较好的分割结果。