计算机视觉与模式识别(1)—— A4纸边缘提取

最近跟着老师学习了计算机视觉领域中，一个教科书般的知识点 —— 图形的边缘提取。随后自己写了一个A4纸边缘提取的小程序（200行），觉得蛮有意思的，测试的效果也不错，于是心血来潮，打算把自己整个设计的感想写下来。代码可能不是很规范，因为精力都侧重在理论知识上了，本章的全部难点几乎都在图形处理的理论知识+_+！

那么首先说说这次的目标是什么：

输入：普通的A4打印纸，上面可能有手写或打印字体，而且拍照时角度不正。

输出：

1. 图像的边缘

2. 计算A4纸边缘的各直线方程

3. 提取A4纸的4个角点

我的环境：

Windows10、C++11、Visual Studio2013

有很多图形处理的相关库文件，这里我使用了老师推荐的CImg库，可以在官网下载：http://www.cimg.eu/，库的使用十分简单，只需将解压文件中的CImg.h头文件加入工程中即可。更多的有关库的使用这里不过多描述，在后面代码讲解中，我会对关键的函数进行解释。事实上，CImg作为专业的图形库，自身已经实现了类似边缘提取的函数，不过我们这里不会使用，因为这样不就没有意义了吗╮(╯▽╰)╭。

话不多说，从输入来看要提取出A4纸的边缘绝非易事。从我们最直观的感觉来看，首先我们要区分颜色，然后把白色的那块圈起来就得到我们的边缘了。说实话，我一开始也是这么想的，典型的小白思想。注意上面的第一张图，A4纸中会有手写体，所以并不是全白的，而且由于拍照时的光线原因，纸上还有阴影。除此以外，图像中还存在其他物体，比如桌子、地板、腿等。第二张图中还被另一张A4纸抢镜，我们不能把它也算进来。

卧cao( ‵o′)，那还怎么做！

莫慌！但不用抱紧我。现在把我的姿势告诉你：

1. 把彩色图像灰度化。

2. 高斯模糊一下。

3. 计算灰度图的梯度。

4. 把梯度高的坐标映射到Hough图坐标上。（Hough转换，一种Voting算法，关键！）

5. 提取Hough图峰值

6. 计算直线方程和角点

经过这一系列的处理，我们就可以得到这样的结果：

效果还不错吧，没有出现上述任何一项问题，A4纸4条边缘和角点的位置也十分贴近真实情况。

我们一步步来。

一、彩色图像灰度化处理。

首先读入BMP格式的原图片：

#include "CImg.h"
using namespace cimg_library;

#define FILE "A4.bmp"

int main {
	CImg<double> image(FILE);     // 原始图
	CImg<double> grayImg(image);  // 灰度图

	image.display;
}

CImg类是CImg.h中的主要图形类，存储了图片的所有数据，包括长宽高，色彩通道和像素值等。它的初始化很简单，可以是文件路径或者是自身对象，也可以通过指定长宽高创建一张默认图片。CImg只是默认支持了BMP格式的图像，如果需要读取JPG等其它格式的图片，要配置GraphicsMagick或者ImageMagick，这里不多讲。

接下来是遍历image的所有像素，将它们变为灰度值。灰度值的计算公式：

Gray = (R * 299 + G * 587 + B * 114 + 500) / 1000;

	// Gray scale;
	cimg_forXY(grayImg, x, y) {
		double R = grayImg(x, y, 0, 0);
		double G = grayImg(x, y, 0, 1);
		double B = grayImg(x, y, 0, 2);
		double Gray = (R * 299 + G * 587 + B * 114 + 500) / 1000;
		grayImg(x, y, 0, 0) = Gray;
		grayImg(x, y, 0, 1) = Gray;
		grayImg(x, y, 0, 2) = Gray;
	}
	grayImg.display;

cimg_forXY函数等价于：

	for (int y = 0; y < grayImg.height; ++y) {
		for (int x = 0; x < grayImg.width; ++x) {
			// do something
		}
	}

由于每个坐标的像素有RGB三个色彩通道，所以需要将所有通道设置成灰度值。

完成这部分代码，运行后应有如下结果：

二、高斯模糊。

