Python - opencv（十）特征匹配和圖像拼接

參考：

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

https://www.bilibili.com/video/BV1PV411774y?p=58

https://docs.opencv.org/3.0-beta/doc/py_tutorials/py_feature2d/py_feature_homography/py_feature_homography.html

手機的全景拍照功能可以將數張照片無縫憑借成一張長照片，用的便是特征匹配和圖像拼接的相關算法。本文介紹基於Python-opencv的實現

 

SIFT特征匹配

理論部分懶得碼字，搬運博客原文：

SIFT（Scale Invariant Feature Transform，尺度不變特征變換匹配算法）是由David G. Lowe教授在1999年（《Object Recognition from Local Scale-Invariant Features》）提出的高效區域檢測算法，在2004年（《Distinctive Image Features from Scale-Invariant Keypoints》）得以完善。

SIFT可以應用到物體辨識、機器人地圖感知與導航、影像縫合、3D模型建立、手勢辨識、影像追蹤和動作比對等方向。

SIFT算法的特點:

穩定性
獨特性
多量性
高速性
可擴展性
SIFT算法可以的解決問題：

目標的旋轉、縮放、平移（RST）
圖像放射/投影變換（視點viewpoint）
光照影響（illumination）
部分目標遮擋（occlusion）
雜物場景（clutter）
噪聲

SIFT原理

1.檢測尺度空間極值

檢測尺度空間極值就是搜索所有尺度上的圖像位置，通過高斯微分函數來識別對於尺度和旋轉不變的興趣點。其主要步驟可以分為建立高斯金字塔、生成DOG高斯差分金字塔和DOG局部極值點檢測。為了讓大家更清楚，我先簡單介紹一下尺度空間，再介紹主要步驟。

（1）尺度空間

一個圖像的尺度空間，定義為一個變化尺度的高斯函數與原圖像的卷積。即：

其中，*表示卷積計算。

其中，m、n表示高斯模版的維度，(x,y)代表圖像像素的位置。為尺度空間因子，值越小表示圖像被平滑的越少，相應的尺度就越小。小尺度對應於圖像的細節特征，大尺度對應於圖像的概貌特征，效果如下圖所示，尺度從左到右，從上到下，依次增大。

（2）建立高斯金字塔

尺度空間在實現時，使用高斯金字塔表示，高斯金字塔的構建分為兩部分：

1.對圖像做不同尺度的高斯模糊

2.對圖像做降采樣（隔點采樣）

圖像的金字塔模型是指，將原始圖像不斷降階采樣，得到一系列大小不一的圖像，由大到小，從下到上構成的塔狀模型。原圖像為金子塔的第一層，每次降采樣所得到的新圖像為金字塔的上一層(每層一張圖像)，每個金字塔共n層。金字塔的層數根據圖像的原始大小和塔頂圖像的大小共同決定。

為了讓尺度體現其連續性，高斯金字塔在簡單降采樣的基礎上加上了高斯濾波。如上圖所示，將圖像金字塔每層的一張圖像使用不同參數做高斯模糊，使得金字塔的每層含有多張高斯模糊圖像，將金字塔每層多張圖像合稱為一組(Octave)，金字塔每層只有一組圖像，組數和金字塔層數相等，每組含有多層Interval圖像。

高斯圖像金字塔共o組、s層， 則有：

其中，σ表示尺度空間坐標，s表示sub-level層坐標，表示初始尺度，S表示每組層數（一般為3~5）

（3）建立DOG高斯差分金字塔

為了有效提取穩定的關鍵點，利用不同尺度的高斯差分核與卷積生成。

DOG函數：

DOG在計算上只需相鄰高斯平滑后圖像相減，因此簡化了計算！

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

（4）DOG局部極值檢測

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

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

2.關鍵點的精確定位

以上方法檢測到的極值點是離散空間的極值點，以下通過擬合三維二次函數來精確確定關鍵點的位置和尺度，同時去除低對比度的關鍵點和不穩定的邊緣響應點(因為DOG算子會產生較強的邊緣響應)，以增強匹配穩定性、提高抗噪聲能力。

