注:
1.特征点必须具有旋转不变性
如:①几何不变性: 位移、旋转、尺度,......
②光度不变性: 光照、 曝光, …...
1.2.特征匹配的基本流程
①根据特征准则提取图像中的特征点
②提取特征点周围的图像块,构造特征描述符
③通过特征描述符对比,实现特征匹配
一.Harris角点检测算法
- 角点:① 局部窗口沿各个方向移动,都产生明显变化的点; ②图像局部曲率突变的点
- 好的角点检测算法:①可以检测出图像中“真实的”角点 ; ② 准确的定位性能; ③很高的稳定性 ;④ 具有对噪声的鲁棒性; ⑤具有较高的计算效率
- Harris步骤:①.读入图像转灰度图;②计算响应函数;③基于响应值选择角点;④.在原图中画出检测的角点
1.Harris角点检测算法的基本思想:①.从图像局部的小窗口观察图像特征; ②.角点定义⬅窗口向任意方向的移动都导致图像灰度的明显变化
2.Harris算法的数学表达:
将图像窗口平移[u,v]产生灰度变化E(u,v),其中w(x,y) 是窗口函数,一般是高斯函数,所以可以把w(x,y)看做是高斯滤波器。I(x,y)是图像灰度, I(x+u,y+v)是平移后的图像灰度。
由:
得到:
对于局部微小的移动量 [u,v],近似表达为:


记M的特征值为,
角点响应函数R:
为了去除加权常数k,我们通常使用商数:
Harris角点检测算法就是对角点响应函数R进行阈值处理:R > threshold,即提取R的局部极大值
3.代码实现
im=array(Image.open('jmu.jpg').convert('L'))
figure()
gray()
subplot(321)
imshow(im)
title('原灰度图像')
axis('off')
axis('equal')
def compute_harris_response(im,sigma=3):
#计算导数
imx=zeros(im.shape)
#构造一个3×3的高斯滤波
gaussian_filter(im,(sigma,sigma),(0,1),imx)
imy=zeros(im.shape)
gaussian_filter(im,(sigma,sigma),(1,0),imy)
#计算Harris矩阵的分量
Wxx=gaussian_filter(imx*imx,sigma)
Wxy=gaussian_filter(imx*imy,sigma)
Wyy=gaussian_filter(imy*imy,sigma)
#计算特征值和迹
Wdet=Wxx*Wyy-Wxy*2
Wtr=Wxx+Wyy
return Wdet/Wtr
#提取角点
def get_harris_points(harrisim,min_dist=10,threshold=0.5):
#选取高于阈值的候选角点
corner_threshold=harrisim.max()*threshold
harrisim_t=(harrisim>corner_threshold)
#得到候选角点
coords=array(harrisim_t.nonzero()).T
#以及它们的harris响应值
candidate_values=[harrisim[c[0],c[1]] for c in coords]
#对harris响应值排序
index=argsort(candidate_values)
#将可行点保存到数组中
allowed_locations=zeros(harrisim.shape)
allowed_locations[min_dist:-min_dist,min_dist:-min_dist]=1
#按照min_dist原则,选择最佳角点
filtered_coords = []
for i in index:
if allowed_locations[coords[i,0],coords[i,1]] == 1:
filtered_coords.append(coords[i])
allowed_locations[(coords[i,0]-min_dist):(coords[i,0]+min_dist),
(coords[i,1]-min_dist):(coords[i,1]+min_dist)]=0
return filtered_coords
#检测harris角点
harrism=compute_harris_response(im)
#响应Harris函数
threshold=[0.01,0.05,0.1,0.5]
for i,thres in enumerate(threshold):
filtered_coords=get_harris_points(harrism,50,thres)
subplot(3,2,i+2)
title(str(threshold[i]))
imshow(im)
plot([p[1] for p in filtered_coords],[p[0] for p in filtered_coords],'*')
axis('off')
show()
3.1实验结果
二.SIFT(尺度不变特征变换)
SIFT特征是基于物体上的一些局部外观的兴趣点而与影像的大小和旋转无关。对于光线、噪声、些微视角改变的容忍度也相当高,所以用来检测要拼接图像的特征及关键点就很有优势
使用开源工具包 VLFeat 提供的二进制文件来计算图像的 SIFT 特征.将3个文件copy(sift.exe,vl.dll,vl.lib)到项目所在文件夹中.
2.1 步骤
- 提取关键点
- 对关键点附加局部特征,即描述符
- 通过特征点(附带上特征向量的关键点)的两两比较找出相互匹配的若干特征点建立景物间的对应关系
2.2代码
检查兴趣点
def process_image(imagename, resultname, params="--edge-thresh 10 --peak-thresh 5"):
""" 处理一幅图像,然后将结果保存在文件中"""
if imagename[-3:] != 'pgm':
#创建一个pgm文件
im = Image.open(imagename).convert('L')
im.save('tmp.pgm')
imagename ='tmp.pgm'
cmmd = str(r"D:\pycharm\实验\vlfeat\sift.exe "+imagename+" --output="+resultname+" "+params)
os.system(cmmd)
print ('processed', imagename, 'to', resultname)
def read_features_from_file(filename):
"""读取特征属性值,然后将其以矩阵的形式返回"""
f = loadtxt(filename)
return f[:,:4], f[:,4:] #特征位置,描述子
由于该二进制文件需要的图像格式为灰度.pgm,所以如果图像为其他各是,我们需要首先将其转换成.pgm格式文件。
转换生成的文本文件如下:
上面的数据前两列代表兴趣点的坐标,第三四列分别是兴趣点的尺度和方向角度。后面紧跟着对应描述符的128维向量。
将上面输出文件中的特征点读取到Numpy数组中的函数
def read_features_from_file(filename):
"""读取特征属性值,然后将其以矩阵的形式返回"""
f = loadtxt(filename)
return f[:,:4], f[:,4:] #特征位置,描述子
读取特征后,通过图像绘制出它们的位置,将其可视化。
def plot_features(im,locs,circle=False):
"""显示带有特征的图像
输入:im(数组图像),locs(每个特征的行、列、尺度和朝向)"""
def draw_circle(c,r):
t = arange(0,1.01,.01)*2*pi
x = r*cos(t) + c[0]
y = r*sin(t) + c[1]
plot(x, y, 'b', linewidth=2)
imshow(im)
if circle:
for p in locs:
draw_circle(p[:2], p[2])
else:
plot(locs[:,0], locs[:,1], 'ob')
axis('off')
2.3实验结果
可以看出,两个算法所选择特征点的位置不同
2.4 匹配描述子
对于将一幅图像中的特征匹配到另一幅图像的特征,一种稳健的准则(同样是由LOWE提出)是使用这两个特征距离和两个最匹配特征距离的比率。相对于图像中的其他特征,该准则保证能够找到足够相似的唯一特征
def plot_features(im,locs,circle=False):
"""显示带有特征的图像
输入:im(数组图像),locs(每个特征的行、列、尺度和朝向)"""
def draw_circle(c,r):
t = arange(0,1.01,.01)*2*pi
x = r*cos(t) + c[0]
y = r*sin(t) + c[1]
plot(x, y, 'b', linewidth=2)
imshow(im)
if circle:
for p in locs:
draw_circle(p[:2], p[2])
else:
plot(locs[:,0], locs[:,1], 'ob')
axis('off')
实验结果:
def match(desc1, desc2):
"""对于第一幅图像的每个描述子,选取其在第二幅图像中的匹配
输入:desc1(第一幅图像中的描述子),desc2(第二幅图像中的描述子)"""
desc1 = array([d/linalg.norm(d) for d in desc1])
desc2 = array([d/linalg.norm(d) for d in desc2])
dist_ratio = 0.6
desc1_size = desc1.shape
matchscores = zeros((desc1_size[0],1), 'int')
desc2t = desc2.T #预先计算矩阵转置
for i in range(desc1_size[0]):
dotprods = dot(desc1[i,:], desc2t) #向量点乘
dotprods = 0.9999*dotprods
# 反余弦和反排序,返回第二幅图像中特征的索引
index = argsort(arccos(dotprods))
# 检查最近邻的角度是否小于dist_ratio乘以第二近邻的角度
if arccos(dotprods)[index[0]] < dist_ratio * arccos(dotprods)[index[1]]:
matchscores[i] = int(index[0])
return matchscores
def match_twosided(desc1,decs2):
"""双向对称版本的match"""
matches_12 = match(desc1, decs2)
matches_21 = match(decs2, decs2)
ndx_12 = matches_12.nonzero()[0]
# 去除不对称匹配
for n in ndx_12:
if matches_21[int(matches_12[n])] != n:
matches_12[n] = 0
return matches_12
def appendimages(im1, im2):
"""返回将两幅图像并排拼接成的一幅新图像"""
# 选取具有最少行数的图像,然后填充足够的空行
row1 = im1.shape[0]
row2 = im2.shape[0]
if row1 < row2:
im1 = concatenate((im1,zeros((row2-row1,im1.shape[1]))), axis=0)
elif row1 > row2:
im2 = concatenate((im2,zeros((row1-row2,im2.shape[1]))), axis=0)
# 如果这些情况都没有,那么他们的行数相同,不需要进行填充
return concatenate((im1,im2), axis=1)
def plot_matches(im1, im2, locs1, locs2, matchscores, show_below=True):
"""显示一幅带有连接匹配之间连线的图片
输入:im1,im2(数组图像),locs1,locs2(特征位置),matchscores(match的输出),
show_below(如果图像应该显示再匹配下方)"""
im3 = appendimages(im1,im2)
if show_below:
im3 = vstack((im3,im3))
imshow(im3)
cols1 = im1.shape[1]
for i in range(len(matchscores)):
if matchscores[i] > 0:
plot([locs1[i, 0], locs2[matchscores[i, 0], 0] + cols1], [locs1[i, 1], locs2[matchscores[i, 0], 1]], 'c')
axis('off')
使用两张图像,匹配得到结果: