Ex2:用CImg重写、封装给定的Canny代码,并测试(Code0)
github地址:https://github.com/linjiafengyang/ComputerVision
代码文件
用CImg重写Code0文件夹中的canny.c和canny.h为canny.cpp和canny.h,并增加一个main.cpp文件作为图像数据输入并测试。
重写过程:首先替换原来的图像输入为CImg图像输入,以及对像素点的for循环更改为CImg版本cimg_forXY(),最后把原来的C语言struct结构体改成C++类。
canny.h文件
#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文件
#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文件
#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个步骤:
- 使用高斯滤波器,以平滑图像,滤除噪声。
- 计算图像中每个像素点的梯度强度和方向。
- 应用非极大值(Non-Maximum Suppression)抑制,以消除边缘检测带来的杂散响应
- 应用双阈值(Double-Threshold)检测来确定真实的和潜在的边缘。
- 通过抑制孤立的弱边缘最终完成边缘检测。