Programming Computer Vision with Python （学习笔记十）

摘要：如下图所示，左右是两张稍有不同的图片，但都包含了广州塔，左图红色框中标出了两个感兴趣的点，我期望找出它们在右图的对应位置即对应点。

现在考虑一个全景图拼接的应用场景，假设现有两张图片需要拼接成一张全景图，这两张图片是通过相机右转一定角度拍摄出来的，两张图片有部分取景是重叠的。如何实现拼接？当然这是一个不简单的问题，我们现在只考虑实现拼接目标的第一步：找出图像中重叠的内容，以及分别在两张图片中的位置。

如下图所示，左右是两张稍有不同的图片，但都包含了广州塔，左图红色框中标出了两个感兴趣的点，我期望找出它们在右图的对应位置（即对应点）。

首先，要确定检测哪些点，即哪些点是我们感兴趣的？这可以使用Harris角检测（见上篇笔记）方法来得到图像的角点集合，然后通过设置合适的阈值和坐标范围来找出我们感兴趣的点。有了两个图像的兴趣点集后，又如何能计算出它们的对应关系呢？这就需要解决两个问题：

兴趣点如何描述

兴趣点之间的对应关系如何计算

兴趣点，也即用Harris角检测出来的结果，它只有坐标和像素值，只有这些信息不足以用于匹配，无法从另一张图像中查找是否包含这个点。所以需要增加点的表征信息，一种方法是使用围绕点周围一小块的图像来描述这个点，如采用上图中所标记的方式，即：以兴趣点为中心划出一个小矩形，将区域内所有像素值以一向量进行存储，用这个向量来描述这个兴趣点，那么此向量称为兴趣点描述符（interest point descriptor，下简称IPD）。

下面实现一个函数，为所有角点生成IPD：

def get_desc(image, filtered_coords, wid = 5):
    #image为原图像，filtered_coords为角点的坐标，wid为矩形的“半径”
    desc = []
    for coords in filtered_coords:
        ipd = image[coords[0] - wid : coords[0] + wid + 1,
                coords[1] - wid : coords[1] + wid + 1].flatten()

        if ipd.shape[0] > 0:
            desc.append(ipd)

    return desc

如何确定左图中的某个兴趣点，对应右图中的某个兴趣点？对应关系，不一定完全是等价关系，即两个点虽然是对应关系，但它们对应的IPD并不完全相同。因为我们这里讨论的找对应点的方法，允许两张图像在亮度、缩放上有一定的区别。所以两个点的对应关系不能用IPD等价来匹配，而是要采用相似度或相关度来计算，相关度越高，它们越可能是对应关系。而相关度，可以使用现成的数学模型——皮尔逊相关系数（Pearson"s r，也被称为皮尔森相关系数r,下简称r系数）来表示。所以，计算两个点的对应关系就转化为计算两个IPD的r系数。

r系数被广泛用于度量两个变量之间的相关（线性相关）程度，它是两个变量之间的协方差和标准差的商，一种等价表达式为标准分的均值：
r =

I1和I2为样本集，μ1为I1的平均值，μ2为I2的平均值，σ1为I1的标准差，σ2为I2的标准差，上式计算结果即为r系数，范围为-1到1。 r系数为正且越大，表示I1和I2同时趋向于它们各自的平均值，变化方向一致，相关度越高。系数为0意味着两个变量之间没有线性关系。

把两个点对应的IPD代入上述公式的I，可得到两个点的相关程度。所以找两个图像之间兴趣点的对应关系，计算步骤是：

分别对两个图像应用Harris角检测，得到图像1的兴趣点集1，和图像2的兴趣点集2

设定IPD的矩形大小，计算所有兴趣点的IPD，得到IPD_SET_1和IPD_SET_2两个集合

设定r系数的阈值，如0.5，即相关度在[0.5,1]之间我们才考虑，那么，对IPD_SET_1中指定的某个IPD，计算它与IPD_SET_2中所有IPD的r系数，若最大的r系数落在[0.5,1]区间，则其对应的IPD是最相关的。

下面实现一个IPD匹配函数，传入两个IPD集合，找出所有r系数符合给定阈值的(即认为有对应关系的)IPD：

