SIFT特征（五）

经过前面四个章节的介绍，我们已经知道了如何从一幅图像中提取稳定的SIFT特征。接下来，我们将针对OpenCV中提供的FeatureDetector函数来说明具体如何实现SIFT特征检测，以及如何使用检测到的SIFT特征点进行目标识别。

对于目标检测问题，通常是给定模板图像(如book.jpg)（数据库，可包含一张或多张模板图像提取的SIFT特征点），希望从场景图片(如scene.jpg)（待检测场景）中检测出模板物体的位置及相应变换矩阵。场景中包含的模板物体通常由于相机参数的改变，光照、噪声、遮挡等因素的影响产生改变。而SIFT特征点对于这些改变恰好具有较好的稳定性。

基于SIFT特征点进行物体检测的基本步骤包括，
a. 分别提取模板图像和场景图像的SIFT特征点，及每个特征点生成的128维度的特征向量。
b. 特征点匹配。 对于场景图像提取的每个特征点，在模板图像中寻找最优匹配和次优匹配特征点，如果最优匹配与次优匹配的比值过大，则认为是错误匹配对。
c. 使用RANSAC计算两个图像平面变换的单应矩阵，筛选不符合当前变换的“局外”匹配点（outliers）。

6.1 特征点匹配
关于如何提取场景图像以模板图像的SIFT特征点及特征向量，这里不再赘述。

对于场景图像中检测的特征点，我们希望从模板图像中寻找可能的匹配点，可以通过特征向量的欧式距离计算最近邻的方法获得。但是对于场景中的背景特征点，极有可能是无法在数据库中找到对应点的，或者部分特征点有可能被错误匹配，那么就必须确定一种筛选方法将错误匹配点删除。

最简单的思路是确定一个距离阈值，凡是大于该阈值的匹配认为是错误匹配，但是所有特征点选择统一的阈值度量并不鲁棒。更优的决策是对于每个特征点，计算在数据库中的最优匹配和次优匹配（欧氏距离最小及排名第二小的匹配点）距离。通常情况下，正确匹配点的最优正确匹配应明显地区别于其他错误匹配，而对于错误匹配点，由于描述特征的多维度，在计算距离排序时可能存在前若干个距离相近且较大的错误匹配。所以可以通过计算每个特征点的最优匹配和次优匹配的距离比值，如果该比值小于某个给定的阈值（最优匹配距离要明显大于次优匹配距离），则认为该匹配点是被正确匹配的，否则被认为是错误匹配的。（详见代码part 2部分）

6.2 RANSAC计算单应矩阵筛选特征点
6.2.1 单应矩阵
在计算机视觉中，平面的单应性被定义为从一个二维平面到另一个二维平面的投影映射，包括平移、放射、尺度变换，试图找到一个投影矩阵，该投影矩阵可以将一个平面的所有点，匹配到另一个平面中。例如，照相机标定可以采用棋盘格平面来辅助完成，在这个过程中，从棋盘格平面到摄像机投影平面的映射反应的是平面的单应性；又比如，空间平面在两个摄像机下的投影图像之间具有点一一对应的关系（见下图），这种对应关系也可以通过单应矩阵来进行表达。

具体描述为，设π 是不通过两摄像机任一光心的空间平面，它在两个摄像机下的图像分别记为I和I′，如上图所示。令X 是平面π 上的任意一点，它在两个摄像机下的像分别为m和m′。空间平面π 到两个像平面之间存在两个单应矩阵H1 和H2，使得m=H1X，m′=H2X，这两个单应矩阵H1和H2 实现了空间平面π 到对应像平面之间的一一变换。同时，空间平面中的任一点投影到两投影面上的两点之间也存在单应关系，他们之间的变换可以表达为m′=Hm，其中H=H2H−11，即平面在两个摄像机下的投影图像之间存在一个二维映射变换，矩阵H 实现了从第一个平面到第二个平面的一一变换。在我们的检测例子中，模板图像和场景图像中的目标区域同样存在着这种单应关系，我们的目标正是通过计算两者的单应矩阵来排除错误匹配点。

单应矩阵H 是大小为3*3的矩阵，对于左右像上的对应点的齐次坐标表示为p1=[x1,y1,1]T和p2=[x2,y2,1]T ，单应矩阵可以将p1变换为点p2，即，

