BimAnt

海量点云的可视化

admin -

数据可视化是一个大问题🌶️：通过使用视觉元素对信息进行图形表示，我们可以最好地呈现和理解数据中的趋势、异常值和模式。你猜对了：使用代表真实世界形状的 3D 点云数据集，这是强制性的 🙂。

本文处理和可视化无人机 3D 点云。你将在实时可视化和创建动画的同时学习特征提取、交互式和自动分割

但是，当从激光扫描仪或摄影测量等 3D 重建技术收集时，点云通常过于密集，无法进行经典渲染。在许多情况下，数据集将远远超过 1000 万大关，这使得它们对于 Matplotlib 等经典可视化库来说不切实际。

你可以注意到，左侧（Open3D）与右侧（PPTK）相比速度有多慢，PPTK 使用八叉树结构来加速可视化。Matplotlib 会更糟😅。

这意味着我们经常需要跳出 Python 脚本（因此使用 I/O 函数将数据写入文件）并在外部进行可视化，这可能会成为一个非常繁琐的过程🤯。我不会撒谎，这几乎就是我在论文第一年所做的，试图猜测特定算法的结果🥴。

直接在脚本中可视化这些点云不是很棒吗？更好的是，将视觉反馈连接到脚本？想象一下，现在 iPhone 12 Pro 配备了 LiDAR；你可以创建一个完整的在线应用程序！好消息，有一种方法可以实现这一点，而无需离开舒适的 Python 环境和 IDE。☕ 准备好了吗？

1、启动你的 Python 环境

在上一篇文章中，我们了解了如何使用 Anaconda 轻松设置环境以及如何使用 IDE Spyder 管理代码。如果你打算成为一名成熟的 Python 应用程序开发人员，我建议您继续以这种方式学习 😆。

如果你使用的是 Jupyter Notebook 或 Google Colab，脚本可能需要进行一些调整才能使可视化后端工作，但性能不稳定。如果你想继续使用这些 IDE，我建议你查看步骤 4 中给出的所选库的替代方案。

2、下载点云数据集

我在之前的教程中说明了使用摄影测量和来自 Open Topography 的航空 LiDAR 获得的 3D 数据集的点云处理和网格划分。我将跳过下文中介绍的有关 LiDAR I/O 的详细信息，直接跳转到使用高效的 .las 文件格式。

只是这一次，我们将使用空中无人机数据集。它是通过摄影测量获得的，让一架小型 DJI Phantom Pro 4 在我们的大学校园里飞行，收集一些图像并运行摄影测量重建，如此处所述。

以下链接提供的 3D 点云来自 DJI Phantom 4 飞行，随后进行了摄影测量重建过程

🤓 注意：对于本操作指南，你可以使用此存储库中的点云，我已经对其进行了过滤和翻译，以便你处于最佳条件。如果你想在不安装任何软件的情况下预先进行可视化和操作，可以查看 webGL 版本。

3、在脚本中加载点云

我们首先在脚本中导入必要的库（NumPy 和 LasPy），然后将 .las 文件加载到名为 point_cloud 的变量中。

import numpy as np
import laspy as lp
input_path="D:/CLOUD/POUX/ALL_DATA/"
dataname="2020_Drone_M"
point_cloud=lp.file.File(input_path+dataname+".las", mode="r")

很好，我们快准备好了！很棒的是，LasPy 库还为 point_cloud 变量提供了一个结构，我们可以使用简单的方法来获取例如 X、Y、Z、红色、蓝色和绿色字段。让我们这样做来将坐标与颜色分开，并将它们放入 NumPy 数组中：

points = np.vstack((point_cloud.x, point_cloud.y, point_cloud.z)).transpose()
colors = np.vstack((point_cloud.red, point_cloud.green, point_cloud.blue)).transpose()

🤓 注意：我们使用 NumPy 的垂直堆栈方法，我们必须将其转置以从 (n x 3) 变为 (3 x n) 点云矩阵。

4、最终预处理，可选

如果你的数据集太重，或者你想尝试子采样版本，我鼓励您可以查看这篇文章，其中提供了几种完成此类任务的方法。