注意到灰度图中有很多噪声点，这不利于待会儿的梯度计算。梯度计算的目的就是把色差较大的像素找出来，通过这种方法找出边缘（注意我没说A4纸的边缘），我们不希望把噪声点也当做边缘。CImg里面有自带的模糊函数，我只需要调用一下就好了，定义一个BLUR宏变量，方便日后修改。

#define BLUR 2	

	grayImg.blur(BLUR);   

三、计算灰度图的梯度。

接下来就涉及专业的理论知识了，如何计算一个像素点的梯度呢？我要介绍一下Prewitt算子。Prewitt算子是一个3*3的矩阵，它长这样：

横向模板 纵向模板：

3*3矩阵的位置描述：9个位置依次是Ipp、Icp、Inp、Ipc、Icc、Inc、Ipn、Icn、Inn。

（p表示previous，c表示current，n表示next，Ipp = img(x-1,y-1)，Icn=img(x,y+1)，...）

计算一个点的梯度首先要找出它的邻域矩阵，此时被计算的点位于邻域矩阵的中心，如果是3*3邻域矩阵，则被计算的点位于Icc位置，其它8个点是被计算点的相邻点。通过Prewitt算子，我们可以计算出点的横向平均梯度是(Inc - Ipc)，纵向平均梯度是(Icp - Icn)。

梯度的值 = sqrt ( (Inc-Ipc)^2 + (Icp-Icn)^2 )

#define GRADLIMIT 20

	CImg<double> gradnum(image.width, image.height, 1, 1, 0);
	// 定义3*3领域矩阵I
	CImg_3x3(I, double);
	// 遍历计算梯度值
	cimg_for3x3(grayImg, x, y, 0, 0, I, double) {
		const double ix = Inc - Ipc;
		const double iy = Icp - Icn;
		double grad = std::sqrt(ix*ix + iy*iy);
		// 梯度大于阈值才赋值
		if (grad > GRADLIMIT) {
			gradnum(x, y) = grad;
		}
	}
	gradnum.display;

CImg_3x3是CImg提供的定义领域矩阵的函数，第一个参数是邻域矩阵的名字，第二个参数是邻域矩阵的数据类型。通过将邻域矩阵传入cimg_for3x3循环函数，数据就会自动保存在 Ipp~Inn 九个变量中。我们只打算找出梯度值较大的像素点，因此设置了一个阈值GRADLIMIT，声明为宏变量，只有比阈值大的点才能显示在新创建的图片gradnum里。注意gradnum初始化时只有一个色彩通道，这样就只用修改一个通道值。

完成这部分代码，运行后应有如下结果：

四、梯度图映射到Hough图（Voting算法）

细心的你可能已经发现，经过梯度处理，我们得到的图片已经得到了A4纸的大概边缘，并且消去了很多字体的边缘和其它物体的边缘。但是，仍有一些顽强的噪声点幸存下来了，它们可能是纸上的浓墨重彩，也可能是另一些边缘明显的物体，比如另一张A4纸的一角。因此我们不能利用最小二乘法或者线性回归方法来拟合A4纸的边缘，这些方法受干扰点的影响很大。

伟大的Paul Hough提出了一种解决的办法，称为Hough变换，中文叫霍夫变换。它的主要思路是通过投票的方式，让原始数据点们自己选出一个可以代表他们的函数表达式，因此是一种投票算法。说得这么玄，不如举个例子来说明：

假设我们现在有这么一些数据点：

我们知道直线的表达式形式是：y = mx + b。通常我们把x和y看成自变量，把m和b看成是常数。

现在我们做一次变换，把x和y看成常数，把m和b看成变量：b = -xm + y

由于我们已经获得了所有的坐标（x，y） ，每个坐标（x，y）代进表达式 b = -xm + y 都表示一条关于m和b的直线，因此可以在mb坐标系上画出一簇直线：

这些直线会近乎相交于一点（注意是近乎不是完全），这一点的坐标（1，1）就是原始点们选举出来的值，这个值能最好地表示原始点的直线表达式，也就是：y = x + 1。

为什么原始图中，右下方的噪声点没有影响到直线表达式呢？我们可以看到，噪声点在mb坐标系上表示的直线偏离了大众交点，它固然也会和其它直线相交，但是那些交点的值（每有一条直线经过就加一），远远比不上大众交点的值，因此被毫不留情地淘汰掉了，这不正体现了投票的精髓吗——少数服从多数。

