多尺度Harris角点函数opencv代码解析

多尺度Harris角点原理：

http://blog.csdn.net/tostq/article/details/49018777

http://www.cnblogs.com/ronny/p/4009425.html

函数代码：

void detectHarrisLaplace(const Mat& imgSrc, Mat& imgDst)
{
	Mat gray;
	if (imgSrc.channels() == 3)
	{
		cvtColor(imgSrc, gray, CV_BGR2GRAY);//彩色图转换为灰度图
	}
	else
	{
		gray = imgSrc.clone();//克隆原图
	}
	gray.convertTo(gray, CV_64F);// 图像像素转换为64位的双精度浮点数

	/* 尺度设置*/
	double dSigmaStart = 1.5;//尺度空间设置
	double dSigmaStep = 1.2;
	int iSigmaNb = 13;

	vector<double> dvecSigma(iSigmaNb);//存放积分尺度
	for (int i = 0; i < iSigmaNb; i++)
	{
		dvecSigma[i] = dSigmaStart + i*dSigmaStep;//尺度空间等差数列生成？
	}
	vector<Mat> harrisArray(iSigmaNb);//存放不同尺度对应的空间位置角点图
	
	for (int i = 0; i < iSigmaNb; i++)
	{
		double iSigmaI = dvecSigma[i];//取积分尺度
		double iSigmaD = 0.7 * iSigmaI;//计算微分尺度

		int iKernelSize = 6*round(iSigmaD) + 1;//round返回四舍五入的整数值,指定一维高斯的尺寸,高斯尺寸一般选取一边要大于3*sigma
		/*微分算子*/
		Mat dx(1, iKernelSize, CV_64F);//定义微分尺度对应的一维高斯模板
		for (int k =0; k < iKernelSize; k++)
		{
			int pCent = (iKernelSize - 1) / 2;//求中心点
			int x = k - pCent;//距离中心点的偏差
			dx.at<double>(0,i) = x * exp(-x*x/(2*iSigmaD*iSigmaD))/(iSigmaD*iSigmaD*iSigmaD*sqrt(2*CV_PI));//求取高斯值，同时取微分
		}
	
		Mat dy = dx.t(); //对dx进行转置
		Mat Ix,Iy;
		/*图像微分*/
		filter2D(gray, Ix, CV_64F, dx);//求偏微分，等于先与高斯一维函数卷积，再求偏微分
		filter2D(gray, Iy, CV_64F, dy);

		Mat Ix2,Iy2,Ixy;
		Ix2 = Ix.mul(Ix);//点乘
		Iy2 = Iy.mul(Iy);
		Ixy = Ix.mul(Iy);

		int gSize = 6*round(iSigmaI) + 1;//积分尺度对应的高斯模板尺寸
		Mat gaussKernel = getGaussianKernel(gSize, iSigmaI);//返回一个一维的高斯数组 
		filter2D(Ix2, Ix2, CV_64F, gaussKernel);//默认填充对称法
		filter2D(Iy2, Iy2, CV_64F, gaussKernel);
		filter2D(Ixy, Ixy, CV_64F, gaussKernel);

		/*自相关矩阵*/
		double alpha = 0.06;
		Mat detM = Ix2.mul(Iy2) - Ixy.mul(Ixy);
		Mat trace = Ix2 + Iy2;
		Mat cornerStrength = detM - alpha * trace.mul(trace);//计算自相关矩阵

		

		double maxStrength;
		minMaxLoc(cornerStrength, NULL, &maxStrength, NULL, NULL);//寻找最小最大值及位置，maxStrength为最大值
		Mat dilated;
		Mat localMax;
		dilate(cornerStrength, dilated, Mat());//膨胀，寻找最大值
		compare(cornerStrength, dilated, localMax, CMP_EQ);//将相关矩阵与最大值矩阵比较相等，相等的位置取255，不等的为0，存入localmax
	

		Mat cornerMap;
		double qualityLevel = 0.2;
		double thresh = qualityLevel * maxStrength;//设定相关矩阵的阈值
		cornerMap = cornerStrength > thresh;//相关矩阵二值化
		bitwise_and(cornerMap, localMax, cornerMap);//获取角点映射图
		harrisArray[i] = cornerMap.clone();	
	}

	/*计算尺度归一化Laplace算子*/
	vector<Mat> laplaceSnlo(iSigmaNb);
	for (int i = 0; i < iSigmaNb; i++)
	{
		double iSigmaL = dvecSigma[i];
		Size kSize = Size(6 * floor(iSigmaL) +1, 6 * floor(iSigmaL) +1);
		Mat hogKernel = getHOGKernel(kSize,iSigmaL);
		filter2D(gray, laplaceSnlo[i], CV_64F, hogKernel);
		laplaceSnlo[i] *= (iSigmaL * iSigmaL);//进行拉普拉斯响应计算
	}
	
	/*检测每个特征点在某一尺度LOG相应是否达到最大*/
	Mat corners(gray.size(), CV_8U, Scalar(0));
	for (int i = 0; i < iSigmaNb; i++)
	{
		for (int r = 0; r < gray.rows; r++)
		{
			for (int c = 0; c < gray.cols; c++)
			{
				if (i ==0)
				{
					if (harrisArray[i].at<uchar>(r,c) > 0 && laplaceSnlo[i].at<double>(r,c) > laplaceSnlo[i + 1].at<double>(r,c))
					{
						corners.at<uchar>(r,c) = 255;
					}
				}
				else if(i == iSigmaNb -1)
				{
					if (harrisArray[i].at<uchar>(r,c) > 0 && laplaceSnlo[i].at<double>(r,c) > laplaceSnlo[i - 1].at<double>(r,c))
					{
						corners.at<uchar>(r,c) = 255;
					}
				}
				else
				{
					if (harrisArray[i].at<uchar>(r,c) > 0 &&
					laplaceSnlo[i].at<double>(r,c) > laplaceSnlo[i + 1].at<double>(r,c) &&
					laplaceSnlo[i].at<double>(r,c) > laplaceSnlo[i - 1].at<double>(r,c))
					{
						corners.at<uchar>(r,c) = 255;
					}
				}
			}
		}
	}
	imgDst = corners.clone();
	
}