【Lidar】基于Python格网法计算点云体积（eg.树木体积）

????????这两天一直不在状态，不是特别想分享文章，所以也没怎么更新。但是代码放在文件里始终不是它的归宿，只有被不断使用它才能进步，才能诠释它的意义。所以今天抽空给大家分享一下如何基于Python利用格网法计算点云的体积，我这里是做林业的点云，所以是按照树木体积编写的代码。

1 代码逻辑

? ? ? ? 逻辑部分其实很简单，就是将三维点云数据投影至二维平面，再通过格网将二维平面切割成无数个小块，统计位于该格网中点云的最高点和最低点差值，用这个高度差值乘以格网大小，即这个格网的体积，再将所有格网体积累加即该物体的体积。类似于积分学，将无数个格网体积近似等同于物体的体积。

2 完整代码

? ? ? ? 这里就不过多解释了，有注释！原理上面已经说明白了。

# -*- coding: utf-8 -*-
"""
@Time ： 2023/9/3 14:37
@Auth ： RS迷途小书童
@File ：Grid method for volume calculation.py
@IDE ：PyCharm
@Purpose：格网法计算点云体积并展示刨面图
"""
import open3d as o3d
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap


def grid_method(point_cloud, grid_size):
    points = np.asarray(point_cloud.points)
    # 将点云数据转换为数组
    points_x = points[:, 0]  # 获取x的数组
    points_y = points[:, 1]  # 获取y的数组
    cols_x = int(np.ceil((np.max(points_x) - np.min(points_x)) / grid_size))+1  # 计算栅格的列数
    cols_y = int(np.ceil((np.max(points_y) - np.min(points_y)) / grid_size))+1  # 计算栅格的行数
    start_x = np.min(points_x)  # 获取x的最小值作为格网的起始x
    start_y = np.min(points_y)  # 获取y的最小值作为格网的起始y
    result = 0  # 初始化体积为0
    array_show = np.zeros((cols_x, cols_y))  # 创建用来显示的数组
    for col_x in range(0, cols_x):  # 遍历格网的列
        end_x = np.min(points_x) + grid_size * (col_x+1)
        # 获取当前格网的范围
        for col_y in range(0, cols_y):  # 遍历格网的行
            end_y = np.min(points_y) + grid_size * (col_y+1)
            # 获取当前格网的范围
            mask = (points[:, 0] >= start_x) & (points[:, 0] <= end_x) & (points[:, 1] >= start_y) & (points[:, 1] <= end_y)
            # 通过格网范围筛选出落在当前网格内的点
            if np.sum(mask) > 0:
                # 如果当前网格内有点，则记录最大高度
                max_z = np.max(points[mask, 2])
                min_z = np.min(points[mask, 2])
                height = max_z - min_z  # 计算格网内的高度插差值
                # print("Max_z:",max_z)
                # print("Min_z:",min_z)
                # print("Height:",height)
                if height <= 0.0005:
                    # 记录高度差非常小的格网，原因：格网范围过小，最高点和最低点没有落在格网内
                    print("当前起始点：", start_x, start_y)
                    print("当前终止点：", end_x, end_y)
                    print("当前格网内的点云为：")
                    print(np.array([points[mask, 0], points[mask, 1], points[mask, 2]]))
                result += height * grid_size ** 2  # 将格网的体积累加
                array_show[col_x, col_y] = height
            start_y = end_y
        start_x = end_x
    return result, array_show


if __name__ == "__main__":
    pcd_path = "彭俊喜/1 - Cloud.pcd"  # 读取点云文件
    pcd = o3d.io.read_point_cloud(pcd_path)
    Grid_size = 0.1  # 设置栅格大小（根据实际情况调整）
    volume, array = grid_method(pcd, Grid_size)  # 计算树冠体积
    print("树冠体积:", volume)
    fig, ax = plt.subplots()  # 创建图形和轴对象
    cmap = LinearSegmentedColormap.from_list('white_to_red', ['white', (0.8, 0.5, 0.2), 'yellow', 'green'])
    im = ax.imshow(array, cmap=cmap)  # 绘制图像'coolwarm'
    cbar = fig.colorbar(im, ax=ax, orientation='vertical')  # 添加颜色拉伸的图例
    cbar.set_label('Colorbar Label')
    plt.savefig(r"G:\Anaconda\ProjectYOLO\yolov5-7.0\Point Cloud/1.png")  # 将图形窗口中的内容保存到指定的路径
    plt.show()  # 显示图形
    

3 效果展示

4 总结

? ? ? ? 对于点云数据体积的计算，格网法是个不错的选择。但是要注意格网大小的选择，如果过大可能导致体积拟合大于实际值过多。如果过小，可能很多格网都在点云的缝隙之间，这就会导致拟合的体积过小。

? ? ? ? 大家如果对点云数据的处理感兴趣可以随时留言交流。如果觉得博主的文章对你有帮助，可以给我点个赞，加个关注，我后面会更新更多点云数据处理的教程。