二、【python计算机视觉编程】局部图像描述子

这篇博客旨在寻找图像间的对应点和对应区域，将介绍两种用于图像匹配的局部描述子算法，它们在很多应用中都有重要作用，比如创建全景图、增强现实（AR）技术和计算图像的三维重建。

（一）Harris角点检测器

概念： Harris角点检测器是一种简单的角点检测算法。其主要思想是，如果像素周围显示存在多于一个方向的边，就可以认为该点为兴趣点，该点就称为角点。

角点可以分为下面几种情况：

1. 一阶导数（即灰度的梯度）的局部最大所对应的像素点；
2. 两条及两条以上边缘的交点；
3. 指物体边缘变化不连续的方向；
4. 一阶导数最大，二阶导数为0；
5. 指图像中梯度值和梯度方向的变化速率都很高的点。
   

Harris角点检测器：
Harris 算子是一种简单的点特征提取算子，这种算子受信号处理中自相关函数的启发，给出与自相关函数相联系的矩阵M。M阵的特征值是自相关函数的一阶曲率，如果两个曲率值都高，那么就认为该点是特征点。为了消除噪声对于角点检测的影响，可以使用一个高斯滤波器来平滑图像。

数学推导：
Harris算法是利用的窗口内图像灰度的自相关性进行的，设定一个窗口，并在图像中移动，计算移动前与移动后窗口所在区域图像的自相关系数。

自相关函数计算如下：(x,y)为窗口中心位置，w(u,v)为权重（一般取高斯函数），L表示窗口，(u,v)表示窗口中的图像位置。

泰勒展开并取近似值：

将近似值代入自相关函数，有：

将平方项展开并写成矩阵形式，有：

回到自相关表达式：

经过上面的数学形式推导，已经得到了自相关函数的表达式。可以看得这也是一个椭圆的矩阵表示形式（非标准椭圆），因此其系数矩阵M的特征值与椭圆的半轴长短有关。
假设M的特征值为λ1、λ2，有三种情况：

• 如果λ1和λ2都是很大的正数，则该x点为角点；
• 如果λ1很大，λ2≈0，则该区域内存在一个边，该区域内的平均M_I的特征值不会变化很大；
• 如果λ1≈λ2≈0，该区域内为空。

通过上面的情况，计算出特征值后就可以判别是否是角点了。
可以看出，这样计算量非常大，因为图像中的几乎每个点都需要进行一次特征值的计算，因此可以由下面的公式简易计算：

detM表示M的行列式，traceM表示M的迹，R表示角点响应值。α为加权常数，一般在0.04至0.06之间取值。

判断准则： 当R超过某个设定的阈值时，可认为是角点；反之，则不是。

数学定义：
把图像域中点x上的对称半正定矩阵M_I​=M_I​(x)定义为：

其中▽I为包含导数I_x和I_y的图像梯度。由于该定义，M_I的秩为1，特征值为λ₁=|▽I|²和λ₂=0。现在对于图像的每一个像素，都可以计算出该矩阵。
选择权重矩阵W（通常为高斯滤波器G_σ），我们可以得到卷积：

该卷积的目的是得到M_I在周围像素上的局部平均。计算出的矩阵称为Harris矩阵。

Harris角点性质：

1. 阈值决定检测点数量： 增大$\alpha$的值，将减小角点响应值$R$，降低角点检测的灵性，减少被检测角点的数量；减小$\alpha$值，将增大角点响应值$R$，增加角点检测的灵敏性，增加被检测角点的数量。
2. 对亮度和对比度的变化不敏感： 由于使用了微分算子对图像进行微分运算，而微分运算对图像密度的拉升或收缩和对亮度的抬高或下降不敏感。
3. 具有旋转不变性： 该算子使用的是角点附近的区域灰度二阶矩矩阵。而二阶矩矩阵可以表示成一个椭圆，椭圆的长短轴正是二阶矩矩阵特征值平方根的倒数。当特征椭圆转动时，特征值并不发生变化，所以判断角点响应值也不发生变化，由此说明Harris角点检测算子具有旋转不变性。
4. 不具有尺度不变性。

