基于python将激光点云.las文件投影为深度图/强度图/反强度图

使用python3.7版本 需安装laspy、numpy、matplotlib、opencv-python python使用laspy详情参见：https://laspy.readthedocs.io/en/latest/installation.html

pip install laspy pip install numpy pip install matplotlib pip install opencv-python 

部分代码参考：点云投影为深度图

# coding=utf-8  import laspy import cv2 import numpy as np import matplotlib.pyplot as plt  #相机内参 CAM_WID, CAM_HGT = 2875, 1875  # 深度图尺寸 CAM_FX, CAM_FY = 2777.75, 2777.75  # fx/fy CAM_CX, CAM_CY = 1437.5, 937.5  # cx/cy   EPS = 1.0e-16#阈值  # 加载.las点云数据 las = laspy.read('E:/peizhun/code/python/Toronto_Strip_05.las') pc=np.vstack((las.x,las.y,las.z,las.intensity)).transpose()   #plt.hist(las.intensity) #plt.hist(pc[:, 2])=plt.hist(las.z) #plt.hist(pc[:, 2]) #plt.title("Histogram of the Depth Dimension") #plt.show()  # 过滤镜头后面的点 valid = pc[:, 2] > EPS z = pc[valid, 2] i = pc[valid, 3]  ''' print(las.header.maxs[0]) print(las.header.maxs[1]) print(las.header.mins[0]) print(las.header.mins[1]) print(las.header.maxs[2]) print(las.header.mins[2]) '''  #坐标系平移-控制点云反投影范围，即反射到像素的坐标范围 #有些坐标系不需要转换，可以根据透视投影直接转换 x=pc[valid, 0]-(las.header.maxs[0]+las.header.mins[0])/2
y=pc[valid, 1]-(las.header.maxs[1]+las.header.mins[1])/2
z=z+1300


''' plt.title("Histogram of the X Dimension") plt.hist(x) plt.show() plt.title("Histogram of the Y Dimension") plt.hist(y) plt.show() plt.title("Histogram of the Z Dimension") plt.hist(z) plt.show() '''


# 点云反向映射到像素坐标位置（透视）
u = np.round(x * CAM_FX / z + CAM_CX).astype(int)
v = np.round(y * CAM_FY / z + CAM_CY).astype(int)

#显示转换后的点分布的直方图
#plt.title("Histogram of the u Dimension")
#plt.hist(u)
#plt.show()
#plt.title("Histogram of the v Dimension")
#plt.hist(v)
#plt.show()


# 滤除超出图像尺寸的无效像素
valid = np.bitwise_and(np.bitwise_and((u >= 0), (u < CAM_WID)),
                       np.bitwise_and((v >= 0), (v < CAM_HGT)))
u, v, z, i = u[valid], v[valid], z[valid], i[valid]


#显示滤波后的点分布的直方图
#plt.title("Histogram of the u Dimension")
#plt.hist(u)
#plt.show()
#plt.title("Histogram of the v Dimension")
#plt.hist(v)
#plt.show()
#plt.title("Histogram of the z Dimension")
#plt.hist(z)
#plt.show()
#plt.title("Histogram of the i Dimension")
#plt.hist(i)
#plt.show()


#将深度/强度/反强度归一化，并转为0——255灰度，0为黑，255为白
z_scaler =256* (z - min(z)) / (max(z) - min(z))
i_scaler =256* (i - min(i)) / (max(i) - min(i))
i_inv_scaler =256* (max(i) - i) / (max(i) - min(i))


''' #统计直方图的时候发现大部分都集中在一个小区域，猜测是个别点位的值与其它相差太大，导致归一化后依旧分布不均匀 #作者自己目视了一下大概没有点的区域，然后对归一化公式进行了更改 z_scaler =256* (z - min(z) - 40) / (max(z) - min(z) - 40) i_scaler =256* (i - min(i)) / (max(i) - 400 - min(i)) i_inv_scaler=256* (max(i) - 400 - i) / (max(i) - 400 - min(i)) plt.title("Histogram of the z_scaler Dimension") plt.hist(z_scaler) plt.show() plt.title("Histogram of the i_scaler Dimension") plt.hist(i_scaler) plt.show() plt.title("Histogram of the i_inv_scaler Dimension") plt.hist(i_inv_scaler) plt.show() '''

# 先生成一个空白图，数组行列数与图像长宽正好相反
# 按距离填充生成深度图，近距离覆盖远距离
img_z = np.full((CAM_HGT,CAM_WID), np.inf)
for ui, vi, zi in zip(u, v, z_scaler):
    img_z[vi, ui] = min(img_z[vi, ui], zi)  # 近距离像素屏蔽远距离像素

# 按强度填充生成强度图
img_i = np.full((CAM_HGT,CAM_WID), np.inf)
for ui, vi, ii in zip(u, v, i_scaler):
    img_i[vi, ui] = min(img_i[vi, ui], ii) 

# 按反强度填充生成强度图
img_i_inv = np.full((CAM_HGT,CAM_WID), np.inf)
for ui, vi, i_inv in zip(u, v, i_inv_scaler):
    img_i_inv[vi, ui] = min(img_i_inv[vi, ui], i_inv) 

	
	
# 小洞和“透射”消除
img_z_shift = np.array([img_z, \
                        np.roll(img_z, 1, axis=0), \
                        np.roll(img_z, -1, axis=0), \
                        np.roll(img_z, 1, axis=1), \
                        np.roll(img_z, -1, axis=1)])
img_z = np.min(img_z_shift, axis=0)

img_i_shift = np.array([img_i, \
                        np.roll(img_i, 1, axis=0), \
                        np.roll(img_i, -1, axis=0), \
                        np.roll(img_i, 1, axis=1), \
                        np.roll(img_i, -1, axis=1)])
img_i = np.min(img_i_shift, axis=0)

img_i_inv_shift = np.array([img_i_inv, \
                        np.roll(img_i_inv, 1, axis=0), \
                        np.roll(img_i_inv, -1, axis=0), \
                        np.roll(img_i_inv, 1, axis=1), \
                        np.roll(img_i_inv, -1, axis=1)])
img_i_inv = np.min(img_i_inv_shift, axis=0)


#镜像翻转，0表示绕x轴，1表示绕y轴，-1为绕x、y轴
#保存并显示最终结果图
img_depth=cv2.flip(img_z,0)
cv2.imwrite("depth.jpg", img_depth)
plt.imshow(img_depth)
plt.show()

img_intensity=cv2.flip(img_i,0)
cv2.imwrite("intensity.jpg", img_intensity)
plt.imshow(img_intensity)
plt.show()

img_intensity_inv=cv2.flip(img_i_inv,0)
cv2.imwrite("intensity_inv.jpg", img_intensity_inv)
plt.imshow(img_intensity_inv)
plt.show()

最终的灰度影像图

百度网盘：https://pan.baidu.com/s/1ra8rgO9nZtihlvO3OpFhfQ 提取码：zwn4