基于聚类分析的管路图像亚像素边缘提取算法

王 骁，刘检华，刘少丽，金 鹏，吴天一

(北京理工大学 机械与车辆学院数字化制造研究所，北京 100081)

1 问题的提出

管路系统广泛应用在航空航天产品中，它是压力系统、动力系统、冷却系统和控制系统等重要的组成部分，各系统的贮箱、阀门和发动机等零件都由管路连接和中转。这些系统的寿命、稳定性和可靠性由管路系统的性能决定[1]。管路具有复杂的外形和走向，为了保证生产符合设计需求的管路，在制造过程中需要进行测量，以验证管路是否合格，从而保证精确装配和无应力装配[2]。在测量管路零件的方法中，视觉测量是一种快速有效的测量方法[3]，其主要思路是重建管路的三维模型，从模型中获取管路有关尺寸。但是管路零件外形均匀、表面缺少特征点，难以利用匹配特征点的方法进行三维重建，因此大量基于视觉的管路重建方法[4-8]都利用了图像中管路边缘特征。

为了凸显管路边缘通常采用背光光源照射，但是由于光源光线不均匀，在提取边缘时容易将光源中不均匀部分误提取。同时，管路为金属材质，表面经过处理具有镜面反射的特征，又由于管路表面为圆柱曲面，导致反光区域不均匀，而且随着管路摆放位置不同，反光区域也会不同。此外自然光干扰和反光区域与背光光源相近等问题都影响了对管路边缘的提取。如图1a所示，由于管路反光区域和背光区域相近，导致边缘不明显；图1b中管路反光区域复杂，管路区域内部明暗不均处容易被当边缘提取。

由于数字图像由像素点构成，利用图像中梯度变化，通过计算梯度的微分算子可以获得图像中边缘特征。但是这种计算只能获得像素精度的边缘，即不连续、不平滑的边缘，精度只到像素精度。应用这样的边缘重建的管路模型精度无法达到测量需求，为此需要利用亚像素精度的边缘。亚像素精度的边缘可以让提取的边缘精度在一个像素以内，在获取亚像素精度时通常利用了图像中的其他信息。提取亚像素边缘算法通常利用了边缘附近的灰度变化[23]，通过拟合灰度变化曲面或对像素灰度变化插值等手段获取具有亚像素精度的边缘。

常见的亚像素边缘提取方法可以分为矩法、插值法和拟合法。矩法是利用模板算子，以边缘点附近灰度值计算求解亚像素边缘。Tabatabai[9]提出一种基于灰度矩的亚像素边缘提取算法。Ghosal[10]提出基于Zernike矩的亚像素边缘提取方法，利用3个算子提取了亚像素边缘。Qu[11]改进了Zernike矩方法，提出一种快速Zernike矩法，提升了计算效率。Bin[12]提出了基于Fourier-Mellin矩的方法，提取了亚像素边缘，使直线上亚像素精度达到0.16个像素尺寸，曲线上亚像素精度达到0.23个像素尺寸，提高了矩法的精度。但是矩法容易受到噪声干扰，在提取管路图像边缘时效果不佳，不适合用于管路图像亚像素边缘提取。

插值法则是在边缘点附近根据灰度值进行插值计算，获取亚像素精度边缘。Devernay[13]利用非极大值抑制实现了亚像素边缘提取，Rockett[14]拟合二次多项式插值求解亚像素边缘，Xie[15]采用Hermite插值求解亚像素边缘，这些方法都采用了Canny检测算子求解图像中的亚像素边缘。此外，Kubota[16]提出一种双线性插值方法提高图像分辨率从而提取亚像素精度边缘。Hermosilla[17]提出了一种根据本质无震荡格式的4阶非线性插值方法求解亚像素边缘。插值求解亚像素边缘运算速度快效率高，但是抗噪性差，在管路图像中噪声干扰多，难以获得准确稳定的结果。

