import cv2
import numpy as np
import math
from tqdm import tqdm # 添加进度条库
class CDomPho:
def __init__(self):
"""初始化摄影测量参数"""
self.m_f = self.m_Pxz = self.m_k1 = self.m_k2 = self.m_p1 = self.m_p2 = 0.0
self.m_X0 = self.m_Y0 = 0.0
self.m_cols = self.m_rows = 0
self.m_aopX = self.m_aopY = self.m_aopZ = self.m_aopP = self.m_aopO = self.m_aopK = 0.0
self.m_rotM = np.zeros(9) # 旋转矩阵
self.m_ref = np.zeros(6) # 参考参数
def POK2M(self, p, o, k, r):
"""将俯仰角(p)、翻滚角(o)、偏航角(k)转换为旋转矩阵"""
cosp = math.cos(p)
cosw = math.cos(o)
cosk = math.cos(k)
sinp = math.sin(p)
sinw = math.sin(o)
sink = math.sin(k)
# 计算旋转矩阵元素
r[0] = cosp * cosk - sinp * sinw * sink
r[1] = -cosp * sink - sinp * sinw * cosk
r[2] = -sinp * cosw
r[3] = cosw * sink
r[4] = cosw * cosk
r[5] = -sinw
r[6] = sinp * cosk + cosp * sinw * sink
r[7] = -sinp * sink + cosp * sinw * cosk
r[8] = cosp * cosw
def SetPara(self, aopX, aopY, aopZ, aopP, aopO, aopK, f, pxz, cols, rows, k1, k2, p1, p2, X0, Y0):
"""设置摄影测量参数"""
self.m_rotM = np.zeros(9)
self.POK2M(aopP, aopO, aopK, self.m_rotM) # 计算旋转矩阵
# 存储参数
self.m_aopX, self.m_aopY, self.m_aopZ = aopX, aopY, aopZ
self.m_aopP, self.m_aopO, self.m_aopK = aopP, aopO, aopK
self.m_f, self.m_Pxz = f, pxz
self.m_k1, self.m_k2, self.m_p1, self.m_p2 = k1, k2, p1, p2
self.m_X0, self.m_Y0 = X0, Y0
self.m_cols, self.m_rows = cols, rows
# 设置参考参数
self.m_ref[0] = -(cols / 2) * pxz
self.m_ref[1] = pxz
self.m_ref[2] = 0
self.m_ref[3] = -(rows / 2) * pxz
self.m_ref[4] = 0
self.m_ref[5] = -pxz
def Scan2Pho(self, scX, scY):
"""从扫描坐标转换到摄影坐标"""
phoX = self.m_ref[0] + self.m_ref[1] * scX - self.m_ref[2] * scY
phoY = self.m_ref[3] + self.m_ref[4] * scX - self.m_ref[5] * scY
# 应用畸变校正
if abs(self.m_k1) > 1e-12:
x, y = phoX, phoY
r2 = x * x + y * y
r4 = r2 * r2
dx = x * (self.m_k1 * r2 + self.m_k2 * r4) + self.m_p1 * (r2 + 2 * x * x) + 2 * self.m_p2 * x * y
dy = y * (self.m_k1 * r2 + self.m_k2 * r4) + self.m_p2 * (r2 + 2 * y * y) + 2 * self.m_p1 * x * y
phoX -= dx
phoY -= dy
return phoX, phoY
def Pho2Scan(self, phoX, phoY):
"""从摄影坐标转换到扫描坐标"""
# 应用畸变校正
if abs(self.m_k1) > 1e-12:
x, y = phoX, phoY
r2 = x * x + y * y
r4 = r2 * r2
dx0 = x * (self.m_k1 * r2 + self.m_k2 * r4) + self.m_p1 * (r2 + 2 * x * x) + 2 * self.m_p2 * x * y
dy0 = y * (self.m_k1 * r2 + self.m_k2 * r4) + self.m_p2 * (r2 + 2 * y * y) + 2 * self.m_p1 * x * y
x = phoX + dx0
y = phoY + dy0
r2 = x * x + y * y
r4 = r2 * r2
dx = x * (self.m_k1 * r2 + self.m_k2 * r4) + self.m_p1 * (r2 + 2 * x * x) + 2 * self.m_p2 * x * y
dy = y * (self.m_k1 * r2 + self.m_k2 * r4) + self.m_p2 * (r2 + 2 * y * y) + 2 * self.m_p1 * x * y
phoX += dx
phoY += dy
phoX -= self.m_ref[0]
phoY -= self.m_ref[3]
tt = self.m_ref[1] * self.m_ref[5] - self.m_ref[2] * self.m_ref[4]
scX = (self.m_ref[5] * phoX - self.m_ref[2] * phoY) / tt
scY = -(-self.m_ref[4] * phoX + self.m_ref[1] * phoY) / tt
return scX, scY
def PhoZ2Grd(self, px, py, gz):
"""从摄影坐标和高程转换到地面坐标"""
A = self.m_rotM[0] * px + self.m_rotM[1] * py - self.m_rotM[2] * self.m_f
B = self.m_rotM[3] * px + self.m_rotM[4] * py - self.m_rotM[5] * self.m_f
C = self.m_rotM[6] * px + self.m_rotM[7] * py - self.m_rotM[8] * self.m_f
if abs(C) > 1e-12:
N = (gz - self.m_aopZ) / C
else:
N = 0.0
gx = self.m_aopX + A * N
gy = self.m_aopY + B * N
return gx, gy
def Grd2Pho_(self, gx, gy, gz):
"""从地面坐标转换到摄影坐标 - 修复数组处理问题"""
# 确保输入是标量
if isinstance(gx, np.ndarray) or isinstance(gy, np.ndarray) or isinstance(gz, np.ndarray):
# 处理数组输入
px = np.zeros_like(gx)
py = np.zeros_like(gy)
for i in range(len(gx)):
dx = gx[i] - self.m_aopX
dy = gy[i] - self.m_aopY
dz = gz[i] - self.m_aopZ
denom = self.m_rotM[2] * dx + self.m_rotM[5] * dy + self.m_rotM[8] * dz
if abs(denom) > 1e-12:
fm = -self.m_f / denom
else:
fm = 0.0
px[i] = (self.m_rotM[0] * dx + self.m_rotM[3] * dy + self.m_rotM[6] * dz) * fm
py[i] = (self.m_rotM[1] * dx + self.m_rotM[4] * dy + self.m_rotM[7] * dz) * fm
return px, py
else:
# 处理标量输入
dx = gx - self.m_aopX
dy = gy - self.m_aopY
dz = gz - self.m_aopZ
denom = self.m_rotM[2] * dx + self.m_rotM[5] * dy + self.m_rotM[8] * dz
if abs(denom) > 1e-12:
fm = -self.m_f / denom
else:
fm = 0.0
px = (self.m_rotM[0] * dx + self.m_rotM[3] * dy + self.m_rotM[6] * dz) * fm
py = (self.m_rotM[1] * dx + self.m_rotM[4] * dy + self.m_rotM[7] * dz) * fm
return px, py
def ImgZ2Grd(self, xyz):
"""从图像坐标和高程转换到地面坐标"""
results = []
for i in range(0, len(xyz), 3):
scX, scY, Z = xyz[i], xyz[i + 1], xyz[i + 2]
phoX, phoY = self.Scan2Pho(scX, scY)
gx, gy = self.PhoZ2Grd(phoX, phoY, Z)
results.extend([gx, gy, Z])
return results
def Grd2Img(self, gx, gy, gz):
"""从地面坐标转换到图像坐标 - 修复数组处理问题"""
# 确保输入是标量
if isinstance(gx, np.ndarray) or isinstance(gy, np.ndarray) or isinstance(gz, np.ndarray):
# 处理数组输入
scX = np.zeros_like(gx)
scY = np.zeros_like(gy)
for i in range(len(gx)):
px, py = self.Grd2Pho_(gx[i], gy[i], gz[i])
scX[i], scY[i] = self.Pho2Scan(px, py)
return scX, scY
else:
# 处理标量输入
px, py = self.Grd2Pho_(gx, gy, gz)
scX, scY = self.Pho2Scan(px, py)
return scX, scY
# 全局变量用于图像交互
scale = 1.0
angle = 0.0
offset_x = 0
offset_y = 0
last_mouse_pos = (0, 0)
dragging = False
domImg = None # 全局正射影像
def mouse_callback(event, x, y, flags, param):
"""鼠标回调函数,处理图像交互"""
global scale, angle, offset_x, offset_y, dragging, last_mouse_pos
if event == cv2.EVENT_MBUTTONDOWN: # 鼠标中键按下
dragging = True
last_mouse_pos = (x, y)
elif event == cv2.EVENT_MOUSEMOVE and dragging: # 拖动图像
dx = x - last_mouse_pos[0]
dy = y - last_mouse_pos[1]
offset_x += dx
offset_y += dy
last_mouse_pos = (x, y)
update_display()
elif event == cv2.EVENT_MBUTTONUP: # 鼠标中键释放
dragging = False
elif event == cv2.EVENT_MOUSEWHEEL: # 鼠标滚轮缩放
if flags > 0:
scale *= 1.1 # 放大
else:
scale /= 1.1 # 缩小
update_display()
elif event == cv2.EVENT_RBUTTONDOWN: # 鼠标右键旋转
angle += 10
if angle >= 360:
angle -= 360
update_display()
elif event == cv2.EVENT_LBUTTONDOWN: # 左键点击显示坐标
print(f"点击位置: x={x}, y={y}")
def update_display():
"""根据当前变换参数更新显示图像"""
global domImg, scale, angle, offset_x, offset_y
if domImg is None:
return
height, width = domImg.shape[:2]
center = (width // 2, height // 2)
# 创建变换矩阵
M_rot = cv2.getRotationMatrix2D(center, angle, scale)
M_translate = np.float32([[1, 0, offset_x], [0, 1, offset_y]])
# 组合变换:先旋转缩放,再平移
M_combined = np.vstack([M_translate, [0, 0, 1]]) @ np.vstack([M_rot, [0, 0, 1]])
M_final = M_combined[:2, :] # 转换为2x3矩阵用于warpAffine
# 应用变换并显示
transformed_img = cv2.warpAffine(domImg, M_final, (width, height))
# 添加坐标信息
cv2.putText(transformed_img, f"Scale: {scale:.2f} | Angle: {angle}°",
(10, 30), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (0, 255, 0), 2)
cv2.putText(transformed_img, f"Offset: ({offset_x}, {offset_y})",
(10, 60), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (0, 255, 0), 2)
cv2.imshow('Orthophoto', transformed_img)
def select_roi_custom(image):
"""自定义ROI选择函数,添加了辅助线"""
print("请用鼠标选择感兴趣区域,按空格键或回车确认,按ESC取消")
clone = image.copy()
roi_selected = False
roi_box = (0, 0, 0, 0)
drawing = False
ix, iy = -1, -1
def mouse_callback(event, x, y, flags, param):
nonlocal ix, iy, drawing, clone, roi_box, roi_selected
if event == cv2.EVENT_LBUTTONDOWN:
drawing = True
ix, iy = x, y
clone = image.copy()
elif event == cv2.EVENT_MOUSEMOVE:
if drawing:
clone = image.copy()
# 绘制矩形
cv2.rectangle(clone, (ix, iy), (x, y), (0, 255, 0), 2)
# 绘制中心十字线
center_x = (ix + x) // 2
center_y = (iy + y) // 2
cv2.line(clone, (center_x, 0), (center_x, clone.shape[0]), (0, 255, 255), 1)
cv2.line(clone, (0, center_y), (clone.shape[1], center_y), (0, 255, 255), 1)
# 显示尺寸信息
width = abs(ix - x)
height = abs(iy - y)
cv2.putText(clone, f"{width}x{height}", (x + 10, y + 20),
cv2.FONT_HERSHEY_SIMPLEX, 0.6, (0, 255, 255), 1)
cv2.imshow('Select ROI', clone)
elif event == cv2.EVENT_LBUTTONUP:
drawing = False
roi_box = (min(ix, x), min(iy, y), abs(ix - x), abs(iy - y))
clone = image.copy()
cv2.rectangle(clone, (roi_box[0], roi_box[1]),
(roi_box[0] + roi_box[2], roi_box[1] + roi_box[3]),
(0, 255, 0), 2)
cv2.imshow('Select ROI', clone)
roi_selected = True
cv2.namedWindow('Select ROI', cv2.WINDOW_NORMAL)
cv2.setMouseCallback('Select ROI', mouse_callback)
cv2.imshow('Select ROI', image)
while True:
key = cv2.waitKey(1) & 0xFF
if key == 13 or key == 32: # 回车或空格确认
break
elif key == 27: # ESC取消
roi_box = (0, 0, 0, 0)
break
cv2.destroyWindow('Select ROI')
if roi_selected:
x, y, w, h = roi_box
cropped = image[y:y + h, x:x + w]
print(f"已选择ROI: x={x}, y={y}, width={w}, height={h}")
return cropped, (x, y, w, h)
else:
print("未选择ROI,将使用整张图像")
return image, (0, 0, image.shape[1], image.shape[0])
def generate_dom(dompho, m_img, mingx, maxgx, mingy, maxgy, resolution, h=64.3):
"""生成正射影像的优化函数"""
dx = maxgx - mingx
dy = maxgy - mingy
domCols = int(dx / resolution) + 1
domRows = int(dy / resolution) + 1
# 创建空图像
domImg = np.zeros((domRows, domCols, 3), dtype=np.uint8)
# 获取ROI尺寸
rows, cols = m_img.shape[:2]
print(f"生成正射影像: {domCols}x{domRows} 像素 (分辨率: {resolution}米/像素)")
# 使用进度条
progress_bar = tqdm(total=domRows, desc="生成DOM")
for r in range(domRows):
gy = mingy + r * resolution
for c in range(domCols):
gx = mingx + c * resolution
# 计算图像坐标
ix, iy = dompho.Grd2Img(gx, gy, h)
i = int(ix + 0.5) # 列坐标
j = int(iy + 0.5) # 行坐标
# 确保坐标在图像范围内
if 0 <= i < cols and 0 <= j < rows:
domImg[r, c] = m_img[j, i] # 使用ROI内的图像数据
progress_bar.update(1)
progress_bar.close()
return domImg
if __name__ == "__main__":
# 读取摄影测量参数
try:
with open(r'D:\shuiwei\para.txt') as file:
file.readline() # 跳过标题行
params = list(map(float, file.readline().split()))
except Exception as e:
print(f"读取参数文件时出错: {e}")
exit(1)
# 解析参数
aopX, aopY, aopZ, aopP, aopO, aopK = params[0:6]
f, Pxz, k1, k2, p1, p2, X0, Y0 = params[6:14]
# 读取原始图像
original_img = cv2.imread(r'D:\shuiwei\8.jpg')
if original_img is None:
raise FileNotFoundError("无法读取图像文件")
# 使用自定义的ROI选择函数
m_img, roi_coords = select_roi_custom(original_img)
x0, y0, w, h = roi_coords
# 调整像主点偏移 - 关键修正!
X0_adjusted = X0 - x0
Y0_adjusted = Y0 - y0
# 获取ROI后的图像尺寸
rows, cols = m_img.shape[:2]
# 初始化摄影测量模型
dompho = CDomPho()
dompho.SetPara(aopX, aopY, aopZ, aopP, aopO, aopK, f, Pxz, cols, rows,
k1, k2, p1, p2, X0_adjusted, Y0_adjusted) # 使用调整后的像主点
# 定义图像四个角点及其高程(基于ROI)
corPt = [
0.0, 0.0, 64.3, # ROI左上角
cols, 0, 64.3, # ROI右上角
0, rows, 64.3, # ROI左下角
cols, rows, 64.3 # ROI右下角
]
# 计算地面坐标范围
grd_coords = dompho.ImgZ2Grd(corPt)
gxs = grd_coords[0::3]
gys = grd_coords[1::3]
mingx, maxgx = min(gxs), max(gxs)
mingy, maxgy = min(gys), max(gys)
# 创建正射影像
resolution = 3 # 地面分辨率(米/像素)
print("开始生成正射影像...")
# 使用优化函数生成DOM
domImg = generate_dom(dompho, m_img, mingx, maxgx, mingy, maxgy, resolution)
# 创建交互窗口
cv2.namedWindow('Orthophoto', cv2.WINDOW_NORMAL)
cv2.setMouseCallback('Orthophoto', mouse_callback)
update_display()
print("操作说明:")
print("- 鼠标中键拖动:平移图像")
print("- 鼠标滚轮:缩放图像")
print("- 鼠标右键:旋转图像")
print("- 鼠标左键:显示点击位置坐标")
print("- ESC键:退出程序")
print("- S键:保存当前视图")
# 主循环
while True:
key = cv2.waitKey(1)
if key == 27: # ESC键退出
break
elif key == ord('s'): # 添加保存功能
cv2.imwrite(r'D:\shuiwei\current_view1.jpg', domImg)
print("当前视图已保存")
cv2.destroyAllWindows()检查一下这个代码,几个坐标系换来换去很容易出错,尤其是opencv图像坐标系的坐标原点在图像左上角,而扫描坐标系坐标原点在左下角
最新发布