为了方便起见，如果你的点云超过 1 亿个点，我们可以使用以下方法快速切分你的数据集：

factor=10
decimated_points_random = points[::factor]

🤓 注意：运行此操作将每 10 行保留 1 行，从而将原始点云的大小除以 10。

5、选择你的可视化策略

现在，让我们选择我们想要如何可视化我们的点云。我会在这里说实话：虽然可视化本身就很好，可以避免繁琐的 I/O 操作，但能够在 Python 中包含一些视觉交互和处理工具是一个很好的补充！因此，我推崇的解决方案是使用一个点云处理工具包，它可以实现这一点甚至更多。如果你想探索其他可能性，我仍然会为你提供替代方案⚖️。

解决方案 A（保留）：PPTK

PPTK 包有一个 3-d 点云查看器，它直接将 3 列 NumPy 数组作为输入，并且可以交互式地可视化 10 到 1 亿个点。它通过使用八叉树剔除视锥体之外的点并将远处点组近似为单个点，减少了每帧中需要渲染的点数。

模拟八叉树结构中的视锥体剔除

首先，你可以使用 Pip 管理器安装库：

pip install pptk

然后，你可以通过键入以下内容从点云中可视化之前创建的点变量：

import pptk
import numpy as np
v = pptk.viewer(points)
在 PPTK 中可视化的点云

启动时，查看器将输入点组织成八叉树。在操纵视点时，八叉树用于将远处的点组近似为单个点并剔除视锥体外的点，从而显著减少要渲染的点数。一旦视点不再发生变化，查看器将继续执行更耗时的详细点渲染。

你不觉得我们缺少一些颜色吗？让我们通过在控制台中输入以下内容来解决这个问题：

v.attributes(colors/65535)
PPTK 查看器中带有颜色信息的 3D 点云

🤓 注意：我们的颜色值是从 .las 文件中以 16 位编码的。我们需要 [0,1] 间隔内的值；因此，我们除以 65535。

这样好多了！但如果我们还想可视化其他属性怎么办？好吧，您只需将属性链接到路径，它就会即时更新。

可视化预先计算的几个属性的示例

💡 提示：不要最大化窗口大小以保持超过 30 FPS 的良好帧速率。目标是在拥有可读脚本的同时拥有最佳的执行运行时

你还可以参数化窗口以显示有关特定颜色渐变的每个属性，管理点大小，将背景设为黑色，而不显示网格和轴信息：

v.color_map('cool')
v.set(point_size=0.001,bg_color=[0,0,0,0],show_axis=0,show_grid=0)

替代方案 B：Open3D

对于任何想知道在 Python 中读取和显示点云的绝佳替代方案的人，我推荐 Open3D。你也可以使用 Pip 包管理器来安装必要的库：

pip install open3d

这将在你的机器上安装 Open3D，然后你将能够通过执行以下脚本读取和显示点云：

import open3d as o3d
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(points)
pcd.colors = o3d.utility.Vector3dVector(colors/65535)
pcd.normals = o3d.utility.Vector3dVector(normals)
o3d.visualization.draw_geometries([pcd])
Open3D 中可视化的 3D 点云。请注意如何很好地使用法线来增强几何图形的视觉效果

Open3D 实际上正在不断发展，你可以通过一些有趣的方式来显示点云，以填补最终的漏洞，例如创建体素结构：

voxel_grid = o3d.geometry.VoxelGrid.
create_from_point_cloud(pcd,voxel_size=0.40)
o3d.visualization.draw_geometries([voxel_grid])
使用 Python 在 open3D 中将点云转换为体素

点云的 3D 体素表示，其中每个体素代表一个 40 x 40 厘米的立方体。

🤓 注意：为什么 Open3d不是此时的选择？如果你处理的数据集少于 5000 万个点，那么这就是我推荐的。如果你需要在此阈值以上进行交互式可视化，我建议你对数据集进行采样以进行可视化，或者使用 PPTK，因为你为此目的创建了八叉树结构，因此可视化效率更高。

其他（Colab 友好型）替代方案：Pyntcloud 和 Pypotree

如果你希望能够简单且交互式地探索点云数据，无论使用哪种传感器来生成它或用例是什么，我建议你研究 Pyntcloud 或 PyPotree。这些将允许你在笔记本中可视化点云，但要注意性能！Pyntcloud 实际上依赖于 Matplotlib，而 PyPotree 需要 I/O 操作；因此，两者实际上都不是超级高效的。尽管如此，我想提到它们，因为对于小点云和 Google Colab 中的简单实验，你可以集成可视化。一些示例：

### PyntCloud ###
conda install pyntcloud -c conda-forge
from pyntcloud import PyntCloud
pointcloud = PyntCloud.from_file("example.ply")
pointcloud.plot()
### PyntCloud ###
pip install pypotree
import pypotree 
import numpy as np
xyz = np.random.random((100000,3))
cloudpath = pypotree.generate_cloud_for_display(xyz)
pypotree.display_cloud_colab(cloudpath)

6、与点云交互

返回 PPTK。要进行交互式选择，比如停车场上的汽车，我将移动相机顶视图（快捷键为 7），然后按住 Ctrl+LMB 拖动矩形选择进行选择。

选择要分割的点云数据

💡 提示：如果你对选择不满意，只需按 RMB 即可清除当前选择。是的，你可以进行多项选择 😀。

做出选择后，你可以返回 Python 控制台，然后获取分配的点标识符。

selection=v.get('selected')

这实际上会返回一个 1D 数组，如下所示：

选择是一个包含每个选定点的索引的数组

你实际上可以扩展该过程以一次选择多个元素（Ctrl+LMB），同时细化选择以删除特定点（Ctrl+Shift+LMB）。

3D 多点云分割。从点云中创建多个选择

在此之后，可以毫不费力地在保存选定点索引的选择变量上交互地应用一系列过程。

让我们复制一个场景，在该场景中，你可以自动在地面和非地面元素之间细化初始选择（汽车）。

7、实现自动分割

在包含完整点云（存储在变量 v 中）的查看器中，我进行以下选择：

selection=v.get('selected')
步骤 1：我们从初始 3D 点云中选择点

然后我计算每个点的法线。为此，我想说明使用 PPTK 的另一个关键要点：函数estimate_normals，可用于根据半径搜索或 k-最近邻获取每个点的法线。别担心，我将在另一本指南中深入说明这些概念，但现在，我将使用 6 个最近邻居来估计我的法线来运行它：

normals=pptk.estimate_normals(points[selection],k=6,r=np.inf)

💡 提示：请记住，selection 变量保存点的索引，即点云中的“线号”，从 0 开始。因此，如果我只想处理这个点子集，我会将其作为 points[selection] 传递。然后，我选择 k-NN 方法，只使用每个点的 6 个最近邻居，同时将 radius 参数设置为 np.inf，以确保我不使用它。我也可以使用这两个约束，或者如果我想进行纯半径搜索，则将 k 设置为 -1。

这基本上会返回以下内容：

每个点的法线样本

然后，我想过滤并返回法线不与 Z 轴共线的原始点的索引。我建议使用以下代码行：

idx_normals=np.where(abs(normals[...,2])<0.9)

🤓 注意：normals[...,2] 是 NumPy 的一种说法，表示我只处理 3 x n 点矩阵的第 3 列，保存法线的 Z 属性。它相当于 normals[:,2]。然后，我将绝对值作为比较点，因为我的法线没有方向（因此可以指向天空或地球中心），并且只会使用函数 np.where() 保留满足条件 <0.9 的那个。

为了直观地展示结果，我创建了一个新的查看器窗口对象：

viewer1=pptk.viewer(points[idx_normals],colors[idx_normals]/65535)
自动法线过滤器后的 3D 点云片段。看看车顶和整个汽车结构上的一些点是如何被丢弃的

如你所见，我们还过滤了汽车的一些点。这不好 🤨。因此，我们应将过滤与另一个过滤器相结合，以确保仅选择靠近地面的点作为法线过滤的主机：

