关于 VLFeat 的说明可参考：[1] - VLFeat指南——SIFT检测器和描述子 - 2014.11.20[2] - 图像检索(1): 再论SIFT-基于vlfeat实现 - 2018.06.13VLFeat 的参数说明可见：https://www.vlfeat.org/api/sift_8h.html1. 命令行实现1.1. 下载编译vlfeat下载： https://www.vlfeat.org/download.htmlGithub：https://www.github.com/vlfeat/vlfeat.git由于 VLFeat 0.9.21 编译后运行 sift 输出的结果为空，这里采用 VLFeat 0.9.20.wget https://www.vlfeat.org/download/vlfeat-0.9.21.tar.gz tar -xzvf vlfeat-0.9.21.tar.gz cd vlfeat-0.9.20 make即可以在 vlfeat-0.9.20/bin 路径生成 sift 工具.1.2. Python 调用命令行Vlfeat sift：https://www.vlfeat.org/man/sift.htmlimport os from PIL import Image def vlfeat_sift_extract(img_file, params="--edge-thresh 10 --peak-thresh 1"): img_pil = Image.open(img_file).convert('L') gray_imgfile = os.path.basename(img_file) + '.pgm' img_pil.save(gray_imgfile) # sift_file = os.path.basename(img_file) + '.sift' cmd_str = "vlfeat-0.9.20/bin/glnxa64/sift {} --output={} {}".format(gray_imgfile, sift_file, params) msg = os.system(cmd_str) print('[INFO]msg: ', msg)2. Cyvlfeat 实现2.1. 安装Github - menpo/cyvlfeatConda 安装(推荐)：#https://anaconda.org/menpo/cyvlfeat conda install -c conda-forge cyvlfeat # 或pip #pip install cyvlfeat (调用时可能会出错)2.2. 示例import numpy as np from cyvlfeat.sift import sift as vlsift def vlfeat_sift_extract(img_file): img_cv2 = cv2.imread(img_file) gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY).astype(np.float32) kp, des = vlsift(gray, peak_thresh=5.0, compute_descriptor=True, float_descriptors=True) sift_vector = kp + des # sift_file = os.path.basename(img_file) + '.sift' np.savetxt(sift_file, np.array(sift_vector), fmt='%0.8f')3. SIFT 可视化import numpy as np import matplotlib.pyplot as plt from PIL import Image def vis_sifts(imgfile, locs, circle=True): def draw_circle(c,r): t = np.arange(0,1.01,.01)*2*np.pi x = r*np.cos(t) + c[0] y = r*np.sin(t) + c[1] plt.plot(x,y,'b',linewidth=2) plt.figure(figsize=(10, 8)) img_pil = Image.open(imgfile) plt.imshow(img_pil) if circle: for p in locs: draw_circle(p[:2],p[2]) else: plt.plot(locs[:,0],locs[:,1],'ob') plt.axis('off') plt.title('{} sift features'.format(len(locs))) plt.show()参考：vlfeat python
import cv2 # def match_snn(desc1, desc2, snn_th = 0.8): index_params = dict(algorithm=1, trees=4) search_params = dict(checks=128) # or pass empty dictionary matcher = cv2.FlannBasedMatcher(index_params, search_params) matches = matcher.knnMatch(desc1.astype(np.float32), desc2.astype(np.float32), k=2) #ratio test good_matches = [] for m,n in matches: if m.distance < snn_th * n.distance: good_matches.append(m) return good_matches取图像1中的某个关键点，并找出其与图像2中距离最近的前两个关键点，在这两个关键点中，若最近的距离除以次近的距离小于某个阈值，则接受这一对匹配点. 由于特征空间的高维性，相近的距离可能有大量的错误匹配，从而它的ratio值比较高. Lowe推荐ratio的取值为0.8，但通过对大量存在尺度、旋转和亮度变化的两幅图片进行匹配，实验结果表明ratio取值在0. 4~0. 6为最佳，小于0. 4会造成非常少的匹配点，大于0. 6会造成大量的错误匹配点.
From: https://github.com/ducha-aiki/affnet/blob/master/examples/SIFT-AffNet-HardNet-kornia-matching.ipynbimport cv2 import numpy as np def extract_sift_keypoints_upright(img, n_feat = 5000): sift = cv2.xfeatures2d.SIFT_create(2 * n_feat, contrastThreshold=-10000, edgeThreshold=-10000) keypoints = sift.detect(img, None) response = np.array([kp.response for kp in keypoints]) respSort = np.argsort(response)[::-1] kpts = [cv2.KeyPoint(keypoints[i].pt[0], keypoints[i].pt[1], keypoints[i].size, 0) for i in respSort] kpts_unique = [] for x in kpts: if x not in kpts_unique: kpts_unique.append(x) return kpts_unique[:n_feat] # def match_snn(desc1, desc2, snn_th = 0.8): index_params = dict(algorithm=1, trees=4) search_params = dict(checks=128) # or pass empty dictionary matcher = cv2.FlannBasedMatcher(index_params, search_params) matches = matcher.knnMatch(desc1.astype(np.float32), desc2.astype(np.float32), k=2) good_matches = [] for m,n in matches: if m.distance < snn_th * n.distance: good_matches.append(m) return good_matches
参考：Feature Matching + Homography to find Objects1. 实现如：import numpy as np import cv2 from matplotlib import pyplot as plt MIN_MATCH_COUNT = 10 img1 = cv2.imread('box.png',0) # queryImage img2 = cv2.imread('box_in_scene.png',0) # trainImage # Initiate SIFT detector sift = cv2.SIFT() # find the keypoints and descriptors with SIFT kp1, des1 = sift.detectAndCompute(img1,None) kp2, des2 = sift.detectAndCompute(img2,None) FLANN_INDEX_KDTREE = 0 index_params = dict(algorithm = FLANN_INDEX_KDTREE, trees = 5) search_params = dict(checks = 50) flann = cv2.FlannBasedMatcher(index_params, search_params) matches = flann.knnMatch(des1,des2,k=2) # store all the good matches as per Lowe's ratio test. good = [] for m,n in matches: if m.distance < 0.7*n.distance: good.append(m) # if len(good)>MIN_MATCH_COUNT: src_pts = np.float32([ kp1[m.queryIdx].pt for m in good ]).reshape(-1,1,2) dst_pts = np.float32([ kp2[m.trainIdx].pt for m in good ]).reshape(-1,1,2) #M: 3x3 变换矩阵. M, mask = cv2.findHomography(src_pts, dst_pts, cv2.RANSAC,5.0) matchesMask = mask.ravel().tolist() h,w,_ = img1.shape pts = np.float32([ [0,0],[0,h-1],[w-1,h-1],[w-1,0] ]).reshape(-1,1,2) dst = cv2.perspectiveTransform(pts,M) img2 = cv2.polylines(img2,[np.int32(dst)],True,255,3, cv2.LINE_AA) else: print "Not enough matches are found - %d/%d" % (len(good),MIN_MATCH_COUNT) matchesMask = None # draw_params = dict(matchColor = (0,255,0), # draw matches in green color singlePointColor = None, matchesMask = matchesMask, # draw only inliers flags = 2) img3 = cv2.drawMatches(img1,kp1,img2,kp2,good,None,**draw_params) plt.imshow(img3, 'gray'),plt.show()如图：2. 参考[1] - 在OpenCV3中进行特征点筛选与优化-Python版 - 2018.01.12[2] - 为何RANSAC算法提纯SIFT算法后仍存在误匹配对？- 知乎[3] - SIFT图像匹配及其python实现 - 2020.07.09[4] - 【Python计算机视觉编程】第二章 SIFT特征提取与检索及RANSAC实例 - 2020.03.08
1. BF 匹配BF 匹配，Brute-Force Matcher，暴力匹配. 其原理比较简单，首先从集合A中选择一个特征的描述子，然后与集合B中所有的其他特征计算某种相似度，进行匹配，并返回最接近的项.OpenCV 中，首先使用 cv2.BFMatcher()创建 BFMatcher 实例，其包含两个可选参数：normType 和 crossCheck.class BFMatcher(DescriptorMatcher): def create(self, normType=None, crossCheck=None): """ normType：如 NORM_L1, NORM_L2, NORM_HAMMING, NORM_HAMMING2. NORM_L1 和 NORM_L2 更适用于 SIFT 和 SURF 描述子; NORM_HAMMING 和 ORB、BRISK、BRIEF 一起使用； NORM_HAMMING2 用于 WTA_K==3或4 的 ORB 描述子. crossCheck：默认为 False，其寻找每个查询描述子的 k 个最近邻. 若值为 True，则 knnMatch() 算法 k=1，仅返回(i, j)的匹配结果， 即集合A中的第 i 个描述子在集合B中的第 j 个描述子是最佳匹配. 也就是说，两个集合中的两个描述子特征是互相匹配的. 其提供连续的结果. 当有足够的匹配项时，其往往能够找到最佳的匹配结果. """ pass实例化 BFMatcher 后，两个重要的方法是 BFMatcher.match() 和 BFMatcher.knnMatch(). 前者仅返回最佳匹配结果，后者返回 k 个最佳匹配结果.类似于 cv2.drawKeypoints() 画出关键点；cv2.drawMatches() 能够画出匹配结果. 其水平的堆叠两幅图像，并画出第一张图像到第二张图像的点连线，以表征最佳匹配. cv2.drawMatchKnn()能够画出 k 个最佳匹配.2. 基于ORB描述子的BF匹配示例import numpy as np import cv2 from matplotlib import pyplot as plt img1 = cv2.imread('box.png',0) # queryImage img2 = cv2.imread('box_in_scene.png',0) # trainImage # Initiate ORB detector orb = cv2.ORB() # find the keypoints and descriptors with SIFT kp1, des1 = orb.detectAndCompute(img1,None) kp2, des2 = orb.detectAndCompute(img2,None) # create BFMatcher object bf = cv2.BFMatcher(cv2.NORM_HAMMING, crossCheck=True) # Match descriptors. matches = bf.match(des1,des2) # Sort them in the order of their distance. matches = sorted(matches, key = lambda x:x.distance) # Draw first 10 matches. img3 = cv2.drawMatches(img1,kp1,img2,kp2,matches[:10], flags=2) plt.imshow(img3) plt.show()如：3. Matcher 对象matches = bf.match(des1,des2) 得到的结果是 DMatch 对象，其具有如下几个属性：DMatch.distance - 描述子间的距离，值越小越好.DMatch.trainIdx - 训练描述子的索引DMatch.queryIdx - 查询描述子的索引DMatch.imgIdx - 训练图像的索引.4. SIFT 描述子和 Ratio Test 的 BF 匹配这里采用 BFMatcher.knnMatch() 得到 k 个最佳匹配. 示例中，设置 k=2，以便于使用 D.Lowe 在论文中的 ratio test.import numpy as np import cv2 from matplotlib import pyplot as plt img1 = cv2.imread('box.png',0) # queryImage img2 = cv2.imread('box_in_scene.png',0) # trainImage # Initiate SIFT detector sift = cv2.SIFT() # find the keypoints and descriptors with SIFT kp1, des1 = sift.detectAndCompute(img1,None) kp2, des2 = sift.detectAndCompute(img2,None) # BFMatcher with default params bf = cv2.BFMatcher() matches = bf.knnMatch(des1,des2, k=2) # Apply ratio test good = [] for m,n in matches: if m.distance < 0.75*n.distance: good.append([m]) # cv2.drawMatchesKnn expects list of lists as matches. img3 = cv2.drawMatchesKnn(img1,kp1,img2,kp2,good,flags=2) plt.imshow(img3) plt.show()如：5. FLANN 匹配FLANN，Fast Library for Approximate Nearest Neighbors. 其是针对大规模高维数据集进行快速最近邻搜索的优化算法库.FLANN Matcher 需要设定两个字典参数，以指定算法和对应的参数，分别为 IndexParams 和 SearchParams.[1] - IndexParams使用 SIFT、SURF等描述子，可以传递如下：index_params = dict(algorithm = FLANN_INDEX_KDTREE, trees = 5)使用 ORB 描述子是，可以传递如下：index_params= dict(algorithm = FLANN_INDEX_LSH, table_number = 6, # 12 key_size = 12, # 20 multi_probe_level = 1) #2[2] - SearchParams参数指定了索引里的树应该被递归遍历的次数，值越高准确度越高，但计算耗时也更久.search_params = dict(checks=100)6. SIFT 描述子和 Ratio Test 的 FLANN 匹配import numpy as np import cv2 from matplotlib import pyplot as plt img1 = cv2.imread('box.png',0) # queryImage img2 = cv2.imread('box_in_scene.png',0) # trainImage # Initiate SIFT detector sift = cv2.SIFT() # find the keypoints and descriptors with SIFT kp1, des1 = sift.detectAndCompute(img1,None) kp2, des2 = sift.detectAndCompute(img2,None) # FLANN parameters FLANN_INDEX_KDTREE = 0 index_params = dict(algorithm = FLANN_INDEX_KDTREE, trees = 5) search_params = dict(checks=50) # or pass empty dictionary flann = cv2.FlannBasedMatcher(index_params,search_params) matches = flann.knnMatch(des1,des2,k=2) # Need to draw only good matches, so create a mask matchesMask = [[0,0] for i in xrange(len(matches))] # ratio test as per Lowe's paper for i,(m,n) in enumerate(matches): if m.distance < 0.7*n.distance: matchesMask[i]=[1,0] draw_params = dict(matchColor = (0,255,0), singlePointColor = (255,0,0), matchesMask = matchesMask, flags = 0) img3 = cv2.drawMatchesKnn(img1,kp1,img2,kp2,matches,None,**draw_params) plt.imshow(img3) plt.show()如：https://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_feature2d/py_matcher/py_matcher.html
RootSift 是论文 Three things everyone should know to improve object retrieval - 2012 所提出的.A Comparative Analysis of RootSIFT and SIFT Methods for Drowsy Features Extraction - 2020当比较直方图时，使用欧氏距离通常比卡方距离或Hellinger核时的性能差，但是在使用 SIFT 特征点为什么一直都使用欧氏距离呢？不论是对 SIFT 特征点进行匹配，还是对 SIFT 特征集合进行聚类得到视觉词汇表，又或者对图像进行BoW编码，都使用的是欧氏距离. 但是 SIFT 特征描述子本质上也是一种直方图，为什么对 SIFT 特征描述子进行比较的时候要使用欧氏距离呢，有没有一种更精确的比较方法呢？SIFT 描述子统计的是关键点邻域的梯度直方图.论文作者认为之所以一直使用欧氏距离来测量 SIFT 特征的相似度，是由于在 SIFT 提出时，使用的是欧氏距离的度量，可以找出一种比较欧氏距离更为精确的度量方法. 故，提出了RootSift 对 SIFT 特征进行扩展.具体操作如下：在提取到 SIFT 描述向量 $x$ 后，进行如下处理，即可得到 RootSIFT:[1] - 对特征向量 $x$ 进行 $l_1$ 的归一化得到 $x'$;[2] - 对 $x'$ 的每一个元素求平方根；[3] - 进行 $l_2$ 归一化.(可选)[3]中，是否进行l2归一化，有些不一致. 在RootSIFT 论文 中并没有指出需要进行 l2 归一化，但是在 presentation, 却有 l2归一化. 也有认为，显式地执行L2规范化是不需要的. 通过采用L1规范，然后是平方根，已经有L2标准化的特征向量，不需要进一步的标准化.参考：图像检索(4):IF-IDF,RootSift,VLAD -- RootSIFT1. RootSIFT 实现Python 实现如：https://www.pyimagesearch.com/2015/04/13/implementing-rootsift-in-python-and-opencv/https://github.com/jrosebr1/imutils/blob/master/imutils/feature/rootsift.pyimport numpy as np import cv2 class RootSIFT: def __init__(self): # initialize the SIFT feature extractor #OpenCV2.4 # self.extractor = cv2.DescriptorExtractor_create("SIFT") #OpenCV3+ self.extractor = cv2.xfeatures2d.SIFT_create() def compute(self, image, kps, eps=1e-7): # compute SIFT descriptors kps, descs = self.extractor.compute(image, kps) # if there are no keypoints or descriptors, return an empty tuple if len(kps) == 0: return ([], None) # apply the Hellinger kernel by first L1-normalizing and taking the # square-root descs /= (descs.sum(axis=1, keepdims=True) + eps) descs = np.sqrt(descs) #descs /= (np.linalg.norm(descs, axis=1, ord=2) + eps) # return a tuple of the keypoints and descriptors return (kps, descs) # image = cv2.imread("test.jpg") gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY) # detect Difference of Gaussian keypoints in the image sift = cv2.xfeatures2d.SIFT_create() kps, descs = sift.detectAndCompute(img1, None) print "SIFT: kps=%d, descriptors=%s " % (len(kps), descs.shape) # extract RootSIFT descriptors root_sift = RootSIFT() (kps, descs) = root_sift.compute(image, kps) print "RootSIFT: kps=%d, descriptors=%s " % (len(kps), descs.shape)C++ 实现：for(int i = 0; i < siftFeature.rows; i ++){ // Conver to float type Mat f; siftFeature.row(i).convertTo(f,CV_32FC1); normalize(f,f,1,0,NORM_L1); // l1 normalize sqrt(f,f); // sqrt-root root-sift rootSiftFeature.push_back(f); }2. 基于 RootSIFT 的相似图搜索
OpenCV 库中 SIFT 特征提取实现，如：import cv2 sift = cv2.xfeatures2d.SIFT_create() img = cv2.imread('test.jpg') kps, des = sift.detectAndCompute(img, None)1. SIFT_create() 参数def SIFT_create(nfeatures=None, nOctaveLayers=None, contrastThreshold=None, edgeThreshold=None, sigma=None): """ SIFT_create([, nfeatures[, nOctaveLayers[, contrastThreshold[, edgeThreshold[, sigma]]]]]) -> retval Use cv.SIFT_create() instead """其中，[1] - nfeatures：特征点数目（算法对检测出的特征点排名，返回最好的nfeatures个特征点）[2] - nOctaveLayers：金字塔中每组的层数（算法中会自己计算这个值）[3] - contrastThreshold：过滤掉较差的特征点的对阈值. contrastThreshold越大，返回的特征点越少.[4] - edgeThreshold：过滤掉边缘效应的阈值. edgeThreshold越大，特征点越多（被过滤掉的越少）.[5] - sigma：金字塔第0层图像高斯滤波系数.注：[1] - 参数nfeatures指定最终返回的特征点数量，并不影响SIFT特征检测的结果.[2] - 参数nOctaveLayers和sigma主要影响图像高斯金字塔的构成. [3] - contrastThreshold和edgeThreshold 会影响在DOG中寻找极值点的过程与结果.在使用OpenCV SIFT特征检测算子时，如果提取效果不够理想，可根据实际情况调整参数后进行重新提取，如：import cv2 sift = cv2.xfeatures2d.SIFT_create( nfeatures=500, contrastThreshold=0.01, edgeThreshold=30, sigma=1.6) kps, des = sift.detectAndCompute(img_cv2, None)OpenCV 默认参数：nfeatures =0 nOctaveLayers =3 contrastThreshold = 0.04 edgeThreshold = 10 sigma =1.6参考：opencv3.1的SIFT特征检测参数图文详解 - 2016.11.042. kps 与 deskps, des = sift.detectAndCompute(img, None)其中，返回值 des为 128 为向量组成的列表； kps 为关键点列表，每个元素为一个 KeyPoint，其包含信息有：[1] - angle：角度，表示关键点的方向. 为了保证方向不变形，SIFT算法通过对关键点周围邻域进行梯度运算，求得该点方向. -1为初值.[2] - class_id：当要对图片进行分类时，可以用 class_id 对每个特征点进行区分. 未设定时为-1，需要自己设定.[3] - octave：代表是从金字塔哪一层提取的得到的数据. [4] - pt：关键点的坐标(x, y).[5] - response：响应程度，代表该点强壮大小，更确切的说，是该点角点的程度.[6] - size：该点直径的大小
摘自：策略算法工程师之路-基于内容的图像检索BoW - Bag of WordsSIFT - Scale Invariant Feature Transform,基于SIFT和CNN特征的通用检索流程 - AIUAIFrom: https://www.robots.ox.ac.uk/~vgg/publications/2012/Arandjelovic12/presentation.pdf1. BoW 模型BoW 模型最初被用在文本分类中，将文档表示成特征矢量. BoW 模型基本思想是，假定对于一个文本，忽略其词序和语法、句法，仅仅将其看做是一些词汇的集合，而文本中的每个词汇都是独立的. 如：文档1：John likes to watch movies. Mary likes too文档2：John also likes to watch football games.对于给定词典:Vocabulary ={ "John": 1, "likes":2, "to": 3, "watch":4, "movies":5, "also":6, "football":7, "games":8, "Mary": 9, "too": 10 } 上述文档可以表示成:文档1： [1, 2, 1, 1, 1, 0, 0, 0, 1, 1]文档2： [1, 1,1, 1, 0, 1, 1, 1, 0, 0]在计算机图像处理中也借鉴了BOW的思想，假设一张图像是有若干视觉词汇组成. 什么样的视觉词典是期望的?光照不敏感尺度不敏感视角变化不敏感旋转不变...2. SIFTSIFT特征之杂七杂八解析 - AIUAISIFT 算法简单来讲，可分为以下几步：2.1. 多分辨率图像金字塔想象人从不同距离看同一个事物，不同尺度的同一个事物从不同视角反映了事物的本质. 距离越远尺度越小更反映事物的全局特性，相反更反映事物的局部特性.2.2. 极值点检测为了寻找尺度空间的极值点，每个像素点要和其图像域(同一尺度空间)和尺度域(相邻的尺度空间)的所有相邻点进行比较，当其大于(或者小于)所有相邻点时，改点就是极值点. 如图所示，中间的检测点要和其所在图像的3×3邻域8个像素点，以及其相邻的上下两层的3×3领域18个像素点，共26个像素点进行比较. 2.3. 特征点梯度直方图生成通过对特征点周围的像素进行分块，计算块内梯度直方图，生成具有独特性的向量，这个向量是该区域图像信息的一种抽象，具有唯一性. 得到的梯度直方图就是最后SIFT特征向量. SIFT特征点的检测实例如图：3. SIFT 特征构建 BoW对提取的SIFT特征点K-means聚类:如下三个Cluster:如此，每个Cluster的中心即是需要的视觉词典(Visual Vocabulary)，如下:对图像检测出的SIFT特征点，指派给离它最近的中心点(将这一步称为hard-assignment)，如此得到图像的视觉词汇分布. 基于BOW的图像检索算法可以描述如下:[1] - 提取图像库中所有图像的局部特征(如SIFT、SURF、ORB等)，得到特征集合 $F$.[2] - 对特征集合 $F$ 进行聚类，得到 $k$ 个聚类中心 $\lbrace C_i | i = 1,...,k \rbrace$，每个聚类中心 $C_i$ 表示一个视觉词汇. 聚类中心的集合就是视觉词典.[3] - 计算特征 $f_i$ 属于哪个词汇 $C_i$(距离该聚类中心的距离最近.)[4] - 统计每个词汇 $C_i$ 在各图像中出现的频次，得到一维的特征向量，即 BoW 表示.基于BOW的方法仅仅关注图像的局部特征而忽略了图像的全局语义，导致会召回不少语义差异较大的图像.4. Python 实现import cv2 import numpy as np from sklearn.cluster import KMeans from scipy.cluster.vq import * from sklearn import preprocessing from PIL import Image #加载数据集 imgFiles = ['1.jpg', '2.jpg', '3.jpg', '4.jpg', '5.jpg'] #聚类中心数/视觉词典数 numWords = 64 #提取图像SIFT特征 sift = cv2.xfeatures2d.SIFT_create() des_list=[] for i in range(0, len(imgFiles)): img=cv2.imread(imgFiles[i]) #转换成灰度图像 gray=cv2.cvtColor(img,cv2.COLOR_RGB2GRAY) #sift kp, des = sift.detectAndCompute(gray,None) #kp, des = sift.detectAndCompute(img1, None) des_list.append((image_path, des)) #视觉词典生成 #生成向量数组(堆叠所有sift特征，以进行 K-means聚类) descriptors = des_list[0][1] for image_path, descriptor in des_list[1:]: descriptors = np.vstack((descriptors, descriptor)) #K-means聚类生成视觉词典 print("开始K-means聚类: %d words, %d key points" %(numWords, descriptors.shape[0])) voc, variance = kmeans(descriptors, numWords, 1) #图像视觉单词抽取 im_features = np.zeros((len(image_paths), numWords), "float32") for i in range(0,len(imgFiles)): #将图像的SIFT特征转换成视觉单词分布 words, distance = vq(des_list[i][1],voc) for w in words: #单张图像单词出现频次统计 im_features[i][w] += 1 #Tf-Idf计算 #统计idf nbr_occurences = np.sum((im_features > 0) * 1, axis = 0) idf = np.array(np.log((1.0*len(imgFiles)+1) / (1.0*nbr_occurences + 1)), 'float32') # tf*idf im_features = im_features*idf # L2 归一化 im_features = preprocessing.normalize(im_features, norm='l2') #查询 #读取图像 query = imgFiles[idx] im = cv2.imread(query) # 提取图像sift特征 des_list = [] gray = cv2.cvtColor(im, cv2.COLOR_RGB2GRAY) kp, des = sift_det.detectAndCompute(gray, None) des_list.append((image_path, des)) descriptors = des_list[0][1] # 抽取图像Tf-Idf特征 test_features = np.zeros((1, numWords), "float32") words, distance = vq(descriptors,voc) for w in words: test_features[0][w] += 1 test_features = test_features*idf # 归一化 test_features = preprocessing.normalize(test_features, norm='l2') # 计算所有图像的相似度 scores = np.dot(test_features, im_features.T) rank_ID = np.argsort(-scores) # 输出Top最相似的图像 for i, ID in enumerate(rank_ID[0][0:20]): img = Image.open(imgFiles[ID]) plt.imshow(img) plt.show()
SIFT- Search Site Local Invariant Feature[转]
原文：Local Invariant Feature Detectors学习、记录、备忘.
SIFT, Scale Invariant Feature Transform，尺度不变特征变换，由加拿大教授David G.Lowe提出的. SIFT特征对旋转、尺度缩放、亮度变化等保持不变性，是一种非常稳定的局部特征.注：该文是学习记录，比较杂七杂八，比较多的复制粘贴，这里提前把相关博文备注如下：[1] - SIFT特征详解-2015.10.03[2] - SIFT - 百度百科[3] - SIFT算法详解-2012.04.28[4] - SIFT特征提取分析 - 2012.06.06[5] - SIFT特征点提取 - 2018.01.15[6] - sift特征提取算法 - 2015.05.23[7] - 【动手学计算机视觉】第七讲：传统目标检测之SIFT特征 - 2019.10.30[8] - python利用opencv实现SIFT特征提取与匹配 - 2020.11.05[9] - 基于SIFT特征的图像检索 vs CNN - 2018.06.12[10] - Implementing sift in python a complete guide-2019[11] - 超强大的SIFT图像匹配技术详细指南（附Python代码）- 2019.11.20[12] - SIFT(尺度不变特征变换) - bilibili1. SIFT 特点图像的局部特征，对旋转、尺度缩放、亮度变化保持不变，对视角变化、仿射变换、噪声也保持一定程度的稳定性.独特性好，信息量丰富，适用于海量特征库进行快速、准确的匹配.多量性，即使是很少几个物体也可以产生大量的SIFT特征高速性，经优化的SIFT匹配算法甚至可以达到实时性扩招性，可以很方便的与其他的特征向量进行联合.SIFT算法可以解决的问题：目标的自身状态、场景所处的环境和成像器材的成像特性等因素影响图像配准、目标识别跟踪的性能. 而SIFT算法在一定程度上可解决：目标的旋转、缩放、平移（RST）图像仿射/投影变换（视点viewpoint）光照影响（illumination）目标遮挡（occlusion）杂物场景（clutter）噪声SIFT算法的实质是在不同的尺度空间上查找关键点(特征点)，并计算出关键点的方向. SIFT所查找到的关键点是一些十分突出，不会因光照，仿射变换和噪音等因素而变化的点，如角点、边缘点、暗区的亮点及亮区的暗点等.2. SIFT 特征检测步骤主要分四步：尺度空间的极值检测 搜索所有尺度空间上的图像，通过高斯微分函数来识别潜在的对尺度和选择不变的兴趣点.特征点定位 在每个候选的位置上，通过一个拟合精细模型来确定位置尺度，关键点的选取依据他们的稳定程度.特征方向赋值 基于图像局部的梯度方向，分配给每个关键点位置一个或多个方向，后续的所有操作都是对于关键点的方向、尺度和位置进行变换，从而提供这些特征的不变性.特征点描述 在每个特征点周围的邻域内，在选定的尺度上测量图像的局部梯度，这些梯度被变换成一种表示，这种表示允许比较大的局部形状的变形和光照变换.3. 尺度空间在一定的范围内，无论物体是大还是小，人眼都可以分辨出来. 然而计算机要有相同的能力却不是那么的容易，在未知的场景中，计算机视觉并不能提供物体的尺度大小，其中的一种方法是把物体不同尺度下的图像都提供给机器，让机器能够对物体在不同的尺度下有一个统一的认知. 在建立统一认知的过程中，要考虑的就是在图像在不同的尺度下都存在的特征点.3.1. 多分辨率图像金字塔在早期图像的多尺度通常使用图像金字塔表示形式. 图像金字塔是同一图像在不同的分辨率下得到的一组结果，其生成过程一般包括两个步骤：对原始图像进行平滑对处理后的图像进行降采样（通常是水平、垂直方向的1/2） 降采样后得到一系列不断尺寸缩小的图像. 显然，一个传统的金字塔中，每一层的图像是其上一层图像长、高的各一半. 多分辨率的图像金字塔虽然生成简单，但其本质是降采样，图像的局部特征则难以保持，也就是无法保持特征的尺度不变形.3.2. 高斯模糊使用高斯模糊技术(Gaussian Blur) 可以用来降低图像中的噪点. 对于图像中的每个像素，高斯模糊技术会基于其相邻像素计算一个值. 以下是应用高斯模糊之前和之后的图像示例. 如图所示，纹理和次要细节将从图像中删除，并且仅保留诸如形状和边缘之类的相关信息.高斯模糊成功地去除了图像中的噪点，强调了图像的重要特征.4. 高斯尺度空间可以通过图像的模糊程度来模拟人在距离物体由远到近时物体在视网膜上成像过程，距离物体越近其尺寸越大图像也越模糊，这就是高斯尺度空间，使用不同的参数模糊图像（分辨率不变），是尺度空间的另一种表现形式. 图像和高斯函数进行卷积运算能够对图像进行模糊，使用不同的“高斯核”可得到不同模糊程度的图像. 一副图像其高斯尺度空间可由其和不同的高斯卷积得到：$$ L(x, y, \sigma) = G(x, y, \sigma) * I(x, y) $$其中，$x$ 和 $y$ 是图像坐标，$G(x, y, \sigma)$ 是高斯函数:$$ G(x, y, \sigma) = \frac{1}{2 \pi \sigma^2} e ^{\frac{x^2 + y^2}{2 \sigma^2}} $$$\sigma$ 称为尺度空间因子，它是高斯正态分布的标准差，反映了图像被模糊的程度，其值越大图像越模糊，对应的尺度也就越大.$L(x, y, \sigma)$ 代表着图像的高斯尺度空间.构建尺度空间的目的是为了检测出在不同的尺度下都存在的特征点，而检测特征点较好的算子是 $\Delta ^2 G$(高斯拉普拉斯， LoG),$$ \Delta ^2 = \frac{\partial ^2}{\partial x^2} + \frac{\partial ^2}{\partial y^2} $$使用LoG虽然能较好的检测到图像中的特征点，但是其运算量过大，通常可使用DoG（差分高斯，Difference of Gaussina) 来近似计算 LoG.设 k 为相邻两个高斯尺度空间的比例因子，则 DoG 定义为：$$ D(x, y, \sigma) = [G(x, y, k \sigma) - G(x, y, \sigma)] * I(x, y) $$$$ D(x, y, \sigma) = L(x, y, k \sigma) - L(x, y, \sigma) $$从上式可以知道，将相邻的两个高斯空间的图像相减就得到了DoG的响应图像. 为了得到DoG图像，先要构建高斯尺度空间，而高斯的尺度空间可以在图像金字塔降采样的基础上加上高斯滤波得到，也就是对图像金字塔的每层图像使用不同的参数 $\sigma$ 进行高斯模糊，使每层金字塔有多张高斯模糊过的图像，这样一组图像就是一个 Octave. 降采样时，金字塔上边一组图像的第一张是由其下面一组图像倒数第三张降采样得到. 尺度空间是从单个图像生成的具有不同比例的图像的集合.易知，高斯金字塔有多组，每组又有多层. 一组中的多个层之间的尺度是不一样的(也就是使用的高斯参数 $\sigma$ 是不同的)，相邻的两层之间的尺度相差一个比例因子 k. 如果每组有 S 层，则$k = 2^{\frac{1}{S}}$. 上一组图像的最底层图像是由下一组中尺度为 $2 \sigma$ 的图像进行因子为 2 的降采样得到的(高斯金字塔)高斯金字塔的组数一般是：$$ o = [log_2 min(m, n)] - a $$其中，o 表示金字塔的层数，m、n 分别为图像的行数和列数. 减去的系数 $a$ 可以在 $[0, log_2(m, n)]$ 区间内任意取值，和具体需要的金字塔的顶层图像的大小有关.高斯模糊参数 $\sigma$ (尺度空间) 可由下面关系式得到：$$ \sigma(o, s) = \sigma_0 \cdot 2 ^{\frac{o+s}{S}} $$其中，o 为所在的组，s 为所在的层，$\sigma_0$ 为初始的尺度，S 为每组的层数.在 Lowe 的算法实现中，$\sigma_0 = 1.6, o_{min} = -1, S=3$，$o_{min}=-1$ 就是首先讲原图像的长和宽各扩展一倍.从上面可知，同一组内相邻层的图像尺度关系为：$$ \sigma _{s+1} = k \cdot \sigma_s = 2^{\frac{1}{S}} \cdot \sigma_s $$相邻组之间的尺度关系为：$$ \sigma_{o+1} = 2 \sigma_{o} $$示例：假设原始图像尺寸为 275x183，缩放图像尺寸为 138x93，对于该两个图像，将创建两个模糊图像，然而，要缩放多少次呢，而且每个缩放图像需要创建多少个后续的模糊图像呢？理想的数量为 4 个，并且对于每个 Octave，模糊图像的数目应该为 5 个.4.1. 高斯金字塔构建示例假设图像 I 为 512x512 图像，构建高斯金字塔步骤为：(从 0 开始计数，倒立的金字塔)[1] - 金字塔的组数，$log_2(512)=9$，减去因子 3，则构建的金字塔的组数为 6. 取每层的层数为 3.[2] - 构建第 0 组，将图像宽和高增加一倍，变成 1024x1024($I_0$)第 0 层 $I_0 * G(x, y, \sigma_0)$第 1 层 $I_1 * G(x, y, k \sigma_0)$第 2 层 $I_0 * G(x, y, k^2 \sigma_0)$[3] - 构建第 1 组，对 $I_0$ 降采样，变成 512x512 ($I_1$)第 0 层 $I_1 * G(x, y, 2 \sigma_0)$第 1 层 $I_1 * G(x, y, 2k \sigma_0)$第 2 层 $I_1 * G(x, y, 2k^2 \sigma_0)$[4] - 类似地，依次进行[5] - 构建第 o 组，第 s 层 $I_o * G(x, y, 2^o k^s \sigma_0)$高斯金字塔构建成功后，将每一组相邻的两层想减，即可得到 DoG 金字塔.相关源码：From: SIFT特征点提取 - 2018.01.15for (o = 0; o < octvs; o++)//金字塔组数为octvs， for (i = 0; i < intvls + 3; i++)//每一组有intvls + 3 层，intvls一般为3 { if (o == 0 && i == 0)//如果是第一组第1层 gauss_pyr[o][i] = cvCloneImage(base);//base 为原始灰度图像经过升采样或降采样得到的图像 /* base of new octvave is halved image from end of previous octave */ else if (i == 0)//建立非第一组的第1层 gauss_pyr[o][i] = downsample(gauss_pyr[o - 1][intvls]);//降采样图像 /* blur the current octave's last image to create the next one */ else//建立非第一组的非第1层 { gauss_pyr[o][i] = cvCreateImage(cvGetSize(gauss_pyr[o][i - 1]),IPL_DEPTH_32F, 1); cvSmooth(gauss_pyr[o][i - 1], gauss_pyr[o][i],CV_GAUSSIAN, 0, 0, sig[i], sig[i]);// sig[i]为模糊系数 }//cvSmooth 为平滑处理函数，也即模糊处理. CV_GAUSSIAN 为选用高斯函数对图像模糊 return gauss_pyr;//返回建好的金字塔4.2. 高斯差分金字塔构建差分金字塔的是在高斯金字塔的基础上操作的，其建立过程是：在高斯金子塔中的每组中相邻两层相减（下一层减上一层）就生成高斯差分金字塔.如图：示例： 以第一个 Octave 为例，图像都具有相同比例，通过在前一个图像上应用高斯模糊，来创建每个后续图像，如图：其源码实现如：for (o = 0; o < octvs; o++)//octvs为高斯金字塔组数 for (i = 0; i < intvls + 2; i++)//因为相减，故高斯金字塔中每组有（intvls + 2）层图像 { dog_pyr[o][i] = cvCreateImage(cvGetSize(gauss_pyr[o][i]),IPL_DEPTH_32F, 1); cvSub(gauss_pyr[o][i + 1], gauss_pyr[o][i], dog_pyr[o][i], NULL);//cvSub为opencv内置相减函数 } return dog_pyr;//返回高斯差分金字塔5. DoG 空间极值点检测为了寻找尺度空间的极值点，每个像素点要和其图像域（同一尺度空间）和尺度域（相邻的尺度空间）的所有相邻点进行比较，当其大于（或者小于）所有相邻点时，改点就是极值点. 极值点(关键点)是由DOG空间的局部极值点组成的，关键点的初步探查是通过同一组内各DoG相邻两层图像之间比较完成的. 为了寻找DoG函数的极值点，每一个像素点要和它所有的相邻点比较，看其是否比它的图像域和尺度域的相邻点大或者小. 如图下图所示，中间的检测点和它同尺度的8个相邻点和上下相邻尺度对应的 9×2 个点共26个点比较，以确保在尺度空间和二维图像空间都检测到极值点. 一个点如果在DOG尺度空间本层以及上下两层的26个邻域中是最大或最小值时，就认为该点是图像在该尺度下的一个特征点. 下图中将叉号点要比较的26个点都标为了绿色. 因为首尾两张图都只有一个相邻图，无法进行极值比较（因为极值定义中，极值点是连续的）. 为了满足尺度变化连续性，要在Octave每一组的顶层使用高斯模糊再生成3幅图像，则高斯金字塔有S+3幅图像，而DOG金字塔有S+2幅图像.5.1. 尺度变化的连续性假设 S=3，也就是每组有 3 层，则 $k = 2 ^{\frac{1}{S}} = 2 ^{\frac{1}{3}}$，也就是高斯金字塔每组有 (S-1)3 层图像，DoG金字塔每组有 (S-2)2 层图像. 在 DoG 金字塔的第一组有两层尺度分别为 $\sigma, k\sigma$，第二组有两层的尺度分别为 $2\sigma, 2k\sigma$，由于只有两项是无法比较取得极值的(只有左右两边均有值才能有极值).由于无法取得极值，那么就需要继续对每组的图像进行高斯模糊，使得尺度行形成 $\sigma, k\sigma, k^2\sigma, k^3 \sigma, k^4 \sigma$，这样就可以选择中间的三项 $k\sigma, k^2\sigma, k^3 \sigma$.对应的下一组由上一组降采样得到的三项是 $2k\sigma, 2k^2\sigma, 2k^3 \sigma$，其首项 $2k \sigma = 2 \cdot 2^{\frac{1}{3}}\sigma=2^{\frac{4}{3}}\sigma$，刚好与上一组的最后一项 $k^3 \sigma = 2^{\frac{3}{3}} \sigma$ 的尺度连续起来.5.2. DoG极值点检测实现参考实现如：if (val > 0)//极大值检测 { for (i = -1; i <= 1; i++) for (j = -1; j <= 1; j++) for (k = -1; k <= 1; k++) if (val < pixval32f(dog_pyr[octv][intvl + i], r + j, c + k))//pixval32f为提取图像像素位置上的灰度值 return 0;} else /* check for minimum */ { for (i = -1; i <= 1; i++) for (j = -1; j <= 1; j++) for (k = -1; k <= 1; k++) if (val > pixval32f(dog_pyr[octv][intvl + i], r + j, c + k))//r c为图像的行数和列数，dog_pyr为高斯差分图 return 0;}6. 删除不好的极值点(特征点)通过比较检测得到的DoG的局部极值点是在离散的空间搜索得到的，由于离散空间是对连续空间采样得到的结果，因此在离散空间找到的极值点不一定是真正意义上的极值点，如图，因此要设法将不满足条件的点剔除掉. 可以通过尺度空间DoG函数进行曲线拟合寻找极值点，这一步的本质是去掉DoG局部曲率非常不对称的点. 要剔除掉的不符合要求的点主要有两种：低对比度的特征点不稳定的边缘响应点6.1. 剔除低对比度的特征点候选特征点$x$，其偏移量定义为 $\Delta x$，其对比度为 $D(x)$的绝对值$|D(x)|$，对$D(x)$ 应用二阶泰勒展开式：$$ D(x) = D + \frac{\partial D^T}{\partial x} \Delta x + \frac{1}{2} \Delta x^T \frac{\partial ^2 D }{\partial x^2} \Delta x $$由于 $x$ 是 $D(x)$ 的极值点，所以对上式求导，并令其为 0，得到：$$ \Delta = - \frac{\partial ^2 D^{-1}}{\partial x^2} \frac{\partial D(x)}{\partial x} $$然后，再把求得的 $\Delta x$ 代入到 $D(x)$ 的泰勒展开式中，有：$$ D(\hat{x}) = D + \frac{1}{2} \frac{\partial D^T}{\partial x} \hat{x} $$设对比度的阈值为 T，若 $|D(\hat{x})| \ge T$，则该特征点保留，否则剔除.注：任意极值点在其 $(x_0, y_0, \sigma_0)$ 处的二阶泰勒展开式如：6.2. 剔除不稳定的边缘响应点在边缘梯度的方向上主曲率值比较大，而沿着边缘方向则主曲率值较小. 候选特征点的DoG函数D(x) 的主曲率与 2x2 Hessian 矩阵 H 的特征值成正比.$$ H= \left[ \begin{matrix} D_{xx} & D_{yx} \\ D_{xy} & D_{yy} \end{matrix} \right] $$其中, $D_{xx}, D_{xy}, D_{yx}$ 是候选点邻域对应位置的差分求得的. 为了避免求具体的值，可以使用 H 特征值的比例. 设 $\alpha = \lambda _{max}$ 为 H 的最大特征值，$\beta = \lambda _{min}$ 为 H 的最小特征值，则：$$ Tr(H) = D_{xx} + D_{yy} = \alpha + \beta $$$$ Det(H) = D_{xx} + D_{yy} - D_{xy}^2 = \alpha \cdot \beta $$其中，$Tr(H)$ 为矩阵 H 的迹，$Det(H)$ 为矩阵 H 的行列式.设 $\gamma = \frac{\alpha}{\beta}$ 表示最大特征值和最小特征值的比值，则：$$ \frac{Tr(H)^2}{Det(H)} = \frac{(\alpha + \beta)^2}{\alpha \beta} = \frac{(\gamma \beta + \beta)^2}{\gamma \beta^2} = \frac{(\gamma + 1)^2}{\gamma} $$上式的结果与两个特征值的比例有关，和具体的大小无关，当两个特征值想等时其值最小，并且随着 $\gamma$ 的增大而增大. 因此为了检测主曲率是否在某个阈值下，只需检测：$$ \frac{Tr(H)^2}{Det(H)} > \frac{(T_{\gamma} + 1)^2}{T_{\gamma}} $$如果上式成立，则剔除该特征点，否则保留. (Lowe 论文中取 $T_{\gamma} =10$).6.3. 泰勒插值的源码实现while (i < SIFT_MAX_INTERP_STEPS)//SIFT_MAX_INTERP_STEPS=5为最大迭代次数，避免长时迭代 { interp_step(dog_pyr, octv, intvl, r, c, &xi, &xr, &xc);// 泰勒展开拟合，xi,xr,xc依次为x、y、σ方向偏移量, if (ABS(xi) < 0.5 && ABS(xr) < 0.5 && ABS(xc) < 0.5)//如果当前偏移量绝对值中的每个值均小于0.5，退出迭代 break; c += cvRound(xc);//计算行坐标，cvRound 为四舍五入. r += cvRound(xr); intvl += cvRound(xi); if (intvl < 1 ||//不在计算的图像层中 intvl > intvls ||//高斯差分每组的层数为intvls c < SIFT_IMG_BORDER ||//靠近图像边缘5个像素的区域不做检测，SIFT_IMG_BORDER=5, r < SIFT_IMG_BORDER || c >= dog_pyr[octv][0]->width - SIFT_IMG_BORDER ||//靠近图像边缘5个像素的区域不做检测 r >= dog_pyr[octv][0]->height - SIFT_IMG_BORDER) { return NULL; } i++;//迭代计数 } # static void interp_step(IplImage*** dog_pyr, int octv, int intvl, int r, int c,double* xi, double* xr, double* xc) { CvMat* dD, *H, *H_inv, X; double x[3] = { 0 }; dD = deriv_3D(dog_pyr, octv, intvl, r, c);//一阶偏导数 H = hessian_3D(dog_pyr, octv, intvl, r, c);//Hessian 矩阵即二阶导数组成的矩阵 H_inv = cvCreateMat(3, 3, CV_64FC1); cvInvert(H, H_inv, CV_SVD);//求Hessian矩阵的逆矩阵 cvInitMatHeader(&X, 3, 1, CV_64FC1, x, CV_AUTOSTEP); cvGEMM(H_inv, dD, -1, NULL, 0, &X, 0); //cvGEMM为矩阵乘法，//第一个矩阵的系数；//H_inv、dD第一二个矩阵//-1矩阵前的常数//X结果矩阵 cvReleaseMat(&dD); cvReleaseMat(&H); cvReleaseMat(&H_inv); *xi = x[2]; *xr = x[1]; *xc = x[0]; }7. 求取特征点的主方向经过上面的步骤已经找到了在不同尺度下都存在的特征点，为了实现图像旋转不变性，需要给特征点的方向进行赋值. 利用特征点邻域像素的梯度分布特性来确定其方向参数，再利用图像的梯度直方图求取关键点局部结构的稳定方向. 找到了特征点，也就可以得到该特征点的尺度$\sigma$，也就可以得到特征点所在的尺度图像$$ L(x, y) = G(x, y, \sigma) * I(x, y) $$计算以特征点为中心，以 $3 \times 1.5\sigma$ 为半径的区域图像的幅角和幅值，每个点 $L(x, y)$ 的梯度的模 $m(x, y)$ 以及方向 $\theta (x, y)$ 可通过如下公式求得：$$ m(x, y) = \sqrt{[L(x+1, y) - L(x-1, y)]^2 + [L(x, y+1) - L(x, y-1)]^2} $$$$ \theta(x, y) = \text{arctan} \frac{L(x, y+1) - L(x, y-1)}{L(x+1, y) - L(x-1, y)} $$计算得到梯度方向后，就要使用直方图统计特征点邻域内像素对应的梯度方向和幅值. 梯度方向的直方图的横轴是梯度方向的角度（梯度方向的范围是0到360度，直方图每36度一个柱共10个柱，或者没45度一个柱共8个柱），纵轴是梯度方向对应梯度幅值的累加，在直方图的峰值就是特征点的主方向. 在完成关键点的梯度计算后，使用直方图统计领域内像素的梯度和方向. 梯度直方图将0~360度的方向范围分为36个柱(bins)，其中每柱10度. 如图所示，直方图的峰值方向代表了关键点的主方向，(为简化，图中只画了八个方向的直方图). 在Lowe的论文还提到了使用高斯函数对直方图进行平滑以增强特征点近的邻域点对关键点方向的作用，并减少突变的影响. 为了得到更精确的方向，通常还可以对离散的梯度直方图进行插值拟合. 对梯度方向直方图的平滑处理：Opencv 采用的平滑公式为：$$ H(i) = \frac{h(i-2) + h(i+2)}{16} + \frac{4 \times (h(i-1) + h(i+1))}{16} + \frac{6 \times h(i)}{16} $$其中，$i \in [0, 35]$，$h$ 和 $H$ 分别表示平滑前和平滑后的直方图. 由于角度是循环的，即 $0^o = 360^o$，如果出现 $h(i)$，$i$ 超出了 [0, 35] 的范围，则可以通过圆周循环的方法找到它所对应的、在 $0^o = 360^o$ 之间的值，如 $h(-1) = h(35)$.具体而言，关键点的方向可以由和主峰值最近的三个柱值通过抛物线插值得到. 梯度直方图抛物线插值:在梯度直方图中，当存在一个相当于主峰值80%能量的柱值时，则可以将这个方向认为是该特征点辅助方向. 所以，一个特征点可能检测到多个方向（也可以理解为，一个特征点可能产生多个坐标、尺度相同，但是方向不同的特征点）. Lowe在论文中指出: 15%的关键点具有多方向，而且这些点对匹配的稳定性很关键. 直方图峰值代表该关键点邻域内图像梯度的主方向，当存在另一个相当于主峰值 80% 能量的峰值时，则认为这个方向是该关键点的辅方向. 所以一个关键点可能检测得到多个方向，这可以增强匹配的鲁棒性.得到特征点的主方向后，对于每个特征点可以得到三个信息$(x, y, \sigma, \theta)$，即位置、尺度和方向. 由此可以确定一个SIFT特征区域，一个SIFT特征区域由三个值表示，中心表示特征点位置，半径表示关键点的尺度，箭头表示主方向. 具有多个方向的关键点可以被复制成多份，然后将方向值分别赋给复制后的特征点，一个特征点就产生了多个坐标、尺度相等，但是方向不同的特征点.8. 生成特征描述通过以上的步骤已经找到了SIFT特征点位置、尺度和方向信息，下面就需要使用一组向量来描述关键点也就是生成特征点描述子，这个描述符不只包含特征点，也含有特征点周围对其有贡献的像素点. 描述子应具有较高的独立性，以保证匹配率. 特征描述符的生成大致有三个步骤：校正旋转主方向，确保旋转不变性.生成描述子，最终形成一个128维的特征向量归一化处理，将特征向量长度进行归一化处理，进一步去除光照的影响.为了保证特征矢量的旋转不变性，要以特征点为中心，在附近邻域内将坐标轴旋转 $\theta$ (特征点的主方向)角度，即将坐标轴旋转为特征点的主方向. 旋转后邻域内像素的新坐标为：$$ \left[ \begin{matrix} x' \\ y' \end{matrix} \right] = \left[ \begin{matrix} cos \theta & -sin \theta \\ sin \theta & cos \theta \end{matrix} \right] \left[ \begin{matrix} x \\ y \end{matrix} \right] $$旋转后以主方向为中心取 8x8 的窗口. 下图所示，左图的中央为当前关键点的位置，每个小格代表为关键点邻域所在尺度空间的一个像素，求取每个像素的梯度幅值与梯度方向，箭头方向代表该像素的梯度方向，长度代表梯度幅值，然后利用高斯窗口对其进行加权运算. 最后在每个 4×4 的小块上绘制8个方向的梯度直方图，计算每个梯度方向的累加值，即可形成一个种子点，如图所示. 每个特征点由4个种子点组成，每个种子点有8个方向的向量信息. 这种邻域方向性信息联合增强了算法的抗噪声能力，同时对于含有定位误差的特征匹配也提供了比较理性的容错性. 与求主方向不同，此时每个种子区域的梯度直方图在0-360之间划分为8个方向区间，每个区间为45度，即每个种子点有8个方向的梯度强度信息. 在实际的计算过程中，为了增强匹配的稳健性，Lowe建议：对每个关键点使用4×4共16个种子点来描述，这样一个关键点就可以产生128维(4x4x8)的SIFT特征向量. 通过对特征点周围的像素进行分块，计算块内梯度直方图，生成具有独特性的向量，这个向量是该区域图像信息的一种抽象，具有唯一性.特征向量形成后，为了去除光照变化的影响，需要对它们进行归一化处理. 主要处理概括如图：9. SIFT 特征匹配在 SIFT 特征提取完成以后，即可生成图像的描述子.假设两幅图像 A、B 的 SIFT 描述子分别为 k1*128 维、k2*128 ，可以将两图的各个尺度的描述子进行匹配，如果两个特征点匹配上 128 维，则表示其匹配上了.当两幅图像的SIFT特征向量生成后，下一步我们采用关键点特征向量的欧式距离来作为两幅图像中关键点的相似性判定度量. 取图像 A 中的某个关键点，并找出其与图像 B 中欧式距离最近的前两个关键点，在这两个关键点中，如果最近的距离除以次近的距离少于某个比例阈值，则接受这一对匹配点. 降低这个比例阈值，SIFT匹配点数目会减少，但更加稳定. 为了排除因为图像遮挡和背景混乱而产生的无匹配关系的关键点, Lowe提出了比较最近邻距离与次近邻距离的方法, 距离比率ratio小于某个阈值的认为是正确匹配. 因为对于错误匹配,由于特征空间的高维性, 相似的距离可能有大量其他的错误匹配,从而它的ratio值比较高. Lowe推荐ratio的阈值为0.8. 但对大量任意存在尺度、旋转和亮度变化的两幅图片进行匹配，结果表明ratio取值在0. 4~0. 6之间最佳，小于0. 4的很少有匹配点，大于0. 6的则存在大量错误匹配点. 故，建议ratio的取值原则如下:ratio=0. 4 对于准确度要求高的匹配；ratio=0. 6 对于匹配点数目要求比较多的匹配； ratio=0. 5 一般情况下. 也可按如下原则:当最近邻距离<200时 ratio=0. 6，反之ratio=0. 4. ratio的取值策略能排分错误匹配点.10. OpenCV 实现#!/usr/bin/python3 #!--*-- coding: utf-8 --*-- import cv2 # imgfile1 = 'test1.jpg' imgfile2 = 'test2.jpg' img1 = cv2.imread(imgfile1) img2 = cv2.imread(imgfile2) # sift = cv2.xfeatures2d.SIFT_create() # 计算关键点和描述子 # kp为关键点keypoints, des为描述子descriptors kp1, des1 = sift.detectAndCompute(img1, None) kp2, des2 = sift.detectAndCompute(img2, None) # 关键点可视化 img3 = cv2.drawKeypoints(img1, kp1, img1, color=(0, 255, 255)) img4 = cv2.drawKeypoints(img2, kp2, img2, color=(0, 255, 255)) # M1：匹配方法一 # 描述子匹配 index_params = dict(algorithm=1, trees=6) search_params = dict(checks=50) flann = cv2.FlannBasedMatcher(index_params, search_params) # knn 计算描述子的匹配 matche = flann.knnMatch(des1, des2, k=2) matchesMask = [[0, 0] for i in range(len(matche))] # 匹配结果可视化 result = [] for m, n in matche: if m.distance < 0.6 * n.distance: result.append([m]) img5 = cv2.drawMatchesKnn(img1, kp1, img2, kp2, matche, None, flags=2) cv2.imshow("SIFTMatchResult", img5) cv2.waitKey(0) ## M2: 匹配方法二 ## K近邻算法求取在空间中距离最近的K个数据点，并将这些数据点归为一类 #matcher = cv2.BFMatcher() #raw_matches = matcher.knnMatch(des1, des2, k = 2) #good_matches = [] #for m1, m2 in raw_matches: # # 如果最接近和次接近的比值大于一个既定的值，那么保留这个最接近的值，认为它和其匹配的点为good_match # if m1.distance < ratio * m2.distance: # good_matches.append([m1]) #img5 = cv2.drawMatchesKnn(img1, kp1, img2, kp2, matche, None, flags=2)11. SIFT 的缺点SIFT在图像的不变特征提取方面拥有无与伦比的优势，但并不完美，仍然存在：实时性不高.有时特征点较少.对边缘光滑的目标无法准确提取特征点，等.如图所示，对模糊的图像和边缘平滑的图像，检测出的特征点过少，对圆更是无能为力. 近来不断有人改进，其中最著名的有SURF和CSIFT. 图像特征提取方法 - AIUAI.1. 尺度变化的连续性假设 S=3，也就是每组有 3 层，则 $k = 2 ^{\frac{1}{S}} = 2 ^{\frac{1}{3}}$，也就是高斯金字塔每组有 (S-1)3 层图像，DoG金字塔每组有 (S-2)2 层图像. 在 DoG 金字塔的第一组有两层尺度分别为 $\sigma, k\sigma$，第二组有两层的尺度分别为 $2\sigma, 2k\sigma$，由于只有两项是无法比较取得极值的(只有左右两边均有值才能有极值).由于无法取得极值，那么就需要继续对每组的图像进行高斯模糊，使得尺度行形成 $\sigma, k\sigma, k^2\sigma, k^3 \sigma, k^4 \sigma$，这样就可以选择中间的三项 $k\sigma, k^2\sigma, k^3 \sigma$.对应的下一组由上一组降采样得到的三项是 $2k\sigma, 2k^2\sigma, 2k^3 \sigma$，其首项 $2k \sigma = 2 \cdot 2^{\frac{1}{3}}\sigma=2^{\frac{4}{3}}\sigma$，刚好与上一组的最后一项 $k^3 \sigma = 2^{\frac{3}{3}} \sigma$ 的尺度连续起来.5.2. DoG极值点检测实现参考实现如：if (val > 0)//极大值检测 { for (i = -1; i <= 1; i++) for (j = -1; j <= 1; j++) for (k = -1; k <= 1; k++) if (val < pixval32f(dog_pyr[octv][intvl + i], r + j, c + k))//pixval32f为提取图像像素位置上的灰度值 return 0;} else /* check for minimum */ { for (i = -1; i <= 1; i++) for (j = -1; j <= 1; j++) for (k = -1; k <= 1; k++) if (val > pixval32f(dog_pyr[octv][intvl + i], r + j, c + k))//r c为图像的行数和列数，dog_pyr为高斯差分图 return 0;}6. 删除不好的极值点(特征点)通过比较检测得到的DoG的局部极值点是在离散的空间搜索得到的，由于离散空间是对连续空间采样得到的结果，因此在离散空间找到的极值点不一定是真正意义上的极值点，如图，因此要设法将不满足条件的点剔除掉. 可以通过尺度空间DoG函数进行曲线拟合寻找极值点，这一步的本质是去掉DoG局部曲率非常不对称的点. 要剔除掉的不符合要求的点主要有两种：低对比度的特征点不稳定的边缘响应点6.1. 剔除低对比度的特征点候选特征点$x$，其偏移量定义为 $\Delta x$，其对比度为 $D(x)$的绝对值$|D(x)|$，对$D(x)$ 应用二阶泰勒展开式：$$ D(x) = D + \frac{\partial D^T}{\partial x} \Delta x + \frac{1}{2} \Delta x^T \frac{\partial ^2 D }{\partial x^2} \Delta x $$由于 $x$ 是 $D(x)$ 的极值点，所以对上式求导，并令其为 0，得到：$$ \Delta = - \frac{\partial ^2 D^{-1}}{\partial x^2} \frac{\partial D(x)}{\partial x} $$然后，再把求得的 $\Delta x$ 代入到 $D(x)$ 的泰勒展开式中，有：$$ D(\hat{x}) = D + \frac{1}{2} \frac{\partial D^T}{\partial x} \hat{x} $$设对比度的阈值为 T，若 $|D(\hat{x})| \ge T$，则该特征点保留，否则剔除.注：任意极值点在其 $(x_0, y_0, \sigma_0)$ 处的二阶泰勒展开式如：6.2. 剔除不稳定的边缘响应点在边缘梯度的方向上主曲率值比较大，而沿着边缘方向则主曲率值较小. 候选特征点的DoG函数D(x) 的主曲率与 2x2 Hessian 矩阵 H 的特征值成正比.$$ H= \left[ \begin{matrix} D_{xx} & D_{yx} \\ D_{xy} & D_{yy} \end{matrix} \right] $$其中, $D_{xx}, D_{xy}, D_{yx}$ 是候选点邻域对应位置的差分求得的. 为了避免求具体的值，可以使用 H 特征值的比例. 设 $\alpha = \lambda _{max}$ 为 H 的最大特征值，$\beta = \lambda _{min}$ 为 H 的最小特征值，则：$$ Tr(H) = D_{xx} + D_{yy} = \alpha + \beta $$$$ Det(H) = D_{xx} + D_{yy} - D_{xy}^2 = \alpha \cdot \beta $$其中，$Tr(H)$ 为矩阵 H 的迹，$Det(H)$ 为矩阵 H 的行列式.设 $\gamma = \frac{\alpha}{\beta}$ 表示最大特征值和最小特征值的比值，则：$$ \frac{Tr(H)^2}{Det(H)} = \frac{(\alpha + \beta)^2}{\alpha \beta} = \frac{(\gamma \beta + \beta)^2}{\gamma \beta^2} = \frac{(\gamma + 1)^2}{\gamma} $$上式的结果与两个特征值的比例有关，和具体的大小无关，当两个特征值想等时其值最小，并且随着 $\gamma$ 的增大而增大. 因此为了检测主曲率是否在某个阈值下，只需检测：$$ \frac{Tr(H)^2}{Det(H)} > \frac{(T_{\gamma} + 1)^2}{T_{\gamma}} $$如果上式成立，则剔除该特征点，否则保留. (Lowe 论文中取 $T_{\gamma} =10$).6.3. 泰勒插值的源码实现while (i < SIFT_MAX_INTERP_STEPS)//SIFT_MAX_INTERP_STEPS=5为最大迭代次数，避免长时迭代 { interp_step(dog_pyr, octv, intvl, r, c, &xi, &xr, &xc);// 泰勒展开拟合，xi,xr,xc依次为x、y、σ方向偏移量, if (ABS(xi) < 0.5 && ABS(xr) < 0.5 && ABS(xc) < 0.5)//如果当前偏移量绝对值中的每个值均小于0.5，退出迭代 break; c += cvRound(xc);//计算行坐标，cvRound 为四舍五入. r += cvRound(xr); intvl += cvRound(xi); if (intvl < 1 ||//不在计算的图像层中 intvl > intvls ||//高斯差分每组的层数为intvls c < SIFT_IMG_BORDER ||//靠近图像边缘5个像素的区域不做检测，SIFT_IMG_BORDER=5, r < SIFT_IMG_BORDER || c >= dog_pyr[octv][0]->width - SIFT_IMG_BORDER ||//靠近图像边缘5个像素的区域不做检测 r >= dog_pyr[octv][0]->height - SIFT_IMG_BORDER) { return NULL; } i++;//迭代计数 } # static void interp_step(IplImage*** dog_pyr, int octv, int intvl, int r, int c,double* xi, double* xr, double* xc) { CvMat* dD, *H, *H_inv, X; double x[3] = { 0 }; dD = deriv_3D(dog_pyr, octv, intvl, r, c);//一阶偏导数 H = hessian_3D(dog_pyr, octv, intvl, r, c);//Hessian 矩阵即二阶导数组成的矩阵 H_inv = cvCreateMat(3, 3, CV_64FC1); cvInvert(H, H_inv, CV_SVD);//求Hessian矩阵的逆矩阵 cvInitMatHeader(&X, 3, 1, CV_64FC1, x, CV_AUTOSTEP); cvGEMM(H_inv, dD, -1, NULL, 0, &X, 0); //cvGEMM为矩阵乘法，//第一个矩阵的系数；//H_inv、dD第一二个矩阵//-1矩阵前的常数//X结果矩阵 cvReleaseMat(&dD); cvReleaseMat(&H); cvReleaseMat(&H_inv); *xi = x[2]; *xr = x[1]; *xc = x[0]; }7. 求取特征点的主方向经过上面的步骤已经找到了在不同尺度下都存在的特征点，为了实现图像旋转不变性，需要给特征点的方向进行赋值. 利用特征点邻域像素的梯度分布特性来确定其方向参数，再利用图像的梯度直方图求取关键点局部结构的稳定方向. 找到了特征点，也就可以得到该特征点的尺度$\sigma$，也就可以得到特征点所在的尺度图像$$ L(x, y) = G(x, y, \sigma) * I(x, y) $$计算以特征点为中心，以 $3 \times 1.5\sigma$ 为半径的区域图像的幅角和幅值，每个点 $L(x, y)$ 的梯度的模 $m(x, y)$ 以及方向 $\theta (x, y)$ 可通过如下公式求得：$$ m(x, y) = \sqrt{[L(x+1, y) - L(x-1, y)]^2 + [L(x, y+1) - L(x, y-1)]^2} $$$$ \theta(x, y) = \text{arctan} \frac{L(x, y+1) - L(x, y-1)}{L(x+1, y) - L(x-1, y)} $$计算得到梯度方向后，就要使用直方图统计特征点邻域内像素对应的梯度方向和幅值. 梯度方向的直方图的横轴是梯度方向的角度（梯度方向的范围是0到360度，直方图每36度一个柱共10个柱，或者没45度一个柱共8个柱），纵轴是梯度方向对应梯度幅值的累加，在直方图的峰值就是特征点的主方向. 在完成关键点的梯度计算后，使用直方图统计领域内像素的梯度和方向. 梯度直方图将0~360度的方向范围分为36个柱(bins)，其中每柱10度. 如图所示，直方图的峰值方向代表了关键点的主方向，(为简化，图中只画了八个方向的直方图). 在Lowe的论文还提到了使用高斯函数对直方图进行平滑以增强特征点近的邻域点对关键点方向的作用，并减少突变的影响. 为了得到更精确的方向，通常还可以对离散的梯度直方图进行插值拟合. 对梯度方向直方图的平滑处理：Opencv 采用的平滑公式为：$$ H(i) = \frac{h(i-2) + h(i+2)}{16} + \frac{4 \times (h(i-1) + h(i+1))}{16} + \frac{6 \times h(i)}{16} $$其中，$i \in [0, 35]$，$h$ 和 $H$ 分别表示平滑前和平滑后的直方图. 由于角度是循环的，即 $0^o = 360^o$，如果出现 $h(i)$，$i$ 超出了 [0, 35] 的范围，则可以通过圆周循环的方法找到它所对应的、在 $0^o = 360^o$ 之间的值，如 $h(-1) = h(35)$.具体而言，关键点的方向可以由和主峰值最近的三个柱值通过抛物线插值得到. 梯度直方图抛物线插值:在梯度直方图中，当存在一个相当于主峰值80%能量的柱值时，则可以将这个方向认为是该特征点辅助方向. 所以，一个特征点可能检测到多个方向（也可以理解为，一个特征点可能产生多个坐标、尺度相同，但是方向不同的特征点）. Lowe在论文中指出: 15%的关键点具有多方向，而且这些点对匹配的稳定性很关键. 直方图峰值代表该关键点邻域内图像梯度的主方向，当存在另一个相当于主峰值 80% 能量的峰值时，则认为这个方向是该关键点的辅方向. 所以一个关键点可能检测得到多个方向，这可以增强匹配的鲁棒性.得到特征点的主方向后，对于每个特征点可以得到三个信息$(x, y, \sigma, \theta)$，即位置、尺度和方向. 由此可以确定一个SIFT特征区域，一个SIFT特征区域由三个值表示，中心表示特征点位置，半径表示关键点的尺度，箭头表示主方向. 具有多个方向的关键点可以被复制成多份，然后将方向值分别赋给复制后的特征点，一个特征点就产生了多个坐标、尺度相等，但是方向不同的特征点.8. 生成特征描述通过以上的步骤已经找到了SIFT特征点位置、尺度和方向信息，下面就需要使用一组向量来描述关键点也就是生成特征点描述子，这个描述符不只包含特征点，也含有特征点周围对其有贡献的像素点. 描述子应具有较高的独立性，以保证匹配率. 特征描述符的生成大致有三个步骤：校正旋转主方向，确保旋转不变性.生成描述子，最终形成一个128维的特征向量归一化处理，将特征向量长度进行归一化处理，进一步去除光照的影响.为了保证特征矢量的旋转不变性，要以特征点为中心，在附近邻域内将坐标轴旋转 $\theta$ (特征点的主方向)角度，即将坐标轴旋转为特征点的主方向. 旋转后邻域内像素的新坐标为：$$ \left[ \begin{matrix} x' \\ y' \end{matrix} \right] = \left[ \begin{matrix} cos \theta & -sin \theta \\ sin \theta & cos \theta \end{matrix} \right] \left[ \begin{matrix} x \\ y \end{matrix} \right] $$旋转后以主方向为中心取 8x8 的窗口. 下图所示，左图的中央为当前关键点的位置，每个小格代表为关键点邻域所在尺度空间的一个像素，求取每个像素的梯度幅值与梯度方向，箭头方向代表该像素的梯度方向，长度代表梯度幅值，然后利用高斯窗口对其进行加权运算. 最后在每个 4×4 的小块上绘制8个方向的梯度直方图，计算每个梯度方向的累加值，即可形成一个种子点，如图所示. 每个特征点由4个种子点组成，每个种子点有8个方向的向量信息. 这种邻域方向性信息联合增强了算法的抗噪声能力，同时对于含有定位误差的特征匹配也提供了比较理性的容错性. 与求主方向不同，此时每个种子区域的梯度直方图在0-360之间划分为8个方向区间，每个区间为45度，即每个种子点有8个方向的梯度强度信息. 在实际的计算过程中，为了增强匹配的稳健性，Lowe建议：对每个关键点使用4×4共16个种子点来描述，这样一个关键点就可以产生128维(4x4x8)的SIFT特征向量. 通过对特征点周围的像素进行分块，计算块内梯度直方图，生成具有独特性的向量，这个向量是该区域图像信息的一种抽象，具有唯一性.特征向量形成后，为了去除光照变化的影响，需要对它们进行归一化处理. 主要处理概括如图：9. SIFT 特征匹配在 SIFT 特征提取完成以后，即可生成图像的描述子.假设两幅图像 A、B 的 SIFT 描述子分别为 k1*128 维、k2*128 ，可以将两图的各个尺度的描述子进行匹配，如果两个特征点匹配上 128 维，则表示其匹配上了.当两幅图像的SIFT特征向量生成后，下一步我们采用关键点特征向量的欧式距离来作为两幅图像中关键点的相似性判定度量. 取图像 A 中的某个关键点，并找出其与图像 B 中欧式距离最近的前两个关键点，在这两个关键点中，如果最近的距离除以次近的距离少于某个比例阈值，则接受这一对匹配点. 降低这个比例阈值，SIFT匹配点数目会减少，但更加稳定. 为了排除因为图像遮挡和背景混乱而产生的无匹配关系的关键点, Lowe提出了比较最近邻距离与次近邻距离的方法, 距离比率ratio小于某个阈值的认为是正确匹配. 因为对于错误匹配,由于特征空间的高维性, 相似的距离可能有大量其他的错误匹配,从而它的ratio值比较高. Lowe推荐ratio的阈值为0.8. 但对大量任意存在尺度、旋转和亮度变化的两幅图片进行匹配，结果表明ratio取值在0. 4~0. 6之间最佳，小于0. 4的很少有匹配点，大于0. 6的则存在大量错误匹配点. 故，建议ratio的取值原则如下:ratio=0. 4 对于准确度要求高的匹配；ratio=0. 6 对于匹配点数目要求比较多的匹配； ratio=0. 5 一般情况下. 也可按如下原则:当最近邻距离<200时 ratio=0. 6，反之ratio=0. 4. ratio的取值策略能排分错误匹配点.10. OpenCV 实现#!/usr/bin/python3 #!--*-- coding: utf-8 --*-- import cv2 # imgfile1 = 'test1.jpg' imgfile2 = 'test2.jpg' img1 = cv2.imread(imgfile1) img2 = cv2.imread(imgfile2) # sift = cv2.xfeatures2d.SIFT_create() # 计算关键点和描述子 # kp为关键点keypoints, des为描述子descriptors kp1, des1 = sift.detectAndCompute(img1, None) kp2, des2 = sift.detectAndCompute(img2, None) # 关键点可视化 img3 = cv2.drawKeypoints(img1, kp1, img1, color=(0, 255, 255)) img4 = cv2.drawKeypoints(img2, kp2, img2, color=(0, 255, 255)) # M1：匹配方法一 # 描述子匹配 index_params = dict(algorithm=1, trees=6) search_params = dict(checks=50) flann = cv2.FlannBasedMatcher(index_params, search_params) # knn 计算描述子的匹配 matche = flann.knnMatch(des1, des2, k=2) matchesMask = [[0, 0] for i in range(len(matche))] # 匹配结果可视化 result = [] for m, n in matche: if m.distance < 0.6 * n.distance: result.append([m]) img5 = cv2.drawMatchesKnn(img1, kp1, img2, kp2, matche, None, flags=2) cv2.imshow("SIFTMatchResult", img5) cv2.waitKey(0) ## M2: 匹配方法二 ## K近邻算法求取在空间中距离最近的K个数据点，并将这些数据点归为一类 #matcher = cv2.BFMatcher() #raw_matches = matcher.knnMatch(des1, des2, k = 2) #good_matches = [] #for m1, m2 in raw_matches: # # 如果最接近和次接近的比值大于一个既定的值，那么保留这个最接近的值，认为它和其匹配的点为good_match # if m1.distance < ratio * m2.distance: # good_matches.append([m1]) #img5 = cv2.drawMatchesKnn(img1, kp1, img2, kp2, matche, None, flags=2)11. SIFT 的缺点SIFT在图像的不变特征提取方面拥有无与伦比的优势，但并不完美，仍然存在：实时性不高.有时特征点较少.对边缘光滑的目标无法准确提取特征点，等.如图所示，对模糊的图像和边缘平滑的图像，检测出的特征点过少，对圆更是无能为力. 近来不断有人改进，其中最著名的有SURF和CSIFT. 图像特征提取方法 - AIUAI