简介：面向点云处理、林业遥感与三维重建方向的学习者，这份小型实践包给出基于PCL欧氏距离聚类完成单木分割的可调用实现，重点解决从无序点云中区分相邻树木、输出独立簇的入门与验证问题。压缩包共3个文件，约93KB，包含cpp源码、xyz示例点云和txt聚类结果，分别用于算法调用、输入测试与结果核对，体量轻便，适合快速查看接口组织与数据流转。已有815人学习，说明该方案在单木分割入门场景中具有一定参考价值。读者可结合源码理解欧氏距离聚类的建簇流程，通过示例数据观察邻域半径、簇点数阈值等参数对分割结果的影响，并利用结果文本对照检查每棵树的点集输出，为后续调参与二次开发提供可复用的起点。

1. 一片林地几百万个点，怎么把每棵树单独剥出来

拿到一份带三维坐标的林地激光点云，最直接的诉求就是数树：这片林子有多少棵、每棵多高、冠幅多大。传统路线是先插值出冠层高度模型，再在栅格上找局部极大值，阈值调起来全是玄学；这个项目换了条道——直接在原始点云上跑欧氏距离聚类，把空间里挨得够近的点归成同一棵，外面套一层 HTTP API，上传文件、传几个参数，拿回每棵树的点集和统计量。它适合手里已经有 LAS/LAZ 数据、又不想为单木分割去啃深度学习框架的开发者，也适合要把分割能力接进自己业务系统的后端同学。算法本身几十行就写完了，真正难的是尺度和参数怎么定。

2. 项目骨架与数据流：一个 LAS 文件进来，N 棵树出去

先把这个项目当成一个黑匣子看：左边是点云文件，右边是每棵树的编号、中心坐标、树高、冠幅和点集。中间夹着的所有东西，本质上是三件事——把文件读成 numpy 数组、把数组切成若干个簇、把簇翻译成业务字段。只要这三件事的边界划清楚了，代码放哪儿、怎么扩，都不会乱。实际拆下来，这个项目的模块切分基本就是这个思路，没有多余的抽象层。

2.1 目录结构：每个文件只干一件事

一份能跑的项目骨架大致长这样，命名可能不同，但职责划分基本一致：

tree-seg-api/ ├── app/ │ ├── main.py # FastAPI 应用入口，挂路由、挂中间件 │ ├── config.py # 参数默认值、上传目录、并发上限 │ ├── api/ │ │ ├── routes.py # /segment、/jobs/{id}、/health 三个接口 │ │ └── schemas.py # Pydantic 请求/响应模型 │ ├── core/ │ │ ├── io_las.py # LAS/LAZ 读取、坐标归一、写回 │ │ ├── preprocess.py # 体素下采样、地面点剔除、ROI 裁剪 │ │ ├── cluster.py # 欧氏距离聚类内核 │ │ └── metrics.py # 树高、冠幅、中心点、点云密度 │ ├── services/ │ │ └── pipeline.py # 把上面几个模块串成一次完整任务 │ └── storage/ │ ├── uploads/ # 原始文件落地 │ └── results/ # 结果 JSON / 分棵树 LAS ├── tests/ │ └── test_cluster.py # 用合成点云做单元测试 └── requirements.txt

这个结构里最关键的一条纪律：core/下面的四个文件不允许 importfastapi。聚类、读文件、算指标这些逻辑必须能在纯 Python 环境里被单元测试直接调用，否则后面想换个调度方式（比如改成命令行批处理）就得连锅端。这个做法并不高级，但我见过太多把业务逻辑写进路由函数的项目，最后连一次回归测试都跑不起来。

services/pipeline.py是唯一知道"顺序"的地方：读文件 → 预处理 → 聚类 → 后处理 → 算指标 → 写结果。它不关心 HTTP，也不关心点云具体怎么聚，只负责编排。想加一步"按坡度过滤"，改这一个文件就够了。

2.2 一次请求的完整数据流

从接口角度看，一次分割请求经过这么几个阶段：