idx_ground=np.where(points[...,2]>np.min(points[...,2]+0.3))
idx_wronglyfiltered=np.setdiff1d(idx_ground, idx_normals)
idx_retained=np.append(idx_normals, idx_wronglyfiltered)
viewer2=pptk.viewer(points[idx_retained],colors[idx_retained]/65535)
过滤垂直法线接近初始段最低 Z 值的点的 3D 点云

这很棒！现在，你可以探索这种强大的思维方式，并结合任何过滤（例如，在 RGB 上玩以摆脱剩余的草……）来创建一个完全交互式的分割应用程序。甚至更好的是，您可以将其与 3D 深度学习分类相结合！呵呵！但那是另一回事了 😉。

8、使用函数打包你的脚本

最后，我建议将你的脚本打包成函数，以便你可以直接将其部分重用为块。我们可以首先定义一个preparedata()，它将以任何.laspoint云作为输入，并对其进行格式化：

def preparedata():
    input_path="D:/CLOUD/OneDrive/ALL_DATA/GEODATA-ACADEMY/"
    dataname="2020_Drone_M_Features"
    point_cloud=lp.file.File(input_path+dataname+".las", mode="r")
    points = np.vstack((point_cloud.x, point_cloud.y, point_cloud.z) 
    ).transpose()
    colors = np.vstack((point_cloud.red, point_cloud.green,
    point_cloud.blue)).transpose()
    normals = np.vstack((point_cloud.normalx, point_cloud.normaly, 
    point_cloud.normalz)).transpose()
    return point_cloud,points,colors,normals

然后，我们编写一个显示函数pptkviz，返回查看器对象：

def pptkviz(points,colors):
    v = pptk.viewer(points)
    v.attributes(colors/65535)
    v.set(point_size=0.001,bg_color= [0,0,0,0],show_axis=0,
    show_grid=0)
    return v

另外，作为奖励，这里是 cameraSelector 函数，用于从打开的查看器中获取相机的当前参数：

def cameraSelector(v):
    camera=[]
    camera.append(v.get('eye'))
    camera.append(v.get('phi'))
    camera.append(v.get('theta'))
    camera.append(v.get('r'))
    return np.concatenate(camera).tolist()

我们定义了 computePCFeatures 函数来自动细化交互式分割：

def computePCFeatures(points, colors, knn=10, radius=np.inf):
    normals=pptk.estimate_normals(points,knn,radius)
    idx_ground=np.where(points[...,2]>np.min(points[...,2]+0.3))
    idx_normals=np.where(abs(normals[...,2])<0.9)
    idx_wronglyfiltered=np.setdiff1d(idx_ground, idx_normals)
    common_filtering=np.append(idx_normals, idx_wronglyfiltered)
    return points[common_filtering],colors[common_filtering]

Et voilà 😁，现在你只需启动包含上述函数的脚本，然后开始使用 computePCFeatures、cameraSelector 和更多功能与你的选择进行交互:

import numpy as np
import laspy as lp
import pptk
#Declare all your functions here
if __name__ == "__main__":
    point_cloud,points,colors,normals=preparedata()
    viewer1=pptkviz(points,colors,normals)

然后可以轻松调用脚本，然后使用控制台作为实验的基准。例如，我可以保存多个相机位置并创建动画：

cam1=cameraSelector(v)
#Change your viewpoint then -->
cam2=cameraSelector(v)
#Change your viewpoint then -->
cam3=cameraSelector(v)
#Change your viewpoint then -->
cam4=cameraSelector(v)
poses = []
poses.append(cam1)
poses.append(cam2)
poses.append(cam3)
poses.append(cam4)
v.play(poses, 2 * np.arange(4), repeat=True, interp='linear')
点云的 PPTK 中 4 个关键帧之间的线性插值

9、结束语

我们刚刚学习了如何导入、可视化和分割由 3000 多万个点组成的点云！做得好！

有趣的是，直接在 GPU 上执行的点云片段和单个点的交互式选择现在可以用于实时点云编辑和分割。但路径并未就此结束，未来的文章将深入探讨点云空间分析、文件格式、数据结构、分割 [2–4]、动画和深度学习 [1]。我们将特别研究如何管理下文定义的大点云数据。

原文链接：Visualise Massive point cloud in Python

汇智网翻译整理，转载请标明出处