图像特征点检测与匹配评价准则——量化

欢迎转载，转载请注明出处，谢谢！

目前图像匹配中，局部特征匹配占据了绝大部分，常用的局部特征匹配方法有Harris、SIFT、SURF、ORB等等，不同的特征点检测和匹配方法尤其独特的优势和不足； 
特征点匹配经过Ransac算法优化后仍存在错误匹配点对，需要优化后的匹配结果进行量化评价； 
特征点检测和匹配评价一般包括两个部分，分别为检测和匹配的评价。

1、特征点检测评价

评价特征点检测方法的优劣常常用到repeatability（重复率）这个概念。 
《a performance evaluation of local descriptors》一文提出了特征点检测器和匹配结果的评价准则，我的理解是：图A、B是两幅待匹配图像，图A映射到图B有一个单应性矩阵H1，图B映射到图A有单应性矩阵H2，图A检测出N1个特征点，图B检测出N2个特征点，因为图像A和B有部分图像不重叠，故将A图检测的特征点坐标由H1算出在B图的坐标，去掉不合格（计算结果超出在B图像坐标）的特征点，剩下的特征点数记为n1；同样，B图的特征点经过处理剩下n2个；分母便是min(n1, n2)。将图A剩下的特征点由H1计算出在图B中的坐标，与图B检测出的特征点的坐标求距离，即dist(h1*a1, b1)，若距离小于阈值ε，则认为是重复的，这么做是因为得到的单应性矩阵不一定完全精确以及一些别的误差原因。

《An affine invariant interest point detector》中对repeatability定义进行了延伸，针对affine invariant的图像特征点。 

repeatability计算式的分母可以和原来定义相同，分子的定义correspondence满足两个条件： 
1、同之前定义，即设定ε=1.5； 
2、regions映射到另一幅图像中重叠误差小于0.2。 
这里第二点是因为特征点检测器得到的特征点实际上是根据特征点附近区域计算梯度来选取特征点，实际上检测的特征区域（region），如Harris-Affine检测的是一个椭圆区域，SIFT检测的是一个圆形区域；当左图一特征点区域在经过乘上单应性矩阵H后实际，该区域会发生变化，所以需要判断变化之后的区域与右图某一特征点区域的重叠误差，误差<0.2，则判断自是一对特征点。 

2、特征点匹配评价

特征点匹配的评价一般涉及到两个概念，即Recall（召回率）和Precision（精确率），我们先来了解一下这两个概念。 
精确率和召回率是广泛用于信息检索和统计学分类领域的两个度量值，用来评价结果的质量。 
——精确率是检索出相关文档数与检索出的文档总数的比率，衡量的是检索系统的查准率； 
——召回率是指检索出的相关文档数和文档库中所有的相关文档数的比率，衡量的是检索系统的查全率。 
在信息检索中：

1. 精确率 = 提取出的正确信息条数 / 提取出的信息条数
2. 召回率 = 提取出的正确信息条数 / 样本中的信息条数

下面是一个非常经典的解释Recall和Precision的图，我做了一点儿说明： 

Precision和Recall指标有时候会出现的矛盾的情况，可以绘制Precision-Recall曲线，曲线越靠上部，结果越好。

对比信息检索的评价，在特征点匹配时： 
A指的是实际为匹配点对，且该匹配算法得到了这些匹配对； 
B指的是实际为匹配点对，但该匹配算法没有得到这些匹配对； 
C指的是实际为错误匹配对，但该算法得到了这些匹配对； 
即Precision为匹配结果中有多少是准确的，Recall就是所有正确的匹配结果有多少通过匹配算法得到了。

在第一部分中已经说明了根据两图像间的单应性矩阵H计算两个特征点是否重复，这里设置一个correspondence变量，作为重复特征点的数量，则匹配结果评价的两个参数Precision和Recall有如下公式： 

 

公式中，#correct matches表示匹配结果中正确匹配点对，#correspondences表示特征点检测时重复特征点对（特征点重复即认为实际为正确匹配点对，但可能被匹配算法匹配上，也可能未匹配上），#false matches表示匹配结果中错误匹配点对。其中，正确与错误匹配点对可以用上面提到的两个判断特征点是否重复的条件来判定。

因此，在评价特征点匹配结果时，首先应当知道两个图像间的单应性矩阵H，然后通过特征点检测算法得到左右图像中的特征点信息，根据单应性矩阵得到重复特征点数，即#correspondences；接着由特征点匹配算法得到的匹配结果计算Precision和Recall，绘制1-precision和recall曲线，曲线靠上方的结果较出色！ 
参考曲线图： 

 
（图片来源——《A Performance Evaluation of Local Descriptors》）

以上内容为我的理解，可能存在理解不当的地方，欢迎指正！

转载请注明出处：图像特征点检测与匹配评价准则——量化

