2.计算梯度 这一步的主要目的是对边缘进行增强,以便进一步进行边缘提 取。这里用常用的Sobel算子来计算梯度,由于后期涉及梯度方向上的 非最大值抑制,需同时计算出梯度的模值和方向。C++代码如下:
/*接着求Sobel边缘*/
int Templet_Y[3][3]= {
-1,-2,-1,
0,0,0,
1,2,1
};
int Templet_X[3][3]= {
-1,0,1,
-2,0,2,-1,0,1
};
double Sobel_Result_X,Sobel_Result_Y,Sobel_Result;
double* p_SobelResult_X = new
double[m_dwWidth*m_dwHeight];
/*X方向梯度幅值*/
double* p_SobelResult_Y = new
double[m_dwWidth*m_dwHeight];
/*Y方向梯度幅值*/
double* p_SobelResult = new
double[m_dwWidth*m_dwHeight];
/*梯度幅值*/
double* p_theta = new double[m_dwWidth*m_dwHeight];
/*梯度方向*/
BYTE *p_Sobel_Disp = new BYTE[m_dwWidth*m_dwHeight];
/*梯度显示缓存*/
for (i = 1; i<m_dwHeight - 1; i++)
for (j = 1; j<m_dwWidth - 1; j++)
{
Sobel_Result_X = 0;
Sobel_Result_Y = 0;
Sobel_Result = 0;
for (m = -1; m <= 1; m++)
for (k = -1; k <= 1; k++)
{
Sobel_Result_X += Templet_X[m+1]
[k+1]*p_GausResult[(i + m)*m_dwWidth + j + k];
Sobel_Result_Y += Templet_Y[m+1]
[k+1]*p_GausResult[(i + m)*m_
dwWidth + j + k];
}
Sobel_Result = ((double)Sobel_Result_X)*
((double)Sobel_Result_X);
Sobel_Result += ((double)Sobel_Result_Y)*
((double)Sobel_Result_Y);
Sobel_Result = sqrt(Sobel_Result);
p_Sobel_Disp[i*m_dwWidth + j] = (Sobel_Result >255)
? 255 :Sobel_Result;
p_SobelResult[i*m_dwWidth + j] = Sobel_Result; //梯
度模值
p_SobelResult_X[i*m_dwWidth + j] = Sobel_Result_X;
//X方向模值
p_SobelResult_Y[i*m_dwWidth + j] = Sobel_Result_Y;
//Y方向模值
p_theta[i*m_dwWidth + j] =
atan2(Sobel_Result_Y,Sobel_
Result_X)* 180/3.1415926; //Y方向,以°为单位
if (p_theta[i*m_dwWidth + j]< 0) //将角度转换到0~
360范围
p_theta[i*m_dwWidth + j]+= 360;
}
Sobel计算的结果如图10-5所示。
3.非最大值抑制(Non-Maximum Suppression)
图像梯度幅值矩阵中的元素值越大,说明图像中该点的梯度值越 大,但这不能说明该点就是边缘(这仅仅是属于图像增强的过程)。 在Canny算法中,非极大值抑制是进行边缘检测的重要步骤,通俗意义 上是指寻找像素点局部最大值,将非极大值点所对应的灰度值置为0, 从而去除潜在的伪边缘。 我们接下来介绍极大值抑制的原理。对当前像素的梯度值(也就 是Sobel计算结果)进行3×3开窗,如图10-6所示。 图10-6中a0 ~a8 为当前窗口的3×3邻域像素点,a4 为当前窗口 中心像素。向量m2 m1 为当前像素点a4 的梯度方向。 所谓极大值抑制就是确定a4 是否是在邻域内最大,图中斜线方向 为a4 点的梯度方向。因此,可以确定其局部的最大值分布在这条线 上,即除了a4 点,梯度方向的交点m1 和m2 这两个点的值也可能会是 局部最大值。 然而m1 和m2 并不是刚好分布在整数邻域上,这时候需要估计m1 和m2 的值。最好的办法当然是对其进行线性插值。在图10-6中,我们用a2 和a5 的值来对m1 进行插值,用a3 和a6 的值来对m1 进行插值。设a4 a5 的距离为x,a5 m1 的距离为y(x与y 均大于0),插值函数为f,则插值结果为
我们把整个坐标轴分为8个象限,如图10-8中的虚线所示,即图中 的象限位置为第一象限,逆时针依次为第2象限第3象限,一直到第8象 限。 可以预见的是在不同的象限,由于梯度方向所落到的邻域位置不 同,插值公式是不一样的。图10-6中展示了其中的两个象限的插值公 式,不妨列出另外三种情况,如图10-7~图10-9所示。
/*下一步进行非最大值抑制*/
BYTE* N = new BYTE[m_dwWidth*m_dwHeight]; //非极大值抑
制结果
int g1 = 0,g2 = 0,g3 = 0,g4 = 0;//用于进行插值,得到亚
像素点坐标值
double dTmp1 = 0.0,dTmp2 = 0.0; //保存两个亚像素点插值
得到的灰度数据
double dWeight = 0.0; //插值的权重
for (i = 1; i< m_dwHeight - 1; i++)
for (j = 1; j< m_dwWidth - 1; j++)
{
int nPointIdx = i*m_dwWidth + j;//当前点在图像数组
中的索引值if (p_SobelResult[nPointIdx]== 0)N[nPointIdx]= 0;
//如果当前梯度幅值为0,则不是局部最大对该点赋为0
else
{
if (((p_theta[nPointIdx]>= 90) &&
(p_theta[nPointIdx]<135)) ||
((p_theta[nPointIdx]>=270)&&
(p_theta[nPointIdx]<315)))
{
//根据斜率和四个中间值进行插值求解
g1 = p_SobelResult[nPointIdx - m_dwWidth -
1];
g2 = p_SobelResult[nPointIdx - m_dwWidth];
g3 = p_SobelResult[nPointIdx + m_dwWidth];
g4 = p_SobelResult[nPointIdx + m_dwWidth +
1];
dWeight = fabs(p_SobelResult_X[nPointIdx])
/ fabs
(p_SobelResult_Y[nPointIdx]); //反正切
dTmp1 = g1*dWeight + g2*(1 - dWeight);
dTmp2 = g4*dWeight + g3*(1 - dWeight);
}
else if(((p_theta[nPointIdx]>=135)&&
(p_theta[nPointIdx]<180))||
((p_theta[nPointIdx]>= 315) &&
(p_theta[nPointIdx]<360)))
{g1 = p_SobelResult[nPointIdx - m_dwWidth -
1];
g2 = p_SobelResult[nPointIdx - 1];
g3 = p_SobelResult[nPointIdx + 1];
g4 = p_SobelResult[nPointIdx + m_dwWidth +
1];
dWeight = fabs(p_SobelResult_Y[nPointIdx])
/ fabs
(p_SobelResult_X[nPointIdx]); //正切
dTmp1 = g1*dWeight + g2*(1 - dWeight);
dTmp2 = g4*dWeight + g3*(1 - dWeight);
}
else if (((p_theta[nPointIdx]>= 45) &&
(p_theta[nPointIdx]<90)) ||
((p_theta[nPointIdx]>= 225) &&
(p_theta[nPointIdx]<270)))
{
g1 = p_SobelResult[nPointIdx - m_dwWidth];
g2 = p_SobelResult[nPointIdx - m_dwWidth +
1];
g3 = p_SobelResult[nPointIdx + m_dwWidth];
g4 = p_SobelResult[nPointIdx + m_dwWidth -
1];
dWeight = fabs(p_SobelResult_X[nPointIdx])
/ fabs
(p_SobelResult_Y[nPointIdx]); //反正切
dTmp1 = g2*dWeight + g1*(1 - dWeight);dTmp2 = g4*dWeight + g3*(1 - dWeight);
}
else if (((p_theta[nPointIdx]>= 0) &&
(p_theta[nPointIdx]<45)) ||
((p_theta[nPointIdx]>= 180) &&
(p_theta[nPointIdx]<225)))
{
g1 = p_SobelResult[nPointIdx - m_dwWidth +
1];
g2 = p_SobelResult[nPointIdx + 1];
g3 = p_SobelResult[nPointIdx + m_dwWidth -
1];
g4 = p_SobelResult[nPointIdx - 1];
dWeight = fabs(p_SobelResult_Y[nPointIdx])
/ fabs
(p_SobelResult_X[nPointIdx]); //正切
dTmp1 = g1*dWeight + g2*(1 - dWeight);
dTmp2 = g3*dWeight + g4*(1 - dWeight);
}
}
// 进 行 局 部 最 大 值 判 断 , 并 写 入 检 测 结
果
if((p_SobelResult[nPointIdx]>=dTmp1)&&
(p_SobelResult[nPointIdx]>=dTmp2))
N[nPointIdx]= 255;
else
N[nPointIdx]= 0;}
图10-10是进行非极大值抑制后的结果。
完成非极大值抑制后,会得到一个二值图像,非边缘的点灰度值 均为0,可能为边缘的局部灰度极大值点可设置其灰度为255。但是这 样一个检测结果还是包含了很多由噪声及假边缘。 4.滞后阈值分割及边缘连接 算子的最后一步是滞后阈值分割与边缘连接。采用阈值分割的主 要目的是消除假边缘。Canny算法中减少假边缘数量的方法是采用双阈 值法。选择两个阈值,根据高阈值得到一个边缘图像,这样一个图像 含有很少的假边缘。但是由于阈值较高,产生的图像边缘可能不闭 合,为解决这样一个问题采用了另外一个低阈值。 具体算法如下: (1)若梯度值大于高阈值,则认定为边缘。 (2)若梯度值小于低阈值,则认定为非边缘。 (3)若梯度值在两个阈值之间,在该像素周围(3×3的邻域)寻 找是否有边缘点(梯度值大于高阈值)。若有,则认定该点为边缘 点。
算法示例如下:
double dThrHigh = 20; //高阈值
double dThrLow = 10; //低阈值
BYTE * m_pTemp = new BYTE[m_dwHeight*m_dwWidth];
memset(m_pTemp,255,m_dwHeight*m_dwWidth);
double sum,temp1,temp2 = 0;
//最后的计算结果
BYTE* p_CannyResult = new BYTE[m_dwWidth*m_dwHeight];//非
极大值抑制结果
for (i = 1; i<m_dwHeight-1; i++)
for (j = 1; j<m_dwWidth-1; j++)
{
if (N[i*m_dwWidth + j]== 255)
{
//小于低阈值的置零
if ((p_SobelResult[i*m_dwWidth+j]<
=dThrLow))p_CannyResult[i*m_
dwWidth+j] = 0;
//大于高阈值的置1
else
if((p_SobelResult[i*m_dwWidth+j]>=dThrHigh))p_CannyResult
[i*m_dwWidth+j]=255;
//中间值进行邻域判决
else
{
int flag = 0;
for (m = -1; m <= 1; m++)for (k = -1; k <= 1; k++)
{
if (p_SobelResult[(i + m)*m_dwWidth + j + k]>=
dThrHigh)
{
flag++;
}
}
if(flag>=1)p_CannyResult[i*m_dwWidth + j]= 255;
else p_CannyResult[i*m_dwWidth + j]= 0;
}
}
else p_CannyResult[i*m_dwWidth + j]= 0;
}
可见,阈值处理消除了大量的假边缘,得到了比较精细的边缘图 像。我们可以把结果与Soble运算的结果进行比较:Canny算子的运算结果中的边缘仅有一个像素单元的宽度,大大减小了边缘的范围。 不足的是,此算法需要两个手动输入的阈值进行分割。当然也有 一些自动阈值计算的算法,例如前面所介绍的大津法阈值(OTSU)计 算,图10-12是以OTSU分割结果作为低阈值并且对高阈值取2倍低阈值 的Canny运算结果。