欧氏距离聚类与KDTree实现LiDAR点云单木分割工程详解

简介：面向点云处理、林业遥感及三维分割方向的学习者与研究人员，该代码包基于PCL欧氏距离聚类实现单木分割，解决从原始点云中自动分离单棵树木的流程搭建与算法复现问题。资源共3个文件，包含1个cpp聚类主程序、1个xyz点云样例数据以及1个txt聚类结果文件，压缩包大小仅93KB，轻量易用，便于快速移植到自己的实验环境。算法主程序完整演示了PCL中欧氏距离聚类的调用流程，配合输入的点云样例可一键运行，生成的txt分割结果可直观对比不同距离阈值下的聚类效果，便于理解参数对单木分割精度的影响。资源包同时适合作为课程设计或科研实验的起点代码，已有813人浏览学习，适合具备一定PCL基础、希望掌握单木分割实现细节的开发者参考。

1. 欧氏距离聚类做单木分割：为什么不用人工标定？

林业调查里最头痛的不是测树高，而是面对一片LiDAR点云，得把每一棵树单独从杂乱的“点云海”里捞出来。人工在软件里一棵棵描框、切割、编号，一片样地少说花半天，而且不同人切出来的树位置还不一样，复现性差。这个项目直接用欧氏距离聚类，把空间上彼此靠近的激光点自动归到同一棵树上，输出每一棵树的点簇和几何中心，同时把整个流程封装成API——你只需把点云文件传进去，就能拿到“树号 + 点云 + 质心坐标”的结构化结果。适合刚接触点云处理的算法工程师，也适合林业信息化系统想快速接入单木分割能力的后端开发者。它不搞语义分割那套复杂模型，只靠一个几何距离阈值，就把单木分割这件事做成了开箱即用的服务。

2. 技术架构与工程结构：从点云输入到单木簇输出的完整链路

2.1 技术选型：为什么用KDTree而不是暴力计算

点云单木分割本质上是一个“空间近邻聚类”问题：如果A点与B点之间欧氏距离小于给定阈值，就把它们视为同一棵树上的点。最朴素的做法是双重循环遍历所有点，但一片10万平方米的森林点云往往有数百万个点，O(n²)的复杂度谁也扛不住。所以这里选择用KDTree做近邻搜索，它能把最近邻查询的平均复杂度降到O(n log n)，这也是点云处理里最常见的加速手段。

底层依赖只用了三个：NumPy管理数组，SciPy的spatial.KDTree做近邻索引，FastAPI提供HTTP接口。为什么不引入PCL或Open3D？因为这里只做空间聚类，不需要复杂的网格重建或可视化，开源点云库反而会让部署体积变大。用SciPy的KDTree已经足够，且它内部用的是C++实现的，查询性能很稳。

2.2 目录结构：一个最小可跑的API工程

整个工程按“核心算法 + API层 + 配置”三层拆分，便于以后把算法单独抽出来做离线批处理或直接内嵌到其他系统。目录结构如下：

pointcloud_treeseg/ ├── app.py # FastAPI 应用，定义HTTP路由 ├── config.yaml # 阈值、采样密度、输出格式配置 ├── core/ │ ├── __init__.py │ ├── preprocess.py # 点云清洗：体素降采样 + 离群点过滤 │ ├── clustering.py # 欧氏距离聚类核心算法 │ └── postprocess.py # 把聚类结果整理成单木点簇和属性表 ├── models/ │ └── schemas.py # API 请求/响应模型定义 ├── tests/ │ └── test_segmentation.py # 用小雨量合成点云做回归测试 └── requirements.txt

config.yaml是调参最常用的入口。下面是一个基本配置示例：

clustering: max_distance: 1.5 # 聚类最大欧氏距离，单位：米 min_points_per_tree: 10 # 少于10个点的簇不认为是树 leaf_size: 16 # KDTree叶子节点大小 preprocess: voxel_size: 0.2 # 体素降采样网格大小 outlier_std_ratio: 2.0 # 离群点过滤标准差倍数