拟合法是根据边缘部分灰度值的变化曲线，拟合参数模型从而定位亚像素精度边缘。可以利用二次曲线[18]，B样条曲线[19]，Sigmoid函数[20]或者能量泛函模型[21]拟合求解亚像素边缘。也有对图像中边缘区域采用高斯模糊，利用边缘区域灰度值拟合高斯模糊曲面求解边缘，达到亚像素精度[22]。这些方法将亚像素边缘求解问题转换为非线性寻优问题，通过拟合边缘区域灰度变化求解边缘区域连续变化函数，从而求解亚像素边缘。但是寻优计算通常为无初值计算，消耗大量计算时间，而且容易得到局部最优解。Agustin[23]提出的一种最新拟合方法，计算效率高，结果准确。但是在求解管路图像亚像素边缘时都不能达到最好的效果，仍然存在提取错误或者提取失败的情况。

但是已有方法在提取管路边缘时容易受到管路光照条件的干扰，导致漏提取、多提取、错提取甚至无法提取的情况。为了提取准确、完整的亚像素边缘，本文提出一种针对管路图像提取亚像素精度边缘的方法，先利用频域滤波滤除图像中低频噪声并凸显边缘等高频信号，然后利用聚类分析准确分割管路图像区域，并根据分割结果准确定位求解亚像素边缘的区域，最后利用求解区域内灰度值拟合灰度连续变化曲面，寻找曲面梯度最大点作为边缘，从而获得亚像素精度边缘。本文最大的优势是利用图像处理等一系列手段自动筛选了管路边缘点，排除了由于噪声等问题引起的错误提取点，为重建管路三维模型提供了基础。经过实验验证，本文方法提取管路亚像素边缘准确率达到99%以上，且达到0.04个像素单位的精度，也可用于管路三维重建。

2 处理方法

本文提出的针对管路图像亚像素精度边缘提取方法流程如图2所示，首先频域滤波与区域粗分割获得滤除低频噪声的管路图像，并将管路图像区域连同周围部分背景分割出来，然后通过聚类分析对粗分割管路图像区域进行细致分割，获取精确的管路区域，再利用图像形态学提取像素精度边缘，用于亚像素精度边缘的初值计算和最后的亚像素精度边缘提取计算。

2.1 频域滤波和区域粗分割

利用多目视觉系统拍摄的管路图像具有多种噪声，这些噪声中存在很多低频噪声。本文采用频域高通滤波器，滤除了图像中不同的低频噪声信号，同时达到了锐化管路边缘的目的，其过程如图3所示。

首先，利用傅里叶变换将图像从空域转换为频域，如式(1)所示，

(1)

其中u和v表示频域空间中点的坐标，F(u,v)表示了对应点在频域空间中的值;x和y表示图像中的点坐标，f(x,y)表示了图像中对应像素点的灰度值。

然后利用高通滤波器滤除图像中的低频噪声。滤波后，还需要通过傅里叶逆变换将频域图像转换回空域图像，如式(2)所示。

(2)

式中：f表示还原后的图像，F表示频域图像，M和N分别表示原图像纵横像素数。

对滤波后的图像进行简单的阈值分割，这样的分割结果容易包含非管路区域，或者在分割的管路区域内出现孔洞。因此，利用图像形态学膨胀管路区域，将孔洞或者非管路区域包含在膨胀区域内，从而获取包含管路区域的粗分割结果，如图3中所示。

2.2 聚类分析分割图像

对管路的粗分割区域可以进行细致分割。如图4所示，根据图像中像素灰度特征，可以将粗分割的管路区域分为背景区域、边缘区域、反光区域和背光区域。各区域之间的界线距离较近，反光区域和背光区域复杂交错，采用阈值分割方法容易产生错误提取，导致结果不够准确。为了准确提取管路区域，本文采用聚类分析细致分割图像，提取完整准确的管路区域，从而获取用于亚像素边缘计算的初值。

利用聚类分析可以自动将管路粗分割结果分为如图4所示的4个区域，本文采用模糊C均值聚类(Fuzzy C-Means, FCM)算法计算。对管路粗分割结果中的像素xi进行初始化，任取点作为聚类中心，记为cj(j=1,2,3,4)。每个xi可计算系数uij，代表评定xi属于第j类的测量值，如式(3)表示。

