Harris Corner Detection Principle and C + + implementation __c++

1. First of all, we can't help but ask what is Harris corner point.

For corner points, there is no clear mathematical definition so far. But you can think of a corner point as the extreme point, that is, the point where the attribute is particularly prominent. General corner Detection is a specific definition, or can be specific detection of points of interest detection. This means that the point of interest can be a corner point, an isolated point with the strongest or smallest strength on some attributes, an endpoint of a line segment, or a point with the largest local curvature on the curve.

Generally speaking, in a pair of images, we can think that the corner point is the object contour line connection point (see Figure 1), when the shooting angle of view changes, these characteristic points still can maintain the stable property.

Figure 1 Corner

While preserving the important features of image graphics, corner points can effectively reduce the amount of information, so that the content of the information is very high, effectively improve the speed of the calculation, is conducive to the reliable image matching, making real-time processing possible. Its various applications, here I will not repeat.

2. How to detect the Harris corner point.

Fig. 2 The basic idea of corner point detection

The original idea of corner detection is to take a neighborhood window of a pixel, and when the window moves in a small range in all directions, observe the change in the average pixel gray value in the window (that is, E (u,v), window-averaged changes of intensity). From the figure above, we can roughly divide an image into three regions (' flat ', ' edge ', ' Corner '), and these three areas vary.

Wherein, the U,V is the window in the horizontal, vertical direction of the offset,

Here's a quick review of Taylor's progression, because it's going to be used right now.

This is a one-dimensional case, and for multivariate functions there is a similar Taylor formula.

In the two-dimensional Taylor series expansion of I (X+U,Y+V), we take the first order approximation and have

The part of the blue coil in the figure we call the structural tensor (structure tensor), expressed in M.

Here, let's be clear about what our purpose is. Our aim is to find such pixel points, it makes our u,v regardless of how the value, E (u,v) are changed relatively large, this pixel is what we are looking for the corner point. It is not difficult to observe that the E (u,v) after the approximate treatment is a two-second type, which is known by the following theorem.

E (U,V) = constant, we can use an ellipse to depict this function.

The length axis of an ellipse is the quantity corresponding to the two eigenvalues of the structure tensor m. By judging the situation we can distinguish between the ' flat ', ' edge ', ' Corner ' of the three areas, because the most intuitive impression:

Corner: In the horizontal, vertical two direction changes are larger points, that is, IX, iy are larger;
Edge: Only in the horizontal, or only in the vertical direction there are larger points, that is, IX and iy only one larger;
Flat: In the horizontal, vertical direction of the changes are small points, that is, IX, iy are smaller;

and the structure tensor m is composed of Ix,iy, its eigenvalues can reflect the situation of ix,iy, and I will explain the physical meaning of the ellipse in a more comprehensible way.

Is this clearer? ^_^ ... so we can draw the conclusion that

Of course, Daniel did not stop here, so it is not convenient to judge the corner point by judging the value of two variables. So they came up with a better way, yes, is to define the corner response function R (Corner response functions),

What is the value of R for three different regions?

At this point, we can judge the value of R to determine whether a dot is a corner.

Corner point: R is a large numeric integer

Edge: R is a large numeric negative number

Flat area: Absolute value R is a decimal value

3. Harris Corner Point Detection algorithm steps

1. Using soble to calculate the gradient value of the XY direction

2. Calculate the Ix^2,iy^2,ix*iy

3. Using Gaussian function to filter the Ix^2,iy^2,ix*iy