min_points_per_tree很关键，LiDAR扫描时地面和低矮灌丛会产生零散点，这些点数量少但也会聚成小簇，不把它过滤掉，会把草丛当成树。

2.3 API接口设计：请求与响应格式

API只暴露两个端点：

/api/segment请求使用multipart/form-data上传点云文件，支持.xyz、.txt、.las（需自行安装laspy扩展）。响应是一个JSON数组，每个元素对应一棵树：

{ "tree_id": 1, "point_count": 356, "centroid": [31.24, 42.87, 8.95], "bbox": [30.1, 41.2, 7.5, 32.3, 44.6, 10.2], "points": "base64_string" }

这里points用Base64编码后的二进制数组传输，避免JSON里嵌套大数组导致网络开销过大。如果只是一次性离线处理，也可以直接调用core.clustering模块，不走HTTP。

3. 核心模块拆解：欧氏距离聚类的完整实现

3.1 点云数据预处理：降噪与下采样

原始点云不能直接聚类。机载LiDAR扫描时，同一位置可能被多次扫描，导致局部点密度极高；同时还有空气中灰尘、飞鸟产生的孤立点。我们的做法是先做体素降采样，让点密度均匀化，再过滤离群点。

import numpy as np from scipy.spatial import KDTree def voxel_downsample(points, voxel_size=0.2): """体素降采样：将点云网格化，每个网格只保留重心点。 Args: points: (N,3) 数组，每行一个点的 xyz 坐标 voxel_size: 体素边长（米） Returns: 降采样后的点云数组 (M,3) """ if voxel_size <= 0: return points coords = np.floor(points / voxel_size).astype(np.int32) # 用字典按键聚合，把同一个体素内的点求平均 voxel_dict = {} for p, c in zip(points, coords): key = tuple(c) if key not in voxel_dict: voxel_dict[key] = [p, 1] else: voxel_dict[key][0] += p voxel_dict[key][1] += 1 return np.array([v[0] / v[1] for v in voxel_dict.values()])

这段代码里，coords是把点坐标除以网格大小后向下取整得到的整数坐标，相同整数坐标的点落在同一个体素里。voxel_dict在统计总和的同时记录点数，最后除以点数得到重心。体素大小直接决定降采样强度：设大了树冠轮廓会变模糊，设小了去重效果不明显，我一般从0.2米试起。

离群点过滤用的是统计法：先计算每个点与它最近32个邻居的平均距离，如果平均距离超出全局均值两倍标准差，就认为它是离群的。这一步能有效剔除高空噪点，避免它们在聚类时形成微型单木簇。

3.2 近邻搜索与聚类核心

聚类算法是整条链路的灵魂。我们采用经典的“种子生长法”：随机挑一个未访问点作为种子，用KDTree找出它半径max_distance内的所有邻居，把这些邻居再作为新种子，继续向外扩展，直到没有新点加入。这本质上是广度优先搜索。

def euclidean_cluster(points, max_distance=1.5, min_points=10): """基于欧氏距离的点云聚类。 Args: points: (N,3) 点云数组 max_distance: 聚类半径（米），两点距离小于此值则属于同一簇 min_points: 簇内最少点数，低于此值视为噪点 Returns: labels: (N,) int 数组，每个点对应的簇标签，-1 表示噪点 """ tree = KDTree(points, leaf_size=16) N = len(points) labels = -np.ones(N, dtype=np.int32) visited = np.zeros(N, dtype=bool) cluster_id = 0 for i in range(N): if visited[i]: continue # BFS 队列 queue = [i] visited[i] = True cluster_points = [] while queue: current = queue.pop() cluster_points.append(current) # 找出当前点半径 max_distance 内的邻居 neighbors = tree.query_ball_point(points[current], r=max_distance) for nb in neighbors: if not visited[nb]: visited[nb] = True queue.append(nb) if len(cluster_points) >= min_points: labels[cluster_points] = cluster_id cluster_id += 1 # else: 小于 min_points 的簇保持为 -1，视为噪点 return labels

