Algorithm-Repo/SourceCode/CVAlgorithmUtils/LesionCompare/SIFTMatch/FeatureMatch.cpp

#include "FeatureMatch.h"

/// 高斯核卷积平滑。img 是输入图像，ikernel 是高斯核
/// 返回值为平滑图像
Mat Conv2(const Mat &img, const Mat &ikernel)
{
	Mat dest;
	Mat kernel;
	cv::flip(ikernel, kernel, -1);
	Mat source = img;
	Point anchor(kernel.cols - kernel.cols / 2 - 1, kernel.rows - kernel.rows / 2 - 1);
	int borderMode = cv::BORDER_CONSTANT;
	cv::filter2D(source, dest, img.depth(), kernel, anchor, 0, borderMode);
	return dest;
}

void GaussianSmooth(Mat src, Mat dst, float sigma)
{
	if (src.type() != CV_32FC1)
	{
		src.convertTo(src, CV_32FC1);
	}
	int kSize = static_cast<int>(round(2 * sigma) * 2) + 1;      // matlab: int kSize = round(2 * sigma) * 2 + 1;
	Mat K = cv::getGaussianKernel(kSize, sigma, CV_32FC1);       // 获取 Gaussian kernel, sigma是标准差
	K = K * K.t();
	Mat dstImg = Conv2(src, K);
	dstImg.convertTo(dst, CV_8U, 1);
}


//速度最快，但是特征点识别能力最弱
void ORBMatch(Mat img1, Mat img2, vector<KeyPoint>& keypoints1, vector<KeyPoint>& keypoints2,
	vector<DMatch>& bestMatches)
{
	const  float ratio = 0.8;
	//1. initinal
	//Ptr<FeatureDetector> detector = ORB::create();
	//Ptr<DescriptorExtractor> extractor = ORB::create();
	Ptr<ORB> orb = ORB::create();

	Mat descriptors1, descriptors2;
	vector<vector<DMatch>> matches12, matches21;
	Ptr<DescriptorMatcher> matcher = DescriptorMatcher::create("BruteForce-Hamming");

	// 2. detect features and extract the descriptors
	//detector->detect(img1, keypoints1);
	//extractor->compute(img1, keypoints1, descriptors1);
	//detector->detect(img2, keypoints2);
	//extractor->compute(img2, keypoints2, descriptors2);
	orb->detect(img1, keypoints1);
	orb->compute(img1, keypoints1, descriptors1);
	orb->detect(img2, keypoints2);
	orb->compute(img2, keypoints2, descriptors2);


#ifdef _DEBUG
	Mat imgKeypoints1, imgKeypoints2;
	drawKeypoints(img1, keypoints1, imgKeypoints1, Scalar::all(-1), DrawMatchesFlags::DEFAULT);
	drawKeypoints(img2, keypoints2, imgKeypoints2, Scalar::all(-1), DrawMatchesFlags::DEFAULT);
#endif

	if (!(keypoints1.size() > 0 && keypoints2.size() > 0))
	{
		return;
	}

	//3.Match the descriptors in two directions...
	matcher->knnMatch(descriptors1, descriptors2, matches12, 2);
	matcher->knnMatch(descriptors2, descriptors1, matches21, 2);

	//4. ratio test proposed by David Lowe's paper 0.8
	vector<DMatch> goodMatches1, goodMatches2;
	goodMatches1 = RatioTest(matches12, ratio);
	goodMatches2 = RatioTest(matches21, ratio);

	// Symmetric Test
	bestMatches = SymmetricTest(goodMatches1, goodMatches2);
#ifdef _DEBUG
	 cout << "Good matches1:" << goodMatches1.size() << endl;
	 cout << "Good matches2:" << goodMatches2.size() << endl;
	 cout << "Best matches:" << bestMatches.size() << endl;
#endif
}


