摘要：如下圖所示，左右是兩張稍有不同的圖片，但都包含了廣州塔，左圖紅色框中標(biāo)出了兩個(gè)感興趣的點(diǎn)，我期望找出它們?cè)谟覉D的對(duì)應(yīng)位置即對(duì)應(yīng)點(diǎn)。

現(xiàn)在考慮一個(gè)全景圖拼接的應(yīng)用場(chǎng)景，假設(shè)現(xiàn)有兩張圖片需要拼接成一張全景圖，這兩張圖片是通過(guò)相機(jī)右轉(zhuǎn)一定角度拍攝出來(lái)的，兩張圖片有部分取景是重疊的。如何實(shí)現(xiàn)拼接？當(dāng)然這是一個(gè)不簡(jiǎn)單的問(wèn)題，我們現(xiàn)在只考慮實(shí)現(xiàn)拼接目標(biāo)的第一步：找出圖像中重疊的內(nèi)容，以及分別在兩張圖片中的位置。

如下圖所示，左右是兩張稍有不同的圖片，但都包含了廣州塔，左圖紅色框中標(biāo)出了兩個(gè)感興趣的點(diǎn)，我期望找出它們?cè)谟覉D的對(duì)應(yīng)位置（即對(duì)應(yīng)點(diǎn)）。

首先，要確定檢測(cè)哪些點(diǎn)，即哪些點(diǎn)是我們感興趣的？這可以使用Harris角檢測(cè)（見(jiàn)上篇筆記）方法來(lái)得到圖像的角點(diǎn)集合，然后通過(guò)設(shè)置合適的閾值和坐標(biāo)范圍來(lái)找出我們感興趣的點(diǎn)。有了兩個(gè)圖像的興趣點(diǎn)集后，又如何能計(jì)算出它們的對(duì)應(yīng)關(guān)系呢？這就需要解決兩個(gè)問(wèn)題：

興趣點(diǎn)如何描述

興趣點(diǎn)之間的對(duì)應(yīng)關(guān)系如何計(jì)算

興趣點(diǎn)，也即用Harris角檢測(cè)出來(lái)的結(jié)果，它只有坐標(biāo)和像素值，只有這些信息不足以用于匹配，無(wú)法從另一張圖像中查找是否包含這個(gè)點(diǎn)。所以需要增加點(diǎn)的表征信息，一種方法是使用圍繞點(diǎn)周圍一小塊的圖像來(lái)描述這個(gè)點(diǎn)，如采用上圖中所標(biāo)記的方式，即：以興趣點(diǎn)為中心劃出一個(gè)小矩形，將區(qū)域內(nèi)所有像素值以一向量進(jìn)行存儲(chǔ)，用這個(gè)向量來(lái)描述這個(gè)興趣點(diǎn)，那么此向量稱為興趣點(diǎn)描述符（interest point descriptor，下簡(jiǎn)稱IPD）。

下面實(shí)現(xiàn)一個(gè)函數(shù)，為所有角點(diǎn)生成IPD：

def get_desc(image, filtered_coords, wid = 5):
    #image為原圖像，filtered_coords為角點(diǎn)的坐標(biāo)，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

如何確定左圖中的某個(gè)興趣點(diǎn)，對(duì)應(yīng)右圖中的某個(gè)興趣點(diǎn)？對(duì)應(yīng)關(guān)系，不一定完全是等價(jià)關(guān)系，即兩個(gè)點(diǎn)雖然是對(duì)應(yīng)關(guān)系，但它們對(duì)應(yīng)的IPD并不完全相同。因?yàn)槲覀冞@里討論的找對(duì)應(yīng)點(diǎn)的方法，允許兩張圖像在亮度、縮放上有一定的區(qū)別。所以兩個(gè)點(diǎn)的對(duì)應(yīng)關(guān)系不能用IPD等價(jià)來(lái)匹配，而是要采用相似度或相關(guān)度來(lái)計(jì)算，相關(guān)度越高，它們?cè)娇赡苁菍?duì)應(yīng)關(guān)系。而相關(guān)度，可以使用現(xiàn)成的數(shù)學(xué)模型——皮爾遜相關(guān)系數(shù)（Pearson"s r，也被稱為皮爾森相關(guān)系數(shù)r,下簡(jiǎn)稱r系數(shù)）來(lái)表示。所以，計(jì)算兩個(gè)點(diǎn)的對(duì)應(yīng)關(guān)系就轉(zhuǎn)化為計(jì)算兩個(gè)IPD的r系數(shù)。

r系數(shù)被廣泛用于度量?jī)蓚€(gè)變量之間的相關(guān)（線性相關(guān)）程度，它是兩個(gè)變量之間的協(xié)方差和標(biāo)準(zhǔn)差的商，一種等價(jià)表達(dá)式為標(biāo)準(zhǔn)分的均值：
r =

I1和I2為樣本集，μ1為I1的平均值，μ2為I2的平均值，σ1為I1的標(biāo)準(zhǔn)差，σ2為I2的標(biāo)準(zhǔn)差，上式計(jì)算結(jié)果即為r系數(shù)，范圍為-1到1。 r系數(shù)為正且越大，表示I1和I2同時(shí)趨向于它們各自的平均值，變化方向一致，相關(guān)度越高。系數(shù)為0意味著兩個(gè)變量之間沒(méi)有線性關(guān)系。

