基于sobel和canny算法的水珠边缘检测
李建国1，李幸汶1，卜祥洲1，陈涛2
（1.河南宏博测控技术有限公司；2.河南工业大学）
摘要：水珠识别分析方法是利用数字图像处理及模式识别技术实现，本文简要说明了水珠通过边缘检测实现模式识别和计算的原理。同时，着重讨论了水珠识别过程中的噪声识别算法，详尽地阐述了识别过程中的阈值选取方法。此应用表明了算法的实时性好，可靠性高，能很好地满足工业生产和检测要求，大大提高了水珠识别的自动化与智能化。
关键词：水珠识别 边缘检测 边缘轮廓
引言
图像识别是计算机视觉领域的一个非常活跃的课题。近几年来人们对图像识别做了大量的研究，提出了很多种算法和方法，但电力行业应用的相对较少。传统的图像识别方法已不能满足高精度，高效率的要求。电力输配电企业对绝缘子的憎水性要求较高，因为其憎水性直接决定着绝缘子绝缘性能的好坏，并可以减少高压电击穿和污闪故障的发生。因绝缘子的憎水性变差而导致泄漏电流的增加，大大提高了输配电线路的安全隐患。因此，迫切需要提高绝缘子憎水性能安全检测的手段，实现检测的自动化和智能化。目前，检测系统都是拆卸绝缘子并通过实验室法进行测试，不能实现在线检测。对于在线喷水并通过CCD取像得到的图像，核心功能是水珠识别算法，算法的不同就决定了检测效果有很大的差别。那么，设计一个既快速又稳定的水珠识别算法，对水珠识别来说具有十分重要的意义。本文就是在如何准确提取水珠面积（或轮廓）的思想下，提出的一种新的识别算法。且通过人机交互接口，还可以满足其它行业（如建材或医疗）不同精度和质量级别需求。
工作原理
其基本工作原理：喷水装置在喷雾过程中，通过取像 CCD（Charge coupled Device）在线拍摄图像数据，图像数据经图像采集卡收集并传入计算机，由系统水珠识别算法设置相应的阈值，检测水珠边缘界线并计算水珠面积和所占用面积比，同时通过水珠的边缘轮廓计算水珠的形状因子。形状因子决定了水珠是否有粘连或流动，而面积反应了水珠的大小。计算这些元素，首先需要得到水珠的外形或形状轮廓，即：水珠边缘的检测。
Sobel（索贝尔）算子：
索贝尔算子（Sobel operator）是图像处理中的常用算子之一，主要用作边缘检测。在
技术上，它是一离散性差分算子，用来运算图像亮度函数的梯度之近似值。在图像的任何一点使用此算子，将会产生对应的梯度矢量或是其法矢量。对梯度幅值进行非极大值抑制，若（x，y）像素点的梯度幅值M（x，y）大于或等于沿梯度方向上，即：平行于H（x，y）方向的两个相邻像素点的梯度幅值，则判定该点为可能的边缘点。[5]p3
该算子包含两组3x3的矩阵，分别为横向及纵向运算模板，如果用Sobel算子检测图像的边缘，可以先分别用水平算子和垂直算子对图像进行平面卷积[2]p1，即可分别得出横向及纵向的亮度差分近似值。如果以A代表原始图像，Gx及Gy分别代表经横向及纵向边缘检测的图像，其公式如下:
??10?1???1?2?1?????00??A Gx=??20?2??A Gy=?0??10?1???1?2?1?????图像的每一个像素的横向及纵向梯度近似值可用以下的公式结合，来计算梯度的大小： Gx = (-1)*f(x-1,y-1)+0*f(x,y-1)+1*f(x+1,y-1)+(-2)*f(x-1,y)+0*f(x,y)+2*f(x+1,y)+
(-1)*f(x-1,y+1)+0*f(x,y+1)+1*f(x+1,y+1)
= [f(x+1,y-1)+2*f(x+1,y)+f(x+1,y+1)]-[f(x-1,y-1)+2*f(x-1,y)+f(x-1,y+1)]Gy = 1*f(x-1,y-1)+2*f(x,y-1)+1*f(x+1,y-1)+0*f(x-1,y)+0*f(x,y)+0*f(x+1,y)+ (-1)*f(x-1,y+1)+(-2)*f(x,y+1)+(-1)*f(x+1,y+1)
= [f(x-1,y-1)+2*f(x,y-1)+f(x+1,y-1)]-[f(x-1,y+1)+2*f(x,y+1)+f(x+1,y+1)] 其中f(a,b),表示图像(a,b)点的灰度值。图像的每一个像素的横向及纵向灰度值通过以下公式结合,来计算该点灰度的大小：
22G?Gx?Gy[3]p2
通常,为了提高效率,使用不开方的近似值：
G?Gx?Gy
如果梯度G大于某一阈值，则认为该点（x，y）为边缘点。 可用以下公式计算梯度方向： ?Gy???arctan?G?x??? ? 上例中，如果角度Θ等于零，即代表图像在该处拥有纵向边缘，左方较右方暗。
在边沿检测中，常用的一种模板是Sobel 算子。Sobel 算子有两个，一个用来检测水平边沿，另一个用来检测垂直边沿。与“和”相比，Sobel算子对于象素的位置的影响做了加权，因此效果更好。
Sobel算子另一种形式是各向同性Sobel(Isotropic Sobel)算子，也由一个水平和垂直边沿检测算子组成。各向同性Sobel算子和普通Sobel算子相比，它的位置加权系数更为准确，在检测不同方向的边沿时梯度的幅度一致。由于水珠图像的特殊性，处理该类型图像轮廓时，并不需要对梯度方向进行运算，所以程序并没有给出各向同性Sobel算子的处理方法。
由于Sobel算子是滤波算子的形式，用于提取边缘，可以利用快速卷积函数。因为其简
单有效，因此应用广泛。美中不足的是Sobel算子并没有将图像的主体与背景严格地区分开来，Sobel算子没有基于图像灰度进行处理，且没有严格地模拟人的视觉生理特征，所以提取的图像轮廓有时并不能令人满意。处理图像的时候关注的是图像与背景不同的部分，并没有把关注点与图像背景区分开来。基于此给出了下面阈值化轮廓提取算法，该算法已在数学上证明当像素点满足正态分布时所求解是最优化的。 Canny算子：
拉普拉斯算子是2阶微分算子，相当于求取2次微分，精度比较高，但对噪声过于敏感是它的重大缺点，所以这种算子并不是特别常用。索贝尔（Sobel）算子是最常用的一种一阶算子，方法简单、效果好，但提取出的边缘比较粗，要进行细化处理。坎尼（Canny）算子在使用时，要提供给一些参数用于控制算法的性能。对于不同的图像或不同的边缘提取要求，应该提供不同的参数，以达到最佳效果。它也有模板运算方法，模板的大小也比较大，与提供的参数有关。标准的大小差不多是17*17，可根据算子的可分离性用快速算法，使用标准算子非常慢。坎尼（Canny）算子的2值化也很有特色，具有一定的智能性。
RGB图像边缘计算在图像中的方向相同，只有两个边缘具有相同的方向。如果问题只是检测边缘，单独分量方法通常可产生可接受的结果。然而，如果把准确度做为重点，明显需要一个适用于向量梯度的新定义。[1]p281
Canny边缘检测算法一般分为四步：
一、用一个高斯滤波器平滑输入图象；[1]p463~466
二、用一阶偏导的有限差分来计算梯度幅值和方向（角度）图像； 三、对梯度幅值图像进行非极大值抑制； 四、用双阈值算法检测和连接边缘。
坎尼方法基于三个基本目标：
1．低错误率。所有边缘都应被找到，并且应该没有伪响应。也就是检测以的边缘必须
可能是真实的边缘。
2．边缘点应被很好地定位。已定位边缘必须尽可能接近真实边缘。也就是由检测标记为边缘的点和真实边缘的中心之间的距离应该最小。
3．单一的边缘点响应。对于真实的边缘点，检测器仅应返回一个点。也就是真实边缘周围的局部最大数应该是最小的。这意味着在仅存一个单一边缘点的位置，检测器不应指出多个边缘像素。
坎尼的工作本质是从数学上表达了上面的三个准则，并试图找到这些表达式的最佳解。通常，寻找一个满足这一目标的、严格形式的解是很困难的（或不可能的）。然而，对由加性高斯白噪声污染的一维台阶边缘使用数字最佳化会得到这样一个结论，即参最佳台阶边缘检测器的一个较好近似是高斯一阶导数：
dedx?x22?2??x?2e?x22?2
使用一个环形二维高斯函数平滑图像，计算结果的梯度，然后使用幅度和方向来估计每一点处的边缘强度和方向。令
f(x,y)表示输入图像， G(x,y)表示高斯函数：
x2?y2G(x,y)?e2?2
利用G(x,y)和f(x,y)的卷积形成一幅平滑的图像fx(x,y)：
fs(x,y)?G(x,y)?f(x,y)
这一操作源自计算梯度和方向：
?gy?22M(x,y)?gx?gy和?(x,y)?arctan ?? 其中，gx??fs?x，gy??fs?y。
?gx?任何滤波器模板对都可以用来得到gx和gy，使用一个大小为n?n的高斯模板来实现，
M(x,y)和?(x,y)与计算出它们的图像尺寸相同的阵列。
Canny算法提取的图像边缘基本上可以满足视觉的要求，但达不到单像素级别，尤其是在边缘的拐角点处，多对一响应现象严重。为获得单像素级精细边缘，我们引入数学形态学算法实现对边缘细化处理，本文采用文献
[4]的算法。
算法的设计
基本工作原理：CCD（Charge coupled Device）拍摄的水珠图像数据经采集卡收集并传入
计算机，由系统降噪算法识别杂质，经过滤波后去除无关信息。经过边缘检测和连接边缘形成完整的水珠轮廓，分别计算闭合区域面积和周长，得到最大水珠的面积比。 2.1 建立憎水性对照表
根据不同的水珠面积和比例，参照相关的周长与形状因子的大小，可以确定介质的憎水性能。所以，在进行测试和计算之前，需要一个对照表。 2.2边缘检测实验程序：
这是一个图像对象提取方案,分别在水平方向和垂直方向检测边缘。
// Sobel算子
/* m_nValues = {-1, 0, 1
-2, 0, 2 -1, 0, 1} */
// 检测垂直平边沿
C3By3Filter Filter1;
Filter1.m_nValues[0][0] = -1; Filter1.m_nValues[0][1] = 0; Filter1.m_nValues[0][2] = 1; Filter1.m_nValues[1][0] = -2; Filter1.m_nValues[1][1] = 0; Filter1.m_nValues[1][2] = 2; Filter1.m_nValues[2][0] = -1; Filter1.m_nValues[2][1] = 0; Filter1.m_nValues[2][2] = 1;
// 检测水平边沿
/* m_nValues = {-1,-2,-1
0, 0, 0 1, 2, 1} */
C3By3Filter Filter2;
Filter2.m_nValues[0][0] = -1; Filter2.m_nValues[0][1] = -2; Filter2.m_nValues[0][2] = -1;