363 lines
13 KiB
Plaintext
363 lines
13 KiB
Plaintext
|
import numpy as np
|
|||
|
|
import cv2
|
|||
|
|
import math
|
|||
|
|
import os
|
|||
|
|
import pyproj
|
|||
|
|
# import gdal
|
|||
|
|
from osgeo import gdal, osr
|
|||
|
|
gdal.UseExceptions() # 显式启用 GDAL 异常
|
|||
|
|
from touying.ImageReproject_python.img_types import Point, POS, ImgEO, CamParas, ImgExif
|
|||
|
|
from touying.ImageReproject_python.dpal import DPAL
|
|||
|
|
|
|||
|
|
class ImageReproject:
|
|||
|
|
"""图像重投影类，实现目标定位、红线重投影和图像重投影功能"""
|
|||
|
|
|
|||
|
|
def __init__(self):
|
|||
|
|
"""初始化"""
|
|||
|
|
pass
|
|||
|
|
|
|||
|
|
def target_positioning(self, u, v, interval, gimbal_pitch, gimbal_yaw, gimbal_roll,
|
|||
|
|
h, l, b, cam_paras, img_width, img_height, dem_file,
|
|||
|
|
x=None, y=None, z=None):
|
|||
|
|
"""目标定位
|
|||
|
|
|
|||
|
|
Args:
|
|||
|
|
u: 图像x坐标
|
|||
|
|
v: 图像y坐标
|
|||
|
|
interval: 间隔
|
|||
|
|
gimbal_pitch: 云台俯仰角(弧度)
|
|||
|
|
gimbal_yaw: 云台偏航角(弧度)
|
|||
|
|
gimbal_roll: 云台横滚角(弧度)
|
|||
|
|
h: 高度
|
|||
|
|
l: 经度(弧度)
|
|||
|
|
b: 纬度(弧度)
|
|||
|
|
cam_paras: 相机参数
|
|||
|
|
img_width: 图像宽度
|
|||
|
|
img_height: 图像高度
|
|||
|
|
dem_file: DEM文件路径
|
|||
|
|
x: 输出经度
|
|||
|
|
y: 输出纬度
|
|||
|
|
z: 输出高度
|
|||
|
|
|
|||
|
|
Returns:
|
|||
|
|
(x, y, z): 目标点的经度、纬度和高度
|
|||
|
|
"""
|
|||
|
|
# 创建ImgExif对象并设置参数
|
|||
|
|
img_exif = ImgExif()
|
|||
|
|
img_exif.width = img_width
|
|||
|
|
img_exif.height = img_height
|
|||
|
|
img_exif.B = b
|
|||
|
|
img_exif.L = l
|
|||
|
|
img_exif.H = h
|
|||
|
|
img_exif.GimbalYaw = gimbal_yaw
|
|||
|
|
img_exif.GimbalPitch = gimbal_pitch
|
|||
|
|
img_exif.GimbalRoll = gimbal_roll
|
|||
|
|
img_exif.camPara = cam_paras
|
|||
|
|
|
|||
|
|
# 计算投影矩阵
|
|||
|
|
img_exif.com_zone()
|
|||
|
|
img_exif.computing_p_mat()
|
|||
|
|
|
|||
|
|
# 去畸变处理
|
|||
|
|
u_undistort, v_undistort = img_exif.undistortion_point(u, v)
|
|||
|
|
|
|||
|
|
# 创建UTM投影
|
|||
|
|
wgs84 = pyproj.CRS('EPSG:4326')
|
|||
|
|
utm_proj = pyproj.CRS(f'EPSG:{32600+img_exif.UTMzone}')
|
|||
|
|
transformer = pyproj.Transformer.from_crs(utm_proj, wgs84, always_xy=True)
|
|||
|
|
# 获取DEM的最小和最大高程
|
|||
|
|
minZ, maxZ = self.get_dem_min_max(dem_file)
|
|||
|
|
|
|||
|
|
planes = []
|
|||
|
|
h = minZ
|
|||
|
|
while h <= maxZ:
|
|||
|
|
planes.append([0.0, 0.0, 1.0, -h])
|
|||
|
|
h += interval
|
|||
|
|
|
|||
|
|
points = []
|
|||
|
|
for plane in planes:
|
|||
|
|
new_row = np.array([plane], dtype=np.float64)
|
|||
|
|
P_ = np.vstack([img_exif.P, new_row])
|
|||
|
|
mat_inv = np.linalg.inv(P_)
|
|||
|
|
img_pt = np.array([[u_undistort], [v_undistort], [1], [0]], dtype=np.float64)
|
|||
|
|
target_ptr = np.dot(mat_inv, img_pt)
|
|||
|
|
point = self.normalized(target_ptr)
|
|||
|
|
# utm转经纬度
|
|||
|
|
point = transformer.transform(point[0], point[1], point[2])
|
|||
|
|
points.append(point)
|
|||
|
|
|
|||
|
|
# 过滤不符合高程的点
|
|||
|
|
filtered_points = []
|
|||
|
|
for pt in points:
|
|||
|
|
projectZ = pt[2]
|
|||
|
|
demZ = self.get_dem_z(dem_file, pt[0], pt[1])
|
|||
|
|
diff = projectZ - demZ
|
|||
|
|
if abs(diff) <= interval:
|
|||
|
|
filtered_points.append(pt)
|
|||
|
|
|
|||
|
|
if len(filtered_points) == 0:
|
|||
|
|
print("未得到采样点！")
|
|||
|
|
return False, None
|
|||
|
|
if len(filtered_points) == 1:
|
|||
|
|
return filtered_points[0][0], filtered_points[0][1], filtered_points[0][2]
|
|||
|
|
else:
|
|||
|
|
cam_center = np.array([img_exif.EO.Xs, img_exif.EO.Ys, img_exif.EO.Zs])
|
|||
|
|
dists = [np.linalg.norm(np.array(pt) - cam_center) for pt in filtered_points]
|
|||
|
|
min_idx = int(np.argmin(dists))
|
|||
|
|
return filtered_points[min_idx][0], filtered_points[min_idx][1], filtered_points[min_idx][2]
|
|||
|
|
|
|||
|
|
def red_line_reproject(self, gimbal_pitch, gimbal_yaw, gimbal_roll, h, l, b,
|
|||
|
|
cam_paras, img_width, img_height, points, img_file, img_save_file):
|
|||
|
|
"""红线重投影
|
|||
|
|
|
|||
|
|
Args:
|
|||
|
|
gimbal_pitch: 云台俯仰角(弧度)
|
|||
|
|
gimbal_yaw: 云台偏航角(弧度)
|
|||
|
|
gimbal_roll: 云台横滚角(弧度)
|
|||
|
|
h: 高度
|
|||
|
|
l: 经度(弧度)
|
|||
|
|
b: 纬度(弧度)
|
|||
|
|
cam_paras: 相机参数
|
|||
|
|
img_width: 图像宽度
|
|||
|
|
img_height: 图像高度
|
|||
|
|
points: 点坐标列表
|
|||
|
|
img_file: 输入图像文件路径
|
|||
|
|
img_save_file: 输出图像文件路径
|
|||
|
|
|
|||
|
|
Returns:
|
|||
|
|
None
|
|||
|
|
"""
|
|||
|
|
# 创建ImgExif对象并设置参数
|
|||
|
|
img_exif = ImgExif()
|
|||
|
|
img_exif.width = img_width
|
|||
|
|
img_exif.height = img_height
|
|||
|
|
img_exif.B = b
|
|||
|
|
img_exif.L = l
|
|||
|
|
img_exif.H = h
|
|||
|
|
img_exif.GimbalYaw = gimbal_yaw
|
|||
|
|
img_exif.GimbalPitch = gimbal_pitch
|
|||
|
|
img_exif.GimbalRoll = gimbal_roll
|
|||
|
|
img_exif.camPara = cam_paras
|
|||
|
|
|
|||
|
|
# 计算投影矩阵
|
|||
|
|
img_exif.com_zone()
|
|||
|
|
img_exif.computing_p_mat()
|
|||
|
|
|
|||
|
|
# 创建UTM投影
|
|||
|
|
wgs84 = pyproj.CRS('EPSG:4326')
|
|||
|
|
utm_proj = pyproj.CRS(f'EPSG:{32600+img_exif.UTMzone}')
|
|||
|
|
transformer = pyproj.Transformer.from_crs(wgs84, utm_proj, always_xy=True)
|
|||
|
|
# 将点坐标转换为UTM坐标
|
|||
|
|
utm_points = []
|
|||
|
|
for point in points:
|
|||
|
|
# 将经纬度转换为UTM坐标
|
|||
|
|
x, y = transformer.transform(point.Y , point.X )
|
|||
|
|
utm_point = Point(x, y, point.Z)
|
|||
|
|
utm_points.append(utm_point)
|
|||
|
|
|
|||
|
|
# 读取图像
|
|||
|
|
img = cv2.imread(img_file)
|
|||
|
|
if img is None:
|
|||
|
|
print(f"无法读取图像文件: {img_file}")
|
|||
|
|
return
|
|||
|
|
pixel_lists=[]
|
|||
|
|
# 投影点到图像上
|
|||
|
|
for point in utm_points:
|
|||
|
|
# 构建三维点的齐次坐标
|
|||
|
|
point_3d = np.array([point.X, point.Y, point.Z, 1.0])
|
|||
|
|
|
|||
|
|
# 投影到图像平面
|
|||
|
|
point_2d = np.matmul(img_exif.P, point_3d)
|
|||
|
|
|
|||
|
|
# 归一化
|
|||
|
|
point_2d = point_2d / point_2d[2]
|
|||
|
|
|
|||
|
|
# 应用畸变
|
|||
|
|
u, v = img_exif.distortion_point(point_2d[0], point_2d[1])
|
|||
|
|
|
|||
|
|
#
|
|||
|
|
# # 检查点是否在图像范围内
|
|||
|
|
if 0 <= u < img_width and 0 <= v < img_height:
|
|||
|
|
pixel_lists.append({
|
|||
|
|
"u": int(u),
|
|||
|
|
"v": int(v)
|
|||
|
|
})
|
|||
|
|
# # 在图像上绘制点
|
|||
|
|
# cv2.circle(img, (int(u), int(v)), 5, (0, 0, 255), -1)
|
|||
|
|
|
|||
|
|
# 保存图像
|
|||
|
|
# cv2.imwrite(img_save_file, img)
|
|||
|
|
return pixel_lists
|
|||
|
|
def img_reproject(self, gimbal_pitch, gimbal_yaw, gimbal_roll, h, l, b,
|
|||
|
|
cam_paras, img_width, img_height, obj_longitude, obj_latitude, obj_height,
|
|||
|
|
pixel_x=None, pixel_y=None):
|
|||
|
|
"""图像重投影
|
|||
|
|
|
|||
|
|
Args:
|
|||
|
|
gimbal_pitch: 云台俯仰角(弧度)
|
|||
|
|
gimbal_yaw: 云台偏航角(弧度)
|
|||
|
|
gimbal_roll: 云台横滚角(弧度)
|
|||
|
|
h: 高度
|
|||
|
|
l: 经度(弧度)
|
|||
|
|
b: 纬度(弧度)
|
|||
|
|
cam_paras: 相机参数
|
|||
|
|
img_width: 图像宽度
|
|||
|
|
img_height: 图像高度
|
|||
|
|
obj_longitude: 目标点经度(度)
|
|||
|
|
obj_latitude: 目标点纬度(度)
|
|||
|
|
obj_height: 目标点高度
|
|||
|
|
pixel_x: 输出像素x坐标
|
|||
|
|
pixel_y: 输出像素y坐标
|
|||
|
|
|
|||
|
|
Returns:
|
|||
|
|
(pixel_x, pixel_y): 目标点在图像上的像素坐标
|
|||
|
|
"""
|
|||
|
|
# 创建ImgExif对象并设置参数
|
|||
|
|
img_exif = ImgExif()
|
|||
|
|
img_exif.width = img_width
|
|||
|
|
img_exif.height = img_height
|
|||
|
|
img_exif.B = b
|
|||
|
|
img_exif.L = l
|
|||
|
|
img_exif.H = h
|
|||
|
|
img_exif.GimbalYaw = gimbal_yaw
|
|||
|
|
img_exif.GimbalPitch = gimbal_pitch
|
|||
|
|
img_exif.GimbalRoll = gimbal_roll
|
|||
|
|
img_exif.camPara = cam_paras
|
|||
|
|
|
|||
|
|
# 计算投影矩阵
|
|||
|
|
img_exif.com_zone()
|
|||
|
|
img_exif.computing_p_mat()
|
|||
|
|
|
|||
|
|
# 创建UTM投影
|
|||
|
|
wgs84 = pyproj.CRS('EPSG:4326')
|
|||
|
|
utm_proj = pyproj.CRS(f'EPSG:{32600+img_exif.UTMzone}')
|
|||
|
|
transformer = pyproj.Transformer.from_crs(wgs84, utm_proj, always_xy=True)
|
|||
|
|
# 将目标点经纬度转换为UTM坐标
|
|||
|
|
obj_x, obj_y = transformer.transform(obj_longitude, obj_latitude)
|
|||
|
|
|
|||
|
|
# 构建三维点的齐次坐标
|
|||
|
|
point_3d = np.array([obj_x, obj_y, obj_height, 1.0])
|
|||
|
|
|
|||
|
|
# 投影到图像平面
|
|||
|
|
point_2d = np.matmul(img_exif.P, point_3d)
|
|||
|
|
|
|||
|
|
# 归一化
|
|||
|
|
point_2d = point_2d / point_2d[2]
|
|||
|
|
|
|||
|
|
# 应用畸变
|
|||
|
|
u, v = img_exif.distortion_point(point_2d[0], point_2d[1])
|
|||
|
|
|
|||
|
|
# 设置输出结果
|
|||
|
|
if pixel_x is not None:
|
|||
|
|
pixel_x = u
|
|||
|
|
if pixel_y is not None:
|
|||
|
|
pixel_y = v
|
|||
|
|
|
|||
|
|
return u, v
|
|||
|
|
|
|||
|
|
def get_dem_min_max(self, dem_file):
|
|||
|
|
"""获取DEM文件的最小和最大高程
|
|||
|
|
|
|||
|
|
Args:
|
|||
|
|
dem_file: DEM文件路径
|
|||
|
|
|
|||
|
|
Returns:
|
|||
|
|
(min_z, max_z): 最小和最大高程
|
|||
|
|
"""
|
|||
|
|
# 打开DEM文件
|
|||
|
|
dataset = gdal.Open(dem_file, gdal.GA_ReadOnly)
|
|||
|
|
if dataset is None:
|
|||
|
|
print(f"无法打开DEM文件: {dem_file}")
|
|||
|
|
return 0, 1000 # 返回默认值
|
|||
|
|
|
|||
|
|
# 获取波段
|
|||
|
|
band = dataset.GetRasterBand(1)
|
|||
|
|
|
|||
|
|
# 获取统计信息
|
|||
|
|
min_z, max_z, mean, std_dev = band.GetStatistics(True, True)
|
|||
|
|
|
|||
|
|
# 关闭数据集
|
|||
|
|
dataset = None
|
|||
|
|
|
|||
|
|
return min_z, max_z
|
|||
|
|
|
|||
|
|
def get_dem_z(self, dem_file, x, y):
|
|||
|
|
"""从DEM文件中获取指定点的高程
|
|||
|
|
|
|||
|
|
Args:
|
|||
|
|
dem_file: DEM文件路径
|
|||
|
|
x: UTM X坐标
|
|||
|
|
y: UTM Y坐标
|
|||
|
|
|
|||
|
|
Returns:
|
|||
|
|
z: 高程值，如果无法获取则返回None
|
|||
|
|
"""
|
|||
|
|
# 打开DEM文件
|
|||
|
|
dataset = gdal.Open(dem_file, gdal.GA_ReadOnly)
|
|||
|
|
if dataset is None:
|
|||
|
|
print(f"无法打开DEM文件: {dem_file}")
|
|||
|
|
return None
|
|||
|
|
|
|||
|
|
# 获取地理变换参数
|
|||
|
|
geotransform = dataset.GetGeoTransform()
|
|||
|
|
|
|||
|
|
# 获取波段
|
|||
|
|
band = dataset.GetRasterBand(1)
|
|||
|
|
|
|||
|
|
# 计算像素坐标
|
|||
|
|
px = int((x - geotransform[0]) / geotransform[1])
|
|||
|
|
py = int((y - geotransform[3]) / geotransform[5])
|
|||
|
|
|
|||
|
|
# 检查像素坐标是否在范围内
|
|||
|
|
if px < 0 or px >= dataset.RasterXSize or py < 0 or py >= dataset.RasterYSize:
|
|||
|
|
dataset = None
|
|||
|
|
return None
|
|||
|
|
|
|||
|
|
# 读取高程值
|
|||
|
|
data = band.ReadAsArray(px, py, 1, 1)
|
|||
|
|
if data is None:
|
|||
|
|
dataset = None
|
|||
|
|
return None
|
|||
|
|
|
|||
|
|
# 获取高程值
|
|||
|
|
z = data[0, 0]
|
|||
|
|
|
|||
|
|
# 检查是否为NoData值
|
|||
|
|
no_data_value = band.GetNoDataValue()
|
|||
|
|
if no_data_value is not None and z == no_data_value:
|
|||
|
|
dataset = None
|
|||
|
|
return None
|
|||
|
|
|
|||
|
|
# 关闭数据集
|
|||
|
|
dataset = None
|
|||
|
|
|
|||
|
|
return z
|
|||
|
|
def normalized(self, mat): # 归一化
|
|||
|
|
return [mat[0,0]/mat[3,0], mat[1,0]/mat[3,0], mat[2,0]/mat[3,0]]
|
|||
|
|
|
|||
|
|
# C风格接口
|
|||
|
|
def ir_target_positioning(u, v, interval, gimbal_pitch, gimbal_yaw, gimbal_roll,
|
|||
|
|
h, l, b, cam_paras, img_width, img_height, dem_file,
|
|||
|
|
x=None, y=None, z=None):
|
|||
|
|
"""目标定位的C风格接口"""
|
|||
|
|
reproject = ImageReproject()
|
|||
|
|
return reproject.target_positioning(u, v, interval, gimbal_pitch, gimbal_yaw, gimbal_roll,
|
|||
|
|
h, l, b, cam_paras, img_width, img_height, dem_file,
|
|||
|
|
x, y, z)
|
|||
|
|
|
|||
|
|
def ir_red_line_reproject(gimbal_pitch, gimbal_yaw, gimbal_roll, h, l, b,
|
|||
|
|
cam_paras, img_width, img_height, points, img_file, img_save_file):
|
|||
|
|
"""红线重投影的C风格接口"""
|
|||
|
|
reproject = ImageReproject()
|
|||
|
|
pixel_list=reproject.red_line_reproject(gimbal_pitch, gimbal_yaw, gimbal_roll, h, l, b,
|
|||
|
|
cam_paras, img_width, img_height, points, img_file, img_save_file)
|
|||
|
|
return pixel_list
|
|||
|
|
|
|||
|
|
def ir_img_reproject(gimbal_pitch, gimbal_yaw, gimbal_roll, h, l, b,
|
|||
|
|
cam_paras, img_width, img_height, obj_longitude, obj_latitude, obj_height,
|
|||
|
|
pixel_x=None, pixel_y=None):
|
|||
|
|
"""图像重投影的C风格接口"""
|
|||
|
|
reproject = ImageReproject()
|
|||
|
|
return reproject.img_reproject(gimbal_pitch, gimbal_yaw, gimbal_roll, h, l, b,
|
|||
|
|
cam_paras, img_width, img_height, obj_longitude, obj_latitude, obj_height,
|
|||
|
|
pixel_x, pixel_y)
|