OpenCV实践之特征点评估代码解析

重复率评估指标简述

　　我的博客写了一些特征匹配（尺度+仿射）相关匹配算法，但是一直都没有写过如何评估特征点的鲁棒性策略。如何进行特征点提取算法的鲁棒性评估是写这篇博客目的所在。重复率指标就是采用最多的特征点检测性能的评估策略：

$$repeatability=\frac{\sum dist(m_{points}^{{}^{\prime }},n_{points})}{min(m_{points}^{{}^{\prime }},n_{points}^{{}^{\prime }})}$$

　　参数解释：$H$为图像变换映射矩阵，$m_{points}$和$n_{points}$分别为算法在两幅图像检测出的点集。$m_{points}^{{}^{\prime }}$和$n_{points}^{{}^{\prime }}$为经过一些约束条件后，剩余的特征点集（下面会详细介绍评估的过程）。

重复率计算流程：

假设$img1$和$img2$图像经过特征检测算法后点集分别为$m_{points}$和$n_{points}$。

　　$Step1$: 首先将$img1$中的$m_{points}$个特征点集分别乘以变换矩阵$H$，生成新的坐标点集$m_{pointsH}$；

　　$Step2$: 新的坐标点集$m_{point{s}_H}$经过$img2$图像范围约束筛选后，剩余点集为$m_{points}^{{}^{\prime }}$； 筛选条件为：$m_{pointsH}$不能小于0大于$img2$图像的宽高。

$$0<x_i<img2.width$$

$$0<y_i<img2.height$$

其中，$m(x_i,y_i)\in m_{pointsH}$

　　$Step3$: 将$img2$中的$n_{points}$个特征点集分别乘以变换矩阵$H^{-1}$，生成新的坐标点集$n_{pointsH}$；

　　$Step4$: 新的坐标点集$n_{pointsH}$经过$img1$图像范围约束筛选后，剩余点集为$n_{points}^{{}^{\prime }}$；筛选条件为：$n_{pointsH}$不能小于0大于$img1$的宽高。

$$0<x_i<img1.width$$

$$0<y_i<img1.height$$

其中，$n(x_i,y_i)\in n_{pointsH}$

　　到此处，我们已经获取重复率计算公式的分母参数$min(m_{points}^{{}^{\prime }},n_{points}^{{}^{\prime }})$。下面我们解释一下分子约束条件：

　　分子中$\sum dist(m_{points}^{{}^{\prime }},n_{points}^{{}^{\prime }})$其实$dist$为计算点之间的欧式距离，将$m_{points}^{{}^{\prime }}$点集分别于$n_{points}$计算欧式距离，如果欧式距离阈值小于设定值$\epsilon$，那么就是一个重复点。最后，对所有满足欧氏距离阈值的点集求和看有多少对。

$$repeatability=\frac{\sum dist(m_{points}^{{}^{\prime }},n_{points})}{min(m_{points}^{{}^{\prime }},n_{points}^{{}^{\prime }})}$$

　　后来论文An Affine invariant interest point detector对重复率评估进行了延伸，针对仿射不变的特征点评估为：

$$repeatability=\frac{\sum correspondence}{min(m_{points}^{{}^{\prime }},n_{points}^{{}^{\prime }})}$$

延伸后，分子$correspondence$满足下述两个条件即可：

　　1 与上述的欧式距离统计筛选原理一致；

　　2 特征点区域映射到另一幅图像中重叠误差小于0.2（默认参数）；

区域overlap重叠比率示意图：

加上限制条件2主要考虑原因如下：

　　仿射不变特征点检测器是基于特征点附近的区域进行综合判断来获取特征点（例如：SIFT、Harris-Affine、MSER等），当特征点与映射矩阵H点乘后，特征点附近区域也会相应发生变化，则需要进一步判断区域映射是否与另一幅特征点区域重复比例。

特征点评估函数解析及其代码

　　评估主函数就是$evaluateFeatureDetector()$函数，相关接口参数解释如下：$img1$与 $mg2$为输入图像，$H1to2$为img1与img2的变换映射矩阵，$keypoints1$与$keypoints2$分别为img1与img2的特征点数据结构（下面具体会解释），$repeatability$与$correspCount$为重复率与一致性点个数，$fdetector$为特征点检测算法声明。

　　当输入$keypoints1$与$keypoints2$为空时候，$fdetector$必须指定某一种特征点检测算法；相反，$keypoints1$与$keypoints2$不为空时候，$fdetector$是否指定都没有意义。具体可以参考$evaluateFeatureDetector()$函数内部如下代码片段：

