python计算机视觉编程第二章——局部图像描述子

一、Harris角点检测器 

Harris角点检测算法（也称Harris & Stephens角点检测器）是一个极为简单的角点检测算法。该算法的主要思想是,如果像素周围显示存在多于一个方向的边，我们认为该点为兴趣点。该点就称为角点。 

from PIL import Image
from matplotlib.pyplot import *
from scipy.ndimage import filters
from numpy import *
from matplotlib import pyplot as plt

def compute_harris_response(im,sigma=3):
 """ 在一幅灰度图像中，对每个像素计算 Harris 角点检测器响应函数 """
 # 计算导数
 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)
 # 计算 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=0.1):
    """ 从一幅 Harris 响应图像中返回角点。min_dist 为分割角点和图像边界的最少像素数目 """
    # 寻找高于阈值的候选角点
    corner_threshold = harrisim.max() * threshold
    harrisim_t = (harrisim > corner_threshold) * 1

    # 得到候选点的坐标
    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_distance 原则，选择最佳 Harris 点
    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

im = array(Image.open('eye.jpg').convert('L'))
harrisim = compute_harris_response(im)
filtered_coords = get_harris_points(harrisim,6,0.1)
filtered_coords1 = get_harris_points(harrisim,6,0.05)

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.subplot(141), plt.imshow(im, plt.cm.gray), plt.title('原始图像'), plt.axis('off')
plt.subplot(142), plt.imshow(harrisim, plt.cm.gray), plt.title('Harris响应函数'), plt.axis('off')
plt.subplot(143), plt.imshow(im, plt.cm.gray), plt.plot([p[1] for p in filtered_coords],
                                                         [p[0] for p in filtered_coords], '.'), plt.title(
    '阈值为0.1'), plt.axis('off')
plt.subplot(144), plt.imshow(im, plt.cm.gray), plt.plot([p[1] for p in filtered_coords1],
                                                         [p[0] for p in filtered_coords1], '.'), plt.title(
    '阈值为0.05'), plt.axis('off')
plt.show()

在图像间寻找对应点

兴趣点描述子是分配给兴趣点的一个向量，描述该点附近的图像的表观信息。描述子越好，寻找到的对应点越好。我们用对应点或者点的对应来描述相同物体和场景点在不同图像上形成的像素点。
Harris 角点的描述子通常是由周围图像像素块的灰度值，以及用于比较的归一化互相关矩阵构成的。图像的像素块由以该像素点为中心的周围矩形部分图像构成。

主要代码：

def get_descriptors(image, filtered_coords, wid=5):
    desc = []
    height, width = image.shape
    for coords in filtered_coords:
        # 但是当特征点位于图像边缘时，提取窗口可能会超出图像边界。所以需要通过添加边界检查来解决
        y, x = coords

        # 计算窗口的边界
        y_min = max(0, y - wid)
        y_max = min(height, y + wid + 1)
        x_min = max(0, x - wid)
        x_max = min(width, x + wid + 1)

        patch = image[y_min:y_max, x_min:x_max].flatten()

        desc.append(patch)
    return desc


def match(desc1,desc2,threshold=0.5):
    """ 对于第一幅图像中的每个角点描述子，使用归一化互相关，选取它在第二幅图像中的匹配角点 """
    n = len(desc1[0])
    # 点对的距离
    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])#mean平均数   std标准差
            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)#argsort函数返回的是数组值从小到大的索引值
    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 = 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):
    """ 返回将两幅图像并排拼接成的一幅新图像 """

    # 选取具有最少行数的图像，然后填充足够的空行
    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)

    # 如果这些情况都没有，那么它们的行数相同，不需要进行填充
    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,m in enumerate(matchscores):
        if m>0:
            plot([locs1[i][1],locs2[m][1]+cols1],[locs1[i][0],locs2[m][0]],'c')
    axis('off')

 

二、SIFT (尺度不变特征变换)