首先介绍几种在harris角点检测中需要用到的函数：

• 用于计算角点检测器响应函数：

from scipy.ndimage import filters
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),(0,1),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

• 显示图像中的角点：

def plot_harris_points(image,filtered_coords):
    """绘制图像中检测到的角点"""
    
    figure()
    gray()
    imshow(image)
    plot([p[1] for p in filtered_coords],[p[0] for p in filtered_coords],'*')
    axis('off')
    show()

编写实验代码：

from pylab import *
from PIL import Image
import harris
from numpy import *

#读取图像
im = array(Image.open('picture\house2.jpg').convert('L'))
#计算harris角点响应函数
harrisim = harris.compute_harris_response(im)

subplot(141)
imshow(im)
axis('off')

#获取角点响应图像并返回角点
#filtered_coords = harris.get_harris_points(harrisim,6)
#绘制图像中的角点
#harris.plot_harris_points(im,filtered_coords)

threshold = [0.01, 0.05, 0.1]
for i, thres in enumerate(threshold):
    filtered_coords = harris.get_harris_points(harrisim, 10,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')
    show()

第一次验证代码发现，得到的图像结果是没有检测角点的。如图所示：

根据这个情况，询问了同学也没有解决这个问题。注意一点是， 上面代码需要在python3才能使用，环境不一样会影响代码运行效果。

修改上面的代码后，得到如下结果：

from pylab import *
from PIL import Image
from PCV.localdescriptors import harris

# 读入图像
im = array(Image.open('picture\jimei1.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')

show()

结果分析：
发现使用不同的阈值，检测出的角点个数会不同。阈值越小，则检测出的角点会越多，符合该条件的角点就会多；增大阈值范围，则检测到的角点会变少，降低角点检测的灵性，减少被检测角点的数量。

想要进行改进harris角点检测器，可以在网站 http://en.wikipedia.org/wiki/Corner_detection 。但是我并不能打开这个网站，感兴趣的读者可以试试看。

图像中寻找对应点：
Harris 角点检测器仅仅能够检测出图像中的兴趣点，但是没有给出通过比较图像间的兴趣点来寻找匹配角点的方法。我们需要在每个点上加入描述子信息，并给出一个比较这些描述子的方法。

兴趣点描述子是分配给兴趣点的一个向量，描述该点附近的图像的表观信息。描述子越好，寻找到的对应点越好。我们用对应点或者点的对应来描述相同物体和场景点在不同图像上形成的像素点。

import harris
from pylab import *
from PIL import Image
 

im1 = array(Image.open('picture\house2.jpg').convert("L"))
im2 = array(Image.open('picture\house2.jpg').convert("L"))

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 mathching')
matches = harris.match_twosided(d1,d2)

figure()
gray()
harris.plot_matches(im1,im2,filtered_coords1,filtered_coords2,matches)
show()

结果分析： 发现该算法存在一些不正确的匹配，是由于图像像素块的互相关矩阵具有较弱的描述性，在实际运用中，可以使用更加稳健的方法来处理这些对应匹配。但是，总的来说，在要求并不是很高的情况下，这种描述点匹配的方法是最简便的一种处理方法。

（二）SIFT（尺度不变特征变换）

概念： 尺度不变特征转换(Scale-invariant feature transform或SIFT)是一种电脑视觉的算法用来侦测与描述影像中的局部性特征，它在空间尺度中寻找极值点，并提取出其位置、尺度、旋转不变量。该描述子具有非常强的稳健性。

应用： 物体辨识、机器人地图感知与导航、影像缝合、3D模型建立、手势辨识、影像追踪和动作比对。

特点：

1. SIFT特征是图像的局部特征，其对旋转、尺度缩放、亮度变化保持不变性，对视角变化、仿射变换、噪声也保持一定程度的稳定性；

2. 独特性好，信息量丰富，适用于在海量特征数据库中进行快速、准确的匹配；

3. 多量性，即使少数的几个物体也可以产生大量的SIFT特征向量；

4. 高速性，经优化的SIFT匹配算法甚至可以达到实时的要求；

5. 可扩展性，可以很方便的与其他形式的特征向量进行联合。

可以解决的问题：

目标的自身状态、场景所处的环境和成像器材的成像特性等因素影响图像配准/目标识别跟踪的性能，一般可以解决下面六种问题：

1. 目标的旋转、缩放、平移（RST）

2. 图像仿射/投影变换（视点viewpoint）

3. 光照影响（illumination）

4. 目标遮挡（occlusion）

5. 杂物场景（clutter）

6. 噪声

SIFT算法的实质是在不同的尺度空间上查找关键点(特征点)，并计算出关键点的方向。SIFT所查找到的关键点是一些十分稳定，不会因光照，仿射变换和噪音等因素而变化的点，如角点、边缘点、暗区的亮点及亮区的暗点等。

SIFT算法分解为以下四步：

1. 尺度空间极值检测：搜索所有尺度上的图像位置。通过高斯微分函数来识别潜在的对于尺度和旋转不变的兴趣点。

2. 关键点定位：在每个候选的位置上，通过一个拟合精细的模型来确定位置和尺度。关键点的选择依据于它们的稳定程度。

3. 方向确定：基于图像局部的梯度方向，分配给每个关键点位置一个或多个方向。所有后面的对图像数据的操作都相对于关键点的方向、尺度和位置进行变换，从而提供对于这些变换的不变性。

4. 关键点描述：在每个关键点周围的邻域内，在选定的尺度上测量图像局部的梯度。这些梯度被变换成一种表示，这种表示允许比较大的局部形状的变形和光照变化。

补充知识：

• 尺度空间理论
  尺度空间(scale space)思想最早是由Iijima于1962年提出的，后经witkin和Koenderink等人的推广逐渐得到关注，在计算机视觉邻域使用广泛。

尺度空间理论的基本思想是：在图像信息处理模型中引入一个被视为尺度的参数，通过连续变化尺度参数获得多尺度下的尺度空间表示序列，对这些序列进行尺度空间主轮廓的提取，并以该主轮廓作为一种特征向量，实现边缘、角点检测和不同分辨率上的特征提取等。

尺度空间方法将传统的单尺度图像信息处理技术纳入尺度不断变化的动态分析框架中，更容易获取图像的本质特征。尺度空间中各尺度图像的模糊程度逐渐变大，能够模拟人在距离目标由近到远时目标在视网膜上的形成过程。

尺度空间满足视觉不变性。该不变性的视觉解释如下：当我们用眼睛观察物体时，一方面当物体所处背景的光照条件变化时，视网膜感知图像的亮度水平和对比度是不同的，因此要求尺度空间算子对图像的分析不受图像的灰度水平和对比度变化的影响，即满足灰度不变性和对比度不变性。另一方面，相对于某一固定坐标系，当观察者和物体之间的相对位置变化时，视网膜所感知的图像的位置、大小、角度和形状是不同的，因此要求尺度空间算子对图像的分析和图像的位置、大小、角度以及仿射变换无关，即满足平移不变性、尺度不变性、欧几里德不变性以及仿射不变性。

• 尺度空间的表示
  一个图像的尺度空间，L(x,y,σ)定义为一个变化尺度的高斯函数G(x,y,σ)与原图像I(x,y)的卷积。公式如下：

其中，
m，n表示高斯模板的维度。(x, y)代表图像的像素位置。σ是尺度空间因子，值越小表示图像被平滑的越少，相应的尺度也就越小。大尺度对应于图像的概貌特征，小尺度对应于图像的细节特征。

• 高斯金字塔的构建：尺度空间在实现时使用高斯金字塔表示，高斯金字塔的构建分为两部分：

1. 对图像做不同尺度的高斯模糊；

2. 对图像做降采样(隔点采样)。

图像的金字塔模型是指，将原始图像不断降阶采样，得到一系列大小不一的图像，由大到小，从下到上构成的塔状模型。原图像为金子塔的第一层，每次降采样所得到的新图像为金字塔的一层(每层一张图像)，每个金字塔共n层。金字塔的层数根据图像的原始大小和塔顶图像的大小共同决定，其计算公式如下：

其中M，N为原图像的大小,t为塔顶图像的最小维数的对数值。

例如，对于大小为512×512的图像，金字塔上各层图像的大小如下表所示，当塔顶图像为4×4时，n=7，当塔顶图像为2 ×2时，n=8。

为了让尺度体现其连续性，高斯金字塔在简单降采样的基础上加上了高斯滤波。如图3.1所示，将图像金字塔每层的一张图像使用不同参数做高斯模糊，使得金字塔的每层含有多张高斯模糊图像，将金字塔每层多张图像合称为一组(Octave)，金字塔每层只有一组图像，组数和金字塔层数相等，使用公式(3-3)计算，每组含有多张(也叫层Interval)图像。另外，降采样时，高斯金字塔上一组图像的初始图像(底层图像)是由前一组图像的倒数第三张图像隔点采样得到的。

• 空间极值点检测(关键点的初步探查)
  关键点是由DoG空间的局部极值点组成的，关键点的初步探查是通过同一组内各DoG相邻两层图像之间比较完成的。为了寻找DoG函数的极值点，每一个像素点要和它所有的相邻点比较，看其是否比它的图像域和尺度域的相邻点大或者小。如下图所示，中间的检测点和它同尺度的8个相邻点和上下相邻尺度对应的9×2个点共26个点比较，以确保在尺度空间和二维图像空间都检测到极值点。

SIFT算法原理可以参见，https://blog.csdn.net/u010440456/article/details/81483145
这位大佬写的很详细了。

（1）兴趣点

SIFT 特征使用高斯差分函数来定位兴趣点：

其中，G_σ 是上一章中介绍的二维高斯核，I_σ 是使用G_σ 模糊的灰度图像，κ 是决定相差尺度的常数。兴趣点是在图像位置和尺度变化下 D(x,σ) 的最大值和最小值点。这些候选位置点通过滤波去除不稳定点。基于一些准则，比如认为低对比度和位于边上的点不是兴趣点，我们可以去除一些候选兴趣点。

（2）描述子

基于上面讨论的兴趣点位置描述子给出了兴趣点的位置和尺度信息。为了实现旋转不变性，基于每个点周围图像梯度的方向和大小，SIFT描述子引入了参考方向。SIFT描述子使用主方向描述参考方向，主方向使用方向直方图（以大小为权重）来度量。

（3）检测兴趣点

使用开源工具包VLFeat提供的二进制文件来计算图像SIFT特征。其中VLFeat工具包可以从 http://www.vlfeat.org/ 下载，二进制文件可以在所有主要的平台上运行。VLFeat库是用C语言写的，但是可以使用该库提供的命令行接口。

首先理解一下代码，把sift算法进行简单的调用：

from PIL import Image
from numpy import *
from pylab import *
import os

if __name__ == '__main__':
    imname = ('picture\jimei1.jpg')     #待处理图像路径
    im=Image.open(imname)
    process_image(imname,'test.sift')
    l1,d1 = read_features_from_file('test.sift')  #l1为兴趣点坐标、尺度和方位角度 l2是对应描述符的128 维向
    figure()
    gray()
    plot_features(im,l1,circle = True)
    title('sift-features')
    show()

实验分析：
简单运用sift算法得到如上结果，从中发现，sift算法角点检测和harris角点检测不同，选取的兴趣点也不一样。SIFT选取的对象会使用DoG检测关键点，并且对每个关键点周围的区域计算特征向量，它主要包括两个操作：检测和计算，操作的返回值是关键点信息和描述符，最后在图像上绘制关键点，并用imshow函数显示这幅图像。

修改并整合实验代码：

from PIL import Image
from numpy import *
from pylab import *
import os
import harris
# 添加中文字体支持
#from matplotlib.font_manager import FontProperties
#font = FontProperties(fname=r"c:\windows\fonts\simsunb.ttc", size=14)

def process_image(imagename,resultname,params="--edge-thresh 10 --peak-thresh 5"):
    """ Process an image and save the results in a file. """

    if imagename[-3:] != 'pgm':
        # create a pgm file
        im = Image.open(imagename).convert('L')  #.convert('L') 将RGB图像转为灰度模式，灰度值范围[0,255]
        im.save('tmp.pgm')                       #将灰度值图像信息保存在.pgm文件中
        imagename = 'tmp.pgm'
   
    cmmd = str("F:\Anaconda\chapter2\sift.exe "+imagename+" --output="+resultname+
                " "+params)
    os.system(cmmd)                              #执行sift可执行程序，生成resultname(test.sift)文件
    print ('processed', imagename, 'to', resultname)


def read_features_from_file(filename):
    """ Read feature properties and return in matrix form. """
    
    f = loadtxt(filename)
    return f[:,:4],f[:,4:] # feature locations, descriptors


def plot_features(im,locs,circle=True):
    """ Show image with features. input: im (image as array), 
        locs (row, col, scale, orientation of each feature). """

    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')

if __name__ == '__main__':
    imname = ('picture\house1.jpg')               #待处理图像路径
    im=Image.open(imname)
    process_image(imname,'test.sift')
    l1,d1 = read_features_from_file('test.sift')           #l1为兴趣点坐标、尺度和方位角度 l2是对应描述符的128 维向
    figure()
    gray()
    
    subplot(131)
    plot_features(im, l1, circle=False)
    title('(a)')
    
    subplot(132)
    plot_features(im,l1,circle = True)
    title('(b)')

实验分析：
可以从上图看出，选择圆圈表示SIFT特征尺度和不用圆圈表示SIFT尺度的检测出的是不一样的，非圆圈表示法没有能够清楚的显示出尺度大小，不能直观表示出该检测结果。

（4）匹配描述子

对于将一幅图像中的特征匹配到另一幅图像的特征，一种稳健的准则（同样是由 Lowe 提出的）是使用这两个特征距离和两个最匹配特征距离的比率。相比于图像中的其他特征，该准则保证能够找到足够相似的唯一特征。使用该方法可以使错误的匹配数降低。

更新：由于某些原因重新学习了以上的内容，因此增加了内容。

# -*- coding: utf-8 -*-
from PIL import Image
from pylab import *
from numpy 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("sift "+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_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 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 match_twosided(desc1, desc2):
    """双向对称版本的match()"""
    matches_12 = match(desc1, desc2)
    matches_21 = match(desc2, desc1)
 
    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):
    """返回将两幅图像并排拼接成的一幅新图像"""
    #选取具有最少行数的图像，然后填充足够的空行
    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 in range(len(matchscores)):
        if matchscores[i]>0:
            plot([locs1[i,0],locs2[matchscores[i,0],0]+cols1], [locs1[i,1],locs2[matchscores[i,0],1]],'c')
    axis('off')
 
# 示例
im1f = r'C:\Users\jiajia520\Pictures\Saved Pictures\house1.jpg'
im2f = r'C:\Users\jiajia520\Pictures\Saved Pictures\house2.jpg'
im1 = array(Image.open(im1f))
im2 = array(Image.open(im2f))
 
process_image(im1f, 'out_sift_1.txt')
l1,d1 = read_features_from_file('out_sift_1.txt')
figure()
gray()
subplot(121)
plot_features(im1, l1, circle=False)
 
process_image(im2f, 'out_sift_2.txt')
l2,d2 = read_featuyudares_from_file('out_sift_2.txt')
subplot(122)
plot_features(im2, l2, circle=False)
 
matches = match_twosided(d1, d2)
print('{} matches'.format(len(matches.nonzero()[0])))
figure()
gray()
plot_matches(im1,im2,l1,l2,matches, show_below=True)
show()

实验结果：



遇到的bug①：

OSError: out_sift_1.txt not found.

解决方案：

添加开源工具包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）进如Anaconda安装目录，找到Lib\site-packages\PCV\localdescriptors中的sift.py文件，打开(notepad或者记事本均可)，修改其中的cmmd路径
cmmd = str(r"C:\Users\XXX\sift.exe “+imagename+” --output="+resultname+ " "+params)（路径是你项目文件夹中的sift.exe的路径）记得在路径前加r，并且路径最后要加一个空格

遇到的bug②：

解决方案：

在这个链接下载缺失的文件：MSVCR100.dll
对应的解决方案，可以参考这个
找不到 vcomp100.dll错误的解决

（三）匹配地理标记图像

（1）从Panoramio下载地理标记图像

可以从谷歌提供的照片共享服务Panoramio（http://www.panoramio.com/ )获取地理标记图像，Panoramio提供一个API接口，方便易用。

# -*- coding: utf-8 -*-
import json
import os
import urllib
import urlparse
from PCV.tools.imtools import get_imlist
from pylab import *
from PIL import Image

# change the longitude and latitude here
# here is the longitude and latitude for Oriental Pearl
minx = '-77.037564'
maxx = '-77.035564'
miny = '38.896662'
maxy = '38.898662'

# number of photos
numfrom = '0'
numto = '20'
url = 'http://www.panoramio.com/map/get_panoramas.php?order=popularity&set=public&from=' + numfrom + '&to=' + numto + '&minx=' + minx + '&miny=' + miny + '&maxx=' + maxx + '&maxy=' + maxy + '&size=medium'

c = urllib.urlopen(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

# 显示下载到的20幅图像
figure()
gray()
filelist = get_imlist('./')
for i, imlist in enumerate(filelist):
    im = Image.open(imlist)
    subplot(4, 5, i + 1)
    imshow(im)
    axis('off')
show()

由于谷歌关闭了该服务，因此这个部分暂时就不进行叙述了。

（2）使用局部描述子匹配

对已经下载的图像，对这些图像提取局部描述子。

# -*- coding: utf-8 -*-
import json
import os
import urllib
#import urlparse
from pylab import *
from PIL import Image
from PCV.tools import imtools
import pydot
import harris

def process_image(imagename,resultname,params="--edge-thresh 10 --peak-thresh 5"):
    """ Process an image and save the results in a file. """

    if imagename[-3:] != 'pgm':
        # create a pgm file
        im = Image.open(imagename).convert('L')  #.convert('L') 将RGB图像转为灰度模式，灰度值范围[0,255]
        im.save('tmp.pgm')                       #将灰度值图像信息保存在.pgm文件中
        imagename = 'tmp.pgm'
   
    cmmd = str("F:\Anaconda\chapter2\sift.exe "+imagename+" --output="+resultname+
                " "+params)
    os.system(cmmd)                              #执行sift可执行程序，生成resultname(test.sift)文件
    print ('processed', imagename, 'to', resultname)


def read_features_from_file(filename):
    """ Read feature properties and return in matrix form. """
    
    f = loadtxt(filename)
    return f[:,:4],f[:,4:] # feature locations, descriptors


def plot_features(im,locs,circle=True):
    """ Show image with features. input: im (image as array), 
        locs (row, col, scale, orientation of each feature). """

    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_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 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])
            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
    
download_path = 'picture'
path = 'picture'

imlist = imtools.get_imlist(download_path)
nbr_images = len(imlist)

featlist = [imname[:-3] + 'sift' for imname in imlist]
for i, imname in enumerate(imlist):
    process_image(imname, featlist[i])

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])
        l1, d1 = read_features_from_file(featlist[i])
        l2, d2 = read_features_from_file(featlist[j])
        matches = 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]