4. Eigenvalues and response functions of the computed matrix M of local characteristic results (c (i,j) =det (m)-K (Trace (m)) ^2 (0.04<=k<=0.06)

5. The value C of the response function is calculated to suppress the maximum value, filter out some points that are not corner points, and satisfy the threshold value greater than the set.

Put the source code below:

#include "opencv2/imgproc/imgproc.hpp" #include "opencv2/highgui/highgui.hpp" #include <iostream> #include
;cmath> using namespace CV;
using namespace Std; /* RGB conversion to grayscale images is a common formula: Gray = r*0.299 + g*0.587 + b*0.114///****************** grayscale conversion function *************************//First parameter
Image input a reference to a color RGB image;//The second parameter imagegray is a reference to the grayscale image output after the conversion;//*******************************************************

void Convertrgb2gray (const Mat &image, Mat &imagegray); Sobel convolution factor calculates x, y direction gradient and gradient direction angle ********************//First parameter IMAGESOURC original grayscale image; The second parameter imagesobelx is the x direction gradient image;//The third parameter imagesobely is the y-directional gradient image;//The fourth parameter pointdrection is the gradient direction angle array pointer//************************** void Sobelgraddirction (Mat &imagesource, Mat &imagesobelx, Mat &

imagesobely); Calculates the sobel of the X-directional gradient amplitude of the square *************************//First parameter IMAGEGRADX is the X-directional gradient image; The second parameter sobelampxx is the output of the X-directional gradient image squared//************************************************void sobelxx (const Mat IMAGEGRADX, mat_<float> &sobelampxx); Calculate the y-direction gradient amplitude of the Sobel squared *************************//The first parameter imagegrady is the y-directional gradient image; The second parameter sobelampxx is the output of the Y-direction gradient image of the square//************************************************************* void Sobelyy (const

Mat Imagegrady, mat_<float> &sobelampyy); The product *************************//The first parameter of the XY-gradient amplitude of the Sobel is computed IMAGEGRADX is the X-directional gradient image; The second parameter imagegrady is the y-directional gradient image;//The second parameter sobelampxy is the output xy direction gradient image//*******************************************************

void Sobelxy (const Mat IMAGEGRADX, const Mat imagegrady, mat_<float> &sobelampxy); Calculate the weights of a Gaussian array *****************//The first parameter size is the size of the side length of the convolution nucleus represented by//The second parameter sigma represents the size of the Sigma//**************

Double *getoneguassionarray (int size, double sigma);
Realization of Gaussian filter function *****************//The first parameter srcimage is represented by the input of the original//The second parameter DST represents the output Graph///The third parameter size represents the size of the side length of the convolution core. //****void Mygaussianblur (mat_<float> &srcimage, mat_<

float> &dst, int size); Compute the eigenvalues of the local special-rise result matrix M and the response function H = (a*b-c)-k* (a+b) ^2******//m//a C//c B//tr (m) =a+b=a+b//det (m) =a*b=a*b-c^2//COMPUTE Output Response Letter The worthy matrix of the number//**************************************************************************** void HarrisResponse (Mat_ <float> &gaussxx, mat_<float> &gaussyy, mat_<float> &gaussxy, Mat_<float> &


Resultdata,float k); Non-maximal suppression and satisfying thresholds and local maxima within a neighborhood are angular points **************//The first parameter is the matrix of the response function///The second parameter is the input grayscale image//The third parameter is the output of the corner point detected by the result graph void

