基于空间标记的岩心孔隙并行分割

陈国军，李 胜

(中国石油大学(华东) 计算机科学与技术学院，山东 青岛 266580)

0 引 言

图像分割是指将图像分成若干具有相似性质的区域的过程，是许多图像处理任务的预处理步骤。现有的图像分割方法主要分以下几类：基于阈值的分割方法、基于区域的分割方法[1]、基于边缘的分割方法[2]以及基于特定理论[3]的分割方法等。

岩心图像中感兴趣的目标是孔隙和骨架，三维岩心模型建立的真实程度和三维分析结果依赖于岩心二维图像分割结果的好坏。目前针对岩心图像分割方法主要有两种，一种是选择其中的一帧图像根据实际孔隙度计算合适的阈值，并以该阈值对所有图像进行阈值分割，再根据实际孔隙度进行阈值优化；另一种是根据实际孔隙度对每一帧图像设定独立的阈值进行分割。文献[4-6]中提出了基于最大类间方差法(Otsu)[7]的不同改进方法，并取得了不错的效果；文献[8]提出了适用于体视图像岩粒分割的算法，采用两次分水岭并结合最小割的方法对岩石图像进行分割；徐永进等[9]提出了利用图像层间相关性实现CT图像半自动分割，原理是利用分割效果好的帧图像修复其余帧图像；文献[10]提出了一种基于模糊距离变换的改进分水岭[11-12]算法，通过对区域生长速度进行控制，优化了过分割问题；文献[13]利用了薄片视域的单偏光图与正交偏光序列图信息对演示薄片骨架孔隙进行分割；文献[14]采用聚类的方法对岩石CT图像分割及量化。不论是将一帧图像的分割阈值用于所有图像分割，还是每一帧图像采用不同的阈值进行分割，虽然也能取得不错的效果，但都只是考虑了岩心截面的二维信息，而忽略了岩心帧图像在三维空间的联系。为解决此问题，该文提出了基于空间标记的岩心孔隙并行分割方法，通过采用连通域标记的方法建立岩心图像在三维空间中的连通关系并实现自动分割。

1 图像预处理

Micro-CT扫描得到的岩心断层图像亮度较低，对于尺寸小于分辨率或尺寸在分辨率附近的特征体很难分辨，尤其是孔隙和骨架的边界区域，通常较为模糊，没有明显的分界，多为过渡灰度值，需对岩心图像进行对比度增强和滤波处理。直方图均衡化实质是图像增强的一种方式，经过直方图均衡化的图像对比度明显增强，使得图像有很强的清晰感。直方图均衡化处理的原理是把原始图像的灰度直方图从比较集中的某个灰度区间变成在全部灰度范围内的均匀分布，对图像进行非线性拉伸，重新分配图像像素值，使一定灰度范围内的像素数量大致相同，就是把给定图像的直方图分布改变成“均匀”分布直方图分布。中值滤波是非线性滤波，避免因线性滤波方法导致的边界模糊问题，并且不受奇异值的影响，可以保护图像细节，故在岩心滤波处理中被广泛采用。经过直方图均衡化和中值滤波后岩心图像骨架孔隙边界模糊问题得到改善，有利于进一步对图像进行处理。

图像二值化是岩心切片处理的中间步骤，通过设置阈值T将图像分为像素值为255或0的两类。岩心图像是在相同的环境中由CT机扫描而来，因此同组切片采用同一阈值分割。首先选取具有代表性的起始、中间和末尾三张切片采用Otsu计算初始分割阈值，并根据理论孔隙度优化分割阈值，最后以该阈值二值化所有切片图像得到包含孤立孔隙的二值图像。

图像二值化公式为：

(1)

式中，f(x,y)为当前坐标点的像素值，g(x,y)为二值化后该点的像素值，T为分割阈值，其中0≤T≤255。

2 二值图像连通域标记