分析： 将每一副图像间的匹配特征数保存在 matchscores 数组中。由于距离度量的对称性，可以在用matchscores矩阵填充。

运行匹配结果：

processed tmp.pgm to picture/mao\1.sift
processed tmp.pgm to picture/mao\2.sift
processed tmp.pgm to picture/mao\3.sift
processed tmp.pgm to picture/mao\4.sift
processed tmp.pgm to picture/mao\5.sift
processed tmp.pgm to picture/mao\6.sift
processed tmp.pgm to picture/mao\7.sift
processed tmp.pgm to picture/mao\8.sift
processed tmp.pgm to picture/mao\9.sift
comparing  picture/mao\1.jpg picture/mao\1.jpg
number of matches =  143
comparing  picture/mao\1.jpg picture/mao\2.jpg
number of matches =  10
comparing  picture/mao\1.jpg picture/mao\3.jpg
number of matches =  3
comparing  picture/mao\1.jpg picture/mao\4.jpg
number of matches =  54
comparing  picture/mao\1.jpg picture/mao\5.jpg
number of matches =  8
comparing  picture/mao\1.jpg picture/mao\6.jpg
number of matches =  10
comparing  picture/mao\1.jpg picture/mao\7.jpg
number of matches =  56
comparing  picture/mao\1.jpg picture/mao\8.jpg
number of matches =  10
comparing  picture/mao\1.jpg picture/mao\9.jpg
number of matches =  2
comparing  picture/mao\2.jpg picture/mao\2.jpg
number of matches =  183
comparing  picture/mao\2.jpg picture/mao\3.jpg
number of matches =  0
comparing  picture/mao\2.jpg picture/mao\4.jpg
number of matches =  11
comparing  picture/mao\2.jpg picture/mao\5.jpg
number of matches =  16
comparing  picture/mao\2.jpg picture/mao\6.jpg
number of matches =  183
comparing  picture/mao\2.jpg picture/mao\7.jpg
number of matches =  9
comparing  picture/mao\2.jpg picture/mao\8.jpg
number of matches =  15
comparing  picture/mao\2.jpg picture/mao\9.jpg
number of matches =  5
comparing  picture/mao\3.jpg picture/mao\3.jpg
number of matches =  127
comparing  picture/mao\3.jpg picture/mao\4.jpg
number of matches =  1
comparing  picture/mao\3.jpg picture/mao\5.jpg
number of matches =  1
comparing  picture/mao\3.jpg picture/mao\6.jpg
number of matches =  0
comparing  picture/mao\3.jpg picture/mao\7.jpg
number of matches =  1
comparing  picture/mao\3.jpg picture/mao\8.jpg
number of matches =  0
comparing  picture/mao\3.jpg picture/mao\9.jpg
number of matches =  0
comparing  picture/mao\4.jpg picture/mao\4.jpg
number of matches =  188
comparing  picture/mao\4.jpg picture/mao\5.jpg
number of matches =  9
comparing  picture/mao\4.jpg picture/mao\6.jpg
number of matches =  10
comparing  picture/mao\4.jpg picture/mao\7.jpg
number of matches =  134
comparing  picture/mao\4.jpg picture/mao\8.jpg
number of matches =  11
comparing  picture/mao\4.jpg picture/mao\9.jpg
number of matches =  0
comparing  picture/mao\5.jpg picture/mao\5.jpg
number of matches =  189
comparing  picture/mao\5.jpg picture/mao\6.jpg
number of matches =  16
comparing  picture/mao\5.jpg picture/mao\7.jpg
number of matches =  12
comparing  picture/mao\5.jpg picture/mao\8.jpg
number of matches =  13
comparing  picture/mao\5.jpg picture/mao\9.jpg
number of matches =  0
comparing  picture/mao\6.jpg picture/mao\6.jpg
number of matches =  183
comparing  picture/mao\6.jpg picture/mao\7.jpg
number of matches =  9
comparing  picture/mao\6.jpg picture/mao\8.jpg
number of matches =  15
comparing  picture/mao\6.jpg picture/mao\9.jpg
number of matches =  5
comparing  picture/mao\7.jpg picture/mao\7.jpg
number of matches =  221
comparing  picture/mao\7.jpg picture/mao\8.jpg
number of matches =  14
comparing  picture/mao\7.jpg picture/mao\9.jpg
number of matches =  1
comparing  picture/mao\8.jpg picture/mao\8.jpg
number of matches =  188
comparing  picture/mao\8.jpg picture/mao\9.jpg
number of matches =  0
comparing  picture/mao\9.jpg picture/mao\9.jpg
number of matches =  158

