基于聚類分析的管路圖像亞像素邊緣提取算法

王 驍，劉檢華，劉少麗，金 鵬，吳天一

(北京理工大學 機械與車輛學院數(shù)字化制造研究所，北京 100081)

1 問題的提出

管路系統(tǒng)廣泛應用在航空航天產(chǎn)品中，它是壓力系統(tǒng)、動力系統(tǒng)、冷卻系統(tǒng)和控制系統(tǒng)等重要的組成部分，各系統(tǒng)的貯箱、閥門和發(fā)動機等零件都由管路連接和中轉。這些系統(tǒng)的壽命、穩(wěn)定性和可靠性由管路系統(tǒng)的性能決定[1]。管路具有復雜的外形和走向，為了保證生產(chǎn)符合設計需求的管路，在制造過程中需要進行測量，以驗證管路是否合格，從而保證精確裝配和無應力裝配[2]。在測量管路零件的方法中，視覺測量是一種快速有效的測量方法[3]，其主要思路是重建管路的三維模型，從模型中獲取管路有關尺寸。但是管路零件外形均勻、表面缺少特征點，難以利用匹配特征點的方法進行三維重建，因此大量基于視覺的管路重建方法[4-8]都利用了圖像中管路邊緣特征。

為了凸顯管路邊緣通常采用背光光源照射，但是由于光源光線不均勻，在提取邊緣時容易將光源中不均勻部分誤提取。同時，管路為金屬材質(zhì)，表面經(jīng)過處理具有鏡面反射的特征，又由于管路表面為圓柱曲面，導致反光區(qū)域不均勻，而且隨著管路擺放位置不同，反光區(qū)域也會不同。此外自然光干擾和反光區(qū)域與背光光源相近等問題都影響了對管路邊緣的提取。如圖1a所示，由于管路反光區(qū)域和背光區(qū)域相近，導致邊緣不明顯；圖1b中管路反光區(qū)域復雜，管路區(qū)域內(nèi)部明暗不均處容易被當邊緣提取。

由于數(shù)字圖像由像素點構成，利用圖像中梯度變化，通過計算梯度的微分算子可以獲得圖像中邊緣特征。但是這種計算只能獲得像素精度的邊緣，即不連續(xù)、不平滑的邊緣，精度只到像素精度。應用這樣的邊緣重建的管路模型精度無法達到測量需求，為此需要利用亞像素精度的邊緣。亞像素精度的邊緣可以讓提取的邊緣精度在一個像素以內(nèi)，在獲取亞像素精度時通常利用了圖像中的其他信息。提取亞像素邊緣算法通常利用了邊緣附近的灰度變化[23]，通過擬合灰度變化曲面或?qū)ο袼鼗叶茸兓逯档仁侄潍@取具有亞像素精度的邊緣。

常見的亞像素邊緣提取方法可以分為矩法、插值法和擬合法。矩法是利用模板算子，以邊緣點附近灰度值計算求解亞像素邊緣。Tabatabai[9]提出一種基于灰度矩的亞像素邊緣提取算法。Ghosal[10]提出基于Zernike矩的亞像素邊緣提取方法，利用3個算子提取了亞像素邊緣。Qu[11]改進了Zernike矩方法，提出一種快速Zernike矩法，提升了計算效率。Bin[12]提出了基于Fourier-Mellin矩的方法，提取了亞像素邊緣，使直線上亞像素精度達到0.16個像素尺寸，曲線上亞像素精度達到0.23個像素尺寸，提高了矩法的精度。但是矩法容易受到噪聲干擾，在提取管路圖像邊緣時效果不佳，不適合用于管路圖像亞像素邊緣提取。

