Moraver算子提取影像特征点

01 Moraver算子

1.基本原理

        Moravec于1977年提出利用灰度方差提取点特征的算子，是一种基于灰度方差的角点检测方法，它通过计算图像中每个像素点沿着水平、垂直、对角线及反对角线的四个方向的灰度方差来检测角点。

2.算法步骤

（1）兴趣值计算

        对于图像中的每个像素点，以其为中心取一个窗口，计算该窗口内像素在四个方向上的灰度差的平方和。如公式（1）所示。（注意在此处c为列，r为行）。

         其中K为窗口大小一半取整，对于整幅图像来说，窗口大小可以忽略不计，因此不考虑边缘像素处理。

        取计算结果的最小值作为兴趣值。

（2）阈值筛选

        给定一个经验阈值，将兴趣值大于该阈值的点作为候选点。

（3）非极大值抑制

        给定另一个窗口，选取兴趣值最大的点，作为特征点。

02 C语言(c++)代码实现

#include <iostream>
#include <stdio.h>
#include <omp.h>
#include <opencv2/opencv.hpp>
#include <iomanip>

using namespace cv;
using namespace std;
//文件读取，灰度图或者RGB图像
void readData(const char* filePath, int r, Mat& img) {
	if (r == 1)
		img = imread(filePath, IMREAD_GRAYSCALE);
	else
		img = imread(filePath);
}
//计算兴趣值
void Moravec_iv(Mat& iv, Mat img, int win_size) {
	int k = win_size / 2;
	for (int i = k; i < img.rows - k; i++) {
		for (int j = k; j < img.cols - k; j++) {
			double v1 = 0, v2 = 0, v3 = 0, v4 = 0;
			for (int l = -k; l <= k - 1; l++) {
				double temp1 = img.at<uchar>(i, j + l) - img.at<uchar>(i, j + l + 1);
				v1 += temp1 * temp1;
				double temp2 = img.at<uchar>(i + l, j + l) - img.at<uchar>(i + l + 1, j + l + 1);
				v2 += temp2 * temp2;
				double temp3 = img.at<uchar>(i + l, j) - img.at<uchar>(i + l + 1, j);
				v3 += temp3 * temp3;
				double temp4 = img.at<uchar>(i - l, j + l) - img.at<uchar>(i - l - 1, j + l + 1);
				v4 += temp4 * temp4;
			}
			iv.at<double>(i, j) = fmin(fmin(v1, v2), fmin(v3, v4));
		}
	}
}
//阈值筛选
void Moravec_threshold_value_points(Mat& iv, Mat img, int win_size, double threshold_value) {
	for (int i = 0; i < img.rows; i++) {
		for (int j = 0; j < img.cols; j++) {
			if (iv.at<double>(i, j) > threshold_value)
				iv.at<double>(i, j) = iv.at<double>(i, j);
			else
				iv.at<double>(i, j) = 0;
		}
	}
}
//寻找抑制窗口内的极大值
void find_max_in_window(Mat img, Mat iv, int window_size, int row, int col, int& max_row, int& max_col) {
	int k = window_size / 2;
	double max_value = 0;
	max_row = -1;
	max_col = -1;

	for (int i = -k; i <=k; i++) {
		for (int j = -k; j <= k; j++) {
			int r = row + i;
			int c = col + j;
			if (r >= 0 && r < img.rows && c >= 0 && c < img.cols) {
				double value = iv.at<double>(r, c);
				if (value > max_value) {
					max_value = value;
					max_row = r;
					max_col = c;
				}
			}
		}
	}
}
//特征点选取
void Moravec_extremum_points(Mat iv, Mat img, Mat& mask, int win_size1, int win_size2, double threshold_value) {
	int k = win_size2 / 2;
    //注意为整个窗口滑动而不是逐个像素滑动
	for (int i = k; i < img.rows - k; i+= win_size2){
		for (int j = k; j < img.cols - k; j+= win_size2){
			int max_row, max_col;
			find_max_in_window(img, iv, win_size2, i, j, max_row, max_col);
			if (max_row >= 0 && max_row < img.rows && max_col >= 0 && max_col < img.cols)
			mask.at<double>(max_row, max_col) = 1;
		}
	}
}

int main() {
	Mat img;
	const char* filePath = "img.jpg";
	readData(filePath, 1, img);
	if (img.empty()) {
		cerr << "无法读取图像文件: " << filePath << endl;
		return -1;
	}
	cout << "文件已经读取" << endl;

	Mat iv = Mat::zeros(img.rows, img.cols, CV_64F);
	Mat mask = Mat::zeros(img.rows, img.cols, CV_64F);
	int win_size = 5;//兴趣值计算窗口
	int win_size2 = 100;//抑制窗口
	int threshold_value = 2000;//经验阈值
	Moravec_iv(iv, img, win_size);
	cout << "兴趣值已经计算" << endl;
	Moravec_threshold_value_points(iv, img, win_size, threshold_value);
	cout << "候选点已经选取" << endl;
	Moravec_extremum_points(iv, img, mask, win_size, win_size2, threshold_value);
	cout << "极值点已经选择" << endl;

	Mat img_with_points;
	readData(filePath, 2, img_with_points);
	int count = 0;
	for (int i = 0; i < mask.rows; i++) {
		for (int j = 0; j < mask.cols; j++) {
			if (mask.at<double>(i, j) == 1) {
				circle(img_with_points, Point(j, i), 2.5, Scalar(0, 0, 255), -1);
				string pointIndex = to_string(count);
				putText(img_with_points, pointIndex, Point(j + 5, i - 5), FONT_HERSHEY_SIMPLEX, 1, Scalar(0, 0, 255), 1);//根据图像大小需要调整字体，确保正确显示
				count++;
			}
		}
	}

	cout << "共有" << count << "个特征点" << endl;
    //调整图像显示比例
	int screenWidth = 1920;
	int screenHeight = 1080;
	int winWidth = img.cols, winHeight = img.rows;
	double imgRatio = static_cast<double>(img.cols) / img.rows;
	
	if (img.cols > screenWidth)
	{
		winWidth = screenWidth;
		winHeight = winWidth / imgRatio;
	}
	else if (img.rows > screenHeight)
	{
		winHeight = screenHeight;
		winWidth = winHeight * imgRatio;
	}

	namedWindow("Image with Points", WINDOW_NORMAL);
	resizeWindow("Image with Points", winWidth, winHeight);
	imshow("Image with Points", img_with_points);
	waitKey(0);
	destroyAllWindows();
	return 0;
}

03 结果展示

        原始图像：

        提取特征点叠加之后的图像：