首先，你需要散货库存的点云。 我将使用 IntelRealSense 捕获的散货库存的 .ply文件。 然而，任何其他产生点云的成像技术都同样有效。
点击这里查看本教程的 Github 上的代码。
NSDT工具推荐： Three.js AI纹理开发包 - YOLO合成数据生成器 - GLTF/GLB在线编辑 - 3D模型格式在线转换 - 可编程3D场景编辑器 - REVIT导出3D模型插件 - 3D模型语义搜索引擎 - Three.js虚拟轴心开发包 - 3D模型在线减面 - STL模型在线切割
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")希望你喜欢这篇文章😃
原文链接：基于点云的散货体积计算 - BimAnt