插值法則是在邊緣點附近根據(jù)灰度值進行插值計算，獲取亞像素精度邊緣。Devernay[13]利用非極大值抑制實現(xiàn)了亞像素邊緣提取，Rockett[14]擬合二次多項式插值求解亞像素邊緣，Xie[15]采用Hermite插值求解亞像素邊緣，這些方法都采用了Canny檢測算子求解圖像中的亞像素邊緣。此外，Kubota[16]提出一種雙線性插值方法提高圖像分辨率從而提取亞像素精度邊緣。Hermosilla[17]提出了一種根據(jù)本質(zhì)無震蕩格式的4階非線性插值方法求解亞像素邊緣。插值求解亞像素邊緣運算速度快效率高，但是抗噪性差，在管路圖像中噪聲干擾多，難以獲得準確穩(wěn)定的結果。

擬合法是根據(jù)邊緣部分灰度值的變化曲線，擬合參數(shù)模型從而定位亞像素精度邊緣。可以利用二次曲線[18]，B樣條曲線[19]，Sigmoid函數(shù)[20]或者能量泛函模型[21]擬合求解亞像素邊緣。也有對圖像中邊緣區(qū)域采用高斯模糊，利用邊緣區(qū)域灰度值擬合高斯模糊曲面求解邊緣，達到亞像素精度[22]。這些方法將亞像素邊緣求解問題轉換為非線性尋優(yōu)問題，通過擬合邊緣區(qū)域灰度變化求解邊緣區(qū)域連續(xù)變化函數(shù)，從而求解亞像素邊緣。但是尋優(yōu)計算通常為無初值計算，消耗大量計算時間，而且容易得到局部最優(yōu)解。Agustin[23]提出的一種最新擬合方法，計算效率高，結果準確。但是在求解管路圖像亞像素邊緣時都不能達到最好的效果，仍然存在提取錯誤或者提取失敗的情況。

但是已有方法在提取管路邊緣時容易受到管路光照條件的干擾，導致漏提取、多提取、錯提取甚至無法提取的情況。為了提取準確、完整的亞像素邊緣，本文提出一種針對管路圖像提取亞像素精度邊緣的方法，先利用頻域濾波濾除圖像中低頻噪聲并凸顯邊緣等高頻信號，然后利用聚類分析準確分割管路圖像區(qū)域，并根據(jù)分割結果準確定位求解亞像素邊緣的區(qū)域，最后利用求解區(qū)域內(nèi)灰度值擬合灰度連續(xù)變化曲面，尋找曲面梯度最大點作為邊緣，從而獲得亞像素精度邊緣。本文最大的優(yōu)勢是利用圖像處理等一系列手段自動篩選了管路邊緣點，排除了由于噪聲等問題引起的錯誤提取點，為重建管路三維模型提供了基礎。經(jīng)過實驗驗證，本文方法提取管路亞像素邊緣準確率達到99%以上，且達到0.04個像素單位的精度，也可用于管路三維重建。

2 處理方法

本文提出的針對管路圖像亞像素精度邊緣提取方法流程如圖2所示，首先頻域濾波與區(qū)域粗分割獲得濾除低頻噪聲的管路圖像，并將管路圖像區(qū)域連同周圍部分背景分割出來，然后通過聚類分析對粗分割管路圖像區(qū)域進行細致分割，獲取精確的管路區(qū)域，再利用圖像形態(tài)學提取像素精度邊緣，用于亞像素精度邊緣的初值計算和最后的亞像素精度邊緣提取計算。

2.1 頻域濾波和區(qū)域粗分割

利用多目視覺系統(tǒng)拍攝的管路圖像具有多種噪聲，這些噪聲中存在很多低頻噪聲。本文采用頻域高通濾波器，濾除了圖像中不同的低頻噪聲信號，同時達到了銳化管路邊緣的目的，其過程如圖3所示。

首先，利用傅里葉變換將圖像從空域轉換為頻域，如式(1)所示，

(1)

其中u和v表示頻域空間中點的坐標，F(xiàn)(u,v)表示了對應點在頻域空間中的值;x和y表示圖像中的點坐標，f(x,y)表示了圖像中對應像素點的灰度值。

然后利用高通濾波器濾除圖像中的低頻噪聲。濾波后，還需要通過傅里葉逆變換將頻域圖像轉換回空域圖像，如式(2)所示。

(2)