（1）關鍵點的精確定位

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

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

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

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

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

（2）去除邊緣響應

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

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

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

該值在兩特征值相等時達最小。Lowe論文中建議閾值T為1.2，即時保留關鍵點，反之剔除。

在Lowe的論文中，取r＝10。下圖右側為消除邊緣響應后的關鍵點分布圖。

3.關鍵點主方向分配

關鍵點主方向分配就是基於圖像局部的梯度方向，分配給每個關鍵點位置一個或多個方向。所有后面的對圖像數據的操作都相對於關鍵點的方向、尺度和位置進行變換，使得描述符具有旋轉不變性。

對於在DOG金字塔中檢測出的關鍵點，采集其所在高斯金字塔圖像3σ鄰域窗口內像素的梯度和方向分布特征。梯度的模值和方向如下：

L為關鍵點所在的尺度空間值，按Lowe的建議，梯度的模值m(x,y)按σ=1.5σ_oct的高斯分布加成，按尺度采樣的3σ原則，鄰域窗口半徑為3*1.5σ_oct。

在完成關鍵點的梯度計算后，使用直方圖統計鄰域內像素的梯度和方向。梯度直方圖將0~360度的方向范圍分為36個柱(bins)，其中每柱10度。如下圖所示，直方圖的峰值方向代表了關鍵點的主方向，(為簡化，圖中只畫了八個方向的直方圖)。

方向直方圖的峰值代表了該特征點處鄰域梯度的方向，以直方圖中最大值作為該關鍵點的主方向。為了增強匹配的魯棒性，只保留峰值大於主方向峰值80％的方向作為該關鍵點的輔方向。Lowe的論文指出大概有15%關鍵點具有多方向，但這些點對匹配的穩定性至為關鍵。檢測結果如下圖：

至此，將檢測出的含有位置、尺度和方向的關鍵點即是該圖像的SIFT特征點。

4.關鍵點的特征描述

通過以上步驟，對於每一個關鍵點，擁有三個信息：位置、尺度以及方向。接下來就是為每個關鍵點建立一個描述符，用一組向量將這個關鍵點描述出來，使其不隨各種變化而改變，比如光照變化、視角變化等。這個描述子不但包括關鍵點，也包含關鍵點周圍對其有貢獻的像素點，並且描述符應該有較高的獨特性，以便於提高特征點正確匹配的概率。 

SIFT描述子是關鍵點鄰域高斯圖像梯度統計結果的一種表示。通過對關鍵點周圍圖像區域分塊，計算塊內梯度直方圖，生成具有獨特性的向量，這個向量是該區域圖像信息的一種抽象，具有唯一性。

Lowe建議描述子使用在關鍵點尺度空間內4*4的窗口中計算的8個方向的梯度信息，共4*4*8=128維向量表征。表示步驟如下：

（1）計算描述子所需的圖像區域

（2）將坐標軸旋轉為關鍵點的方向

（3）將鄰域內的采樣點分配到對應的子區域內

（4）插值計算每個種子八個方向的梯度

（5）描述符向量元素門限化(閾值)

（6）描述符向量元素歸一化

 關鍵點匹配

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

代碼示例：

1. Brute-Force暴力匹配方法 

代碼：

import cv2

img1 = cv2.imread('../pics/part.jpg')
img2 = cv2.imread('../pics/whole.jpg')

sift = cv2.xfeatures2d.SIFT_create()

kp1, des1 = sift.detectAndCompute(img1, None)
kp2, des2 = sift.detectAndCompute(img2, None)

# 1v1 mapping
bf = cv2.BFMatcher(crossCheck=True)

matches = bf.match(des1, des2)
matches = sorted(matches, key=lambda x: x.distance)

img3 = cv2.drawMatches(img1, kp1, img2, kp2, matches[:10], None, flags=2)
cv2.imshow('match', img3)
cv2.waitKey()

# multi mapping
bf = cv2.BFMatcher()
matches = bf.knnMatch(des1, des2, k=2)
good = []
for m,n in matches:
    if m.distance < 0.75 * n.distance:
        good.append([m])