[x2,y2,1]T=H[x1,y1,1]T=⎛⎝⎜h11h21h31h12h22h32h13h231⎞⎠⎟[x1,y1,1]T

注意，这里表示的是齐次坐标，其对应在图像上的坐标分别为[x1,y1]T 和[x2,y2]T。关于其次坐标的表达，可以参考阅读博文，阐述的比较清楚。http://blog.csdn.net/wqvbjhc/article/details/6003112

注意，单应矩阵的最后一个元素为常数值1，原因分析是，对于两点p1 和p2 之间的单应矩阵H ，如果它的每个元素同乘以一个常数a，得到单应矩阵aH和H的作用是相同的，因为它无非是将齐次坐标点由p1变换从 p2变为ap2，对应到图像中的点都同为一点。所以，我们可将单应矩阵的所有元素同除以最后一个元素，进行归一化处理。

可以看出，单应矩阵具有8个未知量，因此需要8个方程求解。而每一对点提供两个方程，所以至少需要4对点求解。每一对匹配点提供的方程是，

x2=x1h11+y1h12+h13x1h31+y1h32+1

y2=x1h21+y1h22+h23x1h31+y1h32+1

将这个方程重新组织一下，得到等价的矩阵形式表达为，

Au=v

其中，

u=(x10y10100x10y101−x1x2−x1y2−x2y1−y1y2)

u=(h11h12h13h21h22h23h31h32)T

v=(x2y2)T

取四对不共线匹配点对，则可以将单应矩阵求解出来，且存在唯一解，而如果点对大于四对，可以通过最小二乘或SVD分解求解H。
在opencv中，提供求解两个点集之间的单应变换矩阵函数findhomography，函数原型为（具体参数描述见文档http://docs.opencv.org/2.4/modules/calib3d/doc/camera_calibration_and_3d_reconstruction.html?highlight=findhomography#cv2.findHomography）

Mat findHomography(InputArray srcPoints, InputArray dstPoints, int method=0, double ransacReprojThreshold=3, OutputArray mask=noArray() )

返回单应矩阵H，大小为3∗3；
1.srcPoints和dstPoints分别表示两个对应点集，大小为N∗2或者齐次坐标表示大小N∗3；
2.method可以选择基于RANSAC方法（method=CV_RANSAC），RANSAC通过多次随机选择四个匹配点来计算单应矩阵，最终返回的单应矩阵所对应误差最小，误差表示为，error=dis(dstPoints,H∗srcPoints)；
3.ransacRojThreshold用在RANSAC方法中，表示被分类为“局内点”的阈值，即满足error=dis(dstPoints,H∗srcPoints)≤ransacRojThreshold的点对被认为“局内点”，否则认为是“局外点”。关于RANSAC方法参考博文 http://grunt1223.iteye.com/blog/961063
4.status表示掩码，长度为输入匹配点对数量。如果为1表示对应点对为“局内点”，否则认为是“局外点”，即错误匹配点。我们可以根据status返回的信息判断匹配点保留与否。

单应矩阵参考链接：http://www.zhihu.com/question/23310855

6.3 C++源代码

