Ex2:用CImg重写、封装给定的Canny代码,并测试

计算机视觉作业二

代码文件

用CImg重写Code0文件夹中的canny.c和canny.h为canny.cpp和canny.h,并增加一个main.cpp文件作为图像数据输入并测试。
重写过程:首先替换原来的图像输入为CImg图像输入,以及对像素点的for循环更改为CImg版本cimg_forXY(),最后把原来的C语言struct结构体改成C++类。

canny.h文件

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
#ifndef canny_h
#define canny_h

#include "CImg.h"
using namespace cimg_library;

class Canny
{
public:
CImg<unsigned char> canny(CImg<unsigned char> grey, int width, int height);
CImg<unsigned char> cannyparam(CImg<unsigned char> grey, int width, int height,
float lowthreshold, float highthreshold,
float gaussiankernelradius, int gaussiankernelwidth,
int contrastnormalised);

Canny *allocatebuffers(const CImg<unsigned char> & grey, int width, int height);
void killbuffers(Canny *can);
int computeGradients(Canny *can, float kernelRadius, int kernelWidth);
void performHysteresis(Canny *can, int low, int high);
void follow(Canny *can, int x1, int y1, int i1, int threshold);
void normalizeContrast(CImg<unsigned char> & data, int width, int height);
float hypotenuse(float x, float y);
float gaussian(float x, float sigma);

private:
CImg<unsigned char> data; /* input image */
int width;
int height;
int *idata; /* output for edges */
int *magnitude; /* edge magnitude as detected by Gaussians */
float *xConv; /* temporary for convolution in x direction */
float *yConv; /* temporary for convolution in y direction */
float *xGradient; /* gradients in x direction, as detected by Gaussians */
float *yGradient; /* gradients in x direction,a s detected by Gaussians */

};

#endif

canny.cpp文件

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
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
#include "canny.h"
#include <math.h>

#define ffabs(x) ( (x) >= 0 ? (x) : -(x) )
#define GAUSSIAN_CUT_OFF 0.005f
#define MAGNITUDE_SCALE 100.0f
#define MAGNITUDE_LIMIT 1000.0f
#define MAGNITUDE_MAX ((int) (MAGNITUDE_SCALE * MAGNITUDE_LIMIT))

/*
Canny edge detection with default parameters
Params: grey - the greyscale image
width, height - image width and height
Returns: binary image with edges as set pixels
*/
CImg<unsigned char> Canny::canny(CImg<unsigned char> grey, int width, int height)
{
return cannyparam(grey, width, height, 2.5f, 7.5f, 2.0f, 16, 0);
}

/*
Canny edge detection with parameters passed in by user
Params: grey - the greyscale image
width, height - image dimensions
lowthreshold - default 2.5
highthreshold - default 7.5
gaussiankernelradius - radius of edge detection Gaussian, in standard deviations
(default 2.0)
gaussiankernelwidth - width of Gaussian kernel, in pixels (default 16)
contrastnormalised - flag to normalise image before edge detection (defualt 0)
Returns: binary image with set pixels as edges
*/
CImg<unsigned char> Canny::cannyparam(CImg<unsigned char> grey, int width, int height,
float lowthreshold, float highthreshold,
float gaussiankernelradius, int gaussiankernelwidth,
int contrastnormalised)
{
Canny *can = 0;
CImg<unsigned char> answer = CImg<unsigned char>(width, height, 1, 1, 0);
int low, high;
int err;
int i;

can = allocatebuffers(grey, width, height);

if (!can)
goto error_exit;
if (contrastnormalised)
normalizeContrast(can->data, width, height);

err = computeGradients(can, gaussiankernelradius, gaussiankernelwidth);

if (err < 0)
goto error_exit;

low = (int)(lowthreshold * MAGNITUDE_SCALE + 0.5f);
high = (int)(highthreshold * MAGNITUDE_SCALE + 0.5f);
performHysteresis(can, low, high);

cimg_forXY(answer, x, y)
{
i = y * width + x;
answer(x, y, 0) = can->idata[i] > 0 ? 255 : 0;
}

killbuffers(can);
return answer;
error_exit:
free(answer);
killbuffers(can);
return CImg<unsigned char>();
}

