陳國(guó)軍，李 勝

(中國(guó)石油大學(xué)(華東) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，山東 青島 266580)

0 引 言

圖像分割是指將圖像分成若干具有相似性質(zhì)的區(qū)域的過(guò)程，是許多圖像處理任務(wù)的預(yù)處理步驟?，F(xiàn)有的圖像分割方法主要分以下幾類(lèi)：基于閾值的分割方法、基于區(qū)域的分割方法[1]、基于邊緣的分割方法[2]以及基于特定理論[3]的分割方法等。

巖心圖像中感興趣的目標(biāo)是孔隙和骨架，三維巖心模型建立的真實(shí)程度和三維分析結(jié)果依賴于巖心二維圖像分割結(jié)果的好壞。目前針對(duì)巖心圖像分割方法主要有兩種，一種是選擇其中的一幀圖像根據(jù)實(shí)際孔隙度計(jì)算合適的閾值，并以該閾值對(duì)所有圖像進(jìn)行閾值分割，再根據(jù)實(shí)際孔隙度進(jìn)行閾值優(yōu)化；另一種是根據(jù)實(shí)際孔隙度對(duì)每一幀圖像設(shè)定獨(dú)立的閾值進(jìn)行分割。文獻(xiàn)[4-6]中提出了基于最大類(lèi)間方差法(Otsu)[7]的不同改進(jìn)方法，并取得了不錯(cuò)的效果；文獻(xiàn)[8]提出了適用于體視圖像巖粒分割的算法，采用兩次分水嶺并結(jié)合最小割的方法對(duì)巖石圖像進(jìn)行分割；徐永進(jìn)等[9]提出了利用圖像層間相關(guān)性實(shí)現(xiàn)CT圖像半自動(dòng)分割，原理是利用分割效果好的幀圖像修復(fù)其余幀圖像；文獻(xiàn)[10]提出了一種基于模糊距離變換的改進(jìn)分水嶺[11-12]算法，通過(guò)對(duì)區(qū)域生長(zhǎng)速度進(jìn)行控制，優(yōu)化了過(guò)分割問(wèn)題；文獻(xiàn)[13]利用了薄片視域的單偏光圖與正交偏光序列圖信息對(duì)演示薄片骨架孔隙進(jìn)行分割；文獻(xiàn)[14]采用聚類(lèi)的方法對(duì)巖石CT圖像分割及量化。不論是將一幀圖像的分割閾值用于所有圖像分割，還是每一幀圖像采用不同的閾值進(jìn)行分割，雖然也能取得不錯(cuò)的效果，但都只是考慮了巖心截面的二維信息，而忽略了巖心幀圖像在三維空間的聯(lián)系。為解決此問(wèn)題，該文提出了基于空間標(biāo)記的巖心孔隙并行分割方法，通過(guò)采用連通域標(biāo)記的方法建立巖心圖像在三維空間中的連通關(guān)系并實(shí)現(xiàn)自動(dòng)分割。

1 圖像預(yù)處理

Micro-CT掃描得到的巖心斷層圖像亮度較低，對(duì)于尺寸小于分辨率或尺寸在分辨率附近的特征體很難分辨，尤其是孔隙和骨架的邊界區(qū)域，通常較為模糊，沒(méi)有明顯的分界，多為過(guò)渡灰度值，需對(duì)巖心圖像進(jìn)行對(duì)比度增強(qiáng)和濾波處理。直方圖均衡化實(shí)質(zhì)是圖像增強(qiáng)的一種方式，經(jīng)過(guò)直方圖均衡化的圖像對(duì)比度明顯增強(qiáng)，使得圖像有很強(qiáng)的清晰感。直方圖均衡化處理的原理是把原始圖像的灰度直方圖從比較集中的某個(gè)灰度區(qū)間變成在全部灰度范圍內(nèi)的均勻分布，對(duì)圖像進(jìn)行非線性拉伸，重新分配圖像像素值，使一定灰度范圍內(nèi)的像素?cái)?shù)量大致相同，就是把給定圖像的直方圖分布改變成“均勻”分布直方圖分布。中值濾波是非線性濾波，避免因線性濾波方法導(dǎo)致的邊界模糊問(wèn)題，并且不受奇異值的影響，可以保護(hù)圖像細(xì)節(jié)，故在巖心濾波處理中被廣泛采用。經(jīng)過(guò)直方圖均衡化和中值濾波后巖心圖像骨架孔隙邊界模糊問(wèn)題得到改善，有利于進(jìn)一步對(duì)圖像進(jìn)行處理。