回到我们的程序，现实是，我们不可能定义一个无限长的m轴和b轴，因此需要把直角坐标系转换为极坐标系。这样就有：p = x*cos(a)+y*sin(a)

极坐标系的好处是，角度a的定义域是 [0, 360] 度，距离p的定义域是 [0, sqrt(x^2+y^2)]。

为了完成hough变换，首先创建一个新的CImg对象hough，用来记录点的值：

	int maxDistance = distance(image.width, image.height);
	CImg<double> hough(360, maxDistance, 1, 1, 0);

distance的函数定义如下：

// 计算两点间距离
double distance(double x, double y) {
	return sqrt(x*x + y*y);
}

修改计算梯度的代码，因为我们只对梯度值大于阈值的像素点感兴趣。找出该点后在hough图上对经过 p = x*cos(a)+y*sin(a) 的点值加1：

	// 定义3*3领域矩阵I
	CImg_3x3(I, double);
	// 遍历计算梯度值
	cimg_for3x3(grayImg, x, y, 0, 0, I, double) {
		const double ix = Inc - Ipc;
		const double iy = Icp - Icn;
		double grad = std::sqrt(ix*ix + iy*iy);
		// 梯度大于阈值才赋值
		if (grad > GRADLIMIT) {
			gradnum(x, y) = grad;
			cimg_forX(hough, angle) {      // 遍历x轴计算y，[x,y]+=1
				double rangle = (double)angle*PI / 180.0;
				int polar = (int)(x*cos(rangle) + y*sin(rangle));
				if (polar >= 0 && polar < hough.height) {      // 确定y在图像高度内
					hough(angle, polar) += 1;
				}
			}
		}
	}

显示hough图：

	hough.display;

完成这部分代码，运行后应有如下结果：（右边是局部放大）

五、提取Hough峰值。

在hough图上，我们可以找到4个亮点，那就是我们希望得到的4条直线。但是每个亮点放大后看，都不是规规矩矩的一个点，而是一簇相邻的亮点，所以要进行提取。Hough峰值的提取有很多种办法，比较理想的是质心算法。我使用了一种较为简单的算法，只是单纯的找出范围内最亮的点，这个范围通过一个阈值DIFF表示：

#define DIFF 200

我希望把找到的峰值点信息都保存起来，所以使用了STL的Vector类，并且定义了一个新结构体Dot：

// 2D Dot(x, y) [value]  
struct Dot {
	int x, y, value;
	Dot(int _x, int _y, int _value)
		:x(_x), y(_y), value(_value){}
};

另外要注意，我们的峰值所表示的关于x和y的直线方程并不一定合法，一般是由极坐标转换回直线坐标的正负误差等造成的。因此首先要进行一番筛选，判断直线是否经过我们的图像。我的判断方法是，如果直线与图像任意一条边有交点，则直线是合法的。该方法涉及两个函数CrossX和CrossY，分别计算直线在指定x或者y上的交点。先看峰值提取代码：

// Find peaks
	std::vector<Dot*> peaks;
	cimg_forXY(hough, angle, polar) {
		if (hough(angle, polar) > THRESHOLD) { // 是否是峰值
			bool flag = false;
			const int ymin = 0;
			const int ymax = image.height - 1;
			const int x0 = CrossY(angle, polar, ymin);
			const int x1 = CrossY(angle, polar, ymax);

			const int xmin = 0;
			const int xmax = image.width - 1;
			const int y0 = CrossX(angle, polar, xmin);
			const int y1 = CrossX(angle, polar, xmax);

			if (x0 >= 0 && x0 <= xmax ||				// 表示的直线是否在图像内
				x1 >= 0 && x1 <= xmax ||
				y0 >= 0 && y0 <= ymax ||
				y1 >= 0 && y1 <= ymax) {
				for (int i = 0; i < peaks.size; ++i) {    // 遍历数组，找相邻峰值
					if (distance(peaks[i]->x - angle, peaks[i]->y - polar) < DIFF) {   // 存在相邻峰值
						flag = true;
						if (peaks[i]->value < hough(angle, polar)) {	// 如果比相邻峰值还大
							peaks[i] = new Dot(angle, polar, hough(angle, polar));   // 替换为当前峰值
						}
					}
				}
				if (flag == false) {		// 当前峰值无相邻峰值
					peaks.push_back(new Dot(angle, polar, hough(angle, polar)));   // 加入新峰值
				}
			}
		}
	}

