colmap笔记：将模型和影像位姿转换到GPS坐标系下

colmap笔记：将模型和影像位姿转换到GPS坐标系下

1.背景和需求

在 COLMAP 中使用 GPS 信息作为先验，可以帮助提高重建的精度和效率，特别是在处理大范围或大尺度的场景时。GPS 信息通常被用来约束相机的位置，使重建过程更加准确。GPS 数据通常以 EXIF 元数据的形式存储在 JPEG 图像中，包含 GPS 坐标（经度、纬度、高度）。如果图像含有 EXIF 数据，COLMAP 会自动识别并读取这些信息。

但在稀疏重建阶段，COLMAP 主要关注于图像之间的相对位置和方向，以及场景中3D点的相对位置。默认情况下，COLMAP 不会使用 EXIF 中的 GPS 数据来定义重建模型的坐标系。因此，重建出的模型是与实际地理位置无关的局部坐标系。

那么现在想要将重建模型转换到地理坐标系（如基于 GPS 的全球坐标系），发现可以在重建完成后使用 model_aligner 工具将重建的模型缩放、平移和旋转，使其与真实世界的地理坐标系对齐。这要求提供包含 GPS 位置信息的参考图像列表。

官方文档的说明如下：

让我来一探究竟。

1.尝试

首先已经用colmap做完了稀疏重建，输出模型，然后进行model_aligner,设置好参数，将新的模型输出到gps文件夹。

colmap model_aligner \
    --input_path ./sparse/0 \
    --output_path ./sparse/gps \
    --database_path  ./database.db \
    --ref_is_gps 1 \
    --alignment_type ecef \
    --robust_alignment 1 \
    --robust_alignment_max_error 3.0

此时gps文件夹下的images.txt里记录的四元数和平移参数就是新的坐标系下的

根据以下代码，可以将其转换到GPS所在wgs84坐标系下：

import numpy as np
from scipy.spatial.transform import Rotation as R
import pyproj
from pyproj import Proj, transform

#0.016938035523210708 0.58455146147157355 -0.48870579156409283 0.64744060819180593 -129342.74756339534 3469822.8668770161 5343696.2087743161
# 四元数 (QW, QX, QY, QZ)
qw, qx, qy, qz = 0.016938035523210708, 0.58455146147157355, -0.48870579156409283, 0.64744060819180593   # 替换为实际值

# 平移向量 (TX, TY, TZ)
tx, ty, tz = -129342.74756339534, 3469822.8668770161, 5343696.2087743161

# 将四元数转换为旋转矩阵
rotation = R.from_quat([qx, qy, qz, qw])
rotation_matrix = rotation.as_matrix()

# 计算相机中心坐标(ECEF坐标系下)
camera_center = -np.dot(rotation_matrix.T, np.array([tx, ty, tz]))
print("投影中心坐标:（ecef坐标系下）", camera_center)

# 定义坐标转换
# 从ECEF转换到WGS84
ecef = Proj(proj='geocent', ellps='WGS84', datum='WGS84')
wgs84 = Proj(proj='latlong', ellps='WGS84', datum='WGS84')

# 执行转换
lon, lat, alt = transform(ecef, wgs84, camera_center[0],camera_center[1], camera_center[2], radians=False)

# 输出WGS84经纬度坐标
print("经度:", lon, "纬度:", lat, "高度:", alt)

# 从WGS84转换到UTM Zone 50N
utm_zone_50n = Proj(proj="utm", zone=50, ellps='WGS84', datum='WGS84', south=False)
utm_x, utm_y = transform(wgs84, utm_zone_50n, lon, lat)

# 输出UTM坐标
print("WGS84坐标系下：UTM Zone 50N X:", utm_x, "Y:", utm_y,"Z:",alt)

# 将四元数转换为旋转矩阵
rotation = R.from_quat([qx, qy, qz, qw])
rotation_matrix = rotation.as_matrix()

# 将旋转矩阵转换为欧拉角 (Omega, Phi, Kappa)
# 摄影测量中通常使用 ZYX 旋转顺序
omega, phi, kappa = rotation.as_euler('ZYX', degrees=True)
print("旋转矩阵:\n", rotation_matrix)

print("Omega:", omega, "Phi:", phi, "Kappa:", kappa)