基于点云的散货体积计算

首先,你需要散货库存的点云。 我将使用 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翻译整理,转载请标明出处