点云单木分割实战：欧氏距离聚类与FastAPI封装指南

简介面向激光雷达点云处理与林业遥感初学者的一份可直接运行的源码包借助PCL库中的欧氏距离聚类算法实现单木分割。压缩包内共3个文件一个C实现的主程序cpp一份用于测试的xyz格式点云数据以及聚类后生成的txt结果文件整个资源包仅93KB轻量且聚焦便于快速复盘整个算法流程。目前已有811人学习下载适合正在入门点云分割、希望深入理解欧氏聚类参数如聚类容差、最小/最大点数如何影响分割效果的研究者与开发者。通过对照配套博文中的原理讲解与示例输出读者能掌握从点云读取、欧氏聚类调用到单木结果导出的完整链路并可根据自身数据微调阈值参数为后续开展株数统计、单木结构参数提取等应用打下基础。资源整体结构紧凑无冗余文件适合作为算法入门或课程实验的参考实现。1. 从点云到单木分割为什么欧氏距离聚类是首选的 API 内核做林业激光扫描数据处理的人几乎都会卡在同一道坎上一片样地的点云扫出来有几十万到几千万个点每一棵树都一堆点但软件里呈现的是一整片没有 ID 的点海。要算平均树高、冠幅、株数密度第一步必须是“单木分割”——把每棵树的点从整片点云里独立拆出来。欧氏距离聚类是这个环节最常用的无监督方案逻辑直白空间距离近的点归成一棵树距离远的点属于另一棵树。它的优势在于不依赖训练数据、不生成栅格、参数少且物理含义明确因此特别适合做成一个可复用的 API 服务。这篇文章就按我从数据读取到接口封装、再到参数调优的实际操作顺序把欧氏距离实现点云单木分割的完整链路拆开讲清楚哪里容易翻车也会一并说明。2. 欧氏距离聚类分割单木的原理核心参数与实现代码2.1 单木分割的三条技术路线欧氏距离为什么值得先用单木分割的主流方案大概分三类各自适用场景差别很大。第一类是 CHM 分水岭先把点云投影成冠层高度模型栅格再对倒置的 CHM 做分水岭分割。它的优点是成熟、计算快缺点是栅格化会丢失冠层内部的垂直结构信息低矮树木容易被高大邻居遮挡树冠粘连时分割边界容易溢出。第二类是区域生长从树顶种子点出发按邻域规则逐步合并点。这个方法需要先可靠地检测种子点种子位置稍有偏差树冠就会长歪。第三类就是欧氏距离聚类直接在三维空间里做连通域分析完全绕开栅格化和种子点问题。从 API 服务的角度看默认算法应该满足一个条件输入数据不确定、场景变化大时仍然稳定。欧氏距离聚类恰好符合这一点它不针对特定树种建模不要求预先训练只依赖“树冠内部点间距显著小于树间间距”这个普遍成立的物理事实。遇到树冠紧密相接的场景它确实容易把两棵树合并到一起但这个短板可以通过后处理弥补后面避坑章节会专门讲。相比之下CHM 分水岭和区域生长在树种变化、点密度变化时调整成本更高。我在实际项目里选择欧氏距离聚类还有一个现实理由它便于逐步调试。整个流程可以拆成“读取点云→体素下采样→聚类→后处理输出”四段每一段都能单独检查中间结果这对接 API 接口排障特别重要。即使后续要换更复杂的语义分割模型欧氏距离聚类的结果也可以作为训练样本的初筛基线并不浪费。2.2 三个核心参数eps、min_points 与体素尺寸的联动关系欧氏距离聚类直接相关的参数只有两个聚类半径阈值 eps 和最小点数 min_points。但点云数据往往有几十万甚至上千万个点直接聚类速度慢且内存压力大所以通常要先做体素下采样这又引入了第三个参数 voxel_size。这三个参数不是独立的它们的联动关系决定了分割成败。先看 eps。它代表“距离多近的点算同一个簇”单位与点云坐标一致通常为米。eps 太小一棵树内部的点因为枝条间有空隙会被拆成几块eps 太大相邻两棵树的树冠会被合并。树的冠层内部点平均间距一般在 0.5 到 1.5 米之间而相邻树冠之间的最小间距经常只有 1 到 3 米这意味着 eps 的理想取值窗口并不宽。我一般会先估算树冠内部点的平均间距再把 eps 取为该间距的 2 到 3 倍这样既保证树冠内部连通又不容易跨过树间间隙。min_points 的作用是过滤小噪点簇。它的合理值取决于体素下采样后的点密度。假设体素尺寸 0.3 米一棵 10 米高的树冠大约能保留几百到上千个点那么 min_points 取 30 到 50 比较合适如果体素放到 1 米一棵树的点可能只剩几十个min_points 就得降到 10 到 20。voxel_size 决定下采样后点云的均匀程度一般取 0.2 到 0.5 米。体素越大数据量越小但树冠形状细节丢失越严重。参数之间的经验关系可以总结成一张表参数含义经验取值调整方向voxel_size体素边长决定下采样后点间距0.2–0.5 m点云越密取越小eps聚类半径阈值体素尺寸的 2–3 倍树冠越密取越小min_points簇最少点数30–500.3 m 体素噪点多时调大调参顺序也有讲究。我习惯先定 voxel_size再根据下采样后的点间距定 eps最后根据单棵树保留的点数定 min_points。反过来调会经常出现“eps 看着合理但分割结果一团糟”的玄学问题。2.3 用 cKDTree 实现欧氏距离聚类代码与复杂度分析欧氏距离聚类在实现层面等价于“半径邻域连通域分析”。对点云中的任意一点半径为 eps 范围内所有点属于同一个簇簇继续向外扩展直到边界处没有新的点加入。这个逻辑直接用双重循环遍历所有点对是行不通的N 个点的距离矩阵是 N 的平方量级百万点会直接卡死。正确做法是构建 KD 树加速最近邻查询代码里使用 scipy.spatial.cKDTree。import numpy as np from scipy.spatial import cKDTree def euclidean_cluster(points, eps, min_points): 基于欧氏距离的连通域聚类。 points: (N, 3) 数组通常为 XY Z 坐标单位为米 eps: 聚类半径阈值 min_points: 聚类最少点数小于该值视为噪声 返回: labels 数组每个点对应一个簇编号-1 表示噪声 n len(points) if n 0: return np.array([], dtypeint) tree cKDTree(points) labels np.full(n, -1, dtypeint) visited np.zeros(n, dtypebool) cluster_id 0 for i in range(n): if visited[i]: continue stack [i] visited[i] True members [] while stack: idx stack.pop() members.append(idx) # 查询半径 eps 内所有邻居索引包含自身 neighbors tree.query_ball_point(points[idx], eps) for j in neighbors: if not visited[j]: visited[j] True stack.append(j) if len(members) min_points: labels[members] cluster_id cluster_id 1 return labels这段代码的核心是tree.query_ball_point它一次返回某个点半径 eps 内的所有邻居索引比逐点计算距离快几个数量级。cKDTree 建树复杂度大约是 O(N log N)每个点的邻域查询平均 O(log N)整体聚类过程接近线性百万点规模在普通笔记本上几秒内能跑完。labels 数组中值为 -1 的是噪点其余簇编号从 0 开始递增。这个实现其实等价于一个 mini 版 DBSCAN只不过固定了欧氏距离度量不需要额外依赖 scikit-learn封装成 API 时依赖更少。要注意一个细节query_ball_point返回的邻居包含点自身但因为当前点已经标记为 visited不会重复入栈。这个写法在面对点云密度不均时表现尚可但遇到两棵树的树冠在空间上交错、点间距小于 eps 的场景仍然会合并这是算法本身的局限后处理环节再解决。3. 从 LAS 到单木点云一条完整的最小实现流程3.1 LAS 读取、体素下采样与地面点过滤林业激光扫描点云最常见的格式是 LAS/LAZ读取用 laspy 库。读完之后拿到的是三维坐标数组后续所有计算都基于 numpy 数组完成。这里有个容易忽略的问题laspy 1.x 和 2.x 的属性接口不一样2.x 使用las.x、las.y、las.z老代码里常见的las.X在新版本已经废弃后面避坑章节会细说。体素下采样这一步我用质心代替体素内的所有点而不是随机取一个点。原因是随机取点会引入抖动导致树冠轮廓不均匀质心则能保留体素内的平均位置。np.add.at 是处理重复索引累加的关键函数比循环快很多。import laspy import numpy as np def read_las_points(file_path): 读取 LAS/LAZ 点云返回 (N, 3) 的 XYZ 数组单位通常为米。 las laspy.read(file_path) points np.vstack((las.x, las.y, las.z)).T return las, points def voxel_downsample(points, voxel_size): 体素下采样每个体素内用质心点代替该体素内所有点。 if voxel_size 0: return points voxel_coords np.floor(points / voxel_size).astype(np.int64) _, inverse_indices np.unique(voxel_coords, axis0, return_inverseTrue) summed np.zeros((inverse_indices.max() 1, 3), dtypenp.float64) np.add.at(summed, inverse_indices, points) counts np.bincount(inverse_indices).reshape(-1, 1) return summed / countsvoxel_size的单位和点云坐标一致一般给 0.3 米。为什么用np.floor而不是np.round因为 floor 能保证同一个体素内的点坐标范围一致左闭右开round 在边界处可能把相邻两个体素的点混到一起。质心下采样后点数大幅减少聚类速度明显提升同时树冠形状的关键特征仍然保留。地面点和低矮植被是单木分割的干扰项。LAS 标准里类别 2 是地面点如果点云已经做过分类直接按类别过滤。很多数据没有分类信息就用高度阈值过滤。def remove_ground_and_low(las, points, min_height2.0): 剔除地面点与低矮点优先用 LAS 分类字段其次用高度阈值。 mask np.ones(len(points), dtypebool) try: cls np.asarray(las.classification).astype(np.int8) mask (cls ! 2) except AttributeError: pass if min_height 0: mask (points[:, 2] min_height) return points[mask]高度阈值 2.0 米是经验值主要目的是排除地面低矮灌木。如果你的目标区域是灌木林或者幼林这个值要下调到 0.5 米左右否则会把真正的目标植株过滤掉。这段顺序上放在体素下采样之后可以少处理一些点体素尺寸不变的情况下不影响过滤效果。3.2 核心聚类把欧氏距离分割函数接入完整流程有了预处理基础核心分割流程可以串起来。输入是 LAS 文件路径输出是每个簇的点云和统计信息。流程顺序是读点云 → 体素下采样 → 过滤地面和低矮点 → 欧氏聚类 → 按簇拆分输出。def run_segmentation(input_path, eps1.2, min_points30, voxel_size0.3, min_height2.0): 执行完整的分割流程返回每个树簇的点云列表和统计信息。 las, points read_las_points(input_path) points voxel_downsample(points, voxel_size) points remove_ground_and_low(las, points, min_height) if len(points) 0: return {status: empty, tree_count: 0, tree_clusters: [], stats: []} labels euclidean_cluster(points, eps, min_points) n_clusters labels.max() 1 tree_clusters [] stats [] for cid in range(n_clusters): cluster_points points[labels cid] if cluster_points.shape[0] min_points: continue max_z cluster_points[:, 2].max() min_z cluster_points[:, 2].min() # 冠幅用 XY 平面上最远两点距离近似 xy cluster_points[:, :2] crown_diameter float(np.linalg.norm(xy.max(axis0) - xy.min(axis0))) tree_clusters.append(cluster_points) stats.append({ tree_id: cid 1, point_count: cluster_points.shape[0], max_height: round(float(max_z), 3), height_range: round(float(max_z - min_z), 3), crown_diameter: round(crown_diameter, 3), }) return {status: ok, tree_count: len(tree_clusters), tree_clusters: tree_clusters, stats: stats}这段流程里eps1.2、min_points30、voxel_size0.3是一组在机载激光雷达点云上常见的起点值。如果你用的是地基扫描仪点密度高很多eps 要降到 0.3 到 0.5 米min_points 可以提高到 50 以上。冠幅这里用 XY 坐标最大值和最小值的差来近似方向固定只适合粗略统计要更准确可以计算凸包或 PCA 主方向直径但对多数林业调查场景这个近似已经够用了。注意epsilon和min_points的匹配关系。如果 eps 设置过大两棵树的树冠直接合并成一个簇stats 里 tree_count 会明显少于实际株数点云密度大时这个现象尤其明显。如果 stats 里出现大量 point_count 只有十几的簇多半是 eps 太小把一棵树切碎了或者灌木没过滤干净。3.3 聚类结果后处理与文件输出写出 LAS 并生成汇总表分割完成只是第一步API 使用者需要拿到可交付的结果。常见做法是把每棵树的点云单独导出为一个 LAS 文件同时生成一个 CSV 汇总表包含每棵树的编号、点数、树高和冠幅。写 LAS 文件时有两个关键设置scale 和 offset。LAS 格式内部用 int32 存储坐标浮点坐标要乘以 scale 并减去 offset 才能转成整数如果 scale 太大或 offset 不设置输出文件的坐标会偏移甚至溢出。import csv def write_tree_las(output_path, points): 把单个树簇的点云写成 LAS 文件。 header laspy.LasHeader(point_format3, version1.2) header.scales np.array([0.001, 0.001, 0.001]) header.offsets points.min(axis0) las laspy.LasData(header) las.x points[:, 0] las.y points[:, 1] las.z points[:, 2] las.write(output_path) def write_summary_csv(stats, output_path): 把每棵树的统计信息写入 CSV。 with open(output_path, w, newline) as f: writer csv.DictWriter(f, fieldnames[tree_id, point_count, max_height, height_range, crown_diameter]) writer.writeheader() writer.writerows(stats)scale 设为 0.001 表示坐标精度到毫米满足林业测量的常规需求。offset 设为当前树簇点云的坐标最小值这样单棵树的局部坐标数值较小不会触发 int32 溢出。如果整个样地一起输出一个文件offset 应该用原始点云的全局最小值而不是单棵树的最小值。文件输出后整个最小实现已经闭环输入一个 LAS得到一组单木 LAS 文件和一张汇总表。到这一步你已经可以在本地脚本里完成单木分割。但脚本每次都要手动改参数、手动跑换成 API 服务后别人就能直接调用你的分割能力这也是标题里 API 的含义。4. 把分割流程封装成 APIFastAPI 接口、异步任务与参数设计4.1 FastAPI 接口设计输入输出与参数模型把分割流程封装成 HTTP 服务常见做法是用 FastAPI配合 python-multipart 处理文件上传。接口接收一个 LAS/LAZ 文件和四个参数返回一个 zip 包和一棵树的 CSV 汇总。下面是一个最小可用实现。from fastapi import FastAPI, File, UploadFile, Form, BackgroundTasks from fastapi.responses import FileResponse, JSONResponse import uuid, zipfile, shutil, os app FastAPI(titlePoint Cloud Tree Segmentation API) JOBS {} BASE_DIR workdir os.makedirs(BASE_DIR, exist_okTrue) app.post(/segment) async def create_segmentation_job( background_tasks: BackgroundTasks, file: UploadFile File(...), eps: float Form(1.2), min_points: int Form(30), voxel_size: float Form(0.3), min_height: float Form(2.0), ): 提交单木分割任务返回任务 ID。 job_id uuid.uuid4().hex input_path os.path.join(BASE_DIR, f{job_id}_input.las) with open(input_path, wb) as f: shutil.copyfileobj(file.file, f) JOBS[job_id] {status: pending, input_path: input_path} background_tasks.add_task(run_job, job_id, input_path, eps, min_points, voxel_size, min_height) return {job_id: job_id, status: pending}这里四个参数用Form(...)声明是因为文件用 multipart 上传时普通字段会以表单形式提交。eps和voxel_size是浮点数min_points是整数min_height用于过滤低矮灌木。上传文件先写到本地磁盘再处理避免大文件长时间占用内存。FastAPI 自动生成接口文档启动后访问/docs就能看到参数说明和测试界面。接口返回的job_id是后续查询任务状态的钥匙不要在前端直接等接口返回结果因为大点云文件分割可能要几十秒到几分钟HTTP 请求很容易超时。4.2 异步任务设计大点云文件不在请求线程里分割同步处理大文件是 API 设计里最容易翻车的地方。如果直接在 POST 请求里执行分割一个百万点的点云可能让请求挂起几分钟前端以为服务挂了。正确做法是后台任务加任务状态查询提交后立刻返回任务 ID前端轮询任务状态完成后从下载接口拿结果。def run_job(job_id, input_path, eps, min_points, voxel_size, min_height): 后台执行分割并打包结果。 try: JOBS[job_id] {status: running, input_path: input_path} result run_segmentation(input_path, eps, min_points, voxel_size, min_height) if result[status] ! ok: JOBS[job_id] {status: failed, error: no points remain after filtering} return output_dir os.path.join(BASE_DIR, job_id) os.makedirs(output_dir, exist_okTrue) zip_path os.path.join(output_dir, trees.zip) with zipfile.ZipFile(zip_path, w) as zf: for i, cluster_points in enumerate(result[tree_clusters]): las_name ftree_{i 1:03d}.las las_path os.path.join(output_dir, las_name) write_tree_las(las_path, cluster_points) zf.write(las_path, arcnamelas_name) write_summary_csv(result[stats], os.path.join(output_dir, summary.csv)) zf_path os.path.join(output_dir, trees.zip) with zipfile.ZipFile(zf_path, a) as zf: zf.write(os.path.join(output_dir, summary.csv), arcnamesummary.csv) JOBS[job_id] { status: done, tree_count: result[tree_count], zip_path: zf_path, input_path: input_path, } except Exception as exc: JOBS[job_id] {status: failed, error: str(exc)} app.get(/jobs/{job_id}) async def get_job(job_id: str): 查询任务状态。 job JOBS.get(job_id) if job is None: return JSONResponse({error: not found}, status_code404) if job[status] done: return {status: done, tree_count: job[tree_count], download_url: f/jobs/{job_id}/download} return {status: job[status]} app.get(/jobs/{job_id}/download) async def download_result(job_id: str): 下载分割结果 zip 包。 job JOBS.get(job_id) if job is None or job[status] ! done: return JSONResponse({error: not ready}, status_code404) return FileResponse(job[zip_path], filenameftrees_{job_id}.zip)这里用一个内存字典JOBS保存任务状态。单进程、单 worker 测试没问题但生产环境如果开了多个 worker每个进程的字典不共享任务状态会互相看不到。常见做法是把任务状态放到 Redis 等外部存储或者用 Celery/ARQ 这类任务队列。对个人工具和小团队内部系统内存字典足够对外提供服务强烈建议换外部存储。轮询间隔建议 2 到 3 秒查一次。任务状态从 pending 到 running 再到 done失败时会有 error 字段说明原因。整个异步模型让大文件分割变得可靠这是 API 化与本地脚本最大的区别。4.3 本地跑通与调用测试uvicorn 启动和 curl 验证本地验证整条链路先启动服务再用 curl 提交一个 LAS 文件。uvicorn 启动命令如下监听 8000 端口。uvicorn main:app --host 0.0.0.0 --port 8000启动后用 curl 提交分割任务参数通过-F传递文件用指定路径。curl -X POST http://127.0.0.1:8000/segment \ -F file./sample.las \ -F eps1.2 \ -F min_points30 \ -F voxel_size0.3 \ -F min_height2.0返回结果是一个 JSON里面包含job_id和statuspending。拿到 job_id 后轮询任务状态curl http://127.0.0.1:8000/jobs/job_id状态变为 done 后下载结果包curl -O -L http://127.0.0.1:8000/jobs/job_id/download下载得到的 zip 包里是每棵树的 LAS 文件和 summary.csv。用 CloudCompare 或 QGIS 打开可以快速检查分割结果是否符合预期。第一次跑通后建议用一个小样地数据几百棵树规模做基准测试记录不同参数下的 tree_count 和运行时间后续调参就有参照了。5. 单木分割避坑指南5 个踩坑点与参数调优方法5.1 参数相关的三个坑过分割、欠分割与灌木误判坑一体素下采样后 eps 没有跟着调一棵树被切成七八块。现象是分割结果里出现大量点数为几十的小簇同一棵树的树冠明显被拆散。原因是体素下采样把点间距拉大了原来的 eps 小于采样后点间距的 2 倍树冠内部的点无法连通成簇。解决办法是调参顺序固定为“先定 voxel_size再算下采样后的点间距最后定 eps”。点间距可以用np.min(pdist)估算也可以直接用经验公式eps 不小于体素尺寸的 2 倍机载数据通常 1.0 到 1.5 米。如果你的样地是密植人工林树间距小eps 要相应下调。坑二相邻树冠交错合并两个树冠输出成一个簇。现象是 stats 里出现点数量特别大的簇空间上明显覆盖两棵树的位置。原因是树冠枝条交错后两棵树之间的最近点距离小于 eps聚类把中间的空隙也连接起来了。单纯调小 eps 能缓解但容易导致过分割。我一般会加一个高度约束的二次处理把簇内点按高度排序只保留上部 70% 的点重新聚类树冠上部的间隙通常比下部明显更容易分开。如果你用的是机载数据还可以尝试对点云做高程归一化后单独在树冠层做一次聚类效果比直接全点聚类要好。坑三低矮灌木被当成树输出。现象是结果里有大量点数量在 20 到 80、高度只有 1 到 3 米的簇。原因是欧氏距离聚类只认距离不区分乔木和灌木。解决办法是组合使用高度阈值和点数阈值min_height过滤低矮点min_points过滤小簇输出阶段再加一道按树高筛选。比如设定只有 max_height 大于 5 米的簇才写入最终结果。矿业和林业项目里这道筛选能滤掉一半以上的假目标。5.2 数据与实现相关的两个坑点密度不均与 LAS 写入异常坑四机载点云密度不均匀稀疏区域的大树被拆分稠密区域的小树被合并。现象是同一片样地里有些树输出为多个簇有些树两棵合并。原因是扫描角度和遮挡导致树冠各部位点密度差异太大稀疏区域树冠内部点间距超过 eps稠密区域树间间隙又被点填平。解决思路是先统计整体点密度再做体素下采样把密度拉均匀。如果单棵树平均点数低于 50说明数据本身太稀欧氏聚类已经不是合适工具可以考虑改用 CHM 分水岭。对密度差异特别大的数据可以按区域分块设置不同 eps但这会让 API 参数设计复杂化我通常优先保证均匀下采样。坑五laspy 版本差异导致读取报错或者写出的 LAS 坐标偏移。现象是老脚本在新环境跑直接报AttributeError或者在 GIS 里打开输出文件位置偏离真实坐标。原因是 laspy 1.x 和 2.x 接口不兼容以及写文件时没有正确设置 scale 和 offset。解决办法是统一用 2.x 语法las.x / las.y / las.z写文件时显式设置header.scales (0.001, 0.001, 0.001)和header.offsets。另外 LAZ 压缩格式需要 laszip 支持没有安装对应扩展会读取失败遇到这种情况可以先转成 LAS 再处理。6. 验证分割效果与批量调参让结果可交付6.1 先做人工对照再用指标说话分割结果可不可靠不能只看 tree_count。我一般会先在点云可视化工具里加载原始点云和分割结果把每个簇用不同颜色渲染快速判断有没有明显欠分割或过分割。然后选一个 100 棵树左右的子区域人工标出树顶位置作为参考再对比算法输出。人工标注只需要在点云缩略图里点出明显的树冠最高点这一步不追求像素级准确足够评估算法效果。6.2 完整率、正确率与 F1 三个指标评估单木分割三个指标最常用。完整率等于正确分割的树的数量除以人工参考树的总数衡量有没有漏树正确率等于正确分割的树的数量除以算法输出的簇总数衡量有没有多分F1 是两者的调和平均综合反映分割质量。判定“一棵树分割正确”的标准是算法输出的某个簇与某棵人工参考树在水平位置上的偏差小于该树冠幅的一半。6.3 批量参数搜索脚本快速找到最优组合参数调优靠手试效率太低写一个简单网格搜索循环把所有组合跑一遍把 tree_count 和指标输出成表一眼就能看出参数趋势。import itertools param_grid { eps: [0.8, 1.0, 1.2, 1.5], min_points: [20, 30, 50], voxel_size: [0.2, 0.3, 0.5], } for eps, min_points, voxel_size in itertools.product( param_grid[eps], param_grid[min_points], param_grid[voxel_size] ): result run_segmentation(sample.las, eps, min_points, voxel_size, min_height2.0) print(eps, min_points, voxel_size, result[tree_count])跑完一遍你会直观看到 eps 从 0.8 升到 1.5 时 tree_count 的变化曲线多数数据会呈现“快速增长后趋稳再下降”的形态拐点附近就是较优参数。我现在的习惯是先跑一组小体素网格试探趋势再在拐点附近加密搜索。算法本身不复杂但参数调优的流程和验证方法才是 API 稳定输出的保障。希望这篇文章能帮你在自己的点云数据上少走几步弯路。本文还有配套的精品资源点击获取