圖像二值化是巖心切片處理的中間步驟，通過(guò)設(shè)置閾值T將圖像分為像素值為255或0的兩類(lèi)。巖心圖像是在相同的環(huán)境中由CT機(jī)掃描而來(lái)，因此同組切片采用同一閾值分割。首先選取具有代表性的起始、中間和末尾三張切片采用Otsu計(jì)算初始分割閾值，并根據(jù)理論孔隙度優(yōu)化分割閾值，最后以該閾值二值化所有切片圖像得到包含孤立孔隙的二值圖像。

圖像二值化公式為：

(1)

式中，f(x,y)為當(dāng)前坐標(biāo)點(diǎn)的像素值，g(x,y)為二值化后該點(diǎn)的像素值，T為分割閾值，其中0≤T≤255。

2 二值圖像連通域標(biāo)記

二值圖像連通域標(biāo)記依據(jù)掃描圖像次數(shù)分為一趟掃描算法[15]和兩趟掃描算法[16]，一趟掃描算法主要有種子填充算法和游程標(biāo)記算法；兩趟掃描算法典型代表為基于快速連通域標(biāo)記的Two-Pass算法。該文采用基于并查集的兩遍掃描算法標(biāo)記單幀圖像連通域，通過(guò)結(jié)合多線程技術(shù)對(duì)多幀連通區(qū)域并行標(biāo)記，使標(biāo)記時(shí)間僅為串行標(biāo)記時(shí)間的三分之一。

初始化與當(dāng)前圖像大小相同且全為0的標(biāo)簽矩陣，設(shè)(x,y)為當(dāng)前處理的像素點(diǎn)，g(x,y)為該點(diǎn)像素值，l(x,y)為該點(diǎn)在標(biāo)簽矩陣中的標(biāo)簽值，label為標(biāo)簽矩陣當(dāng)前的最大標(biāo)簽值，其中l(wèi)abel≥2，基于快速連通域標(biāo)記的Two-Pass算法的流程如下所述。

2.1 第一趟掃描

由于第一趟掃描時(shí)標(biāo)簽矩陣全部初始化為0，當(dāng)前像素點(diǎn)后的像素點(diǎn)尚未處理，因此在處理時(shí)只需要考慮當(dāng)前像素點(diǎn)左上方、上方、右上方和左邊像素點(diǎn)的標(biāo)簽值，如圖1所示，深色像素為當(dāng)前像素點(diǎn)。

圖1 八鄰接像素點(diǎn)標(biāo)簽值示意圖

掃描當(dāng)前像素點(diǎn)(x,y)，如果g(x,y)=0，繼續(xù)下一步，否則處理下一個(gè)像素點(diǎn)。

(1)如果l(x-1,y-1)=l(x-1,y)=l(x-1,y+1)=l(x,y-1)=0，即八鄰接像素點(diǎn)標(biāo)簽值全為0，當(dāng)前像素點(diǎn)和已被標(biāo)記的像素點(diǎn)不鄰接(新的連通區(qū)域)，則當(dāng)前最大標(biāo)簽值加1，作為當(dāng)前像素點(diǎn)的標(biāo)簽值，即：

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，八鄰接點(diǎn)像素值相同且不為0，即只與一個(gè)連通域鄰接，則將該點(diǎn)的標(biāo)簽值設(shè)置為鄰接點(diǎn)的標(biāo)簽值，表示如下：

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)中至少存在兩個(gè)相等且不為0的值，即八鄰接點(diǎn)像素值不全相同，當(dāng)前像素點(diǎn)與多個(gè)連通域相連，則將該點(diǎn)標(biāo)簽值設(shè)置為多個(gè)連通域標(biāo)簽值的最小值，內(nèi)容如下所示：

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)作為等價(jià)對(duì)存入映射集合Map中，供第二趟掃描使用，Map內(nèi)容如下所示：

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

2.2 第二趟掃描