式中：f表示還原后的圖像，F(xiàn)表示頻域圖像，M和N分別表示原圖像縱橫像素數(shù)。

對濾波后的圖像進行簡單的閾值分割，這樣的分割結果容易包含非管路區(qū)域，或者在分割的管路區(qū)域內(nèi)出現(xiàn)孔洞。因此，利用圖像形態(tài)學膨脹管路區(qū)域，將孔洞或者非管路區(qū)域包含在膨脹區(qū)域內(nèi)，從而獲取包含管路區(qū)域的粗分割結果，如圖3中所示。

2.2 聚類分析分割圖像

對管路的粗分割區(qū)域可以進行細致分割。如圖4所示，根據(jù)圖像中像素灰度特征，可以將粗分割的管路區(qū)域分為背景區(qū)域、邊緣區(qū)域、反光區(qū)域和背光區(qū)域。各區(qū)域之間的界線距離較近，反光區(qū)域和背光區(qū)域復雜交錯，采用閾值分割方法容易產(chǎn)生錯誤提取，導致結果不夠準確。為了準確提取管路區(qū)域，本文采用聚類分析細致分割圖像，提取完整準確的管路區(qū)域，從而獲取用于亞像素邊緣計算的初值。

利用聚類分析可以自動將管路粗分割結果分為如圖4所示的4個區(qū)域，本文采用模糊C均值聚類(Fuzzy C-Means, FCM)算法計算。對管路粗分割結果中的像素xi進行初始化，任取點作為聚類中心，記為cj(j=1,2,3,4)。每個xi可計算系數(shù)uij，代表評定xi屬于第j類的測量值，如式(3)表示。

(3)

式中：m表示大于1的實數(shù)作為計算的因子。聚類中心cj隨著uij變化而變化，如式(4)所示。

(4)

構造目標函數(shù)Jm如式(5)所示，以Jm的極小值為優(yōu)化目標，通過迭代計算完成聚類計算。迭代計算的終止條件如式(6)所示。

(5)

‖U(k+1)-U(k)‖<ε。

(6)

式(6)中ε∈(0,1)，表示終止條件，當兩次迭代結果中U的差值小于ε時可以認為達到迭代終止條件。FCM算法步驟如下：

步驟1初始化U，標記為U(0)，并開始迭代計算。

步驟2計算k+1次迭代計算中的各聚類中心，記為C(k+1)。

步驟3根據(jù)中心C(k+1)，計算目標函數(shù)極值解，計算U(k+1)。

步驟4比較當前U(k)和上一輪迭代U(k-1)的差值，判斷終止還是進入下次迭代。

通過計算，圖像被劃分為4個區(qū)域，如圖5所示。

采用FCM算法可以對管路圖像進行細致分割，提取了完整的管路區(qū)域，管路區(qū)域包含了反光區(qū)域、背光區(qū)域和邊緣區(qū)域。實際上FCM根據(jù)灰度值計算了4個聚類中心，將圖像分割為4個區(qū)域，其中背景區(qū)域和邊緣區(qū)域是進行亞像素邊緣提取的主要計算區(qū)域。

2.3 提取像素精度邊緣

在提取亞像素邊緣前，準確定位像素邊緣可以有效降低漏提取的情況，也可以為提取亞像素邊緣提供計算初值。對于管路圖像，提取像素邊緣需要滿足如下兩個需求：①單像素邊緣；②分割區(qū)域的外層邊緣。經(jīng)典的像素邊緣提取算法主要是構建微分算子，與圖像卷積運算，設置梯度閾值求取極大閾值檢測邊緣。此外，還有根據(jù)數(shù)學形態(tài)學，移除區(qū)域內(nèi)部像素提取邊緣的算法。通過驗證，數(shù)學形態(tài)學方法更適合本文應用場合。

2.3.1 微分算子提取邊緣

檢測圖像邊緣就是檢測灰度變化劇烈的部分，通常采用圖像灰度梯度變化檢測邊緣。對于圖像中各像素點的梯度如式(7)所示：

(7)

