opencv学习_12 (harris角点检测)

素颜马尾好姑娘i 2022-09-18 06:47 260阅读 0赞

原理：

Harris角点检测最直观的解释是:在任意两个相互垂直的方向上，都有较大变化的点。\---harris在A combined corner and edge detector这篇文章中提出来的。

![SouthEast][]

在moravec角点检测中，w(x,y)的取值是二元的，在窗口内部就取值为1，在窗口外部就取值为0，在harris的角点检测中，使用的是高斯窗口，所以w(x,y)表示的是高斯窗口中的权重。此时 当u和v取两组相互垂直的值时，E(u,v)都有较大值的点。

公式推导：

![SouthEast 1][]

根据二阶taylor展开式得出：

![SouthEast 2][]

![SouthEast 3][]

![SouthEast 4][]

![SouthEast 5][]

如何考虑E(u,v)在两个相互垂直的方向上都取得最大值呢？我们可以先求出E(u,v)在一个方向上的最大值为r1，即E(u1,v1) = r1,然后求其垂直方向上的值r2，我们意外的发现这里的r1和r2分别对应M的两个特征值，对应的向量(u1,v1),(u2,v2)为M的特征向量，此处的M被称为Hessian矩阵。

判断是否为角点的标准:

r1,r2都很小,对应于图像中的平滑区域

r1,r2都很大，对应于图像中的角点

r1，r2一个很大，一个很小，对应于图像中的边缘

Harris角点量计算公式: R = det(M) – a\*trace(M)^2,其中a为经验值，取值范围介于\[0.04, 0.06\],另一种角点量计算方法:R = det(M)/trace(M)—这个公式自由发挥，只要能反应角点的特征就可以了

其中det(M) = r1\*r2 = AB –C^2

Trace(M) = r1+r2 =A+B （其中A=Ix^2;; B = Iy^2;;; C=IxIy）

角点检测的步骤

根据上述的讨论，harris角点检测的步骤可以总结为:

<1>计算图像I(x,y)在x和y两个方向的梯度Ix,Iy：

![SouthEast 6][]

<2>计算梯度方向的乘积

![SouthEast 7][]

<3>使用高斯核对Ix^2,Iy^2,IxIy进行加权，计算矩阵M的元素A，B，C

![SouthEast 8][]

<4>根据计算角点量，并设定阀值，当角点量小于阀值，不是候选角点

<5>进行局部极大值抑制

代码：