/*
buffer allocation
*/
Canny * Canny::allocatebuffers(const CImg<unsigned char> & grey, int width, int height)
{
Canny *answer;

answer = new Canny;
if (!answer)
goto error_exit;
answer->data = CImg<unsigned char>(width, height, 1, 1, 0);
answer->idata = new int[width * height];
answer->magnitude = new int[width * height];
answer->xConv = new float[width * height];
answer->yConv = new float[width * height];
answer->xGradient = new float[width * height];
answer->yGradient = new float[width * height];

if (!answer->data || !answer->idata || !answer->magnitude ||
!answer->xConv || !answer->yConv ||
!answer->xGradient || !answer->yGradient)
goto error_exit;

cimg_forXY(grey, x, y)
{
answer->data(x, y, 0) = grey(x, y, 0);
}
answer->width = width;
answer->height = height;

return answer;
error_exit:
killbuffers(answer);
return 0;
}

/*
buffers destructor
*/
void Canny::killbuffers(Canny *can) {
if (can)
{
delete(can->idata);
delete(can->magnitude);
delete(can->xConv);
delete(can->yConv);
delete(can->xGradient);
delete(can->yGradient);
}
}

/* NOTE: The elements of the method below (specifically the technique for
non-maximal suppression and the technique for gradient computation)
are derived from an implementation posted in the following forum (with the
clear intent of others using the code):
http://forum.java.sun.com/thread.jspa?threadID=546211&start=45&tstart=0
My code effectively mimics the algorithm exhibited above.
Since I don't know the providence of the code that was posted it is a
possibility (though I think a very remote one) that this code violates
someone's intellectual property rights. If this concerns you feel free to
contact me for an alternative, though less efficient, implementation.
*/
int Canny::computeGradients(Canny *can, float kernelRadius, int kernelWidth)
{
float *kernel;
float *diffKernel;
int kwidth;

int width, height;

int initX;
int maxX;
int initY;
int maxY;

int x, y;
int i;
int flag;

width = can->width;
height = can->height;

kernel = new float[kernelWidth];
diffKernel = new float[kernelWidth];

if (!kernel || !diffKernel)
goto error_exit;

/* initialise the Gaussian kernel */
for (kwidth = 0; kwidth < kernelWidth; kwidth++)
{
float g1, g2, g3;
g1 = gaussian((float)kwidth, kernelRadius);
if (g1 <= GAUSSIAN_CUT_OFF && kwidth >= 2)
break;
g2 = gaussian(kwidth - 0.5f, kernelRadius);
g3 = gaussian(kwidth + 0.5f, kernelRadius);
kernel[kwidth] = (g1 + g2 + g3) / 3.0f / (2.0f * (float) 3.14 * kernelRadius * kernelRadius);
diffKernel[kwidth] = g3 - g2;
}

initX = kwidth - 1;
maxX = width - (kwidth - 1);
initY = width * (kwidth - 1);
maxY = width * (height - (kwidth - 1));

/* perform convolution in x and y directions */
for (x = initX; x < maxX; x++)
{
for (y = initY; y < maxY; y += width)
{
int index = x + y;
float sumX = can->data[index] * kernel[0];
float sumY = sumX;
int xOffset = 1;
int yOffset = width;
while (xOffset < kwidth)
{
sumY += kernel[xOffset] * (can->data[index - yOffset] + can->data[index + yOffset]);
sumX += kernel[xOffset] * (can->data[index - xOffset] + can->data[index + xOffset]);
yOffset += width;
xOffset++;
}

can->yConv[index] = sumY;
can->xConv[index] = sumX;
}
}

for (x = initX; x < maxX; x++)
{
for (y = initY; y < maxY; y += width)
{
float sum = 0.0f;
int index = x + y;
for (i = 1; i < kwidth; i++)
sum += diffKernel[i] * (can->yConv[index - i] - can->yConv[index + i]);

can->xGradient[index] = sum;
}
}

for (x = kwidth; x < width - kwidth; x++)
{
for (y = initY; y < maxY; y += width)
{
float sum = 0.0f;
int index = x + y;
int yOffset = width;
for (i = 1; i < kwidth; i++)
{
sum += diffKernel[i] * (can->xConv[index - yOffset] - can->xConv[index + yOffset]);
yOffset += width;
}

can->yGradient[index] = sum;
}
}

initX = kwidth;
maxX = width - kwidth;
initY = width * kwidth;
maxY = width * (height - kwidth);
for (x = initX; x < maxX; x++)
{
for (y = initY; y < maxY; y += width)
{
int index = x + y;
int indexN = index - width;
int indexS = index + width;
int indexW = index - 1;
int indexE = index + 1;
int indexNW = indexN - 1;
int indexNE = indexN + 1;
int indexSW = indexS - 1;
int indexSE = indexS + 1;

float xGrad = can->xGradient[index];
float yGrad = can->yGradient[index];
float gradMag = hypotenuse(xGrad, yGrad);

/* perform non-maximal supression */
float nMag = hypotenuse(can->xGradient[indexN], can->yGradient[indexN]);
float sMag = hypotenuse(can->xGradient[indexS], can->yGradient[indexS]);
float wMag = hypotenuse(can->xGradient[indexW], can->yGradient[indexW]);
float eMag = hypotenuse(can->xGradient[indexE], can->yGradient[indexE]);
float neMag = hypotenuse(can->xGradient[indexNE], can->yGradient[indexNE]);
float seMag = hypotenuse(can->xGradient[indexSE], can->yGradient[indexSE]);
float swMag = hypotenuse(can->xGradient[indexSW], can->yGradient[indexSW]);
float nwMag = hypotenuse(can->xGradient[indexNW], can->yGradient[indexNW]);
float tmp;
/*
* An explanation of what's happening here, for those who want
* to understand the source: This performs the "non-maximal
* supression" phase of the Canny edge detection in which we
* need to compare the gradient magnitude to that in the
* direction of the gradient; only if the value is a local
* maximum do we consider the point as an edge candidate.
*
* We need to break the comparison into a number of different
* cases depending on the gradient direction so that the
* appropriate values can be used. To avoid computing the
* gradient direction, we use two simple comparisons: first we
* check that the partial derivatives have the same sign (1)
* and then we check which is larger (2). As a consequence, we
* have reduced the problem to one of four identical cases that
* each test the central gradient magnitude against the values at
* two points with 'identical support'; what this means is that
* the geometry required to accurately interpolate the magnitude
* of gradient function at those points has an identical
* geometry (upto right-angled-rotation/reflection).
*
* When comparing the central gradient to the two interpolated
* values, we avoid performing any divisions by multiplying both
* sides of each inequality by the greater of the two partial
* derivatives. The common comparand is stored in a temporary
* variable (3) and reused in the mirror case (4).
*
*/
flag = ((xGrad * yGrad <= 0.0f) /*(1)*/
? ffabs(xGrad) >= ffabs(yGrad) /*(2)*/
? (tmp = ffabs(xGrad * gradMag)) >= ffabs(yGrad * neMag - (xGrad + yGrad) * eMag) /*(3)*/
&& tmp > fabs(yGrad * swMag - (xGrad + yGrad) * wMag) /*(4)*/
: (tmp = ffabs(yGrad * gradMag)) >= ffabs(xGrad * neMag - (yGrad + xGrad) * nMag) /*(3)*/
&& tmp > ffabs(xGrad * swMag - (yGrad + xGrad) * sMag) /*(4)*/
: ffabs(xGrad) >= ffabs(yGrad) /*(2)*/
? (tmp = ffabs(xGrad * gradMag)) >= ffabs(yGrad * seMag + (xGrad - yGrad) * eMag) /*(3)*/
&& tmp > ffabs(yGrad * nwMag + (xGrad - yGrad) * wMag) /*(4)*/
: (tmp = ffabs(yGrad * gradMag)) >= ffabs(xGrad * seMag + (yGrad - xGrad) * sMag) /*(3)*/
&& tmp > ffabs(xGrad * nwMag + (yGrad - xGrad) * nMag) /*(4)*/
);
if (flag)
{
can->magnitude[index] = (gradMag >= MAGNITUDE_LIMIT) ? MAGNITUDE_MAX : (int)(MAGNITUDE_SCALE * gradMag);
/*NOTE: The orientation of the edge is not employed by this
implementation. It is a simple matter to compute it at
this point as: Math.atan2(yGrad, xGrad); */
}
else
{
can->magnitude[index] = 0;
}
}
}
free(kernel);
free(diffKernel);
return 0;
error_exit:
free(kernel);
free(diffKernel);
return -1;
}

