基于特征检测（SURF，SIFT方法）与特征匹配（Feature Matching）（FLANN方法）_surf.detectandcomput()返回值

本文主要讲述三个部分：
Feature extraction, Feature matching, Feature tracking.
另外还提到了透视变换：
Perspective transform.

内容的最后有我的完整代码实现。

特征获取（Feature extraction）

很多的低级特征，例如边，角，团，脊会比一个像素的灰度值所带有的信息多的多。在不同的应用，一些特征会比其它特征更加的有用。一旦想好我们想要的特征的构成，我们就要想办法在图片里找到我们想要的特征。

特征检测（Feature detection）

在图片里找到我们感兴趣的区域的过程就叫做特征检测。OpenCV中提供了多个特征检测算法：

• Harris corner detection: 角点时在所有方向像素变化剧烈的点，Harris和Stephens提出了检测这样区域的快速的方法。opencv:cv2.cornerHarris
• Shi-Tomasi corner detection:通常比Harris方法更优，他们查找N个最强的角点。opencv:cv2.goodFeaturesToTrack
• Scale-Invariant Feature Transform(SIFT):在图像大小改变时角点检测的效果就不好了，Lowe提出了一个描述图像里与角度大小无关的关键点的方法。在opencv3中，SIFT在contrib模块里，ubuntu环境安装opencv_contrib的方法见我写的教程。windows的话，opencv3.2中集成了contrib，可以在这里找到对应的whl文件通过pip安装即可。opencv2:cv2.SIFT，opencv3:cv2.xfeatures2d.SIFT_create()
• Speeded-Up Robust Features(SURF):SIFT是一个很好的方法，但是对于大部分应用来说，它不够快。SURF将SIFT中的Laplacian of a Gaussian（LOG）用一个方框滤波（box filter）代替。opencv2:cv2.SURF，opencv3:cv2.xfeatures2d.SURF_create()
• OpenCV支持很多的特征描述，包括Features fromAccelerated Segment Test (FAST), Binary Robust IndependentElementary Features (BRIEF), Oriented FAST，Rotated BRIEF(ORB)。

使用SURF在图片里检测特征

SURF算法可以粗略分成两个步骤：检测兴趣点，描述描述符。SURF依赖于Hessian角点检测方法对于兴趣点的探测，因此需要设置一个min_hessian的阈值。这个阈值决定了一个点要称为兴趣点，它对应的Hessian filter输出至少要有多大。大的值输出的数量比较少但是它们更为突出，相比之下输出较小的值虽然多但是不够突出（就是与普通差别不够大）。文中代码阈值设置为400：

    def extract_features(self):
        self.min_hessian = 400
        self.SURF = cv2.xfeatures2d.SURF_create(self.min_hessian)

特征和描述符只需要一步就能获得：

        #  detectAndCompute函数返回关键点和描述符，mask为None
        # 注意，书中的query和train图片和opencv官方的turorials里的是相反的
        key_query, desc_query = self.SURF.detectAndCompute(self.img_query, None)

通过以下函数就能简单的将关键点画出：
imgOut = cv2.drawKeypoints(self.img, self.key_train, None, (255, 0, 0), 4)
注意：在获得特征点后，要先检查一下特征点的数量:len(key_query)，以避免返回过多的特征点（太多则修改min_hessian）。

特征匹配

通过Fast Library for Approximate Nearest Neighbors(FLANN)方法将当前帧中像我们感兴趣的对象给找出来：
good_matches = self.matchfeatures(desc_query)
找到帧和帧之间的一致性的过程就是在一个描述符集合（询问集）中找另一个集合（相当于训练集）的最近邻。

通过FLANN方法进行特征匹配

可选的方法是利用近似k近邻算法去寻找一致性，FLANN方法比BF（Brute-Force）方法快的多：

    def matchfeatures(self, desc_frame):
        # 函数返回一个训练集和询问集的一致性列表
        matches = self.flann.knnMatch(self.desc_train, desc_frame, k=2)

用比值判别法（ratio test）删除离群点

检测出的匹配点可能有一些是错误正例（false positives）。

以为这里使用过的kNN匹配的k值为2（在训练集中找两个点），第一个匹配的是最近邻，第二个匹配的是次近邻。直觉上，一个正确的匹配会更接近第一个邻居。换句话说，一个不正确的匹配，两个邻居的距离是相似的。因此，我们可以通过查看二者距离的不同来评判距匹配程度的好坏。比值检测认为第一个匹配和第二个匹配的比值小于一个给定的值（一般是0.5），这里是0.7：

        # 丢弃坏的匹配
        good_matches = filter(lambda x: x[0].distance < 0.7*x[1].distance, matches)

