Exercise7--边缘检测

发表于 | Visitors
Words count in article 3.6k | Reading time 17

源码在Github

通常来说，边缘所在的位置，都是与其四周的像素颜色区别较大的地点的。

因此，边缘检测是基于灰度突变来分割图像的常用方法。

Prewitt

prewitt算子是一种比较简单的，基于一阶微分算子的边缘检测。

其原理即是，利用像素点的上下，左右邻点之间的灰度值之差，从而计算出他们的一阶导数，当导数达到极值点，即为边缘。

最终计算出来的值可以是

也可以是

在我们计算出他们之间的一阶导数值之后，只需要指定一个阈值，既可以简单的分割出一个图片的边缘了。

这种方法的有点在于简单。而缺点也很明显，那就是对于噪点的抵挡能力比较弱，很容易就会被噪点所干扰，因为许多噪点的灰度值也是比较大的。

代码的简单实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
ImageUtil::IMGDATA prewitt(ImageUtil::IMGDATA data, const int threadhold)
{
	BYTE *img = new BYTE[data.width * data.height];
	memset(img, 0, data.width * data.height);
	for (int i = 1; i < data.height - 1; i++)
	{
		for (int j = 1; j < data.width - 1; j++)
		{
			const int dx = data.pImg[(i + 1)*data.width + j - 1] + data.pImg[(i + 1)*data.width + j] + data.pImg[(i + 1)*data.width + j + 1]
				- (data.pImg[(i - 1)*data.width + j - 1] + data.pImg[(i - 1)*data.width + j] + data.pImg[(i - 1)*data.width + j + 1]);

			const int dy = data.pImg[(i - 1) *data.width + j + 1] + data.pImg[i*data.width + j + 1] + data.pImg[(i - 1)*data.width + j + 1]
				- (data.pImg[(i + 1) *data.width + j - 1] + data.pImg[i*data.width + j - 1] + data.pImg[(i - 1)*data.width + j - 1]);

			if (std::_Max_value(dx, dy) > threadhold)
			{
				img[i * data.width + j] = 1;
			}

		}
	}
	data.pImg = img;
	return data;
}

Sobel

那么有没有比Prewitt更好的算子吗？

那当然有，那就是Sobel算子。

这个算子同样是基于一阶偏导实现的，与Prewitt不同的是，Sobel对中间区域的像素加上了权值，从而降低了边缘模糊程度，从而达到了更好的效果

而最终的梯度大小为

同样是，根据阈值来确定边缘与非边缘

简单的实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
ImageUtil::IMGDATA sobel(ImageUtil::IMGDATA data, const int threadhold)
{
	BYTE *img = new BYTE[data.width * data.height];
	memset(img, 0, data.width * data.height);
	for (int i = 1; i < data.height - 1; i++)
	{
		for (int j = 1; j < data.width - 1; j++)
		{
			const int gx = data[i - 1][j - 1] * -1 + data[i][j - 1] * -2 + data[i + 1][j - 1] * -1 +
				data[i - 1][j + 1] * 1 + data[i][j + 1] * 2 + data[i + 1][j + 1] * 1;

			const int gy = data[i + 1][j - 1] * -1 + data[i + 1][j] * -2 + data[i + 1][j + 1] * -1 +
				data[i - 1][j - 1] * 1 + data[i - 1][j] * 2 + data[i - 1][j + 1] * 1;


			const double g = std::sqrt(gx*gx + gy * gy);
			if (g > threadhold)
				img[i*data.width + j] = 1;
		}
	}
	data.pImg = img;
	return data;
}

LoG

不过，Sobel算子虽然对于图像的噪声有所采取防范措施，但是，总的来说，他还是没有对图像进行滤波的操作。

因此，我们这里采用的方法则是更加高级的方法。

因此，LoG就应运诞生了。

LoG算法是指的是高斯拉普拉斯。简单来说就是，先对图像进行高斯滤波，从而计算出平滑的图像，然后再对图像进行拉普拉斯边缘检测，就可以得出一个比较好的图像了。

这个算法一个比较难的难点就是对图像取二阶高斯滤波。

而结合拉普拉斯算子则是

只要我们根据输入的σ，以及卷积模板的大小，那么就可以计算出LoG算子的模板了。比较常用的近似的是这个

当然，常用的是这个的负模板。

LoG算法的详细步骤可以为

1，对图像取一个$n*n$的高斯低通滤波器对输入图像进行滤波