if( (keypoints1->empty() || keypoints2->empty()) && !fdetector )
        CV_Error( Error::StsBadArg, "fdetector must not be empty when keypoints1 or keypoints2 is empty" );

    if( keypoints1->empty() )
        fdetector->detect( img1, *keypoints1 );
    if( keypoints2->empty() )
        fdetector->detect( img2, *keypoints2 );

下面参考一下评估简短代码：

#include <opencv2/features2d.hpp>
#include <opencv2/opencv.hpp>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/xfeatures2d.hpp>
#include <iostream>

using namespace std;
using namespace cv;

int main(void)
{	
	// Read img1 and img2 from path 
	string filename1 = "..\\featuresDetectorEval\\image\\img1.jpg";
	string filename2 = "..\\featuresDetectorEval\\image\\img3.jpg";
	Mat img1 = imread(filename1);
	Mat img2 = imread(filename2);

	// keypoints1 keypoints2 is initializing 
	vector<cv::KeyPoint> *keypoints1 = NULL, *keypoints2 = NULL;

	// Homography MatRix
	string xml_path = "..\\featuresDetectorEval\\image\\H1to3p.xml";
	FileStorage fs(xml_path, FileStorage::READ);
	// H1to2 matrix get datat from H1to3p.xml 
	Mat H1to2;
	fs.getFirstTopLevelNode() >> H1to2;

	float repeatability = 0.0;
	int correspCount = 0;

	//Ptr<FeatureDetector> _fdetector = FastFeatureDetector::create();
	Ptr<Feature2D> _fdetector = xfeatures2d::MSDDetector::create();
	//Ptr<Feature2D> _fdetector = xfeatures2d::StarDetector::create();
	//Ptr<Feature2D> _fdetector = xfeatures2d::HarrisLaplaceFeatureDetector::create();
	//Ptr<Feature2D> _fdetector = xfeatures2d::SIFT::create();
	//Ptr<Feature2D> _fdetector = xfeatures2d::SURF::create();
	//Ptr<Feature2D> _fdetector = ORB::create();
	//Ptr<Feature2D> _fdetector = MSER::create();
	//Ptr<Feature2D> _fdetector = AgastFeatureDetector::create();
	//Ptr<Feature2D> _fdetector = GFTTDetector::create();
	//Ptr<Feature2D> _fdetector = SimpleBlobDetector::create();

	evaluateFeatureDetector(img1, img2, H1to2, keypoints1, keypoints2, repeatability, correspCount, _fdetector);

	cout << "repeatability is " << repeatability << "  correspCount is : " << correspCount << endl;

	system("pause");
	return 0;
}

各个特征点检测算法重复率与一致性点数实验对比：

　　根据上述对比有点小诧异，最大相似性算法$MSDDetector$重复率最高，$ORB$、$SIFT$、$SURF$紧随其后，$MSER$算法的性能中等。按照算法设计上来看，$MSER$应该性能较佳，对比结果来看并不是这样。当然，可能原因在于我只是简单测试一幅图像，要想系统了解各个算法的检测性能，还是需要多测试一下不同类别的图像。

下面是源码解析，如果只是想调用的话，就可以跳过这部分。
---------------------------------------源码解析-------------------------------------------

评估主函数源码：

void cv::evaluateFeatureDetector(const Mat& img1, const Mat& img2, const Mat& H1to2,
	std::vector<KeyPoint>* _keypoints1, std::vector<KeyPoint>* _keypoints2,
	float& repeatability, int& correspCount,
	const Ptr<FeatureDetector>& _fdetector)
{
	CV_INSTRUMENT_REGION();
	//-- 特征点检测模式
	Ptr<FeatureDetector> fdetector(_fdetector);
	//-- 特征点数据结构初始化
	std::vector<KeyPoint> *keypoints1, *keypoints2, buf1, buf2;
	keypoints1 = _keypoints1 != 0 ? _keypoints1 : &buf1;
	keypoints2 = _keypoints2 != 0 ? _keypoints2 : &buf2;
	// 判断输入的keypoints1与keypoints2是否是空  同时fdetector是否初始化
	if ((keypoints1->empty() || keypoints2->empty()) && !fdetector)
		CV_Error(Error::StsBadArg, "fdetector must not be empty when keypoints1 or keypoints2 is empty");
	//-- 如果keypoints1为空的话，进行检测
	if (keypoints1->empty())
		fdetector->detect(img1, *keypoints1);
	//-- 如果keypoints2为空的话，进行检测
	if (keypoints2->empty())
		fdetector->detect(img2, *keypoints2);
	//-- 计算特征点的重复比例与正确的点数
	calculateRepeatability(img1, img2, H1to2, *keypoints1, *keypoints2, repeatability, correspCount);
}