CrossX和CrossY的定义如下：

// Polar coordinate intersection at x
const int CrossX(int theta, int distance, int x) {
	double angle = (double)theta*PI / 180.0;
	double m = -cos(angle) / sin(angle);
	double b = (double)distance / sin(angle);
	return m*x + b;
}

// Polar coordinate intersection at y
const int CrossY(int theta, int distance, int y) {
	double angle = (double)theta*PI / 180.0;
	double m = -cos(angle) / sin(angle);
	double b = (double)distance / sin(angle);
	return ((double)(y - b) / m);
}

如果你的阈值设置得当，那么你会得到一个含有4个值的峰值数组。

六、计算直线方程和角点坐标。

有了峰值数组，实际上我们已经得到了A4纸4条边缘的极坐标直线方程，我们要把它们转换成直角坐标方程，并且计算直线的交点：

// Transform polar coordinates to rectangular coordinates
	std::vector<Line*> lines;
	for (int i = 0; i < peaks.size; ++i) {
		double angle = (double)peaks[i]->x*PI / 180.0;
		double m = -cos(angle) / sin(angle);
		double b = (double)peaks[i]->y / sin(angle);
		lines.push_back(new Line(m, b));
		std::cout << "y = (" << m << ") x + (" << b << ")" << std::endl;
	}

	std::cout << std::endl << "intersections: " << std::endl;

	// Calculate line intersections
	std::vector<Dot*> intersections;
	for (int i = 0; i < lines.size; ++i) {
		for (int j = i + 1; j < lines.size; ++j) {
			double m0 = lines[i]->m;
			double m1 = lines[j]->m;
			double b0 = lines[i]->b;
			double b1 = lines[j]->b;
			double x = (b1 - b0) / (m0 - m1);
			double y = (m0*b1 - m1*b0) / (m0 - m1);
			if (x >= 0 && x < result.width && y >= 0 && y < result.height) {
				intersections.push_back(new Dot(x, y, 0));
				std::cout << "(" << x << ", " << y << ")" << std::endl;
			}
		}
	}

其中Line是新定义的结构体：

// 2D Straight Line(x, y) [y = m*x + b]
struct Line {
	double m, b;
	Line(double _m, double _b)
		:m(_m), b(_b){}
};

完成这部分代码，运行后可以在控制台看到直线方程和角点坐标：

Oops！看来之前只找出了3条直线，不过不要紧，只需改一下阈值就好，后面会详细讲解如何调阈值。

七、绘图。

有了直线方程和角点坐标，我想把他们绘制在原图上，这样可以有直观的感受，知道是否出现了问题。直线和点的绘制可以直接使用CImg的draw_line和draw_circle函数：

	CImg<double> result(image);

	// Draw lines
	for (int i = 0; i < lines.size; ++i) {
		const int ymin = 0;
		const int ymax = result.height - 1;
		const int x0 = (double)(ymin - lines[i]->b) / lines[i]->m;
		const int x1 = (double)(ymax - lines[i]->b) / lines[i]->m;

		const int xmin = 0;
		const int xmax = result.width - 1;
		const int y0 = xmin*lines[i]->m + lines[i]->b;
		const int y1 = xmax*lines[i]->m + lines[i]->b;

		const double color = { 255, 0, 255 };

		if (abs(lines[i]->m) > SLOPE_FLAG) {
			result.draw_line(x0, ymin, x1, ymax, color);
		}
		else {
			result.draw_line(xmin, y0, xmax, y1, color);
		}
	}

	// Draw intersections
	for (int i = 0; i < intersections.size; ++i) {
		const double color = { 255, 100, 255 };
		result.draw_circle(intersections[i]->x, intersections[i]->y, 50, color);
	}

	result.display;

注意特殊情况，要考虑直线是水平或者垂直情况下的绘制。

完成绘制后，我们终于可以看见结果了：