1. 客户端POST /segment上传.las或.laz，同时带上eps、min_points、voxel_size三个参数。
2. 服务端校验后缀、大小、参数范围，把文件落到storage/uploads/，生成一个job_id。
3. 任务进入执行队列（小体量项目常见做法是BackgroundTasks或一个ThreadPoolExecutor，体量大就上独立 worker）。
4. io_las.read_points()把文件读成(N, 3)的float64数组，同时把intensity、classification这些属性暂存下来。
5. preprocess按voxel_size做体素下采样，把点数量降一个量级，顺便剔除明显的离群点。
6. cluster.euclidean_cluster()返回每个点的簇标签，-1表示噪声。
7. metrics按簇统计点数、最高点、中心 xy、冠幅（用 xy 方向的外接圆或分位数直径估算）。
8. 结果写storage/results/{job_id}.json，同时可选地按簇导出成多个小 LAS。
9. 客户端拿job_id去GET /jobs/{job_id}轮询，status从pending变成done后取result_url。

这条链路里，第 4 到第 7 步是纯计算，第 2、3、8 步是 IO 和调度。分开之后，最耗时的聚类就可以单独压测：同样的点云，改eps跑十次，看耗时和簇数量怎么变，不用起服务。

2.3 依赖与配置项

依赖不多，都是点云和 Web 方向的常规选择：

配置项集中在config.py，一般来说这几项要暴露成环境变量，因为不同数据源差别很大：

最后两项是血泪经验换来的。点云文件不像图片，一个 300MB 的 LAS 解压后可能占好几个 G 的内存，如果不设上限，两个请求同时进来就能把容器顶到 OOM。默认值不要照抄，eps必须按你自己的数据密度重新标定，后面避坑那章会细说。

3. 欧氏距离聚类内核：KDTree 半径搜索怎么把点连成树

算法名字听着唬人，其实一句话能说完：把所有点看成图上的节点，两个点的距离小于eps就连一条边，最后这个图里有几个连通分量，就有几棵树。真正决定成败的有三个细节——用什么数据结构找邻居、连通分量怎么扩、扩完之后怎么收拾残局。这三个细节处理不好，跑出来的结果不是把整片林子连成一坨，就是被砍成一地碎渣。

3.1 体素下采样：先把百万级点降到十万级

原始点云的密度往往远超需要。同一棵树冠上有几万个点，每个点都参与半径搜索，cKDTree也扛不住。常见做法是先做体素栅格下采样：把空间切成边长voxel_size的小立方体，每个立方体里的所有点用它们的质心代替。

import numpy as np def voxel_downsample(points: np.ndarray, voxel_size: float) -> np.ndarray: """体素下采样：每个体素内的点取质心，保留空间结构的同时降点数""" if voxel_size <= 0: return points # 关键：先减去最小值再取整，避免负坐标除法截断方向不一致 shifted = points - points.min(axis=0) voxel_idx = np.floor(shifted / voxel_size).astype(np.int64) # 把三维体素索引压成一个整数键，方便做唯一分组 keys = (voxel_idx[:, 0] * 73856093) ^ \ (voxel_idx[:, 1] * 19349663) ^ \ (voxel_idx[:, 2] * 83492791) # 按 key 排序后找分组边界，比 Python 字典快一个数量级 order = np.argsort(keys, kind="stable") sorted_keys = keys[order] boundaries = np.flatnonzero(np.diff(sorted_keys)) + 1 groups = np.split(order, boundaries) centroids = np.array([points[g].mean(axis=0) for g in groups]) return centroids

这段代码的逻辑是：先把坐标平移到非负区间，除以体素边长取整得到每个点所属的体素编号；然后用三个大质数异或的方式把三维编号压成一维键；接着排序、找相邻差异不为零的位置作为分组边界；最后每组求质心。参数上voxel_size一般取eps的 1/5 到 1/10，比如eps=0.4时用0.05。如果你的数据是单木级别的背包激光雷达，点密度本来就高，这个参数直接决定了后面聚类跑 3 秒还是 3 分钟。

提示：下采样会改变点云的绝对数量，min_points必须按降采样后的规模重新设。很多人把原始点云上标定的min_points=200直接搬过来，结果一棵树都留不下。

下采样有一个副作用要注意：质心会往点密集的方向偏，树冠边缘变稀疏。如果后续要算树高，最高点必须从原始点云里取，不能从降采样结果里取，否则树高会系统性偏低十几厘米。

3.2 半径搜索 + 广度优先：连通分量怎么扩出来