这里tree.query_ball_point是SciPy KDTree的方法，返回当前点周围r米内所有点的索引列表。它的内部实现用了递归的盒子检测，比手动遍历邻居高效得多。leaf_size=16控制KDTree的叶子节点大小，太小则树太深，查询反而慢；太大则叶子内线性扫描开销增大，16是经验默认值。

聚类半径max_distance的取值直接决定分割效果：森林中一棵树冠的相邻点间距通常在0.1~0.5米，树冠边缘与相邻树木的间隙往往大于0.5米。所以0.3~1.5米都是合理区间，具体要看点云密度。

3.3 单木点簇的输出与后处理

聚类得到的是每个点的标签数组，还需要把它们整理成树级别的结构。后处理做了三件事：按标签分组提取点云，计算每棵树的质心和包围盒，过滤不符合树形特征的簇。

def postprocess(points, labels, min_height=2.0): """把labels转为单木列表。 Args: points: 原始点云 labels: 聚类标签 min_height: 点簇的Z轴高度差，低于此值可能是灌木或地面段 Returns: trees: list of dict，每个dict包含points, centroid, bbox """ unique_labels = set(labels[labels >= 0]) trees = [] for label in unique_labels: mask = labels == label cluster_pts = points[mask] z_diff = cluster_pts[:, 2].max() - cluster_pts[:, 2].min() if z_diff < min_height: continue # 高度不够，跳过 centroid = cluster_pts.mean(axis=0) bbox = [ cluster_pts[:, 0].min(), cluster_pts[:, 1].min(), cluster_pts[:, 2].min(), cluster_pts[:, 0].max(), cluster_pts[:, 1].max(), cluster_pts[:, 2].max() ] trees.append({ "points": cluster_pts.tolist(), "centroid": centroid.tolist(), "bbox": bbox }) return trees

这里min_height=2.0是一个启发式阈值：成熟的乔木树冠垂直落差通常超过2米，而低矮灌木和草丛的Z范围很小。通过这一步能进一步去掉误聚类的非目标物。但注意，如果测区存在幼树，高度不足但确实是目标，这个阈值需要调低。

4. 常见问题与避坑：欧氏距离聚类跑不动的那些事

4.1 聚类阈值设多大？——尺度决定一切

• 现象：设max_distance=0.5时，一棵树的点云被切成了七八个碎片；设成2.0后，相邻两棵树又被并成了一簇。
• 原因：不同测区的点云密度差异很大。机载LiDAR点密度可能是每平米几到十几点，而地面扫描仪每平米能到上千点。固定阈值没有适配点密度。
• 解决：先估算点云平均点间距。常见做法是计算每个点与最近邻的欧氏距离，取平均值，然后设置阈值为点间距的2~3倍。这样密度高的区域阈值小，密度低的区域阈值大，自适应效果好。

4.2 内存爆掉怎么办？——分块与降采样

• 现象：处理一片500万点的LAS文件，程序直接“MemoryError”或者卡死。
• 原因：KDTree在构建时要维护一个索引树，占的内存约为点数的8倍。500万点就是4000万条索引记录，加上原始数组和BFS队列，很容易超出8GB内存。
• 解决：先体素降采样把点降到100万以内，再用分块策略。把点云切成256米×256米的方块，每个方块单独聚类，最后把跨块边界邻近的簇合并。合并逻辑是：如果两个簇的质心距离小于两倍max_distance，就合并。注意方块之间要保留5米的overlap，防止树正好处于切割线上。

4.3 离群点被误判成单木？——加密度过滤

• 现象：聚类结果里出现大量只有一二十个点的“树”，而且位置在空中，明显是飞鸟或噪点。
• 原因：噪声点虽然离散，但在KDTree的query_ball_point里，只要某两三个噪点彼此距离足够近，也能聚成一个簇。而min_points设得太低，比如10，就容易放过它们。
• 解决：把min_points提高到50以上，同时增加高度约束：树簇的Z轴均值不能低于地面高度。更稳妥的方案是在预处理阶段用统计离群点过滤，把每个点的近邻平均距离大于全局均值2.5倍的点直接删掉。这样噪点数量会大幅减少。