(3)

式中：m表示大于1的实数作为计算的因子。聚类中心cj随着uij变化而变化，如式(4)所示。

(4)

构造目标函数Jm如式(5)所示，以Jm的极小值为优化目标，通过迭代计算完成聚类计算。迭代计算的终止条件如式(6)所示。

(5)

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

(6)

式(6)中ε∈(0,1)，表示终止条件，当两次迭代结果中U的差值小于ε时可以认为达到迭代终止条件。FCM算法步骤如下：

步骤1初始化U，标记为U(0)，并开始迭代计算。

步骤2计算k+1次迭代计算中的各聚类中心，记为C(k+1)。

步骤3根据中心C(k+1)，计算目标函数极值解，计算U(k+1)。

步骤4比较当前U(k)和上一轮迭代U(k-1)的差值，判断终止还是进入下次迭代。

通过计算，图像被划分为4个区域，如图5所示。

采用FCM算法可以对管路图像进行细致分割，提取了完整的管路区域，管路区域包含了反光区域、背光区域和边缘区域。实际上FCM根据灰度值计算了4个聚类中心，将图像分割为4个区域，其中背景区域和边缘区域是进行亚像素边缘提取的主要计算区域。

2.3 提取像素精度边缘

在提取亚像素边缘前，准确定位像素边缘可以有效降低漏提取的情况，也可以为提取亚像素边缘提供计算初值。对于管路图像，提取像素边缘需要满足如下两个需求：①单像素边缘；②分割区域的外层边缘。经典的像素边缘提取算法主要是构建微分算子，与图像卷积运算，设置梯度阈值求取极大阈值检测边缘。此外，还有根据数学形态学，移除区域内部像素提取边缘的算法。通过验证，数学形态学方法更适合本文应用场合。

2.3.1 微分算子提取边缘

检测图像边缘就是检测灰度变化剧烈的部分，通常采用图像灰度梯度变化检测边缘。对于图像中各像素点的梯度如式(7)所示：

(7)

梯度幅值变化为：

‖

(8)

边缘切线方向斜率可以表示为：

(9)

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

根据式(7)构造微分算子检测边缘，常见的微分算子有Prewitt算子、Roberts算子和Sobel算子[24]，如图6所示。Canny算子[24]也是一种检测边缘的算法，相较于其他算子，Canny算子抗噪性突出，利用高斯核函数构造算子滤除噪声，计算梯度极值区域后进行非极大抑制处理，降低非极值当做边缘提取的情况，再利用之后阈值处理，利用双阈值求取边缘并细化获得最终值。Canny算子提取边缘考虑了更多因素，应用更加广泛。LoG算子提取边缘也考虑抗噪问题，通过高斯核函数滤除了大部分噪声。

2.3.2 数学形态学提取边缘

利用形态学将目标区域像素分为边缘像素和内部像素，移除内部像素保留边缘像素提取边缘。通常将目标区域内的像素置1，区域外像素置0，提取区域内3×3区域的像素模板，当中心像素和与其相邻的4个像素都为1时说明中心像素为内部像素，将其置0即完成移除，对目标区域内所有像素都进行处理后即完成边缘提取，其主要流程如图7所示。

为了提取符合前文所述的两个标准，利用不同的方法提取像素精度边缘，得到如图8所示的结果。

图中白色区域表示管路区域，浅灰色区域表示超出管路区域提取的边缘，深灰色表示在管路区域内提取的边缘。对比结果，各方法都成功提取了单像素边缘，但是利用微分算子提取的边缘容易超出管路区域，这是因为利用微分算子卷积计算时按照图像从左到右从上到下的顺序，只根据梯度提取像素精度边缘，忽略了边缘所在区域。而根据图像形态学提取的边缘都是在管路区域内，保证了边缘从区域内到区域外的梯度方向。从提取精度而言，出现超出管路区域边缘并不是提取错误，但是为了结合后续操作提取亚像素精度边缘，推荐采用形态学的方法提取边缘，并计算梯度法线方向作为亚像素精度边缘计算初值。

