CN104376556B - 一种岩石ct图像目标分割方法 - Google Patents

Abstract

本发明公开了一种岩石CT图像分割方法，主要分为5个步骤：(1)对岩石CT图像进行自适应阈值分割，所谓自适应阈值分割是根据CT图像中相邻图像相关性，在分割过程中自动调节分割阈值；(2)对分割结果进行视觉判断，判断是否满意，是否需要进一步修复；(3)如果需要修复，则对分割结果自动判别分割优劣，进行自动区域生长修复；(4)对自动修复结果再此进行人眼判断，判断是否有个别感兴趣的目标需要进一步手工修复；(5)如果需要进一步手工修复，则进行手工选取特定目标的自动修复。通过本发明分割岩石CT图像的效率很高，整个过程中几乎没有人为参与，对大规模的CT图像分割效果很好。

Description

一种岩石CT图像目标分割方法

技术领域

本发明设计一种CT图像分割方法，尤其涉及一种岩石CT图像目标分割方法，属于CT图像分割技术领域。

技术背景

在石油地质分析过程中，CT成像技术(Computed Tomography，电子计算机X射线断层扫描技术)得到越来越广泛的应用，为了分析演示孔隙和颗粒结构，得到连通性、渗透率和形状因子等微观信息，通常采用对岩石进行CT扫描得到二维序列图，再根据二维序列图进行相应目标(孔隙和颗粒)分割提取，然后重建岩石的三维立体模型的方法。二维图像分割结果的好坏，直接关系到重建之后三维模型的真实程度。

图像分割是目标识别和特征提取的基础，是近年来计算机视觉研究的热点之一。图像分割技术发展到现在，产生了很多从经典到前沿的分割方法。主要分为以下几种：阈值分割，基于区域的分割，基于边缘的分割，基于特定数学理论的分割(包括数学形态学，模糊理论，神经网络等)。其中经典的分割技术当属阈值分割。根据特定的阈值来判断目标像素，从而得到感兴趣的目标。而根据阈值的选取的个数不同，阈值分割又分为单阈值，双阈值，多阈值等分割方法。

双阈值分割方法简单有效，速度较快。在岩石图像中，感兴趣的目标往往是孔隙和颗粒，其边缘比较模糊，需要通过预览手动调节合适的阈值来进行目标分割。但对于作为三维重建的岩石CT图像，图像的数量巨大，手动对每一张图像设定合适阈值进行分割的工作量非常大。针对这一问题，目前的一种解决方案是选择其中一幅图像，设定合适的阈值，其他所有图像都根据这个阈值进行分割，或者将CT图像分几个小组，每个小组选择一个合适的阈值进行分割。另一种解决方案是针对每一幅图像使用合适的分割方法，每一幅图像独立分割。第一种解决方案中分组分割是在分割效果和效率中间取一个平衡点。这种方法一方面平衡点的寻找需要多次的试探，工作量本身就很大，另一方面如果CT图像亮度变化比较剧烈，这种解决方案可能需要对每一幅图进行单独分割，工作量不会得到有效的降低。第二种解决方案理论上是最好的方法，整个过程几乎不需要进行人为干预。但是分割技术发展到现在，还没有一种分割方法对亮度不均匀的所有CT图像都能根据人的视觉判断进行理想的分割。况且如果分割方法过度复杂，对于大规模的CT图像进行处理的效率将会很低。而且在图像分割过程中，经常会有对分割结果进行局部手工修复的需求。并且期望手工修复的时候可以尽量降低工作量。因为当规模较大的时候，对每一幅图像进行手工修复的工作量是巨大的。然而对所有图像采用同一种固定修复方法，又很难达到预期的修复效果。针对这种需求，目前的研究情况还没有提出一种可以对岩石CT图像进行自适应，高效的分割和修复方法。本发明针对这种问题和需求提出并实现了一种可以根据人眼视觉判断，进行的自适应，并且高效的分割方法，并且利用序列图的相关性提出并实现了一种自适应的修复方法。可以在不加大工作量的前提下，进而对分割效果进一步提高。

发明内容

本发明的目的就是为了满足上述需求而提出并实现一种具有实用性的，高效的岩石CT图像目标分割方法。