二值图像连通域标记依据扫描图像次数分为一趟扫描算法[15]和两趟扫描算法[16]，一趟扫描算法主要有种子填充算法和游程标记算法；两趟扫描算法典型代表为基于快速连通域标记的Two-Pass算法。该文采用基于并查集的两遍扫描算法标记单帧图像连通域，通过结合多线程技术对多帧连通区域并行标记，使标记时间仅为串行标记时间的三分之一。

初始化与当前图像大小相同且全为0的标签矩阵，设(x,y)为当前处理的像素点，g(x,y)为该点像素值，l(x,y)为该点在标签矩阵中的标签值，label为标签矩阵当前的最大标签值，其中label≥2，基于快速连通域标记的Two-Pass算法的流程如下所述。

2.1 第一趟扫描

由于第一趟扫描时标签矩阵全部初始化为0，当前像素点后的像素点尚未处理，因此在处理时只需要考虑当前像素点左上方、上方、右上方和左边像素点的标签值，如图1所示，深色像素为当前像素点。

图1 八邻接像素点标签值示意图

扫描当前像素点(x,y)，如果g(x,y)=0，继续下一步，否则处理下一个像素点。

(1)如果l(x-1,y-1)=l(x-1,y)=l(x-1,y+1)=l(x,y-1)=0，即八邻接像素点标签值全为0，当前像素点和已被标记的像素点不邻接(新的连通区域)，则当前最大标签值加1，作为当前像素点的标签值，即：

label=label+1；l(x,y)=label

(2)如果l(x-1,y-1)=l(x-1,y)=l(x-1,y+1)=l(x,y-1)≠0，八邻接点像素值相同且不为0，即只与一个连通域邻接，则将该点的标签值设置为邻接点的标签值，表示如下：

l(x,y)=l(xi,yj),i,j=0,±1

(3)若l(x-1,y-1)、l(x-1,y)、l(x-1,y+1)、l(x,y-1)中至少存在两个相等且不为0的值，即八邻接点像素值不全相同，当前像素点与多个连通域相连，则将该点标签值设置为多个连通域标签值的最小值，内容如下所示：

l(x,y)=min{l(x-1,y-1),l(x-1,y),l(x-1,y+1),l(x,y-1)}

并将l(x-1,y-1)、l(x-1,y)、l(x-1,y+1)、l(x,y-1)与l(x,y)作为等价对存入映射集合Map中，供第二趟扫描使用，Map内容如下所示：

Map={l(x+i,y+j),l(x,y)|i,j=0,±1}

2.2 第二趟扫描

在第一趟扫描(3)步骤中，得出等价关系集合Map，根据Map中的等价关系，将当前像素点标签值更新为l(x,y)，即：

l(x+i,y+j)=l(x,y),i,j=0,±1

由于标签合并，导致标签号不连续，重新编号使标签号连续。扫描结束后，图像中具有相同标签号的像素就构成一个连通区域。

3 岩心图像空间标记

三维空间标记以二维连通域合并为基础，将多帧相邻图像连通域合并，形成多个空间连通体，根据岩心孔隙连通性规则，去除孤立的连通体(孤立孔隙)，由于岩心图像分辨率高且数量较多，采用多线程技术加速图像连通域合并。

3.1 标签合并

设定Si为上层帧图像，Si+1为下层帧图像，Si.l(x,y)为当前连通域标签值，下层帧图像对应像素点的连通域标签值为Si+1.l(x,y)，标签的合并采取小标签值为主原则即存在多个标签值时统一用最小的标签值替换其余标签值。

3.1.1 更新下层图像标签

在对单帧图像标记时，设定上层图像标签值小于下层标签值。如果Si.l(x,y)不为0(0为背景像素点直接跳过)且Si.l(x,y)所属的连通域只与上层图像的一个连通域相邻，则将Si+1.l(x,y)所在连通域的值替换为Si.l(x,y)的值，即：