核心聚类函数就这一段。用cKDTree.query_ball_point找邻域，用队列做广度优先扩展，一个簇扩完再开下一个。

from collections import deque import numpy as np from scipy.spatial import cKDTree def euclidean_cluster(points: np.ndarray, eps: float = 0.4, min_points: int = 30): """欧氏距离聚类，返回 labels（-1 表示噪声）和簇数量""" if len(points) == 0: return np.array([], dtype=int), 0 tree = cKDTree(points) labels = np.full(len(points), -1, dtype=int) cluster_id = 0 for seed in range(len(points)): if labels[seed] != -1: continue # 从种子点开始做 BFS，把能连上的点全部收进同一个簇 queue = deque([seed]) labels[seed] = cluster_id members = 1 while queue: cur = queue.popleft() # 查询半径 eps 内的所有邻居索引 neighbors = tree.query_ball_point(points[cur], r=eps) for nb in neighbors: if labels[nb] == -1: labels[nb] = cluster_id queue.append(nb) members += 1 # 人数不够的簇整体回滚成噪声 if members < min_points: labels[labels == cluster_id] = -1 else: cluster_id += 1 return labels, cluster_id

逐段说明。cKDTree建树是一次性开销，复杂度大约O(N log N)，之后每次query_ball_point是O(log N + k)，k是邻域内点的数量。BFS 用deque而不是 list，popleft()是O(1)。labels数组一次开好，用-1表示未访问，避免额外维护一个 visited 集合。

参数上，eps是唯一真正影响分割结果的量。它的物理含义是"同一棵树内部相邻点之间允许的最大间距"，所以它和点云密度强相关：密度 1000 点/㎡ 的数据，点间距大约 3 厘米，eps取 0.3～0.5 米足够把树冠连起来又不会跨到邻树；密度 50 点/㎡ 的数据，点间距 15 厘米，eps就得放大到 0.8 米以上，否则一棵树内部自己就断成好几块。

min_points是第二道闸门，用来清掉灌木、杂草和飞点。它的取值依赖下采样后的密度，eps越大，一个簇能圈进来的点越多，min_points也得相应放大。

这里有个性能陷阱：labels[labels == cluster_id] = -1这一句在簇数量多的时候是O(N)的全数组扫描，如果碎簇很多，整体会退化。数据量上千万点时，我一般改成先记录members的索引列表，最后一次性回滚，代价是内存多占一点。

3.3 后处理：过滤、合并、算出业务字段

聚类出来的原始结果不能直接交给业务。三个问题必须处理：碎簇、过大的簇、以及每棵树到底怎么描述。

碎簇是指同一棵树被切成两三个小簇，通常发生在树冠稀疏或者有遮挡的位置。简单的处理规则是：两个簇的中心水平距离小于某个阈值（比如0.5 * eps + 0.3），且高度范围有重叠，就合并。这个规则不是万能的，但比不做强得多。

过大的簇通常是几棵树连在一起，尤其在密林里。如果簇的 xy 方向直径超过某个上限（比如 8 米），可以尝试按高度做二次切分，或者直接把整个簇标记为merged，让业务侧决定要不要人工复核。

业务字段的计算在metrics.py里：

import numpy as np def cluster_metrics(points: np.ndarray, labels: np.ndarray): """按簇统计树高、中心、冠幅，输出可直接入库的字典列表""" results = [] for cid in np.unique(labels): if cid < 0: continue pts = points[labels == cid] # 树高：该簇内最高点相对最低点的高差 height = float(pts[:, 2].max() - pts[:, 2].min()) # 中心：用 xy 中位数比均值更抗边缘飞点 cx, cy = float(np.median(pts[:, 0])), float(np.median(pts[:, 1])) # 冠幅：xy 方向 95 分位跨度，避免个别离群点把直径撑大 dx = np.percentile(pts[:, 0], 95) - np.percentile(pts[:, 0], 5) dy = np.percentile(pts[:, 1], 95) - np.percentile(pts[:, 1], 5) crown = float(max(dx, dy)) results.append({ "tree_id": int(cid), "center_x": round(cx, 3), "center_y": round(cy, 3), "height_m": round(height, 2), "crown_m": round(crown, 2), "point_count": int(len(pts)), }) return results