参考资料： 
介绍Precision和Recall的博客： 
http://blog.csdn.net/pirage/article/details/9851339 
http://bookshadow.com/weblog/2014/06/10/precision-recall-f-measure/ 
文章《A Performance Evaluation of Local Descriptors》 
https://www.robots.ox.ac.uk/~vgg/research/affine/det_eval_files/mikolajczyk_pami2004.pdf 
K Mikolajczyk大神的特征点重复率计算相关代码： 
http://www.robots.ox.ac.uk/~vgg/research/affine/descriptors.html#binaries 
特征点检测匹配标准数据库： 

http://www.robots.ox.ac.uk/~vgg/data/data-aff.html

研究特征提取的方法，找到提取出来的特征精度最高，这里的最高是指在多帧运动变化图片里应用此算法都能找到固定位置的特征，而且这个特征不限，可以是边缘、角点或区域块等，然后根据不同特征、其运动变化及图像成像质量提供不同的算法，以达到最优。这个方向刚开始的确十分头大，完全不知道从什么着手，是边缘提取算法呢？还是先作目标分割，大概想一两天，决定先从特征点（包括角点）提取算法做起，先找到在特征点算法里非常优秀的，然后结合各种背景（噪声，运动变化）测试其变化。

这里主要的思路来自于一个仿射不变特征提取方法的检测：http://www.robots.ox.ac.uk/~vgg/research/affine/

这个网站上面主要提供了一个仿射不变特征提取方法的评估，它提供了八组图片，分别从噪声、模糊、视点变化、缩放、旋转、亮度变化来评估其仿射不变特征提取方法的效果。不过它主要是对相似区域进行提取，涉及到椭圆区域的仿射变换，而我主要进行点特征点提取的评估，要稍微简单点，当然也借鉴了其代码（这里面只有评估代码，没有特征提取算法的代码，不过提供了一系列算法的执行文件）

我首先要测试的sift算法对视点变化的评估，从上面的网站里下载了几个图片【img1.ppm、img2.ppm】、变化矩阵【H1to2p】及特征点提取文件【img1.haraff.sift】，我主要从特征向量匹配程度和实际匹配的位置进行比较。我的代码如下：