Localmaxvalue (mat_<float> &resultdata, Mat &srcgray, Mat &resultimage,int);
	int main () {const Mat srcimage = Imread ("3.jpg");
		if (!srcimage.data) {printf ("Could not load image...\n");
	return-1;
	} imshow ("Srcimage", srcimage);
	Mat Srcgray;
	Convertrgb2gray (Srcimage, Srcgray);
	Mat Imagesobelx;
	Mat imagesobely;
	Mat Resultimage; MaT_<float> IMAGESOBELXX;
	Mat_<float> Imagesobelyy;
	Mat_<float> Imagesobelxy;
	Mat_<float> gaussianxx;
	Mat_<float> Gaussianyy;
	Mat_<float> Gaussianxy;
	Mat_<float> Harrisrespond;
	Calculate soble xy Gradient sobelgraddirction (Srcgray, Imagesobelx, imagesobely);
	Calculates the square sobelxx of the gradient in the x direction (Imagesobelx, IMAGESOBELXX);
	Sobelyy (imagesobely, imagesobelyy);
	Sobelxy (Imagesobelx, imagesobely, Imagesobelxy);
	Compute gauss Blur xx YY XY mygaussianblur (IMAGESOBELXX, gaussianxx,3);
	Mygaussianblur (Imagesobelyy, Gaussianyy, 3);
	Mygaussianblur (Imagesobelxy, Gaussianxy, 3);
	Harrisresponse (Gaussianxx, Gaussianyy, Gaussianxy, Harrisrespond, 0.05);
	Localmaxvalue (Harrisrespond, Srcgray, Resultimage, 3);
	Imshow ("Imagesobelx", Imagesobelx);
	Imshow ("imagesobely", imagesobely);
	Imshow ("Resultimage", resultimage);
	Waitkey (0);
return 0; 
	} void Convertrgb2gray (const Mat &image, Mat &imagegray) {if (!image.data | | image.channels ()!= 3) {return; //Create a single channel grayscale image imagEgray = Mat::zeros (Image.size (), CV_8UC1);
	Remove the pointer to the array that stores the pixel of the image uchar *pointimage = Image.data;
	Uchar *pointimagegray = Imagegray.data;
	The number of bytes taken out of the image per line size_t stepimage = image.step;
	size_t Stepimagegray = Imagegray.step; for (int i = 0; i < imagegray.rows. i++) {for (int j = 0; J < Imagegray.cols; J +) {Pointimagegray[i*stepi Magegray + j] = (Uchar) (0.114*pointimage[i*stepimage + 3 * j] + 0.587*pointimage[i*stepimage + 3 * j + 1] + 0.299*pointima
		Ge[i*stepimage + 3 * j + 2]); }}///storage gradient membrane long void sobelgraddirction (Mat &imagesource, Mat &imagesobelx, Mat &imagesobely) {Imagesobelx
	= Mat::zeros (Imagesource.size (), CV_32SC1);
	imagesobely = Mat::zeros (Imagesource.size (), CV_32SC1);
	Remove the first address of an array of original and x and y gradient graphs uchar *p = imagesource.data;
	Uchar *px = Imagesobelx.data;

	Uchar *py = Imagesobely.data;
	Remove the number of bytes occupied by each row int step = Imagesource.step;

	int stepxy = Imagesobelx.step; int index = 0;//gradient Orientation angle index for (int i = 1; i < imagesource.rows-1; ++i) {foR (Int j = 1; j < imagesource.cols-1 ++j) {//traversal of each pixel on the image by pointer double GradY = p[(i + 1) *step + j-1] + p[( i + 1) *step + j] * 2 + p[(i + 1) *step + j + 1]-p[(i-1) *step + j-1]-p[(i-1) *step + j] * 2-p[(i-1) *step + j +
			1];

			Py[i*stepxy + j* (stepxy/step)] = ABS (GradY); Double GRADX = p[(i-1) *step + j + 1] + P[i*step + j + 1] * 2 + p[(i + 1) *step + j + 1]-p[(i-1) *step + j-1]-p[i*s
			TEP + j-1] * 2-p[(i + 1) *step + j-1];
		Px[i*stepxy + j* (stepxy/step)] = ABS (GRADX);
	}///convert gradient array to 8-bit unsigned integer convertscaleabs (imagesobelx, IMAGESOBELX);