主函数内部调用计算重复率的函数：

static void calculateRepeatability(const Mat& img1, const Mat& img2, const Mat& H1to2,
	const std::vector<KeyPoint>& _keypoints1, const std::vector<KeyPoint>& _keypoints2,
	float& repeatability, int& correspondencesCount,
	Mat* thresholdedOverlapMask = 0)
{
	std::vector<EllipticKeyPoint> keypoints1, keypoints2, keypoints1t, keypoints2t;
	//-- 检测特征点进行椭圆初始化格式
	EllipticKeyPoint::convert(_keypoints1, keypoints1);
	EllipticKeyPoint::convert(_keypoints2, keypoints2);

	//-- 根据变换矩阵求出坐标映射
	EllipticKeyPoint::calcProjection(keypoints1, H1to2, keypoints1t);
	Mat H2to1; invert(H1to2, H2to1);
	EllipticKeyPoint::calcProjection(keypoints2, H2to1, keypoints2t);

	float overlapThreshold;
	bool ifEvaluateDetectors = thresholdedOverlapMask == 0;
	if (ifEvaluateDetectors)
	{
		overlapThreshold = 1.f - 0.4f;

		//-- 移除映射后特征点超出图像边界范围
		Size sz1 = img1.size(), sz2 = img2.size();
		//-- 根据特征点过滤掉一些限制条件的点，具体见filterEllipticKeyPointsByImageSize()解析
		filterEllipticKeyPointsByImageSize(keypoints1, sz1);
		filterEllipticKeyPointsByImageSize(keypoints1t, sz2);
		filterEllipticKeyPointsByImageSize(keypoints2, sz2);
		filterEllipticKeyPointsByImageSize(keypoints2t, sz1);
	}
	else
	{
		overlapThreshold = 1.f - 0.5f;

		thresholdedOverlapMask->create((int)keypoints1.size(), (int)keypoints2t.size(), CV_8UC1);
		thresholdedOverlapMask->setTo(Scalar::all(0));
	}
	size_t size1 = keypoints1.size(), size2 = keypoints2t.size();
	//-- 映射特征点进行约束后，重复率分母因子为min二者较小的特征点数
	size_t minCount = MIN(size1, size2);

	//-- 计算区域重复比率
	std::vector<SIdx> overlaps;
	//-- 计算特征点集之间的椭圆交叉比率
	computeOneToOneMatchedOverlaps(keypoints1, keypoints2t, ifEvaluateDetectors, overlaps, overlapThreshold);

	correspondencesCount = -1;
	repeatability = -1.f;
	if (overlaps.empty())
		return;

	if (ifEvaluateDetectors)
	{
		//-- regions one-to-one matching
		correspondencesCount = (int)overlaps.size();
		repeatability = minCount ? (float)correspondencesCount / minCount : -1;
	}
	else
	{
		for (size_t i = 0; i < overlaps.size(); i++)
		{
			int y = overlaps[i].i1;
			int x = overlaps[i].i2;
			thresholdedOverlapMask->at<uchar>(y, x) = 1;
		}
	}
}

区域overlaps计算函数：