这里有两个刻意的选择。中心点用中位数而不是均值，因为树冠边缘常有零星飞点，均值会被拽偏。冠幅用 5%～95% 分位跨度而不是最大最小值，同样是抗离群点。树高用簇内高差而不是绝对高程，这样不需要知道当地高程基准，输出的是相对高度，业务侧如果要绝对高程再自己加地面点高程。

最后输出的 JSON 大概是这个样子，一棵树一条记录，前端拿来直接画点、画圈都行：

{ "job_id": "a1b2c3", "tree_count": 47, "params": {"eps": 0.4, "min_points": 30, "voxel_size": 0.05}, "trees": [ {"tree_id": 0, "center_x": 512345.12, "center_y": 3102345.67, "height_m": 12.4, "crown_m": 3.8, "point_count": 1287} ] }

4. API 层封装：把聚类算法包成能并发调用的 HTTP 接口

算法跑通只是第一步。真正给别的系统用的时候，暴露的是接口，不是函数。这一层要解决的问题很具体：文件怎么收、参数怎么校验、长任务怎么不阻塞、失败怎么让调用方看得懂。做得糙一点也能用，但一旦并发上来或者输入不合规，问题全在这层冒出来。

4.1 路由设计与请求参数校验

接口一共三个，不多：

参数校验用 Pydantic 做，比在函数里写一堆if强：

from fastapi import APIRouter, UploadFile, File, Form, HTTPException from pydantic import BaseModel, Field import uuid, shutil, os router = APIRouter() class SegmentParams(BaseModel): """分割参数，范围限制死，超出直接 422""" eps: float = Field(0.4, gt=0.05, lt=3.0, description="邻域半径，米") min_points: int = Field(30, gt=1, lt=10000, description="最小簇点数") voxel_size: float = Field(0.05, gt=0.0, lt=1.0, description="下采样体素，米") @router.post("/segment") async def segment( file: UploadFile = File(...), eps: float = Form(0.4), min_points: int = Form(30), voxel_size: float = Form(0.05), ): # 1) 后缀白名单，别指望 content-type，浏览器和脚本给的值不可靠 ext = os.path.splitext(file.filename)[-1].lower() if ext not in {".las", ".laz"}: raise HTTPException(400, detail=f"不支持的文件类型: {ext}") # 2) 参数走同一套校验逻辑，非法值在这里被拦下 try: params = SegmentParams(eps=eps, min_points=min_points, voxel_size=voxel_size) except ValueError as e: raise HTTPException(422, detail=str(e)) job_id = uuid.uuid4().hex[:12] dst = f"storage/uploads/{job_id}{ext}" # 3) 分块写盘，避免一次性 read() 把大文件全塞进内存 with open(dst, "wb") as f: while chunk := await file.read(1024 * 1024): f.write(chunk) # 4) 投递后台任务，立即返回 job_id submit_job(job_id, dst, params.dict()) return {"job_id": job_id, "status": "pending"}

这段代码有三个点是踩过坑才加的。第一，文件类型判断只看后缀不看content-type，因为不同客户端给的值差别太大，靠不住。第二，写盘用分块read，await file.read()不带参数会一次性把整个文件读进内存，512MB 的文件直接就是 512MB 内存。第三，参数校验走 Pydantic 模型，和内部调用共用同一份约束，避免命令行能跑、接口拒绝的不一致。

eps的上限设成 3.0 米是故意的。超过 3 米，在绝大多数林地数据里已经足以把相邻的树连成一片，与其让用户等十分钟拿到一坨结果，不如在入口就拦掉，并在错误信息里写清楚建议范围。

4.2 任务调度与状态机

长任务不能同步等着。小体量项目最省事的做法是一个固定大小的线程池，submit_job把任务丢进去，维护一个内存字典记录状态：

from concurrent.futures import ThreadPoolExecutor from threading import Lock from datetime import datetime _executor = ThreadPoolExecutor(max_workers=2) _jobs = {} _lock = Lock() def submit_job(job_id: str, path: str, params: dict): with _lock: _jobs[job_id] = { "status": "pending", "created_at": datetime.utcnow().isoformat(), "result": None, "error": None, } _executor.submit(_run_job, job_id, path, params) def _run_job(job_id: str, path: str, params: dict): with _lock: _jobs[job_id]["status"] = "running" try: result = run_pipeline(path, params) # 读文件→聚类→指标 _save_result(job_id, result) with _lock: _jobs[job_id].update(status="done", result=result) except MemoryError: with _lock: _jobs[job_id].update(status="failed", error="点云规模超出内存上限") except Exception as e: with _lock: _jobs[job_id].update(status="failed", error=f"{type(e).__name__}: {e}")