通过cv2.drawMatchesKnn画出匹配的特征点，再将好的匹配返回：

        return good_matches

在复杂的环境中，FLANN算法不容易将对象混淆，而像素级算法则容易混淆。以下是书中的结果：

单应性估计

由于我们的对象是平面且固定的，所以我们就可以找到两幅图片特征点的单应性变换。得到单应性变换的矩阵后就可以计算对应的目标角点：

# 代码与书中不同，做了一些修改
def detectcorner_points(self, key_frame, good_matches):
        # 将所有好的匹配的对应点的坐标存储下来
        src_points = np.float32([self.key_train[m.trainIdx].pt for m in good_matches]).reshape(-1, 1, 2)
        dst_points = np.float32([keyQuery[m.queryIdx].pt for m in good_matches]).reshape(-1, 1, 2)
        H, mask= cv2.findHomography(src_points, dst_points, cv2.RANSAC, 5.0)
        matches_mask = mask.ravel().tolist()
        # 有了H单应性矩阵，我们可以查看源点被映射到query image中的位置
        self.sh_train = self.img_train.shape[:2] # rows, cols
        src_corners = np.float32([(0, 0), (self.sh_train[1], 0), (self.sh_train[1], self.sh_train[0]), (0,self.sh_train[0])]).reshape(-1, 1, 2)
        # perspectiveTransform返回点的列表
        dst_corners = cv2.perspectiveTransform(src_corners, H)
        dst_corners = map(tuple, dst_corners[0])
        # 将点向右移动img_train的宽度大小，方便我们同时显示两张图片
        dst_corners = [(np.int(dst_corners[i][0]+self.sh_train[1]), np.int(dst_corners[i][1])
        for i in range(0,len(dst_corners)):
            cv2.line(img_flann, dst_corners[i], dst_corners[(i+1) % 4],(0, 255, 0), 3)

结果图：

弯曲图片

我们可以将场景改变，使得看上去像正对着这本书。我们可以简单的将单应性矩阵取逆：
Hinv = cv2.linalg.inverse(H)
但是，这会将书左上角的点变成新图片的原点，书本左边和上面的部分都会被剪断。我们试图仅仅大概的将书本放在图片的中间。因此我们计算一个新的单应性矩阵。将场景点作为输入，输出的图片中的书本要和模板图片里的一样大：

   def warp_keypoints(self):
       dst_size = img_in.shape[:2]

将书本大小缩小到dst_size大小的1/2，同时还要移动1/4的距离：

        scale_row = 1./src_size[0]*dst_size[0]/2.
        bias_row = dst_size[0]/4.
        scale_col = 1./src_size[1]*dst_size[1]/2.
        bias_col = dst_size[1]/4.
        # 将每个点应用这样的变换
        src_points = [key_frame[good_matches[i].trainIdx].pt for i in range(len(good_matches))]
        dst_points = [self.key_train[good_matches[i].queryIdx].pt for i in range(len(good_matches))]
        dst_points = [[x*scale_row+bias_row, y*scale_col+bias_col] for x, y in dst_points]
        Hinv, = cv2.findHomography(np.array(srcpoints), np.array(dst_points), cv2.RANSAC)
        img_warp = cv2.warpPerspective(img_query, Hinv, dst_size)

特征跟踪

如何保证一个帧里找到的图在下一帧里再被找到。

在FearureMatching类的构造函数中，创建了一些记录的变量。主要的想法是从一帧跑到下一帧时要加强一些连贯性。因此我们抓取了大约每秒10帧的图，虽然上一帧里的图和下一帧变化并不会太大，但是也不能因此而把新的一帧里的一些离群点认为是正确的。为了解决这个问题，我们保存了我们没有找到合适结果的帧数量：self.num_frames_no_success，如果这个数量小于self.max_frames_no_success，我们将这些帧进行比较。如果大于阈值，我们就假定距离最后一次在帧中获取结果的时间已经过去了很久，这种情况下就不需要再帧中比较结果了。

早期的离群点检测与排除

我们可以将离群点的排除放在每次计算的步骤中，尽量保证获取好的匹配。

def match(self, frame):
        img_query = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
        sh_query = img_query.shape[:2] # rows,cols
        # 获得好的matches
        key_query, desc_query = self._extract_features(img_query)
        good_matches = self._match_features(descQuery)
        # 为了让RANSAC方法可以尽快工作，至少需要4个好的匹配，否则视为匹配失败
        if len(good_matches) < 4:
            self.num_frames_no_success=self.num_frames_no_success + 1
            return False, frame
        # 在query_image中找到对应的角点
        dst_corners = self._detect_corner_points(key_query, good_matches)
        # 如果这些点位置距离图片内太远（至少20像素），那么意味着我们没有找到我们感兴趣
        # 的目标或者说是目标没有完整的出现在图片内，对于这两种情况，我们都视为False
        if np.any(filter(lambda x: x[0] < -20 or x[1] < -20
            or x[0] > sh_query[1] + 20 or x[1] > sh_query[0] + 20, dst_corners)):
            self.num_frames_no_success =
            self.num_frames_no_success + 1
            return False, frame
        # 如果4个角点没有围出一个合理的四边形，意味着我们可能没有找到我们的目标。
        # 计算面积
        area = 0
        for i in range(0, 4):
            next_i = (i + 1) % 4
            area = area + (dst_corners[i][0]*dst_corners[next_i]
            [1]- dst_corners[i][1]*dst_corners[next_i][0])/2.
        # 如果面积太大或太小，将它排除
        if area < np.prod(sh_query)/16. or area > np.prod(sh_query)/2.:
            self.num_frames_no_success=self.num_frames_no_success + 1
            return False, frame
        # 如果我们此时发现的单应性矩阵和上一次发现的单应性矩阵变化太大，意味着我们可能找到了
        # 另一个对象，这种情况我们丢弃这个帧并返回False
        np.linalg.norm(Hinv – self.last_hinv)
        # 这里要用到self.max_frames_no_success的，作用就是距离上一次发现的单应性矩阵
        # 不能太久时间，如果时间过长的话，完全可以将上一次的hinv抛弃，使用当前计算得到
        # 的Hinv
        recent = self.num_frames_no_success < self.max_frames_no_success
        similar = np.linalg.norm(Hinv - self.last_hinv) < self.max_error_hinv
        if recent and not similar:
            self.num_frames_no_success = self.num_frames_no_success + 1
            return False, frame
        self.num_frames_no_success = 0
        self.last_hinv = Hinv   
        img_out = cv2.warpPerspective(img_query, Hinv, dst_size)
        img_out = cv2.cvtColor(img_out, cv2.COLOR_GRAY2RGB)
            return True, imgOut

书中运行起来的效果如下：

warp image是变化不大的，近乎静止：

如果算出了错误的单应性矩阵，因为self.max_frames_no_success的设置也能很快的恢复，使得算法精确且有效率。

整个程序的流程总结如下：

1. 抽取每个视频帧的感兴趣的特征：_extract_features
2. 将模板和视频帧的特征点进行匹配，没有匹配到则跳过：_match_features
3. 检测视频帧中的对应模板图片的角点，如果一些重要的角点在该帧之外则跳过：_detect_corner_points
4. 计算四个角点构成的四边形的面积，如果太小或太大则跳过
5. 在当前帧大概画出模板图片的角点
6. 计算将当前物体从当前帧移动到前额平行平面的透视变换，如果当前的结果与之前帧的结果差别很大，则跳过：_warp_keypoints
7. 扭曲当前帧使得目标物体出现在中央并且正面朝上

实际运行代码及结果

书里的代码组合起来有些错误，下面是我修改过的整个FeatureMatching类的实现：

（注意几个变动：

• 其中的train image和query image我是按照opencv的官方python教程来设置的，本书的train image和query image与之是相反的；
• 其中求四边形的面积我改用的是行列式的方式；
• 书中使用的是SURF方法求特征点，但是就我的两张图片SURF得不到好的结果，故改用了SIFT方法；
  ）

import cv2
import numpy as np
from matplotlib import pyplot as plt
 
class FeatureMatching:
    # 官方教程的目标图片是query image
    def __init__(self, query_image='data/query.jpg'):
        # 创建SURF探测器，并设置Hessian阈值，由于效果不好，我改成了SIFT方法
        # self.min_hessian = 400（surf方法使用）
        # self.surf = cv2.xfeatures2d.SURF_create(min_hessian)
        self.sift = cv2.xfeatures2d.SIFT_create()
        self.img_query = cv2.imread(query_image, 0)
        # 读取一个目标模板
        if self.img_query is None:
            print("Could not find train image " + query_image)
            raise SystemExit
        self.shape_query = self.img_query.shape[:2]  # 注意，rows，cols，对应的是y和x，后面的角点坐标的x,y要搞清楚
        #  detectAndCompute函数返回关键点和描述符
        self.key_query, self.desc_query = self.sift.detectAndCompute(self.img_query, None)
        # 设置FLANN对象
        FLANN_INDEX_KDTREE = 0
        index_params = dict(algorithm=FLANN_INDEX_KDTREE, trees=5)
        search_params = dict(checks=50)
        self.flann = cv2.FlannBasedMatcher(index_params, search_params)
        # 保存最后一次计算的单应矩阵
        self.last_hinv = np.zeros((3, 3))
        # 保存没有找到目标的帧的数量
        self.num_frames_no_success = 0
        # 最大连续没有找到目标的帧的次数
        self.max_frames_no_success = 5
        self.max_error_hinv = 50.
        # 防止第一次检测到时由于单应矩阵变化过大而退出
        self.first_frame = True
 
    def _extract_features(self, frame):
        # self.min_hessian = 400
        # sift = cv2.xfeatures2d.SURF_create(self.min_hessian)
        sift = cv2.xfeatures2d.SIFT_create()
        #  detectAndCompute函数返回关键点和描述符，mask为None
        key_train, desc_train = sift.detectAndCompute(frame, None)
        return key_train, desc_train
 
    def _match_features(self, desc_frame):
        # 函数返回一个训练集和询问集的一致性列表
        matches = self.flann.knnMatch(self.desc_query, desc_frame, k=2)
        # 丢弃坏的匹配
        good_matches = []
        # matches中每个元素是两个对象，分别是与测试的点距离最近的两个点的信息
        # 留下距离更近的那个匹配点
        for m, n in matches:
            if m.distance < 0.7 * n.distance:
                good_matches.append(m)
        return good_matches
 
    def _detect_corner_points(self, key_frame, good_matches):
        # 将所有好的匹配的对应点的坐标存储下来
        src_points = np.float32([self.key_query[m.queryIdx].pt for m in good_matches]).reshape(-1, 1, 2)
        dst_points = np.float32([key_frame[m.trainIdx].pt for m in good_matches]).reshape(-1, 1, 2)
        H, mask = cv2.findHomography(src_points, dst_points, cv2.RANSAC, 5.0)
        matchesMask = mask.ravel().tolist()
        # 有了H单应性矩阵，我们可以查看源点被映射到img_query中的位置
        # src_corners = np.float32([(0, 0), (self.shape_train[1], 0), (self.shape_train[1], self.shape_train[0]),
        #                           (0, self.shape_train[0])]).reshape(-1, 1, 2)
        h, w = self.img_query.shape[:2]
        src_corners = np.float32([[0, 0], [0, h - 1], [w - 1, h - 1], [w - 1, 0]]).reshape(-1, 1, 2)
        # perspectiveTransform返回点的列表
        dst_corners = cv2.perspectiveTransform(src_corners, H)
        return dst_corners, H, matchesMask
 
    def _center_keypoints(self, frame, key_frame, good_matches):
        dst_size = frame.shape[:2]
        # 将图片的对象大小缩小到query image的1/2（书里是train image，和官方命名相反而已）
        scale_row = 1. / self.shape_query[0] * dst_size[0] / 2.
        bias_row = dst_size[0] / 4.
        scale_col = 1. / self.shape_query[1] * dst_size[1] / 2.
        bias_col = dst_size[1] / 4.
        # 将每个点应用这样的变换
        src_points = [self.key_query[m.queryIdx].pt for m in good_matches]
        dst_points = [key_frame[m.trainIdx].pt for m in good_matches]
        dst_points = [[x * scale_row + bias_row, y * scale_col + bias_col] for x, y in dst_points]
        Hinv, _ = cv2.findHomography(np.array(src_points), np.array(dst_points), cv2.RANSAC, 5.0)
        img_center = cv2.warpPerspective(frame, Hinv, dst_size, flags=2)
        return img_center
 
    def _frontal_keypoints(self, frame, H):
        Hinv = np.linalg.inv(H)
        dst_size = frame.shape[:2]
        img_front = cv2.warpPerspective(frame, Hinv, dst_size, flags=2)
        return img_front
 
    def match(self, frame):
        img_train = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
        cv2.waitKey(0)
        shape_train = img_train.shape[:2]  # rows,cols
 
        # 获得好的matches
        key_train, desc_train = self._extract_features(img_train)
        good_matches = self._match_features(desc_train)
 
        # 为了让RANSAC方法可以尽快工作，至少需要4个好的匹配，否则视为匹配失败
        if len(good_matches) < 4:
            self.num_frames_no_success += 1
            return False, frame
        # 画出匹配的点
        img_match = cv2.drawMatchesKnn(self.img_query, self.key_query, img_train, key_train, [good_matches], None,
                                       flags=2)
        plt.imshow(img_match), plt.show()
 
        # 在query_image中找到对应的角点
        dst_corners, Hinv, matchesMask = self._detect_corner_points(key_train, good_matches)
        # 如果这些点位置距离图片内太远（至少20像素），那么意味着我们没有找到我们感兴趣
        # 的目标或者说是目标没有完整的出现在图片内，对于这两种情况，我们都视为False
        dst_ravel = dst_corners.ravel()
        if (dst_ravel > shape_train[0] + 20).any() and (dst_ravel > -20).any() \
                and (dst_ravel > shape_train[1] + 20).any():
            self.num_frames_no_success += 1
            return False, frame
 
        # 如果4个角点没有围出一个合理的四边形，意味着我们可能没有找到我们的目标。
        # 通过行列式计算四边形面积
        area = 0.
        for i in range(0, 4):
            D = np.array([[1., 1., 1.],
                          [dst_corners[i][0][0], dst_corners[(i + 1) % 4][0][0], dst_corners[(i + 2) % 4][0][0]],
                          [dst_corners[i][0][1], dst_corners[(i + 1) % 4][0][1], dst_corners[(i + 2) % 4][0][1]]])
            area += abs(np.linalg.det(D)) / 2.
        area /= 2.
        # 以下注释部分是书中的计算方式，我使用时是错误的
        # for i in range(0, 4):
        #     next_i = (i + 1) % 4
        #     print(dst_corners[i][0][0])
        #     print(dst_corners[i][0][1])
        #     area += (dst_corners[i][0][0] * dst_corners[next_i][0][1] - dst_corners[i][0][1] * dst_corners[next_i][0][
        #         0]) / 2.
        # 如果面积太大或太小，将它排除
        if area < np.prod(shape_train) / 16. or area > np.prod(shape_train) / 2.:
            self.num_frames_no_success += 1
            return False, frame
 
        # 如果我们此时发现的单应性矩阵和上一次发现的单应性矩阵变化太大，意味着我们可能找到了
        # 另一个对象，这种情况我们丢弃这个帧并返回False
        # 这里要用到self.max_frames_no_success的，作用就是距离上一次发现的单应性矩阵
        # 不能太久时间，如果时间过长的话，完全可以将上一次的hinv抛弃，使用当前计算得到
        # 的Hinv
        recent = self.num_frames_no_success < self.max_frames_no_success
        similar = np.linalg.norm(Hinv - self.last_hinv) < self.max_error_hinv
        if recent and not similar and not self.first_frame:
            self.num_frames_no_success += self.num_frames_no_success
            return False, frame
        # 第一次检测标志置否
        self.first_frame = False
        self.num_frames_no_success = 0
        self.last_hinv = Hinv
 
        draw_params = dict(matchColor=(0, 255, 0),  # draw matches in green color
                           singlePointColor=None,
                           matchesMask=matchesMask,  # draw only inliers
                           flags=2)
        img_dst = cv2.polylines(img_train, [np.int32(dst_corners)], True, (0, 255, 255), 5, cv2.LINE_AA)
        img_dst = cv2.drawMatches(self.img_query, self.key_query, img_dst, key_train, good_matches, None,
                                  **draw_params)
        plt.imshow(img_dst)
        plt.show()
 
        img_center = self._center_keypoints(frame, key_train, good_matches)
        plt.imshow(img_center)
        plt.show()
 
        # 转换成正面视角
        img_front = self._frontal_keypoints(frame, Hinv)
        plt.imshow(img_front)
        plt.show()
        return True, img_dst

测试：

import cv2
from feature_matching import FeatureMatching
 
img_train = cv2.imread('data/BM_left1.jpg')
 
matching = FeatureMatching(query_image='data/query.jpg')
flag = matching.match(img_train)