在第一趟掃描(3)步驟中，得出等價(jià)關(guān)系集合Map，根據(jù)Map中的等價(jià)關(guān)系，將當(dāng)前像素點(diǎn)標(biāo)簽值更新為l(x,y)，即：

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

由于標(biāo)簽合并，導(dǎo)致標(biāo)簽號(hào)不連續(xù)，重新編號(hào)使標(biāo)簽號(hào)連續(xù)。掃描結(jié)束后，圖像中具有相同標(biāo)簽號(hào)的像素就構(gòu)成一個(gè)連通區(qū)域。

3 巖心圖像空間標(biāo)記

三維空間標(biāo)記以二維連通域合并為基礎(chǔ)，將多幀相鄰圖像連通域合并，形成多個(gè)空間連通體，根據(jù)巖心孔隙連通性規(guī)則，去除孤立的連通體(孤立孔隙)，由于巖心圖像分辨率高且數(shù)量較多，采用多線程技術(shù)加速圖像連通域合并。

3.1 標(biāo)簽合并

設(shè)定Si為上層幀圖像，Si+1為下層幀圖像，Si.l(x,y)為當(dāng)前連通域標(biāo)簽值，下層幀圖像對(duì)應(yīng)像素點(diǎn)的連通域標(biāo)簽值為Si+1.l(x,y)，標(biāo)簽的合并采取小標(biāo)簽值為主原則即存在多個(gè)標(biāo)簽值時(shí)統(tǒng)一用最小的標(biāo)簽值替換其余標(biāo)簽值。

3.1.1 更新下層圖像標(biāo)簽

在對(duì)單幀圖像標(biāo)記時(shí)，設(shè)定上層圖像標(biāo)簽值小于下層標(biāo)簽值。如果Si.l(x,y)不為0(0為背景像素點(diǎn)直接跳過(guò))且Si.l(x,y)所屬的連通域只與上層圖像的一個(gè)連通域相鄰，則將Si+1.l(x,y)所在連通域的值替換為Si.l(x,y)的值，即：

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

若Si+1.l(x,y)所屬的連通域與上層圖像中的多個(gè)連通域相鄰，則將Si+1.l(x,y)所屬的連通域的標(biāo)簽值替換為上層多個(gè)連通域標(biāo)簽值的最小值，即：

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

式中，Si.l(xj,yj),…,Si.l(xk,yk)為上層圖像不同連通域的標(biāo)簽值，如圖2(左)所示。

3.1.2 更新上層圖像標(biāo)簽

若連通域Si+1.l(x,y)與上層圖像的多個(gè)連通域相鄰，即通過(guò)Si+1.l(x,y)使Si中原本不連通的兩個(gè)或多個(gè)連通域連通，則合并這兩個(gè)或多個(gè)連通域，用多個(gè)連通域的最小標(biāo)簽值代替其他標(biāo)簽值。

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

式中，Si+1.l(x,y)為下層連通域標(biāo)簽值，{Si.l(xj,yj)}為上層圖像中與Si+1.l(x,y)相鄰的連通域的標(biāo)簽值集合，用上層連通域中標(biāo)簽值的最小值替換其余連通域標(biāo)簽值。如圖2(右)為更新后上下層圖像連通域的標(biāo)簽值，更新后Si中標(biāo)簽值為3的連通域合并到標(biāo)簽值為2的連通域中。

圖2 標(biāo)簽標(biāo)記(左)與合并圖(右)

3.2 去除孤立孔隙

數(shù)字巖心中有效孔隙為相互連通的，在一般壓力條件下，可允許液體在其中流動(dòng)的孔隙，除此之外為孤立孔隙(封閉孔隙)。即不屬于起始和末尾幀圖像且不與起始、末尾幀及多幀圖像邊界孔隙像素相連的孔隙為孤立孔隙。設(shè)Ubound為圖像邊界孔隙像素點(diǎn)的標(biāo)簽值和起始、末尾幀圖像孔隙像素點(diǎn)的標(biāo)簽值集合，集合U為多幀圖像連通體標(biāo)簽值集合，Uremove為孤立連通體標(biāo)簽值集合。

3.3 連通域合并優(yōu)化

