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

• 1 原理
• 2 代码

1 原理

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

1. 获取点云数据边界(Xmax,Ymax,Xmin,Ymin)
2. 计算投影格网行列总数(Row,Col)
   需要注意点云坐标系，当坐标系的水平轴为X，则投影的列于X对应；同理，坐标系的纵轴为Y，则投影的行于Y对应。
   R o w = c e i l ( Y m a x − Y m i n s t e p ) C o l = c e i l ( X m a x − X m i n s t e p ) Row=ceil(frac{Y_{max}-Y_{min}}{step}) Col=ceil(frac{X_{max}-X_{min}}{step}) Row=ceil(stepYmax​−Ymin​​)Col=ceil(stepXmax​−Xmin​​)
   step: 格网间隔
3. 计算当前点 P i ( x i , y i ) P_{i}(x_{i},y_{i}) Pi​(xi​,yi​)所在行列索引 ( r o w s i , c o l s i ) (rows_{i},cols_{i}) (rowsi​,colsi​)
   r o w s i = c e i l ( y i − Y m i n s t e p ) c o l s i = c e i l ( x i − X m i n s t e p ) rows_{i}=ceil(frac{y_{i} - Ymin} {step}) cols_{i}=ceil(frac{x_{i} - Xmin}{step}) rowsi​=ceil(stepyi​−Ymin​)colsi​=ceil(stepxi​−Xmin​)
4. 为了方便准确获取当前点的格网位置，可以建立数组，每一行为[x,y,z,行索引,列索引]，当然也可通过格网的行列索引为格网编号number，不同格网之间编号不同。上一步计算出每个点的行列索引，为其编号，能够快速定位格网序列，编号后建立数组，每一行为[x,y,z,编号]。
   n u m b e r i = ( r o w s i − 1 ) ∗ C o l + c o l s i number_{i}=(rows_{i} - 1) * Col + cols_{i} numberi​=(rowsi​−1)∗Col+colsi​

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

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