1. function [matchFeaIndex, matchFeaRatio,matchLocal1]=matchPointFea(file1,file2,H,image1,image2)
2. % 主要用特征点匹配效果测试：
3. % 用法：[matchFeaIndex, matchFeaRatio,matchLocal1]=matchPointFea('img1.haraff.sift','img2.haraff.sift','H1to2p','img1.ppm','img2.ppm');
4. % 输入:
5. % file1 - 图像1特征点检测文件
6. % file2 - 图像2特征点检测文件
7. % H - 图像1到图像2的仿射变换矩阵
8. % image1 - 图像1
9. % image2 - 图像2
10. % 输出:
11. % matchFeaIndex - 保存图像1中同图像2匹配的特征点在file2中的标号，若不匹配则为0。
12. % matchFeaRatio - 对应于matchIndex，保存特征点匹配程度，不匹配则为0。
13. % matchLocal1 - 保存图像2中特征点位置同经过仿射变换1之间的位置差别
14. % matchLocal2 - 保存图像1中特征点位置同经过仿射变换2之间的位置差别（保留）
15. % 特征点检测文件格式：
16. % file1, file2:
17. % x1 y1 a1 b1 c1 c2 c3 ...
18. % x2 y2 a2 b2 c1 c2 c3 ...
19. %---------------------
20. % 第一行 保存描述子数目及大小
21. % x, y - 特征点位置
22. % a, b - 特征点大小及方向（主要针对于sift变换）
23. % d1 d2 d3 ... - 特征描述子向量（如果小于，则无效）
25. close all;
26. % loc存放位置及大小角度，des存放特征描述子向量,dimdesc存放描述子大小及数目
27. [des1, loc1, dimdesc1]=loadFeatures(file1);
28. [des2, loc2, dimdesc2]=loadFeatures(file2);
29. % 导入图像
30. Im1=imread(image1);
31. Im2=imread(image2);
32. % 导入仿射矩阵文件
33. H=load(H);
35. fprintf(1,'numbers and length of descriptor in file1 %d %d\n',dimdesc1(1),dimdesc1(2));
36. fprintf(1,'numbers and length of descriptor in file2 %d %d\n',dimdesc2(1),dimdesc2(2));
38. if dimdesc1(2)&gt;1 &amp;&amp; dimdesc1(2)==dimdesc2(2)
39. fprintf(1,'%s, %s look like files with descriptors...\n',file1,file2);
40. else
41. error('Different descriptor dimension in %s or %s files.',file1,file2);
42. end
44. % 计算两个特征向量的匹配程度可以通过向量空间余弦相似度来衡量.
45. % 设置比值distRatio，保证所匹配的特征点具有显著相似度，即与第二相似度有较大差异.
46. distRatio = 0.6;
48. % 在图像2中找到图像1中每一个匹配点.
49. des2t = des2'; % 转置
50. desNum = dimdesc1(1); % 图像1特征点个数
51. matchFeaIndex=zeros(desNum,1);
52. matchFeaRatio=zeros(desNum,1);
53. matchNum=0;
54. for i = 1 : desNum
55. dotprods = des1(i,:) * des2t; % 计算乘积项
56. [vals,indx] = sort(acos(dotprods)); % 排序余弦相似度
57. % 找到最大余弦相似度.
58. if (vals(1) &lt; distRatio * vals(2))
59. matchFeaIndex(i) = indx(1);
60. matchFeaRatio(i) = vals(1);
61. matchNum=matchNum+1;
62. else
63. matchFeaIndex(i) = 0;
64. end
65. end
67. fprintf(1,'numbers of match descriptor is %d \n',matchNum);
68. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
69. % 确定图像1中特征点变换到图像2中的位置，及图像2变换到图像1中的位置
70. % 要进行仿射变换测试，要求图像1和2必须大小相同
71. if size(Im1)==size(Im2)
72. HI=H(:, 1:3); %图1到图2
73. H=inv(HI); %图2到图1
74. fprintf(1,'Projecting 1 to 2...\n');
75. loc1t=projectPoints(loc1',HI);
76. loc1t=loc1t';
77. fprintf(1,'and 2 to 1...\n');
78. loc2t=projectPoints(loc2',H);
79. loc2t=loc2t';
80. ImageSize=size(Im1);
81. matchLocal1=matchLocal(loc1t,loc2,ImageSize,matchFeaIndex);
82. end
84. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
85. % 画图
86. % showkeys(image1,loc1); % 特征点显示
87. % showkeys(image2,loc2);
89. % 显示匹配图像过程
90. Im=appendimages(Im1,Im2);
91. % Show a figure with lines joining the accepted matches.
92. figure('Position', [0 0 size(Im,2) size(Im,1)]);
93. colormap('gray');
94. imagesc(Im);
95. hold on;
96. cols1 = size(Im1,2);
97. for i = 1: size(des1,1)
98. if (matchFeaIndex(i) &gt; 0 &amp;&amp; matchLocal1(i)&lt;5)
99. line([loc1(i,1) loc2(matchFeaIndex(i),1)+cols1], ...
100. [loc1(i,2) loc2(matchFeaIndex(i),2)], 'Color', 'c'); %画对应线条
101. plot(loc1(i,1),loc1(i,2),'r.'); %画特征点
102. plot(loc2(matchFeaIndex(i),1)+cols1,loc2(matchFeaIndex(i),2),'g.'); %画特征点
103. end
104. end
105. hold off;
107. end
109. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
110. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
111. % 点的投影运算
112. function loct=projectPoints(loc,H)
113. num=size(loc(1)); %特征点数目
114. loct=loc;
115. for i=1:num
116. l1=[loc(i,1),loc(i,2),1];
117. l1_2=H*l1';
118. l1_2=l1_2/l1_2(3);
119. loct(i,1)=l1_2(1);
120. loct(i,2)=l1_2(2);
121. end
122. end
124. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
125. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
126. % 加载描述子文件
127. function [des, loc, dimdes]=loadFeatures(file)
128. fid = fopen(file, 'r');
129. [dimdes, count]=fscanf(fid, '%d %d',[1 2]);
130. if count ~= 2
131. error('Invalid keypoint file beginning.');
132. end
133. num = dimdes(1); %描述子数目
134. len = dimdes(2); %描述子长度
135. loc = double(zeros(num, 4));
136. des = double(zeros(num, len));
137. % 将描述子导入向量
138. for i = 1:num
139. [vector, count] = fscanf(fid, '%f %f %f %f', [1 len+5]); %row col scale ori
140. if count ~= (5+len)
141. error('Invalid keypoint file format');
142. end
143. loc(i, :) = vector(1, 1:4);
144. descrip(1, :) = vector(1, 6:len+5);
145. descrip = descrip / sqrt(sum(descrip.^2));
146. des(i, :) = descrip(1, :);
147. end
148. fclose(fid);
149. end
151. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
152. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
153. % 计算实际特征匹配的偏差位置
154. function matchLoc=matchLocal(loct,loc,ImageSize,matchFeaIndex)
155. matchLoc=ones(size(matchFeaIndex));
156. for i=1:size(loct)
157. if matchFeaIndex(i)~=0 %找到匹配的点
158. %注意仿射变换，位置可能超过图像范围
159. if loct(i,1)&gt;0&amp;&amp;loct(i,1)&lt;ImageSize(1)&amp;&amp;loct(i,2)&gt;0&amp;&amp;loct(i,2)&lt;ImageSize(2)...
160. &amp;&amp;loc(matchFeaIndex(i),1)&gt;0&amp;&amp;loc(matchFeaIndex(i),1)&lt;ImageSize(1)&amp;&amp;loc(matchFeaIndex(i),2)&gt;0&amp;&amp;loc(matchFeaIndex(i),2)&lt;ImageSize(2)
161. matchLoc(i)=(loct(i,1)-loc(matchFeaIndex(i),1))^2+(loct(i,2)-loc(matchFeaIndex(i),2))^2;
162. matchLoc(i)=sqrt(matchLoc(i)); %这个偏差即是欧氏距离
163. end
164. end
165. end
166. matchLoc=100.*matchLoc./max(matchLoc); %计算百分比
167. end