2，根据第一步取得的图像进行拉普拉斯算子卷积运算

3，根据得出的结果，使用阈值去取得边缘

实现

这里没有直接使用模板去计算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
//NxN高斯卷积模板
double** getGaussianKernel(const int size,const double sqrSigma)
{
	double** gaus = new double*[size];
	for (int i = 0; i < size; i++)
	{
		gaus[i] = new double[size];
	}
	const double pi = 4.0 * std::atan(1.0);
	const int center = size / 2;
	double sum = 0;
	for (int i = 0; i < size; i++)
	{
		for (int j = 0; j < size; j++)
		{
			gaus[i][j] = (1 / (2 * pi*sqrSigma))*exp(-((1 - center)*(1 - center) + (j - center)*(j - center)) / (2 * sqrSigma));
			sum += gaus[i][j];
		}
        
	}

	for (int i = 0; i < size; i++)
	{
		for (int j = 0; j < size; j++)
		{
			gaus[i][j] /= sum;
		}
	}

	return gaus;
}

//LoG
ImageUtil::IMGDATA LOG(ImageUtil::IMGDATA data,double sqrSigma, const int threadhold)
{
	BYTE *img = new BYTE[data.width * data.height];
	memset(img, 0, data.width * data.height);
	double** gaus = getGaussianKernel(5, sqrSigma);
	for (int i = 2; i < data.height - 2; i++)
	{
		for (int j = 2; j < data.width - 2; j++)
		{
			int sum = 0;
			for (int x = -2; x <= 2; x++)
			{
				for (int y = -2; y <= 2; y++)
				{
					sum += data[i + x][j + y] * gaus[4 - (x + 2)][y + 2];
				}
			}
            img[i * data.width + j] = sum;

		}
	}

	for (int i = 0; i < data.height; i++)
	{
		for (int j = 0; j < data.width; j++)
		{
			int up, down, left, right;
			if (i == 0)
				up = 0;
			else
				up = i - 1;

			if (i == data.height - 1)
				down = data.height - 1;
			else
				down = i + 1;

			if (j == 0)
				left = 0;
			else
				left = j - 1;

			if (j == data.width - 1)
				right = data.width - 1;
			else
				right = j + 1;

			img[i * data.width + j] = ImageUtil::clamp(
				1 * data.pImg[up * data.width + left] + 1 * data.pImg[up * data.width + j] + 1 * data.pImg[up * data.width + right] +
				1 * data.pImg[i * data.width + left] + -8 * data.pImg[i * data.width + j] + 1 * data.pImg[i * data.width + right] +
				1 * data.pImg[down * data.width + left] + 1 * data.pImg[down * data.width + j] + 1 * data.pImg[down * data.width + right]);
		}

	}

	for(int i = 0;i < data.width * data.height;i++)
	{
		if (img[i] > threadhold)
			img[i] = 1;
		else
			img[i] = 0;
	}

	data.pImg = img;
	return data;
}

Canny

最后要介绍的，就是这个Canny检测器了。这是这四个检测器当中最优秀的一个。Canny方法有三个最基本的目标

1，低错误率。所有边缘都被找到，并且没有伪相应。

2，边缘点被很好地定位。已定位的边缘尽可能的接近真实边缘。

3，单一响应点。对于真实的边缘点，探测器只返回一个点。

Canny的本质工作就是对以上三个准则进行数学公式化并求其最优解。这通常来说是很困难的。然而，对由加性高斯白噪音污染的一阶台阶边缘使用数字最佳化很好近似高斯一阶导数的。

而把结果推广到二维也同样适用。

那么，我们开始实现Canny方法吧。

降噪

第一步，自然是绕不开降噪的。使用如同LoG算子当中的高斯滤波即可实现降噪的效果了。

寻找图像当中的亮度梯度

第二步，分为两个步骤，首先是求得降噪后的图像的梯度幅度。由于，我们需要知道图像的梯度方向，所以，我们必须将图像的x与y方向的梯度都求出来

在这里，我是用的是sobel算子去求得图像的x方向与y方向的梯度。

第二个步骤，就是对梯度求出来的值进行抑制了。

由于，使用梯度产生的边缘，$M(x,y)$通常在局部最大值当中包含了一个更宽的范围（图像的一阶导数并非只有突变点才是突变值）

因此，我们需要对其中的非最大值进行抑制。

那么我们需要做的就是

1，寻找点$\theta(i,j)$的方向$d_k$

2，若$M(x,y)$的值是沿$d_k$方向上的两个邻居当中的最大值，那就不进行操作。否则，置为0

这样，每一个边缘点就都只剩下最大值的那个了，也符合了Canny的目的之一，单一响应点

在图像当中追踪边缘

在得出了非最大值抑制的边缘之后，我们就需要对图像进行阈值处理，来减少其中的伪边缘点了。

在之前，我们的算法都是使用单阈值点去处理这个问题。但是，这会导致一个问题，那就是，阈值过低会使得伪边缘点过多，而阈值过高却会导致有效边缘点被误删。

因此，我们这里通过使用滞后阈值来试图改变这个状况。

这里使用了两个阈值$T_L,T_H$，根据Canny的建议，这两个阈值的比应该是1：2或者1：3。

首先我们分别对图像进行阈值处理

在阈值处理之后，我们就得出了由两个阈值所得出的两张处理过后的图片。而此时，$G_L$是包含着我们使用较高阈值处理的图片的，因此，我们再对$G_L$进行处理

那么得出得就是一个是强阈值的图案，一个是弱阈值的图案。此时，我们可以以此为基础检测真正的边缘了。

1，首先将所有的$G_H$设定为有效的边缘边

2，从$G_H$当中取得一个还没有访问过的像素点p

3，从$G_L$当中与p相邻的点也被纳入到有效边缘边当中,并从这个点出发，用8联通的方式将其周围的点连接在一起,最后将这个点从$G_L$加入到$G_H$当中

4，若$G_H$遍历完成，结束遍历

那么此时标记出来的所有有效边就是真正的边缘了。

整个Canny方法可以总结为如下的操作

1，使用高斯滤波器平滑图像

2，计算图像的梯度值和角度

3，对梯度值使用非最大值抑制

4，使用双阈值处理以及连接分析来检测并连接边缘

代码的实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
//只用于记录像素点的位置
struct Pixel
{
	int x = 0, y = 0;
	BYTE pix = 0;

	Pixel(const int x, const int y, const BYTE p)
		:x(x),y(y),pix(p)
	{
	}

	Pixel() = default;
};