2.4 亚像素边缘提取

为了得到准确的亚像素边缘，需要根据边缘附近局部灰度值，拟合灰度连续变化曲面，根据局部灰度阈值求解亚像素边缘。以像素精度边缘上一点为中心，其n×n区域内灰度变曲面作为亚像素边缘求解区域，灰度值变化如图9所示。这样的变化趋势与高斯概率密度函数相似，因此本文考虑利用高斯概率密度函数作为边缘连续变化曲面，拟合求解亚像素边缘。采用高斯概率密度函数有两个优势：①相较于高斯模糊函数拟合[21]，不需要对管路图像进行高斯模糊处理，避免了因为高斯模糊导致的图像边缘区域丢失的情况。如图10所示，管路图像边缘区域与背光区域边界距离较近，利用高斯模糊处理后容易将两个区域并在一起，这容易导致提取的亚像素边缘不准确，以背光区域的边缘作为管路边缘错误提取。②拟合初值容易计算，利用计算区域内各像素点灰度值统计结果即可计算大部分初值用于拟合寻优计算，提高了计算效率。

根据图8反映的求解区域灰度变化离散值，拟合高斯概率密度函数，求解亚像素精度边缘，其解析式如式(10)所示，

(10)

式中G(x,y)表示图像中点(x,y)对应的灰度值，a、b、μ、σ和k为需要拟合求解的参数。a表示求解区域灰度最小值，b表示求解区域灰度最大值与最小值的差值，μ和σ为高斯概率密度函数参数，k表示边缘点在图像中的切线方向。以Is为求解区域，拟合目标函数如式(11)所示。

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

(11)

式中G(x,y)表示拟合曲面计算值，Gr(x,y)表示实际灰度值。

拟合灰度连续变化曲面是非线性寻优问题，利用Levenberg-Marquardt法可求解目标函数全局最优解[25]。为了获得稳定收敛的求解结果，需要选取合适的初值进行寻优计算。利用本文提出的亚像素边缘计算方法，可以根据计算区域的灰度值统计结果求解可用于寻优计算的初值。选取n×n的计算区域Is，各点灰度值为Gr(x,y)，其中x,y∈[1,n]。边缘切线斜率k初值通过式(9)求解，其余参数可以通过统计求解区域的灰度值设置初值，各参数初值求解如式(12)～式(15)所示：

a=minGr;

(12)

b=maxGr-a;

(13)

(14)

(15)

为了获得准确的拟合结果，需选取大小适中的计算区域。导管边缘区域比较狭长，为了保证精度，推荐选取7×7或5×5的局部区域用于拟合计算。当选取的计算区域过大时，会受到反光、背光等区域的干扰，导致无法获得收敛结果。而拟合的区域过小，会受到噪声影响，导致拟合结果不准确。通过多次实验，选取7×7或5×5的计算区域能获得准确稳定的结果，当管路截面直径较小时(如5～8 mm)，由于管路边缘区域太小，选取5×5的计算区域更合适。

利用灰度连续变化函数，可以求解亚像素精度边缘以及其对应的梯度方向。设置一个计算区域内的灰度阈值即可定位准确的亚像素边缘点，这个阈值可通过聚类分析的结果获得，记为Gt。计算如式(16)所示，cb表示背景区灰度值中心，ce表示边缘区域灰度值中心。ω为权值参数，经过计算和实验一般可设置在0.2～0.5之间。

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

(16)

将阈值Gt代入式(10)即可得到亚像素边缘的切线方程，如式(17)所示。

(17)

(18)

根据梯度方向[-A-B]T，可以求得亚像素边缘点Esub坐标为

(19)

3 管路亚像素边缘提取实验

