基于空間標記的巖心孔隙并行分割

陳國軍，李 勝

(中國石油大學(華東) 計算機科學與技術(shù)學院，山東 青島 266580)

0 引 言

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

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

1 圖像預(yù)處理

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

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

圖像二值化公式為：

(1)

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

2 二值圖像連通域標記

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

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

2.1 第一趟掃描

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

圖1 八鄰接像素點標簽值示意圖

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

(1)如果l(x-1,y-1)=l(x-1,y)=l(x-1,y+1)=l(x,y-1)=0，即八鄰接像素點標簽值全為0，當前像素點和已被標記的像素點不鄰接(新的連通區(qū)域)，則當前最大標簽值加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，即只與一個連通域鄰接，則將該點的標簽值設(shè)置為鄰接點的標簽值，表示如下：

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的值，即八鄰接點像素值不全相同，當前像素點與多個連通域相連，則將該點標簽值設(shè)置為多個連通域標簽值的最小值，內(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)作為等價對存入映射集合Map中，供第二趟掃描使用，Map內(nèi)容如下所示：

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

2.2 第二趟掃描

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

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

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

3 巖心圖像空間標記

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

3.1 標簽合并

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

3.1.1 更新下層圖像標簽

在對單幀圖像標記時，設(shè)定上層圖像標簽值小于下層標簽值。如果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 去除孤立孔隙

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

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

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

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

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

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

4 閾值優(yōu)化

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

(2)

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

(3)

5 實驗分析

該文采用的實驗平臺為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，實驗開發(fā)環(huán)境為Visual Studio 2015。實驗圖像為650張分辨率為968*995吉林砂巖CT圖像。

5.1 多層切片標記效率分析

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

表1 圖像標記效率對比

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

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

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

表2 孔隙度和閾值變化

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

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

圖4 二值圖像對比圖

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

為了驗證圖像合并的有效性和準確性，選取跟蹤了中間位置幀圖像0240.tiff圖像的特定標簽值觀察連通域合并過程。

圖5 連通域合并跟蹤圖

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

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

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

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

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

6 結(jié)束語

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