#include <iostream>
    #include "cv.h"
    #include "cxcore.h"
    #include "highgui.h"
    using namespace std;
    
    
    
    //harris 角点检测需要的参数
    typedef struct HARRISPARAMS
    {
    	int gaussSize; //高斯窗口
    	float gaussSigma; //高斯方差
    	double threshold; //对角点量设定的阈值
    	int maximumSize; //局部极大值抑制时的窗口大小
    	
    }HARRISPARAMS,*PHARRISPARAMS;
    
    
    /*******************************
    *对源图像进行卷积运算
    *输入项
    *srcFloat    源图像
    *Ix          卷积的结果
    *dxTemplate  卷积模板
    *widthTemplate 模板的宽度
    *heightTemplate 模板的高度
    ********************************/
    void convolution(IplImage* srcFloat,IplImage* Ix,double* dxTemplate , int widthTemplate,int heightTemplate);
    
    
    /***************************
    *harris角点检测函数
    *输入项
    *srcIn   源图像
    *params  harris角点检测需要的参数
    *corners 存放harris角点坐标
    ****************************/
    void getHarrisPoints(IplImage* srcIn,PHARRISPARAMS params,CvSeq* corners);
    
    //主函数
    int main(int argc, char* argv[])
    {
    
    	//相关变量
    	IplImage* src,*srcGray;
    	CvMemStorage* mem = cvCreateMemStorage(0);
    	CvSeq* harrisPoints;
    	HARRISPARAMS harrisParams;
    
    	src = cvLoadImage("E:\\study_opencv_video\\lesson17_2\\2.jpg");//源图像
    	
    	srcGray = cvCreateImage(cvGetSize(src),8,1);        // 灰度图像
    	
    	if(!src)
    	{
    		cout << " src is null";
    		return 0;
    	}
    	
    	cvCvtColor(src,srcGray,CV_BGR2GRAY);
    
    	//harris角点保存的空间  角点坐标保存在一个序列中
    	harrisPoints = cvCreateSeq(0,sizeof(CvSeq),sizeof(CvPoint),mem);
    	
    	//设置相关参数
    	harrisParams.gaussSize = 5;   // 高斯窗口的大小
    	harrisParams.gaussSigma = 0.8;
    	harrisParams.threshold = 1000;
    	harrisParams.maximumSize = 21;
    
    	//进行harris角点检测
    	getHarrisPoints(srcGray,&harrisParams,harrisPoints);
    	
    	//获取每一个角点的坐标
    	for(int x=0;x<harrisPoints->total;x++)
    	{
    		//获取第x个角点的坐标
    		CvPoint* pt = (CvPoint*)cvGetSeqElem(harrisPoints,x);
    
    		//以角点坐标为中心  绘制一个半径为5的圆
    		cvCircle(src,*pt,2,cvScalar(255,0,255,0));
    	}
    
    	//cvSaveImage("dst.jpg",src);
    
    	//显示图像
    	cvNamedWindow("dst");
    	cvShowImage("dst",src);
    	cvWaitKey(0);
    
    	//释放资源
    	cvReleaseImage(&src);
    	cvReleaseImage(&srcGray);
    	cvReleaseMemStorage(&mem);
    	return 0;
    }
    
    /***************************
    *harris角点检测函数
    *输入项
    *srcIn   源图像
    *params  harris角点检测需要的参数
    *corners 存放harris角点坐标
    ****************************/
    void getHarrisPoints(IplImage* srcIn,PHARRISPARAMS params,CvSeq* corners)
    {
    	int x,y;
    
    	IplImage* srcFloat;
    	
    	srcFloat = cvCreateImage(cvGetSize(srcIn),32,1);
    	
    	cvConvertScale(srcIn,srcFloat);
    	
    	IplImage *Ix,*Iy,*IxIx,*IyIy,*IxIy,*A,*B,*C,*cornerness;
    	double *gaussWindow = new double[sizeof(double)*params->gaussSize*params->gaussSize];
    	//水平方向差分算子并求Ix
         double dxTemplate[9]={-1,0,1,
    		                   -1,0,1,
    						   -1,0,1};
    
    	//垂直方向差分算子并求Iy
         double dyTemplate[9]={-1,-1,-1,
    		                    0, 0, 0,
    							1, 1, 1};
    	
    	 //此处内存用得有点奢侈 请大家自行修改 节省内存
    	 //下面变量的含义与 第二十四集harris角点检测数学原理 的文档资料中定义的变量含义是一样的  请参阅
    	Ix = cvCreateImage(cvGetSize(srcFloat),32,1);
    	Iy = cvCreateImage(cvGetSize(srcFloat),32,1);
    	IxIx = cvCreateImage(cvGetSize(srcFloat),32,1);
    	IyIy = cvCreateImage(cvGetSize(srcFloat),32,1);
    	IxIy = cvCreateImage(cvGetSize(srcFloat),32,1);
    	A = cvCreateImage(cvGetSize(srcFloat),32,1);
    	B = cvCreateImage(cvGetSize(srcFloat),32,1);
    	C = cvCreateImage(cvGetSize(srcFloat),32,1);
    	cornerness = cvCreateImage(cvGetSize(srcFloat),32,1); //保存角点量
    
        convolution(srcFloat,Ix,dxTemplate,3,3); //计算Ix
    	convolution(srcFloat,Iy,dyTemplate,3,3); //计算Iy 
    
    	//计算Ix2、Iy2、IxIy
    	for(y=0;y<srcFloat->height;y++)
    	{
    		for(x=0;x<srcFloat->width;x++)
    		{
    			float IxValue,IyValue;
    			
    			IxValue = cvGetReal2D(Ix,y,x);
    			IyValue = cvGetReal2D(Iy,y,x);
    
    			cvSetReal2D(IxIx,y,x,IxValue*IxValue);
    			cvSetReal2D(IyIy,y,x,IyValue*IyValue);
    			cvSetReal2D(IxIy,y,x,IxValue*IyValue);
    
    			
    		}
    	}
    
    	//计算高斯窗口
    	for( y=0;y<params->gaussSize;y++)
    	{
    		for( x=0;x<params->gaussSize;x++)
    		{
    			
    			float dis,weight;
    			
    			dis = (y-params->gaussSize/2)*(y-params->gaussSize/2)+(x-params->gaussSize/2)*(x-params->gaussSize/2);
    			weight = exp(-dis/(2.0*params->gaussSigma));
    			*(gaussWindow+y*params->gaussSize+x)=weight;		
    			
    		}
    	}
    	 
    	convolution(IxIx,A,gaussWindow,params->gaussSize,params->gaussSize);//计算IxIx与高斯的卷积
    	convolution(IyIy,B,gaussWindow,params->gaussSize,params->gaussSize);//计算IyIy与高斯的卷积
    	convolution(IxIy,C,gaussWindow,params->gaussSize,params->gaussSize);//计算IxIy与高斯的卷积
    
    	//计算角点量
    	for(y=0;y<srcFloat->height;y++)
    	{
    		for(x=0;x<srcFloat->width;x++)
    		{
    			double cornernessValue,Avalue,Bvalue,Cvalue;
    			Avalue = cvGetReal2D(A,y,x);
    			Bvalue = cvGetReal2D(B,y,x);
    			Cvalue = cvGetReal2D(C,y,x);
    			
    			cornernessValue = (Avalue*Bvalue-Cvalue*Cvalue)/(Avalue+Bvalue+0.0000001);
    			
    			cvSetReal2D(cornerness,y,x,fabs(cornernessValue));	
    		}
    	
    	}
    
    	//计算局部极大值 及 极大值是否大于阈值
    	int beginY,endY,beginX,endX;
    	int halfWinSize = params->maximumSize/2;
    
    	beginY = halfWinSize;
    	endY = cornerness->height - halfWinSize;
    
    	beginX = halfWinSize;
    	endX = cornerness->width - halfWinSize;
    	
    	for(y=beginY;y<endY;)
    	{
    		for(x=beginX;x<endX;)
    		{
    			//寻找局部极大值 及其位置信息
    			float maxValue=0;
    			int flag = 0 ;
    			CvPoint maxLoc;
    			maxLoc.x = -1;
    			maxLoc.y = -1;
    
    			//首先计算以点(x,y)位中心的maximumSize*maximumSize的窗口内部的局部极大值
    			for(int winy=-halfWinSize;winy<=halfWinSize;winy++)
    			{
    				for(int winx=-halfWinSize;winx<=halfWinSize;winx++)
    				{
    					float value ;
    					value = cvGetReal2D(cornerness,y+winy,x+winx);
    					
    					//计算该窗口内 最大值 保存到max 并保存其坐标到maxLoc
    					if(value>maxValue)
    					{
    						maxValue = value;
    						maxLoc.x = x+winx;
    						maxLoc.y = y+winy;
    						flag = 1;
    					}
    				}
    			}
    
    			
    			//如果找到局部极大值 并且该值大于预先设定的阈值 则认为是角点
    			if(flag==1 && maxValue>params->threshold)
    			{
    				cvSeqPush(corners,&maxLoc);
    			
    			}
    					
    			x = x+params->maximumSize;
    
    		}
    		
    		y = y + params->maximumSize;
    
    	}
    	
    	delete []gaussWindow;
    	cvReleaseImage(&Ix);
    	cvReleaseImage(&Iy);
    	cvReleaseImage(&IxIx);
    	cvReleaseImage(&IyIy);
    	cvReleaseImage(&IxIy);
    	cvReleaseImage(&A);
    	cvReleaseImage(&B);
    	cvReleaseImage(&C);
    	cvReleaseImage(&cornerness);
    	cvReleaseImage(&srcFloat);
    }
    /*******************************
    *对源图像进行卷积运算
    *输入项
    *srcFloat    源图像
    *Ix          卷积的结果
    *dxTemplate  卷积模板
    *widthTemplate 模板的宽度
    *heightTemplate 模板的高度
    ********************************/
    void convolution(IplImage* srcFloat,IplImage* Ix,double* dxTemplate , int widthTemplate,int heightTemplate)
    {
    	int x,y,beginY,endY,beginX,endX;
    	
    	beginY = heightTemplate/2;
    	endY = srcFloat->height - heightTemplate/2;
    	
    	beginX = widthTemplate/2;
    	endX = srcFloat->width - widthTemplate/2;
    	
    	for(y=beginY;y<endY;y++)
    	{
    		for(x=beginX;x<endX;x++)
    		{
    	
    			int i,j;
    			double curDx=0;
    			
    
    			for(i=0;i<heightTemplate;i++)
    			{
    				for(j=0;j<widthTemplate;j++)
    				{
    					curDx += cvGetReal2D(srcFloat,y+i-heightTemplate/2,x+j-widthTemplate/2)**(dxTemplate+i*widthTemplate+j);
    				}
    			}		
    			cvSetReal2D(Ix,y,x,curDx);
    		
    		}
    		
    	}
    }

作者：小村长  出处：http://blog.csdn.net/lu597203933 欢迎转载或分享，但请务必声明文章出处。 （新浪微博：小村长zack, 欢迎交流！）

[SouthEast]: /images/20220828/6d30a489007d41a0979e611320ef9ef1.png
[SouthEast 1]: /images/20220828/394466728c254b3bb48c5da2f3db75cc.png
[SouthEast 2]: /images/20220828/9173f5807d9c4bafa6cebf7969d5d426.png
[SouthEast 3]: /images/20220828/1ff83db342004d3e9113d42c01839aff.png
[SouthEast 4]: /images/20220828/5ee69e865115463ea863141ad3945ca9.png
[SouthEast 5]: /images/20220828/8d00ba2aef524c7c8d8342c095b41f90.png
[SouthEast 6]: /images/20220828/8fb29a18e65b4f3aa6e4b6fa09eafff0.png
[SouthEast 7]: /images/20220828/7207a148f4944162858cbee4abba069d.png
[SouthEast 8]: /images/20220828/97955b2cc24e4a49a30fdb9acf1e1809.png