基于Python的GDAL高分变率遥感影像配准算法

#!/usr/bin/env python # -*- encoding: utf-8 -*- ''' @File : utility.py @Time : 2022/05/25 18:17:10 @Author : qingchengboy @Version : 1.0 @Contact : 1441210996@qq.com @License : (C)Copyright 2021-2024 @Desc : None ''' # here put the import lib import numpy as np from sift_method import sift_matching from osgeo import gdal def block_matching(img1,img2,m,n,overlay=16): r,c,d = img1.shape b_r,b_c=r // m,c // n total_sift = np.array([]) for i in range(0, m): for j in range(0,n): #Use Rec's coordinates to capture images and match them start_rows = i *b_r - overlay start_cols = j*b_c - overlay end_rows = (i+1)*b_r + overlay end_cols = (j+1)*b_c + overlay if i == 0: start_rows = start_rows + overlay if i == m - 1: end_rows = end_rows - overlay if j== 0: start_cols = start_cols + overlay if j== n- 1: end_cols = end_cols - overlay #block matching l_img = img1[start_rows:end_rows, start_cols:end_cols, :] r_img = img2[start_rows:end_rows, start_cols:end_cols,:] if np.all(sift_matching(l_img,r_img) != 0): list_sift = sift_matching(l_img,r_img) #Combine n rows and 4 columns of the returned coordinate point pairs temp_sift = np.empty((list_sift.shape[0],4)) temp_sift[:,0] = list_sift[:,0] + start_cols temp_sift[:,1] = list_sift[:,1] + start_rows temp_sift[:,2] = list_sift[:,2] + start_cols temp_sift[:,3] = list_sift[:,3] + start_rows # Delete points within overlapping areas condition = (temp_sift[:, 1] > i * b_r) & (temp_sift[:, 1] <= (i + 1)* b_r) &(temp_sift[:,3]>i*b_r) & ( temp_sift[:,3]<= (i+1)*b_r) & (temp_sift[:,0]>j*b_c) & (temp_sift[:, 0] <= (j+1) *b_c) &( temp_sift[:,2]>j*b_c) & (temp_sift[:, 2] <= (j+ 1) * b_c) sift_new = temp_sift[np.where(condition)] total_sift=np.append(total_sift,sift_new) return total_sift def cps_write(pts,filename, ref_path, flag ='IMG2IMG'): if flag =='IMG2IMG': np.savetxt(filename,pts.reshape(-1,4),fmt="%8.6f") else: dataset = gdal.Open(ref_path) GeoTrans = dataset.GetGeoTransform() gcp_list =[] for i in range(pts.shape[0]): Xpixel = pts[i,0] Yline = pts[i,1] Xp = GeoTrans[0] + Xpixel * GeoTrans[1] + Yline * GeoTrans[2] Yp = GeoTrans[3] + Xpixel * GeoTrans[4] + Yline * GeoTrans[5] gcp_list.append([Xp,Yp,pts[i,2],pts[i,3]]) np.savetxt(filename, np.array(gcp_list).reshape(-1, 4), fmt="%8.9f",delimiter=',')