Open3D处理点云数据（三）点云单位方格内取平均高度信息处理

之所以写本文是因为：在过程中陷入了一个坑中：

背景应用是在点云单位方格内（比如1000mm*1000mm内）取平均高度传输给指定服务需求方。通过降低密度来提高输出效率。（虽然我们知道可以通过提高角度分辨率来降低采样率，减少点的个数，但是业务场景如此，不能较真吧。干吧）

结果我想检测我们自己的处理是否到位嘛，所以通过Open3d显示一下，陷入一个坑中：我只对高度信息取了平均值，但是呢，其实X轴和Y轴变换1000倍，这样导致的后果，图像是一个在Z轴上拉的特别长的图像，根本无法看出还是原来的物体。我还一直以为是不是我写错了呀。结果写了三种计算方式，才发现问题在第39行代码上：

z = point[2] / 250   # 不统一缩放一个比例的话，显示时看不出原图了。

 这个低级的错误，记录一下，也是初入点云世界的第一次坑。

全码如下：

import time
import numpy as np
import pandas as pd
import open3d as o3d


pcd = o3d.geometry.PointCloud()    # pcd类型的数据。

with open("0.csv", encoding="utf-8") as f:
    data = pd.read_csv(f, header=None).values.tolist()
    # print(data[:5])
    # print(data[10][2])
    # np_points = np.array(data)
    np_points = np.array(data)[:, 1:4]
    time1 = time.time()
    # np.savetxt('test0.csv', np_points, delimiter=',')

    # print(np_points)
    # print(np_points.shape)
    Mx = max(np_points[:, 0])
    # print("Max of x is ", Mx)
    mx = np_points[:, 0].min()
    # print(mx)
    My = max(np_points[:, 1])
    # print("Max of y is ", My)
    my = np_points[:, 1].min()
    # print(my)

    MX = int((Mx - mx) / 250) + 1
    MY = int((My - my) / 250) + 1

    # print("MX:", MX)
    # print("MY:", MY)
    sum_z = np.zeros((MX, MY))
    count_z = np.zeros((MX, MY))
    for point in np_points:
        x = int((point[0] - mx)/250)
        y = int((point[1]-my)/250)
        z = point[2] / 250   # 不统一缩放一个比例的话，显示时看不出原图了。
        # print(x, y, z)
        sum_z[x][y] += z
        count_z[x][y] += 1

    # print(np.where(count_z > 0)[0].shape)
    # print(np.where(count_z > 0)[1].shape)

    # print(np.where(count_z > 0)[0])  # x
    # print(np.where(count_z > 0)[1])  # y

    # print("点数:", count_z[np.where(count_z > 0)])
    # print("高度：", sum_z[np.where(count_z > 0)])
    X = np.where(count_z > 0)[0]  # x
    # print(type(X))
    Y = np.where(count_z > 0)[1]  # y
    # 第一种验证：
    # Z = np.zeros(len(X))
    # Z = [int(sum_z[np.where(count_z > 0)][i]/count_z[np.where(count_z > 0)][i]) for i in range(len(X))]  # z

    # 第二种验证。
    Z = np.zeros(len(X))
    for i in range(len(X)):
        Z[i] = sum_z[X[i]][Y[i]] / count_z[X[i]][Y[i]]
    Z = np.array(Z)
    print(type(Z))
    print(X.shape, Y.shape, Z.shape)

    # 第三种验证：
    # X = []
    # Y = []
    # Z = []
    # for x in range(MX):
    #     for y in range(MY):
    #         count = count_z[x][y]
    #         if count > 0:
    #             X.append(x)
    #             Y.append(y)
    #             Z.append(sum_z[x][y]/count)
    # print(len(X), len(Y), len(Z))
    # print(max(X), max(Y), min(X), min(Y))

    time2 = time.time()
    print("运行时间：", time2 - time1)
    m3d = np.array([X, Y, Z]).T
    # print(m3d.shape)
    np.savetxt('test1.csv', m3d, delimiter=',')

    pcd.points = o3d.utility.Vector3dVector(m3d)
    o3d.visualization.draw_geometries([pcd])

如果要输出Z轴的话，需要Z轴整体乘以250。

上面有一个坑：在计算网格数目的时候。如果出现整除现象就会多出一个方格。为此：解决方案：

第29行和30行处理如下：

MX = int((Mx - mx - 1) / 250) + 1
MY = int((My - my - 1) / 250) + 1

在计算除号之前，先减一，然后加一，原因是：假设5和2的情形，5减去1是4，除以2是2，加1是3成立。4减去1是3，3除以2是1，加1是2成立。

或者使用天花板除

MX = math.ceil((Mx - mx) / 250)
MY = math.ceil((My - my) / 250)