Si+1.l(x,y)=Si.l(x,y)

若Si+1.l(x,y)所属的连通域与上层图像中的多个连通域相邻，则将Si+1.l(x,y)所属的连通域的标签值替换为上层多个连通域标签值的最小值，即：

Si+1.l(x,y)=Min{Si.l(xj,yj),…,Si.l(xk,yk)}

式中，Si.l(xj,yj),…,Si.l(xk,yk)为上层图像不同连通域的标签值，如图2(左)所示。

3.1.2 更新上层图像标签

若连通域Si+1.l(x,y)与上层图像的多个连通域相邻，即通过Si+1.l(x,y)使Si中原本不连通的两个或多个连通域连通，则合并这两个或多个连通域，用多个连通域的最小标签值代替其他标签值。

Si.l(xj,yj)=Si+1.l(x,y)

式中，Si+1.l(x,y)为下层连通域标签值，{Si.l(xj,yj)}为上层图像中与Si+1.l(x,y)相邻的连通域的标签值集合，用上层连通域中标签值的最小值替换其余连通域标签值。如图2(右)为更新后上下层图像连通域的标签值，更新后Si中标签值为3的连通域合并到标签值为2的连通域中。

图2 标签标记(左)与合并图(右)

3.2 去除孤立孔隙

数字岩心中有效孔隙为相互连通的，在一般压力条件下，可允许液体在其中流动的孔隙，除此之外为孤立孔隙(封闭孔隙)。即不属于起始和末尾帧图像且不与起始、末尾帧及多帧图像边界孔隙像素相连的孔隙为孤立孔隙。设Ubound为图像边界孔隙像素点的标签值和起始、末尾帧图像孔隙像素点的标签值集合，集合U为多帧图像连通体标签值集合，Uremove为孤立连通体标签值集合。

3.3 连通域合并优化

为简化图像标签合并，在对单帧图像标签时保证Si+1帧的标签值大于Si帧的标签值，Si帧标签的起始值为图像编号乘以大于所有单帧图像连通域个数的定值N，即Si×N。采用多线程并行合并技术，开始合并前先将所有帧图像按照层次分组，并设计数据结构SliceGroup存储合并所需的参数。

SlicesGroup:{slices,{G1.Si+slices,G2.Si}}

图像数量slices表明合并的层次，规模为2n(n=0,1,2…)，若共有2n帧图像，共需要进行n层合并。G1、G2为同一层中待合并的两组帧图像，G1.Si+slices为G1中的最后一帧图像，G2.Si为G2中的第一帧图像，其中G1.(Si+slices).l(x,y)

图3 合并前(左)与合并后(右)

4 阈值优化

去除孤立孔隙后，重新计算图像的孔隙度并与实验孔隙度对比，若误差超出限定范围，则调整阈值并二值化重复步骤2、3，直至误差缩小到限定范围以内。设f(x)为孔隙度计算函数，porosity为实验孔隙度，T为复杂度分割阈值。孔隙度f(x)计算公式为：

(2)

式中，AllHolePixel为总的孔隙像素个数，AllPixel为总的像素个数。阈值调整方法为：

(3)

5 实验分析

该文采用的实验平台为Intel(R) Xeon(R) CPU E5-1620，主频为3.70 GHz，系统内存为16 GB，显卡为NVIDIA Quadro K2000，显卡内存为2 GB，操作系统为Windows 10，实验开发环境为Visual Studio 2015。实验图像为650张分辨率为968*995吉林砂岩CT图像。

5.1 多层切片标记效率分析

为了验证基于多线程加速技术的有效性，任意选取不同编号的图像，采用八邻接搜索方式进行判别，分别测试普通Two-Pass扫描法、种子填充算法、基于并查集的Two-Pass扫描法的运行时间，并对650帧图像在多线程下标记时间与串行标记时间进行对比，对比结果如表1所示。

表1 图像标记效率对比