Convertscaleabs (imagesobely, imagesobely); } void Sobelxx (const Mat IMAGEGRADX, mat_<float> &sobelampxx) {sobelampxx = mat_<float> (imageGradX.siz
	E (), CV_32FC1); for (int i = 0; i < sobelampxx.rows. i++) {for (int j = 0; J < Sobelampxx.cols; J +) {Sobelampxx.at<flo
		At> (i, j) = Imagegradx.at<uchar> (i, J) *imagegradx.at<uchar> (I, j); }}//convertscaleabs(Sobelampxx, SOBELAMPXX); } void Sobelyy (const Mat imagegrady, mat_<float> &sobelampyy) {sobelampyy = mat_<float> (imageGradY.size
	(), CV_32FC1); for (int i = 0; i < sobelampyy.rows. i++) {for (int j = 0; J < Sobelampyy.cols; J +) {Sobelampyy.at<flo
		At> (i, j) = Imagegrady.at<uchar> (i, J) *imagegrady.at<uchar> (I, j);
}//convertscaleabs (Sobelampyy, SOBELAMPYY); } void Sobelxy (const Mat IMAGEGRADX, const Mat imagegrady, mat_<float> &sobelampxy) {sobelampxy = Mat_<flo
	At> (Imagegradx.size (), CV_32FC1); for (int i = 0; i < sobelampxy.rows. i++) {for (int j = 0; J < Sobelampxy.cols; J +) {Sobelampxy.at<flo
		At> (i, j) = Imagegradx.at<uchar> (i, J) *imagegrady.at<uchar> (I, j);
}//convertscaleabs (Sobelampxy, SOBELAMPXY);
	///computes a Gaussian array of weights double *getoneguassionarray (int size, double sigma) {double = 0.0;

	Define Gauss kernel radius int kerR = SIZE/2; Create a dynamic one-dimensional array of size sizes double *arr = newDouble[size]; for (int i = 0; i < size; i++) {//Gaussian function constants can not be computed, will be eliminated in the process of normalization arr[i] = exp (-(I-KERR) * (I-kerr))/(2 * sigm
		A*sigma));
		Sum + + arr[i];//All values are added}//normalized for (int i = 0; i < size; i++) {arr[i]/= sum;
	cout << Arr[i] << Endl;
return arr; } void Mygaussianblur (mat_<float> &srcimage, mat_<float> &dst, int size) {Cv_assert (SrcImage.chann Els () = = 1 | | Srcimage.channels () = = 3);
	Only process single channel or three channel image int kerR = SIZE/2;
	DST = Srcimage.clone ();
	int channels = Dst.channels ();
	double* arr;
	arr = Getoneguassionarray (size, 1); First, the convolution for the horizontal direction of the Gaussian array//traversal image is obtained (int i = KerR; i < Dst.rows-kerr; i++)
			{for (int j = KerR; J < Dst.cols-kerr; J + +) {float Guassionsum[3] = {0};  The sliding window search completes the Gaussian kernel smoothing for (int k =-kerr k <= KerR; k++) {if (channels = 1)//If only single channel {guassionsum[0] = Arr[kerr + K] * dst.at<float> (i, j + k);/rows unchanged, column transforms, first horizontal convolution} else if (ChannEls = = 3)//If it is a three-channel situation {vec3f BGR = dst.at<vec3f> (i, j + K);
					Auto A = Arr[kerr + K];
					Guassionsum[0] + = a*bgr[0];
					GUASSIONSUM[1] + = a*bgr[1];
				GUASSIONSUM[2] + = a*bgr[2];
				for (int k = 0; k < channels; k++) {if (guassionsum[k) < 0) Guassionsum[k] = 0;
			else if (Guassionsum[k] > 255) guassionsum[k] = 255;
			if (channels = = 1) dst.at<float> (i, j) = Static_cast<float> (guassionsum[0)); else if (channels = 3) {vec3f BGR = {static_cast<float> (guassionsum[0)), static_cast<float> (guassion
				SUM[1]), static_cast<float> (guassionsum[2])};
			Dst.at<vec3f> (i, j) = BGR;
		}}//vertical for (int i = Kerr; i < Dst.rows-kerr; i++) {for (int j = Kerr; J < Dst.cols-kerr; J +)
			{float Guassionsum[3] = {0};  The sliding window search completes the Gaussian kernel smoothing for (int k =-kerr k <= KerR; k++) {if (channels = 1)//If only single channel {guassionsum[0] = Arr[kerr + K] *Dst.at<float> (i + K, j),//row change, column does not change, then do the vertical convolution} else if (channels = 3)//If it is a three-channel situation {vec3f BGR = dst.a
					T<vec3f> (i + K, j);
					Auto A = Arr[kerr + K];
					Guassionsum[0] + = a*bgr[0];
					GUASSIONSUM[1] + = a*bgr[1];
				GUASSIONSUM[2] + = a*bgr[2];
				for (int k = 0; k < channels; k++) {if (guassionsum[k) < 0) Guassionsum[k] = 0;
			else if (Guassionsum[k] > 255) guassionsum[k] = 255;
			if (channels = = 1) dst.at<float> (i, j) = Static_cast<float> (guassionsum[0)); else if (channels = 3) {vec3f BGR = {static_cast<float> (guassionsum[0)), static_cast<float> (guassion
				SUM[1]), static_cast<float> (guassionsum[2])};
			Dst.at<vec3f> (i, j) = BGR;
}} delete[] arr; } void Harrisresponse (mat_<float> &gaussxx, mat_<float> &gaussyy, mat_<float> &GaussXY, Mat_<float> &resultdata,float k) {//Create a matrix of response function output Resultd