图像边缘信息主要集中在高频段，通常说图像锐化或检测边缘，实质就是高频滤波。我们知道微分运算是求信号的变化率，具有加强高频分量的作用。在空域运算中来说，对图像的锐化就是计算微分。由于数字图像的离散信号，微分运算就变成计算差分或梯度。图像处理中有多种边缘检测（梯度）算子，常用的包括普通一阶差分，Robert算子（交叉差分），Sobel算子等等，是基于寻找梯度强度。拉普拉斯算子（二阶差分）是基于过零点检测。通过计算梯度，设置阈值，得到边缘图像。
何为边缘？图象局部区域亮度变化显著的部分，对于灰度图像来说，也就是灰度值有一个明显变化，既从一个灰度值在很小的缓冲区域内急剧变化到另一个灰度相差较大的灰度值。在提高对景物边缘的敏感性的同时，可以抑制噪声的方法才是好的边缘提取方法。
Canny算子求边缘点具体算法步骤如下：
1. 用高斯滤波器平滑图像．
2. 用一阶偏导有限差分计算梯度幅值和方向.
3. 对梯度幅值进行非极大值抑制.
4. 用双阈值算法检测和连接边缘．
第一步：灰度化
把彩色图像变成灰度图像，该部分是按照Canny算法通常处理的图像为灰度图，如果获取的彩色图像，那首先就得进行灰度化。以RGB格式的彩图为例，通常灰度化采用的公式是：
Gray=0.299R+0.587G+0.114B;
第二步：高斯滤波
对图像高斯滤波，图像高斯滤波的实现可以用两个一维高斯核分别两次加权实现，也就是先一维X方向卷积，得到的结果再一维Y方向卷积。当然也可以直接通过一个二维高斯核一次卷积实现。也就是二维卷积模板，由于水平有限，只说二维卷积模板怎么算。
首先，一维高斯函数：
二维高斯函数：
模板中每一个点的高斯系数可以由上面的公式计算，这样得到的是不是最终的模板呢？答案不是，需要归一化，也即是每一个点的系数要除以所有系数之和，这样才是最终的二维高斯模板。
这个里面有个小知识点，要想计算上面的系数，需要知道高斯函数的标准差σ （sigma），还需要知道选3*3还是5*5的模板，也就是模板要多大，实际应用的时候，这两者是有关系的，根据数理统计的知识，高斯分布的特点就是数值分布在（μ—3σ,μ+3σ)中的概率为0.9974，也就是模板的大小其实就是6σ这么大就OK了，但是6σ可能不是奇数，因为我们一定要保证有核心。所以模板窗口的大小一般采用1+2*ceil(3*nSigma) ceil是向上取整函数，例如ceil(0.6)=1。
计算得到模板，那就是直接卷积就OK，卷积的意思就是图像中的点附近的模板大小区域乘以高斯模板区域，得到的结果就是该点卷积后的结果。卷积的核心意义就是获取原始图像中像模板特征的性质。
第三步：计算梯度值和方向
图像的边缘可以指向不同方向，因此经典Canny算法用了四个梯度算子来分别计算水平，垂直和对角线方向的梯度。但是通常都不用四个梯度算子来分别计算四个方向。常用的边缘差分算子（如Rober，Prewitt，Sobel）计算水平和垂直方向的差分Gx和Gy。这样就可以如下计算梯度模和方向：
梯度角度θ范围从弧度-π到π，然后把它近似到四个方向，分别代表水平，垂直和两个对角线方向（0°,45°,90°,135°）。可以以±iπ/8（i=1,3,5,7）分割，落在每个区域的梯度角给一个特定值，代表四个方向之一。
这里选择Sobel算子计算梯度，相对于其他边缘算子，Sobel算子得出来的边缘粗大明亮。
Sobel算子，无阈值
第四步：非极大值抑制
非极大值抑制是进行边缘检测的一个重要步骤，通俗意义上是指寻找像素点局部最大值。沿着梯度方向，比较它前面和后面的梯度值进行了。见下图。
上图中左右图：g1、g2、g3、g4都代表像素点，很明显它们是c的八领域中的4个，左图中c点是我们需要判断的点，蓝色的直线是它的梯度方向，也就是说c要是局部极大值，它的梯度幅值M需要大于直线与g1g2和g2g3的交点，dtmp1和dtmp2处的梯度幅值。但是dtmp1和dtmp2不是整像素，而是亚像素，也就是坐标是浮点的，那怎么求它们的梯度幅值呢？线性插值，例如dtmp1在g1、g2之间，g1、g2的幅值都知道，我们只要知道dtmp1在g1、g2之间的比例，就能得到它的梯度幅值，而比例是可以靠夹角计算出来的，夹角又是梯度的方向。
写个线性插值的公式：设g1的幅值M(g1)，g2的幅值M(g2),则dtmp1可以很得到：
M(dtmp1)=w*M(g2)+(1-w)*M(g1)
其中w=distance(dtmp1,g2)/distance(g1,g2)
distance(g1,g2) 表示两点之间的距离。实际上w是一个比例系数，这个比例系数可以通过梯度方向(幅角的正切和余切)得到。
右边图中的4个直线就是4个不同的情况，情况不同，g1、g2、g3、g4代表c的八领域中的4个坐标会有所差异，但是线性插值的原理都是一致的。
下图是非最大值抑制的结果。可见边缘宽度已经大大减小。但是这个图像中因为没有应用任何阈值，还含有大量小梯度模值的点，也就是图中很暗的地方。下面，阈值要上场了。
第五步：双阈值的选取
一般的边缘检测算法用一个阈值来滤除噪声或颜色变化引起的小的梯度值，而保留大的梯度值。Canny算法应用双阈值，即一个高阈值和一个低阈值来区分边缘像素。如果边缘像素点梯度值大于高阈值，则被认为是强边缘点。如果边缘梯度值小于高阈值，大于低阈值，则标记为弱边缘点。小于低阈值的点则被抑制掉。
第六步：滞后边界跟踪
强边缘点可以认为是真的边缘。弱边缘点则可能是真的边缘，也可能是噪声或颜色变化引起的。为得到精确的结果，后者引起的弱边缘点应该去掉。通常认为真实边缘引起的弱边缘点和强边缘点是连通的，而又噪声引起的弱边缘点则不会。所谓的滞后边界跟踪算法检查一个弱边缘点的8连通领域像素，只要有强边缘点存在，那么这个弱边缘点被认为是真是边缘保留下来。
这个算法搜索所有连通的弱边缘，如果一条连通的弱边缘的任何一个点和强边缘点连通，则保留这条弱边缘，否则抑制这条弱边缘。搜索时可以用广度优先或者深度优先算法，我在这里实现了应该是最容易的深度优先算法。一次连通一条边缘的深度优先算法如下：
- 准备一个栈s，一个队列q，设联通指示变量connected为假。从图像的第一个点开始，进入2。
- 如果这个点是弱边界点并且没有被标记，把它标记，并把它作为第一个元素放入栈s中，同时把它放入记录连通曲线的队列q，进入3。如果这个点不是弱边界或者已经被标记过，到图像的下一个点，重复2。
- 从栈s中取出一个元素，查找它的8像素领域。如果一个领域像素是弱边界并且没有被标记过，把这个领域像素标记，并加入栈s中，同时加入队列q。同时查找领域对应的强边界图，如果有一个像素是强边界，表示这条弱边界曲线和强边界联通，设置connected为真。重复3直到栈中没有元素了。如果connected为假，则依次从队列q中取出每个元素，清空标记。如果connected为真，保留标记。
- 清空队列q，设置connected为假，移动到图像的下一个点，到2。
第七步：结果查看
下面是对Lena图计算Canny边缘检测的梯度模图和二值化图，高斯半径2，高阈值100，低阈值50。
作为对比，下面是用一阶差分和Sobel算子对原图计算的结果，阀值100。由于一阶差分的梯度值相对较小，我对一阶差分的梯度值放大了一定倍数，使得它和Sobel的梯度值保持同样的水平。
很明显，Canny边缘检测的效果是很显著的。相比普通的梯度算法大大抑制了噪声引起的伪边缘，而且是边缘细化，易于后续处理。对于对比度较低的图像，通过调节参数，Canny算法也能有很好的效果。
public Bitmap CANNY(Bitmap oldBitmap, int Width, int Height) //定义canny函数
{
Bitmap newBitmap = new Bitmap(Width, Height);
cannyClass2 cannyclass3 = new cannyClass2(); //生成类的实例
float gradientMax = 0;
float[,] gradient = new float[Width, Height];
byte[,] degree = new byte[Width, Height];
int[,] srcBytes = new int[Width, Height];
int highThreshould = 24; //高阀值
int lowThreshould = 10; //低阀值
int ret = 0;
for (int j = 0; j < Height; j++)
{
for (int i = 0; i < Width; i++)
{
srcBytes[i, j] = oldBitmap.GetPixel(i, j).R;
}
}
cannyclass3.GetGradientDegree(srcBytes, ref gradient, ref degree, ref gradientMax, Width, Height); //调用类中 求梯度函数
cannyclass3.NonMaxMini(gradient, ref srcBytes, gradientMax, Width, Height, degree); //非极大值抑制
cannyclass3.TwoThreshouldJudge(highThreshould, lowThreshould, ref srcBytes, Width, Height); //双阈值边缘判断
for (int j = 0; j < Height; j++)
{
for (int i = 0; i < Width; i++)
{
ret = srcBytes[i, j];
newBitmap.SetPixel(i, j, Color.FromArgb(ret, ret, ret));
}
}
return newBitmap;
}
高斯滤波
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Drawing;
namespace WindowsFormsApplication1
{
class gsClass1
{
public int BlurRadius = 1; //高斯函数宽度
public List<double> BlurArray = new List<double>();
public int MaxWidth;
public int MaxHeight;
private double GetWeighing(int x, int y) //定义高斯公式
{
double q = (BlurRadius * 2 + 1) / 2;
return 1 / (2 * Math.PI * Math.Pow(q, 2)) * Math.Exp(-(x * x + y * y) / (2 * q * q));
}
public Color GetBlurColor(Bitmap oldBitmap, int x, int y)
{
double r = 0, g = 0, b = 0;
int index = 0;
for (var t = y - this.BlurRadius; t <= y + this.BlurRadius; t++)
{
for (var l = x - this.BlurRadius; l <= x + this.BlurRadius; l++)
{
var color = GetDefautColor(oldBitmap, l, t);
var weighValue = BlurArray[index];
r += color.R * weighValue;
g += color.G * weighValue;
b += color.B * weighValue;
index++;
}
}
return Color.FromArgb((byte)r, (byte)g, (byte)b);
}
public Color GetDefautColor(Bitmap oldBitmap, int x, int y)
{
if (x < 0 && y < 0)
return oldBitmap.GetPixel(0, 0);
else if (x < 0)
return oldBitmap.GetPixel(0, Math.Min(MaxHeight, y));
else if (y < 0)
return oldBitmap.GetPixel(Math.Min(MaxWidth, x), 0);
else
return oldBitmap.GetPixel(Math.Min(MaxWidth, x), Math.Min(MaxHeight, y));
}
public void SetBlurArray()
{
int blur = BlurRadius;
double sum = 0;
for (var y = blur; y >= blur * -1; y--)
{
for (var x = blur * -1; x <= blur; x++)
{
var d = GetWeighing(x, y);
BlurArray.Add(d);
sum += d;
}
}
for (var i = 0; i < BlurArray.Count; i++)
BlurArray[i] = BlurArray[i] / sum;
//sum = 0;
//foreach (var item in this.BlurArray)
// sum += item;
}
}
}
CANNY
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Drawing;
namespace WindowsFormsApplication1
{
class cannyClass2
{
public void NonMaxMini(float[,] gradient, ref int[,] srcBytes, float GradientMax, int x, int y, byte[,] degree) //非极大值抑制
{
float leftPixel = 0, rightPixel = 0;
for (int j = 1; j < y - 1; j++)
{
for (int i = 1; i < x - 1; i++)
{
switch (degree[i, j])
{
case 0:
leftPixel = gradient[i - 1, j];
rightPixel = gradient[i + 1, j];
break;
case 45:
leftPixel = gradient[i - 1, j + 1];
rightPixel = gradient[i + 1, j - 1];
break;
case 90:
leftPixel = gradient[i, j + 1];
rightPixel = gradient[i, j - 1];
break;
case 135:
leftPixel = gradient[i + 1, j + 1];
rightPixel = gradient[i - 1, j - 1];
break;
default:
break;
}
if ((gradient[i, j] < leftPixel) || (gradient[i, j] < rightPixel))
{
srcBytes[i, j] = 0;
}
else
{
srcBytes[i, j] = (int)(255 * gradient[i, j] / GradientMax);
}
}
}
}
//双阈值边缘判断
public void TwoThreshouldJudge(int highThreshold, int lowThreshould, ref int[,] srcBytes, int x, int y)
{
for (int j = 1; j < y - 1; j++)
{
for (int i = 1; i < x - 1; i++)
{
if (srcBytes[i, j] > highThreshold)
{
srcBytes[i, j] = 255;
}
else if (srcBytes[i, j] < lowThreshould)
{
srcBytes[i, j] = 0;
}
else
{
if (srcBytes[i - 1, j - 1] < highThreshold && srcBytes[i, j - 1] < highThreshold && srcBytes[i + 1, j - 1] < highThreshold && srcBytes[i - 1, j] < highThreshold
&& srcBytes[i + 1, j] < highThreshold && srcBytes[i - 1, j + 1] < highThreshold && srcBytes[i, j + 1] < highThreshold && srcBytes[i + 1, j + 1] < highThreshold)
{
srcBytes[i, j] = 0;
}
else
srcBytes[i, j] = 255;
}
}
}
}
public void GetGradientDegree(int[,] srcBytes, ref float[,] gradient, ref byte[,] degree, ref float GradientMax, int x, int y) //梯度公式
{
gradient = new float[x, y];
degree = new byte[x, y];
int gx, gy;
int temp;
double div;
for (int j = 1; j < y - 1; j++)
{
for (int i = 1; i < x - 1; i++)
{
gx = srcBytes[i + 1, j - 1] + 2 * srcBytes[i + 1, j] + srcBytes[i + 1, j + 1] - srcBytes[i - 1, j - 1] - 2 * srcBytes[i - 1, j] - srcBytes[i - 1, j + 1];
gy = srcBytes[i - 1, j - 1] + 2 * srcBytes[i, j - 1] + srcBytes[i + 1, j - 1] - srcBytes[i - 1, j + 1] - 2 * srcBytes[i, j + 1] - srcBytes[i + 1, j + 1];
gradient[i, j] = (float)Math.Sqrt((double)(gx * gx + gy * gy)); //梯度公式
if (GradientMax < gradient[i, j])
{
GradientMax = gradient[i, j];
}
if (gx == 0)
{
temp = (gy == 0) ? 0 : 90;
}
else
{
div = (double)gy / (double)gx;
if (div < 0)
{
temp = (int)(180 - Math.Atan(-div) * 180 / Math.PI);
}
else
{
temp = (int)(Math.Atan(div) * 180 / Math.PI);
}
if (temp < 22.5)
{
temp = 0;
}
else if (temp < 67.5)
{
temp = 45;
}
else if (temp < 112.5)
{
temp = 90;
}
else if (temp < 157.5)
{
temp = 135;
}
else
temp = 0;
}
degree[i, j] = (byte)temp;
}
}
}
}
}
