Python SIFT特征匹配算法、匹配地理標記圖像

一、SIFT算法的原理

1、檢測適度空間極值

　　檢測尺度空間極值就是搜索所有尺度上的圖像位置，通過高斯微分函數來識別對於尺寸和旋轉不變的興趣點。其主要步驟可以分為建立高斯金字塔、生成DOG高斯差分金字塔和DOG局部極值點檢測。

　　1.1 尺度空間的構建

　　圖像的尺度空間是這幅圖像在不同解析度下的表示。一幅圖像可以產生幾組（octave）圖像，一組圖像包括幾層圖像。構造尺度空間傳統的方法即構造一個高斯金字塔，原始圖像作為最底層，然后對圖像進行高斯模糊再降采樣（2倍）作為下一層圖像（即尺度越大，圖像越模糊），循環迭代下去。

對圖像進行尺度變換，以滿足特征點的尺度不變性，保留圖像輪廓和細節。

 

　　1.2 DOG算子

　　DoG（Difference of Gaussian）函數：
　　
　　
　　該函數在計算上只需相鄰高斯平滑后圖像相減，因此簡化了計算。

 

　　1.2.1 DoG高斯差分金字塔

　　對應DOG算子，需構建DOG金字塔。
　　
　　可以通過高斯差分圖像看出圖像上的像素值變化情況。（如果沒有變化，也就沒有特征。特征必須是變化盡可能多的點。）DOG圖像描繪的是目標的輪廓。

 

　　1.3 DOG局部極值檢測

　　特征點是由DOG空間的局部極值點組成的。為了尋找DOG函數的極值點，每一個像素點要和它所有的相鄰點比較，看其是否比它的圖像域和尺度域 的相鄰點大或者小。

　　

　　中間的檢測點和它同尺度的8個相鄰點和上下相鄰尺度對應的9×2個 點共26個點比較，以確保在尺度空間和二維圖像空間都檢測到極值點。

 

2、關鍵點的精確定位

　　2.1 關鍵點精確定位

　　利用已知的離散空間點插值得到的連續空間極值點的方法叫做子像素插值（Sub-pixel Interpolation）。

 

 　　

　　為了提高關鍵點的穩定性，需要對尺度空間DOG函數進行曲線擬合。利用DOG函數在尺度空間的Taylor展開式(擬合函數)為：

 　　

　　其中，求導並讓方程等於零，可以得到極值點的偏移量為：

　　對應極值點，方程的值為：

　　其中, 代表相對插值中心的偏移量，當它在任一維度上的偏移量大於0.5時（即x或y或），意味着插值中心已經偏移到它的鄰近點上，所以必須改變當前關鍵點的位置。同時在新的位置上反復插值直到收斂；也有可能超出所設定的迭代次數或者超出圖像邊界的范圍，此時這樣的點應該刪除，在Lowe中進行了5次迭代。另外，|D(x)| 過小的點易受噪聲的干擾而變得不穩定，所以將小於某個經驗值(Lowe論文中使用0.03，Rob Hess等人實現時使用0.04/S)的極值點刪除。同時，在此過程中獲取特征點的精確位置(原位置加上擬合的偏移量)以及尺度()。

 

　　2.2 去除邊緣響應

　　由於DOG函數在圖像邊緣有較強的邊緣響應，因此需要排除邊緣響應DOG函數的峰值點在邊緣方向有較大的主曲率，而在垂直邊緣的方向有較小的主曲率。主曲率可以通過計算在該點位置尺度的2×2的Hessian矩陣得到，導數由采樣點相鄰差來估計：

 　　

　　Dxx表示DOG金字塔中某一尺度的圖像x方向求導兩次。

　　D的主曲率和H的特征值成正比。令 α ，β為特征值，則

 　　 

　　該值在兩特征值相等時達最小。

 

3、關鍵點的主方向分配

　　通過尺度不變性求極值點，可以使其具有縮放不變的性質。而利用關鍵點鄰域像素的梯度方向分布特性，可以為每個關鍵點指定方向參數方向，從而使描述子對圖像旋轉具有不變性通過求每個極值點的梯度來為極值點賦予方向。
像素點的梯度表示：
　　
梯度幅值：
　　

梯度方向：
　　

 

 

 

 

4、關鍵點的特征描述

　　下圖是一個SIFT描述子事例。其中描述子由2×2×8維向量表征，也即是2×2個8方向的方向直方圖組成。左圖的種子點由8×8單元組成。每一個小格都代表了特征點鄰域所在的尺度空間的一個像素，箭頭方向代表了像素梯度方向，箭頭長度代表該像素的幅值。然后在4×4的窗口內計算8個方向的梯度方向直方圖。繪制每個梯度方向的累加可形成一個種子點，如右圖所示：一個特征點由4個種子點的信息所組成。
　　
　　

二、SIFT算法關鍵點的匹配

　　分別對模板圖（參考圖，reference image）和實時圖（觀測圖， observation image）建立關鍵點描述子集合。目標的識別是通過兩點 集內關鍵點描述子的比對來完成。具有128維的關鍵點描述子的相似 性度量采用歐式距離。

　　模板圖中關鍵點描述子：

　　實時圖中關鍵點描述子：

　　任意兩描述子相似性度量：

　　要得到配對的關鍵點描述子， 需滿足 

　　關鍵點的匹配可以采用窮舉法來完成，但是這樣耗費的時間太多，一 般都采用kd樹的數據結構來完成搜索。搜索的內容是以目標圖像的關鍵點為基准，搜索與目標圖像的特征點最鄰近的原圖像特征點和次鄰近的原圖像特征點。

　　

 

三、SIFT算法的實現