#include "opencv2/highgui/highgui.hpp"#include "opencv2/imgproc/imgproc.hpp"#include "opencv2/nonfree/nonfree.hpp"#include "opencv2/nonfree/features2d.hpp"#include <opencv2/calib3d/calib3d.hpp>#include <iostream>using namespace cv;using namespace std;int main(int argc, char* argv[]){    //Part 1. 检测关键点，提取特征向量    initModule_nonfree();//初始化    Ptr<FeatureDetector> detector = FeatureDetector::create("SIFT");//create sift feature dector    Ptr<DescriptorExtractor> descriptor_extractor = DescriptorExtractor::create("SIFT");//create feature vector extractor    Ptr<DescriptorMatcher> descriptor_matcher = DescriptorMatcher::create("BruteForce");//create feature matcher    if (detector.empty() || descriptor_extractor.empty())        return - 1;    //读取图片    Mat img1 = imread("scene.jpg");    Mat img2 = imread("book.jpg");    //检测keypoints    vector<KeyPoint> keypoint1, keypoint2;    detector->detect(img1, keypoint1);    detector->detect(img2, keypoint2);    //根据keypoints计算128维度描述符    Mat descriptor1, descriptor2;    descriptor_extractor->compute(img1, keypoint1, descriptor1);    descriptor_extractor->compute(img2, keypoint2, descriptor2);    //绘制keypoints    Mat img_keypoints1, img_keypoints2;    drawKeypoints(img1, keypoint1, img_keypoints1);    drawKeypoints(img2, keypoint2, img_keypoints2);    imshow("keypoints1", img_keypoints1);    imshow("keypoints2", img_keypoints2);    //Part 2. 根据sift描述符匹配，并初步筛选    vector<vector<DMatch>> knn_match;    descriptor_matcher->knnMatch(descriptor1, descriptor2, knn_match, 2);//返回最优匹配和次优匹配    //正确匹配保证最优匹配和次优匹配之间具有较大差距，如果两者之间差距较小，则认为是错误匹配    const float minRatio = 1.f / 1.5f;    vector<DMatch> matches;    for (int i = 0; i < knn_match.size(); i++){        const DMatch& bestMatch = knn_match[i][0];        const DMatch& betterMatch = knn_match[i][1];        float r = bestMatch.distance / betterMatch.distance;        if (r < minRatio)            matches.push_back(bestMatch);    }    //绘制匹配结果    Mat img_knnMatches;    drawMatches(img1, keypoint1, img2, keypoint2, matches, img_knnMatches);    imshow("KNN matches", img_knnMatches);    // Part 3:RANSAC计算单应变换矩阵排除错误匹配点对    vector<KeyPoint> alignedKps1, alignedKps2;    for (int i = 0; i < matches.size(); i++){        alignedKps1.push_back(keypoint1[matches[i].queryIdx]);        alignedKps2.push_back(keypoint2[matches[i].trainIdx]);    }    vector<Point2f> ps1, ps2;    for (int i = 0; i < alignedKps1.size(); i++)        ps1.push_back(alignedKps1[i].pt);    for (int i = 0; i < alignedKps2.size(); i++)        ps2.push_back(alignedKps2[i].pt);    if (matches.size() < 4){        printf("Matches number is not enough");        return -1;    }    Mat status = Mat::zeros(matches.size(), 1, CV_8UC1);    Mat H = findHomography(ps2, ps1, CV_FM_RANSAC, 3, status);    vector<DMatch> RANSAC_matches;    uchar *status_p;    vector<DMatch>::const_iterator it_match = matches.begin();    for (int i = 0; i < matches.size(); i++){        status_p = status.ptr<uchar>(i);        if (*status_p)            RANSAC_matches.push_back(it_match[i]);    }    //绘制RANSAC筛选后的匹配结果    Mat img_RANSAC_matches;    drawMatches(img1, keypoint1, img2, keypoint2, RANSAC_matches, img_RANSAC_matches);    imshow("RANSAC_matches", img_RANSAC_matches);    //检测目标物体在场景图像中的位置    Mat img_detect = img1.clone();;    vector<Point2f> model_corners(4);//存储模板图像的四个角点    vector<Point2f> scene_corners(4);    //Part 4:确定模板在场景中的位置    model_corners[0] = cvPoint(0, 0);//左上角为原点(0,0)；x轴指向→；y轴指向↓    model_corners[1] = cvPoint(img2.cols, 0);    model_corners[2] = cvPoint(img2.cols, img2.rows);    model_corners[3] = cvPoint(0, img2.rows);    perspectiveTransform(model_corners, scene_corners, H);    line(img_detect, scene_corners[0], scene_corners[1], Scalar(0, 255, 0), 4);    line(img_detect, scene_corners[1], scene_corners[2], Scalar(0, 255, 0), 4);    line(img_detect, scene_corners[2], scene_corners[3], Scalar(0, 255, 0), 4);    line(img_detect, scene_corners[3], scene_corners[0], Scalar(0, 255, 0), 4);    imshow("Object detection", img_detect);    cvWaitKey(0);    return 0;}

6.3 结果展示
scene.jpg/ book.jpg SIFT特征点

KNN-match 筛选后匹配关系

RANSAC 筛选后匹配关系

模板位置检测，其中绿色框框出的表示被检测物体在场景中的位置

根据knn最优/次优比例筛选后，特征点匹配数目将至41对，RANSAC求解单应矩阵筛选特征点对后剩余38对特征点对。