img4 = cv2.drawMatchesKnn(img1, kp1, img2, kp2, good, None, flags=2)
cv2.imshow('match k', img4)
cv2.waitKey()

(別人做的效果，本人實驗效果很差)

2. ransec 算法進行圖像拼接

因為沒有找到很好的素材做，本人自己拍了幾個樣張做拼接效果都很差，搬運一篇別人做好的：

https://www.jb51.net/article/183678.htm

代碼：

import cv2
import numpy as np
from matplotlib import pyplot as plt
import time
MIN = 10
starttime=time.time()
img1 = cv2.imread('1.jpg') #query
img2 = cv2.imread('2.jpg') #train

#img1gray=cv2.cvtColor(img1,cv2.COLOR_BGR2GRAY)
#img2gray=cv2.cvtColor(img2,cv2.COLOR_BGR2GRAY)
surf=cv2.xfeatures2d.SURF_create(10000,nOctaves=4,extended=False,upright=True)
#surf=cv2.xfeatures2d.SIFT_create()#可以改為SIFT
kp1,descrip1=surf.detectAndCompute(img1,None)
kp2,descrip2=surf.detectAndCompute(img2,None)

FLANN_INDEX_KDTREE = 0
indexParams = dict(algorithm = FLANN_INDEX_KDTREE, trees = 5)
searchParams = dict(checks=50)

flann=cv2.FlannBasedMatcher(indexParams,searchParams)
match=flann.knnMatch(descrip1,descrip2,k=2)

good=[]
for i,(m,n) in enumerate(match):
    if(m.distance<0.75*n.distance):
        good.append(m)

if len(good)>MIN:
    src_pts = np.float32([kp1[m.queryIdx].pt for m in good]).reshape(-1,1,2)
    ano_pts = np.float32([kp2[m.trainIdx].pt for m in good]).reshape(-1,1,2)
    M,mask=cv2.findHomography(src_pts,ano_pts,cv2.RANSAC,5.0)
    warpImg = cv2.warpPerspective(img2, np.linalg.inv(M), (img1.shape[1]+img2.shape[1], img2.shape[0]))
    direct=warpImg.copy()
    direct[0:img1.shape[0], 0:img1.shape[1]] =img1
    simple=time.time()

    #cv2.namedWindow("Result", cv2.WINDOW_NORMAL)
    #cv2.imshow("Result",warpImg)
    rows,cols=img1.shape[:2]
    left, right = (0,0)

    for col in range(0,cols):
        if img1[:, col].any() and warpImg[:, col].any():#開始重疊的最左端
            left = col
            break
    for col in range(cols-1, 0, -1):
        if img1[:, col].any() and warpImg[:, col].any():#重疊的最右一列
            right = col
            break

    res = np.zeros([rows, cols, 3], np.uint8)
    for row in range(0, rows):
        for col in range(0, cols):
            if not img1[row, col].any():#如果沒有原圖，用旋轉的填充
                res[row, col] = warpImg[row, col]
            elif not warpImg[row, col].any():
                res[row, col] = img1[row, col]
            else:
                srcImgLen = float(abs(col - left))
                testImgLen = float(abs(col - right))
                alpha = srcImgLen / (srcImgLen + testImgLen)
                res[row, col] = np.clip(img1[row, col] * (1-alpha) + warpImg[row, col] * alpha, 0, 255)

                warpImg[0:img1.shape[0], 0:img1.shape[1]]=res
                final=time.time()
                img3=cv2.cvtColor(direct,cv2.COLOR_BGR2RGB)
                plt.imshow(img3,),plt.show()
                img4=cv2.cvtColor(warpImg,cv2.COLOR_BGR2RGB)
                plt.imshow(img4,),plt.show()
                print("simple stich cost %f"%(simple-starttime))
                print("\ntotal cost %f"%(final-starttime))
                cv2.imwrite("simplepanorma.png",direct)
                cv2.imwrite("bestpanorma.png",warpImg)
else:
    print("not enough matches!")

效果：

運行結果

原圖1.jpg

原圖2.jpg

特征點匹配

直接拼接和平滑對比

效果