用Python处理图像自学（三）
这一篇用Harris和SIFT两种方法检测角点并实现图像间的匹配
所谓角点，举个例子，桌子拐角
Harris
检测点
# --------------------------Harris检测角点-------------------------------
from PIL import Image
import numpy as np
from pylab import *
from scipy.ndimage import filters
import matplotlib.pylab as plt
import matplotlib as mpl#添加中文字体
mpl.rcParams['font.family']='SimHei'
mpl.rcParams['font.sans-serif']=['SimHei']
mpl.rcParams['axes.unicode_minus']=False # 正常显示负号
#-------------------------对每一个像素计算Harris角点检测器-------------------------------------------------------------
def compute_harris_response(im,sigma=3):
#计算导数
imx=np.zeros(im.shape)
filters.gaussian_filter(im,(sigma,sigma),(0,1),imx)
imy = np.zeros(im.shape)
filters.gaussian_filter(im, (sigma, sigma), (1, 0), imy)
#计算Harris矩阵的分量
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):#min_dist值小则稀疏，threshold值大则点少
conner_threshold = harrisim.max() * threshold
harrisim_t = (harrisim > conner_threshold) * 1
coords = np.array(harrisim_t.nonzero()).T
candidate_values = [harrisim[c[0], c[1]] for c in coords]
index = np.argsort(candidate_values)
allowed_locations = np.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 # 此处保证min_dist*min_dist仅仅有一个harris特征点
return filtered_coords
im = np.array(Image.open('test.jpg').convert('L'))
harrisim = compute_harris_response(im)
plt.figure(num='Harris')
plt.gray()
plt.subplot(141)
plt.imshow(255-harrisim)
plt.title('响应函数')
plt.axis('off')
threshold=[0.01,0.1,0.5]#这里只是为了有个对比，不想的话在前面的threshold就可以写等于几了
for i,items in enumerate(threshold):
filtered_coords = get_harris_points(harrisim,10,items)
plt.subplot(1,4,i+2)
plt.imshow(im)
plt.plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '+')
plt.title(str(threshold[i]))
plt.axis('off')
plt.show()
效果如下：
图像间匹配
在匹配图这一点上，选图很重要，一个是图片太大容易要很久才出结果，二是我之前出现过im1=图1和im2=图1就报错，交换又正常的情况的情况。当然，图位置不一肯定会有不同，具体见效果图
#---------------------------图像间匹配------------------------------------------------------
def get_descriptors(image, filtered_coords, wid=5):#添加到描述子列表中
"""对于每个返回的点，返回点周围2*wid+1个像素的值（假设选取点的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):#描述子匹配另一图的最优候选点
"""对于第一幅图像中的每个角点描述子，使用归一化互相关，选取它在第二幅图像中的匹配角点"""
n = len(desc1[0])
# 点对的距离
d = -np.ones((len(desc1), len(desc2)))
for i in range(len(desc1)):
for j in range(len(desc2)):
d1 = (desc1[i] - np.mean(desc1[i]))/np.std(desc1[i])
d2 = (desc2[j] - np.mean(desc2[j]))/np.std(desc2[j])
# np.seterr(divide='ignore', invalid='ignore')
ncc_value = sum(d1 * d2)/(n - 1)
if ncc_value > threshold:
d[i, j] = ncc_value
ndx = np.argsort(-d)
matchscores = ndx[:, 0]
return matchscores
def match_twosided(desc1, desc2, threshold=0.5):#图一配图二或者反过来，取最好的
"""两边对称版本的match（）"""
matches_12 = match(desc1, desc2, threshold)
matches_21 = match(desc2, desc1, threshold)
ndx_12 = np.where(matches_12 >= 0)[0]
# 去除非对称的匹配
for n in ndx_12:
if matches_21[matches_12[n]] != n:
matches_12[n] = -1
return matches_12
def appendimages(im1, im2):
"""返回将两幅图像并排拼接成的一幅新图像"""
# 选取具有最少行数的图像，然后填充足够的空行
row1 = im1.shape[0]
row2 = im2.shape[0]
if row1 < row2:
im1 = np.concatenate((im1, np.zeros((row2 - row1, im1.shape[1]))), axis=0)
elif row1 > row2:
im2 = np.concatenate((im2, np.zeros((row1 - row2, im2.shape[1]))), axis=0)
# 如果这些情况都没有，那么他们的行数相同，不需要进行填充
return np.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 = np.vstack((im3, im3))
plt.imshow(im3)
cols1 = im1.shape[1]
for i, m in enumerate(matchscores):
if m > 0:
plt.plot([locs1[i][1], locs2[m][1] + cols1], [locs1[i][0], locs2[m][0]], 'c')
plt.axis('off')
plt.show()
def imresize(im, sz):
pil_im = Image.fromarray(uint8(im))
return np.array(pil_im.resize(sz))
im1 = np.array(Image.open(r'2.jpg').convert('L'))
im2 = np.array(Image.open(r'1.jpg').convert('L'))
# resize加快匹配速度
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
harrisim1 = compute_harris_response(im1,5)
filtered_coords1 = get_harris_points(harrisim1, wid+1)
d1 = get_descriptors(im1, filtered_coords1, wid)
harrisim2 = compute_harris_response(im2,5)
filtered_coords2 = get_harris_points(harrisim2, wid + 1)
d2 = get_descriptors(im1, filtered_coords2, wid)
print('starting matching')
matches = match_twosided(d1,d2)
plt.figure(num='Harris图像间匹配')
# plt.gray()
plot_matches(im1, im2, filtered_coords1, filtered_coords2, matches)
plt.show()
效果如下：
同样的图片，交换下位置，匹配就会有不同喽
SIFT
检测点
注意：cmmd那一行的地址很重要，如果你直接复制代码，但电脑里没有安装或者路径不对，一定会报错
所以请先做好以下准备：
1.下载开源工具包VLFeat，下载地址，下载vlfeat-0.9.20-bin.tar.gz，感觉0.9.21有点问题，也有可能是我不会搞。
2.把vlfeat-0.9.20\bin\win64文件夹下的sift.exe、vl.dll和vl.lib这三个文件复制到项目的文件夹中。
3.在下图的位置写上sift.exe的文件所在。
import os
import matplotlib.pylab as plt
import numpy as np
from PIL import Image
import sys
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"！！！很重要！！！！！！ "+imagename+" --output="+resultname+" "+params)
os.system(cmmd)#r后面必须下载安装才有效果
print('processed',imagename,'to',resultname)
def read_features_from_file(filename):
"""读取特征属性值，然后将其以矩阵的形式返回"""
f =np.loadtxt(filename)
return f[:,:4], f[:,4:] #特征位置，描述子
def write_featrues_to_file(filename, locs, desc):
"""将特征位置和描述子保存到文件中"""
np.savetxt(filename, np.hstack((locs,desc)))
def plot_features(im, locs, circle=False):
"""显示带有特征的图像
输入：im（数组图像），locs（每个特征的行、列、尺度和朝向）"""
def draw_circle(c,r):
t = np.arange(0,1.01,.01)*2*np.pi
x = r*np.cos(t) + c[0]
y = r*np.sin(t) + c[1]
plt.plot(x, y, 'b', linewidth=2)
plt.imshow(im)
if circle:
for p in locs:
draw_circle(p[:2], p[2])
else:
plt.plot(locs[:,0], locs[:,1], 'ob')
plt.axis('off')
imname ='1.jpg'
im1 = np.array(Image.open(imname).convert('L'))
process_image(imname, 'empire.sift')
l1,d1 = read_features_from_file('empire.sift')
plt.figure(num='sift检测点')
plt.gray()
plot_features(im1, l1, circle=True)
plt.show()
效果如下：
图像间匹配
这个匹配真的是搞了我很久，总是报错，最后终于弄出来了
#--------------------------------------sift图像间----------------------------------
def match(desc1, desc2):
"""对于第一幅图像的每个描述子，选取其在第二幅图像中的匹配
输入：desc1（第一幅图像中的描述子），desc2（第二幅图像中的描述子）"""
desc1 = np.array([d/np.linalg.norm(d) for d in desc1])
desc2 = np.array([d/np.linalg.norm(d) for d in desc2])
dist_ratio = 0.6
desc1_size = desc1.shape
matchscores = np.zeros((desc1_size[0],1), 'int')
desc2t = desc2.T #预先计算矩阵转置
for i in range(desc1_size[0]):
dotprods = np.dot(desc1[i,:], desc2t) #向量点乘
dotprods = 0.9999*dotprods
# 反余弦和反排序，返回第二幅图像中特征的索引
index = np.argsort(np.arccos(dotprods))
# 检查最近邻的角度是否小于dist_ratio乘以第二近邻的角度
if np.arccos(dotprods)[index[0]] < dist_ratio * np.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 = np.concatenate((im1,np.zeros((row2-row1,im1.shape[1]))), axis=0)
elif row1 > row2:
im2 = np.concatenate((im2,np.zeros((row1-row2,im2.shape[1]))), axis=0)
# 如果这些情况都没有，那么他们的行数相同，不需要进行填充
return np.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 = np.vstack((im3,im3))
plt.imshow(im3)
cols1 = im1.shape[1]
for i in range(len(matchscores)):
if matchscores[i] > 0:
plt.plot([locs1[i, 0], locs2[matchscores[i, 0], 0] + cols1], [locs1[i, 1], locs2[matchscores[i, 0], 1]], 'c')
plt.axis('off')
plt.show()#你可以试试注释这一行，你会有新的收获
im1f = r'1.jpg'
im2f = r'2.jpg'
im1 = np.array(Image.open(im1f))
im2 = np.array(Image.open(im2f))
plt.figure(num='SIFT匹配')
plt.gray()
process_image(imname, 'out_sift_1.txt')
l1, d1 = read_features_from_file('out_sift_1.txt')
process_image(imname, 'out_sift_2.txt')
l2, d2 = read_features_from_file('out_sift_2.txt')
matches = match_twosided(d1, d2)
print('{} matches'.format(len(matches.nonzero()[0])))
plot_matches(im1, im2, l1, l2, matches, show_below=True)
# plt.show()
效果如下：
用了我母校的风景呢~~
对比了之后，很明显，SIFT完胜
当然也有人专门用了同一照片的不同角度，不同明暗，旋转等之后进行过对比的，但是我忘了地址了，不然放上个链接还挺有学习意义
扫一扫
7万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