//canny算法
ImageUtil::IMGDATA canny(ImageUtil::IMGDATA data, const int minVal,const int maxVal)
{

    //步骤一，高斯滤波
	double** gaus = getGaussianKernel(5,0.01);
	for (int i = 2; i < data.height - 2; i++)
	{
		for (int j = 2; j < data.width - 2; j++)
		{
			int sum = 0;
			for (int x = -2; x <= 2; x++)
			{
				for (int y = -2; y <= 2; y++)
				{
					sum += data[i + x][j + y] * gaus[4 - (x + 2)][y + 2];
				}
			}
            data[i][j] = sum;
			
		}
	}

	//步骤二，计算梯度与角度
	BYTE *sobelImg = new BYTE[data.width * data.height];
	int *gxArr = new int[data.width * data.height];
	int *gyArr = new int[data.width * data.height];
	memset(sobelImg, 0, data.width * data.height);
	memset(gxArr, 0, data.width * data.height);
	memset(gyArr, 0, data.width * data.height);

	for (int i = 1; i < data.height - 1; i++)
	{
		for (int j = 1; j < data.width - 1; j++)
		{
			const int gx = data[i - 1][j - 1] * -1 + data[i][j - 1] * -2 + data[i + 1][j - 1] * -1 +
				data[i - 1][j + 1] * 1 + data[i][j + 1] * 2 + data[i + 1][j + 1] * 1;

			const int gy = data[i + 1][j - 1] * -1 + data[i + 1][j] * -2 + data[i + 1][j + 1] * -1 +
				data[i - 1][j - 1] * 1 + data[i - 1][j] * 2 + data[i - 1][j + 1] * 1;

			gxArr[i*data.width + j] = gx;
			gyArr[i*data.width + j] = gy;

			const double g = std::sqrt(gx*gx + gy * gy);
			sobelImg[i*data.width + j] = g;
		}
	}

    //步骤三：非最大值抑制(从计算出来的梯度方向的法线方向才是边缘宽的方向)
	BYTE *temp = new BYTE[data.width * data.height];
	for(int i = 0;i < data.width * data.height;i++)
	{
		temp[i] = sobelImg[i];
	}

	for (int i = 1; i < data.height - 1; i++)
	{
		for (int j = 1; j < data.width - 1; j++)
		{
			double dir;
			if (gxArr[i*data.width + j] == 0)
				dir = 90;
			else
				dir = (std::atan(gyArr[i*data.width + j] / gxArr[i*data.width + j])) * 180 / ImageUtil::pi;
			//水平
			if ((dir >= 157.5 || dir <= -157.5) || (dir <= 22.5 && dir >= -22.5))
			{
				if (sobelImg[i * data.width + j] < sobelImg[i * data.width + j + 1] ||
					sobelImg[i * data.width + j] < sobelImg[i * data.width + j - 1])
				{
					temp[i * data.width + j] = 0;
				}
			}
			//-45度
			else if ((dir <= -112.5 && dir >= -157.5) || (dir <= 67.5 && dir >= 22.5))
			{
				if (sobelImg[i * data.width + j] < sobelImg[(i + 1) * data.width + j - 1] ||
					sobelImg[i * data.width + j] < sobelImg[(i - 1) * data.width + j + 1])
				{
					temp[i * data.width + j] = 0;
				}
			}
			//垂直
			else if ((dir >= -112.5 && dir <= -67.5) || (dir <= 112.5 && dir >= 67.5))
			{
				if (sobelImg[i * data.width + j] < sobelImg[(i + 1) * data.width + j] ||
					sobelImg[i * data.width + j] < sobelImg[(i - 1) * data.width + j])
				{
					temp[i * data.width + j] = 0;
				}
			}
			//+45度
			else if ((dir >= -67.5 && dir <= -22.5) || (dir <= 157.5 && dir >= 112.5))
			{
				if (sobelImg[i * data.width + j] < sobelImg[(i + 1) * data.width + j + 1] ||
					sobelImg[i * data.width + j] < sobelImg[(i - 1) * data.width + j - 1])
				{
					temp[i * data.width + j] = 0;
				}
			}
		}
	}	
	
	for (int i = 0; i < data.width * data.height; i++)
	{
		sobelImg[i] = temp[i];
	}
	delete[] temp;
	delete[] gxArr;
	delete[] gyArr;

	std::queue<Pixel> highPixQue;
	BYTE *lowPix = new BYTE[data.width * data.height];
	//步骤四，双阈值处理
	for (int i = 0; i < data.height; i++)
	{
		for (int j = 0; j < data.width; j++) {
			if (sobelImg[i * data.width + j] > maxVal)
			{
				const Pixel p(j, i, 1);
				highPixQue.push(p);
			}

			if (sobelImg[i * data.width + j] > minVal && sobelImg[i * data.width + j] <= maxVal)
			{
				lowPix[i * data.width + j] = 1;
			}
			else
			{
				lowPix[i * data.width + j] = 0;
			}
		}
	}

    //步骤五，分析并联通边缘
	memset(sobelImg, 0, data.width*data.height);
	while (!highPixQue.empty())
	{
		const Pixel p = highPixQue.front();
		highPixQue.pop();

		sobelImg[p.y * data.width + p.x] = 1;
		for(int i = -1;i <= 1;i++)
		{
			for(int j = -1;j <= 1;j++)
			{
				if(p.y + i < 0 || p.y + i >= data.height || p.x + j < 0 || p.x + j >= data.width)
					continue;

				if (i == 0 && j == 0)
					continue;

				if(lowPix[(p.y + i) * data.width + p.x + j] == 1)
				{
					//8联通合并
					for (int x = -1; x <= 1; x++)
					{
						for (int y = -1; y <= 1; y++)
						{
							if (p.y + i + y < 0 || 
                                				p.y + i + y >= data.height || 
                                				p.x + j + x < 0 || 
                                				p.x + j + x >= data.width)
								continue;

							if (lowPix[(p.y + i + y) * data.width + p.x + j + x] == 1)
							{
								highPixQue.push(Pixel(p.x + j, p.y + i, 1));
							}

								
						}
					}
					lowPix[(p.y + i) * data.width + p.x + j] = 0;
				}

			}
		}
	}

	delete[] lowPix;


	data.pImg = sobelImg;

	for (int i = 0; i < 5; i++)
		delete[] gaus[i];

	delete[] gaus;

	return data;
}
  • 本文作者： ShinyGX
  • 本文链接： https://ShinyGX.github.io/posts/49b27d17/
  • 版权声明： 本博客所有文章除特别声明外，均采用 https://creativecommons.org/licenses/by-nc-sa/3.0/ 许可协议。转载请注明出处！
0%