CS143-project2 局部特征匹配 Local Feature Matching

本项目通过harris算法对图像进行兴趣点的检测和提取。然后计算兴趣点周围4*4区域的8个方向特征向量，最后对获取到的特征向量对进行匹配。本算法对旋转，缩放，亮度变化保持一定程度的稳定性。

环境：

操作系统：WIN7

实现平台：MATLAB2014a

 

 本项目主要包括以下4部分：

1. 特征点检测和提取(harris).
2. 生成特征算子(sift).
3. 图像匹配(match).
4. 实验DEMO

一.特征点检测和提取(harris)

Harris 角点检测算法是由Chris Harris 和MikeStephens 在1988 年提出， 该算法是在Moravec 算法的基础上发展起来的。Moravec 算法是研究图像中一个局部窗口在不同方向进行少量的偏移后， 考察窗口内图像亮度值的平均变化。

需要考虑下面三种情况：
（1）如果窗口内区域图像的亮度值恒定，那么所有不同方向的偏移几乎不发生变化；
（2）如果窗口跨越一条边，那么沿着这条边的偏移几乎不发生变化， 但是与边垂直的偏移会发生很大的变化；
（3）如果窗口包含一个孤立的点或者角点，那么所有不同方向的偏移会发生很大的变化。
Harris检测方法考虑的是用一个高斯窗或矩形窗在图像上移动,由模板窗口取得原图像衍生出2 ×2的局部结构M 矩阵.对该模板矩阵求取特征值λ1 和λ2 , 建立度量函数R = detM - k(traceM)2,detM =λ1λ2 ,traceM =λ1+λ2, 
其中，Det(M)是矩阵M 的行列式值，Trace(M)是矩阵M 的迹，k=0.04。根据R是否大于0即可判断该点是否是角点。值得注意的是该方法具有旋转不变性
,但检测的角点有较大的冗余,需要根据实际经验来确定R的阈值。

自相关矩阵的计算和R值计算部分代码如下：

for i = 1:row
 for j = 1:col
  M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)]; %自相关矩阵
  R(i,j) = det(M) - k*(trace(M))^2; %根据公式计算R值
 end;
end;
Rmax = max(max(R));

实验得出，获取到的特征点与设定的搜索窗口密切相关。窗口越小，获取到的特征点越多，反之越少。当搜索窗口feature_width =16时，获取到的角点如下：

 

二.生成特征算子(sift)

SIFT算法是本项目的重点及难点，它在空间尺度中寻找极值点，并提取出其位置、尺度、旋转不变量。
Lowe 将SIFT 算法分解为如下四步：
1. 尺度空间极值检测：搜索所有尺度上的图像位置。通过高斯微分函数来识别潜在的对于尺度和旋转不变的兴趣点。
2. 关键点定位：在每个候选的位置上，通过一个拟合精细的模型来确定位置和尺度。关键点的选择依据于它们的稳定程度。
3. 方向确定：基于图像局部的梯度方向，分配给每个关键点位置一个或多个方向。所有后面的对图像数据的操作都相对于关键点的方向、尺度和位置进行变换，从而提供对于这些变换的不变性。
4. 关键点描述：在每个关键点周围的邻域内，在选定的尺度上测量图像局部的梯度。这些梯度被变换成一种表示，这种表示允许比较大的局部形状的变形和光照变化。

特殊说明：

由于之前已经使用了harris对关键点进行了提取，所以对SIFT算法进行简化：直接使用harris获取到的关键点进行第3、4步计算。不再构建尺度空间检测和定位关键点。当然，这对尺度不变性的稳定性会有一些影响。

下面是具体是做法
二.1使用梯度直方图统计每个关键点的主方向

为关键点分配方向，目的是为了保证图像的旋转不变性。 具体实现是围绕关键点，选取一个窗口，窗口内计算各采样点的梯度和幅值，并统计360度（每10度一个直方图）的峰值，作为主方向。

窗口大小与关键点所在的空间尺度有关系，由于本项目没有构建尺度空间，所以窗口大小由多次试验后取效果较好的大小。本项目也设定了一个比例阀值，来表示达到直方图峰值的此比例，就列为关键点的辅方向，以便提高匹配时的鲁棒性。阀值默认是80%，即有一个方向超过峰值的80%的时候，就复制此关键点，并赋予此方向作为主方向参与后续的匹配。

以下是部分实现代码：