梯度幅值變化為：

‖

(8)

邊緣切線方向斜率可以表示為：

(9)

其中I(x,y)表示像素點(x,y)的灰度值，對像素點各方向求偏導值求取梯度I(x,y)，利用[Ix,Iy]表示，最大梯度方向的幅值為‖I(x,y)‖。

根據(jù)式(7)構造微分算子檢測邊緣，常見的微分算子有Prewitt算子、Roberts算子和Sobel算子[24]，如圖6所示。Canny算子[24]也是一種檢測邊緣的算法，相較于其他算子，Canny算子抗噪性突出，利用高斯核函數(shù)構造算子濾除噪聲，計算梯度極值區(qū)域后進行非極大抑制處理，降低非極值當做邊緣提取的情況，再利用之后閾值處理，利用雙閾值求取邊緣并細化獲得最終值。Canny算子提取邊緣考慮了更多因素，應用更加廣泛。LoG算子提取邊緣也考慮抗噪問題，通過高斯核函數(shù)濾除了大部分噪聲。

2.3.2 數(shù)學形態(tài)學提取邊緣

利用形態(tài)學將目標區(qū)域像素分為邊緣像素和內(nèi)部像素，移除內(nèi)部像素保留邊緣像素提取邊緣。通常將目標區(qū)域內(nèi)的像素置1，區(qū)域外像素置0，提取區(qū)域內(nèi)3×3區(qū)域的像素模板，當中心像素和與其相鄰的4個像素都為1時說明中心像素為內(nèi)部像素，將其置0即完成移除，對目標區(qū)域內(nèi)所有像素都進行處理后即完成邊緣提取，其主要流程如圖7所示。

為了提取符合前文所述的兩個標準，利用不同的方法提取像素精度邊緣，得到如圖8所示的結果。

圖中白色區(qū)域表示管路區(qū)域，淺灰色區(qū)域表示超出管路區(qū)域提取的邊緣，深灰色表示在管路區(qū)域內(nèi)提取的邊緣。對比結果，各方法都成功提取了單像素邊緣，但是利用微分算子提取的邊緣容易超出管路區(qū)域，這是因為利用微分算子卷積計算時按照圖像從左到右從上到下的順序，只根據(jù)梯度提取像素精度邊緣，忽略了邊緣所在區(qū)域。而根據(jù)圖像形態(tài)學提取的邊緣都是在管路區(qū)域內(nèi)，保證了邊緣從區(qū)域內(nèi)到區(qū)域外的梯度方向。從提取精度而言，出現(xiàn)超出管路區(qū)域邊緣并不是提取錯誤，但是為了結合后續(xù)操作提取亞像素精度邊緣，推薦采用形態(tài)學的方法提取邊緣，并計算梯度法線方向作為亞像素精度邊緣計算初值。

2.4 亞像素邊緣提取

為了得到準確的亞像素邊緣，需要根據(jù)邊緣附近局部灰度值，擬合灰度連續(xù)變化曲面，根據(jù)局部灰度閾值求解亞像素邊緣。以像素精度邊緣上一點為中心，其n×n區(qū)域內(nèi)灰度變曲面作為亞像素邊緣求解區(qū)域，灰度值變化如圖9所示。這樣的變化趨勢與高斯概率密度函數(shù)相似，因此本文考慮利用高斯概率密度函數(shù)作為邊緣連續(xù)變化曲面，擬合求解亞像素邊緣。采用高斯概率密度函數(shù)有兩個優(yōu)勢：①相較于高斯模糊函數(shù)擬合[21]，不需要對管路圖像進行高斯模糊處理，避免了因為高斯模糊導致的圖像邊緣區(qū)域丟失的情況。如圖10所示，管路圖像邊緣區(qū)域與背光區(qū)域邊界距離較近，利用高斯模糊處理后容易將兩個區(qū)域并在一起，這容易導致提取的亞像素邊緣不準確，以背光區(qū)域的邊緣作為管路邊緣錯誤提取。②擬合初值容易計算，利用計算區(qū)域內(nèi)各像素點灰度值統(tǒng)計結果即可計算大部分初值用于擬合尋優(yōu)計算，提高了計算效率。

