import cv2
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
# 设置Matplotlib中文字体
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
# 相机内参和畸变系数
K = np.array([
[1881.4773 , 0 ,1016.363887],
[0 , 1931.988488 , 547.1170723],
[0 , 0 , 1]], dtype=np.float32)
dist = np.array([-0.927214665,0.875249537,0.001301429,-0.006101325, 0], dtype=np.float32)
# 控制点坐标
world_points_3d = np.array([
[4262845.751, 644463.054, 64.289],
[4262844.149, 644468.755, 63.701],
[4262846.106, 644470.001, 63.621],
[4262847.099, 644471.268, 63.765],
[4262853.252, 644455.049, 64.177],
[4262846.305, 644431.715, 66.42],
[4262851.769, 644434.948, 64.503],
[4262852.526, 644431.694, 64.784],
[4262858.526, 644426.641, 64.81],
[4262865.608, 644419.362, 66.008],
[4262858.634, 644475.401, 63.64],
[4262868.267, 644480.64, 63.645],
[4262871.129, 644477.184, 63.719]
], dtype=np.float32)
image_points = np.array([
[205, 357],
[31, 462],
[74, 485],
[88, 501],
[593, 288],
[525, 145],
[644, 196],
[674, 178],
[815, 164],
[968, 130],
[748, 599],
[1697, 763],
[1739, 616]
], dtype=np.float32)
# 读取图像
raw_img = cv2.imread(r"D:\shuiwei\1.jpg")
if raw_img is None:
raise FileNotFoundError("无法加载图像,请检查路径: D:\\shuiwei\\1.jpg")
# 畸变校正
undistorted_img = cv2.undistort(raw_img, K, dist)
# 修正后的旋转矩阵和平移向量
R = np.array([
[-0.99562414, -0.08557365, -0.03754638],
[-0.01233105, -0.27796685, 0.96051152],
[-0.09263113, 0.95677144, 0.27569529]
], dtype=np.float32)
tvec = np.array([[4299360.374, 231647.0009, -221780.8828]], dtype=np.float32)
tvec = tvec.reshape(3, 1) # 确保为3x1向量
# 计算水面平均高程
Z0 = np.mean(world_points_3d[:, 2])
print(f"水面平均高程: {Z0:.3f} 米")
# 计算反投影函数
def inverse_project(points_2d, Z0):
R_inv = np.linalg.inv(R)
K_inv = np.linalg.inv(K)
points_3d = []
for pt in points_2d:
# 归一化坐标
uv = np.array([pt[0], pt[1], 1.0], dtype=np.float32)
xyz_cam = K_inv @ uv
# 计算射线方向
ray_dir = R_inv @ xyz_cam
# 计算与高程平面相交的参数t
# 方程: t * ray_dir[2] + (R_inv @ (-tvec))[2] = Z0
t = (Z0 - (R_inv @ (-tvec))[2, 0]) / ray_dir[2]
# 计算世界坐标
P_c = t * xyz_cam
P_w = R_inv @ (P_c - tvec)
points_3d.append(P_w.flatten())
return np.array(points_3d)
# 使用图像边界反投影
img_height, img_width = undistorted_img.shape[:2]
img_corners = np.array([
[0, 0], # 左上
[0, img_height], # 左下
[img_width, img_height], # 右下
[img_width, 0] # 右上
], dtype=np.float32)
# 反投影得到世界坐标
world_corners = inverse_project(img_corners, Z0)
# 提取XY范围
X_min, X_max = np.min(world_corners[:, 0]), np.max(world_corners[:, 0])
Y_min, Y_max = np.min(world_corners[:, 1]), np.max(world_corners[:, 1])
# 扩展5%边界缓冲
x_range = X_max - X_min
y_range = Y_max - Y_min
X_min -= 0.05 * x_range
X_max += 0.05 * x_range
Y_min -= 0.05 * y_range
Y_max += 0.05 * y_range
print(f"正射影像范围: X: {X_min:.2f}-{X_max:.2f}, Y: {Y_min:.2f}-{Y_max:.2f}")
# 设置分辨率
pixels_per_meter = 5 # 降低分辨率以加速处理
width = int((X_max - X_min) * pixels_per_meter)
height = int((Y_max - Y_min) * pixels_per_meter)
print(f"正射影像尺寸: {width}x{height} 像素 (分辨率: {pixels_per_meter} 像素/米)")
# 创建正射影像网格
ortho_x = np.linspace(X_min, X_max, width)
ortho_y = np.linspace(Y_max, Y_min, height) # Y轴反转
ortho_xx, ortho_yy = np.meshgrid(ortho_x, ortho_y)
ortho_zz = np.ones_like(ortho_xx) * Z0
# 将3D点展平
points_3d = np.dstack([ortho_xx, ortho_yy, ortho_zz]).reshape(-1, 3)
# 世界坐标→图像坐标
projected_points, _ = cv2.projectPoints(
points_3d,
R,
tvec,
K,
None
)
# 重塑为网格坐标
map_x = projected_points[:, 0, 0].reshape(ortho_xx.shape)
map_y = projected_points[:, 0, 1].reshape(ortho_yy.shape)
# 检查坐标范围有效性
valid_mask = (map_x >= 0) & (map_x < img_width) & \
(map_y >= 0) & (map_y < img_height)
valid_percentage = np.mean(valid_mask) * 100
print(f"有效投影点比例: {valid_percentage:.2f}%")
# 重映射生成正射影像
ortho_img = cv2.remap(
undistorted_img,
map_x.astype(np.float32),
map_y.astype(np.float32),
cv2.INTER_LINEAR,
borderMode=cv2.BORDER_CONSTANT,
borderValue=(0, 0, 0)
)
# 创建带透明通道的图像 (RGBA)
b, g, r = cv2.split(ortho_img)
alpha = np.uint8(valid_mask * 255)
ortho_img_rgba = cv2.merge([b, g, r, alpha])
# 保存结果
cv2.imwrite("river_ortho_corrected.png", ortho_img_rgba)
# 可视化比较
plt.figure(figsize=(15, 10))
# 原始图像
plt.subplot(131)
plt.imshow(cv2.cvtColor(raw_img, cv2.COLOR_BGR2RGB))
plt.title("原始图像")
plt.axis('off')
# 畸变校正后的图像
plt.subplot(132)
plt.imshow(cv2.cvtColor(undistorted_img, cv2.COLOR_BGR2RGB))
plt.title("畸变校正后图像")
plt.axis('off')
# 正射影像(带透明通道)
plt.subplot(133)
# 将RGBA转换为RGB显示,透明部分显示为白色背景
ortho_rgb = cv2.cvtColor(ortho_img_rgba, cv2.COLOR_BGRA2RGBA)
plt.imshow(ortho_rgb)
plt.title(f"正射影像 ({width}x{height})\n有效区域: {valid_percentage:.1f}%")
plt.axis('off')
plt.tight_layout()
plt.savefig("ortho_comparison.png", dpi=300, bbox_inches='tight')
plt.show()
# 保存调试图像(投影点可视化)
debug_img = undistorted_img.copy()
for i, pt in enumerate(projected_points.reshape(-1, 2)):
x, y = pt
if 0 <= x < img_width and 0 <= y < img_height:
color = (0, 255, 0) if i % 100 == 0 else (0, 0, 255) # 每隔100个点用绿色,其他用红色
cv2.circle(debug_img, (int(x), int(y)), 2, color, -1)
cv2.imwrite("projection_debug.jpg", debug_img)
# 计算重投影误差
reprojected, _ = cv2.projectPoints(world_points_3d, R, tvec, K, None)
errors = np.linalg.norm(reprojected - image_points, axis=2)
print(f"控制点重投影误差统计:")
print(f" 最小: {np.min(errors):.2f} 像素")
print(f" 最大: {np.max(errors):.2f} 像素")
print(f" 平均: {np.mean(errors):.2f} 像素")
print(f" 中位数: {np.median(errors):.2f} 像素")D:\Anaconda\python.exe D:\水位检测\正射校正\基于控制点的单应性变换.py
水面平均高程: 64.391 米
正射影像范围: X: 4262651.85-4262875.15, Y: 4262780.67-4262985.83
正射影像尺寸: 1116x1025 像素 (分辨率: 5 像素/米)
有效投影点比例: 0.00%
控制点重投影误差统计:
最小: 32.22 像素
最大: 2048.67 像素
平均: 737.62 像素
中位数: 632.97 像素