void SURFMatch(Mat img1, Mat img2, vector<KeyPoint>& keypoints1, vector<KeyPoint>& keypoints2,
	           vector<DMatch>& bestMatches)
{
	// Step 1: Detect the keypoints using SURF Detector
	double minHessian = 400.0 /*800*/;
	//SurfFeatureDetector detector(minHessian);
	Ptr<SurfFeatureDetector> detector = SurfFeatureDetector::create();

	detector->setHessianThreshold(minHessian);

	detector->detect(img1, keypoints1);
	detector->detect(img2, keypoints2);

	if (!(keypoints1.size()>0 && keypoints2.size()>0))
	{
		return;
	}

	// Step 2: Calculate descriptors (feature vectors)
	//SurfDescriptorExtractor extractor;
	Ptr<SurfDescriptorExtractor> extractor = SurfFeatureDetector::create();
	Mat descriptors1, descriptors2;
	extractor->compute(img1, keypoints1, descriptors1);
	extractor->compute(img2, keypoints2, descriptors2);

#ifdef _DEBUG
	Mat imgKeypoints1, imgKeypoints2;
	drawKeypoints(img1, keypoints1, imgKeypoints1, Scalar::all(-1), DrawMatchesFlags::DEFAULT);
	drawKeypoints(img2, keypoints2, imgKeypoints2, Scalar::all(-1), DrawMatchesFlags::DEFAULT);
#endif

	// Step 3: Matching descriptor vectors using FLANN matcher
	FlannBasedMatcher matcher;
	//Ptr<DescriptorMatcher> matcher = DescriptorMatcher::create("FlannBased"); 
	//Ptr<DescriptorMatcher> matcher = DescriptorMatcher::create(DescriptorMatcher::FLANNBASED);
	// BFMatcher matcher(NORM_L2, true);

	vector<DMatch> matches1, matches2;
	matcher.match(descriptors1, descriptors2, matches1);
	matcher.match(descriptors2, descriptors1, matches2);

	vector<DMatch> goodMatches1, goodMatches2;
	goodMatches1 = FilterDistance(descriptors1, matches1);
	goodMatches2 = FilterDistance(descriptors2, matches2);

	for (int i = 0; i < goodMatches1.size(); i++)
	{
		DMatch temp1 = goodMatches1[i];
		for (int j = 0; j < goodMatches2.size(); j++)
		{
			DMatch temp2 = goodMatches2[j];
			if (temp1.queryIdx == temp2.trainIdx && temp2.queryIdx == temp1.trainIdx)
			{
				bestMatches.push_back(temp1);
				break;
			}
		}
	}

#ifdef _DEBUG
	// Draw only "good" matches
	Mat img_matches;
	cv::drawMatches(img1, keypoints1, img2, keypoints2, bestMatches, img_matches, Scalar::all(-1), Scalar::all(-1),
		vector<char>(), DrawMatchesFlags::NOT_DRAW_SINGLE_POINTS);
	//for (int i = 0; i < (int)bestMatches.size(); i++)
	//{
	//	printf("-- Good Match [%d] Keypoint 1: %d  -- Keypoint 2: %d  \n", i, bestMatches[i].queryIdx, bestMatches[i].trainIdx);
    //}
#endif 
}