David Lowe提出的SIFT (Scale-Invariant Feature Transform，尺度不变特征变换）是过去十年中最成功的图像局部描述子之一。SIFT特征包括兴趣点检测器和描述子。SIFT描述子具有非常强的稳健性，这在很大程度上也是SIFT特征能够成功和流行的主要原因。SIFT特征对于尺度、旋转和亮度都具有不变性，因此，它可以用于三维视角和噪声的可靠匹配。

sift.py：

from PIL import Image
from pylab import *
import os

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"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:] # 特征位置，描述子


def write_features_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 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
        # 反余弦和反排序，返回第二幅图像中特征的索引
        indx = argsort(arccos(dotprods))
        # 检查最近邻的角度是否小于 dist_ratio 乘以第二近邻的角度
        if arccos(dotprods)[indx[0]] < dist_ratio * arccos(dotprods)[indx[1]]:
            matchscores[i] = int(indx[0])

    return matchscores


def appendimages(im1,im2):
    """ 返回将两幅图像并排拼接成的一幅新图像 """

    # 选取具有最少行数的图像，然后填充足够的空行
    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)

    # 如果这些情况都没有，那么它们的行数相同，不需要进行填充
    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,m in enumerate(matchscores):
        if m>0:
            plot([locs1[i][1],locs2[m][1]+cols1],[locs1[i][0],locs2[m][0]],'c')
    axis('off')

运行代码：

import sift
from PIL import Image
from pylab import *


imname = 'eye.jpg'
im1 = array(Image.open(imname).convert('L'))
sift.process_image(imname,'eye.sift')
l1,d1 = sift.read_features_from_file('eye.sift')
figure()
gray()
sift.plot_features(im1,l1,circle=True)
show()

其中： 

添加开源工具包VLFeat，下载链接：http://www.vlfeat.org/download/，最好下载vlfeat-0.9.20-bin.tar.gz。
1）把vlfeat-0.9.20\bin\win64文件夹下的sift.exe、vl.dll和vl.lib这三个文件复制到项目的文件夹中。
2）打开sift.py文件，修改其中的cmmd路径

 但

原因不详。

三、匹配地理标记图像

import os
import urllib, urlparse
import simplejson as json
# 查询图像
url = 'http://www.panoramio.com/map/get_panoramas.php?order=popularity&\
    set=public&from=0&to=20&minx=-77.037564&miny=38.896662&\
    maxx=-77.035564&maxy=38.898662&size=medium'
c = urllib.urlopen(url)

# 从 JSON 中获得每个图像的 url
j = json.loads(c.read())
imurls = []
for im in j['photos']:
    imurls.append(im['photo_file_url'])

# 下载图像
for url in imurls:
    image = urllib.URLopener()
    image.retrieve(url, os.path.basename(urlparse.urlparse(url).path))
    print 'downloading:', url

使用局部描述子匹配 

对所有组合图像对进行逐个匹配，如下: 

import sift
nbr_images = len(imlist)
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]
        l1,d1 = sift.read_features_from_file(featlist[i])
        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]

我们首先通过图像间是否具有匹配的局部描述子来定义图像间的连接，然后可视化这些连接情况。为了完成可视化，我们可以在图中显示这些图像，图的边代表连接。我们将会使用 pydot 工具包（http://code.google.com/p/pydot/），该工具包是功能强大的 GraphViz 图形库的 Python 接口。 Pydot 使用 Pyparsing（http://pyparsing.wikispaces.com/）和 GraphViz（Graphviz）；不用担心，这些都非常容易安装，只需要几分钟就可以安装成功。

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('graph.jpg',prog='neato')

我们接下来继续探讨地理标记图像处理的例子。为了创建显示可能图像组的图，如果匹配的数目高于一个阈值，我们使用边来连接相应的图像节点。为了得到图中的图像，需要使用图像的全路径（在下面例子中，使用 path 变量表示）。为了使图像看起来漂亮，我们需要将每幅图像尺度化为缩略图形式，缩略图的最大边为 100 像素。下面是具体实现代码： 

import pydot
threshold = 2 # 创建关联需要的最小匹配数目
g = pydot.Dot(graph_type='graph') # 不使用默认的有向图

for i in range(nbr_images):
    for j in range(i+1,nbr_images):
        if matchscores[i,j] > threshold:
            # 图像对中的第一幅图像
            im = Image.open(imlist[i])
            im.thumbnail((100,100))
            filename = str(i)+'.png'
            im.save(filename) # 需要一定大小的临时文件
            g.add_node(pydot.Node(str(i),fontcolor='transparent',shape='rectangle',image=path+filename))
            # 图像对中的第二幅图像
            im = Image.open(imlist[j])
            im.thumbnail((100,100))
            filename = str(j)+'.png'
            im.save(filename) # 需要一定大小的临时文件
            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')