def match(desc1, desc2, threshold = 0.5):
    n = len(desc1[0])
    count1 = len(desc1)
    count2 = len(desc2)

    d = -np.ones((count1, count2)) #每个图1的IPD，其对应的力2的IPD下标初始化为-1
    for i in range(count1):
        ipd1 = desc1[i]
        d1 = (ipd1 - np.mean(ipd1)) / np.std(ipd1)
        for j in range(count2):
            ipd2 = desc2[j]
            if ipd1.shape[0] == ipd2.shape[0]: #忽略位于边缘的IPD
                d2 = (ipd2 - np.mean(ipd2)) / np.std(ipd2)
                r = np.sum(d1 * d2) / (n - 1)
                if r > threshold:
                    d[i, j] = r #i为图像1角点坐标, j为符合阈值的图像2角点坐标

    ndx = np.argsort(-d)  #将d的列降序排列，第0列即为r系数最大的
    match_index_array = ndx[:, 0] #只保留第0列
    return match_index_array

上述的函数为图1的每个IPD，从右边找到最佳的匹配（如果存在），但这还不够，因为这不代表对右边的这个IPD来说，左边的的这个IPD是它的最佳匹配，所以，如果使用两向匹配，互相认为是最佳的，我们才认为是对应关系，这样效果会更好一些，双向匹配的函数实现：

def match_twosided(desc1, desc2, threshold = 0.5):
    m_12 = match(desc1, desc2, threshold)
    m_21 = match(desc2, desc1, threshold)

    for i,j in enumerate(m_12):
        if j >= 0 and m_21[j] != i:
            m_12[i] = -1  #非双向匹配的，置为-1，上层应该忽略之

    return m_12

下面代码使用以上的两向匹配方法找出两张图像的对应点，并用白色线连接起来，看一下效果，两张图像是并排显示的：

从图中可以看出，两个图像中的广州塔上的关键角点基本能找到对应的位置，但图像的底部即建筑物的角点，与右图的建筑物连接起来，即使它们不是相同的建筑物，这是因为这些角点看起来很像，准确点讲，相关度（r系数）很高，所以被认为是对应点。

从这个例子也可以看出，要准确的找到对象在图像间的对应点，还需要考虑一些因素，来使效果更佳：

为兴趣点定义一个范围，比如上面例子，如果只关注塔尖的兴趣点，得出的效果令人满意

在寻找对应关系时，可限定对应点的y坐标的距离不能超过一定范围（如50个像素，根据实际应用而定），这样可以有效排除一些虽然r系数高，但事实上不对应的点。

代码如下，注意点的疏密可以通过参数微调：

from PIL import Image
import matplotlib.pyplot as plt
import numpy as np
from scipy.ndimage import filters
from skimage.feature import corner_peaks

def harris_eps(im, sigma=3): #harris角检测，见上个笔记
    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)
    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 * 2 / (Wtr + 1e-06)

# def get_desc(image, filtered_coords, wid = 5):
# 省略，见上文

# def match_twosided(desc1, desc2, threshold = 0.5):
# 省略，见上文

im1 = np.array(Image.open("tower-left.jpg").convert("L"))
im2 = np.array(Image.open("tower-right.jpg").convert("L"))

coords_1 = corner_peaks(harris_eps(im1, sigma=1), min_distance=3, threshold_abs=0, threshold_rel=0.1)
coords_2 = corner_peaks(harris_eps(im2, sigma=1), min_distance=3, threshold_abs=0, threshold_rel=0.1)

desc1 = get_desc(im1, coords_1, wid=6)
desc2 = get_desc(im2, coords_2, wid=6)

match_index_array = match_twosided(desc1, desc2, threshold=0.5)

im3 = np.concatenate((im1, im2), axis=1) #将两个图像左右合并成一个，以便显示
plt.gray()
plt.imshow(im3)

for ipd_index_1,ipd_index_2 in enumerate(match_index_array):
    if ipd_index_2 != -1:
        x = [coords_1[ipd_index_1][4], coords_2[ipd_index_2][5] + im1.shape[1]]
        y = [coords_1[ipd_index_1][0], coords_2[ipd_index_2][0]]
        if np.abs(y[0] - y[1]) < 50: #这里限制了对应点之间的y坐标距离
            plt.plot(x, y, "w", alpha=0.5) #连接两个对应点

plt.plot(coords_1[:, 1], coords_1[:, 0], "+r", markersize=5) #画图1角点坐标
plt.plot(coords_2[:, 1] + im1.shape[1], coords_2[:, 0], "+r", markersize=5)  #画图2角点坐标
plt.axis("off")
plt.show()

从实例中可以看到，本文使用的描述点的和匹配的方法，存在误配的情况，矩形大小的设置也会影响匹配的结果，而且它也不适用于在图像被旋转和缩放的情况下使用，近年，关于这方面的研究也在不断取得进步，下一笔记将介绍一种称为尺度不变特征转换(Scale-invariant feature transform 或 SIFT)的算法，此算法应用非常广。

你还可以查看我的其它笔记

wiki 皮尔逊积矩相关系数