把兩個(gè)點(diǎn)對(duì)應(yīng)的IPD代入上述公式的I，可得到兩個(gè)點(diǎn)的相關(guān)程度。所以找兩個(gè)圖像之間興趣點(diǎn)的對(duì)應(yīng)關(guān)系，計(jì)算步驟是：

分別對(duì)兩個(gè)圖像應(yīng)用Harris角檢測(cè)，得到圖像1的興趣點(diǎn)集1，和圖像2的興趣點(diǎn)集2

設(shè)定IPD的矩形大小，計(jì)算所有興趣點(diǎn)的IPD，得到IPD_SET_1和IPD_SET_2兩個(gè)集合

設(shè)定r系數(shù)的閾值，如0.5，即相關(guān)度在[0.5,1]之間我們才考慮，那么，對(duì)IPD_SET_1中指定的某個(gè)IPD，計(jì)算它與IPD_SET_2中所有IPD的r系數(shù)，若最大的r系數(shù)落在[0.5,1]區(qū)間，則其對(duì)應(yīng)的IPD是最相關(guān)的。

下面實(shí)現(xiàn)一個(gè)IPD匹配函數(shù)，傳入兩個(gè)IPD集合，找出所有r系數(shù)符合給定閾值的(即認(rèn)為有對(duì)應(yīng)關(guān)系的)IPD：

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

    d = -np.ones((count1, count2)) #每個(gè)圖1的IPD，其對(duì)應(yīng)的力2的IPD下標(biāo)初始化為-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角點(diǎn)坐標(biāo), j為符合閾值的圖像2角點(diǎn)坐標(biāo)

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

上述的函數(shù)為圖1的每個(gè)IPD，從右邊找到最佳的匹配（如果存在），但這還不夠，因?yàn)檫@不代表對(duì)右邊的這個(gè)IPD來(lái)說(shuō)，左邊的的這個(gè)IPD是它的最佳匹配，所以，如果使用兩向匹配，互相認(rèn)為是最佳的，我們才認(rèn)為是對(duì)應(yīng)關(guān)系，這樣效果會(huì)更好一些，雙向匹配的函數(shù)實(shí)現(xiàn)：

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，上層應(yīng)該忽略之

    return m_12

下面代碼使用以上的兩向匹配方法找出兩張圖像的對(duì)應(yīng)點(diǎn)，并用白色線連接起來(lái)，看一下效果，兩張圖像是并排顯示的：

從圖中可以看出，兩個(gè)圖像中的廣州塔上的關(guān)鍵角點(diǎn)基本能找到對(duì)應(yīng)的位置，但圖像的底部即建筑物的角點(diǎn)，與右圖的建筑物連接起來(lái)，即使它們不是相同的建筑物，這是因?yàn)檫@些角點(diǎn)看起來(lái)很像，準(zhǔn)確點(diǎn)講，相關(guān)度（r系數(shù)）很高，所以被認(rèn)為是對(duì)應(yīng)點(diǎn)。

從這個(gè)例子也可以看出，要準(zhǔn)確的找到對(duì)象在圖像間的對(duì)應(yīng)點(diǎn)，還需要考慮一些因素，來(lái)使效果更佳：

為興趣點(diǎn)定義一個(gè)范圍，比如上面例子，如果只關(guān)注塔尖的興趣點(diǎn)，得出的效果令人滿意

在尋找對(duì)應(yīng)關(guān)系時(shí)，可限定對(duì)應(yīng)點(diǎn)的y坐標(biāo)的距離不能超過(guò)一定范圍（如50個(gè)像素，根據(jù)實(shí)際應(yīng)用而定），這樣可以有效排除一些雖然r系數(shù)高，但事實(shí)上不對(duì)應(yīng)的點(diǎn)。

代碼如下，注意點(diǎn)的疏密可以通過(guò)參數(shù)微調(diào)：

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角檢測(cè)，見(jiàn)上個(gè)筆記
    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):
# 省略，見(jiàn)上文

# def match_twosided(desc1, desc2, threshold = 0.5):
# 省略，見(jiàn)上文

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) #將兩個(gè)圖像左右合并成一個(gè)，以便顯示
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: #這里限制了對(duì)應(yīng)點(diǎn)之間的y坐標(biāo)距離
            plt.plot(x, y, "w", alpha=0.5) #連接兩個(gè)對(duì)應(yīng)點(diǎn)

plt.plot(coords_1[:, 1], coords_1[:, 0], "+r", markersize=5) #畫(huà)圖1角點(diǎn)坐標(biāo)
plt.plot(coords_2[:, 1] + im1.shape[1], coords_2[:, 0], "+r", markersize=5)  #畫(huà)圖2角點(diǎn)坐標(biāo)
plt.axis("off")
plt.show()

從實(shí)例中可以看到，本文使用的描述點(diǎn)的和匹配的方法，存在誤配的情況，矩形大小的設(shè)置也會(huì)影響匹配的結(jié)果，而且它也不適用于在圖像被旋轉(zhuǎn)和縮放的情況下使用，近年，關(guān)于這方面的研究也在不斷取得進(jìn)步，下一筆記將介紹一種稱為尺度不變特征轉(zhuǎn)換(Scale-invariant feature transform 或 SIFT)的算法，此算法應(yīng)用非常廣。

你還可以查看我的其它筆記

wiki 皮爾遜積矩相關(guān)系數(shù)