// 对单服图像检测到的特征点较多，优于SURF, 但对同2幅图重复检测时匹配的特征点会跳动
void SIFTMatch(Mat img1, Mat img2, vector<KeyPoint>& keypoints1, vector<KeyPoint>& keypoints2,
	           vector<DMatch>& bestMatches)
{
	const float ratio = 0.8;
	Ptr<SiftFeatureDetector> detector = SiftFeatureDetector::create(0, 3, 0.04, 10, 1.2);
	Ptr<SiftDescriptorExtractor> extractor = SiftDescriptorExtractor::create();
	detector->detect(img1, keypoints1);
	detector->detect(img2, keypoints2);
#ifdef _DEBUG
cout << "keypoints1 的数目是：" << keypoints1.size() << ",  keypoints2 的数目是：" << keypoints2.size() << endl;
#endif
	Mat descriptors1, descriptors2;
	extractor->compute(img1, keypoints1, descriptors1);
	extractor->compute(img2, keypoints2, descriptors2);
#ifdef _DEBUG
	Mat imgKeypoints1, imgKeypoints2;
	drawKeypoints(img1, keypoints1, imgKeypoints1, Scalar::all(-1), DrawMatchesFlags::DEFAULT);
	drawKeypoints(img2, keypoints2, imgKeypoints2, Scalar::all(-1), DrawMatchesFlags::DEFAULT);
#endif
	// 当 keypoints1.size() == keypoints2.size() == 1 时，触发 OpenCV Error:
	// Assertion failed ((globalDescIdx>=0) && (globalDescIdx < size()))
	// in cv::DescriptorMatcher::DescriptorCollection::getLocalIdx,
	// .\opencv\sources\modules\features2d\src\matchers.cpp, line 488
	// 当 keypoints1.size() 与 keypoints2.size()任意一个等于0时，触发OpenCV Error:
	// Unsupported format or combination of formats (type=0) in cv::flann::buildIndex_,
	// .\opencv\sources\modules\flann\src\miniflann.cpp, line 315
	if ((keypoints2.size() <= 1) || (keypoints1.size() <= 1))
	{
		return;
	}
	vector<vector<DMatch>> matches12, matches21;
	Ptr<DescriptorMatcher> matcher = DescriptorMatcher::create("FlannBased");
	
	// knnMatch()的第二个参数对应的 keypoints.size()==0 时，FlannBasedMatcher::train() 报错
	// knnMatch()的第二个参数对应的 keypoints.size()==1 时, flannIndex->knnSearch() 报错
	matcher->knnMatch(descriptors1, descriptors2, matches12, 2);  
	matcher->knnMatch(descriptors2, descriptors1, matches21, 2);

	// ratio test proposed by David Lowe paper = 0.8
	vector<DMatch> goodMatches1, goodMatches2;
	goodMatches1 = RatioTest(matches12, ratio);
	goodMatches2 = RatioTest(matches21, ratio);

	bestMatches = SymmetricTest(goodMatches1, goodMatches2);
#ifdef _DEBUG
	Mat output;
	cv::drawMatches(img1, keypoints1, img2, keypoints2, bestMatches, output);
#endif 
}

//ratio test nearest/second nearest < ratio
vector<DMatch> RatioTest(vector<vector<DMatch>> matches12, double ratio)
{
	vector<DMatch> goodMatches;
	for (int i = 0; i < (int)matches12.size(); i++)
	{
		if (matches12[i].size() == 2)
		{
			if (matches12[i][0].distance < ratio * matches12[i][1].distance)
				goodMatches.push_back(matches12[i][0]);
		}
		else
		{
			break;
		}	
	}
	return goodMatches;
}

// Symmetric Test
vector<DMatch> SymmetricTest(vector<DMatch> goodMatches1, vector<DMatch> goodMatches2)
{
	vector<DMatch> betterMatches;
	for (int i = 0; i < goodMatches1.size(); i++)
	{
		for (int j = 0; j < goodMatches2.size(); j++)
		{
			if (goodMatches1[i].queryIdx == goodMatches2[j].trainIdx && goodMatches2[j].queryIdx == goodMatches1[i].trainIdx)
			{
				betterMatches.push_back(DMatch(goodMatches1[i].queryIdx, goodMatches1[i].trainIdx, goodMatches1[i].distance));
				break;
			}
		}
	}

	return betterMatches;
}

vector<DMatch> FilterDistance(Mat descriptors, vector<DMatch> matches)
{
	double max_dist = 0; double min_dist = 100;

	// Quick calculation of max and min distances between keypoints
	for (int i = 0; i < descriptors.rows; i++)
	{
		double dist = matches[i].distance;
		if (dist < min_dist)
			min_dist = dist;
		if (dist > max_dist)
			max_dist = dist;
	}

	// Draw only "good" matches (i.e. whose distance is less than 2*min_dist,
	// or a small arbitary value ( 0.02 ) in the event that min_dist is very small)
	// PS.- radiusMatch can also be used here.
	std::vector< DMatch > good_matches;

	for (int i = 0; i < descriptors.rows; i++)
	{
		if (matches[i].distance <= 5 * min_dist)
		{
			good_matches.push_back(matches[i]);
		}
	}
	return good_matches;
}