Moravec角点检测

(1)Moravec角点检测算法原理
Moravec角点检测算法是最早的角点检测算法之一。该算法将角点定义为具有低“自相关性”的点。算法会检测图像的每一个像素，将像素周边的一个邻域作为一个patch，并检测这个patch和周围其他patch的相关性。这种相关性通过两个patch间的平方差之和(SSD)来衡量，SSD值越小则相似性越高。
如果像素位于平滑图像区域内，周围的patch都会非常相似。如果像素在边缘上，则周围的patch在与边缘正交的方向上会有很大差异，在与边缘平行的方向上则较为相似。而如果像素是各个方向上都有变化的特征点，则周围所有的patch都不会很相似。
Moravec会计算每个像素patch和周围patch的SSD最小值作为强度值，取局部强度最大的点作为特征点。

(2)Moravec角点检测示意图：


(3)将moravec角点检测方法公式化：
Moravec 在1981年提出Moravec角点检测算子，并将它应用于立体匹配。
首先, 计算每个像素点的兴趣值, 即以该像素点为中心, 取一个w*w(如:5x5)的方形窗口, 计算0度、45度、90度、135度四个方向灰度差的平方和, 取其中的最小值作为该像素点的兴趣值。


如图以3x3为例　黑色窗口为Ｉ(x,y)　红色窗口为I(x+u,y+v)
其中四种移位 (u,v) = (1,0), (1,1), (0,1), (-1, 1).w(x,y)为方形二值窗口,若像素点在窗口内，则取值为1, 否则为0。