為簡(jiǎn)化圖像標(biāo)簽合并，在對(duì)單幀圖像標(biāo)簽時(shí)保證Si+1幀的標(biāo)簽值大于Si幀的標(biāo)簽值，Si幀標(biāo)簽的起始值為圖像編號(hào)乘以大于所有單幀圖像連通域個(gè)數(shù)的定值N，即Si×N。采用多線程并行合并技術(shù)，開(kāi)始合并前先將所有幀圖像按照層次分組，并設(shè)計(jì)數(shù)據(jù)結(jié)構(gòu)SliceGroup存儲(chǔ)合并所需的參數(shù)。

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

圖像數(shù)量slices表明合并的層次，規(guī)模為2n(n=0,1,2…)，若共有2n幀圖像，共需要進(jìn)行n層合并。G1、G2為同一層中待合并的兩組幀圖像，G1.Si+slices為G1中的最后一幀圖像，G2.Si為G2中的第一幀圖像，其中G1.(Si+slices).l(x,y)

圖3 合并前(左)與合并后(右)

4 閾值優(yōu)化

去除孤立孔隙后，重新計(jì)算圖像的孔隙度并與實(shí)驗(yàn)孔隙度對(duì)比，若誤差超出限定范圍，則調(diào)整閾值并二值化重復(fù)步驟2、3，直至誤差縮小到限定范圍以內(nèi)。設(shè)f(x)為孔隙度計(jì)算函數(shù)，porosity為實(shí)驗(yàn)孔隙度，T為復(fù)雜度分割閾值?？紫抖萬(wàn)(x)計(jì)算公式為：

(2)

式中，AllHolePixel為總的孔隙像素個(gè)數(shù)，AllPixel為總的像素個(gè)數(shù)。閾值調(diào)整方法為：

(3)

5 實(shí)驗(yàn)分析

該文采用的實(shí)驗(yàn)平臺(tái)為Intel(R) Xeon(R) CPU E5-1620，主頻為3.70 GHz，系統(tǒng)內(nèi)存為16 GB，顯卡為NVIDIA Quadro K2000，顯卡內(nèi)存為2 GB，操作系統(tǒng)為Windows 10，實(shí)驗(yàn)開(kāi)發(fā)環(huán)境為Visual Studio 2015。實(shí)驗(yàn)圖像為650張分辨率為968*995吉林砂巖CT圖像。

5.1 多層切片標(biāo)記效率分析

為了驗(yàn)證基于多線程加速技術(shù)的有效性，任意選取不同編號(hào)的圖像，采用八鄰接搜索方式進(jìn)行判別，分別測(cè)試普通Two-Pass掃描法、種子填充算法、基于并查集的Two-Pass掃描法的運(yùn)行時(shí)間，并對(duì)650幀圖像在多線程下標(biāo)記時(shí)間與串行標(biāo)記時(shí)間進(jìn)行對(duì)比，對(duì)比結(jié)果如表1所示。

表1 圖像標(biāo)記效率對(duì)比

通過(guò)以上數(shù)據(jù)分析可知，三種算法在處理單張圖像時(shí)相差不明顯，基于并查集的Two-Pass算法速度最快，種子填充算法的速度最慢。對(duì)多幀圖像進(jìn)行標(biāo)記時(shí)，通過(guò)多線程技術(shù)并行處理圖像，時(shí)間為串行處理時(shí)間的三分之一左右。說(shuō)明通過(guò)增加并查集和多線程技術(shù)，有效提高了圖像的標(biāo)記效率。

5.2 多層切片連通域優(yōu)化分析

為了驗(yàn)證空間標(biāo)記分割的有效性，選取同一組圖像，根據(jù)孔隙度進(jìn)行優(yōu)化分割，選取多張圖像對(duì)比分割效果如圖4所示。通過(guò)對(duì)多幀圖像連通域標(biāo)記、合并、去除孤立孔隙，根據(jù)理論孔隙度優(yōu)化分割閾值，孔隙度和閾值隨著迭代次數(shù)增加變化如表2所示。

表2 孔隙度和閾值變化

由表2可知，未經(jīng)過(guò)連通域優(yōu)化(迭代次數(shù)為0)與經(jīng)過(guò)連通域優(yōu)化的圖像雖然具有相同的分割閾值，但未優(yōu)化的孔隙度大于優(yōu)化后的孔隙度，表明非連通孔隙點(diǎn)對(duì)孔隙度有一定的影響。隨著迭代次數(shù)增加分割閾值增大，孔隙度逐漸增大，但始終圍繞著0.12上下波動(dòng)，閾值在103～104之間變動(dòng)，通過(guò)設(shè)置最小誤差，最終確定分割閾值。