验证实验采用AVT F146b型相机，Computar 8mm镜头拍摄，图像像素为1280×960，采用背光光源照射管路。为了说明本文方法能在提取管路亚像素边缘的同时自动筛选提取结果，排除错误的边缘点，利用5个不同的管路进行实验，并通过统计结果说明本文方法准确有效。此外，为了说明本文方法精度，设计了精度实验，求解像素标准圆的圆心坐标偏差和标准方形的边长偏差。

本文方法提取的亚像素边缘结果如图12a和图13a所示，能筛选出需要的管路边缘点，而不受到复杂光照环境的影响。为了显示光照对提取边缘的影响，对比了直接使用的Agustin方法提取的边缘结果[11]。从结果中可知，图12中被拍摄管路的边缘不明显，容易提取失败，因此利用Agustin方法提取的边缘不完整；图13中被拍摄管路的边缘与背光部分边缘较近，容易受到干扰，因此使用Agustin方法容易造成提取错误。而本文方法因为设计了能自动筛选干扰点和错误点的方法，因此可以克服光照环境中的噪声点。

通过对比实验可知，本文针对管路图像的亚像素边缘提取方法能消除管路中复杂光照环境带来的噪声干扰，对比Agustin方法不需要额外的筛选手段，即可自动筛选管路边缘亚像素边缘点，并消除反光位置的边界干扰。为了进一步量化实验结果，以提取的像素精度边缘点数作为提取结果的标准值评判提取结果的准确程度。统计提取的亚像素精度边缘总数，再减去提取错误的情况。与Agustin方法对比提取结果，并通过其他手段对Agustin方法提取的结果进行筛选，选取实际管路图像进行实验，结果如表1所示。

表1 管路边缘提取对比

表1中，提取边缘表示算法直接提取的亚像素边缘点数，通过判断筛选出正确的提取结果，再以筛选结果和标准值计算准确率。对比提取结果，本文方法能不受干扰提取准确的管路亚像素边缘点，而且具有更高的准确率。而Agustin的方法容易受到管路图像中噪声影响，导致提取错误。通过对比可知，本文方法相较于Agustin方法在管路亚像素边缘提取时具有更高的准确率，在计算亚像素边缘点前的一系列处理可以保证获得用于重建计算的亚像素边缘点，而不受到管路复杂光照环境的影响。

为了验证本文亚像素边缘提取精度，采用圆形和方向边缘提取算例，提取结果如图14所示。图14a为直径20个像素的圆形，图14b为边长20像素的正方形。通过本文方法提取的亚像素边缘，计算图14a圆心偏差以及图14b方形边长误差，结果如表2所示。

表2 亚像素边缘提取误差对比

圆心误差XY边长误差0.030 40.029 30.046 1

从提取结果和误差对比结果中可以看出，本文方法提取的亚像素边缘精度达到0.04个像素尺寸，可满足于管路重建与测量。

4 结束语

本文利用提取的管路亚像素精度边缘重建了管路三维模型，测量了管路尺寸。为了保证测量结果准确，需要高精度的亚像素边缘信息。而获取亚像素边缘时容易受到噪声干扰，造成边缘缺失甚至提取错误边缘。本文提出了一种针对管路图像，提取亚像素精度边缘的方法，能在对光照环境复杂的管路图像进行处理，获取亚像素精度边缘，用于管路的三维重建和测量。

为了提取精准的管路图像亚像素边缘，本文首先采用高通滤波器滤除图像中低频信号，保留如边缘高频信号。然后，利用聚类分析划分管路区域，提取完整且准确的管路区域。再利用图像形态学方法获取亚像素边缘求解区域。进一步，根据求解区域的局部灰度变化拟合高斯概率密度函数作为灰度连续变化参数模型。最后，根据之前聚类分析结果求解计算区域内阈值，定位亚像素精度边缘点，完成边缘提取。

通过实验验证，本文方法能提取完整、准确和高精度的亚像素边缘，为基于图像的管路模型重建提供了精度基础。今后会根据提取的亚像素精度边缘重建准确的管路三维模型，用于航空航天管路的快速、高精度测量。