/*
we follow edges. high gives the parameter for starting an edge,
how the parameter for continuing it.
*/
void Canny::performHysteresis(Canny *can, int low, int high)
{
int offset = 0;
int x, y;

memset(can->idata, 0, can->width * can->height * sizeof(int));

for (y = 0; y < can->height; y++)
{
for (x = 0; x < can->width; x++) {
if (can->idata[offset] == 0 && can->magnitude[offset] >= high)
follow(can, x, y, offset, low);
offset++;
}
}
}

/*
recursive portion of edge follower
*/
void Canny::follow(Canny *can, int x1, int y1, int i1, int threshold)
{
int x, y;
int x0 = x1 == 0 ? x1 : x1 - 1;
int x2 = x1 == can->width - 1 ? x1 : x1 + 1;
int y0 = y1 == 0 ? y1 : y1 - 1;
int y2 = y1 == can->height - 1 ? y1 : y1 + 1;

can->idata[i1] = can->magnitude[i1];
for (x = x0; x <= x2; x++)
{
for (y = y0; y <= y2; y++)
{
int i2 = x + y * can->width;
if ((y != y1 || x != x1) && can->idata[i2] == 0 && can->magnitude[i2] >= threshold)
follow(can, x, y, i2, threshold);
}
}
}

void Canny::normalizeContrast(CImg<unsigned char> & data, int width, int height)
{
int histogram[256] = { 0 };
int remap[256];
int sum = 0;
int j = 0;
int k;
int target;
int i;

cimg_forXY(data, x, y)
{
histogram[data(x, y, 0)]++;
}

for (i = 0; i < 256; i++)
{
sum += histogram[i];
target = (sum * 255) / (width * height);
for (k = j + 1; k <= target; k++)
remap[k] = i;
j = target;
}

cimg_forXY(data, x, y)
{
data(x, y, 0) = remap[data(x, y, 0)];
}
}