根據(jù)圖8反映的求解區(qū)域灰度變化離散值，擬合高斯概率密度函數(shù)，求解亞像素精度邊緣，其解析式如式(10)所示，

(10)

式中G(x,y)表示圖像中點(x,y)對應的灰度值，a、b、μ、σ和k為需要擬合求解的參數(shù)。a表示求解區(qū)域灰度最小值，b表示求解區(qū)域灰度最大值與最小值的差值，μ和σ為高斯概率密度函數(shù)參數(shù)，k表示邊緣點在圖像中的切線方向。以Is為求解區(qū)域，擬合目標函數(shù)如式(11)所示。

minx,y∈Is∑‖G(x,y)-Gr(x,y)‖2。

(11)

式中G(x,y)表示擬合曲面計算值，Gr(x,y)表示實際灰度值。

擬合灰度連續(xù)變化曲面是非線性尋優(yōu)問題，利用Levenberg-Marquardt法可求解目標函數(shù)全局最優(yōu)解[25]。為了獲得穩(wěn)定收斂的求解結果，需要選取合適的初值進行尋優(yōu)計算。利用本文提出的亞像素邊緣計算方法，可以根據(jù)計算區(qū)域的灰度值統(tǒng)計結果求解可用于尋優(yōu)計算的初值。選取n×n的計算區(qū)域Is，各點灰度值為Gr(x,y)，其中x,y∈[1,n]。邊緣切線斜率k初值通過式(9)求解，其余參數(shù)可以通過統(tǒng)計求解區(qū)域的灰度值設置初值，各參數(shù)初值求解如式(12)～式(15)所示：

a=minGr;

(12)

b=maxGr-a;

(13)

(14)

(15)

為了獲得準確的擬合結果，需選取大小適中的計算區(qū)域。導管邊緣區(qū)域比較狹長，為了保證精度，推薦選取7×7或5×5的局部區(qū)域用于擬合計算。當選取的計算區(qū)域過大時，會受到反光、背光等區(qū)域的干擾，導致無法獲得收斂結果。而擬合的區(qū)域過小，會受到噪聲影響，導致擬合結果不準確。通過多次實驗，選取7×7或5×5的計算區(qū)域能獲得準確穩(wěn)定的結果，當管路截面直徑較小時(如5～8 mm)，由于管路邊緣區(qū)域太小，選取5×5的計算區(qū)域更合適。

利用灰度連續(xù)變化函數(shù)，可以求解亞像素精度邊緣以及其對應的梯度方向。設置一個計算區(qū)域內(nèi)的灰度閾值即可定位準確的亞像素邊緣點，這個閾值可通過聚類分析的結果獲得，記為Gt。計算如式(16)所示，cb表示背景區(qū)灰度值中心，ce表示邊緣區(qū)域灰度值中心。ω為權值參數(shù)，經(jīng)過計算和實驗一般可設置在0.2～0.5之間。

Gt=ωcb+(1-ω)ce。

(16)

將閾值Gt代入式(10)即可得到亞像素邊緣的切線方程，如式(17)所示。

(17)

(18)

根據(jù)梯度方向[-A-B]T，可以求得亞像素邊緣點Esub坐標為

(19)

3 管路亞像素邊緣提取實驗

驗證實驗采用AVT F146b型相機，Computar 8mm鏡頭拍攝，圖像像素為1280×960，采用背光光源照射管路。為了說明本文方法能在提取管路亞像素邊緣的同時自動篩選提取結果，排除錯誤的邊緣點，利用5個不同的管路進行實驗，并通過統(tǒng)計結果說明本文方法準確有效。此外，為了說明本文方法精度，設計了精度實驗，求解像素標準圓的圓心坐標偏差和標準方形的邊長偏差。

