使用C++的OpenCV实现Halcon的PointsHarris

许久没有更新了

今天就来更新下关于使用OpenCV实现Halcon的PointsHarris

具体公式参考Halcon的官方文档

Halcon points_harris算法介绍

当然这个算法运行速度会很慢 因为我没做任何优化
只是提供一种实现的思路

优化的话 可以考虑把求导和求特征值合并计算
用AVX2或者 AVX512对单精度浮点进行SIMD加速

void derivateGauss(Mat &Input, Mat &Output, const std::string &mode, double Sigma);

void caclCornerSubPix(Mat &Input,double SigmaGrad,double SigmaSmooth,std::vector<Point2d> &inoutPoints);

void points_harris(Mat &Input, Mat &Output, double SigmaGrad, double SigmaSmooth, double Alpha, double Threshold,
	std::vector<Point2d>&harrisPoints
) {

	///彩色通道可自行根据下列算法扩展
	if (Input.channels() != 1) {
		return;
	}

	///计算关于x和y方向的高斯一阶偏导数
	Mat Ix, Iy;
	derivateGauss(Input, Ix, "x", SigmaGrad);
	derivateGauss(Input, Iy, "y", SigmaGrad);

	int rows = Input.rows;
	int cols = Input.cols;

	///协方差矩阵
	Mat Cov(rows,cols,CV_32FC3);
	for (int i = 0; i < rows; i++) {

		float *IxPtr = Ix.ptr<float>(i);
		float *IyPtr = Iy.ptr<float>(i);

		float *CovPtr = Cov.ptr<float>(i);
		for (int j = 0; j < cols; j++) {

			float IxValue = IxPtr[j];
			float IyValue = IyPtr[j];

			CovPtr[j * 3] = IxValue * IxValue;
			CovPtr[j * 3+1] = IxValue * IyValue;
			CovPtr[j * 3+2] = IyValue * IyValue;
		}
	}

	Mat CovFilter;
	derivateGauss(Cov,CovFilter,"none",SigmaSmooth);

	///DetM = randa1 *randa2
	///TraceM= randa1+randa2;
	for (int i = 0; i < rows; i++) {

		float *CovFilterPtr = CovFilter.ptr<float>(i);

		for (int j = 0; j < cols; j++) {

			float IxIx = CovFilterPtr[j * 3];
			float IxIy = CovFilterPtr[j * 3+1];
			float IyIy = CovFilterPtr[j * 3+2];

			Mat M(2, 2, CV_32F);
			M.at<float>(0, 0) = IxIx;
			M.at<float>(0, 1) = IxIy;
			M.at<float>(1, 0) = IxIy;
			M.at<float>(1, 1) = IyIy;

			///特征值
			Mat eigenValue;
			eigen(M,eigenValue);

			float randa1 = eigenValue.ptr<float>(0)[0];
			float randa2 = eigenValue.ptr<float>(1)[0];

			///计算矩阵行列式值
			float detM = randa1 * randa2;

			///计算矩阵迹
			float traceM = randa1 + randa2;

			float R = detM - Alpha * traceM*traceM;


			if (R > Threshold) {

				harrisPoints.push_back(Point2d(j,i));
			}

		}
	}

	///求二阶导数 计算亚像素点
	caclCornerSubPix(Input, SigmaGrad, SigmaSmooth, harrisPoints);


}