本章旨在寻找图像间的对应点和对应区域。本章将介绍用于图像匹配的两种局部描 述子算法。本书的很多内容中都会用到这些局部特征,它们在很多应用中都有重要 作用,比如创建全景图、增强现实技术以及计算图像的三维重建。
1、Harris角点检测器
Harris 角点检测算法(也称 Harris & Stephens 角点检测器)是一个极为简单的角点 检测算法。该算法的主要思想是,如果像素周围显示存在多于一个方向的边,我们 认为该点为兴趣点。该点就称为角点。我们把图像域中点 x 上的对称半正定矩阵定义为:
选择权重矩阵 W(通常为高斯滤波器 Gσ),我们可以得到卷积:
Harris 矩阵的特征值有三种情况:
-
如果 λ1 和 λ2 都是很大的正数,则该 x 点为角点;
-
如果 λ1 很大,λ2 ≈ 0,则该区域内存在一个边,该区域内的平均的特征值不
会变化太大;
-
如果 λ1≈ λ2 ≈ 0,该区域内为空。
在不需要实际计算特征值的情况下,为了把重要的情况和其他情况分开,Harris 和 Stephens 引入了指示函数:
为了去除加权常数 κ,我们通常使用商数:
作为指示器。
下面我们写出 Harris 角点检测程序。代码实现:
# -*- coding: utf-8 -*-
from pylab import *
from PIL import Image
# from PCV.localdescriptors import harris
import harris
"""
Example of detecting Harris corner points (Figure 2-1 in the book).
"""
# 读入图像
im = array(Image.open('../data/empire.jpg').convert('L'))
# 检测harris角点
harrisim = harris.compute_harris_response(im)
# Harris响应函数
harrisim1 = 255 - harrisim
figure()
gray()
#画出Harris响应图
subplot(141)
imshow(harrisim1)
print(harrisim1.shape)
axis('off')
axis('equal')
threshold = [0.01, 0.05, 0.1]
for i, thres in enumerate(threshold):
filtered_coords = harris.get_harris_points(harrisim, 6, thres)
subplot(1, 4, i+2)
imshow(im)
print(im.shape)
plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '*')
axis('off')
#原书采用的PCV中PCV harris模块
# filtered_coords = harris.get_harris_points(harrisim,6)
# harris.plot_harris_points(im, filtered_coords)
# plot only 200 strongest
# harris.plot_harris_points(im, filtered_coords[:200])
show()
其中compute_harris_response函数、 get_harris_points函数和 plot_harris_points函数代码实现详见harris.py文件。
运行结果:
在图像间寻找对应点,代码实现:
# -*- coding: utf-8 -*-
from pylab import *
from PIL import Image
# from PCV.localdescriptors import harris
# from PCV.tools.imtools import imresize
import harris
from imtools import imresize
"""
This is the Harris point matching example in Figure 2-2.
"""
# Figure 2-2上面的图
#im1 = array(Image.open("../data/crans_1_small.jpg").convert("L"))
#im2 = array(Image.open("../data/crans_2_small.jpg").convert("L"))
# Figure 2-2下面的图
im1 = array(Image.open("../data/sf_view1.jpg").convert("L"))
im2 = array(Image.open("../data/sf_view2.jpg").convert("L"))
# resize to make matching faster
im1 = imresize(im1, (im1.shape[1]//2, im1.shape[0]//2))
im2 = imresize(im2, (im2.shape[1]//2, 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()
运行结果:
2、SIFT(尺度不变特征变换)
David Lowe 在文献 [17] 中提出的 SIFT(Scale-Invariant Feature Transform,尺度不 变特征变换)是过去十年中最成功的图像局部描述子之一。SIFT 特征后来在文献 [18] 中得到精炼并详述,经受住了时间的考验。SIFT 特征包括兴趣点检测器和描述 子。SIFT 描述子具有非常强的稳健性,这在很大程度上也是 SIFT 特征能够成功和 流行的主要原因。自从 SIFT 特征的出现,许多其他本质上使用相同描述子的方法 也相继出现。现在,SIFT 描述符经常和许多不同的兴趣点检测器相结合使用(有些 情况下是区域检测器),有时甚至在整幅图像上密集地使用。SIFT 特征对于尺度、 旋转和亮度都具有不变性,因此,它可以用于三维视角和噪声的可靠匹配。你可以在 http://en.wikipedia.org/wiki/Scale-invariant_feature_transform 获得 SIFT 特征的简 要介绍。
(1)兴趣点
SIFT 特征使用高斯差分函数来定位兴趣点:
其中,Gσ 是上一章中介绍的二维高斯核,Iσ 是使用 Gσ 模糊的灰度图像,κ 是决定相 差尺度的常数。
(2)描述子
为了实现旋转不变性,基于每个点周围图像梯度的方向和大小,SIFT 描述子又引入了参考方向。SIFT 描述子使用主方向描述参考方向。主方向使用方向直方图(以大小为权 重)来度量。
构造 SIFT 描述子特征向量的图解:(a)一个围绕兴趣点的网格结构,其中该网格已 经按照梯度主方向进行了旋转;(b)在网格的一个子区域内构造梯度方向的 8-bin 直方图;(c) 在网格的每个子区域内提取直方图;(d)拼接直方图,得到一个长的特征向量。
(3)检测兴趣点
代码实现:
# -*- coding: utf-8 -*-
from PIL import Image
from pylab import *
# from PCV.localdescriptors import sift
# from PCV.localdescriptors import harris
import sift
import harris
# 添加中文字体支持
# from matplotlib.font_manager import FontProperties
# font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)
from matplotlib.font_manager import FontProperties
font = FontProperties(fname="/System/Library/Fonts/PingFang.ttc", size=14)
imname = '../data/empire.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()
运行结果:
(4)匹配描述子
代码实现:
from PIL import Image
from pylab import *
import sys
# from PCV.localdescriptors import sift
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 = '../data/crans_1_small.jpg'
# im2f = '../data/crans_2_small.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()
运行结果:
3、匹配地理标记图像
(1)从Panoramio下载地理标记图像
2016年Google已经关闭Panoramio(http://www.panoramio.com/)网站,所以只能从Programming Computer Vision with Python书中代码对应中文主页Python计算机视觉编程(http://yongyuan.name/pcvwithpython/),去下载相应的部分图片了。
(2)使用局部描述子匹配
代码实现:
import sift
import imtools
from pylab import *
#获取图像列表
imlist = imtools.get_imlist('../data/panoimages/')
nbr_images = len(imlist)
print("imlist\n")
print(imlist)
nbr_images = len(imlist)
print(nbr_images)
# extract features
featlist = [imname[:-3] + 'sift' for imname in imlist]
print("featlist\n")
print(featlist)
matchscores = zeros((nbr_images,nbr_images))
for i in range(nbr_images):
for j in range(i,nbr_images): # 仅仅计算上三角
print('comparing ', imlist[i], imlist[j])
# process and save features to file
sift.process_image(imlist[i],featlist[i])
l1,d1 = sift.read_features_from_file(featlist[i])
sift.process_image(imlist[j],featlist[j])
l2,d2 = sift.read_features_from_file(featlist[j])
matches = sift.match_twosided(d1,d2)
nbr_matches = sum(matches > 0)
print('number of matches = ', nbr_matches)
matchscores[i,j] = nbr_matches
# 复制值
for i in range(nbr_images):
for j in range(i+1,nbr_images): # 不需要复制对角线
matchscores[j,i] = matchscores[i,j]
print("matchscores:\n",matchscores)
运行结果:
matchscores:
[[7.060e+02 3.800e+01 0.000e+00 0.000e+00 1.300e+01 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 5.000e+00
0.000e+00 0.000e+00 1.000e+00 0.000e+00 0.000e+00 0.000e+00]
[3.800e+01 5.950e+02 0.000e+00 0.000e+00 2.000e+00 0.000e+00 1.000e+00
0.000e+00 1.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 3.000e+00
1.000e+00 0.000e+00 2.000e+00 0.000e+00 0.000e+00 0.000e+00]
[0.000e+00 0.000e+00 8.510e+02 0.000e+00 0.000e+00 0.000e+00 1.000e+00
4.000e+00 0.000e+00 0.000e+00 1.000e+00 1.000e+00 0.000e+00 0.000e+00
1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.600e+01]
[0.000e+00 0.000e+00 0.000e+00 1.685e+03 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00]
[1.300e+01 2.000e+00 0.000e+00 0.000e+00 4.160e+02 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00]
[0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 1.032e+03 0.000e+00
0.000e+00 2.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 1.000e+00 0.000e+00 0.000e+00 0.000e+00]
[0.000e+00 1.000e+00 1.000e+00 0.000e+00 0.000e+00 0.000e+00 8.100e+02
1.000e+00 0.000e+00 0.000e+00 1.000e+00 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00]
[0.000e+00 0.000e+00 4.000e+00 0.000e+00 0.000e+00 0.000e+00 1.000e+00
5.490e+02 0.000e+00 0.000e+00 1.000e+00 0.000e+00 0.000e+00 0.000e+00
3.000e+00 9.000e+00 0.000e+00 0.000e+00 1.000e+00 6.000e+00]
[0.000e+00 1.000e+00 0.000e+00 0.000e+00 0.000e+00 2.000e+00 0.000e+00
0.000e+00 5.220e+02 0.000e+00 2.000e+00 1.000e+00 0.000e+00 0.000e+00
0.000e+00 2.000e+00 0.000e+00 1.000e+00 0.000e+00 0.000e+00]
[0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 1.849e+03 0.000e+00 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00]
[0.000e+00 0.000e+00 1.000e+00 0.000e+00 0.000e+00 0.000e+00 1.000e+00
1.000e+00 2.000e+00 0.000e+00 1.154e+03 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 1.000e+00 0.000e+00 0.000e+00 3.000e+00]
[0.000e+00 0.000e+00 1.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
0.000e+00 1.000e+00 0.000e+00 0.000e+00 1.190e+02 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00]
[0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.716e+03 0.000e+00
0.000e+00 0.000e+00 1.300e+01 0.000e+00 0.000e+00 0.000e+00]
[5.000e+00 3.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 6.680e+02
0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00]
[0.000e+00 1.000e+00 1.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
3.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
1.464e+03 0.000e+00 0.000e+00 1.000e+00 0.000e+00 2.000e+00]
[0.000e+00 0.000e+00 1.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
9.000e+00 2.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
0.000e+00 5.570e+02 2.000e+00 7.000e+00 0.000e+00 1.700e+01]
[1.000e+00 2.000e+00 1.000e+00 0.000e+00 0.000e+00 1.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 1.000e+00 0.000e+00 1.300e+01 0.000e+00
0.000e+00 2.000e+00 1.391e+03 0.000e+00 1.000e+00 1.000e+00]
[0.000e+00 0.000e+00 1.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
0.000e+00 1.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
1.000e+00 7.000e+00 0.000e+00 5.400e+02 0.000e+00 3.000e+00]
[0.000e+00 0.000e+00 1.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
1.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 1.000e+00 0.000e+00 1.664e+03 1.000e+00]
[0.000e+00 0.000e+00 1.600e+01 0.000e+00 0.000e+00 0.000e+00 0.000e+00
6.000e+00 0.000e+00 0.000e+00 3.000e+00 0.000e+00 0.000e+00 0.000e+00
2.000e+00 1.700e+01 1.000e+00 3.000e+00 1.000e+00 8.800e+02]]
[Finished in 285.9s]
(3)可视化连接的图像
pydot例子代码实现:
import pydot
g = pydot.Dot(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) + '_' + str(i + 1)))
g.add_edge(pydot.Edge(str(j + 1) + '_' + str(i + 1), str(j + 1)))
# g.write_png('../images/ch02/ch02_fig2-9_graph.png', prog='neato')
# 源码中write_png函数已经改成write函数,默认生产图片format=raw所要重新设置图片格式format,不然普通图片工具打不开
g.write('graph.jpg', prog='neato', format='png', encoding=None)
from PIL import Image
from pylab import *
# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname="/System/Library/Fonts/PingFang.ttc", size=14)
figure()
im = array(Image.open('graph.jpg'))
imshow(im)
title(u'Pydot Sample', fontproperties=font)
axis('off')
show()
pydot例子运行结果:
白宫图片匹配代码实现:
# -*- coding: utf-8 -*-
from pylab import *
from PIL import Image
# from PCV.localdescriptors import sift
# from PCV.tools import imtools
import sift
import imtools
import pydot
import os
""" This is the example graph illustration of matching images from Figure 2-10.
To download the images, see ch2_download_panoramio.py."""
parent_path = os.path.abspath(os.path.join(os.path.dirname("__file__"),os.path.pardir))
download_path = parent_path+'/data/panoimages/' # set this to the path where you downloaded the panoramio images
path = parent_path+'/ch02/' # path to save thumbnails (pydot needs the full system path)
# list of downloaded filenames
imlist = imtools.get_imlist(download_path)
nbr_images = len(imlist)
# extract features
featlist = [imname[:-3] + 'sift' for imname in imlist]
matchscores = zeros((nbr_images,nbr_images))
for i in range(nbr_images):
for j in range(i, nbr_images): # only compute upper triangle
# print('comparing ', imlist[i],imlist[j])
# process and save features to file
sift.process_image(imlist[i],featlist[i])
l1, d1 = sift.read_features_from_file(featlist[i])
sift.process_image(imlist[j],featlist[j])
l2, d2 = sift.read_features_from_file(featlist[j])
matches = sift.match_twosided(d1, d2)
nbr_matches = sum(matches > 0)
# print('number of matches = ', nbr_matches)
matchscores[i, j] = nbr_matches
# copy values
for i in range(nbr_images):
for j in range(i + 1, nbr_images): # no need to copy diagonal
matchscores[j, i] = matchscores[i, j]
print("The match scores is: %d",matchscores)
#savetxt(("../data/panoimages/panoramio_matches.txt",matchscores)
threshold = 2 # min number of matches needed to create link
g = pydot.Dot(graph_type='graph') # don't want the default directed graph
for i in range(nbr_images):
for j in range(i + 1, nbr_images):
if matchscores[i, j] > threshold:
# first image in pair
im = Image.open(imlist[i])
im.thumbnail((100, 100))
filename = str(i) + '.png'
im.save(filename) # need temporary files of the right size
g.add_node(pydot.Node(str(i), fontcolor='transparent', \
shape='rectangle', image=path+filename))
# second image in pair
im = Image.open(imlist[j])
im.thumbnail((100, 100))
filename = str(j) + '.png'
im.save(filename) # need temporary files of the right size
g.add_node(pydot.Node(str(j), fontcolor='transparent', \
shape='rectangle', image=path+filename))
g.add_edge(pydot.Edge(str(i), str(j)))
# g.write_png('whitehouse.png')
g.write('whitehouse.png',format='png')
from PIL import Image
from pylab import *
# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname="/System/Library/Fonts/PingFang.ttc", size=14)
figure()
im = array(Image.open('whitehouse.png'))
imshow(im)
title(u'MatchGraph', fontproperties=font)
axis('off')
show()
白宫图片匹配运行结果:
4、其他
(1)开发环境
本代码是在mac电脑sublimetxt编辑器python3.7.8下调试出来的,如果你要在windows/linux下编辑器/IDE的其他python版本运行的话,请对代码做相应的调整。
(2)Panoramio和vlfeat
1)像前文中提到的,Panoramio(http://www.panoramio.com/)网站已经关闭,所以只能从Programming Computer Vision with Python书中代码对应中文主页http://yongyuan.name/pcvwithpython/,去下载相应的部分图片了(具体见(3)源码和图片)。下载的图片中,png格式是代码生成的,jpg格式才是原图片。
2)vlfeat最新版vlfeat-0.9.21有问题(老是报错内存不够分配),回退到vlfeat-0.9.20即可,特此说明。
(3)源码和图片
已经调试过的源码(含harris.py和sift.py文件)和图片详见:
https://github.com/Abonaventure/pcv-book-code.git 或 https://gitlab.com/Abonaventure/pcv-book-code.git
部分图片可能未上传可在https://pan.baidu.com/share/link?shareid=4059473268&uk=3440544532 或者原书目录http://programmingcomputervision.com下载