如何使用 python 从 2d 坐标中获取 3d 空间坐标？

Posted 2023-03-23

【中文标题】如何使用 python 从 2d 坐标中获取 3d 空间坐标？

【英文标题】：how to get 3d space coordinates from 2d coordinates using python?

【发布时间】：2022-01-07 18:28:23

【问题描述】：

我有一个 2d 正射镶嵌图像及其点云。如果我在图像上的对象周围绘制一个多边形，那么在输出中它应该从 3d 点云返回该对象的坐标。一旦我从 3d 云中获得对象的坐标，我就想计算该对象的体积。 任何人都可以帮助我吗？

【问题讨论】：

假设点的坐标与马赛克（加Z）在同一坐标系中可用，您需要找到XY投影落在多边形内的所有点。为了提高效率，您可以使用加速技术，例如边界框层次结构。计算体积具有挑战性，因为您需要重建表面。您可以尝试对 2D 点进行三角剖分并对相应 3D 棱镜的体积求和。

【参考方案1】：

从点云中检索正射镶嵌图像的对应点并非易事，因为在数字图像中，您所拥有的是具有由行和列索引表示的“位置”的像素，而在点云数据中, 每个点都有一个 (X, Y, Z) 值对，是在现实世界坐标系中测量的坐标。

因此，您需要统一图像和点云之间的坐标。对于正射镶嵌图像（假设它是由一些常用软件生成的），它应该是一个tiff文件，其元数据中定义了转换规则。该变换规则可以将图像像素坐标i, j变换为世界坐标X, Y。

这些名为“Geotransform”的转换规则通常由转换矩阵定义，并且可以使用 GDAL 等库从元数据中检索： （确保坐标系地理变换目标是您的点云数据的确切系统，否则，需要进行点云注册）

from osgeo import GDAL

def get_geotransform(path: str):
    img = gdal.Open(path, 0)    # orthomosaic image path
    return img.GetGeoTransform()

上述方法返回一个(6, )大小的数组，包含六个变换参数。

GDAL 官方文档有a good article 说明如何使用这些参数。如果您想转换多个点，我还实现了矢量化方法（使用风险自负）：

geotransform = get_geotransform(sample_img)

trans_matrix = np.array([[geotransform[2], geotransform[5]], [geotransform[1], geotransform[4]]])

XY = xy @ trans_matrix + np.array([geotransform[0], geotransform[3]])

这里，xy 是一个 (N, 2) ndarray，每行包含 (row, col) 索引。

所以我们从上面的过程中得到了世界坐标X，Y，现在我们需要在点云数据中搜索离X，Y最近的点，因为不能保证转换后，你能找到准确的点云数据中的相同点条目。

您可以通过多种方式实现此目标，即使是蛮力搜索也可以完成这项工作，尽管考虑到整体效率，我建议使用KD-Tree 数据结构。假设您的点云数据是一个(N, 3) 大小的ndarray，您可以使用scikit-learn 拟合一个KDTree：

from sklearn.neighbors import NearestNeighbors

kd_tree = NearestNeighbors(n_neighbors=1, n_jobs=-1)
kd_tree.fit(cloud[:, 0: 2])

并获取最近的点云入口索引：

nearest_point_index = kd_tree.kneighbors(XY, return_distance=False)

使用此索引访问点云数据应该会为您提供所选图像点的最近似坐标。

总结，需要两个步骤：

将图像点索引转换为世界坐标系 通过点云执行最小距离搜索以找到对应的坐标

另外，以上代码纯属演示，请根据自己的使用场景进行调整。例如，如果您只需要查询少量点，而您的点云数据非常大，KD-Tree 方法可能整体效率较低，需要更多细节才能更具体地回答这个问题。

【讨论】：

我宁愿使用矩形查询而不是最近邻，使用多边形的轴对齐边界框。然后使用完整的多边形测试过滤候选点。

虽然我自己没有做过测试，但@YvesDaoust 的点查询解决方案确实应该更适合这个特定任务。