本发明提供的岩石CT图像目标分割方法，是基于以下思想提出的：岩石CT图像具有规模大，灰度对比度低、纹理单一、噪声较多、且结构不规则等特征。由于CT扫描技术的原因，通常会还会造成CT图像前后亮度不均匀，亮度差别较大。岩石内部结构比较复杂，直接进行分割难度较大。但是CT图像有一个很重要的特征，就是相邻CT图的特征信息基本保持一致，包括亮度，结构，统计特性等。而且同一地质环境下的CT图像，孔隙和颗粒在岩石中的分布有很大相关性，特征变化缓慢。合理利用这些特性可以解决很多目前比较难以解决的问题。而且大规模图像处理过程中，图像处理的效果和效率，以及人所参与的工作量，是判断整个图像处理方法优劣的重要指标。均衡这些指标也正是图像处理方法研究的关键所在。图像分割作为图像处理中的代表应用技术不可能避开这些指标。在效果方面，优秀的图像分割技术要保证能将所有感兴趣的目标分割出来，并且非感兴趣目标要被抑制。由于CT图像虽然具有很大相关性，但是图像之间亮度有一定差别，采用固定阈值分割的时候，无法保证每一幅CT图像都达到预期的分割效果。而根据相关性信息，自动调节分割阈值，可以达到预期的分割效果。但即便阈值自动调节，分割效果得到很大改善，还会难免有少量图像达不到人眼视觉判断的分割效果，存在一定的欠分割现象。所谓欠分割就是没有分割彻底，较大的目标只被分割出少量区域。以当前的分割结果，自动进行区域生长，将目标区域进一步修复完善，便可以达到较理想的分割效果。在效率方面，工作量(包括人和计算机资源)要在可接受范围内尽可能降低。自动阈值调节和自适应区域生长就是在改善分割效果的同时兼顾效率的方法。即便在最终对个别目标进行人工修复的过程中，针对相关目标自适应进行区域生长修复也是将效率作为重要指标。只需对一幅图像中的目标进行手工修复，其他CT图中的目标将自动进行判别修复。综合要求分割方法简单有效，分割方法不能太复杂，需尽量使整个过程中没有人参与或者人参与成份很少。

本发明提供岩石CT图像目标分割方法，包括以下步骤：

步骤1：对岩石CT图像进行自动调节阈值分割，把感兴趣的目标提取出来；

步骤2：对步骤1目标提取的结果进行人眼视觉判断，分割效果是否达到要求，如果未达到要求，需要进一步对分割结果修复，进行步骤3，否则分割结束；

步骤3：对步骤1分割结果进行自动区域生长修复，使提取的目标更加完善；

步骤4：对步骤3目标修复的结果进行人眼视觉判断，目标修复效果是否达到要求，如果没有达到要求，需要进一步对分割修复的结果进行手工修复，进行步骤5，否则分割结束；

步骤5：对步骤3分割修复结果进行手工选取特定目标自动修复，以修复的结果作为最终提取的目标。

在本发明上述技术方案中，所述步骤1对岩石CT图像进行自动调节阈值分割，优先采用由下述方法进行分割：

(1)选定CT图像中的第一幅图像，通过双阈值二值化方法，手动调阈值，根据人眼视觉对感兴趣目标区域的判断，选定合适的阈值(th_L，th_U),其中灰度值在[th_l,th_U]范围内的像素点属于分割区域，使当前图像的分割效果较为理想。