填充完整后的matchscores矩阵如下所示：

（3）可视化连接的图像

首先通过图像间是否具有匹配的局部描述子来定义图像间的连接，然后可视化这些连接情况。为了完成可视化，可以在图中显示这些图像，图的边代表连接。实验中将会使用pydot工具包 http://code.google.com/p/pydot/ ,该工具包是功能强大的GraphViz图形库的python接口。pydot使用pyparsing http://pyparsing.wikispaces.com/ 和GraphViz http://www.graphviz.org/ 。
为了创建显示可能图像组的图，如果匹配的数目高于一个阈值，我们使用边来连接相应的图像节点。为了得到图中的 图像，需要使用图像的全路径（在下面例子中，使用 path 变量表示）。为了使图像看起来漂亮，我们需要将每幅图像尺度化为缩略图形式，缩略图的最大边为 100 像素。

注意：在进行实验前要先安装完GraphViz，否则会出现错误！

graphviz实际上是一个绘图工具，可以根据dot脚本画出树形图等，十分方便。我们利用它可以轻松完成树形图等图案的绘制工作。原理其实很简单，利用python代码生成dot脚本，然后调用graphviz软件解析，生成一张图片。

# -*- coding: utf-8 -*-
import json
import os
import urllib
#import urlparse
from pylab import *
from PIL import Image
#from PCV.localdescriptors import sift
from PCV.tools import imtools
import pydot
import sift

