实现Harris Corner Detector，输出结果以及中间过程。

实验过程

实现了一个HarrisCornerDetector函数，函数原型如下：

1
void HarrisCornerDetector(Mat& src, Mat& R, int aperture_size, double k);

整体实现过程按照HarrisCornerDetector的运算过程来。

先将彩色图片转为单通道的灰度图，便于计算。

1
cvtColor(src, src_gray, cv::COLOR_BGR2GRAY);

然后计算X和Y方向的导数，本质上是用一个算子做一下卷积，我这里使用了Sobel算子，会根据aperture_size生成对应的模板来计算近似的导数。

1
2
Sobel(src_gray, Ix, CV_32FC1, 1, 0, aperture_size);
Sobel(src_gray, Iy, CV_32FC1, 0, 1, aperture_size);

然后就对整张导数的图来计算IxIx, IxIy, IyIy:

1
2
3
4
5
6
7
8
for (int i = 0; i < size.height; ++i) {
    for (int j = 0; j < size.width; ++j) {
        IxIx.at<float>(i,j) = Ix.at<float>(i, j) * Ix.at<float>(i, j);
        IxIy.at<float>(i,j) = Ix.at<float>(i, j) * Iy.at<float>(i, j);
        IyIy.at<float>(i,j) = Iy.at<float>(i, j) * Iy.at<float>(i, j);
    }
}

如此就得到了这个矩阵的值：

$$
\begin{vmatrix}
I_xI_x & I_xI_y \
I_xI_y & I_yI_y
\end{vmatrix}
$$

然后就是用W[x,y] 来求和

$$
\sum W(x,y) \begin{vmatrix}
I_xI_x & I_xI_y \
I_xI_y & I_yI_y
\end{vmatrix}
$$
W(x,y)可以取高斯滤波函数。

1
2
3
4
Size block(3,3);
GaussianBlur(IxIx, IxIx, block, 0);
GaussianBlur(IxIy, IxIy, block, 0);
GaussiaBlur(IyIy, IyIy, block, 0);

然后就可以求得到的矩阵的特征值。

因为是2x2的矩阵，特征方程就是 λ^2-(a+d)λ+ad-bc=0, 直接使用求根公式来求特征值，用韦达定理可以得到R的值。
$$
R = det(H) - k \times trace(H) ^2 \
$$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
for (int i = 0; i < size.height; ++i) {
		for (int j = 0; j < size.width; ++j) {
			float a = IxIx.at<float>(i, j);
			float b = IxIy.at<float>(i, j);
			float c = b;
			float d = IyIy.at<float>(i, j);
			// 2-D mat a b c d
			// λ^2-(a+d)λ+ad-bc=0
			// λ1 + λ2 = a+d
			// λ1 * λ2 = ad-bc
			R.at<float>(i,j) = (a*d - b*c) - k*(a + d)*(a + d);
			largeEigen.at<float>(i, j) = ((a + d) + sqrt((a + d)*(a + d) - 4 * (a*d - b*c))) / 2;
			smallEigen.at<float>(i, j) = ((a + d) - sqrt((a + d)*(a + d) - 4 * (a*d - b*c))) / 2;
		}
	}

然后就可以根据R矩阵的值来画出检测出是角的点。

画的时候为避免多个点聚集，使用了Non Maximum Suppression, 只取一定范围内R值最大的点作为角的特征点。

1
2
3
4
5
6
7
8
9
10
for (int i = 0; i < size.height; ++i) {
    for (int j = 0; j < size.width; ++j) {
        if ((int)R.at<uchar>(i, j) > threshold) {
            // Non Maximum Suppression
            if (R.at<uchar>(i, j) == maxValue(R, NMS_size, i, j)) {
              	circle(src, Point(j, i), 5, Scalar(0, 0, 255.0), 2, 8, 0);
            }
        }
    }
}

最后展示结果并存储结果

1
2
3
4
5
6
7
8
9
10
11
imshow("LargeEigen", largeEigen);
imshow("SmallEigen", smallEigen);
imshow("R", R);
imshow("result", src);

waitKey(0);

imwrite("LargeEigen.png", largeEigen);
imwrite("SmallEigen.png", smallEigen);
imwrite("R.png", R);
imwrite("result.png", src);

实验结果

略。

心得体会

实验中实现了 Harris Corner Detector，充分体会了这样的角点检测的运算过程。在实现的过程中，深刻理解了Harris的这种检测的原理，推导了计算的公式，加深了理解。实验的结果也符合预期。

在代码的编写过程中也进一步熟悉了OpenCV的使用，可以熟练地使用OpenCV进行图片的处理和卷积运算，提高了使用的熟练度。

附：源代码

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
#include <opencv2/opencv.hpp>
#include <iostream>
#include <string>
#include <cmath>
using namespace std;
using namespace cv;

void HarrisCornerDetector(Mat& src, Mat& R, int aperture_size, double k);