通过以上数据分析可知，三种算法在处理单张图像时相差不明显，基于并查集的Two-Pass算法速度最快，种子填充算法的速度最慢。对多帧图像进行标记时，通过多线程技术并行处理图像，时间为串行处理时间的三分之一左右。说明通过增加并查集和多线程技术，有效提高了图像的标记效率。

5.2 多层切片连通域优化分析

为了验证空间标记分割的有效性，选取同一组图像，根据孔隙度进行优化分割，选取多张图像对比分割效果如图4所示。通过对多帧图像连通域标记、合并、去除孤立孔隙，根据理论孔隙度优化分割阈值，孔隙度和阈值随着迭代次数增加变化如表2所示。

表2 孔隙度和阈值变化

由表2可知，未经过连通域优化(迭代次数为0)与经过连通域优化的图像虽然具有相同的分割阈值，但未优化的孔隙度大于优化后的孔隙度，表明非连通孔隙点对孔隙度有一定的影响。随着迭代次数增加分割阈值增大，孔隙度逐渐增大，但始终围绕着0.12上下波动，阈值在103～104之间变动，通过设置最小误差，最终确定分割阈值。

如图4所示，第一行为原始帧图像，图中圆域内像素值较小(颜色较深)为孔隙，像素值较大(颜色较浅)为骨架，图中第二行为Otsu分割的二值图像，没有考虑图像间上下层像素点的邻接关系，第三行为考虑了帧与帧之间的联系。通过对比可知，未考虑帧与帧之间像素关联二值图像有更多的孤立点和噪声，去除了孤立孔隙的二值图像黑色像素点明显增多，孔隙轮廓更加明显。由此可知，通过增加了帧与帧之间关联的空间标记方法有效去除图像的孤立孔隙点。

图4 二值图像对比图

5.3 多层切片连通域合并分析

为了验证图像合并的有效性和准确性，选取跟踪了中间位置帧图像0240.tiff图像的特定标签值观察连通域合并过程。

图5 连通域合并跟踪图

如图5所示，未进行合并时标签值连通区域面积较小，随着合并迭代次数的增加图像中标签值所代表的连通区域面积逐渐增加，表明单帧图像中不连通的区域，随着多帧图像标签的合并最终连接在一起，形成一个连通体。剔除孤立孔隙后图像中只剩与外界相同的连通体，如图5中最后一张图片所示。

5.4 数字岩心建模对比分析

为了验证通过连通域标记算法对岩心图像进行三维分割的有效性，分别选取了10、100、200、500张图像生成数字岩心模型，并比较效果，如图6所示。

图6 数字岩心孔隙模型对比图

文献[9]提出的半自动分割方法，需要手动对图像分割结果进行修复，由于不同的实验者具有不同的主观判断无法做到完全一致，造成实验结果具有不一致性。其他文献方法皆为对单帧图像的处理与文中方法有本质区别。文中提出的方法完全实现自动化，避免了人为因素的影响且考虑了岩心三维空间联系，因此不具有对比参考性，故不进行对照实验。

6 结束语

针对数字岩心三维建模分离骨架孔隙只考虑图像二维信息而忽略图像间的空间联系，导致生成的三维岩心模型不准确，孔隙连通情况不精确的问题，提出了基于二维连通域分析并结合多线程技术的空间标记并行分割方法。采用Otsu算法预分割岩心图像，通过对二值图像连通域合并，生成多个互不相通的连通体，根据边界限定准则去除孤立的孔隙，并根据理论孔隙度优化分割阈值。文中使用了650张吉林砂岩的CT图像，分别使用简单阈值分割方法和基于空间标记的串行分割方法作为对比。实验表明，基于空间标记的并行分割方法有效去除了图像中孤立的孔隙像素和噪声像素，提高了岩心模型的准确度和分割速度，为后续建立数字岩心三维模型奠定了基础。下一步将继续优化标签合并速度。