(2)计算图像的直方图H₀(x)，根据H₀(x)记录用于分割的下阈值点(th_L，H₀(th_L))和与其相邻的左右各5个点组成下阈值匹配模版T_l。同样记录用于分割的上阈值点(th_U，H₀(th_U)和与其相邻的左右各5个点组成上阈值匹配模版T_u，T_l和T_u的中心点分别为(th_L，H₀(th_L))和(th_U，H₀(th_U))。

(3)计算当前图像的直方图H(x),使用T_l和T_u在H(x)上进行匹配。计算T_l在H(x)上对应位置差值绝对值τ，找到使τ最小的中心点thl_i，将thl_i用于当前图像分割的下阈值点。使用同样的方法找到T_u最匹配的中心点thu_i，将thu_i用于当前图像分割的上阈值点。

(4)对当前图像使用(thl_i，thu_i)阈值进行分割，并用thl_i更新th_L，用thu_i更新th_U，用H(x)更新H₀(x)。

(5)重复步骤(2)，(3)，(4)，直到所有图像分割完毕。

在本发明上述技术方案中，所述步骤3对分割后的岩石CT图像进行自动区域生长修复，优先采用由下述方法进行修复：

(1)根据人眼视觉判断，选择一幅分割比较理想的图片(通常直接取第一幅)。计算分割区域的覆盖率δ₀，以δ₀作为分割优劣的标准。

(2)计算CT图中每一幅图已分割的覆盖率δ，以及和当前图像覆盖率和δ₀的比值ω，根据ω将当前图像的分割程度归类到A₁,A₂,A₃,A₄。

(3)对CT图像进行扫描，通过区域生长方法对已分割区域进行识别，得到已分割目标区域的像素点集Ω，并计算目标区域Ω的平均灰度值Ψ。

(4)计算判别像素点的8领域灰度均值和Ψ的差值△以及8领域灰度值与当前像素点灰度的平均差值ε。将△和ε作为区域生长的判别准则,并根据A₁,A₂,A₃,A₄不同类别自动设定特定的△和ε判别值，对目标区域进行区域生长。

(5)在Ω的基础上根据步骤(4)中的判别规则对目标区域进行改进的区域生长方法修复，使目标区域分割更加完整。

(6)重复步骤(4)和(5)直到图像中所有目标都修复完毕。

在本发明上述技术方案中，所述步骤5对分割修复结果进行手工选取特定目标自动修复，优先采用由下述方法进行自动修复：

(1)对CT图中个别分割不理想的目标进行手工修复，得到修复区域的点集Φ。

(2)将Φ所有像素点作为种子点，计算Φ区域的加权平均灰度值φ。判别领域像素点P，计算P的8领域灰度均值与φ的差值△。并计算P的8领域灰度值λ_i与当前像素点的灰度值λ₄的平均差值ε。这两个值在一定范围内的像素点将归并到当前区域中。最终的到区域A。

(3)计算区域A的面积s₀和像素加权平均值g₀，并运用拉普拉斯变换在图形层上提取A的边缘D。

(4)将区域像素点集A投影到下一幅相邻图像，通过对应位置图像灰度值的比较，进行像素点取舍，得到投影点集B。

(5)以点集B代替Φ作为种子点，重复步骤(3)。

(6)将当前区域生长完成的区域C与A做比较，将C中存在而在A中不存在的像素点集标记为区域F。

(7)计算边缘D的总线长d₀，以及区域F与边缘D的交线线长d₁，并计算d₁和d₀的比值μ_d，根据μ₀判断F是从中C删除还是保留。

(8)计算区域C的面积s和其在原始图像的加权平均灰度值g。计算s与s₀进行比值μ_s，以及g和g₀的差值绝对值μ_g。用μ_s和μ_g的取值作为方法的终止条件。

(9)用C更新A，s更新s₀的值，保持g₀不变，重复步骤(3)，(4)，(5)，(6)，(7)，(8)。

在所述对岩石CT图像进行自动调节阈值分割的方法步骤(1)中，通过人眼的判断来决定分割阈值范围[th_L，th_U]。保证分割效果较佳。

在所述对岩石CT图像进行自动调节阈值分割的方法步骤(2)中，计算图像的直方图H₀(x)，计算公式为：

其中N(x)为灰度级x在图像中的像素点数，x∈[0,255]。N为图像总共像素点数。根据H₀(x)记录用于分割的下阈值点(th_L，H₀(th_L))和与其相邻的左右各5个点组成下阈值匹配模版T_l。同样记录用于分割的上阈值点(th_U，H₀(th_U)和与其相邻的左右各5个点组成上阈值匹配模版T_u，T_l和T_u的中心点分别为(th_L，H₀(th_L))和(th_U，H₀(th_U))。

在所述对岩石CT图像进行自动调节阈值分割的方法步骤(3)中，使用公式(1)计算当前图像的直方图H(x),使用T_l和T_u在H(x)上进行匹配，匹配方式是在H(x)上定位到[th_L-5，th_L+5]区间，使T_l中心点在[th_L-5，th_L+5]内从中间到两边交替移动，每移动一个位置都会有一个与T_l相对应的模版，其中心点为thl_i，范围区间为[thl_i-5，thl_i+5]。计算T_l在H(x)上对应位置差值绝对值和τ，计算公式为：

找到使τ最小的中心点thl_i，将thl_i用于当前图像分割的下阈值点。使用同样的方法找到T_u最匹配的中心点thu_i，将thu_i用于当前图像分割的上阈值点。

在所述对岩石CT图像进行自动调节阈值分割的方法步骤(4)中，对当前图像使用(thl_i，thu_i)阈值进行分割，并用thl_i更新th_L，用thu_i更新th_U，用H(x)更新H₀(x)。

在所述对岩石CT图像进行自动调节阈值分割的方法步骤(5)中，步骤(2)，(3)，(4)，可以完成根据上一幅图像来完成当前图像的分割，重复步骤(2)，(3)，(4)，可以完成所有CT图的分割。

在所述对分割后的岩石CT图像进行自动区域生长修复的方法步骤(1)中，根据人视觉判断，选择一幅分割较理想的图片(通常直接取第一幅),设图像长和宽分别为w和h,已分割出来的区域像素点数为l。计算分割区域的覆盖率δ₀，计算公式如下：

以δ₀作为分割优劣的标准。如果δ₀比较小表明分割区域较少，如果δ₀比较大代表分割区域较大。一定程度上可以表征当前图像的分割情况。

在所述对分割后的岩石CT图像进行自动区域生长修复的方法步骤(2)中，计算CT图中每一幅图已分割的覆盖率δ和当前图像覆盖率和δ₀的比值计ω，算公式如下：

根据ω将当前图像的分割程度归类到A₁,A₂,A₃,A₄。具体分类情况如表1-1：

表1-1根据ω将图像归类情况

ω	类别	描述
			(0.0，0.3]		过度欠分割
(0.3，0.5]		欠分割
			(0.5，0.8]		轻微欠分割
(0.8，1.0]		分割效果较好

原理是δ₀是人为判断比较满意的分割区域的覆盖率，在同一地质层中岩石的孔隙覆盖率不会相差很大。用当前图像分割区域的覆盖率δ和δ₀进行比较，便可以得出当前图像的分割情况。

在所述对分割后的岩石CT图像进行自动区域生长修复的方法步骤(3)中，对已分割区域进行识别，其实这一步是在图形层上进行的。图像分割的过程中，把分割出来的区域单独生成一个图形层，保持原始图像数据不变，显示过程中使图形层和原始图像叠加显示。图形层是一个和原始图像素点位置对应的等大二值图像。通过对图形层扫描和区域生长，可以把当前图像分割出来的每一个目标识别出来。一旦识别到一个目标区域Ω，就进行步骤(4)和(5)。Ψ是原始图像在Ω相应的像素点的加权平均值。计算公式如下：

其中N为区域Ω所有的像素点数，H_i为i处像素点的灰度值，P(i)为H_i的概率。

在所述对分割后的岩石CT图像进行自动区域生长修复的方法步骤(4)中，像素8领域示意图如图7-1，计算判别像素点P的8领域灰度均值和Ψ的差值△以及8领域灰度值与P灰度的平均差值ε，计算公式如下：

其中λ₄为P的灰度值，λ_i为P的8领域像素灰度值。设初始设定领域均值与加权均值差值为△₀和领域差值平均为ε₀。根据表1-2自动设定判别准则，根据不同的类型设定特定的判别阈值可以有效满足各类图像都达到预期的分割效果。如果都使用统一判别阈值，将无法避免有些图像过生长和有些图像欠生长的情况。而且由于自动调节阈值分割所得到的大部分分割效果都较好，因此分割后的CT图像中大部分图像都将属于A₄类。A₄类在这里是不需要区域生长修复的，这大大提升了区域生长分割修复的效率。以△和ε作为区域生长的判别准则。前者有效的抑制了噪声，后者可以一定程度上检测到边缘，防止过生长现象。

表1-2区域生长自动判别准则设置标准

类别	领域均值与加权均值差值	领域差值平均
	无	无

在所述对分割后的岩石CT图像进行自动区域生长修复的方法步骤(5)中，为进一步改善已分割区域Ω的分割效果，根据步骤(4)中的判别规则对目标区域Ω进行区域生长方法修复。

在所述对分割后的岩石CT图像进行自动区域生长修复的方法步骤(6)中，步骤(4)和(5)所处理的是一幅图像中的一个目标，重复上述两个步骤，就可以把每一幅图像中每一个欠分割的目标区域进行修复。

在所述对分割修复结果进行手工选取特定目标自动修复的方法步骤(1)中，在自动化分割完全结束之后，如果需要还可以根据需求进一步手工修复。但即便是手工修复，也要尽可能减少工作量。于是实现了在一幅图像中手工修复一个目标区域，其他图像中的所有相应的目标区域都会自动得到修复的方法。通过手动画取或者手动设定种子点自动区域生长得到区域Φ。

在所述对分割修复结果进行手工选取特定目标自动修复的方法步骤(2)中，将Φ所有像素点作为种子点，计算Φ区域的加权平均灰度值φ，计算公式如下：

其中N为区域Φ所有的像素点数，H_j为j处像素点的灰度值，P(j)为H_j的概率值。判别领域像素点P，计算P的8领域灰度均值与φ的差值△，计算公式如下：

其中λ_i为P的8领域的灰度值，λ₄为当前像素点的灰度值。并计算P的8领域灰度值λ_i与当前像素点的灰度值λ₄的平均差值ε，计算公式如下：

以△≤m(m∈[0,50])和ε≤n(n∈[0,50])(根据需要可自行调节)作为区域生长的判别准则，得到区域生长完成的区域A。

在所述对分割修复结果进行手工选取特定目标自动修复的方法步骤(3)中，计算区域A的面积s₀(即区域A像素点数)和A在原始图像上的灰度加权平均值g₀，计算公式如下：

其中N为区域A所有的像素点数，H_j为j处像素点的灰度值，P(j)为H_j的概率值。在A上运用拉布拉斯变换，提取A的边缘D。变换使用的模版为图7-2。用此模版遍历区域A中的每一个像素点，如果模版结果为0将此像素点删除，否则保留，便提取出区域A的边缘像素点集D。

在所述对分割修复结果进行手工选取特定目标自动修复的方法步骤(4)中，将区域像素点集A投影到下一幅相邻图像，在原始图像上，比较两幅图像在区域A位置上的像素灰度值(设分别为P_i,P_j)，使用公式|P_i-P_j|≤k(k∈[0,50])筛选出A在下一幅相邻图像上的投影点集B。

在所述对分割修复结果进行手工选取特定目标自动修复的方法步骤(5)中，以点集B代替Φ为种子点，重复步骤5.2，得到当前图像区域生长的区域C。

在所述对分割修复结果进行手工选取特定目标自动修复的方法步骤(6)中，将区域C和区域A做比较，将在C中存在而在A中不存在的像素点集标记为区域F。表示为下：

F＝C∩C_SA,S为整幅图像区域 (12)

在所述对分割修复结果进行手工选取特定目标自动修复的方法步骤(7)中，计算边缘D的总线长d₀(即D像素点总数)以及区域F与边缘D的交线线长d₁。即d₁为F∩D的像素点集中像素点个数。计算d₁和边缘D的总线长d₀的比值μ_d，计算公式如下：

根据μ_d判断F是从中C删除还是保留。如果μ_d≤f(f∈[0,1]),表示区域F不属于感兴趣的区域的，将F从C中删除。否则，保留F。这样通过相邻图像的矫正，有效的避免了过分割现象的发生。

在所述对分割修复结果进行手工选取特定目标自动修复的方法步骤(8)中，计算区域C的面积s和其在原始图像的加权平均灰度值g。用s和s₀进行比较，g和g₀比较。μ_g＝|g-g₀|。若μ_s≤t₁(t₁∈[0,1])或者μ_g≥t₂(t₂∈[0,50]),表明目标区域消失，方法终止。

在所述对分割修复结果进行手工选取特定目标自动修复的方法步骤(9)中，用C更新A，s更新s₀的值，保持g₀不变，重复步骤(3)，(4)，(5)，(6)，(7)，(8)将CT图中同一个目标分割出来。

本发明提供的岩石CT图像目标分割方法，较之现有技术具有以下十分突出的技术效果：

本发明以岩石CT图像的分割效果为基础，兼顾工作量和分割效率，能够快速有效的批量大规模分割出人眼感兴趣的目标区域。而且整个分割过程中不仅能够全部自动化完成，而且在需要的时候可以加入极少量的人工调节，让分割效果得到进一步改善。自动化过程是对CT图中所有相关目标进行自适应分割，但本发明也可以根据人感兴趣的一个目标或几个目标进行选择性分割，并且分割效率和效果都能达到要求，极大的增强了分割的灵活性和方便性。

附图说明

图1-1是本发明实施例中岩石CT原始图像之一(CT图第一张)；

图1-2是本发明实施例中岩石CT原始图像之二(与图1-1相隔200张)；

图1-3是本发明实施例中岩石CT原始图像之三(与图1-2相隔200张)；

图1-4是本发明实施例中岩石CT原始图像之四(与图1-3相隔100张)；

图1-5是本发明实施例中岩石CT原始图像之五(与图1-4相隔3张)；

图1-6是本发明实施例中岩石CT原始图像之六(需手工修复的第一张)；

图1-7是本发明实施例中岩石CT原始图像之七(与图1-6相隔3张)；

图1-8是本发明实施例中岩石CT原始图像之八(与图1-7相隔3张)；

图2-1是本发明实施例中图1-1经过预处理后图像；

图2-2是本发明实施例中图1-2经过预处理后图像；

图2-3是本发明实施例中图1-3经过预处理后图像；

图2-4是本发明实施例中图1-4经过预处理后图像；

图2-5是本发明实施例中图1-5经过预处理后图像；

图2-6是本发明实施例中图1-6经过预处理后图像；

图2-7是本发明实施例中图1-7经过预处理后图像；

图2-8是本发明实施例中图1-8经过预处理后图像；

图3-1是本发明实施例中图2-1使用固定双阈值分割结果；

图3-2是本发明实施例中图2-2使用固定双阈值分割结果；

图3-3是本发明实施例中图2-3使用固定双阈值分割结果；

图3-4是本发明实施例中图2-4使用固定双阈值分割结果；

图3-5是本发明实施例中图2-5使用固定双阈值分割结果；

图4-1是本发明实施例中图2-1使用自动调节双阈值分割结果；

图4-2是本发明实施例中图2-2使用自动调节双阈值分割结果；

图4-3是本发明实施例中图2-3使用自动调节双阈值分割结果；

图4-4是本发明实施例中图2-4使用自动调节双阈值分割结果；

图4-5是本发明实施例中图2-5使用自动调节双阈值分割结果；

图4-6是本发明实施例中图2-6使用自动调节双阈值分割结果；

图4-7是本发明实施例中图2-7使用自动调节双阈值分割结果；

图4-8是本发明实施例中图2-8使用自动调节双阈值分割结果；

图5-1是本发明实施例中图4-4自动区域生长修复后的结果；

图5-2是本发明实施例中图4-5自动区域生长修复后的结果；

图6-1是本发明实施例中图4-6手动修复后的结果；

图6-2是本发明实施例中图4-7手动修复后的结果；

图6-3是本发明实施例中图4-8手动修复后的结果；

图7-1是本发明中所使用的像素点的8领域示意图；

图7-2是本发明中所使用的拉普拉斯变换模版示意图；

图8是本发明的方法流程图。

具体实施方式

下面结合具体实施例和附图对本发明作进一步说明：

实施例：

为了使本发明所述分割方法更加便于理解和接近于真实应用，下面从原始岩石CT图像开始到图像分割完成为止进行整个流程的整体说明，其中包括本发明的核心分割方法,具体步骤如下：

(1)准备一组待分割的岩石CT序列图，图片数量为1000，格式为BMP，尺寸988X1013。

(2)使用图像批处理软件打开岩石CT图像。其中代表如图1-1，图1-2，图1-3，图1-4，图1-5，图1-6，图1-7，图1-8。

(3)对CT图像选定处理区，并对处理区内进行自动白平衡、自动光源矫正、中值滤波等图像预处理，增加图像的对比度，并去除图像中噪声。结果如图2-1，图2-2，图2-3，图2-4，图2-5，图2-6，图2-7，图2-8。

(4)选定步骤(3)中经预处理后的CT图中的一张图像(默认为第一张)即图1-1，使用双阈值分割，手动调节上下阈值，通过预览查看，对感兴趣的目标区域得到比较满意的分割效果，并记录上下阈值(th_L，th_U)。

(5)计算图像的直方图H₀(x)，根据H₀(x)记录用于分割的下阈值点(th_L，H₀(th_L))和与其相邻的左右各5个点组成下阈值匹配模版T_l。同样记录用于分割的上阈值点(th_U，H₀(th_U)和与其相邻的左右各5个点组成上阈值匹配模版T_u，T_l和T_u的中心点分别为(th_L，H₀(th_L))和(th_U，H₀(th_U))。

(6)计算当前图像的直方图H(x)，使用上幅图像分割后的直方图模版T_l和T_u在当前图像直方图上进行匹配,找到T_l最匹配的中心点将此中心点thl_i。使用同样的方法找到T_u最匹配的中心点thu_i。

(7)对当前图像使用(thl_i，thu_i)进行分割，用thl_i更新th_L，用thu_i更新th_U。(8)重复步骤(5)，(6)，(7)，直到所有图像分割完毕。结果如图4-1，图4-2，图4-3，图4-4，图4-5，图4-6，图4-7，图4-8。

(9)根据人眼视觉判断，选择一幅分割较理想的图片(默认第一幅)即图4-1,计算分割区域的覆盖率δ₀作为分割优劣的标准。

(10)计算CT图中每一幅图已分割的覆盖率δ，以及和当前图像覆盖率和δ₀的比值ω，通过ω的值将当前图像的分割程度归类为A₁,A₂,A₃,A₄。

(11)对CT图像进行扫描，通过区域生长方法对已分割区域进行识别，得到已分割目标区域的像素点集Ω，并计算目标区域Ω的加权平均灰度值Ψ。

(12)设初始设定领域均值与加权均值差值为△₀和领域差值平均为ε₀。并根据表1-2自动根据类别设定判别准则。

(13)在Ω的基础上根据步骤(12)中的判别规则对目标区域进行改进的区域生长方法修复，使目标区域分割更加完整。

(14)重复步骤(11)和(12)直到图像中所有目标都修复完毕。代表结果如图5-1，图5-2可与图4-4和图4-5进行对比，从而观察修复效果。

(15)根据需求对CT图(图4-6)中个别需要手工修复的目标进行手工修复，得到修复区域的点集Φ。

(16)将Φ所有像素点作为种子点，用公式(5)、(6)、(7)分别计算φ、△、ε。以△≤10和ε≤6(根据需要可自行调节)作为区域生长的判别准则，得到区域生长完成的区域A。

(17)计算区域A的面积s₀(即像素点数)和A在原始图像上的灰度加权平均值g₀。对A运用拉普拉斯变换，提取A的边缘D。

(18)将区域像素点集A投影到下一幅相邻图像，通过对应位置图像灰度值的比较，进行像素点取舍，得到投影点集B。

(19)以点集B代替Φ为种子点，重复步骤(16)，得到当前图像区域生长的区域C。

(20)将C与A做比较，将C中存在而在A中不存在的像素点集标记为区域F。

(21)计算边缘D的总线长d₀(即D像素点总数)以及区域F与边缘D的交线线长d₁。即d₁为F∩D的像素点集中像素点个数。

(22)计算d₁和d₀的比值μ_d，根据μ_d判断F是从C中删除还是保留。

(23)计算区域C的面积s和其在原始图像的加权平均灰度值g。计算s与s₀进行比值μ_s，以及g和g₀的差值绝对值μ_g。用μ_s≤0.5或者μ_g≥10作为方法的终止条件。

(24)用C更新A，s更新s₀的值，保持g₀不变，重复步骤(17)，(18)，(19)，(20)，(21)，(22)直到修复完毕，结果如图6-1，图6-2，图6-3，可与图4-6，图4-7，图4-8进行对比，观察效果。

本实施例中，以自动调节阈值方法分割岩石CT序列图像，并根据分割效果自动局部调节和修复。为达到更好的分割效果，也可以适量加入少量的人工干预。保证分割效果的同时，提高分割效率，减少工作量。并其一方面将分割结果与原始岩石CT序列图像进行直观地对比观察，另一方面与固定双阈值分割方法结果进行直观地对比观察，从而验证了本发明的可靠性与实用性。

上述实施例只是本发明的较佳实施例，并不是对本发明技术方案的限制，只要是不经过创造性劳动即可在上述实施例的基础上实现的技术方案，均应视为落入本发明专利的权利保护范围内。

Claims (9)

1.一种岩石CT图像目标分割方法，包括以下步骤：

步骤1：对岩石CT图像进行自动调节阈值分割，把感兴趣的目标提取出来；

步骤2：对步骤1目标提取的结果通过人眼视觉进行判断，判断分割效果是否达到要求，如果未达到要求，需要进一步对分割结果修复，进行步骤3，否则分割结束；

步骤3：对步骤1分割结果进行自动区域生长修复，使提取的目标更加完善；

步骤4：对步骤3目标修复的结果通过人眼视觉进行判断，判断目标修复效果是否达到要求，如果没有达到要求，需要进一步对分割修复的结果进行手工修复，进行步骤5，否则分割结束；

步骤5：对步骤3分割修复结果进行手工选取特定目标自动修复，以修复的结果作为最终提取的目标；

步骤1所述的对岩石CT图像进行自动调节阈值分割，采取下述方法进行：

(1)选定CT图像中的一幅图像，作为分割开始的基准图像，通过双阈值二值化分割方法，根据人眼视觉对感兴趣目标区域的判断，手动选定阈值区间[th_L，th_U],其中灰度值在[th_L，th_U]范围内的像素点属于分割区域；

(2)计算基准图像的直方图H₀(x)，根据H₀(x)记录用于分割的下阈值点(th_L，H₀(th_L))和与其相邻的左右各5个点组成下阈值匹配模版T_l，同样记录用于分割的上阈值点(th_U，H₀(th_U)和与其相邻的左右各5个点组成上阈值匹配模版T_u，T_l和T_u的中心点分别为(th_L，H₀(th_L))和(th_U，H₀(th_U))；

(3)计算当前图像的直方图H(x),使用T_l和T_u在H(x)上进行匹配，计算T_l在H(x)上对应位置差值绝对值和τ，找到使τ最小的中心点thl_i，将thl_i用于当前图像分割的下阈值点，使用同样的方法找到T_u最匹配的中心点thu_i，将thu_i用于当前图像分割的上阈值点；

(4)对当前图像使用(thl_i，thu_i)阈值进行分割，并用thl_i更新th_L，用thu_i更新th_U，用H(x)更新H₀(x)；

(5)重复步骤(2)，(3)，(4)，直到所有图像分割完毕。

2.根据权利要求1所述的岩石CT图像目标分割方法，其特征在于，步骤3对自动调节阈值分割结果进行自动区域生长修复，采取下述方法进行：

(1)根据人眼视觉判断，选择一幅分割较理想的图像作为基准图像，计算其分割区域的覆盖率δ₀，以δ₀作为分割优劣的标准；

(2)计算CT图像中每一幅图已分割的覆盖率δ，以及δ和δ₀的比值ω，根据ω将当前图像的分割程度归类到A₁,A₂,A₃,A₄；

(3)对CT图像进行扫描，通过区域生长方法对已分割区域进行识别，得到已分割目标区域的像素点集Ω，并计算Ω的平均灰度值Ψ；

(4)计算判别像素点P的8邻域灰度均值和Ψ的差值Δ以及8邻域灰度值与判别像素点P灰度的平均差值ε，将Δ和ε作为区域生长的判别准则,并根据A₁,A₂,A₃,A₄不同类别自动设定特定的Δ和ε判别值，对目标区域进行区域生长；

(5)在Ω的基础上根据步骤(4)中的判别规则对目标区域进行改进的区域生长方法修复，使目标区域分割更加完整；

(6)重复步骤(4)和(5)直到图像中所有目标都修复完毕。

3.根据权利要求1所述的岩石CT图像目标分割方法，其特征在于，步骤5对分割结果进行手工选取特定目标自动修复，采取下述方法进行：

(1)对CT图像中分割不理想的目标进行手工修复，得到修复区域的点集Φ；

(2)将Φ所有像素点作为种子点，计算Φ区域的加权平均灰度值φ和计算判别像素点P的8邻域灰度均值与φ的差值Δ，并计算P的8邻域灰度值λ_i与判别像素点P的灰度值λ₄的平均差值ε，这两个值在一定范围内的像素点将归并到当前区域中，最终得到区域A；

(3)计算区域A的面积s₀和像素加权平均值g₀，并运用拉普拉斯变换在图形层上提取A的边缘D；

(4)将区域A投影到下一幅相邻图像，通过对应位置图像灰度值的比较，进行像素点取舍，得到投影点集B；以点集B代替Φ作为种子点，重复步骤(2)得到生长完成的区域C；

(5)将当前区域生长完成的区域C与区域A做比较，将C中存在而在A中不存在的像素点集标记为区域F；

(6)计算边缘D的总线长d₀，以及区域F与边缘D的交线线长d₁，计算d₁和d₀的比值μ_d，根据μ_d判断F是从C中删除还是保留；

(7)计算区域C的面积s和其加权平均灰度值g，计算s与s₀的比值μ_s，以及g和g₀的差值绝对值μ_g，用μ_s和μ_g的值作为方法的终止条件；

(8)用C更新A，s更新s₀的值，保持g₀不变，重复步骤(3)，(4)，(5)，(6)，(7)直到修复完毕。

4.根据权利要求1所述的岩石CT图像目标分割方法，其特征在于，在对岩石CT图像进行自动调节阈值分割的步骤(3)中，计算当前图像的直方图H(x),使用T_l和T_u在H(x)上进行匹配，匹配方式为在H(x)上定位到[th_L-5，th_L+5]区间，使T_l中心点在[th_L-5，th_L+5]内从中间到两边交替移动，每移动一个位置都会有一个与T_l相对应的模版，其中心点为thl_i，范围区间为[thl_i-5，thl_i+5]。

5.根据权利要求1所述的岩石CT图像目标分割方法，其特征在于，在对岩石CT图像进行自动调节阈值分割的步骤(3)中计算T_l在H(x)上对应位置差值绝对值和τ的计算公式如下：

$\mathrm{\τ}=\underset{th={\mathrm{thl}}_i+j}{\overset{5}{\underset{{\mathrm{th}}_0={\mathrm{th}}_L+j}{\underset{j=-5}{\Σ}}}}|H\left(th\right)-H_0\left({\mathrm{th}}_0\right)|---\left(1\right).$

6.根据权利要求2所述的岩石CT图像目标分割方法，其特征在于，在对自动调节阈值分割结果进行自动区域生长修复的步骤(2)中，计算CT图像中每一幅图已分割的覆盖率δ以及δ和δ₀的比值ω，根据ω将当前图像的分割程度归类到A₁,A₂,A₃,A₄；具体分类情况为当ω∈(0.0，0.3]时，图像过度欠分割，属于A₁类；当ω∈(0.3，0.5]时，图像欠分割，属于A₂类；当ω∈(0.5，0.8]时，图像轻微欠分割，属于A₃类；当ω∈(0.8，1.0]时，图像分割效果较好，属于A₄类。

7.根据权利要求2所述的岩石CT图像目标分割方法，其特征在于，在对自动调节阈值分割结果进行自动区域生长修复的步骤(4)中，计算判别像素点P的8邻域灰度均值和Ψ的差值Δ以及8邻域灰度值与判别像素点P灰度的平均差值ε，Δ和ε的计算公式如下：

$\mathrm{\Δ}=|\underset{i=0}{\overset{8}{\Σ}}\frac{{\mathrm{\λ}}_i}{9}-\mathrm{\Ψ}|,i\∈\[0,8\]---\left(2\right)$

$\mathrm{\ϵ}=\frac{\underset{i=0}{\overset{8}{\Σ}}|{\mathrm{\λ}}_i-{\mathrm{\λ}}_4|}{8},i\∈\[0,8\]---\left(3\right)$

其中λ₄为P的灰度值，λ_i为P的8邻域像素灰度值，设定初始邻域灰度均值和Ψ的差值为Δ₀，初始邻域灰度值与P灰度的平均差值为差值ε₀，根据A₁,A₂,A₃,A₄自动设定特定的Δ和ε判别值。

8.根据权利要求3所述的岩石CT图像目标分割方法，其特征在于，在对分割结果进行手工选取特定目标自动修复的步骤(2)中，计算Φ区域的加权平均灰度值φ，计算判别像素点P的8邻域灰度均值与φ的差值Δ，其中Δ的计算公式如下所示：

$\mathrm{\Δ}=|\underset{i=0}{\overset{8}{\Σ}}\frac{{\mathrm{\λ}}_i}{9}-\mathrm{\φ}|,i\∈\[0,8\]---\left(4\right)$

其中λ_i为P的8邻域的灰度值，λ₄为判别像素点P的灰度值，并计算P的8邻域灰度值λ_i与判别像素点P的灰度值λ₄的平均差值ε，其中ε的计算公式如下所示：

$\mathrm{\ϵ}=\frac{\underset{i=0}{\overset{8}{\Σ}}|{\mathrm{\λ}}_i-{\mathrm{\λ}}_4|}{8},i\∈\[0,8\]---\left(5\right)$

以Δ≤m和ε≤n作为区域生长的判别准则,其中m∈[0，50],n∈[0，50]。

9.根据权利要求3所述的岩石CT图像目标分割方法，其特征在于，在对分割结果进行手工选取特定目标自动修复的步骤(6)中，边缘D的总线长d₀，区域F与边缘D的交线线长d₁为F∩D的像素点集中像素点个数，计算d₁和d₀的比值μ_d，用μ_d判断F是从C中删除还是保留，如果μ_d≤f，将F从C中删除，否则保留F，其中f∈[0，1]。