状态机只有四个值：pending、running、done、failed。别看简陋，它能覆盖 90% 的实际需求：调用方拿job_id轮询，done之后取结果，failed之后看error字段。

max_workers必须显式设小，别用默认值。线程池默认大小和 CPU 核数挂钩，看起来合理，但欧氏聚类是内存密集型不是 CPU 密集型，8 个任务同时跑，每个占 2GB，机器直接躺平。经验值是max_workers取min(2, 内存GB / 单任务峰值GB)，并且配合前面的MAX_POINTS做双重保险。

内存里存_jobs有个明显缺陷：进程重启就全没了。对于需要持久化的场景，把状态写进 SQLite 或 Redis，键是job_id，值是 JSON，改动量不大。结果本身写文件，_jobs里只存路径和状态，这样内存占用可控。

4.3 错误码与返回体约定

接口返回要让人看得懂，尤其是失败的时候。约定一套固定的错误结构：

{ "detail": { "code": "INVALID_PARAM", "message": "eps 超出允许范围 [0.05, 3.0]", "field": "eps" } }

常见的错误码就这几个：INVALID_FILE（后缀不对或文件损坏）、INVALID_PARAM（参数越界）、TOO_LARGE（超出点数或体积上限）、NOT_FOUND（job_id 不存在）、PROCESS_FAILED（计算异常，message里带原始异常类型）。

有一点容易忽略：run_pipeline里的异常必须被捕获并翻译成上面这套结构，不能直接往上抛。否则 FastAPI 默认会返回 500 加一段堆栈，调用方拿到一坨 HTML 完全不知道发生了什么。我在_run_job里统一包了 try/except，把异常类型和消息塞进error字段，GET /jobs/{id}再把它映射成对应的错误码。

5. 避坑与排查：调参调到怀疑人生的五个现场

这个算法最大的问题是"看起来能跑，结果不对"。它不会报错，只会给你一个荒唐的树数量。下面五条是我在实际数据上反复撞过的墙，按出现频率排。

5.1 整片林子被聚成一棵树

现象：跑完返回tree_count: 1，point_count等于全部点数。

原因：eps太大。欧氏距离聚类的连通性是有传递性的——A 连 B、B 连 C，A 和 C 就算隔着老远也在同一个簇里。只要eps超过了相邻两棵树之间最近点的间距，整个林分就会通过一条"点桥"串成一片。

解决：先算数据的平均最近邻距离，把eps设在它的 2～3 倍。快速估算的方法是对下采样后的点求每个点的最近邻距离，取中位数：

import numpy as np from scipy.spatial import cKDTree tree = cKDTree(points) dists, _ = tree.query(points, k=2) # 第 0 列是自己，距离 0 nn = dists[:, 1] print("最近邻中位数:", np.median(nn), "95分位:", np.percentile(nn, 95))

eps取中位数的 2～3 倍是个不错的起点，然后在这个值附近做小范围扫描，看树数量随eps变化的曲线，找那段平台区。

5.2 一棵大树被切成一堆碎块

现象：树数量明显偏多，每棵树的point_count都很小，可视化时看到一堆小色块挤在一起。

原因：eps太小，或者下采样太狠。树冠外缘的点本来就稀疏，如果voxel_size设得比树冠间隙还大，下采样后冠层会变成一个个孤岛，半径搜索根本连不起来。

解决：两个方向同时调。先把voxel_size降到eps的 1/10 以内，看碎块是否消失；如果还碎，说明eps确实偏小。另外检查min_points，它太小会把噪声当成树，太大又会把真正的小树干掉，一般设成下采样后单棵平均点数的 20%～30%。

5.3 下采样后树高系统性偏低

现象：算出来的树高比实测矮了 0.5～1.5 米，且偏差方向一致。

原因：用了下采样后的质心坐标算树高。体素质心必然落在原始点的凸包内部，树顶那个尖点被平均掉了。