# pydot需要绝对路径，路径分隔符为/而非\
download_path = 'picture/mao'
path = 'picture/mao'
 
# list of downloaded filenames
imlist = imtools.get_imlist(download_path)
nbr_images = len(imlist)

# extract features
featlist = [imname[:-3] + 'sift' for imname in imlist]
for i, imname in enumerate(imlist):
    sift.process_image(imname, featlist[i])

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])
        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
print ("The match scores is: \n", matchscores)

# 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]

# 可视化

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 = path + str(i) + '.jpg'
            im.save(filename)  # need temporary files of the right size
            g.add_node(pydot.Node(str(i), fontcolor='transparent', shape='rectangle', image=filename))

            # second image in pair
            im = Image.open(imlist[j])
            im.thumbnail((100, 100))
            filename = path + str(j) + '.jpg'
            im.save(filename)  # need temporary files of the right size需要一定大小的临时文件
            g.add_node(pydot.Node(str(j), fontcolor='transparent', shape='rectangle',image=filename)) 
            g.add_edge(pydot.Edge(str(i), str(j)))
            g.write_jpg('mao.jpg')

运行之后，所得到的图像下面这样的：

修改代码后，还是没有能够显示出图像：

由于不熟练，这部分学习进行了很长的时间，这部分询问了同学这个问题还是没有能够解决。待我解决了之后，再来完善该部分。

代码运行结果：