目录
1、SIFT（Scale Invariant Feature Transform）简介
一、SIFT原理及简介
1、SIFT（Scale Invariant Feature Transform）简介
传统的图像匹配算法一般是利用图像的角点与边缘进行匹配，鲁棒性较差，在光照变换，角度变换的时候容易出现误差。
SIFT（尺度不变特征变换）算法是有David G.Lowe教授在1999年提出，2004年加以完善。SIFT算法在光照不同，角度不同，尺度不同的情况下应用良好。
1.1、SIFT算法的操作步骤
兴趣点的检测：
在尺度空间（第二节有详细介绍）上，利用高斯微分函数，识别尺度和旋转不变的兴趣点。
特征方向的赋值：
基于特征点的梯度方向，给特征点分配8个方向，选取主方向，与辅方向（在直方图里等于主方向80%）
特征点描述：
把特征点的梯度方向转换成特定表示。
1.2、SIFT算法的适用范围
1.目标图像的旋转，平移，尺度变换。
2.反光图像
3.目标有障碍物
2、尺度空间
在不同的尺度下的同一物体计算机不会认为是同一物体，因此，需要将物体的不同尺度提供给计算机。
2.1、传统图像金字塔
传统的图像金字塔是对图像进行简单的缩放，放大等操作。这种做法简单粗暴，但容易造成采样信息的丢失，降低匹配结果的精确度。
2.2、高斯金字塔
1.对图像做高斯平滑
2.对平滑后得图像做采样
2.3、高斯模糊
高斯模糊的目的是为了拟合人类视觉中对远近不同图像所产生的尺度差异与清晰度差异。
经过高斯模糊的图像会体现出总体模糊化，边缘清晰化的特征
高斯模糊的具体做法是：
L(x,y,σ)=G(x,y,σ)∗I(x,y)
L(x,y,σ)=G(x,y,σ)∗I(x,y)
其中，GG是高斯函数：
G(x,y,σ)=12πσ2ex2+y22σ2
G(x,y,σ)=12πσ2ex2+y22σ2
其中，σσ是尺度空间因子，是高斯正态分布的标准差，反映了图像被模糊的程度，其值越大图像越模糊，对应的尺度也就越大，L(x,y,σ)L(x,y,σ)对应高斯尺度空间。
---------------------
作者：呆呆的猫
来源：优快云
原文：https://blog.youkuaiyun.com/jiaoyangwm/article/details/79986729
3.DOG空间极值检测
DOG层是由两个相邻的高斯空间图像层相减得到的，其目的是简化计算复杂度。
4.删除边缘响应点
一些噪声点会产生边缘响应，会影响匹配结果。可以通过尺度空间DoG函数进行曲线拟合寻找极值点，去除边缘响应点。
二、SIFT算法的简单应用与和Harris角点匹配的对比
这次采用的实验原图片是集美大学的中山纪念馆。
首先我们任意选取两张图片分别进行Harris特征匹配处理与SIFT特征匹配处理的对比
1.Harris特征匹配处理
Harris 可以表示出图像的角点，具体实现代码如下
from scipy.ndimage import filters
from numpy import *
from PIL import Image
import harris
from pylab import *
import os
im=array(Image.open('1.jpg').convert('L'))
harrisim=harris.compute_harris_response(im)
filtered_coords=harris.get_harris_points(harrisim,6)
harris.plot_harris_points(im,filtered_coords)
其中用到的harris函数实现方法如下：
from scipy.ndimage import filters
from numpy import *
from pylab import *
def compute_harris_response(im,sigma=3):
# derivatives
imx = zeros(im.shape)
filters.gaussian_filter(im, (sigma,sigma), (0,1), imx)
imy = zeros(im.shape)
filters.gaussian_filter(im, (sigma,sigma), (1,0), imy)
# compute components of the Harris matrix
Wxx = filters.gaussian_filter(imx*imx,sigma)
Wxy = filters.gaussian_filter(imx*imy,sigma)
Wyy = filters.gaussian_filter(imy*imy,sigma)
# determinant and trace
Wdet = Wxx*Wyy - Wxy**2
Wtr = Wxx + Wyy
return Wdet / Wtr
def get_harris_points(harrisim,min_dist=10,threshold=0.1):
""" Return corners from a Harris response image
min_dist is the minimum number of pixels separating
corners and image boundary. """
# find top corner candidates above a threshold
corner_threshold = harrisim.max() * threshold
harrisim_t = (harrisim > corner_threshold) * 1
# get coordinates of candidates
coords = array(harrisim_t.nonzero()).T
# ...and their values
candidate_values = [harrisim[c[0],c[1]] for c in coords]
# sort candidates (reverse to get descending order)
index = argsort(candidate_values)[::-1]
# store allowed point locations in array
allowed_locations = zeros(harrisim.shape)
allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1
# select the best points taking min_distance into account
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
def plot_harris_points(image,filtered_coords):
""" Plots corners found in image. """
figure()
gray()
imshow(image)
plot([p[1] for p in filtered_coords],[p[0] for p in filtered_coords],'*')
axis('off')
show()
def get_descriptors(image,filtered_coords,wid=5):
""" For each point return pixel values around the point
using a neighbourhood of width 2*wid+1. (Assume points are
extracted with min_distance > wid). """
desc = []
for coords in filtered_coords:
patch = image[coords[0]-wid:coords[0]+wid+1,
coords[1]-wid:coords[1]+wid+1].flatten()
desc.append(patch)
return desc
def match(desc1,desc2,threshold=0.5):
""" For each corner point descriptor in the first image,
select its match to second image using
normalized cross correlation. """
n = len(desc1[0])
# pair-wise distances
d = -ones((len(desc1),len(desc2)))
for i in range(len(desc1)):
for j in range(len(desc2)):
d1 = (desc1[i] - mean(desc1[i])) / std(desc1[i])
d2 = (desc2[j] - mean(desc2[j])) / std(desc2[j])
ncc_value = sum(d1 * d2) / (n-1)
if ncc_value > threshold:
d[i,j] = ncc_value
ndx = argsort(-d)
matchscores = ndx[:,0]
return matchscores
def match_twosided(desc1,desc2,threshold=0.5):
""" Two-sided symmetric version of match(). """
matches_12 = match(desc1,desc2,threshold)
matches_21 = match(desc2,desc1,threshold)
ndx_12 = where(matches_12 >= 0)[0]
# remove matches that are not symmetric
for n in ndx_12:
if matches_21[matches_12[n]] != n:
matches_12[n] = -1
return matches_12
def appendimages(im1,im2):
""" Return a new image that appends the two images side-by-side. """
# select the image with the fewest rows and fill in enough empty rows
rows1 = im1.shape[0]
rows2 = im2.shape[0]
if rows1 < rows2:
im1 = concatenate((im1,zeros((rows2-rows1,im1.shape[1]))),axis=0)
elif rows1 > rows2:
im2 = concatenate((im2,zeros((rows1-rows2,im2.shape[1]))),axis=0)
# if none of these cases they are equal, no filling needed.
return concatenate((im1,im2), axis=1)
def plot_matches(im1,im2,locs1,locs2,matchscores,show_below=True):
""" Show a figure with lines joining the accepted matches
input: im1,im2 (images as arrays), locs1,locs2 (feature locations),
matchscores (as output from 'match()'),
show_below (if images should be shown below matches). """
im3 = appendimages(im1,im2)
if show_below:
im3 = vstack((im3,im3))
imshow(im3)
cols1 = im1.shape[1]
for i,m in enumerate(matchscores):
if m>0:
plot([locs1[i][1],locs2[m][1]+cols1],[locs1[i][0],locs2[m][0]],'c')
axis('off')
实验原图与结果如下：
2.Harris算法在图像间寻找对应点
利用Harris算法比较一下两张图像之间的特征点
代码实现部分：
from scipy.ndimage import filters
from numpy import *
from PIL import Image
import harris
from pylab import *
import os
def imresize(im,sz):
""" Resize an image array using PIL. """
pil_im = Image.fromarray(uint8(im))
return array(pil_im.resize(sz))
im1 = array(Image.open('2.jpg').convert('L'))
im2 = array(Image.open('3.jpg').convert('L'))
im1 = imresize(im1, (int(im1.shape[1]/2), int(im1.shape[0]/2)))
im2 = imresize(im2, (int(im2.shape[1]/2), int(im2.shape[0]/2)))
wid=5
harrisim=harris.compute_harris_response(im1,5)
filtered_coords1=harris.get_harris_points(harrisim,wid+1)
d1=harris.get_descriptors(im1,filtered_coords1,wid)
harrisim=harris.compute_harris_response(im2,5)
filtered_coords2=harris.get_harris_points(harrisim,wid+1)
d2=harris.get_descriptors(im2,filtered_coords2,wid)
print ('starting matching')
matches=harris.match_twosided(d1,d2)
figure()
gray()
harris.plot_matches(im1,im2,filtered_coords1,filtered_coords2,matches)
show()实验结果：
3.SIFT算法检测与Harris角点检测对比
这里要强调一个很重要、重要、重要的问题，就是VLfeat的配置以及PCV应用VLfeat工具包的问题。
由于我第一次使用VLfeat工具包与PCV包，产生了非常麻烦的问题，差不多从困扰了我几个小时，在这里特别提出，希望读者在进行VLfeat工具包以及PCV包的配置时不会出现这个让人崩溃的问题。在PCV的sift函数里有一个函数，它导入了VLfeat包，在路径导入的时候，需要在前面加一个r，不然会出现各种各样莫名其妙的报错。
下面这个是sift函数：
from PIL import Image
import os
from numpy import *
from pylab import *
from scipy.ndimage import filters
def imresize(im,sz):
pil_im = Image.fromarray(uint8(im))
return array(pil_im.resize(sz))
def compute_harris_response(im,sigma=3):
imx = zeros(im.shape)
filters.gaussian_filter(im, (sigma,sigma), (0,1), imx)
imy = zeros(im.shape)
filters.gaussian_filter(im, (sigma,sigma), (1,0), imy)
Wxx = filters.gaussian_filter(imx*imx,sigma)
Wxy = filters.gaussian_filter(imx*imy,sigma)
Wyy = filters.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.1):
corner_threshold = harrisim.max() * threshold
harrisim_t = (harrisim > corner_threshold) * 1
coords = array(harrisim_t.nonzero()).T
candidate_values = [harrisim[c[0],c[1]] for c in coords]
index = argsort(candidate_values)[::-1]
allowed_locations = zeros(harrisim.shape)
allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1
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
def plot_harris_points(image,filtered_coords):
""" Plots corners found in image. """
figure()
gray()
imshow(image)
plot([p[1] for p in filtered_coords],[p[0] for p in filtered_coords],'*')
axis('off')
show()
def get_descriptors(image,filtered_coords,wid=5):
""" For each point return pixel values around the point
using a neighbourhood of width 2*wid+1. (Assume points are
extracted with min_distance > wid). """
desc = []
for coords in filtered_coords:
patch = image[coords[0]-wid:coords[0]+wid+1,
coords[1]-wid:coords[1]+wid+1].flatten()
desc.append(patch)
return desc
def match(desc1,desc2,threshold=0.5):
""" For each corner point descriptor in the first image,
select its match to second image using
normalized cross correlation. """
n = len(desc1[0])
# pair-wise distances
d = -ones((len(desc1),len(desc2)))
for i in range(len(desc1)):
for j in range(len(desc2)):
d1 = (desc1[i] - mean(desc1[i])) / std(desc1[i])
d2 = (desc2[j] - mean(desc2[j])) / std(desc2[j])
ncc_value = sum(d1 * d2) / (n-1)
if ncc_value > threshold:
d[i,j] = ncc_value
ndx = argsort(-d)
matchscores = ndx[:,0]
return matchscores
def match_twosided(desc1,desc2,threshold=0.5):
""" Two-sided symmetric version of match(). """
matches_12 = match(desc1,desc2,threshold)
matches_21 = match(desc2,desc1,threshold)
ndx_12 = where(matches_12 >= 0)[0]
# remove matches that are not symmetric
for n in ndx_12:
if matches_21[matches_12[n]] != n:
matches_12[n] = -1
return matches_12
def appendimages(im1,im2):
""" Return a new image that appends the two images side-by-side. """
# select the image with the fewest rows and fill in enough empty rows
rows1 = im1.shape[0]
rows2 = im2.shape[0]
if rows1 < rows2:
im1 = concatenate((im1,zeros((rows2-rows1,im1.shape[1]))),axis=0)
elif rows1 > rows2:
im2 = concatenate((im2,zeros((rows1-rows2,im2.shape[1]))),axis=0)
# if none of these cases they are equal, no filling needed.
return concatenate((im1,im2), axis=1)
def plot_matches(im1,im2,locs1,locs2,matchscores,show_below=True):
""" Show a figure with lines joining the accepted matches
input: im1,im2 (images as arrays), locs1,locs2 (feature locations),
matchscores (as output from 'match()'),
show_below (if images should be shown below matches). """
im3 = appendimages(im1,im2)
if show_below:
im3 = vstack((im3,im3))
imshow(im3)
cols1 = im1.shape[1]
for i,m in enumerate(matchscores):
if m>0:
plot([locs1[i][1],locs2[m][1]+cols1],[locs1[i][0],locs2[m][0]],'c')
axis('off')
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:\VL20\vlfeat-0.9.20\bin\win64\sift.exe "+imagename+" --output="+resultname+" "+params)
os.system(cmmd)
def read_features_from_file(filename):
"""读取特征属性值，然后将其以矩阵的形式返回"""
f = loadtxt(filename)
return f[:,:4], f[:,4:] #特征位置，描述子
def write_featrues_to_file(filename, locs, desc):
"""将特征位置和描述子保存到文件中"""
savetxt(filename, hstack((locs,desc)))
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 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:\VL20\vlfeat-0.9.20\bin\win64\sift.exe "+imagename+" --output="+resultname+" "+params)
os.system(cmmd)下面是实现SIFT匹配和Harris角点匹配算法的对比：
from PIL import Image
from pylab import *
from PCV.localdescriptors import sift
from PCV.localdescriptors import harris
# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)
imname = '2.jpg'
im = array(Image.open(imname).convert('L'))
sift.process_image(imname, 'empire.sift')
l1, d1 = sift.read_features_from_file('empire.sift')
figure()
gray()
subplot(131)
sift.plot_features(im, l1, circle=False)
title(u'SIFT特征',fontproperties=font)
subplot(132)
sift.plot_features(im, l1, circle=True)
title(u'用圆圈表示SIFT特征尺度',fontproperties=font)
# 检测harris角点
harrisim = harris.compute_harris_response(im)
subplot(133)
filtered_coords = harris.get_harris_points(harrisim, 6, 0.1)
imshow(im)
plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '*')
axis('off')
title(u'Harris角点',fontproperties=font)
show()
以下是实验结果：
可以明显看到SIFT所确定的特征逼Harris角点要优秀，并且在图片大小较大的时候，SIFT的鲁棒性要优于Harris角点：
可以看到在这个时候，SIFT依然很正常得计算出了结果，然而Harris角点已然失效。
4.SIFT算法匹配
以下是实现算法
from PIL import Image
from pylab import *
import sys
from PCV.localdescriptors import sift
if len(sys.argv) >= 3:
im1f, im2f = sys.argv[1], sys.argv[2]
else:
# im1f = '../data/sf_view1.jpg'
# im2f = '../data/sf_view2.jpg'
im1f = '1.jpg'
im2f = '2.jpg'
# im1f = '../data/climbing_1_small.jpg'
# im2f = '../data/climbing_2_small.jpg'
im1 = array(Image.open(im1f))
im2 = array(Image.open(im2f))
sift.process_image(im1f, 'out_sift_1.txt')
l1, d1 = sift.read_features_from_file('out_sift_1.txt')
figure()
gray()
subplot(121)
sift.plot_features(im1, l1, circle=False)
sift.process_image(im2f, 'out_sift_2.txt')
l2, d2 = sift.read_features_from_file('out_sift_2.txt')
subplot(122)
sift.plot_features(im2, l2, circle=False)
#matches = sift.match(d1, d2)
matches = sift.match_twosided(d1, d2)
print ('{} matches').format(len(matches.nonzero()[0]))
figure()
gray()
sift.plot_matches(im1, im2, l1, l2, matches, show_below=True)
show()
以下是上诉代码的匹配结果：
5.pylot的应用
(graph_type='graph')
g.add_node(pydot.Node(str(0), fontcolor='transparent'))
for i in range(5):
g.add_node(pydot.Node(str(i + 1)))
g.add_edge(pydot.Edge(str(0), str(i + 1)))
for j in range(5):
g.add_node(pydot.Node(str(j + 1) + '0' + str(i + 1)))
g.add_edge(pydot.Edge(str(j + 1) + '0' + str(i + 1), str(j + 1)))
g.write_png('C:/Users/ASUS/Desktop/Ss/computerview/FIRST Demo/image/ a.png ', prog='neato')
要注意这里保存文件的斜杠要用好，方向错了会十分麻烦。
未完待续。。。（代码跑得有点慢）
扫一扫
57万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