正如我们第6步结束后看到的，我们只找到了A4纸的3条边缘，左侧的边缘似乎被遗弃掉了。这是由于阈值设置不得当引起的，如果足够细心的话，会发现原图A4纸左侧是有阴影的，因此得到的梯度图中，左侧的点比较少，所以在hough图上亮点就不明显。

目前为止，整个程序的阈值有：

#define BLUR 3
#define GRADLIMIT 20
#define THRESHOLD 800
#define DIFF 200
#define SLOPE_FLAG 1

其中SLOPE_FLAG是用来判断直线是否水平的，因此可以不用管它。我们来讨论剩下的4个阈值。

BLUR：模糊参数，通常在原图噪声点较多的时候适当增大。

GRADLIMIT：梯度阈值，如果边缘与周围物体的区分不明显，则应适当减小，反之则适当增大。

THRESHOLD：Hough峰值阈值，这里应通过调整阈值来获取到刚好4个峰值。

DIFF：峰值提取时判断相邻峰值的距离阈值，如果一条边提取后有两个斜率十分相近的直线方程，那么应增大该阈值。

这里，我们应该减小THRESHOLD阈值，从原先的800降到650，这样就能得到正确的结果：

直线的方程和角点坐标显示在命令行窗口：

结果图的左侧边缘似乎并没有完全贴合我们的A4纸，只是由于我的峰值提取算法导致的，如果采用质心算法效果会更好。

Hough的强大之处是，如果我们把A4纸的一角盖住了，我们也能预测出A4纸原来的样子，仔细想想，是不是这样？

好了，今天给出的是满满的干货，要吸收还得多打打代码。我写得太急，代码一点都不优雅，大家千万别学我，最好把各个功能都结构化，方便日后重用。

最后贴出我的完整代码吧：

#include "CImg.h"
using namespace cimg_library;
#include <vector>
#include <iostream>

#define FILE "A4.bmp"

#define BLUR 3
#define GRADLIMIT 20
#define THRESHOLD 650 
#define DIFF 200
#define SLOPE_FLAG 1
#define PI 3.14159265358979323846

// 2D Dot(x, y) [value]  
struct Dot {
	int x, y, value;
	Dot(int _x, int _y, int _value)
		:x(_x), y(_y), value(_value){}
};

// 2D Straight Line(x, y) [y = m*x + b]
struct Line {
	double m, b;
	Line(double _m, double _b)
		:m(_m), b(_b){}
};

// Calculate the distance
double distance(double x, double y) {
	return sqrt(x*x + y*y);
}

// Polar coordinate intersection at x
const int CrossX(int theta, int distance, int x) {
	double angle = (double)theta*PI / 180.0;
	double m = -cos(angle) / sin(angle);
	double b = (double)distance / sin(angle);
	return m*x + b;
}

// Polar coordinate intersection at y
const int CrossY(int theta, int distance, int y) {
	double angle = (double)theta*PI / 180.0;
	double m = -cos(angle) / sin(angle);
	double b = (double)distance / sin(angle);
	return ((double)(y - b) / m);
}

int main {
	CImg<double> image(FILE);
	CImg<double> grayImg(image);
	CImg<double> result(image);

	// Gray scale;
	cimg_forXY(grayImg, x, y) {
		double R = grayImg(x, y, 0, 0);
		double G = grayImg(x, y, 0, 1);
		double B = grayImg(x, y, 0, 2);
		double Gray = (R * 299 + G * 587 + B * 114 + 500) / 1000;
		grayImg(x, y, 0, 0) = Gray;
		grayImg(x, y, 0, 1) = Gray;
		grayImg(x, y, 0, 2) = Gray;
	}

	// Blur
	grayImg.blur(BLUR);

	CImg<double> gradnum(image.width, image.height, 1, 1, 0);
	int maxDistance = distance(image.width, image.height);
	CImg<double> hough(360, maxDistance, 1, 1, 0);
	// 定义3*3领域矩阵I
	CImg_3x3(I, double);
	// 遍历计算梯度值
	cimg_for3x3(grayImg, x, y, 0, 0, I, double) {
		const double ix = Inc - Ipc;
		const double iy = Icp - Icn;
		double grad = std::sqrt(ix*ix + iy*iy);
		// 梯度大于阈值才赋值
		if (grad > GRADLIMIT) {
			gradnum(x, y) = grad;
			cimg_forX(hough, angle) {
				double rangle = (double)angle*PI / 180.0;
				int polar = (int)(x*cos(rangle) + y*sin(rangle));
				if (polar >= 0 && polar < hough.height) {
					hough(angle, polar) += 1;
				}
			}
		}
	}

	// Find peaks
	std::vector<Dot*> peaks;
	cimg_forXY(hough, angle, polar) {
		if (hough(angle, polar) > THRESHOLD) { // 是否是峰值
			bool flag = false;
			const int ymin = 0;
			const int ymax = image.height - 1;
			const int x0 = CrossY(angle, polar, ymin);
			const int x1 = CrossY(angle, polar, ymax);

			const int xmin = 0;
			const int xmax = image.width - 1;
			const int y0 = CrossX(angle, polar, xmin);
			const int y1 = CrossX(angle, polar, xmax);

			if (x0 >= 0 && x0 <= xmax ||				// 表示的直线是否在图像内
				x1 >= 0 && x1 <= xmax ||
				y0 >= 0 && y0 <= ymax ||
				y1 >= 0 && y1 <= ymax) {
				for (int i = 0; i < peaks.size; ++i) {    // 遍历数组，找相邻峰值
					if (distance(peaks[i]->x - angle, peaks[i]->y - polar) < DIFF) {   // 存在相邻峰值
						flag = true;
						if (peaks[i]->value < hough(angle, polar)) {	// 如果比相邻峰值还大
							peaks[i] = new Dot(angle, polar, hough(angle, polar));   // 替换为当前峰值
						}
					}
				}
				if (flag == false) {		// 当前峰值无相邻峰值
					peaks.push_back(new Dot(angle, polar, hough(angle, polar)));   // 加入新峰值
				}
			}
		}
	}
	std::cout << "lines: " << std::endl;

	// Transform polar coordinates to rectangular coordinates
	std::vector<Line*> lines;
	for (int i = 0; i < peaks.size; ++i) {
		double angle = (double)peaks[i]->x*PI / 180.0;
		double m = -cos(angle) / sin(angle);
		double b = (double)peaks[i]->y / sin(angle);
		lines.push_back(new Line(m, b));
		std::cout << "y = (" << m << ") x + (" << b << ")" << std::endl;
	}

	std::cout << std::endl << "intersections: " << std::endl;

	// Calculate line intersections
	std::vector<Dot*> intersections;
	for (int i = 0; i < lines.size; ++i) {
		for (int j = i + 1; j < lines.size; ++j) {
			double m0 = lines[i]->m;
			double m1 = lines[j]->m;
			double b0 = lines[i]->b;
			double b1 = lines[j]->b;
			double x = (b1 - b0) / (m0 - m1);
			double y = (m0*b1 - m1*b0) / (m0 - m1);
			if (x >= 0 && x < result.width && y >= 0 && y < result.height) {
				intersections.push_back(new Dot(x, y, 0));
				std::cout << "(" << x << ", " << y << ")" << std::endl;
			}
		}
	}

	std::cout << std::endl;

	// Draw lines
	for (int i = 0; i < lines.size; ++i) {
		const int ymin = 0;
		const int ymax = result.height - 1;
		const int x0 = (double)(ymin - lines[i]->b) / lines[i]->m;
		const int x1 = (double)(ymax - lines[i]->b) / lines[i]->m;

		const int xmin = 0;
		const int xmax = result.width - 1;
		const int y0 = xmin*lines[i]->m + lines[i]->b;
		const int y1 = xmax*lines[i]->m + lines[i]->b;

		const double color = { 255, 0, 255 };

		if (abs(lines[i]->m) > SLOPE_FLAG) {
			result.draw_line(x0, ymin, x1, ymax, color);
		}
		else {
			result.draw_line(xmin, y0, xmax, y1, color);
		}
	}

	// Draw intersections
	for (int i = 0; i < intersections.size; ++i) {
		const double color = { 255, 100, 255 };
		result.draw_circle(intersections[i]->x, intersections[i]->y, 50, color);
	}

	result.display;

	// Display
	//grayimg.display;
	//gradnum.display;
	//hough.display;
	result.display;

}