static void computeOneToOneMatchedOverlaps(const std::vector<EllipticKeyPoint>& keypoints1, const std::vector<EllipticKeyPoint>& keypoints2t,
	bool commonPart, std::vector<SIdx>& overlaps, float minOverlap)
{
	CV_Assert(minOverlap >= 0.f);
	overlaps.clear();
	if (keypoints1.empty() || keypoints2t.empty())
		return;

	overlaps.clear();
	overlaps.reserve(cvRound(keypoints1.size() * keypoints2t.size() * 0.01));

	for (size_t i1 = 0; i1 < keypoints1.size(); i1++)
	{
		EllipticKeyPoint kp1 = keypoints1[i1];
		float maxDist = sqrt(kp1.axes.width*kp1.axes.height),
			fac = 30.f / maxDist;
		if (!commonPart)
			fac = 3;

		maxDist = maxDist * 4;
		fac = 1.f / (fac*fac);

		EllipticKeyPoint keypoint1a = EllipticKeyPoint(kp1.center, Scalar(fac*kp1.ellipse[0], fac*kp1.ellipse[1], fac*kp1.ellipse[2]));

		for (size_t i2 = 0; i2 < keypoints2t.size(); i2++)
		{
			EllipticKeyPoint kp2 = keypoints2t[i2];
			Point2f diff = kp2.center - kp1.center;

			if (norm(diff) < maxDist)
			{
				EllipticKeyPoint keypoint2a = EllipticKeyPoint(kp2.center, Scalar(fac*kp2.ellipse[0], fac*kp2.ellipse[1], fac*kp2.ellipse[2]));
				//-- find the largest eigenvalue
				int maxx = (int)ceil((keypoint1a.boundingBox.width > (diff.x + keypoint2a.boundingBox.width)) ?
					keypoint1a.boundingBox.width : (diff.x + keypoint2a.boundingBox.width));
				int minx = (int)floor((-keypoint1a.boundingBox.width < (diff.x - keypoint2a.boundingBox.width)) ?
					-keypoint1a.boundingBox.width : (diff.x - keypoint2a.boundingBox.width));

				int maxy = (int)ceil((keypoint1a.boundingBox.height >(diff.y + keypoint2a.boundingBox.height)) ?
					keypoint1a.boundingBox.height : (diff.y + keypoint2a.boundingBox.height));
				int miny = (int)floor((-keypoint1a.boundingBox.height < (diff.y - keypoint2a.boundingBox.height)) ?
					-keypoint1a.boundingBox.height : (diff.y - keypoint2a.boundingBox.height));
				int mina = (maxx - minx) < (maxy - miny) ? (maxx - minx) : (maxy - miny);

				//--compute the area
				float dr = (float)mina / 50.f;
				int N = (int)floor((float)(maxx - minx) / dr);
				//-- 交叉区域统计
				IntersectAreaCounter ac(dr, minx, miny, maxy, diff, keypoint1a.ellipse, keypoint2a.ellipse);
				parallel_reduce(BlockedRange(0, N + 1), ac);
				if (ac.bna > 0)
				{
					float ov = (float)ac.bna / (float)ac.bua;
					if (ov >= minOverlap)
						overlaps.push_back(SIdx(ov, (int)i1, (int)i2));
				}
			}
		}
	}

	std::sort(overlaps.begin(), overlaps.end());

	typedef std::vector<SIdx>::iterator It;

	It pos = overlaps.begin();
	It end = overlaps.end();

	while (pos != end)
	{
		It prev = pos++;
		end = std::remove_if(pos, end, SIdx::UsedFinder(*prev));
	}
	overlaps.erase(pos, overlaps.end());
}

过滤不符合要求的特征点：

static void filterEllipticKeyPointsByImageSize(std::vector<EllipticKeyPoint>& keypoints, const Size& imgSize)
{
	if (!keypoints.empty())
	{
		std::vector<EllipticKeyPoint> filtered;
		filtered.reserve(keypoints.size());
		std::vector<EllipticKeyPoint>::const_iterator it = keypoints.begin();
		for (int i = 0; it != keypoints.end(); ++it, i++)
		{   //-- 中心点坐标加或者减去椭圆的宽、高，是否越界
			if (it->center.x + it->boundingBox.width < imgSize.width &&
				it->center.x - it->boundingBox.width > 0 &&
				it->center.y + it->boundingBox.height < imgSize.height &&
				it->center.y - it->boundingBox.height > 0)
				filtered.push_back(*it);
		}
		keypoints.assign(filtered.begin(), filtered.end());
	}
}

小结

　　特征点检测算法的评估一般是重复率指标，其实各个算法针对不同的应用场景，例如像基础的特征点检测算法：$FAST$、$Harris$、$Hessian$等并不支持尺度的变换，如果图像对存在一定的尺度变换，那么该类别算子重复率可能就会不佳。还有$SIFT$、$SURF$等具备尺度变换的性能，但是并不具备仿射的因素，$HarrisAffine$与$HessianAffine$与$MSER$算子针对仿射的因子，可能就会对仿射变换的图像性能体现较佳的效果，但是也带来极大的计算时间。介绍完特征点检测算法的性能评估，那么描述符该用什么指标来进行评估呢？下篇博客将会介绍一下描述符的评估策略。最后，如有错误，还请批评指正！

参考

https://docs.opencv.org/4.0.0-beta/da/d9b/group__features2d.html