Python計算機視覺-第2章

本章旨在尋找影像間的對應點和對應區域。本章將介紹用於影像匹配的兩種區域性描 述子演算法。本書的很多內容中都會用到這些區域性特徵，它們在很多應用中都有重要 作用，比如建立全景圖、擴增實境技術以及計算影像的三維重建。

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計算機視覺程式設計，去下載相應的部分圖片了。

（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()

白宮圖片匹配執行結果：

vlfeat最新版vlfeat-0.9.21有問題（老是報錯記憶體不夠分配），回退到vlfeat-0.9.20即可，特此說明。