特征检测与匹配
1 角点检测算法实验
1.1 实验目的与要求
（1）了解及掌握角点检测算法原理。
（2）掌握在MATLAB中角点算法的编程。
（3）掌握Moravec，Harris与SUSAN算法的差异。
1.2 实验原理及知识点
（1）Harris算法原理
Harris角点检测算法的基本原理是取以目标像素点为中心的一个小窗口，计算窗口沿任何方向移动后的灰度变化，并用解析形式表达。设以像素点(x,y)为中心的小窗口在X方向上移动u，y方向上移动v，Harris给出了灰度变化度量的解析表达式：
其中，R为旋转因子，对角化处理后并不改变以u,v为坐标参数的空间曲面的形状，其特征值反应了两个主轴方向的图像表面曲率。当两个特征值均较小时，表明目标点附近区域为“平坦区域”；特征值一大一小时，表明特征点位于“边缘”上；只有当两个特征值均比较大时，沿任何方向的移动均将导致灰度的剧烈变化。Harris的角点响应函数(CRF)表达式由此而得到：
其中：det(M)表示矩阵M的行列式，trace(M)表示矩阵的迹。当目标像素点的CRF值大于给定的阈值时，该像素点即为角点。
（2）SUSAN算法原理
该算法直接利用像素的灰度进行角点检测，而不考虑曲率等复杂的角点特征。SUSAN检测算子的基本原理是通过统计某一像素局部区域内与该像素等灰度值近似的的点的个数，实现角点检测。
该算法一般利用一个37像素的圆形模板来实现的，下图所示。图1中圆形模板e的圆心称为核心，假如模板中的某些像素的亮度与核心相同或相似，就定义这些像素组成的区域为USAN（核值相似区）区域。图2显示出了不同位置的USAN区域面积大小。USAN区域包含了图像结构的以下信息：在a位置，核心点在角点上，USAN面积达到最小:在b位置，核心点在边缘线上时，USAN区域面积接近最大值的一半；在 c、d位置，核心点处于黑色矩形区域之内，USAN区域面积接近最大值。因此，可以根据USAN区的面积大小检测出角点。
具体检测时，是用圆形模板扫描整个图像，比较模板内每一像素与中心像素的灰度值。并给定阈值来判别该像素是否属于USAN区域，式2.1是SUSAN算法的原始相似比较函数。式2.2是在实际应用中比较常用的相似比较函数 ：
用于计算以每个像素点为核心的USAN区的像素个数；是模板中心像素(核)的灰度值；I(r)为模板内其他任意像素的灰度值；t是区分特征目标与的一个重要阈值。
图像中某一点USAN区域大小可由下式2.3表示：
式中，g是抑制噪声的几何阈值 ，它决定了输出角点的USAN区域的最大值。同时它还决定了所检测到的角点的尖锐程度。g取得越小，所检测到的角点越尖锐。用这种原理，取不同的几何门限，不但能检测角点．还可以检测交点、边缘等特征。
1.3 实验内容及步骤
1、参考以下代码，也可自行编程使用Harris算法对lena.png和cameraman.tif图片(或其它图片)进行角点检测。
参考代码：
%%%% Harris角点检测算法 Matlab code
clear all; clc ;tic;
ori_im = imread('C:\Users\33841\Desktop\image\lena.png'); % 读取图像
if(size(ori_im,3)==3)ori_im = rgb2gray(uint8(ori_im)); %转为灰度图像
end
% fx = [5 0 -5;8 0 -8;5 0 -5]; % 高斯函数一阶微分，x方向(用于改进的Harris角点提取算法)
fx = [-2 -1 0 1 2]; % x方向梯度算子(用于Harris角点提取算法)
Ix = filter2(fx,ori_im); % x方向滤波
% fy = [5 8 5;0 0 0;-5 -8 -5]; % 高斯函数一阶微分，y方向(用于改进的Harris角点提取算法)
fy = [-2;-1;0;1;2]; % y方向梯度算子(用于Harris角点提取算法)
Iy = filter2(fy,ori_im); % y方向滤波
Ix2 = Ix.^2;
Iy2 = Iy.^2;
Ixy = Ix.*Iy;
clear Ix;
clear Iy;
h= fspecial('gaussian',[7 7],2); % 产生7*7的高斯窗函数，sigma=2
Ix2 = filter2(h,Ix2);
Iy2 = filter2(h,Iy2);
Ixy = filter2(h,Ixy);
height = size(ori_im,1);
width = size(ori_im,2);
result = zeros(height,width); % 纪录角点位置，角点处值为1
R = zeros(height,width);
for i = 1:heightfor j = 1:widthM = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)]; % auto correlation matrixR(i,j) = det(M)-0.06*(trace(M))^2; end
end
cnt = 0;
for i = 2:height-1for j = 2:width-1% 进行非极大抑制，窗口大小3*3if R(i,j) > R(i-1,j-1) && R(i,j) > R(i-1,j) && R(i,j) > R(i-1,j+1) && R(i,j) > R(i,j-1) && R(i,j) > R(i,j+1) && R(i,j) > R(i+1,j-1) && R(i,j) > R(i+1,j) && R(i,j) > R(i+1,j+1)result(i,j) = 1;cnt = cnt+1;endend
end
Rsort=zeros(cnt,1);
[posr, posc] = find(result == 1);
for i=1:cntRsort(i)=R(posr(i),posc(i));
end
[Rsort,ix]=sort(Rsort,1);
Rsort=flipud(Rsort);
ix=flipud(ix);
ps=100;
posr2=zeros(ps,1);
posc2=zeros(ps,1);
for i=1:psposr2(i)=posr(ix(i));posc2(i)=posc(ix(i));
end
imshow(ori_im);
hold on;
plot(posc2,posr2,'g+');
（2）参考以下代码，编程使用SUSAN算法对lena.png和 cameraman.tif图片(或其它图片)进行角点检测,对结果和harris算法实现的结果进行对比,分析它们之间的不同之处。 进一步自行实现Moravec算法，比较三种角点检测算法的差异。
% detectSUSANFeayures.m Author：Ephemeroptera Date:2018\11\12
%7x7 圆周template的局部地址
% [-1,-3] [0,-3] [1,-3]
% [-2,-2] [-1,-2] [0,-2] [1,-2] [2,-2]
% [-3,-1] [-2,-1] [-1,-1] [0,-1] [1,-1] [2,-1] [3, -1]
% [-3, 0] [-2, 0] [-1, 0] [0, 0] [1, 0] [2, 0] [3 , 0]
% [-3 ,1] [-2, 1] [-1, 1] [0, 1] [1, 1] [2 ,1] [3 , 1]
% [-2, 2] [-1, 2] [0, 2] [1 ,2] [2 ,2]
% [-1 ,3] [0, 3] [1, 3]
%
%formula：c(r,r0)= 1 if |I(r)-I(r0)|<=t;
% 0 if |I(r)-I(r0)|>t;
% USAN(r0)=Σc(r,r0);%设置圆周模板半径和滑动窗口的步长
radius=3;Xstep=1;Ystep=1;
template=fspecial('disk',radius);
template(template>0.01)=1; %模板二值化
template(template<0.01)=0;%提取圆周模板的逻辑地址
[tem_x,tem_y]=find(template==1);
tem_x=tem_x-radius-1;
tem_y=tem_y-radius-1;t=45; %USAN判定阈值
I = imread('cameraman.tif');
W=size(I,2);H=size(I,1); %图像大小
nucleas_X=radius+1:Xstep:W-radius; %模板圆心即nucleas运动范围
nucleas_Y=radius+1:Ystep:H-radius;
USAN=zeros(size(I,1),size(I,2)); %初始化USAN累加器
%网格遍历
tic;
for y=nucleas_Yfor x=nucleas_Xfor e=1:length(tem_x) %圆周模板上进行判定delta=I(y+tem_y(e),x+tem_x(e))-I(y,x);if delta<tUSAN(y,x)=USAN(y,x)+1; %低于阈值则收纳endend fprintf(strcat('已处理第','(',num2str(y),',',num2str(x),')','像素点..\n'));end
end
toc;
%边缘检测
%
%formula： R(r0) = g-USAN(r0) , if USAN(r0)<g
% 0 , otherwise
% where g=3/4*max(USAN(:)), USAN越小边缘响应越强,其中角点极小，对应R极大
g=1/2*max(USAN(:));
R=zeros(H,W);
for i=1:size(USAN,1)for j=1:size(USAN,2)if USAN(i,j)<gR(i,j)=g-USAN(i,j);elseR(i,j)=0;endend
end
BIN=zeros(H,W);
BIN(R>0)=1;
figure(1),imshow(BIN,[]),title('SUSAN 边缘检测');%角点检测（非极大值抑制）
corners=[];
for i=2:H-1for j=2:W-1 if R(i,j)>max([max(R(i-1,j-1:j+1)),R(i,j-1) ,R(i,j+1), max(R(i+1,j-1:j+1))])corners=[corners;[i,j]];endend
end
% I = insertMarker(I,corners);
figure(2);imshow(I);hold on;
set(gca,'xaxislocation','top','yaxislocation','left','ydir','reverse');
scatter(corners(:,2),corners(:,1),'x','g'),title('SUSAN 角点检测');
hold off;
Harris角点检测算法：
Harris算法通过计算像素周围窗口的自相关矩阵来检测角点。它考虑了像素周围的小窗口，通过计算窗口内像素的梯度和方向来确定角点的位置。
Harris算法具有旋转不变性和尺度不变性，可以在不同尺度和方向上检测角点。
该算法对于纹理丰富和噪声较多的图像具有较好的性能，并且可以检测到较为精细的角点。
SUSAN算法：
SUSAN算法基于圆形模板进行角点检测。它比较模板中心像素与模板内其他像素的灰度差异，通过统计差异值来确定角点的位置。
SUSAN算法具有较快的计算速度，并且对于边缘和噪声的干扰较小。
该算法可以检测到较为明显的角点，但对于一些细微的角点可能无法准确检测。
Moravec算法：
Moravec算法是最早的角点检测算法之一。它通过计算像素周围窗口在不同方向上的强度变化来检测角点。
Moravec算法对于明显的角点具有较好的检测效果，但对于噪声和纹理丰富的图像可能会出现误检和漏检的情况。
该算法的计算复杂度较高，对于一些实时性要求较高的应用场景可能不够适用。
Harris算法、SUSAN算法和Moravec算法在角点检测中各具特点。Harris算法具有较好的旋转不变性和尺度不变性，适用于纹理丰富和噪声较多的图像；SUSAN算法具有较快的计算速度和较好的边缘抑制能力，适用于实时性要求较高的场景；而Moravec算法作为最早的角点检测算法之一，对于一些明显的角点具有较好的检测效果，但对于复杂的图像可能会出现一些问题。在具体应用中，需要根据图像的特点和需求来选择合适的算法进行角点检测。
1.4 实验报告要求
描述实验的基本步骤，用数据和图片给出各个步骤中取得的实验结果和源代码，并进行必要的讨论，必须包括原始图像及其计算/处理后的图像。
可使用不同平台实现同样的功能，matlab，Python+Opencv 或者 C + Opencv都可以实现类似的读取图片和显示操作。
2 SIFT特征提取与匹配实验（选做）
2.1 实验目的
掌握SIFT特征提取与匹配的相关原理，懂的如何利用相关编程语言实现，并进行相关分析和调整参数。
2.2 实验原理
SIFT算法是一种提取局部特征的算法，在尺度空间寻找极值点,并确定关键点(Key points)的位置和关键点所处的尺度;然后使用关键点邻域梯度的主方向作为该点的方向特征，以实现算子对尺度和方向的无关性。
2.3 实验内容
准备两张有部分区域重叠的照片，使用opencv的特征提取、匹配函数实现两张图片的SIFT特征提取和匹配，并绘制特征点在匹配前和匹配后的连线。(扩展:将两幅图片进行全景拼接)
参考代码:
SIFT特征提取与匹配
#include "highgui/highgui.hpp"
#include "opencv2/nonfree/nonfree.hpp"
#include "opencv2/legacy/legacy.hpp" using namespace cv;int main(int argc, char *argv[])
{Mat image01 = imread("1.jpg");Mat image02 = imread("2.jpg");Mat image1, image2;GaussianBlur(image01, image1, Size(3, 4), 0.5);GaussianBlur(image02, image2, Size(3, 4), 0.5);SiftFeatureDetector siftDetector(30); vector<KeyPoint> keyPoint1, keyPoint2;siftDetector.detect(image1, keyPoint1); siftDetector.detect(image2, keyPoint2);drawKeypoints(image1, keyPoint1, image1, Scalar::all(-1), DrawMatchesFlags::DRAW_RICH_KEYPOINTS);drawKeypoints(image2, keyPoint2, image2, Scalar::all(-1), DrawMatchesFlags::DRAW_RICH_KEYPOINTS);namedWindow("KeyPoints of image1", 0);namedWindow("KeyPoints of image2", 0);imshow("KeyPoints of image1", image1);imshow("KeyPoints of image2", image2);SiftDescriptorExtractor siftDescriptor;Mat imageDesc1, imageDesc2;siftDescriptor.compute(image1, keyPoint1, imageDesc1);siftDescriptor.compute(image2, keyPoint2, imageDesc2);BruteForceMatcher<L2<float>> matcher;vector<DMatch> matchePoints;matcher.match(imageDesc1, imageDesc2, matchePoints, Mat());Mat imageOutput;drawMatches(image01, keyPoint1, image02, keyPoint2, matchePoints, imageOutput);namedWindow("out", 0);imshow("out", imageOutput);waitKey();return 0;
}3.Harris特征检测与匹配实验的Python实现 (选做)
特征检测与匹配的目标是识别一个图像中的关键点与另一个图像中的对应点之间的配对。在此实验中，你将编写代码以检测图像中的特征点（对于平移、旋转和照明具有一定的不变性），并在另一个图像中找到最佳匹配特征。
为了帮你可视化结果并调试程序，我们提供了一个用户界面，可以显示检测到的特征和最佳匹配。我们还提供了一个示例ORB特征检测器，用于结果比较。
3.1. 实施细节
该实验有三个部分：特征检测、特征描述和特征匹配。您所需要实现的所有代码都在features.py中。
3.1.1 特征检测
你将使用Harris角点检测方法识别图像中的关键点，Harris角点检测的详细信息请参阅课程ppt和论文（harris.pdf）。对于图像中的每个点，考虑该点周围的像素窗口，计算该点的Harris矩阵H，定义为
我们还需要每个像素的度数方向。以梯度的方向作为近似方向。零角度指向右侧，正角度为逆时针方向。注意：不要通过结构张量的特征向量分析来计算方向。
我们将根据c(H)选择最强的关键点，它们是7×7邻域中的局部最大值。
你需要补充完成类HarrisKeypointDetector中的三个函数，computeHarrisValues（TODO 1，为图像中每个像素计算Harris强度函数与方向）、computeLocalMaxima（TODO 2，计算布尔矩阵，指示每个像素是否是局部最大值）、detectKeypoints（TODO 3，根据像素的Harris强度与是否局部最大值生成特征点集合）。
可以使用以下库函数：
- scipy.ndimage.sobel：使用Sobel算子对输入图像滤波。
- scipy.ndimage.gaussian_filter：使用高斯卷积核对输入图像滤波。
- np.arctan2
- scipy.ndimage.filters.maximum_filter：使用最大过滤器过滤输入图像。
- scipy.ndimage.filters.convolve：使用选定的滤波器对输入图像滤波。
3.1.2 特征描述
下一步是为每个特征点提供一个描述符，此描述符将用于比较不同图像中的要素以查看它们是否匹配。
你将实现两个特征描述符，SimpleFeatureDescriptor类和MOPSFeaturesDescriptor类。这些类的describeFeatures函数获取一组特征点的位置和方向信息，并计算这些特征点的描述符，然后将这些描述符存储在一个numpy二维数组中，该二维数组的行数为特征点的个数，列数为该特征描述符的大小（例如，对于5×5简单特征描述符，为25）。
作为初学者，你将首先实现一个简单的描述符SimpleFeatureDescriptor（TODO 4），它是5×5邻域中的像素强度值（即灰度值）。当所比较的图像是平移关系时应该可以正常工作。
其次，你将实现MOPS描述符的简化版本MOPSFeaturesDescriptor（TODO 5）。MOPS的详情可参考课程ppt和论文（MOPS.pdf）。你将计算从特征点周围的40×40像素区域子采样的8×8定向图像块。你必须提供一个变换矩阵，它将围绕特征点的40×40窗口变换为8×8图像块，使其特征点方向指向右侧。你将使用cv2.warpAffine函数实现变换。warpAffine采用2×3前向卷绕的仿射矩阵，它从左边乘以原坐标使得变换后的坐标为列向量。生成2×3变换矩阵的最简单方法是组合多个转换：平移（T1）、旋转（R）、缩放（S）和平移（T2）。变换矩阵是矩阵乘积T2×S×R×T1，下图说明了该序列。变换与卷绕的详情可参考论文，注意这些变换采用齐次坐标。你可能会发现transformations.py中的函数很有用，它们实现了3D仿射变换（你需要手动将它们转换为2D变换）。
你还应该将8×8图像块规范化为零均值和单位方差（TODO 6）。如果方差非常接近零（幅度小于10 -5），那么你应该返回一个全零向量以避免除零错误。你可能会用到np.std、np.mean函数。
3.1.3 特征匹配
下一步是编写匹配它们的代码（即，在一个图像中给定一个特征，在另一个图像中找到最佳匹配特征）。最简单的方法如下：比较两个特征并计算它们之间的标量距离，最佳匹配是距离最小的特征。你将实现两个距离函数：
- 平方差之和（SSD）：这是两个特征向量之间的欧氏距离平方。
- 比率测试：通过SSD距离查找最近和次近的两个特征，比率测试距离是它们的比率（即，最接近的特征匹配的SSD距离除以第二接近的特征匹配的SSD距离）。
你将实现SSDFeatureMatcher（TODO 7）和RatioFeatureMatcher（TODO 8）的matchFeatures函数，该函数返回cv2.DMatch对象的列表。你应该将queryIdx属性设置为第一个图像中特征的索引，将trainIdx属性设置为第二个图像中特征的索引，将距离属性设置为两个特征之间的距离（例如，SSD或比率测试）。你可能会用到scipy.spatial.distance.cdist和numpy.argmin两个库函数。
3.2 编码规则
可以使用NumPy、SciPy和OpenCV2函数来实现数学、滤波和变换操作。不要直接使用关键点检测或特征匹配的函数。
使用Sobel算子或高斯滤波器时，应使用“reflect”模式，该模式在边缘处给出零梯度。
3.3 测试和可视化
可以通过运行tests.py文件来测试你的TODO代码。这将加载一个小图像、运行你的代码并与正确的输出进行比较。这将允许你以递增的方式测试代码，而无需完成所有TODO块。
我们已经为TODO 1-6提供了测试代码。你应该为TODO 7和8设计自己的测试。最后，请注意提供的测试非常简单——它们可以帮助你入门。但是通过提供的测试用例并不意味着将通过评分测试用例。
通过运行featuresUI.py，你将看到一个UI，你可以进行特征点检测、特征匹配和基准测试。UI是一种工具，可帮助你可视化特征点检测和匹配结果。请记住，你的代码将以数字方式进行评分，而不是以可视方式进行评分。
我们提供了一组基准图像，用于测试算法的性能，作为不同类型的受控变化（即旋转，缩放，照明，透视，模糊）的函数。对于这些图像中的每一个，我们都知道正确的变换，因此可以衡量每个特征匹配的准确性。这是使用我们在框架代码中提供的例程来完成的。你还应该出去拍摄自己的照片，看看你的方法在更有趣的数据集上的效果如何。例如，你可以拍摄几个不同对象（例如，书籍，办公室，建筑物等）的图像，并查看它的工作情况。
3.4 所提供的文件
features.py：需要你实现的函数都在这个文件中。
tests.py：可以用该文件测试TODO1-6。
featuresUI.py：提供了一个GUI，可用来调试算法，并获得benchmark图像。
transformations.py：提供了3D仿射变换的矩阵获取函数。
resources目录下是可用的测试图片。
harris.pdf：Harris角点检测算法论文。
MOPS.pdf：MOPS算法论文。
3.5 实验环境配置
先在Linux或Windows上安装Python3，之后安装Numpy、Scipy、OpenCV for python。如果要运行featuresUI.py，还要安装PIL（Pillow）。
3.6 提交文件
features.py：内含需要完成的TODO1-8。
harris.png：每次计算Harris角点强度，都会自动生成一个harris.png文件，用resources目录下的yosemite1.jpg来生成此harris.png文件。
report(.docx,.pdf)：报告中应包括resources目录下的yosemite数据集的基准测试结果中的ROC曲线和AUC。你可以通过运行featuresUI.py，切换到“Benchmark”选项卡，按“Run Benchmark”，选择目录“resources/ yosemite”并等待一会儿来获得ROC曲线，同时AUC将显示在屏幕底部，可以通过按“Screenshot”保存ROC曲线。记录参数，如Threshold，并使用ROC曲线和AUC来比较几种特征检测、描述及匹配方法。
3.7 额外加分
三种方式可以获得额外加分：
- 在实现MOPSFeatureDescriptor类的describeFeatures方法时，加入尺度不变特性。一种方法可以参考MOPS论文（MOPS.pdf），你也可以设计其它方法。
- 实现MOPS论文（MOPS.pdf）中的自适应非最大抑制。
- 设计并实现其它的特征检测与描述方法，或对已有算法进行改进。报告你算法的性能，并分析提高性能的原因。如果最终的比率测试AUC提高15%（AUC提高15%的意思是(1-AUC)减小15%）以上，将被认可为实质性改进。
如果你希望有额外加分，除提交代码外，在报告中必须详细描述你的改进算法及性能提升结果。可以在5个数据集（bikes、graf、leuven、wall、yosemite）上运行UI基准测量你的平均AUC。
简单的超参数更改将获得很小的额外加分，实质性的算法改进更容易获得认可。
在尝试额外加分之前，你首先应该完成基本任务。基本任务比额外加分更容易得分。
3.8 调试中的注意事项
工程由python2.7转化而来，并在python3.6下经过测试。
以下三个问题是opencv的问题，不一定会被触发，哪怕大家用的都是python3也取决于你具体安装的opencv版本，解决方法如下：
1、所有路径中均不允许存在中文，否则会报出NoneType错误；
2、如果你的ORB不起作用并且报错——Incorrect type of self (must be 'Feature2D' or its derivative)，请将features.py文件中的orb = cv2.ORB()替换为：orb = cv2.ORB_create()；反之亦然。
3、如果features.py的以下代码报出错误——OpenCV error: (-215) K == 1 && update == 0 && mask.empty() in function cv::batchDistance，则进行调整：
class ORBFeatureMatcher(FeatureMatcher):
## ...(略)...
def matchFeatures(self, desc1, desc2):
return self.bf.match(desc1.astype(np.uint8), desc2.astype(np.uint8), k=1) # 移除或修改这个k的值，让它直接通过断言。具体值取决于你的OpenCV和底层C语言库的版本。