（Python）Gdal与Opencv遥感影像Sift匹配+RANSAC筛选并计算匹配RMSE精度

简要介绍

SIFT，即尺度不变特征变换（Scale-invariant feature transform，SIFT），是用于图像处理领域的一种描述。这种描述具有尺度不变性，可在图像中检测出关键点，是一种局部特征描述子。 该方法于1999年由David Lowe 首先发表于计算机视觉国际会议（International Conference on Computer Vision，ICCV），2004年再次经David Lowe整理完善后发表于International journal of computer vision（IJCV）。
不仅在近景影像，SIFT在遥感影像也获得了极大的成功，本文旨在提供一个Demo，使用Python版本的Gdal读取遥感影像然后使用Opencv中的Sift等算法进行处理，还有给出一个计算匹精度RMSE的小工具。

代码

调用Opencv的SIFT进行匹配，使用单应性的RANSAC进行误匹配筛除。

from osgeo import gdal
import numpy as np
import cv2
import cmath
from matplotlib import pyplot as plt
from time import *
def readTIFF(tifpath, bandnum):
    """
    Use GDAL to read data and transform them into arrays.
    :param tifpath:tif文件的路径
    :param bandnum:需要读取的波段
    :return:该波段的数据，narray格式。len(narray)是行数，len(narray[0])列数
    """
    image = gdal.Open(tifpath)  # 打开影像
    if image == None:
        print(tifpath + "该tif不能打开!")
        return
    im_width = image.RasterXSize  # 栅格矩阵的列数
    im_height = image.RasterYSize  # 栅格矩阵的行数
    im_bands = image.RasterCount  # 波段数
    im_proj = image.GetProjection()  # 获取投影信息坐标系
    im_geotrans = image.GetGeoTransform()  # 仿射矩阵
    print('tif数据:{}个行，{}个列，{}层波段, 取出第{}层.'.format(im_width, im_height, im_bands, bandnum))
    im_data = image.ReadAsArray(0, 0,  im_width, im_height)
    del image  # 减少冗余
    return im_data,im_proj, im_geotrans
def normalization(data):
    _range = np.max(data) - np.min(data)
    return (data - np.min(data)) / _range
def Tiff16to8bit(img_16):
    if (np.max(img_16) - np.min(img_16) != 0):
        # img_nrm = (img_16 - np.min(img_16)) / (np.max(img_16) - np.min(img_16)) #计算灰度范围,归一化
        img_nrm = normalization(img_16)
        img_8 = np.uint8(255 * img_nrm)
    return img_8

def imagexy2geo(trans, row, col):
    px = trans[0] + col * trans[1] + row * trans[2]
    py = trans[3] + col * trans[4] + row * trans[5]
    return px, py
def geo2imagexy(trans, x, y):
    a = np.array([[trans[1], trans[2]], [trans[4], trans[5]]])
    b = np.array([x - trans[0], y - trans[3]])
    return np.linalg.solve(a, b)  # 使用numpy的linalg.solve进行二元一次方程的求解

def SIFT(img_l, img_r):
    sift = cv2.SIFT_create()
    kp1, des1 = sift.detectAndCompute(img_l, None)
    kp2, des2 = sift.detectAndCompute(img_r, None)
    # 创建设置FLANN匹配
    FLANN_INDEX_KDTREE = 0
    index_params = dict(algorithm=FLANN_INDEX_KDTREE, trees=5)
    search_params = dict(checks=50)
    flann = cv2.FlannBasedMatcher(index_params, search_params)
    matches = flann.knnMatch(des1, des2, k=2)
    # store all the good matches as per Lowe's ratio test.
    good = []
    # 舍弃大于0.7的匹配,初步筛除
    for m, n in matches:
        if m.distance < 0.7 * n.distance:
            good.append(m)
    MIN_MATCH_COUNT = 10  # 设置最低特征点匹配数量为10
    if len(good) > MIN_MATCH_COUNT:
        # 获取关键点的坐标
        src_pts = np.float32([kp1[m.queryIdx].pt for m in good]).reshape(-1, 1, 2)
        dst_pts = np.float32([kp2[m.trainIdx].pt for m in good]).reshape(-1, 1, 2)
        # 计算变换矩阵和MASK
        M, mask = cv2.findHomography(src_pts, dst_pts, cv2.RANSAC, 3.0)
        #M, mask = cv2.findFundamentalMat(src_pts, dst_pts, cv2.FM_RANSAC, 1.0)
        matchesMask = mask.ravel().tolist()
        calRMSE(src_pts, dst_pts, M, mask)#计算匹配算的精度)
    else:
        print("Not enough matches are found - %d/%d" % (len(good), MIN_MATCH_COUNT))
        matchesMask = None
    draw_params = dict(matchColor=(0, 255, 255),
                       singlePointColor=None,
                       matchesMask=matchesMask,
                       flags=2)
    row_l, col_l = img_l.shape
    row_r, col_r = img_r.shape
    img_show = np.empty((max(row_l, row_r), col_l + col_r))
    img_show = cv2.drawMatches(img_l,kp1,img_r,kp2,good,None,**draw_params)
    #plt.imshow(img_show), plt.show()
def calRMSE(src_pts,dst_pts,M,mask):
    # 求残差
    sum_H = 0 #残差和
    num = 0 #参与统计的总个数
    for i, j, m in zip(src_pts, dst_pts, mask):
        P_src = np.float32([i[0][0],i[0][1],1]).reshape((-1, 1))
        P = np.matmul(M, P_src) #通过计算出的矩阵预测点
        p = np.float32([P[0] / P[2], P[1] / P[2]]) #从齐次矩阵变为2维点
        j = j.T
        distance = np.linalg.norm(p - j)
        if (m == True):
            sum_H += distance
            num += 1
    rmse = cmath.sqrt(sum_H/num)
    print("rmse : ",rmse)
    return rmse
B,B_proj, B_geotrans = readTIFF('D:\\SongshanZY3\\ortho\\BWDSC.tif', 1)
F,F_proj, F_geotrans = readTIFF('D:\\SongshanZY3\\ortho\\FWDSC.tif', 1)
x_offset,y_offset = 5000,6000
block_size_x = 3000
block_size_y = 3000
block_B = B[y_offset:y_offset + block_size_y,x_offset:x_offset + block_size_x]
block_F = F[y_offset:y_offset + block_size_y,x_offset:x_offset + block_size_x]
block_B = Tiff16to8bit(block_B)
block_F = Tiff16to8bit(block_F)
begin_time = time()
SIFT(block_B,block_F)
end_time = time()
run_time = end_time-begin_time
print ('匹配耗时运行时间：',run_time,'s') #该循环程序运行时间： 1.4201874732

运行结果

裁出了3000*3000的小块，匹配精度RMSE为1.16像素，thinkPadX1轻薄本耗时29s。