用python写tif遥感影像立体匹配算法并输出DEM

实现遥感影像立体匹配算法，可以采用以下步骤：

1. 读入左右两幅遥感影像，可以采用 Python 的 GDAL 库进行读取。

from osgeo import gdal

# 读取左右两幅遥感影像
left_ds = gdal.Open('left.tif')
right_ds = gdal.Open('right.tif')

# 获取影像的宽度和高度
width = left_ds.RasterXSize
height = left_ds.RasterYSize

2. 提取左右两幅影像的 SIFT 特征点，可以采用 OpenCV 库进行实现。

import cv2

# 读取左右两幅影像
left_img = cv2.imread('left.tif', cv2.IMREAD_GRAYSCALE)
right_img = cv2.imread('right.tif', cv2.IMREAD_GRAYSCALE)

# 创建 SIFT 特征提取器
sift = cv2.SIFT_create()

# 提取左右两幅影像的 SIFT 特征点
left_kp, left_desc = sift.detectAndCompute(left_img, None)
right_kp, right_desc = sift.detectAndCompute(right_img, None)

3. 对左右两幅影像的特征点进行匹配，可以采用 OpenCV 库的 `FlannBasedMatcher` 进行实现。

# 创建 Flann 匹配器
matcher = cv2.FlannBasedMatcher()

# 对左右两幅影像的特征点进行匹配
matches = matcher.knnMatch(left_desc, right_desc, k=2)

4. 进行 RANSAC 算法进行离群点去除和求解基础矩阵，可以采用 OpenCV 库的 `findFundamentalMat` 函数进行实现。

# 进行 RANSAC 算法进行离群点去除和求解基础矩阵
left_pts = []
right_pts = []
for m1, m2 in matches:
    if m1.distance < 0.75 * m2.distance:
        left_pts.append(left_kp[m1.queryIdx].pt)
        right_pts.append(right_kp[m1.trainIdx].pt)

left_pts = np.int32(left_pts)
right_pts = np.int32(right_pts)

F, mask = cv2.findFundamentalMat(left_pts, right_pts, cv2.FM_RANSAC, 0.1, 0.99)

5. 根据基础矩阵和左右两幅影像的特征点，进行立体匹配，可以采用 OpenCV 库的 `cv2.reprojectImageTo3D` 函数进行实现。

# 根据基础矩阵和左右两幅影像的特征点，进行立体匹配
disp = np.zeros_like(left_img)
disp[left_pts[:, 1], left_pts[:, 0]] = cv2.computeDisparity(left_img, right_img)

Q = np.float32([[1, 0, 0, -width/2],
                [0, -1, 0, height/2],
                [0, 0, 0, -focal_length],
                [0, 0, 1, 0]])

points_3D = cv2.reprojectImageTo3D(disp, Q)

6. 将匹配得到的三维点云转换为 DEM，可以采用 Python 的 GDAL 库进行实现。

# 将匹配得到的三维点云转换为 DEM
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create('output_dem.tif', width, height, 1, gdal.GDT_Float32)
out_ds.SetGeoTransform(left_ds.GetGeoTransform())
out_ds.SetProjection(left_ds.GetProjection())

out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(points_3D[:, :, 2])
out_band.SetNoDataValue(-9999)

out_ds.FlushCache()
out_ds = None

完整代码如下：