解决：聚类和统计用两套点云。聚类在下采样数据上做，拿到标签后，用 KDTree 把原始点云按最近邻映射回簇标签，再从原始点算树高。多写十几行代码，换来的是厘米级的树高精度。

# labels_ds 是下采样点的标签，points_ds 是下采样坐标 tree = cKDTree(points_ds) _, idx = tree.query(points_raw) # 原始点找最近的下采样点 labels_raw = labels_ds[idx] # 继承标签

注意这一步前提是voxel_size不能过大，否则原始点和下采样点对不上号，一棵树的点可能被映射到邻树上。

5.4 服务跑两个请求就 OOM

现象：单请求正常，并发两个就容器被杀，日志里连异常都没有。

原因：cKDTree建树本身要占内存（大约是点数的若干倍），加上原始数组、下采样数组、标签数组、query_ball_point返回的邻居列表，一个百万级点的任务峰值可能到 GB 级别。线程池默认大小又偏大，几个任务一起跑必然崩。

解决：三重限制同时上——MAX_POINTS在入口拒绝超大文件；线程池max_workers手动设成 1 或 2；_run_job结束后显式del大数组并gc.collect()。如果业务量确实大，别硬扛，改成外部 worker 进程 + 消息队列，每个 worker 限制内存配额。

5.5 换了一批数据，参数全废

现象：同一个服务，A 区的数据调好的参数，到 B 区跑出来结果完全不对。

原因：eps是绝对距离，和点云密度、扫描方式、树木种类都绑定。机载数据和背包数据密度差两个数量级，针叶林和阔叶林的冠层形态也不一样，参数当然不能通用。

解决：把参数按"数据源"分组管理，配置里存一组 profile，请求时指定 profile 名称而不是裸传eps。同时提供一个/calibrate接口，输入样本点云，自动扫eps的几个候选值，返回树数量曲线，让人工挑。这一步看着麻烦，但比每次上线都靠猜要省事得多。

6. 二次开发：换算法、做标定、验结果

拿到代码之后最常做的一类改动，是把欧氏聚类换成能利用形态信息的算法。欧氏距离只看"近不近"，不看"像不像一棵树"，所以在密度不均的数据上天生吃亏。区域生长是个平滑的替代方案：从局部最高点当种子，按高度差和水平距离双重阈值往外扩，扩到不平缓的位置就停。

def region_growing(points, labels_init, max_dz=1.5, max_dxy=2.0): """在欧氏聚类结果上做区域生长修正： 把标签不一致但高度连续、水平邻近的点并到同一棵树""" from scipy.spatial import cKDTree tree = cKDTree(points[:, :2]) # 只看平面距离 merged = labels_init.copy() for i, p in enumerate(points): cid = merged[i] if cid < 0: continue for j in tree.query_ball_point(p[:2], r=max_dxy): # 相邻但标签不同，且高差在允许范围内，视为同一棵树 if merged[j] != cid and abs(points[j, 2] - p[2]) < max_dz: merged[merged == merged[j]] = cid return merged

这段代码的用途是在原有聚类结果上打补丁，不做全量重算，改动风险小。逻辑是用平面 KDTree 找水平邻近点，如果两个点标签不同但高差小于max_dz（默认 1.5 米，相当于一层枝叶的高度跨度），就把它们并到同一个标签。max_dxy控制合并的激进程度，设大了会把邻树吞并，设小了没效果。

改造之后必须验证，不能只看树数量对不对。我的习惯是做两层校验。第一层是数量校验：拿一块有人工标注样地（或者航拍影像目视解译）的数据，比对分割结果和真值的树数量、位置匹配率。位置匹配用简单的距离阈值匹配，分割中心到真值中心 2 米内就算命中。第二层是稳定性校验：同一块数据，把点云随机抽稀 80%，重跑一遍，看树数量和每棵树的高度变化幅度。变化大的话，说明参数贴着临界值，不可靠。

这套流程我后来做每个点云项目都强制走一遍，哪怕只是改了一个默认参数。有次图省事跳过了抽稀校验，上线之后换了个扫描角度，同一片林子的树数量从 47 变成 63，回头排查才发现是min_points卡在了临界值上。从那以后我每次改聚类参数都强制走一遍「抽稀 + 重跑 + 比对」这三步，宁可多花二十分钟，也不想到线上去猜。希望帮到你。

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