如圖4所示，第一行為原始幀圖像，圖中圓域內(nèi)像素值較小(顏色較深)為孔隙，像素值較大(顏色較淺)為骨架，圖中第二行為Otsu分割的二值圖像，沒(méi)有考慮圖像間上下層像素點(diǎn)的鄰接關(guān)系，第三行為考慮了幀與幀之間的聯(lián)系。通過(guò)對(duì)比可知，未考慮幀與幀之間像素關(guān)聯(lián)二值圖像有更多的孤立點(diǎn)和噪聲，去除了孤立孔隙的二值圖像黑色像素點(diǎn)明顯增多，孔隙輪廓更加明顯。由此可知，通過(guò)增加了幀與幀之間關(guān)聯(lián)的空間標(biāo)記方法有效去除圖像的孤立孔隙點(diǎn)。

圖4 二值圖像對(duì)比圖

5.3 多層切片連通域合并分析

為了驗(yàn)證圖像合并的有效性和準(zhǔn)確性，選取跟蹤了中間位置幀圖像0240.tiff圖像的特定標(biāo)簽值觀察連通域合并過(guò)程。

圖5 連通域合并跟蹤圖

如圖5所示，未進(jìn)行合并時(shí)標(biāo)簽值連通區(qū)域面積較小，隨著合并迭代次數(shù)的增加圖像中標(biāo)簽值所代表的連通區(qū)域面積逐漸增加，表明單幀圖像中不連通的區(qū)域，隨著多幀圖像標(biāo)簽的合并最終連接在一起，形成一個(gè)連通體。剔除孤立孔隙后圖像中只剩與外界相同的連通體，如圖5中最后一張圖片所示。

5.4 數(shù)字巖心建模對(duì)比分析

為了驗(yàn)證通過(guò)連通域標(biāo)記算法對(duì)巖心圖像進(jìn)行三維分割的有效性，分別選取了10、100、200、500張圖像生成數(shù)字巖心模型，并比較效果，如圖6所示。

圖6 數(shù)字巖心孔隙模型對(duì)比圖

文獻(xiàn)[9]提出的半自動(dòng)分割方法，需要手動(dòng)對(duì)圖像分割結(jié)果進(jìn)行修復(fù)，由于不同的實(shí)驗(yàn)者具有不同的主觀判斷無(wú)法做到完全一致，造成實(shí)驗(yàn)結(jié)果具有不一致性。其他文獻(xiàn)方法皆為對(duì)單幀圖像的處理與文中方法有本質(zhì)區(qū)別。文中提出的方法完全實(shí)現(xiàn)自動(dòng)化，避免了人為因素的影響且考慮了巖心三維空間聯(lián)系，因此不具有對(duì)比參考性，故不進(jìn)行對(duì)照實(shí)驗(yàn)。

6 結(jié)束語(yǔ)

針對(duì)數(shù)字巖心三維建模分離骨架孔隙只考慮圖像二維信息而忽略圖像間的空間聯(lián)系，導(dǎo)致生成的三維巖心模型不準(zhǔn)確，孔隙連通情況不精確的問(wèn)題，提出了基于二維連通域分析并結(jié)合多線程技術(shù)的空間標(biāo)記并行分割方法。采用Otsu算法預(yù)分割巖心圖像，通過(guò)對(duì)二值圖像連通域合并，生成多個(gè)互不相通的連通體，根據(jù)邊界限定準(zhǔn)則去除孤立的孔隙，并根據(jù)理論孔隙度優(yōu)化分割閾值。文中使用了650張吉林砂巖的CT圖像，分別使用簡(jiǎn)單閾值分割方法和基于空間標(biāo)記的串行分割方法作為對(duì)比。實(shí)驗(yàn)表明，基于空間標(biāo)記的并行分割方法有效去除了圖像中孤立的孔隙像素和噪聲像素，提高了巖心模型的準(zhǔn)確度和分割速度，為后續(xù)建立數(shù)字巖心三維模型奠定了基礎(chǔ)。下一步將繼續(xù)優(yōu)化標(biāo)簽合并速度。