【问题标题】：Algorithm or command line tool to decimate point cloud of terrain points?抽取地形点的点云的算法或命令行工具？
【发布时间】：2017-09-04 00:12:48
【问题描述】：

我需要获取比所需的激光雷达测量点列表（经度、纬度和海拔）更大（更密集）的地形定义，并根据二维网格对其进行抽取。我们的想法是最终使用经度、纬度 (x,y) 值基于 NxN（即 1 米 x 1 米）维度网格的点，从而消除超出需要的点。目标是确定抽取后网格中每个点的高程是多少，而不是将高程用作抽取规则本身的一部分。

一个实际的或精确结构化的网格不是必需的，也不是这里的目标，我只使用网格术语来最好地近似我设想的点云的剩余部分，然后以我们总是有一个点的方式减少它在一定半径内（即1米）。可能有比网格更好的术语。

如果我可以从抽取算法开始或使用可能已经存在的可以在 Ubuntu 上运行的项目中执行此操作的命令行工具，我想自己用脚本或编程语言编写代码/编写脚本从我们的应用程序作为系统调用调用。该方法不应该要求使用基于 GUI 类型的软件或工具来解决这个问题。它需要成为一组自动化步骤的一部分。

数据当前存在于制表符分隔值文件中，但如果使用数据库/sql 查询驱动算法会更好/更快，我可以将数据加载到 sqlite 数据库文件中。理想的脚本语言是 ruby 或 python，但实际上可以是任何一种，如果已经存在 C/C++/C# 库，那么我们可以根据需要包装它们。

想法？

更新澄清这个抽取列表的结果的用途：给定用户的位置（通过纬度和经度知道），列表中最近的点是什么，反过来它的海拔是多少？我们现在当然可以这样做，但是我们拥有的数据比需要的多，所以我们只想放宽数据的密度，这样如果我们可以在允许的距离（即 1 米）内找到最近的点，如果能够使用抽取列表与完整列表。列表中的纬度、经度值是十进制 GPS（即 38.68616190027656、-121.11013105991036）

【问题讨论】：

可以分享数据吗？？？为了正确理解您，在构建 Grid 之后，您要为每个网格单元保留哪个点？？
@DaviddelaIglesia，我可以尝试创建一些示例数据。精确的网格不是目标，只是列表/云中最接近实际网格可能基于纬度和经度的点。我刚刚注意到您的回答，并会在一分钟内对其进行审核。
我会更新很多答案以满足您的更新
据我所知，您的答案可能是我已经需要的，只是我还没有时间测试它。可能要到星期一才能做到。我更新了我的帖子，以更清楚地说明我这样做的原因。
我认为我的新答案应该是寻找新用户高度的更有效方式

标签： algorithm point-cloud-library point-clouds lidar

【解决方案1】：

第 1 部分：抽取版本

加载数据

从表格文件中加载数据（根据您使用的分隔符更改sep）：

# installed as dependency
import pandas as pd  
# https://github.com/daavoo/pyntcloud
from pyntcloud import PyntCloud  

dense = PyntCloud(pd.read_csv("example.tsv",
                     sep='\t',
                     names=["x","y","z"]))

这是我创建的示例的外观：

构建体素网格

假设文件中的纬度和经度以米为单位，您可以生成如下网格：

grid_id = dense.add_structure("voxelgrid",
                           sizes=[1, 1,None],
                           bb_cuboid=False)
voxelgrid = dense.voxelgrids[grid_id]

这个体素网格在x（纬度）和y（经度）维度上的大小为1。

构建抽取版本

decimated = dense.get_sample("voxelgrid_centroids", voxelgrid=grid_id)

decimated 是一个 numpy (N,3) 数组。您可以将其存储以供以后在 SQL 数据库等中使用。

第 2 部分：查询

选项 A：查询体素网格

获取每个网格单元的平均高度

您可以知道为网格中的每个单元格获取一个具有平均z（高度）值的向量：

z_mean = voxelgrid.get_feature_vector(mode="z_mean")

用用户的位置查询网格：

users_location = np.random.rand(100000, 2)

添加一列零，因为查询需要 3D（这不会影响结果）：

users_location = np.c_[ users_location, np.zeros(users_location.shape[0]) ]

进入每个用户所在的单元格：

users_cell = voxelgrid.query(users_location)

最后，得到每个用户对应的海拔高度：

users_altitude = z_mean[users_cell]

选项 B：使用抽取版本进行查询

构建一个抽取的 KDTree：

从 scipy.spatial 导入 cKDTree kdt = cKDTree(抽取)
使用用户位置查询 KDTree：

users_location = np.random.rand(100000, 2) users_location = np.c_[ users_location, np.zeros(users_location.shape[0])

距离，索引 = kdt.query(user_locations, k=1, n_jobs=-1)

额外，你可以用pickle保存和加载voxelgrid：

pickle.dump(voxelgrid, open("voxelgrid.pkl", "wb"))

voxelgrid = pickle.load(open("voxelgrid.pkl", "rb"))

【讨论】：

感谢您的回答。乍一看，至少从方法名称的概念上看，编码似乎与我正在寻找的东西保持一致。我需要更仔细地测试和审查以提供更多反馈。我可以说您嵌入的视频对我来说很难看到与云中抽取的网格状数据集有任何相似之处，但这可能只是视频的质量/大小。如果列表中的纬度、经度值是十进制 GPS 格式 -（即 38.68616190027656、-121.11013105991036）会改变它的工作方式吗？
@Streamline 我认为我有更好的方法来做到这一点。但是我还没有推送到github。如果您不介意等待，我会尽快推送。
这个更新的答案是否允许我在某个时间点拆分生成抽取列表的任务，然后再尝试用户位置查找？（请参阅我对主要问题的评论）。这部分任务是减少和消除不需要的数据。
这都是关于体素网格的。要抽取点云，您可以将体素网格质心作为新的抽取版本。要查询用户位置，您可以使用我在答案中显示的方法，也可以对抽取的版本使用任何邻居搜索。
我认为我需要 python3 并且在查看 pyntcloud.readthedocs.io/en/latest/installation.html 后确定了依赖项。我已经第一次尝试了，在 python 控制台的上述步骤 1 中的这一行 voxelgrid = dense.voxelgrids[grid_id]，我得到了 AttributeError: 'PyntCloud' object has no attribute 'voxelgrids'。如果我输入dense，我确实会看到1 voxelgrids。建议？自 4 月以来图书馆是否发生了变化，我应该修改这些步骤？

【解决方案2】：

如果您有一个点云作为文本文件 (.xyz)，一个简单快速的解决方案是使用 shuf 从文件中随机抽取样本。

xyz 文件中的 1000 万个点等于 1000 万行文本。你可以运行：

shuf input.xyz -l 5000000 -o out.xyz

您已将文件缩减为原始大小的一半。

【讨论】：