　　計算圖像的SIFT特征，需要用到開源工具包VLFeat。下載地址：www.vlfeat.org。下載后將vl.dll和sift.exe放到項目的目錄下。

　　SIFT算法：

from pylab import *
from PIL import Image
import numpy as np
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 = np.loadtxt(filename)
    return f[:, :4], f[:, 4:]  # 特征位置，描述子


def plot_features(im, locs, circle=False):
    """顯示帶有特征的圖像
        輸入：im（數組圖像），locs（每個特征的行、列、尺度和方向角度）"""

    def draw_circle(c, r):
        t = np.arange(0, 1.01, .01) * 2 * np.pi
        x = r * np.cos(t) + c[0]
        y = r * np.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')
    return


def match(desc1, desc2):
    """對於第一幅圖像的每個描述子，選取其在第二幅圖像中的匹配
        輸入：desc1（第一幅圖像中的描述子），desc2（第二幅圖像中的描述子）"""

    desc1 = np.array([d / np.linalg.norm(d) for d in desc1])
    desc2 = np.array([d / np.linalg.norm(d) for d in desc2])

    dist_ratio = 0.6
    desc1_size = desc1.shape

    matchscores = np.zeros((desc1_size[0], 1), 'int')
    desc2t = desc2.T    # 預先計算矩陣轉置
    for i in range(desc1_size[0]):
        dotprods = np.dot(desc1[i, :], desc2t) # 向量點乘
        dotprods = 0.9999*dotprods
        # 反余弦和反排序，返回第二幅圖像中特征的索引
        index = np.argsort(np.arccos(dotprods))

        # 檢查最近鄰的角度是否小於dist_ratio乘以第二近鄰的角度
        if np.arccos(dotprods)[index[0]] < dist_ratio * np.arccos(dotprods)[index[1]]:
            matchscores[i] = int(index[0])

    return matchscores


def match_twosided(desc1, decs2):
    """雙向對稱版本的match"""

    matches_12 = match(desc1, decs2)
    matches_21 = match(decs2, 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):
    """返回將兩幅圖像並排拼接成的一幅新圖像"""

    # 選取具有最少行數的圖像，然后填充足夠的空行
    row1 = im1.shape[0]
    row2 = im2.shape[0]

    if row1 < row2:
        im1 = np.concatenate((im1, np.zeros((row2 - row1, im1.shape[1]))), axis=0)
    elif row1 > row2:
        im2 = np.concatenate((im2, np.zeros((row1 - row2, im2.shape[1]))), axis=0)

    # 如果這些情況都沒有，那么他們的行數相同，不需要進行填充
    return np.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 = np.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')


if __name__ == '__main__':

    im1f = r'IMG_01.jpg'
    im2f = r'IMG_03.jpg'
    im1 = np.array(Image.open(im1f))
    im2 = np.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_features_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()

　　運行結果：

　　

 

 

 

四、SIFT算法的優缺點

4.1 SIFT算法的優點：

　　圖像的局部特征，對旋轉、尺度縮放、亮度變化保持不變，對視角變化、仿射變換、噪聲也保持一定程度的穩定性。

　　獨特性好，信息量豐富，適用於海量特征庫進行快速、准確的匹配。

　　多量性，即使是很少幾個物體也可以產生大量的SIFT特征

　　高速性，經優化的SIFT匹配算法甚至可以達到實時性

　　擴招性，可以很方便的與其他的特征向量進行聯合

 

4.2 SIFT算法的缺點：

　　因為是通過對特征點構造128維的向量,然后對向量進行匹配,這樣圖像就得滿足足夠多的紋理,否則構造出的128維向量區別性就不是太大,容易造成誤匹配,極限情況如指紋圖像的匹配,星圖識別等這類圖像特征點周圍根本沒有什么紋理,這時SIFT算法就完全失效了。

 

五、地理標記圖像匹配

　　由於最后需要對匹配后的圖像進行連接可視化，所以需要用到pydot工具包。這個工具包需要先安裝配置Graphviz環境才能使用。

　　相關下載地址和安裝教程：

　　下載地址：https://graphviz.gitlab.io/_pages/Download/windows/graphviz-2.38.msi

　　教程：https://blog.csdn.net/lanchunhui/article/details/49472949

　　注：最好將圖片像素調整到100K左右或者更小，否則可能會出現 xxx.sift not found 或者關於數組維度的報錯。像素調小也會讓代碼運行時間縮短很多。

　　代碼實現：

from pylab import *
from PIL import Image
from PCV.localdescriptors import sift
from PCV.tools import imtools
import pydot

# 圖像文件夾地址
images_path = "D:/LearnSrc/PythonSrc/Homework/Lecture03/images"
path = "D:/LearnSrc/PythonSrc/Homework/Lecture03/images"

# 獲取圖像名 形成存放到列表中 計算列表長度
imlist = imtools.get_imlist(images_path)
nbr_images = len(imlist)

# 提取特征值
featlist = [imname[:-3] + 'sift' for imname in imlist]
for i, imname in enumerate(imlist):
    sift.process_image(imname, featlist[i])

matchscores = np.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)

# 復制值
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  # 創建鏈接所需的最小匹配數

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 = path + str(i) + '.png'
            im.save(filename)  # 轉成大小合適的臨時文件
            g.add_node(pydot.Node(str(i), fontcolor='transparent', shape='rectangle', image=filename))

            # 成對的第二個圖像
            im = Image.open(imlist[j])
            im.thumbnail((100, 100))
            filename = path + 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=filename))

            g.add_edge(pydot.Edge(str(i), str(j)))
g.write_png('connect-images.png')

　　運行結果：

 

　　

 

 

內容參考：

https://blog.csdn.net/qq_40369926/article/details/88597406

https://www.pianshen.com/article/1543280501/