图像拼接—SIFT+Flann匹配+估计单应矩阵—（全景图panorama）

开发环境：python+opencv3
笔者拼接图像的步骤如下：
step1: 利用特征算子检测并描述,常见的特征算子在cv2.xfeatures2中都有，比如角点harris,斑点surf,sift等，这些算子的原理不再过多阐述。
step1:对于描述子，要进行匹配，在cv2中匹配可以是暴力匹配也可以选择flann，这是优化过的匹配算法，至于暴力匹配就是sift 的作者lowe 提出的NNDR方法，遍历所有点。
step3: 对于匹配点，我们首先要检查是否存在误点，可以利用基本矩阵（即极线约束的方式来剔除），在笔者源码中，没有这一步。对于提纯后的同名点，我们可以估计一个单应矩阵，即点对点的变换方式（如果h33=1),那么单应矩阵也可称为“投影矩阵”，根据《计算机视觉在多视图几何中的应用》那本书可以知道，投影矩阵具有保线性，即可以消除射影失真（影像上的平行线在看起来不平行，它们都有消影点（vanishing point))。单应矩阵常用语拼接。->基本矩阵和单应矩阵的区别见笔者之前的文章。
step 4: 利用opencv 中的warpPerspective 函数来实现全景图。
下面附上代码和效果图：（代码都加了注释）

# -*- coding: utf-8 -*-
#Panorama Stitcher 
import cv2
import numpy as np
class macthing(object):
    def matchIMG(self,im1,im2,kp1,kp2,des1,des2):
        FLANN_INDEX_KDTREE=0
        index_p=dict(algorithm = FLANN_INDEX_KDTREE, trees = 5)
        searth_p=dict(checks=50)
        flann=cv2.FlannBasedMatcher(index_p,searth_p)
        matches=flann.knnMatch(des1,des2,k=2)
        good =[]
        pts1=[]
        pts2=[]
        for i,(m,n) in enumerate(matches):
            if m.distance<0.6*n.distance:
                good.append(m)
                pts1.append(kp1[m.queryIdx].pt)
                pts2.append(kp2[m.trainIdx].pt)
        pts1=np.float32(pts1)
        pts2=np.float32(pts2)
        H,mask=cv2.findHomography(pts1,pts2,cv2.RANSAC,0.01)
        pts1_1=pts1[mask.ravel()==1]
        pts2_2=pts2[mask.ravel()==1]
        return pts1_1,pts2_2,H
    def appendimage(self,im1,im2):
        r1=im1.shape[0]
        r2=im2.shape[0]
        if r1<r2:
            img=np.zeros((r2-r1,im1.shape[1]),np.uint8)
            im1_1=np.vstack((im1,img))
            im3=np.hstack((im1_1,im2))
        else:
            img=np.zeros((r1-r2,im2.shape[1]),np.uint8)
            im2_2=np.vstack((im2,img))
            im3=np.hstack((im1,im2_2))
        return im3 
#    def panorama_get(self,im1,im2,H):
#        if im1.shape[0]>=im2.shape[0]:
#            result=cv2.warpPerspective(im1,H,(im1.shape[1]+im2.shape[1],im1.shape[0]))
#            result[0:im2.shape[0],0:im2.shape[1]]=im2
#        else:
#            result=cv2.warpPerspective(im1,H,(im1.shape[1]+im2.shape[1],im2.shape[0]))
#            result[0:im2.shape[0],0:im2.shape[1]]=im2
#        return result
    def panorama_get(self,im1,im2,H):
        h1,w1=im1.shape[:2]
        h2,w2=im2.shape[:2]
        pts1=np.float32([[0,0],[0,h1],[w1,h1],[w1,0]]).reshape(-1,1,2)#转为3维坐标
        pts2=np.float32([[0,0],[0,h2],[w2,h2],[w2,0]]).reshape(-1,1,2)
        pts1_=cv2.perspectiveTransform(pts1,H)#H：3*3 矩阵，所以pts1也该为3维坐标
        pts=np.concatenate((pts1_,pts2),axis=0)#列连接
          #np.min 是行优先
        [xmin,ymin]=np.int32(pts.min(axis=0).ravel()-0.5)
        [xmax,ymax]=np.int32(pts.max(axis=0).ravel()+0.5)
        t=[-xmin,-ymin]# 左加右减
        Ht = np.array([[1,0,t[0]],[0,1,t[1]],[0,0,1]]) #相当于一个向右平移
        result = cv2.warpPerspective(im1, Ht.dot(H), (xmax-xmin, ymax-ymin))#最后一个参数是输出图像的宽、高
        result[t[1]:h2+t[1],t[0]:w2+t[0]] = im2
        return result
if __name__=="__main__":
    M=macthing()
    im1_=cv2.imread(r"C:\Users\Y\Desktop\3.jpg")
    im2_=cv2.imread(r"C:\Users\Y\Desktop\4.jpg")
    im1=cv2.cvtColor(im1_,cv2.COLOR_BGR2GRAY)
    im2=cv2.cvtColor(im2_,cv2.COLOR_BGR2GRAY)
    sift=cv2.xfeatures2d.SIFT_create()
    kp1,des1=sift.detectAndCompute(im1,None)
    kp2,des2=sift.detectAndCompute(im2,None)
    pts1_1,pts2_2,H=M.matchIMG(im1,im2,kp1,kp2,des1,des2)
    im3=M.appendimage(im1,im2)
    pts2_new=pts2_2.copy()
    for i in range(len(pts2_2)):
        pts2_new[i,0]=pts2_new[i,0]+np.float32(im1.shape[1])
    for i in range(len(pts1_1)):
        cv2.line(im3,tuple(pts1_1[i]),tuple(pts2_new[i]),255,2)
#    cv  cv2\Y\Desktop\45.jpg",result)
    cv2.namedWindow("panorma ,cv2.WINDOW_NORMAL)
    cv2.imshow("panorama",result),cv2.waitKey(0)

_]UI