基于点云的散货体积计算

首先，你需要散货库存的点云。我将使用 IntelRealSense 捕获的散货库存的 .ply文件。然而，任何其他产生点云的成像技术都同样有效。

点击这里查看本教程的 Github 上的代码。

1、定位点云

让我们打开点云并检查它在笛卡尔 3D 空间中的位置。 Open3D 拥有现成的点云可视化工具。我们还将在可视化中添加轴，以便我们可以找到库存：

pcd = o3d.io.read_point_cloud("data/stockpile.ply")
axes = o3d.geometry.TriangleMesh.create_coordinate_frame()
o3d.visualization.draw_geometries([pcd, axes])

我们可以看到库存和地板。我们还可以看到它相对于轴是无方向的。为了计算体积，我们需要使地板平行于 XY 平面并将其平移到原点。为此，我们需要对地板进行分割。再次幸运的是，Open3D 为我们提供了帮助：

plane_model, inliers = pcd.segment_plane(distance_threshold=0.01,
                                         ransac_n=3,
                                         num_iterations=10000)
[a, b, c, d] = plane_model
plane_pcd = pcd.select_by_index(inliers)
plane_pcd.paint_uniform_color([1.0, 0, 0])
stockpile_pcd = pcd.select_by_index(inliers, invert=True)
stockpile_pcd.paint_uniform_color([0, 0, 1.0])
o3d.visualization.draw_geometries([plane_pcd, stockpile_pcd, axes])

现在我们已经分割了地板平面，我们可以使用一些数学方法通过旋转和平移点云将地板带到 XY 平面：

plane_pcd = plane_pcd.translate((0,0,d/c))
stockpile_pcd = stockpile_pcd.translate((0,0,d/c))
cos_theta = c / math.sqrt(a**2 + b**2 + c**2)
sin_theta = math.sqrt((a**2+b**2)/(a**2 + b**2 + c**2))
u_1 = b / math.sqrt(a**2 + b**2 )
u_2 = -a / math.sqrt(a**2 + b**2)
rotation_matrix = np.array([[cos_theta + u_1**2 * (1-cos_theta), u_1*u_2*(1-cos_theta), u_2*sin_theta],
                            [u_1*u_2*(1-cos_theta), cos_theta + u_2**2*(1- cos_theta), -u_1*sin_theta],
                            [-u_2*sin_theta, u_1*sin_theta, cos_theta]])
plane_pcd.rotate(rotation_matrix)
stockpile_pcd.rotate(rotation_matrix)
o3d.visualization.draw_geometries([plane_pcd, stockpile_pcd, axes])

我们可以丢弃地板，因为我们只关心库存：

o3d.visualization.draw_geometries([stockpile_pcd])

正如我们所看到的，我们有一些错误的分割点，特别是在 3D 图像的边界处。我们需要删除这些异常值。你知道这是怎么回事，Open3D 有一个专门用于此目的的函数：

cl, ind = stockpile_pcd.remove_statistical_outlier(nb_neighbors=30,
                                                    std_ratio=2.0)
stockpile_pcd = stockpile_pcd.select_by_index(ind)
o3d.visualization.draw_geometries([stockpile_pcd])

好的，我们终于对库存点云进行了分割并正确定向。我们现在可以跳到体积计算部分。耶！

2、计算体积

我们无法从点云计算体积，我们需要一个网格。为了从点云到网格，我们可以使用表面重建算法。表面重建是一个不适定（ill-posed）问题，这意味着没有完美的解决方案，算法基于启发式。然而，我们可以通过利用库存（或实际上任何地形）的一个特殊特征来避免这些算法，即表面受 XY 平面约束。没有 2 个点具有相同的 x 和 y，并且从 XY 平面上方观察该表面，我们将看到一个“平面”表面，也就是说，如果在任何点垂直穿过 XY 平面，只会穿过我们的表面一次。

我们将在 2D 空间中构建用于表面重建的三角化。我们将通过删除 Z 值将 PCD 投影到 XY 平面。我们将使用 Delaunay 算法获得 2D 点的三角测量。然后我们将保留 Delaunay 输出的三角形，但我们将使用 3D 顶点而不是 2D 顶点。这种方法有时称为“Delaunay 2.5D”

这听起来可能很复杂，但在代码中很容易看出：

downpdc = stockpile_pcd.voxel_down_sample(voxel_size=0.05)
xyz = np.asarray(downpdc.points)
xy_catalog = []
for point in xyz:
    xy_catalog.append([point[0], point[1]])
tri = Delaunay(np.array(xy_catalog)

tri 将为我们提供三角形的索引。在 2D 空间中，它看起来像：

使用带有 3D 点的三角形将为我们提供库存的表面：

surface = o3d.geometry.TriangleMesh()
surface.vertices = o3d.utility.Vector3dVector(xyz)
surface.triangles = o3d.utility.Vector3iVector(tri.simplices)
o3d.visualization.draw_geometries([surface], mesh_show_wireframe=True)

一旦有了表面，我们就准备好获得体积了。可以采用不同的方法。我将计算表面每个三角形相对于 XY 平面（高度 0）的体积。然后将所有三角形的体积相加。

Open3D 网格中的三角形定义顶点的索引，而不是顶点，因此我们需要平移三角形，以便它们的值实际上是 3D 点，而不是顶点列表的索引：

def get_triangles_vertices(triangles, vertices):
    triangles_vertices = []
    for triangle in triangles:
        new_triangles_vertices = [vertices[triangle[0]], vertices[triangle[1]], vertices[triangle[2]]]
        triangles_vertices.append(new_triangles_vertices)
    return np.array(triangles_vertices)

然后我们需要一个函数来计算每个三角形下的体积。根据布尔巴克数学家凯文·布朗的文章，我们可以将这样的函数定义为：

def volume_under_triangle(triangle):
    p1, p2, p3 = triangle
    x1, y1, z1 = p1
    x2, y2, z2 = p2
    x3, y3, z3 = p3
    return abs((z1+z2+z3)*(x1*y2-x2*y1+x2*y3-x3*y2+x3*y1-x1*y3)/6)

然后，获取库存的体积就像添加所有三角形的体积一样简单：

volume = reduce(lambda a, b:  a + volume_under_triangle(b), get_triangles_vertices(surface.triangles, surface.vertices), 0)
print(f"The volume of the stockpile is: {round(volume, 4)} m3")

希望你喜欢这篇文章😃

原文链接：Stockpile volume computation with Open3D

BimAnt翻译整理，转载请标明出处