int main(int argc, const char** argv) {
	char* filename = new char[100];
	//system("dir");
	int apertureSize = 3;
	double k = 0.04;

	if (argc != 4) {
		cout << "Illegal Input." << endl;
		cout << "HarrisDetector.exe $path $k $aperture_size." << endl;
		//cout << "using default settings." << endl;
		//filename = "sample.png";
		filename = "Sydney_Opera_House_Sails_edit02_adj.jpg";
	}
	else {
		 // parse commandline args
		sscanf_s(argv[1], "%s", filename, 99);
		sscanf_s(argv[2], "%lf", &k);
		sscanf_s(argv[3], "%d", &apertureSize);
	}

	Mat src = imread(filename);
	if (!src.data) {
		cout << "imread failed" << endl;
		return 0;
	}
	Mat dst;

	cout << "path: " << filename << endl;
	cout << "k: " << k << endl;
	cout << "aperture_size: " << apertureSize << endl;

	HarrisCornerDetector(src, dst, apertureSize, k);
	
	return 0;
}

int maxValue(Mat& img, int size, int y, int x)
{
	uchar maxval = 0;
	for (int i = y-size; i<y + size; ++i)
	{
		if (i < 0 || i >= img.rows)continue;
		for (int j = x-size; j< x + size; ++j)
		{
			if (j<0 || j >= img.cols)continue;
			if (img.at<uchar>(i,j) > maxval)
			{
				maxval = img.at<uchar>(i, j);
			}
		}
	}
	//cout << "maxval" << (int)maxval << endl;
	return maxval;
}

void HarrisCornerDetector(Mat& src, Mat& R, int aperture_size, double k)
{
	Mat src_gray;
	// convert to gray
	cvtColor(src, src_gray, cv::COLOR_BGR2GRAY);
	// normalize src
	normalize(src_gray, src_gray, 0, 255, NORM_MINMAX);
	convertScaleAbs(src_gray, src_gray);

	R.create(src_gray.size(), CV_32FC1);
	Mat Ix, Iy;

	//sobel operation get Ix, Iy 
	Sobel(src_gray, Ix, CV_32FC1, 1, 0, aperture_size);
	Sobel(src_gray, Iy, CV_32FC1, 0, 1, aperture_size);
	//cout << Ix.size() << " " << Ix.channels() << " " << Ix.depth() << endl;

	// prepare Mat to store info
	Mat IxIx(src_gray.size(), CV_32FC1);
	Mat IxIy(src_gray.size(), CV_32FC1);
	Mat IyIy(src_gray.size(), CV_32FC1);
	Mat largeEigen(src_gray.size(), CV_32FC1);
	Mat smallEigen(src_gray.size(), CV_32FC1);
	//Mat heat(src_gray.size(), CV_8SC3);

	Size size = src_gray.size();
	for (int i = 0; i < size.height; ++i) {
		for (int j = 0; j < size.width; ++j) {
			IxIx.at<float>(i,j) = Ix.at<float>(i, j) * Ix.at<float>(i, j);
			IxIy.at<float>(i,j) = Ix.at<float>(i, j) * Iy.at<float>(i, j);
			IyIy.at<float>(i,j) = Iy.at<float>(i, j) * Iy.at<float>(i, j);
		}
	}

	// W[x,y] * I
	// W[x,y] is Gaussian Filter
	Size block(3,3);
	GaussianBlur(IxIx, IxIx, block, 0);
	GaussianBlur(IxIy, IxIy, block, 0);
	GaussianBlur(IyIy, IyIy, block, 0);

	for (int i = 0; i < size.height; ++i) {
		for (int j = 0; j < size.width; ++j) {
			float a = IxIx.at<float>(i, j);
			float b = IxIy.at<float>(i, j);
			float c = b;
			float d = IyIy.at<float>(i, j);
			// 2-D mat a b c d
			// λ^2-(a+d)λ+ad-bc=0
			// λ1 + λ2 = a+d
			// λ1 * λ2 = ad-bc
			R.at<float>(i,j) = (a*d - b*c) - k*(a + d)*(a + d);
			largeEigen.at<float>(i, j) = ((a + d) + sqrt((a + d)*(a + d) - 4 * (a*d - b*c))) / 2;
			smallEigen.at<float>(i, j) = ((a + d) - sqrt((a + d)*(a + d) - 4 * (a*d - b*c))) / 2;
		}
	}
	cout << endl;
	normalize(R, R, 0, 255, NORM_MINMAX, CV_32FC1, Mat());
	convertScaleAbs(R, R);
	int threshold = 50;
	int NMS_size = 15;
	// Draw circles with NMS
	for (int i = 0; i < size.height; ++i) {
		for (int j = 0; j < size.width; ++j) {
			if ((int)R.at<uchar>(i, j) > threshold) {
				// Non Maximum Suppression
				if (R.at<uchar>(i, j) == maxValue(R, NMS_size, i, j)) {
					circle(src, Point(j, i), 5, Scalar(0, 0, 255.0), 2, 8, 0);
				}
			}
		}
	}

	// Show the result
	imshow("LargeEigen", largeEigen);
	imshow("SmallEigen", smallEigen);
	//imshow("heat", heat);
	imshow("R", R);
	imshow("result", src);

	waitKey(0);

	imwrite("LargeEigen.png", largeEigen);
	imwrite("SmallEigen.png", smallEigen);
	imwrite("R.png", R);
	imwrite("result.png", src);
	destroyAllWindows();
}