C++实现特征点提取、相关系数影像匹配、最小二乘影像匹配

写在前面

本学期数字摄影测量课程学习中，老师让同学实现书上提到的特征点提取、相关系数影像匹配、最小二乘影像匹配的相关算法，特征点提取确实有很多参考的，但是影像匹配的代码几乎都是付费下载的，自己也去找过，大多还是Matlab的代码。最后还是跟同学交流，然后用C++造轮子（因为选修的计算机视觉是用C++做的，所以这个就也用C++了），希望对有用到的同学有所帮助。

特征点提取

特征点提取有很多算子，例如Moravec、Harris、SIFT等等，毕竟造轮子的能力有限，就选择最简单的Moravec算子实现，下面介绍相关代码，因为在写代码的时候就打算发博客了，注释比较详细，就不啰嗦了。

1. 求兴趣值

//寻找数组最大值、最小值
void find(float a[], int m, float &max, float &min)
{
   
	min = a[0];
	max = a[0];
	for (int i = 0; i < m; i++)
	{
   
		if (a[i] > max)
		{
   
			max = a[i];
			continue;
		}
		else if (a[i] < min)
		{
   
			min = a[i];
			continue;
		}
	}
}
//求兴趣值，输入(图像矩阵、Moravec窗口大小、该像素的x,y)
float getInterestValue(cv::Mat m_srcimg, int Moravecsize, int i_x, int j_y)
{
   
	int halfsize = (Moravecsize) / 2;//定义小半个窗口大小
	float temp4[4];//创立四个方向的数组来进行"*"型的平方和
	//数组初始化
	for (int i = 0; i < 4; i++)
	{
   
		temp4[i] = 0;
	}
	//累加四个方向求平方和
	for (int i = 0; i < Moravecsize; i++)
	{
   
		float l = m_srcimg.at<uchar>(i_x - halfsize + i, j_y);
		temp4[0] += pow(m_srcimg.at<uchar>(i_x - halfsize + i, j_y) - m_srcimg.at<uchar>(i_x - halfsize + i + 1, j_y), 2);
		temp4[1] += pow(m_srcimg.at<uchar>(i_x, j_y - halfsize + i) - m_srcimg.at<uchar>(i_x, j_y - halfsize + i + 1), 2);
		temp4[2] += pow(m_srcimg.at<uchar>(i_x - halfsize + i, j_y - halfsize + i) - m_srcimg.at<uchar>(i_x - halfsize + i + 1, j_y - halfsize + i + 1), 2);
		temp4[3] += pow(m_srcimg.at<uchar>(i_x - halfsize + i, j_y + halfsize - i) - m_srcimg.at<uchar>(i_x - halfsize + i + 1, j_y + halfsize - i - 1), 2);
	}
	float min, max;
	find(temp4, 4, max, min);//给极小值赋值
	return min;//返回极小值，即该像素点的兴趣值
}

2. 特征点提取

int  Moravec(std::string path,std::vector<cv::Point3f> &featurePointLeft)
{
   
	cv::Mat imageRGB = cv::imread(path, cv::IMREAD_COLOR);
	if (imageRGB.empty())
	{
   
		std::cout << "Fail to read the image!" << path << std::endl;
		return -1;
	}
	cv::imshow("image", imageRGB);//显示原图
	cv::waitKey(0);
	cv::Mat imageGray;
	cv::cvtColor(imageRGB, imageGray, cv::COLOR_RGB2GRAY);
	std::vector<cv::Point3f> f;
	GaussianBlur(imageGray, imageGray, cv::Size(5, 5), 0, 0);//使用opencv自带高斯滤波预处理
	cv::Mat result/*结果32bit*/, resultNorm, resultNormUInt8;
	result = cv::Mat::zeros(imageGray.size(), CV_32FC1);

	int Moravecsize = 5;//计算Moravec兴趣值的窗口大小
	int Inhibitionsize = 9;//抑制局部最大的窗口大小

	float sum = 0;
	for (int i = 5; i < imageGray.rows - 5; i++)
	{
   
		for (int j = 5; j < imageGray.cols - 5; j++)
		{
   
			float min = getInterestValue(imageGray, Moravecsize, i, j);
			result.at<float>(i, j) = min;
			sum += min;
		}
	}
	//对小于阈值的该像素兴趣值置为零
	float mean = sum / (result.rows*result.cols);//经验阈值设置
	if (mean < 60) {
    mean = 300; }
	mean = 700;
	for (int i = 0; i <result.rows; i++)
	{
   
		for (int j = 0; j < result.cols; j++)
		{
   
			if (result.at<float>(i, j) < mean)
			{
   
				result.at<float>(i, j) = 0;
			}
		}
	}
	int halfInhibitionsize = Inhibitionsize / 2;//定义小半个窗口大小
	//抑制局部最大
	for (int i = halfInhibitionsize; i < result.rows - 1 - halfInhibitionsize; i++)
	{
   
		for (int j = halfInhibitionsize; j < result.cols - 1 - halfInhibitionsize; j++)
		{
   
			float temp1 = result.at<float>(i, j);
			for (int m = 0; m < Inhibitionsize; m++)
			{
   
				for (int n = 0; n < Inhibitionsize; n++)
				{
   
					float temp2 = result.at<float>(i - halfInhibitionsize + m, j - halfInhibitionsize + n);
					if (temp1 < temp2)
					{
   
						result.at<float>(i, j) = 0;
						n = Inhibitionsize;
						m = Inhibitionsize;
					}
				}
			}
		}
	}

	//存储特征点
	for (int i = halfInhibitionsize; i < result.rows - halfInhibitionsize - 1; i++)
	{
   
		for (int j = halfInhibitionsize; j < result.cols - halfInhibitionsize - 1; j++)
		{
   

			if (result.at<float>(i, j) > 0)
			{
   
				cv::Point3f temp;
				temp.x = i;
				temp.y = j;
				temp.z = result.at<float>(i, j);
				featurePointLeft.push_back(temp);
			}
		}
	}
	//开始删除同意一窗口值一样的重复点，尽管出现概率较小，但较大的图像往往某些窗口中会存在好几个数值相等的极大值
	for (int i = 0; i < featurePointLeft.size() - 1; i++)
	{
   
		for (int j = i + 1; j < featurePointLeft.size(); j++)
		{
   
			if ((featurePointLeft.at(i).z == featurePointLeft.at(j).z))
			{
   
				if (abs(featurePointLeft.at(i).x - featurePointLeft.at(j).x) < Inhibitionsize || abs(featurePointLeft.at(i).y - featurePointLeft.at(j).y) < Inhibitionsize)
				{
   
					featurePointLeft.erase(featurePointLeft.begin() + j);
					i = 0;
					break;
				}
			}
		}
	}
	//Normalizing image to 0-255
	cv::normalize(result, resultNorm, 0, 255, cv::NORM_MINMAX, CV_32FC1, cv::Mat());
	cv::convertScaleAbs(resultNorm, resultNormUInt8);
	//画图展示出来
	int radius = 5;
	for (size_t n = 0; n < featurePointLeft.size(); n++)
	{
   
		int i = int(featurePointLeft.at(n).y);
		int j = int(featurePointLeft.at(n).x);
		cv::circle(resultNormUInt8, cv::Point(i, j), radius, cv::Scalar(255), 1, 8, 0);
		cv::circle(imageRGB, cv::