(4)moravec角点检测步骤：
(1)对于每一个像素点，计算在E(u,v),在我们的算法中，(u,v)的取值是((1,0), (1,1), (0,1), (-1, 1).当然，你自己可以改成(1,0),(1,1),(0,1),(-1,1),(-1,0),(-1,-1),(0,-1),(1,-1) 8种情况
(2)计算最小值对每个位置minValue = min{E(u,v)},其中(u,v) = (1,0), (1,1), (0,1), (-1, 1).
(3)对每个位置minValue 进行判断，是不是大于设定阈值，如果是大于设定阈值，接着判断是不是局部极大值

x-2,y-2	x-1,y-2	x,y-2	x+1,y-2	x+2,y-2
x-2,y-2	x-1,y-1	x,y-1	x+1,y-1	x+2,y-1
x-2,y	x-1,y	x,y	x+1,y	x+2,y
x-2,y+1	x-1,y+1	x,y+1	x+1,y	x+2,y
x-2,y+2	x-1,y+2	x,y+2	x+1,y+2	x+2,y+2

(A1-B1)+(A2-B2)+(A3-B3)+(A4-B4)
A1到A4的 x坐标变化范围：[x-2,x+1]
A1到A4的y坐标不变



(A1-B1)+(A2-B2)+(A3-B3)+(A4-B4)
A1到A4的x坐标变化范围：[x-2,x+2]
A1到A4的y坐标变化范围：[y-2,y+2]



A1到A4的y坐标变化范围：[y-2,y+1]



A1到A4的x坐标变化范围：[x+2,x-1]
A1到A4的y坐标变化范围: [y-2,y+1]

注意在实际使用时还要要进行非局部极大值抑制

下面给出完整代码：

#include "stdafx.h"
#include "cv.h"
#include "cxcore.h"
#include "highgui.h"

/**********************************************************************************
*函数 int getMoravec(IplImage* src,CvSeq* corners)
*输入:
*src : 单通道图像
*corners : 用来保存提取到的角点
*threshold : 角点量的阈值 (具体函数 请看视频)
*输出
*corners : 用来保存提取到的角点
*返回值
*角点的个数
***************************************************************************************/
int getMoravec(IplImage* src,CvSeq* corners , float threshold)
{
	//窗口大小
	const int winSize=5; 
	int x,y,halfWinSize=winSize/2;

	//保存最小的变化量
	IplImage* diffDst = cvCreateImage(cvGetSize(src),32,1);
	cvZero(diffDst);

	//保存角点个数
	int cornersCount=0;

	//1.计算图像上每一个点上的 在4个方向上的变化量 并计算出最小值 保存在矩阵diffDst中
	for(y=halfWinSize;y<src->height-halfWinSize;y++)
	{
		for(x=halfWinSize;x<src->width-halfWinSize;x++)
		{
			//compute the reaction in the four directions(0,45,90,135)
			int winx;

			//数组reaction[4] 用于保持 在四个方向上的灰度值变化量 
			//minValue用于保存4个变化量中的最小值
			float reaction[4],minValue;
			reaction[0]=0;
			reaction[1]=0;
			reaction[2]=0;
			reaction[3]=0;

			//提示  下面的4个循环 可以综合成一个循环

			//0 度方向的变化量保存在reaction[0]中
			for( winx=-halfWinSize;winx<halfWinSize;winx++)
			{
				reaction[0] = reaction[0] + pow( cvGetReal2D(src,y,x+winx)-cvGetReal2D(src,y,x+winx+1), 2 );
			}

			//45 度方向的变化量保存在reaction[1]中
			for( winx=-halfWinSize;winx<halfWinSize;winx++)
			{
				reaction[1] = reaction[1]+pow(cvGetReal2D(src,y+winx,x+winx)-cvGetReal2D(src,y+winx+1,x+winx+1),2);
			}

			//90 度方向的变化量保存在reaction[2]中
			for( winx=-halfWinSize;winx<halfWinSize;winx++)
			{
				reaction[2] = reaction[2]+pow(cvGetReal2D(src,y+winx,x)-cvGetReal2D(src,y+winx+1,x),2);
			}

			//135 度方向的变化量保存在reaction[3]中
			for( winx=-halfWinSize;winx<halfWinSize;winx++)
			{
				reaction[3] = reaction[3]+pow(cvGetReal2D(src,y+winx,x-winx)-cvGetReal2D(src,y+winx+1,x-winx-1),2);
			}

			//计算4个量中最小值  保存到minValue
			minValue = reaction[0];
			minValue = minValue > reaction[1] ? reaction[1] : minValue;
			minValue = minValue > reaction[2] ? reaction[2] : minValue;
			minValue = minValue > reaction[3] ? reaction[3] : minValue;

			//将最小的变化量保存到矩阵
			cvSetReal2D(diffDst,y,x,minValue);
		}
	}

	//2.获取角点坐标
	for(y=halfWinSize;y<src->height-halfWinSize;)
	{
		for(x=halfWinSize;x<src->width-halfWinSize;)
		{
			float max=0;
			int flag = 0 ;
			CvPoint maxLoc;
			maxLoc.x = -1;
			maxLoc.y = -1;

			//首先计算以点(x,y)位中心的winSize*winSize的窗口内部的局部极大值
			for(int winy=-halfWinSize;winy<=halfWinSize;winy++)
			{
				for(int winx=-halfWinSize;winx<=halfWinSize;winx++)
				{
					float value ;
					value = cvGetReal2D(diffDst,y+winy,x+winx);

					//计算该窗口内 最大值 保存到max 并保存其坐标到maxLoc
					if(value>max)
					{
						max = value;
						maxLoc.x = x+winx;
						maxLoc.y = y+winy;
						flag = 1;
					}
				}
			}

			//如果找到局部极大值 并且该值大于预先设定的阈值 则认为是角点
			if(flag==1 && max>threshold)
			{
				cvSeqPush(corners,&maxLoc);
				cornersCount++;	
			}

			//下一个窗口
			x=x+halfWinSize;

		}

		//下一行的第一个窗口
		y=y+halfWinSize;
	}

	cvReleaseImage(&diffDst);

	return cornersCount;
}


int main(int argc, char* argv[])
{
	IplImage* src;
	//加载源图像
	src = cvLoadImage("images/1.BMP",CV_LOAD_IMAGE_GRAYSCALE);
	if(!src)
	{
		printf("图像加载失败");
	}

	//用于保存最终角点的空间
	CvMemStorage* mem = cvCreateMemStorage(0);

	//角点将会保存在一个CvSeq中
	CvSeq* corners;
	corners = cvCreateSeq(0,sizeof(CvSeq),sizeof(CvPoint),mem);

	//角点的个数
	int cornersCount;

	//调用函数getMoravec计算角点
	cornersCount = getMoravec(src,corners,29500);

	//图像show用于显示角的提取结果
	IplImage* show= cvCreateImage(cvGetSize(src),8,3);
	cvCvtColor(src,show,CV_GRAY2BGR);

	//获取每一个角点的坐标
	for(int x=0;x<cornersCount;x++)
	{
		CvPoint* pt = (CvPoint*)cvGetSeqElem(corners,x);

		//以角点坐标为中心  绘制一个半径为5的圆
		cvCircle(show,*pt,5,cvScalar(255,0,255,0));
	}

	//显示结果
	cvNamedWindow("dst");
	cvShowImage("dst",show);
	cvWaitKey(0);

	cvReleaseImage(&src);
	cvReleaseImage(&show);
	cvReleaseMemStorage(&mem);

	return 0;
}