float Canny::hypotenuse(float x, float y)
{
return (float)sqrt(x*x + y*y);
}

float Canny::gaussian(float x, float sigma)
{
return (float)exp(-(x * x) / (2.0f * sigma * sigma));
}

main.cpp文件

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#include "canny.cpp"

int main() {
CImg<unsigned char> bigben_grey("../test_Data/bigben.bmp");
Canny c;
CImg<unsigned char> bigben_canny = c.canny(bigben_grey, bigben_grey.width(), bigben_grey.height());
bigben_canny.save("../result_Data/bigben_canny.bmp");

CImg<unsigned char> lena_grey("../test_Data/lena.bmp");
CImg<unsigned char> lena_canny = c.canny(lena_grey, lena_grey.width(), lena_grey.height());
lena_canny.save("../result_Data/lena_canny.bmp");
//lena_canny.display("lena_canny");

CImg<unsigned char> stpietro_grey("../test_Data/stpietro.bmp");
CImg<unsigned char> stpietro_canny = c.canny(stpietro_grey, stpietro_grey.width(), stpietro_grey.height());
stpietro_canny.save("../result_Data/stpietro_canny.bmp");

CImg<unsigned char> twows_grey("../test_Data/twows.bmp");
CImg<unsigned char> twows_canny = c.canny(twows_grey, twows_grey.width(), twows_grey.height());
twows_canny.save("../result_Data/twows_canny.bmp");
return 0;
}

测试数据

作业中有四张bmp图像,但这里限于篇幅,我以女神lena为例:
这里写图片描述

测试结果

这里写图片描述


修改算法参数后的测试结果