W = 8; %floor(3.0* sigmaw); %高斯加权窗口直径,取半径为4 调整
histo = zeros(1, NBINS); %累积之前初始化直方图
for xs = max(xp - W,1): min((N - 1), xp + W)
for ys = max(yp - W,1) : min((M-1), yp + W)
dx = (xs - xp);
dy = (ys - yp);
if dx^2 + dy^2 <= W^2 % 点在高斯加权圆内
wincoef = exp(-(dx^2 + dy^2)/(2*sigmaw^2)); %高斯窗口函数
bin = round( NBINS * angles(ys, xs)/(2*pi) + 0.5); %确定所在方向直方图的柱

histo(bin) = histo(bin) + wincoef * magnitudes(ys, xs); %用高斯加权后的向量幅度做直方图值累加
end

end
end
theta_max = max(histo); %找到直方图峰值
theta_indx = find(histo == theta_max); %关键点方向索引，大于80%峰值的角度

for i = 1: size(theta_indx, 2)
theta = 2*pi * theta_indx(i) / NBINS;
key_points = cat(1,key_points,[xp,yp,theta]); %幅值峰值大于80%的关键点，提高匹配的鲁棒性 
end

当阀值为1时（即只取最大峰值的方向），得到的中间结果是:

 

二.2构建4*4*8的128维特征向量

通过以上步骤，对于每一个关键点，拥有两个信息：位置、方向。接下来就是为每个关键点建立一个描述符，用一组向量将这个关键点描述出来，使其不随各种变化而改变，比如光照变化、视角变化等等。这个描述子不但包括关键点，也包含关键点周围对其有贡献的像素点，并且描述符应该有较高的独特性，以便于提高特征点正确匹配的概率。
具体实现是将关键点附近的领域划分为d ×d (Lowe 建议d=4)个子区域，每个子区域做为一个种子点，每个种子点有8 个方向。每个方向45度。每个区域大小选取与分配方向时一样，通过多次试验，取效果较好的尺寸。本项目是取的16。

计算时必须将坐标轴旋转为关键点的方向，以确保旋转不变性，如图：

旋转后领域内采样点的新坐标为：

 

其中为关键点的方向。

然后将领域内的采样点分配到对应的子区域内，将子区域内的梯度值分配到8 个方向上，计算其权值。遍历所有采样点，得到4*4*8维的一个特征算子。为使得其对亮度变化不敏感，对其进行归一化处理。设定个分量的阀值为0.2，若大于此阀值，再进行归一化。综上，即可得到描述关键点是描述子。

代码详见get_features.m

 

三.图像匹配(match).

一般采用各种距离函数作为特征的相似性度量,如欧氏距离、马氏距离等。通过相似性度量得到图像间的潜在匹配。本项目采用欧氏距离来计算两特征向量的相似度。取图像1中的某个关键点J，并找出其与图像2中欧式距离最近的关键点I和次近的关键点。如果最近的距离除以次近距离小于比例阈值k(经验值在0.4~0.6内效果较好)，则接受这一对匹配点J和I，否则J没有匹配点。降低这个比例阈值，SIFT匹配点数目会减少，但更加稳定。 
以下是调用evaluate_correspondence对Notre Dame图像对的试验结果

序号k正确匹配数量错误匹配数量10.44020.515030.639240.7791250.813751

可观察到，当k值越大，得到的匹配数量对越多，但是稳定性也更差。最顶端的图像对，就是当k=0.6时的结果

　

四.实验DEMO.

1.尺度不变性测试（美的集团总部大楼）：选取两幅角度不同，尺寸不同的美的总部大楼图像，匹配阀值k取0.4，进行测试，匹配效果较好

 

2.亮度、旋转不变性测试（书本）：用手机拍了两幅背景、亮度不同，且旋转了大概180度的书本图片，匹配阀值k取0.5，进行测试，匹配效果较好

 

 

提示与说明

1.对SIFT算法进行简化：直接使用harris获取到的关键点进行计算。不构建尺度空间检测和定位关键点

2.高斯模糊参数通过试验取效果最好的值

3.对于其他图像对，由于没有评价代码，我自己写了个test函数，用来连接两两特征匹配对，以观察其匹配效果

4.计算图像梯度方向时公式是使用反正切arctan函数，但是在实现中此函数的值域是(-pi/2,pi/2)，无法满足构建36个bin的要求，导致在测试旋转不变性时效果很差。

所以改使用atan2函数

5.get_features函数增加输出参数key_points，以方便在图像中显示观察