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文件

#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个步骤：

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