Harris角点检测（1）

最新推荐文章于 2022-08-08 14:34:10 发布

zizi7

最新推荐文章于 2022-08-08 14:34:10 发布

阅读量1.5k

点赞数 2

分类专栏：图像处理

本文链接：https://blog.csdn.net/zizi7/article/details/50112169

版权

图像处理专栏收录该内容

47 篇文章 17 订阅

订阅专栏

1. 理论

角点是一个局部概念。窗口在局部范围移动，如果采集到的灰度值没有显著变化，那么该区域不存在角点；如果在某一方向上有显著变化，那么该区域可能存在一条直线；如果在各个方向都有变化，那么该区域可能存在角点。

Harris以当前待检点为原点，调整x和y方向的偏移值获取相邻像素与该像素的相似度：，其中w(u,v)反映了每个像素的权重，Harris使用高斯分布对其赋值。

根据泰勒公式，有，其中I_x和I_y分别为x和y方向的偏导。这样，相似度公式可以整理为，

其中。显然，这是个二次项函数，而二次项函数从几何学上来说是一个椭圆函数，由M(x,y)决定的椭圆函数的两个特征值反映了其扁率和尺寸。Harris算法就是根据这两个特征值来确定当前点是否为角点的（图1-1）。

图1-1. 椭圆特征值与图像特征的关系

但是，计算矩阵特征值比较麻烦，在这里用矩阵M的行列式和迹来反映特征值的情况：，其中，，显然，当两个特征值较大时R取较大值。

至此，可以归纳出Harris检测角点的算法步骤：

（1）计算图像在X和Y方向上的梯度

（2）计算图像两个方向的点乘结果

（3）高斯函数对其进行加权，得到矩阵M

（4）计算每个像素的Harris响应值R，使用阈值t的置为0

（5）在3*3或5*5邻域内进行非最大值抑制（取窗口内最大值对应的点作为角点）

2. Harris角点的性质

（1）由于微分算子的特性（邻域灰度值分布特点），图像整体上的亮度、对比度变化并不会改变Harris角点的位置，但由于阈值 t 的关系，可能会影响Harris角点的检测数量

（2）由于二阶矩阵的特性（椭圆方向变化并不影响特征值本身），Harris角点具有旋转不变性

（3）因为窗口尺寸是一定的，因此Harris角点不具备尺度不变性（角点位置放大可能被检测为边缘或曲线）

3. 代码

#define PI 3.14159
#define MAXPOINTS 50000

typedef struct myPoint
{
	int x;
	int y;
}_point;

typedef struct mySize
{
	int height;
	int width;
}_size;

template<typename T1, typename T2>
void convolution(T1* img, T2** out, _size imgSize, T2* kernel, _size kerSize);
void dotMultiply(int* src1, int* src2, int** dst, _size srcSize);
void gaussian(int* src, float** out, _size imgSize, _size gaussSize, float sita);
void strengthR(float* A, float* B, float* C, float** R, float* maxR, _size imgSize, float alpha);
void depressR(float* srcR, _point* dstR, int* numDst, _size imgSize, _size depressWin, float threshold);
void _cornerHarris(unsigned char* src, _size srcSize, _point* corners, int* numCorner, float alpha);

void main()
{
	cv::Mat img = cv::imread("../file/butterfly.jpg", 0);
	unsigned char* imgBuffer = new unsigned char[img.rows*img.cols];
	for(int i=0; i<img.rows; i++)
	{
		uchar* ptr = img.ptr<uchar>(i);
		for(int j=0; j<img.cols; j++)
			imgBuffer[i*img.cols+j] = ptr[j];
	}

	_size srcSize;
	srcSize.height = img.rows;
	srcSize.width = img.cols;
	_point corners[MAXPOINTS];
	int numCorners = 0;
	_cornerHarris(imgBuffer, srcSize, corners, &numCorners, 0.04);

	cv::Mat colorImg;
	cv::cvtColor(img, colorImg, CV_GRAY2BGR);
	for(int i=0; i<numCorners; i++)
	{
		cv::Point center(corners[i].x, corners[i].y);
		cv::circle(colorImg, center, 2, cv::Scalar(0,0,255));
	}
	cv::namedWindow("show");
	cv::imshow("show", colorImg);
	cv::waitKey(0);
}

template<typename T1, typename T2>
void convolution(T1* img, T2** out, _size imgSize, T2* kernel, _size kerSize)
{
	for(int i=0; i<imgSize.height; i++)
	{
		for(int j=0; j<imgSize.width; j++)
		{
			int count = 0;
			T2 sumValue = 0;
			for(int m=i-kerSize.height; m<=i+kerSize.height; m++)
			{
				for(int n=j-kerSize.width; n<=j+kerSize.width; n++)
				{
					if(m>=0 && m<imgSize.height && n>=0 && n<imgSize.width)
						sumValue += T2(img[m*imgSize.width+n])*kernel[count];
					count++;
				}
			}
			(*out)[i*imgSize.width+j] = sumValue;
		}
	}
}