注意到程序中这个函数接收Canny边缘检测算法的相关参数,也就是调参主要是针对这几个参数:其中width和height为输入图像的宽度和高度,无需调整,所以一共有五个参数需要调整,分别为lowthreshold低阈值(默认2.5)、highthreshold高阈值(默认7.5)、gaussiankernelradius标准差(默认2.0)、gaussiankernelwidth高斯卷积核大小(默认16)、contrastnormalised图像是否需要归一化(默认0值表示不需要)。

lowthreshold低阈值

保持其他参数不变,修改低阈值后的测试结果如下:默认2.5。
这里写图片描述

分析:

可以看出,我们修改原程序中的lowthreshold低阈值后,比如改成0.1,低阈值减小,将会增加很多噪声;而改成6.0,低阈值增大,又会丢失很多强边缘像素。这与已知结论吻合:对于弱边缘像素(在低阈值和高阈值之间),这些像素可以从真实边缘提取也可以是因噪声或颜色变化引起的。为了获得准确的结果,应该抑制由后者引起的弱边缘像素。因此需要调整低阈值参数。


highthreshold高阈值

保持其他参数不变,修改高阈值后的测试结果如下:默认7.5。
这里写图片描述

分析:

可以看出,我们修改原程序中的highthreshold高阈值后,比如改成3.0,高阈值减小,此时将会增加很多强边缘像素,因为大于3的像素点将被确定为边缘;而改成10.0,高阈值增大,原来的强边缘像素大部分会转化为弱边缘像素,丢失一部分边缘像素点。这与已知结论吻合:被划分为强边缘的像素点已经被确定为边缘,因为它们是从图像中的真实边缘中提取出来的。


gaussiankernelradius标准差

保持其他参数不变,修改标准差后的测试结果如下:默认2.0。
这里写图片描述

分析:

可以看出,我们修改原程序中的gaussiankernelradius标准差后,比如改成1.0,标准差减小,为标准正态分布,此时将会增加很多噪声,高斯滤波平滑效果不好,表现为噪声过多;而改成4.0,标准差增大,原来的强边缘像素大部分被误认为噪声,因此将会被丢弃,造成边缘缺失,高斯滤波平滑效果也不好。这与已知结论吻合:为了平滑图像,使用高斯滤波器与图像进行卷积,该步骤将平滑图像,以减少边缘检测器上明显的噪声影响。


gaussiankernelwidth高斯卷积核大小

保持其他参数不变,修改高斯卷积核大小后的测试结果如下:默认16。
这里写图片描述
这里写图片描述

分析:

可以看出,我们修改原程序中的gaussiankernelwidth高斯卷积核大小后,比如改成7,发现两张边缘检测图其实没变化,说明原程序中高斯卷积核大小16已经过大,将会浪费内存空间,其实为7就有同样的效果了;而改成5,高斯卷积核大小为5也是比较推荐的折中,再看改成3,此时丢失了很多边缘像素,即误认为噪声,对噪声的敏感度较大。这与已知结论吻合:高斯卷积核大小越大,检测器对噪声的敏感度越低,但是边缘检测的定位误差也将略有增加。


contrastnormalised图像是否需要归一化

保持其他参数不变,修改该参数后的测试结果如下:默认0。
这里写图片描述

分析:

可以看出,我们修改原程序中的contrastnormalised为1后,此时程序将会先对图像矩阵进行归一化处理,对照原图,我发现归一化会丢失一些比较暗、比较模糊的像素,在这种情况下做边缘检测时这些像素点都会被忽略,不会被认为是边缘。


Canny边缘检测算法总结

Canny边缘检测算法可以分为以下5个步骤:

  1. 使用高斯滤波器,以平滑图像,滤除噪声。
  2. 计算图像中每个像素点的梯度强度和方向。
  3. 应用非极大值(Non-Maximum Suppression)抑制,以消除边缘检测带来的杂散响应
  4. 应用双阈值(Double-Threshold)检测来确定真实的和潜在的边缘。
  5. 通过抑制孤立的弱边缘最终完成边缘检测。
------本文结束感谢阅读------