点云二维高程属性格网投影python

1 原理

关于点云投影其实三维点云投影成二维的栅格图像过程，传统的图像在python中通过固定行列的数组表示，数组中的每个数值为格网属性。其投影过程主要分为以下几步：

1. 获取点云数据边界(Xmax,Ymax,Xmin,Ymin)
2. 计算投影格网行列总数(Row,Col)
   需要注意点云坐标系，当坐标系的水平轴为X，则投影的列于X对应；同理，坐标系的纵轴为Y，则投影的行于Y对应。
   $$\begin{gathered} Row=ceil(\frac{Y_{max}-Y_{min}}{step}) \\ Col=ceil(\frac{X_{max}-X_{min}}{step}) \end{gathered}$$
   step: 格网间隔
3. 计算当前点$P_i(x_i,y_i)$所在行列索引$(row{s}_i,col{s}_i)$
   $$\begin{gathered} row{s}_i=ceil(\frac{y_i-Ymin}{step}) \\ col{s}_i=ceil(\frac{x_i-Xmin}{step}) \end{gathered}$$
4. 为了方便准确获取当前点的格网位置，可以建立数组，每一行为[x,y,z,行索引,列索引]，当然也可通过格网的行列索引为格网编号number，不同格网之间编号不同。上一步计算出每个点的行列索引，为其编号，能够快速定位格网序列，编号后建立数组，每一行为[x,y,z,编号]。
   $$numbe{r}_i=(row{s}_i-1)\ast Col+col{s}_i$$

2 代码

将格网内点的平均高程作为格网属性值，生成均值高程属性图像。

def MeanAltitudeProjection(pcd,step):
    points=np.asarray(pcd.points)
    total=points.shape[0]

    # 获取点云数据边界
    x_max = max(points[:, 0])
    x_min = min(points[:, 0])
    y_max = max(points[:, 1])
    y_min = min(points[:, 1])

    # 计算格网行列数
    Row = int(np.ceil((y_max - y_min) / step)) + 1
    Col = int(np.ceil((x_max - x_min) / step)) + 1
    img_meanH=np.zeros((Row,Col))

    # 计算点所在行列
    number = np.zeros(total)  # 该点对应的格网编号
    for i in range(total):
        rows = (int(np.ceil((points[i, 1] - y_min) / step)))
        cols = (int(np.ceil((points[i, 0] - x_min) / step)))
        number[i] = (rows - 1) * Col + cols

    number_sort=sorted(number)
    number_argSort=np.argsort(number)
    begin=0
    for i in range(total-1):
        if number_sort[i]==number_sort[i+1]:  # 若前后点编号相同，视为同一格网
            continue
        else:
            grids_points=points[number_argSort[begin:i+1]]
            rows = int(np.ceil((grids_points[0, 1] - y_min) / step))
            cols = int(np.ceil((grids_points[0, 0] - x_min) / step))
            img_meanH[rows,cols]=np.mean(grids_points[:,2])
        begin=i+1
    return img_meanH

格网投影参考：马新江《一种基于几何特征的高速公路路面点云提取方法》