void dotMultiply(int* src1, int* src2, int** dst, _size srcSize)
{
	for(int i=0; i<srcSize.height; i++)
	{
		for(int j=0; j<srcSize.width; j++)
			(*dst)[i*srcSize.width+j] = src1[i*srcSize.width+j]*src2[i*srcSize.width+j];
	}
}

void gaussian(int* src, float** out, _size imgSize, _size gaussSize, float sita)
{
	int sizeKernel = (2*gaussSize.height+1) * (2*gaussSize.width+1);
	float* gaussKernel = new float[sizeKernel];
	for(int i=-gaussSize.height; i<=gaussSize.height; i++)
	{
		for(int j=-gaussSize.width; j<=gaussSize.width; j++)
		{
			float tmp = -1*(i*i+j*j)/(2*sita*sita);
			gaussKernel[(i+gaussSize.height)*(2*gaussSize.width+1)+(j+gaussSize.width)] =
				exp(tmp)/(2*PI*sita*sita);
		}
	}

	convolution(src, out, imgSize, gaussKernel, gaussSize);
}

void strengthR(float* A, float* B, float* C, float** R, float* maxR, _size imgSize, float alpha)
{
	*maxR = 0;
	for(int i=0; i<imgSize.height; i++)
	{
		for(int j=0; j<imgSize.width; j++)
		{
			float detM = A[i*imgSize.width+j]*B[i*imgSize.width+j] - C[i*imgSize.width+j]*C[i*imgSize.width+j];
			float traceM = A[i*imgSize.width+j] + B[i*imgSize.width+j];
			float result = detM - alpha*(traceM*traceM);
			if(result > *maxR)
				*maxR = result;
			(*R)[i*imgSize.width+j] = result;
		}
	}
}

void depressR(float* srcR, _point* dstR, int* numDst, _size imgSize, _size depressWin, float threshold)
{
	*numDst = 0;
	for(int i=0; i<imgSize.height; i++)
	{
		for(int j=0; j<imgSize.width; j++)
		{
			float flagValue = srcR[i*imgSize.width+j]<threshold?0:srcR[i*imgSize.width+j];
			int numPoint = 0, numPass = 0;
			for(int m=i-depressWin.height; m<=i+depressWin.height; m++)
			{
				for(int n=j-depressWin.width; n<=j+depressWin.width; n++)
				{
					if(m>=0 && m<imgSize.height && n>=0 && n<imgSize.width)
					{
						float compareValue = srcR[m*imgSize.width+n]<threshold?0:srcR[m*imgSize.width+n];
						if(flagValue > compareValue)
							numPass ++;
						numPoint ++;
					}
				}
			}
			if(numPoint == numPass+1)
			{
				_point corner;
				corner.x = j; 
				corner.y = i;
				dstR[(*numDst)++] = corner;
			}
		}
	}
}

void _cornerHarris(unsigned char* src, _size srcSize, _point* corners, int* numCorner, float alpha)
{
	int kernel[3] = {-1, 0, 1};
	_size sizeX, sizeY;
	sizeX.height = 0; sizeX.width = 1;
	sizeY.height = 1; sizeY.width = 0;
	int* gradientX = new int[srcSize.height*srcSize.width];
	int* gradientY = new int[srcSize.height*srcSize.width];
	convolution(src, &gradientX, srcSize, kernel, sizeX);
	convolution(src, &gradientY, srcSize, kernel, sizeY);

	int* gradientXX = new int[srcSize.height*srcSize.width];
	int* gradientYY = new int[srcSize.height*srcSize.width];
	int* gradientXY = new int[srcSize.height*srcSize.width];
	dotMultiply(gradientX, gradientX, &gradientXX, srcSize);
	dotMultiply(gradientY, gradientY, &gradientYY, srcSize);
	dotMultiply(gradientX, gradientY, &gradientXY, srcSize);

	delete[] gradientX;
	delete[] gradientY;

	_size gaussSize;
	gaussSize.height = 1;  gaussSize.width = 1;
	float sita = 1;
	float* M_A = new float[srcSize.height*srcSize.width];
	float* M_B = new float[srcSize.height*srcSize.width];
	float* M_C = new float[srcSize.height*srcSize.width];
	gaussian(gradientXX, &M_A, srcSize, gaussSize, sita);
	gaussian(gradientYY, &M_B, srcSize, gaussSize, sita);
	gaussian(gradientXY, &M_C, srcSize, gaussSize, sita);

	delete[] gradientXX;
	delete[] gradientYY;
	delete[] gradientXY;

	float maxR = 0;
	float* R = new float[srcSize.height*srcSize.width];
	strengthR(M_A, M_B, M_C, &R, &maxR, srcSize, alpha);

	delete[] M_A;
	delete[] M_B;
	delete[] M_C;

	float threshold = 0.1;
	threshold *= maxR;
	_size depressSize;
	depressSize.height = 1; depressSize.width = 1;
	depressR(R, corners, numCorner, srcSize, depressSize, threshold);

	delete[] R;
}

效果图（左：原图，中：旋转30度，右：放大2倍）