4.4 结果空转：为什么API返回空数组？

• 现象：POST点云文件后返回[]，但点云在软件里看着明显有树。
• 原因：最常见的是坐标单位不一致。有的LAS文件是度（经纬度），有的是米。如果坐标是度，max_distance=1.5相当于1.5度的距离，聚类半径大到把所有树都并成一个簇；反之如果坐标是投影米，却用了0.001的阈值，则什么都聚不起来。
• 解决：在预处理里增加一步坐标检查：计算点云坐标差分的百分位数，如果95百分位的两点距离小于0.01，说明数据可能是度单位，需要先投影转换。或者让用户显式在config.yaml里声明coordinate_unit，默认按米处理。

5. 二次开发：把分割结果接到你的业务系统里

拿到单木点簇列表只是第一步，实际业务往往需要把这些结果接进WebGIS或统计报表。这一章给几个具体的改造方向。

输出GeoJSON——林业调查通常要求可视化。把postprocess返回的点簇转换成GeoJSON多边形，可以直接在Leaflet或Cesium里显示。做法是先把点簇的X、Y坐标投影成平面坐标，用Alpha Shape算法提取树冠边界多边形。如果不想额外引库，可以简化成外包圆：以质心为圆心，以到最远点的距离为半径画圆，作为属性标签。对多数统计场景，圆形已经够用。

并发处理大场景——API服务可以用FastAPI的异步特性加线程池。因为聚类是CPU密集操作，GIL会限制多线程，所以用concurrent.futures.ProcessPoolExecutor把不同分块交给独立进程运行。下面是一个简单的改造示意：

from fastapi import UploadFile from concurrent.futures import ProcessPoolExecutor import numpy as np executor = ProcessPoolExecutor(max_workers=4) async def segment_upload(file: UploadFile): data = np.loadtxt(file.file) # 假设是xyz文本 # 按X轴分两个块，注意overlap half_x = (data[:, 0].max() + data[:, 0].min()) / 2 left = data[data[:, 0] < half_x + 1.5] right = data[data[:, 0] > half_x - 1.5] # 提交到进程池 left_future = executor.submit(run_cluster_pipeline, left) right_future = executor.submit(run_cluster_pipeline, right) left_trees, right_trees = left_future.result(), right_future.result() # 合并相邻树簇 return merge_trees(left_trees, right_trees)

这段代码里，run_cluster_pipeline是预处理到后处理的封装函数。注意分块overlap设为1.5米（即max_distance），否则边界上的树会被切成两半。合并簇时判断两个簇的质心距离是否小于max_distance，是则合并。

换用其他距离度量——欧氏距离假设点云空间各向同性，但树木在山坡上时，坡向方向点距变化更剧烈。如果你发现分割结果里山坡上的树不完整，可以考虑把Z轴距离权重放大。实现上只需要在预处理里给Z坐标乘一个权重系数，比如points[:, 2] *= 1.2，再送进KDTree和聚类。这个技巧在数据集含明显地形起伏时非常有效，代价是树高数据会变形，输出质心和包围盒前记得除回来。

验证修改是否正确，我一般用合成点云做回归测试：生成几棵椭球形状的点簇，设置不同间距，跑一遍分割，断言每棵树被完整分开、位置误差小于0.1米。测试脚本放在tests/test_segmentation.py里，每次改代码后先跑一遍，避免越改越偏。

从那以后，我每次调整阈值或数据结构，都会强制走一遍合成数据回归测试，确认分割结果没有退化，再投到真实数据集上跑。这个习惯帮我省掉了不少在客户现场才发现问题的尴尬。希望这个项目能帮你把单木分割的坑少踩几个，让你的点云处理链路早日稳定跑起来。

本文还有配套的精品资源，点击获取