KITTI数据集实战：从零开始处理3D点云数据（附Python可视化代码）
·
KITTI数据集实战：从零开始处理3D点云数据（附Python可视化代码）
第一次拿到KITTI数据集的.bin文件时，我盯着那个几十兆的二进制文件愣了半天。身边有经验的同事告诉我，这里面藏着几十万个三维空间点，是自动驾驶汽车“看到”的真实世界。但怎么把这些冰冷的二进制数据变成屏幕上能旋转、能缩放、能直观感受的3D场景？这成了我入门3D视觉的第一个拦路虎。如果你也正站在这个起点，面对一堆陌生的文件格式和坐标系转换公式感到无从下手，那么这篇文章就是为你准备的。我们将抛开复杂的理论推导，直接从最实际的代码操作入手，一步步拆解KITTI点云数据的处理流程，并对比两种主流的可视化工具，让你在半小时内就能看到自己的第一个3D点云场景。
1. 理解KITTI：不只是数据，更是多传感器融合的典范
很多人把KITTI简单理解为一个“3D目标检测数据集”，这其实大大低估了它的价值。它本质上是一个多传感器同步采集的时空对齐数据宝库。2012年由德国卡尔斯鲁厄理工学院和丰田研究院联合发布时，就旨在解决车载环境下计算机视觉算法的标准化评测问题。我最初接触时，最震撼的是它的“真实性”——所有数据都来自真实的城市、乡村和高速公路驾驶场景，包含了各种天气、光照条件和复杂的交通参与者交互。
KITTI的核心在于其精密的传感器套件和严格的时空同步。想象一下，一辆装备了多种传感器的数据采集车在路上行驶：
- 激光雷达（LiDAR）：Velodyne HDL-64E，每秒旋转10次，产生约130万个点，精确测量周围120米范围内的物体距离。
- 相机系统：两对立体相机（一灰一彩），以54厘米的基线距离安装，用于获取高分辨率的图像信息。
- 定位系统：高精度的GPS/IMU组合导航单元，提供车辆自身的位姿信息。
这些传感器并非独立工作。它们在硬件层面被同步到同一个时钟源，确保每一帧点云、每一张图像都对应着同一个瞬间的世界状态。这种设计使得KITTI成为了研究传感器融合算法的绝佳平台——你可以用图像丰富的纹理信息来补充点云的几何信息，也可以用点云精确的距离信息来辅助图像中的目标定位。
提示：虽然KITTI官方提供了完整的数据包，但对于初学者或实验性项目，我强烈建议先从社区维护的“精简版”或“子集”开始。直接处理近200GB的原始数据，对硬件和耐心都是不小的考验。
在实际研究或工程中，我们通常只关注其中的一个子集。对于3D点云处理入门，最常用的是“3D Object Detection”任务的数据，它包含了：
- 7481个训练样本（带标注）
- 7518个测试样本（标注不可见，用于提交在线评测）
每个样本都包含：
- 一张彩色图像（来自左彩色相机）
- 对应的点云数据（
.bin文件） - 相机标定参数（
.txt文件） - 3D边界框标注（仅训练集有，
.txt文件）
文件夹的组织结构清晰，但初次接触时仍需注意：
kitti/
├── training/
│ ├── calib/ # 标定文件，每个样本一个.txt
│ ├── image_2/ # 左彩色相机图像，.png格式
│ ├── label_2/ # 3D标注，每个样本一个.txt
│ └── velodyne/ # 点云数据，每个样本一个.bin
└── testing/ # 测试集，结构类似但没有label_2
这种结构虽然标准，但在实际编程中，我们往往需要根据自己项目的需求重新组织或创建数据加载管道。接下来，我们就从最核心的点云数据开始。
2. 点云数据解码：从二进制文件到三维坐标
KITTI的点云文件以.bin后缀存储，这是一种紧凑的二进制格式。每个点用4个float32数字表示，分别对应(x, y, z, reflectance)。前三个是点在激光雷达坐标系下的三维坐标，单位是米；第四个是反射强度，表示激光脉冲被物体表面反射回来的强度，这个值对于区分不同材质（如金属、玻璃、植被）很有帮助。
2.1 基础读取与理解坐标系
用Python读取这些数据非常简单，但理解数据背后的坐标系是关键。激光雷达的坐标系定义是：
- x轴：指向车辆前进方向
- y轴：指向车辆左侧
- z轴：指向上方
这符合右手坐标系规则：伸出右手，食指指向x轴正方向，中指指向y轴正方向，那么拇指的方向就是z轴正方向。
import numpy as np
import struct
def read_bin_file(file_path):
"""
读取KITTI点云.bin文件
参数:
file_path: .bin文件的路径
返回:
points: numpy数组，形状为(N, 4)，每行是[x, y, z, reflectance]
"""
# 一次性读取整个文件
with open(file_path, 'rb') as f:
data = f.read()
# 每个点4个float32，共16字节
point_size = 4 * 4 # 4个float32，每个4字节
num_points = len(data) // point_size
# 使用struct模块解析二进制数据
# 'f'表示float32，'<'表示小端字节序
points = np.zeros((num_points, 4), dtype=np.float32)
for i in range(num_points):
start = i * point_size
end = start + point_size
# 解析4个float32
x, y, z, reflectance = struct.unpack('<ffff', data[start:end])
points[i] = [x, y, z, reflectance]
return points
# 更高效的numpy直接读取方式
def read_bin_file_fast(file_path):
"""
使用numpy的fromfile函数快速读取，效率更高
"""
# 直接读取为float32数组，然后重塑形状
points = np.fromfile(file_path, dtype=np.float32).reshape(-1, 4)
return points
# 使用示例
if __name__ == "__main__":
# 替换为你的.bin文件路径
sample_file = "./data/kitti/training/velodyne/000000.bin"
# 方法1：使用自定义解析（更可控）
points_custom = read_bin_file(sample_file)
print(f"自定义方法读取到 {points_custom.shape[0]} 个点")
# 方法2：使用numpy快速读取（推荐）
points_fast = read_bin_file_fast(sample_file)
print(f"快速方法读取到 {points_fast.shape[0]} 个点")
# 查看前5个点
print("前5个点的坐标和反射强度:")
print(points_fast[:5])
# 基本统计信息
print(f"\n数据统计:")
print(f"x范围: [{points_fast[:, 0].min():.2f}, {points_fast[:, 0].max():.2f}]")
print(f"y范围: [{points_fast[:, 1].min():.2f}, {points_fast[:, 1].max():.2f}]")
print(f"z范围: [{points_fast[:, 2].min():.2f}, {points_fast[:, 2].max():.2f}]")
print(f"反射强度范围: [{points_fast[:, 3].min():.2f}, {points_fast[:, 3].max():.2f}]")
运行这段代码，你会看到类似这样的输出：
快速方法读取到 115384 个点
前5个点的坐标和反射强度:
[[ 6.27570009 2.26649976 -1.43900001 0. ]
[ 6.28249979 2.28530002 -1.42900002 0. ]
[ 6.28929996 2.30410004 -1.41900003 0. ]
[ 6.29610014 2.32290006 -1.40900004 0. ]
[ 6.30289984 2.34169984 -1.39900005 0. ]]
数据统计:
x范围: [-79.99, 79.99]
y范围: [-79.99, 79.99]
z范围: [-4.99, 3.49]
反射强度范围: [0.00, 1.00]
注意z轴的数值范围：大部分点在-5米到3.5米之间。这是因为激光雷达安装在车顶，扫描线有一定的俯仰角，能够看到地面（负z值）和上方建筑物（正z值）。
2.2 数据预处理：过滤与降采样
原始点云数据通常包含十几万甚至几十万个点，直接可视化或处理计算量很大。在实际应用中，我们经常需要进行预处理。以下是一些常见的操作：
def preprocess_point_cloud(points, ground_height_threshold=-1.5, max_range=70.0):
"""
点云预处理：过滤地面点、限制范围、降采样
参数:
points: 原始点云，形状(N, 4)
ground_height_threshold: 地面高度阈值，低于此值的点被认为是地面
max_range: 最大距离，超出此距离的点被过滤
返回:
filtered_points: 处理后的点云
"""
# 1. 移除过远的点（基于距离）
distances = np.sqrt(points[:, 0]**2 + points[:, 1]**2)
range_mask = distances < max_range
# 2. 移除地面点（简单高度阈值法）
height_mask = points[:, 2] > ground_height_threshold
# 3. 组合掩码
mask = range_mask & height_mask
filtered_points = points[mask]
# 4. 降采样（随机采样）
if len(filtered_points) > 50000:
# 如果点数太多，随机采样到50000个点
indices = np.random.choice(len(filtered_points), 50000, replace=False)
filtered_points = filtered_points[indices]
return filtered_points
def calculate_point_intensity_stats(points):
"""
计算反射强度的统计特征，可用于点云分割
"""
reflectance = points[:, 3]
stats = {
'mean': np.mean(reflectance),
'std': np.std(reflectance),
'min': np.min(reflectance),
'max': np.max(reflectance),
'percentile_25': np.percentile(reflectance, 25),
'percentile_50': np.percentile(reflectance, 50), # 中位数
'percentile_75': np.percentile(reflectance, 75),
}
# 反射强度直方图（用于分析材质分布）
hist, bins = np.histogram(reflectance, bins=20, range=(0, 1))
return stats, hist, bins
# 预处理示例
if __name__ == "__main__":
# 读取原始数据
points = read_bin_file_fast("./data/kitti/training/velodyne/000000.bin")
# 预处理
processed_points = preprocess_point_cloud(points)
print(f"原始点数: {len(points)}")
print(f"处理后点数: {len(processed_points)}")
print(f"过滤比例: {(1 - len(processed_points)/len(points))*100:.1f}%")
# 反射强度分析
stats, hist, bins = calculate_point_intensity_stats(processed_points)
print(f"\n反射强度统计:")
for key, value in stats.items():
print(f"{key}: {value:.4f}")
预处理后的点云数据量会显著减少，但保留了大部分有意义的信息。在实际的自动驾驶系统中，这种预处理是必不可少的，它可以：
- 减少计算负担
- 移除噪声和无关点（如过远的地面点）
- 为后续的目标检测、分割算法提供更干净的数据
3. 坐标系转换：连接激光雷达与相机视角
KITTI数据集最强大的特性之一就是多传感器数据的对齐。这意味着我们可以将激光雷达点云投影到相机图像上，或者将图像中的像素信息关联到3D点云中。要实现这一点，我们需要理解并应用坐标系转换。
3.1 理解标定文件
每个样本都有一个对应的标定文件（在calib文件夹中），这个文件包含了所有必要的转换矩阵。让我们仔细解析一个典型的标定文件：
P0: 7.215377e+02 0.000000e+00 6.095593e+02 0.000000e+00 0.000000e+00 7.215377e+02 1.728540e+02 0.000000e+00 0.000000e+00 0.000000e+00 1.000000e+00 0.000000e+00
P1: 7.215377e+02 0.000000e+00 6.095593e+02 -3.875744e+02 0.000000e+00 7.215377e+02 1.728540e+02 0.000000e+00 0.000000e+00 0.000000e+00 1.000000e+00 0.000000e+00
P2: 7.070493e+02 0.000000e+00 6.040814e+02 4.575831e+01 0.000000e+00 7.070493e+02 1.805066e+02 -3.454157e-01 0.000000e+00 0.000000e+00 1.000000e+00 4.981016e-03
P3: 7.215377e+02 0.000000e+00 6.095593e+02 -3.395242e+02 0.000000e+00 7.215377e+02 1.728540e+02 2.199936e+00 0.000000e+00 0.000000e+00 1.000000e+00 2.729905e-03
R0_rect: 9.999128e-01 1.009263e-02 -8.511932e-03 -1.012729e-02 9.999406e-01 -4.037671e-03 8.470675e-03 4.123522e-03 9.999556e-01
Tr_velo_to_cam: 6.927964e-03 -9.999722e-01 -2.757829e-03 -2.457729e-02 -1.162982e-03 2.749836e-03 -9.999955e-01 -6.127237e-02 9.999753e-01 6.931141e-03 -1.143899e-03 -3.321029e-01
Tr_imu_to_velo: 9.999976e-01 7.553071e-04 -2.035826e-03 -8.086759e-01 -7.854027e-04 9.998898e-01 -1.482298e-02 3.195559e-01 2.024406e-03 1.482454e-02 9.998881e-01 -7.997231e-01
这些矩阵各自承担着不同的作用：
| 矩阵名称 | 维度 | 作用描述 | 关键用途 |
|---|---|---|---|
| P0-P3 | 3×4 | 相机投影矩阵 | 将3D相机坐标投影到2D图像平面 |
| R0_rect | 3×3 | 相机0的矫正旋转矩阵 | 校正相机镜头的畸变 |
| Tr_velo_to_cam | 3×4 | 激光雷达到相机的变换矩阵 | 将点云从激光雷达坐标系转换到相机坐标系 |
| Tr_imu_to_velo | 3×4 | IMU到激光雷达的变换矩阵 | 用于定位和里程计任务 |
对于3D点云处理，最重要的是Tr_velo_to_cam矩阵，它包含了旋转和平移两部分，可以将激光雷达坐标系下的点转换到相机坐标系。
3.2 实现点云到图像的投影
让我们通过代码实现完整的坐标转换流程：
import cv2
import matplotlib.pyplot as plt
def load_calibration(calib_file):
"""
加载并解析标定文件
返回:
calib_dict: 包含所有标定矩阵的字典
"""
calib_dict = {}
with open(calib_file, 'r') as f:
for line in f:
if ':' in line:
key, values = line.split(':', 1)
key = key.strip()
# 将字符串转换为浮点数列表
values = [float(v) for v in values.strip().split()]
# 根据键名确定矩阵形状
if key.startswith('P'):
# 投影矩阵是3x4
calib_dict[key] = np.array(values).reshape(3, 4)
elif key.startswith('R'):
# 旋转矩阵是3x3
calib_dict[key] = np.array(values).reshape(3, 3)
elif key.startswith('Tr'):
# 变换矩阵是3x4
calib_dict[key] = np.array(values).reshape(3, 4)
return calib_dict
def project_velo_to_image(points_velo, calib_dict, cam_id=2):
"""
将激光雷达点云投影到相机图像平面
参数:
points_velo: 激光雷达坐标系下的点云 (N, 3) 或 (N, 4)
calib_dict: 标定参数字典
cam_id: 相机ID (0:左灰度, 1:右灰度, 2:左彩色, 3:右彩色)
返回:
points_img: 图像坐标系下的2D坐标 (N, 2)
depth: 每个点的深度值 (到相机的距离)
valid_mask: 有效的投影点掩码
"""
# 确保点云是齐次坐标 (N, 4)
if points_velo.shape[1] == 3:
# 添加齐次坐标的1
points_velo_homo = np.hstack([points_velo, np.ones((points_velo.shape[0], 1))])
else:
# 只取前3个坐标，忽略反射强度
points_velo_homo = np.hstack([points_velo[:, :3], np.ones((points_velo.shape[0], 1))])
# 1. 激光雷达到相机坐标系的转换
# Tr_velo_to_cam是3x4矩阵，需要扩展为4x4齐次变换矩阵
Tr_velo_to_cam = calib_dict[f'Tr_velo_to_cam']
Tr_velo_to_cam_homo = np.vstack([Tr_velo_to_cam, [0, 0, 0, 1]])
# 2. 应用相机矫正 (R0_rect)
R0_rect = calib_dict['R0_rect']
# 将R0_rect扩展为4x4齐次矩阵
R0_rect_homo = np.eye(4)
R0_rect_homo[:3, :3] = R0_rect
# 3. 获取指定相机的投影矩阵
P = calib_dict[f'P{cam_id}']
# 将P扩展为3x4齐次矩阵（已经是正确的形状）
# 完整的变换链: 点云 -> 相机坐标系 -> 矫正 -> 投影
# 注意矩阵乘法的顺序
points_cam = Tr_velo_to_cam_homo @ points_velo_homo.T
points_cam_rect = R0_rect_homo @ points_cam
points_img_homo = P @ points_cam_rect
# 转换为非齐次坐标
points_img = points_img_homo[:2, :] / points_img_homo[2, :]
points_img = points_img.T # 转置为 (N, 2)
# 计算深度（相机坐标系下的z值）
depth = points_cam_rect[2, :]
# 创建有效点掩码（点在相机前方且在图像范围内）
valid_mask = (depth > 0) # 深度为正（点在相机前方）
return points_img, depth, valid_mask
def visualize_projection(image_path, points_velo, calib_dict, cam_id=2):
"""
可视化点云在图像上的投影
"""
# 读取图像
image = cv2.imread(image_path)
image_rgb = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)
# 投影点云到图像
points_img, depth, valid_mask = project_velo_to_image(
points_velo, calib_dict, cam_id
)
# 只保留有效点
points_img_valid = points_img[valid_mask]
depth_valid = depth[valid_mask]
# 创建图像副本用于绘制
fig, axes = plt.subplots(1, 2, figsize=(15, 6))
# 左侧：原始图像
axes[0].imshow(image_rgb)
axes[0].set_title('原始图像')
axes[0].axis('off')
# 右侧：带点云投影的图像
axes[1].imshow(image_rgb)
# 根据深度着色（近处红色，远处蓝色）
scatter = axes[1].scatter(
points_img_valid[:, 0], points_img_valid[:, 1],
c=depth_valid, cmap='jet', s=1, alpha=0.6,
vmin=0, vmax=50 # 限制深度范围以便更好显示
)
# 添加颜色条
plt.colorbar(scatter, ax=axes[1], label='深度 (米)')
axes[1].set_title('点云投影（颜色表示深度）')
axes[1].axis('off')
plt.tight_layout()
plt.show()
# 打印统计信息
print(f"总点数: {len(points_velo)}")
print(f"有效投影点数: {np.sum(valid_mask)}")
print(f"投影比例: {np.sum(valid_mask)/len(points_velo)*100:.1f}%")
print(f"深度范围: [{depth_valid.min():.1f}, {depth_valid.max():.1f}] 米")
# 使用示例
if __name__ == "__main__":
# 加载数据
calib_file = "./data/kitti/training/calib/000000.txt"
image_file = "./data/kitti/training/image_2/000000.png"
velodyne_file = "./data/kitti/training/velodyne/000000.bin"
# 读取点云（只取前10000个点以加快演示）
points = read_bin_file_fast(velodyne_file)[:10000, :3]
# 加载标定参数
calib_dict = load_calibration(calib_file)
# 可视化投影
visualize_projection(image_file, points, calib_dict)
这个投影过程是理解多传感器数据融合的关键。通过它，我们可以：
- 验证数据对齐质量：检查点云是否准确地投影到对应的物体上
- 创建深度图像：将点云的深度信息与图像像素对应
- 数据增强：在图像上应用变换（如裁剪、旋转）时，可以同步更新点云的投影
- 多模态训练：为基于图像的检测器提供3D监督信号，或为基于点云的检测器提供纹理信息
在实际项目中，我经常使用这种投影来快速验证数据预处理是否正确。如果发现点云投影明显偏离图像中的物体，那很可能是标定参数加载错误或坐标系理解有误。
4. 可视化对比：Mayavi与Open3D的深度解析
可视化是理解3D点云数据最重要的手段之一。Python中有两个主流的3D可视化库：Mayavi和Open3D。它们各有特点，适用于不同的场景。让我分享一下在实际使用中的对比体验。
4.1 Mayavi：科研场景的传统选择
Mayavi基于VTK，是科学计算领域的经典工具。它的优势在于高度可定制性和强大的渲染能力，特别适合需要精细调整可视化效果的科研场景。
import numpy as np
from mayavi import mlab
def visualize_mayavi_advanced(points,
color_by='height',
save_path=None,
show_axes=True,
bgcolor=(0, 0, 0)):
"""
使用Mayavi进行高级点云可视化
参数:
points: 点云数据，形状(N, 3)或(N, 4)
color_by: 着色方式，可选 'height', 'distance', 'intensity'
save_path: 保存图像的路径，如果为None则不保存
show_axes: 是否显示坐标轴
bgcolor: 背景颜色，RGB元组，范围0-1
"""
# 提取坐标
x = points[:, 0]
y = points[:, 1]
z = points[:, 2]
# 根据选择计算颜色值
if color_by == 'height':
# 根据高度着色（z坐标）
color_values = z
colormap = 'coolwarm' # 冷暖色系，适合高度
title = "高度着色 (蓝色低处，红色高处)"
elif color_by == 'distance':
# 根据到原点的距离着色
distances = np.sqrt(x**2 + y**2)
color_values = distances
colormap = 'viridis' # 渐变色系
title = "距离着色 (近处紫色，远处黄色)"
elif color_by == 'intensity' and points.shape[1] >= 4:
# 根据反射强度着色
color_values = points[:, 3]
colormap = 'hot' # 热力图色系
title = "反射强度着色 (暗色弱反射，亮色强反射)"
else:
# 默认使用高度
color_values = z
colormap = 'coolwarm'
title = "高度着色"
# 创建图形
fig = mlab.figure(
size=(1200, 800),
bgcolor=bgcolor,
fgcolor=(1, 1, 1) # 前景色（文字颜色）
)
# 创建点云可视化
pts = mlab.points3d(
x, y, z, color_values,
mode='point',
colormap=colormap,
scale_factor=0.05, # 点的大小
opacity=0.8,
figure=fig
)
# 添加颜色条
mlab.colorbar(
pts,
title=color_by.capitalize(),
orientation='vertical',
label_fmt='%.1f'
)
# 设置标题
mlab.title(title, height=0.95, size=0.3)
# 设置视角
mlab.view(azimuth=45, elevation=60, distance=80, focalpoint=(0, 0, 0))
# 是否显示坐标轴
if show_axes:
mlab.axes(
xlabel='X (前进方向)',
ylabel='Y (左侧方向)',
zlabel='Z (高度方向)',
ranges=[x.min(), x.max(), y.min(), y.max(), z.min(), z.max()],
x_axis_visibility=True,
y_axis_visibility=True,
z_axis_visibility=True
)
# 添加网格平面（地面）
# 创建网格
xx, yy = np.mgrid[
x.min():x.max():50j,
y.min():y.max():50j
]
zz = np.ones_like(xx) * z.min()
# 绘制地面网格
mlab.mesh(xx, yy, zz, color=(0.3, 0.3, 0.3), opacity=0.2)
# 保存图像
if save_path:
mlab.savefig(save_path, magnification=2)
print(f"图像已保存到: {save_path}")
# 显示
mlab.show()
# Mayavi的交互功能示例
def interactive_mayavi_demo(points):
"""
展示Mayavi的交互功能
"""
# 创建基础可视化
x, y, z = points[:, 0], points[:, 1], points[:, 2]
fig = mlab.figure(size=(1000, 700))
# 按高度着色的点云
pts = mlab.points3d(x, y, z, z,
mode='point',
colormap='spectral',
scale_factor=0.1)
# 添加一些交互控件（通过注释说明，实际需要GUI）
print("Mayavi交互功能说明:")
print("1. 鼠标左键拖动: 旋转视角")
print("2. 鼠标右键拖动: 缩放")
print("3. 鼠标中键拖动: 平移")
print("4. 'r'键: 重置视角")
print("5. 's'键: 保存当前视图")
mlab.show()
# 使用示例
if __name__ == "__main__":
# 读取并预处理点云
points = read_bin_file_fast("./data/kitti/training/velodyne/000000.bin")
processed_points = preprocess_point_cloud(points)
print("Mayavi可视化演示")
print("=" * 50)
# 演示不同着色方式
print("\n1. 按高度着色（默认）")
visualize_mayavi_advanced(
processed_points[:, :3],
color_by='height',
save_path='mayavi_height.png'
)
print("\n2. 按距离着色")
visualize_mayavi_advanced(
processed_points[:, :3],
color_by='distance',
save_path='mayavi_distance.png'
)
print("\n3. 按反射强度着色")
visualize_mayavi_advanced(
processed_points,
color_by='intensity',
save_path='mayavi_intensity.png'
)
Mayavi的强大之处在于其灵活性。你可以通过调整各种参数实现高度定制化的可视化效果。但它的缺点也很明显：安装相对复杂（特别是Windows系统），启动较慢，且交互界面相对陈旧。
4.2 Open3D：现代3D处理的瑞士军刀
Open3D是一个专门为3D数据处理设计的现代库，它提供了更简洁的API和更好的性能。对于大多数应用场景，特别是需要快速原型开发时，我倾向于使用Open3D。
import open3d as o3d
import numpy as np
from matplotlib import cm
def visualize_open3d_advanced(points,
color_by='height',
point_size=2.0,
background_color=(0, 0, 0),
show_coordinate_frame=True):
"""
使用Open3D进行高级点云可视化
参数:
points: 点云数据，形状(N, 3)或(N, 4)
color_by: 着色方式，可选 'height', 'distance', 'intensity', 'rgb'
point_size: 点的大小
background_color: 背景颜色，RGB元组，范围0-1
show_coordinate_frame: 是否显示坐标系
"""
# 创建Open3D点云对象
pcd = o3d.geometry.PointCloud()
# 设置点坐标
pcd.points = o3d.utility.Vector3dVector(points[:, :3])
# 根据着色方式设置颜色
if color_by == 'height':
# 高度着色
z_values = points[:, 2]
z_min, z_max = z_values.min(), z_values.max()
# 归一化到0-1
z_normalized = (z_values - z_min) / (z_max - z_min + 1e-8)
# 使用viridis色图
colors = cm.viridis(z_normalized)[:, :3]
elif color_by == 'distance':
# 距离着色
distances = np.sqrt(points[:, 0]**2 + points[:, 1]**2)
dist_min, dist_max = distances.min(), distances.max()
dist_normalized = (distances - dist_min) / (dist_max - dist_min + 1e-8)
colors = cm.plasma(dist_normalized)[:, :3]
elif color_by == 'intensity' and points.shape[1] >= 4:
# 反射强度着色
intensity = points[:, 3]
int_min, int_max = intensity.min(), intensity.max()
int_normalized = (intensity - int_min) / (int_max - int_min + 1e-8)
colors = cm.hot(int_normalized)[:, :3]
elif color_by == 'rgb' and points.shape[1] >= 6:
# RGB颜色（如果有的话）
colors = points[:, 3:6] / 255.0
else:
# 默认使用高度着色
z_values = points[:, 2]
z_min, z_max = z_values.min(), z_values.max()
z_normalized = (z_values - z_min) / (z_max - z_min + 1e-8)
colors = cm.viridis(z_normalized)[:, :3]
pcd.colors = o3d.utility.Vector3dVector(colors)
# 创建可视化窗口
vis = o3d.visualization.Visualizer()
vis.create_window(
window_name=f"Open3D点云可视化 - {color_by}着色",
width=1200,
height=800,
left=50,
top=50
)
# 设置渲染选项
render_option = vis.get_render_option()
render_option.background_color = np.array(background_color)
render_option.point_size = point_size
render_option.light_on = True
render_option.show_coordinate_frame = show_coordinate_frame
# 添加几何体
vis.add_geometry(pcd)
# 设置视角
view_control = vis.get_view_control()
# 根据点云范围自动调整视角
bbox = pcd.get_axis_aligned_bounding_box()
center = bbox.get_center()
max_extent = bbox.get_max_extent()
# 设置合适的视角
view_control.set_front([0, -1, -0.5]) # 从斜上方看
view_control.set_lookat(center) # 看向点云中心
view_control.set_up([0, 0, 1]) # Z轴向上
view_control.set_zoom(0.8 * max_extent) # 根据范围自动缩放
# 添加坐标系（如果需要）
if show_coordinate_frame:
coordinate_frame = o3d.geometry.TriangleMesh.create_coordinate_frame(
size=2.0, # 坐标系大小
origin=[0, 0, 0] # 坐标系原点
)
vis.add_geometry(coordinate_frame)
# 运行可视化
vis.run()
vis.destroy_window()
def open3d_interactive_features():
"""
演示Open3D的交互功能
"""
print("Open3D交互功能说明:")
print("=" * 50)
print("\n鼠标控制:")
print("• 左键拖动: 旋转视角")
print("• 右键拖动: 平移场景")
print("• 滚轮滚动: 缩放")
print("• Ctrl+左键拖动: 模型旋转")
print("• Shift+左键拖动: 滚动")
print("\n键盘快捷键:")
print("• R: 重置视角")
print("• C: 切换点/线框/面片渲染")
print("• S: 截图保存")
print("• H: 显示帮助")
print("• +/-: 调整点大小")
print("\n高级功能:")
print("• 支持选择点（Shift+左键框选）")
print("• 支持测量距离（工具菜单）")
print("• 支持点云配准（编程接口）")
print("• 实时点云更新")
def compare_visualization_libraries():
"""
对比Mayavi和Open3D的特性
"""
comparison_data = [
["特性", "Mayavi", "Open3D", "推荐场景"],
["安装难度", "较高（依赖VTK）", "中等（pip直接安装）", "Open3D更友好"],
["启动速度", "较慢", "较快", "Open3D胜出"],
["API简洁性", "较复杂", "非常简洁", "Open3D更适合快速开发"],
["渲染质量", "优秀，高度可定制", "良好，现代风格", "Mayavi适合出版物"],
["交互体验", "传统，功能丰富", "现代，直观", "Open3D更易上手"],
["3D处理功能", "基础", "丰富（配准、分割等）", "Open3D是完整工具箱"],
["社区活跃度", "稳定但增长慢", "非常活跃", "Open3D有更好的支持"],
["多平台支持", "良好", "优秀", "两者都不错"],
["实时渲染", "一般", "优秀", "Open3D适合交互应用"]
]
print("Mayavi vs Open3D 详细对比")
print("=" * 60)
# 打印对比表格
for row in comparison_data:
print(f"{row[0]:<15} {row[1]:<20} {row[2]:<20} {row[3]:<20}")
print("\n" + "=" * 60)
print("个人使用建议:")
print("1. 如果是科研论文，需要高质量的可视化图，选择Mayavi")
print("2. 如果是工程应用或快速原型开发，选择Open3D")
print("3. 如果需要复杂的3D数据处理（如配准、重建），选择Open3D")
print("4. 如果项目已经大量使用VTK生态，选择Mayavi")
print("5. 对于KITTI数据可视化入门，建议从Open3D开始")
# 使用示例
if __name__ == "__main__":
# 读取数据
points = read_bin_file_fast("./data/kitti/training/velodyne/000000.bin")
processed_points = preprocess_point_cloud(points)
print("Open3D高级可视化演示")
print("=" * 50)
# 显示交互功能说明
open3d_interactive_features()
# 演示不同着色方式
print("\n\n1. 高度着色演示")
visualize_open3d_advanced(
processed_points,
color_by='height',
point_size=3.0,
background_color=(0.1, 0.1, 0.1)
)
print("\n2. 距离着色演示")
visualize_open3d_advanced(
processed_points,
color_by='distance',
point_size=2.5,
background_color=(0, 0, 0.1) # 深蓝色背景
)
print("\n3. 反射强度着色演示")
visualize_open3d_advanced(
processed_points,
color_by='intensity',
point_size=2.0,
background_color=(0.1, 0, 0) # 深红色背景
)
# 显示对比分析
compare_visualization_libraries()
Open3D的API设计更加现代化和直观。它的PointCloud类封装了丰富的功能，包括点云滤波、配准、分割等。在实际项目中，我特别喜欢Open3D的以下几点：
- 简洁的API：几行代码就能实现复杂可视化
- 丰富的3D算法：内置了常用的点云处理算法
- 良好的性能：支持大规模点云的实时渲染
- 活跃的社区：更新频繁，bug修复及时
4.3 实战：创建交互式点云查看器
结合两者的优点，我们可以创建一个更实用的点云查看器。以下是一个结合了Open3D交互性和Mayavi渲染质量的示例：
import numpy as np
import open3d as o3d
from matplotlib import cm
import threading
import time
class InteractivePointCloudViewer:
"""
交互式点云查看器类
"""
def __init__(self, points, window_name="点云查看器"):
"""
初始化查看器
参数:
points: 点云数据，形状(N, 3)或(N, 4)
window_name: 窗口名称
"""
self.points = points
self.window_name = window_name
self.current_colormap = 'viridis'
self.current_point_size = 2.0
self.is_running = False
# 创建Open3D点云对象
self.pcd = o3d.geometry.PointCloud()
self.pcd.points = o3d.utility.Vector3dVector(points[:, :3])
# 初始颜色（按高度）
self.update_colors('height')
def update_colors(self, color_by='height'):
"""
更新点云颜色
参数:
color_by: 着色方式
"""
if color_by == 'height':
z_values = self.points[:, 2]
z_min, z_max = z_values.min(), z_values.max()
z_normalized = (z_values - z_min) / (z_max - z_min + 1e-8)
colors = cm.viridis(z_normalized)[:, :3]
elif color_by == 'distance':
distances = np.sqrt(self.points[:, 0]**2 + self.points[:, 1]**2)
dist_min, dist_max = distances.min(), distances.max()
dist_normalized = (distances - dist_min) / (dist_max - dist_min + 1e-8)
colors = cm.plasma(dist_normalized)[:, :3]
elif color_by == 'intensity' and self.points.shape[1] >= 4:
intensity = self.points[:, 3]
int_min, int_max = intensity.min(), intensity.max()
int_normalized = (intensity - int_min) / (int_max - int_min + 1e-8)
colors = cm.hot(int_normalized)[:, :3]
else:
# 默认颜色
colors = np.ones((len(self.points), 3)) * 0.7
self.pcd.colors = o3d.utility.Vector3dVector(colors)
self.current_colormap = color_by
def update_point_size(self, size):
"""
更新点的大小
"""
self.current_point_size = max(0.1, min(size, 10.0))
def create_custom_ui(self):
"""
创建自定义UI控件
"""
# 创建可视化器
self.vis = o3d.visualization.VisualizerWithKeyCallback()
self.vis.create_window(
window_name=self.window_name,
width=1400,
height=900
)
# 添加点云
self.vis.add_geometry(self.pcd)
# 设置渲染选项
render_option = self.vis.get_render_option()
render_option.background_color = np.array([0.1, 0.1, 0.1])
render_option.point_size = self.current_point_size
render_option.light_on = True
# 添加坐标系
coordinate_frame = o3d.geometry.TriangleMesh.create_coordinate_frame(
size=3.0,
origin=[0, 0, 0]
)
self.vis.add_geometry(coordinate_frame)
# 注册键盘回调
self.register_key_callbacks()
# 设置初始视角
self.reset_view()
def register_key_callbacks(self):
"""
注册键盘回调函数
"""
# 'H'键：高度着色
self.vis.register_key_callback(ord("H"), self.key_callback_height)
# 'D'键：距离着色
self.vis.register_key_callback(ord("D"), self.key_callback_distance)
# 'I'键：强度着色
self.vis.register_key_callback(ord("I"), self.key_callback_intensity)
# '+'键：增大点大小
self.vis.register_key_callback(ord("+"), self.key_callback_increase_size)
# '-'键：减小点大小
self.vis.register_key_callback(ord("-"), self.key_callback_decrease_size)
# 'R'键：重置视角
self.vis.register_key_callback(ord("R"), self.key_callback_reset)
# 'S'键：保存截图
self.vis.register_key_callback(ord("S"), self.key_callback_screenshot)
# 'Q'键：退出
self.vis.register_key_callback(ord("Q"), self.key_callback_quit)
def key_callback_height(self, vis):
"""高度着色回调"""
self.update_colors('height')
vis.update_geometry(self.pcd)
print("切换到高度着色模式")
return False
def key_callback_distance(self, vis):
"""距离着色回调"""
self.update_colors('distance')
vis.update_geometry(self.pcd)
print("切换到距离着色模式")
return False
def key_callback_intensity(self, vis):
"""强度着色回调"""
if self.points.shape[1] >= 4:
self.update_colors('intensity')
vis.update_geometry(self.pcd)
print("切换到反射强度着色模式")
else:
print("当前点云没有反射强度信息")
return False
def key_callback_increase_size(self, vis):
"""增大点大小"""
self.current_point_size = min(10.0, self.current_point_size + 0.5)
vis.get_render_option().point_size = self.current_point_size
print(f"点大小增加到: {self.current_point_size}")
return False
def key_callback_decrease_size(self, vis):
"""减小点大小"""
self.current_point_size = max(0.5, self.current_point_size - 0.5)
vis.get_render_option().point_size = self.current_point_size
print(f"点大小减小到: {self.current_point_size}")
return False
def key_callback_reset(self, vis):
"""重置视角"""
self.reset_view()
print("视角已重置")
return False
def key_callback_screenshot(self, vis):
"""保存截图"""
timestamp = time.strftime("%Y%m%d_%H%M%S")
filename = f"screenshot_{timestamp}.png"
vis.capture_screen_image(filename)
print(f"截图已保存: {filename}")
return False
def key_callback_quit(self, vis):
"""退出程序"""
print("退出查看器")
self.is_running = False
vis.destroy_window()
return False
def reset_view(self):
"""重置视角到默认位置"""
ctr = self.vis.get_view_control()
ctr.set_front([0, -1, -0.5])
ctr.set_lookat([0, 0, 0])
ctr.set_up([0, 0, 1])
ctr.set_zoom(0.8)
def print_help(self):
"""打印帮助信息"""
help_text = """
=== 点云查看器帮助 ===
鼠标控制:
左键拖动: 旋转视角
右键拖动: 平移场景
滚轮: 缩放
键盘快捷键:
H: 高度着色
D: 距离着色
I: 反射强度着色（如果可用）
+: 增大点大小
-: 减小点大小
R: 重置视角
S: 保存截图
Q: 退出查看器
当前状态:
着色模式: {}
点大小: {:.1f}
点数: {}
""".format(
self.current_colormap,
self.current_point_size,
len(self.points)
)
print(help_text)
def run(self):
"""运行查看器"""
self.is_running = True
self.create_custom_ui()
self.print_help()
# 在单独线程中运行，以便可以处理键盘输入
def run_visualizer():
self.vis.run()
self.is_running = False
# 启动可视化线程
vis_thread = threading.Thread(target=run_visualizer)
vis_thread.daemon = True
vis_thread.start()
# 主线程等待
try:
while self.is_running:
time.sleep(0.1)
except KeyboardInterrupt:
print("\n用户中断")
finally:
if hasattr(self, 'vis'):
self.vis.destroy_window()
# 使用示例
if __name__ == "__main__":
print("启动交互式点云查看器")
print("=" * 50)
# 加载并预处理点云
points = read_bin_file_fast("./data/kitti/training/velodyne/000000.bin")
processed_points = preprocess_point_cloud(points, max_range=50.0)
print(f"加载点云: {len(processed_points)} 个点")
print("按任意键开始可视化，按Q键退出")
# 创建并运行查看器
viewer = InteractivePointCloudViewer(
processed_points,
window_name="KITTI点云交互查看器"
)
viewer.run()
print("查看器已关闭")
这个交互式查看器提供了实时的点云探索功能，你可以：
- 用鼠标旋转、平移、缩放视角
- 用快捷键切换不同的着色模式
- 调整点的大小以获得最佳视觉效果
- 保存当前视图的截图
- 实时查看点云的统计信息
在实际工作中，我经常使用这样的工具来快速检查数据质量、验证预处理效果，或者向非技术同事展示点云数据的特性。这种交互性对于理解3D数据的空间关系非常有帮助。
5. 从可视化到分析：点云数据处理实战技巧
掌握了基础的可视化之后，我们可以进一步探索点云数据的分析技巧。这些技巧在实际的3D感知项目中非常有用。
5.1 点云统计分析
理解点云的统计特性对于设计算法和调试模型至关重要。以下是一些实用的统计分析函数：
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
def analyze_point_cloud_statistics(points, sample_name="点云样本"):
"""
对点云进行全面的统计分析
返回:
stats_dict: 包含各种统计指标的字典
"""
# 基本统计
x, y, z = points[:, 0], points[:, 1], points[:, 2]
stats_dict = {
'sample_name': sample_name,
'total_points': len(points),
'density': len(points) / (x.ptp() * y.ptp()), # 点密度（点/平方米）
# 坐标范围
'x_range': [x.min(), x.max(), x.ptp()],
'y_range': [y.min(), y.max(), y.ptp()],
'z_range': [z.min(), z.max(), z.ptp()],
# 中心统计
'center': [x.mean(), y.mean(), z.mean()],
'median': [np.median(x), np.median(y), np.median(z)],
# 分布统计
'x_std': x.std(),
'y_std': y.std(),
'z_std': z.std(),
# 偏度和峰度
'x_skewness': stats.skew(x),
'y_skewness': stats.skew(y),
'z_skewness': stats.skew(z),
'x_kurtosis': stats.kurtosis(x),
'y_kurtosis': stats.kurtosis(y),
'z_kurtosis': stats.kurtosis(z),
}
# 如果有反射强度信息
if points.shape[1] >= 4:
intensity = points[:, 3]
stats_dict.update({
'intensity_mean': intensity.mean(),
'intensity_std': intensity.std(),
'intensity_range': [intensity.min(), intensity.max()],
'intensity_skewness': stats.skew(intensity),
'intensity_kurtosis': stats.kurtosis(intensity),
})
# 距离分布
distances = np.sqrt(x**2 + y**2)
stats_dict.update({
'mean_distance': distances.mean(),
'max_distance': distances.max(),
'distance_std': distances.std(),
# 距离分位数
'distance_percentiles': {
'25%': np.percentile(distances, 25),
'50%': np.percentile(distances, 50),
'75%': np.percentile(distances, 75),
'90%': np.percentile(distances, 90),
'95%': np.percentile(distances, 95),
}
})
return stats_dict
def visualize_statistics(points, stats_dict):
"""
可视化点云统计信息
"""
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
fig.suptitle(f"点云统计分析 - {stats_dict['sample_name']}", fontsize=16)
x, y, z = points[:, 0], points[:, 1], points[:, 2]
# 1. 3D散点图（颜色表示高度）
ax = axes[0, 0]
scatter = ax.scatter(x, y, c=z, cmap='viridis', s=1, alpha=0.6)
ax.set_xlabel('X (前进方向)')
ax.set_ylabel('Y (左侧方向)')
ax.set_title('3D点云投影（XY平面）')
plt.colorbar(scatter, ax=ax, label='高度 Z')
ax.grid(True, alpha=0.3)
# 2. 高度分布直方图
ax = axes[0, 1]
ax.hist(z, bins=50, alpha=0.7, color='skyblue', edgecolor='black')
ax.axvline(z.mean(), color='red', linestyle='--', label=f'均值: {z.mean():.2f}')
ax.axvline(np.median(z), color='green', linestyle=':', label=f'中位数: {np.median(z):.2f}')
ax.set_xlabel('高度 Z')
ax.set_ylabel('频数')
ax.set_title('高度分布')
ax.legend()
ax.grid(True, alpha=0.3)
# 3. 距离分布
ax = axes[0, 2]
distances = np.sqrt(x**2 + y**2)
ax.hist(distances, bins=50, alpha=0.7, color='lightcoral', edgecolor='black')
ax.axvline(distances.mean(), color='red', linestyle='--',
label=f'平均距离: {distances.mean():.2f}m')
ax.set_xlabel('到原点的距离 (米)')
ax.set_ylabel('频数')
ax.set_title('距离分布')
ax.legend()
ax.grid(True, alpha=0.3)
# 4. 2D密度图
ax = axes[1, 0]
hb = ax.hexbin(x, y, gridsize=50, cmap='inferno', bins='log')
ax.set_xlabel('X (前进方向)')
ax.set_ylabel('Y (左侧方向)')
ax.set_title('点云密度分布')
plt.colorbar(hb, ax=ax, label='对数频数')
ax.grid(True, alpha=0.3)
# 5. 反射强度分析（如果有）
if points.shape[1] >= 4:
ax = axes[1, 1]
intensity = points[:, 3]
ax.hist(intensity, bins=50, alpha=0.7, color='gold', edgecolor='black')
ax.set_xlabel('反射强度')
ax.set_ylabel('频数')
ax.set_title('反射强度分布')
ax.grid(True, alpha=0.3)
else:
axes[1, 1].axis('off')
axes[1, 1].text(0.5, 0.5, '无反射强度数据',
ha='center', va='center', transform=axes[1, 1].transAxes)
# 6. 统计摘要
ax = axes[1, 2]
ax.axis('off')
summary_text = f"""
统计摘要:
================
总点数: {stats_dict['total_points']:,}
点密度: {stats_dict['density']:.2f} 点/平方米
坐标范围:
X: [{stats_dict['x_range'][0]:.1f}, {stats_dict['x_range'][1]:.1f}]
(跨度: {stats_dict['x_range'][2]:.1f}m)
Y: [{stats_dict['y_range'][0]:.1f}, {stats_dict['y_range'][1]:.1f}]
(跨度: {stats_dict['y_range'][2]:.1f}m)
Z: [{stats_dict['z_range'][0]:.1f}, {stats_dict['z_range'][1]:.1f}]
(跨度: {stats_dict['z_range'][2]:.1f}m)
中心位置:
均值: ({stats_dict['center'][0]:.1f},
{stats_dict['center'][1]:.1f},
{stats_dict['center'][2]:.1f})
中位数: ({stats_dict['median'][0]:.1f},
{stats_dict['median'][1]:.1f},
{stats_dict['median'][2]:.1f})
距离统计:
平均距离: {stats_dict['mean_distance']:.1f}m
最远点: {stats_dict['max_distance']:.1f}m
距离标准差: {stats_dict['distance_std']:.1f}m
"""
ax.text(0, 1, summary_text, fontfamily='monospace',
verticalalignment='top', fontsize=9)
plt.tight_layout()
plt.show()
# 打印详细统计
print(f"\n详细统计信息 - {stats_dict['sample_name']}")
print("=" * 60)
for key, value in stats_dict.items():
if key not in ['sample_name', 'distance_percentiles']:
if isinstance(value, list):
print(f"{key:20}: {value}")
else:
print(f"{key:20}: {value:.4f}")
print(f"\n距离百分位数:")
for percentile, dist in stats_dict['distance_percentiles'].items():
print(f" {percentile:5}: {dist:.2f}m")
# 使用示例
if __name__ == "__main__":
# 加载数据
points = read_bin_file_fast("./data/kitti/training/velodyne/000000.bin")
# 预处理
processed_points = preprocess_point_cloud(points, max_range=60.0)
# 统计分析
print("开始点云统计分析...")
stats = analyze_point_cloud_statistics(processed_points, "KITTI样本000000")
# 可视化
visualize_statistics(processed_points, stats)
# 比较多个样本
print("\n" + "="*60)
print("多样本统计对比")
print("="*60)
sample_files = [
"./data/kitti/training/velodyne/000000.bin",
"./data/kitti/training/velodyne/000001.bin",
"./data/kitti/training/velodyne/000002.bin"
]
comparison_results = []
for i, file_path in enumerate(sample_files[:3]): # 只分析前3个样本
try:
points_sample = read_bin_file_fast(file_path)
points_processed = preprocess_point_cloud(points_sample, max_range=60.0)
stats_sample = analyze_point_cloud_statistics(
points_processed,
f"样本{i:06d}"
)
comparison_results.append(stats_sample)
print(f"\n样本 {i}: {file_path}")
print(f" 点数: {stats_sample['total_points']:,}")
print(f" 平均距离: {stats_sample['mean_distance']:.1f}m")
print(f" 高度范围: [{stats_sample['z_range'][0]:.1f}, "
f"{stats_sample['z_range'][1]:.1f}]")
except Exception as e:
print(f"处理样本 {file_path} 时出错: {e}")
这种统计分析不仅有助于理解数据特性，还能为后续的算法设计提供重要参考。比如，通过分析点云密度，我们可以确定合适的体素化分辨率；通过分析距离分布，我们可以设计有效的距离相关特征。
5.2 点云分割与聚类
在实际的3D感知任务中，我们经常需要将点云分割成不同的物体或区域。以下是一个基于简单几何特征的点云分割示例：
from sklearn.cluster import DBSCAN
from scipy.spatial import KDTree
def segment_ground_plane(points, height_threshold=-1.2, distance_threshold=0.3):
"""
使用高度阈值法分割地面点
参数:
points: 输入点云 (N, 3)
height_threshold: 地面高度阈值
distance_threshold: 地面点最大距离阈值
返回:
ground_indices: 地面点的索引
non_ground_indices: 非地面点的索引
"""
# 简单的高度阈值分割
ground_mask = points[:, 2] < height_threshold
# 进一步过滤：地面点应该相对平坦
if np.sum(ground_mask) > 0:
ground_points = points[ground_mask]
# 使用DBSCAN聚类移除离群点
if len(ground_points) > 100:
clustering = DBSCAN(eps=distance_threshold, min_samples=10).fit(ground_points[:, :2])
largest_cluster = np.bincount(clustering.labels_[clustering.labels_ >= 0]).argmax()
refined_ground_mask = clustering.labels_ == largest_cluster
# 更新地面点索引
ground_indices = np.where(ground_mask)[0][refined_ground_mask]
else:
ground_indices = np.where(ground_mask)[0]
else:
ground_indices = np.array([], dtype=int)
# 非地面点
all_indices = np.arange(len(points))
non_ground_indices = np.setdiff1d(all_indices, ground_indices)
return ground_indices, non_ground_indices
def cluster_objects(points, eps=0.5, min_samples=10):
"""
使用DBSCAN聚类分割物体
参数:
points: 输入点云 (N, 3)，建议先移除地面点
eps: 邻域搜索半径
min_samples: 形成核心点所需的最小点数
返回:
labels: 每个点的聚类标签，-1表示噪声点
n_clusters: 发现的聚类数量
"""
# 使用DBSCAN进行聚类
clustering = DBSCAN(eps=eps, min_samples=min_samples).fit(points)
labels = clustering.labels_
# 统计聚类结果
unique_labels = np.unique(labels)
n_clusters = len(unique_labels) - (1 if -1 in unique_labels else 0)
n_noise = np.sum(labels == -1)
print(f"发现 {n_clusters} 个聚类")
print(f"噪声点数量: {n_noise}")
print(f"噪声点比例: {n_noise/len(points)*100:.1f}%")
return labels, n_clusters
def visualize_segmentation(points, ground_indices, object_labels):
"""
可视化分割结果
"""
import matplotlib.pyplot as plt
from matplotlib import cm
fig = plt.figure(figsize=(15, 5))
# 1. 原始点云
ax1 = fig.add_subplot(131, projection='3d')
ax1.scatter(points[:, 0], points[:, 1], points[:, 2],
c='blue', s=1, alpha=0.5)
ax1.set_xlabel('X')
ax1.set_ylabel('Y')
ax1.set_zlabel('Z')
ax1.set_title('原始点云')
ax1.view_init(elev=30, azim=45)
# 2. 地面分割结果
ax2 = fig.add_subplot(132, projection='3d')
# 地面点（绿色）
if len(ground_indices) > 0:
ground_points = points[ground_indices]
ax2.scatter(ground_points[:, 0], ground_points[:, 1], ground_points[:, 2],
c='green', s=1, alpha=0.7, label='地面')
# 非地面点（红色）
all_indices = np.arange(len(points))
non_ground_indices = np.setdiff1d(all_indices, ground_indices)
if len(non_ground_indices) > 0:
non_ground_points = points[non_ground_indices]
ax2.scatter(non_ground_points[:, 0], non_ground_points[:, 1], non_ground_points[:, 2],
c='red', s=1, alpha=0.5, label='非地面')
ax2.set_xlabel('X')
ax2.set_ylabel('Y')
ax2.set_zlabel('Z')
ax2.set_title('地面分割结果')
ax2.legend()
ax2.view_init(elev=30, azim=45)
# 3. 物体聚类结果
ax3 = fig.add_subplot(133, projection='3d')
# 为每个聚类分配颜色
unique_labels = np.unique(object_labels)
colors = cm.rainbow(np.linspace(0, 1, len(unique_labels)))
for label, color in zip(unique_labels, colors):
if label == -1:
# 噪声点用灰色
cluster_mask = object_labels == label
if np.any(cluster_mask):
cluster_points = points[cluster_mask]
ax3.scatter(cluster_points[:, 0], cluster_points[:, 1], cluster_points[:, 2],
c='gray', s=1, alpha=0.3, label='噪声')
else:
# 聚类点用彩色
cluster_mask = object_labels == label
cluster_points = points[cluster_mask]
ax3.scatter(cluster_points[:, 0], cluster_points[:, 1], cluster_points[:, 2],
color=color, s=2, alpha=0.7, label=f'聚类{label}')
ax3.set_xlabel('X')
ax3.set_ylabel('Y')
ax3.set_zlabel('Z')
ax3.set_title('物体聚类结果')
ax3.legend(loc='upper right', fontsize='small')
ax3.view_init(elev=30, azim=45)
plt.tight_layout()
plt.show()
def extract_cluster_features(points, labels):
"""
提取每个聚类的特征
返回:
features_dict: 每个聚类的特征字典
"""
unique_labels = np.unique(labels)
features_dict = {}
for label in unique_labels:
if label == -1:
continue # 跳过噪声点
cluster_mask = labels == label
cluster_points = points[cluster_mask]
if len(cluster_points) < 10:
continue # 跳过太小的聚类
# 计算基本特征
centroid = np.mean(cluster_points, axis=0)
bbox_min = np.min(cluster_points, axis=0)
bbox_max = np.max(cluster_points, axis=0)
bbox_size = bbox_max - bbox_min
# 计算方向（主成分分析）
centered = cluster_points - centroid
cov_matrix = np.cov(centered.T)
eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)
# 最大特征值对应的特征向量是主要方向
main_direction = eigenvectors[:, np.argmax(eigenvalues)]
features_dict[label] = {
'n_points': len(cluster_points),
'centroid': centroid,
'bbox_min': bbox_min,
'bbox_max': bbox_max,
'bbox_size': bbox_size,
'volume': np.prod(bbox_size),
'main_direction': main_direction,
'eigenvalues': eigenvalues,
'compactness': eigenvalues[0] / (np.sum(eigenvalues) + 1e-8), # 紧凑性
'height': bbox_size[2], # 高度
'width': bbox_size[1], # 宽度
'length': bbox_size[0], # 长度
}
return features_dict
# 使用示例
if __name__ == "__main__":
print("点云分割与聚类演示")
print("=" * 50)
# 加载数据
points = read_bin_file_fast("./data/kitti/training/velodyne/000000.bin")
# 预处理：移除过远的点
distances = np.sqrt(points[:, 0]**2 + points[:, 1]**2)
near_mask = distances < 50.0
这里是“一人公司”的成长家园。我们提供从产品曝光、技术变现到法律财税的全栈内容，并连接云服务、办公空间等稀缺资源，助你专注创造，无忧运营。
更多推荐
11
0- 0
已为社区贡献3条内容
所有评论(0)