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

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

1. 透视投影原理简述

针孔相机成像模型：

简单解释来源：

2. 前期准备

使用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

3. 代码

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

# 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()

4. 结果

最终的灰度影像图

5. 点云数据

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