本文方法提取的亞像素邊緣結果如圖12a和圖13a所示，能篩選出需要的管路邊緣點，而不受到復雜光照環(huán)境的影響。為了顯示光照對提取邊緣的影響，對比了直接使用的Agustin方法提取的邊緣結果[11]。從結果中可知，圖12中被拍攝管路的邊緣不明顯，容易提取失敗，因此利用Agustin方法提取的邊緣不完整；圖13中被拍攝管路的邊緣與背光部分邊緣較近，容易受到干擾，因此使用Agustin方法容易造成提取錯誤。而本文方法因為設計了能自動篩選干擾點和錯誤點的方法，因此可以克服光照環(huán)境中的噪聲點。

通過對比實驗可知，本文針對管路圖像的亞像素邊緣提取方法能消除管路中復雜光照環(huán)境帶來的噪聲干擾，對比Agustin方法不需要額外的篩選手段，即可自動篩選管路邊緣亞像素邊緣點，并消除反光位置的邊界干擾。為了進一步量化實驗結果，以提取的像素精度邊緣點數(shù)作為提取結果的標準值評判提取結果的準確程度。統(tǒng)計提取的亞像素精度邊緣總數(shù)，再減去提取錯誤的情況。與Agustin方法對比提取結果，并通過其他手段對Agustin方法提取的結果進行篩選，選取實際管路圖像進行實驗，結果如表1所示。

表1 管路邊緣提取對比

表1中，提取邊緣表示算法直接提取的亞像素邊緣點數(shù)，通過判斷篩選出正確的提取結果，再以篩選結果和標準值計算準確率。對比提取結果，本文方法能不受干擾提取準確的管路亞像素邊緣點，而且具有更高的準確率。而Agustin的方法容易受到管路圖像中噪聲影響，導致提取錯誤。通過對比可知，本文方法相較于Agustin方法在管路亞像素邊緣提取時具有更高的準確率，在計算亞像素邊緣點前的一系列處理可以保證獲得用于重建計算的亞像素邊緣點，而不受到管路復雜光照環(huán)境的影響。

為了驗證本文亞像素邊緣提取精度，采用圓形和方向邊緣提取算例，提取結果如圖14所示。圖14a為直徑20個像素的圓形，圖14b為邊長20像素的正方形。通過本文方法提取的亞像素邊緣，計算圖14a圓心偏差以及圖14b方形邊長誤差，結果如表2所示。

表2 亞像素邊緣提取誤差對比

圓心誤差XY邊長誤差0.030 40.029 30.046 1

從提取結果和誤差對比結果中可以看出，本文方法提取的亞像素邊緣精度達到0.04個像素尺寸，可滿足于管路重建與測量。

4 結束語

本文利用提取的管路亞像素精度邊緣重建了管路三維模型，測量了管路尺寸。為了保證測量結果準確，需要高精度的亞像素邊緣信息。而獲取亞像素邊緣時容易受到噪聲干擾，造成邊緣缺失甚至提取錯誤邊緣。本文提出了一種針對管路圖像，提取亞像素精度邊緣的方法，能在對光照環(huán)境復雜的管路圖像進行處理，獲取亞像素精度邊緣，用于管路的三維重建和測量。

為了提取精準的管路圖像亞像素邊緣，本文首先采用高通濾波器濾除圖像中低頻信號，保留如邊緣高頻信號。然后，利用聚類分析劃分管路區(qū)域，提取完整且準確的管路區(qū)域。再利用圖像形態(tài)學方法獲取亞像素邊緣求解區(qū)域。進一步，根據(jù)求解區(qū)域的局部灰度變化擬合高斯概率密度函數(shù)作為灰度連續(xù)變化參數(shù)模型。最后，根據(jù)之前聚類分析結果求解計算區(qū)域內(nèi)閾值，定位亞像素精度邊緣點，完成邊緣提取。

通過實驗驗證，本文方法能提取完整、準確和高精度的亞像素邊緣，為基于圖像的管路模型重建提供了精度基礎。今后會根據(jù)提取的亞像素精度邊緣重建準確的管路三維模型，用于航空航天管路的快速、高精度測量。