简介面向点云处理、林业遥感与三维重建方向的学习者这份小型实践包给出基于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 一次请求的完整数据流从接口角度看一次分割请求经过这么几个阶段客户端POST /segment上传.las或.laz同时带上eps、min_points、voxel_size三个参数。服务端校验后缀、大小、参数范围把文件落到storage/uploads/生成一个job_id。任务进入执行队列小体量项目常见做法是BackgroundTasks或一个ThreadPoolExecutor体量大就上独立 worker。io_las.read_points()把文件读成(N, 3)的float64数组同时把intensity、classification这些属性暂存下来。preprocess按voxel_size做体素下采样把点数量降一个量级顺便剔除明显的离群点。cluster.euclidean_cluster()返回每个点的簇标签-1表示噪声。metrics按簇统计点数、最高点、中心 xy、冠幅用 xy 方向的外接圆或分位数直径估算。结果写storage/results/{job_id}.json同时可选地按簇导出成多个小 LAS。客户端拿job_id去GET /jobs/{job_id}轮询status从pending变成done后取result_url。这条链路里第 4 到第 7 步是纯计算第 2、3、8 步是 IO 和调度。分开之后最耗时的聚类就可以单独压测同样的点云改eps跑十次看耗时和簇数量怎么变不用起服务。2.3 依赖与配置项依赖不多都是点云和 Web 方向的常规选择依赖用途备注laspy读写 LAS/LAZ读 LAZ 需要额外的解压后端numpy数组运算全程主数据结构scipycKDTree半径搜索聚类性能的关键open3d可选可视化、下采样纯算法可以不用fastapiuvicornHTTP 服务自带 OpenAPI 文档pydantic参数校验请求体模型配置项集中在config.py一般来说这几项要暴露成环境变量因为不同数据源差别很大配置项默认值说明DEFAULT_EPS0.4聚类的邻域半径单位米DEFAULT_MIN_POINTS30一个簇最少多少个点小于它当噪声丢掉DEFAULT_VOXEL_SIZE0.05下采样体素边长单位米MAX_UPLOAD_MB512单文件上限MAX_POINTS20000000硬上限超了直接拒绝防止内存打爆WORKER_CONCURRENCY2同时跑几个分割任务最后两项是血泪经验换来的。点云文件不像图片一个 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(axis0) 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, kindstable) sorted_keys keys[order] boundaries np.flatnonzero(np.diff(sorted_keys)) 1 groups np.split(order, boundaries) centroids np.array([points[g].mean(axis0) for g in groups]) return centroids这段代码的逻辑是先把坐标平移到非负区间除以体素边长取整得到每个点所属的体素编号然后用三个大质数异或的方式把三维编号压成一维键接着排序、找相邻差异不为零的位置作为分组边界最后每组求质心。参数上voxel_size一般取eps的 1/5 到 1/10比如eps0.4时用0.05。如果你的数据是单木级别的背包激光雷达点密度本来就高这个参数直接决定了后面聚类跑 3 秒还是 3 分钟。提示下采样会改变点云的绝对数量min_points必须按降采样后的规模重新设。很多人把原始点云上标定的min_points200直接搬过来结果一棵树都留不下。下采样有一个副作用要注意质心会往点密集的方向偏树冠边缘变稀疏。如果后续要算树高最高点必须从原始点云里取不能从降采样结果里取否则树高会系统性偏低十几厘米。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([], dtypeint), 0 tree cKDTree(points) labels np.full(len(points), -1, dtypeint) 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], reps) 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而不是 listpopleft()是O(1)。labels数组一次开好用-1表示未访问避免额外维护一个 visited 集合。参数上eps是唯一真正影响分割结果的量。它的物理含义是同一棵树内部相邻点之间允许的最大间距所以它和点云密度强相关密度 1000 点/㎡ 的数据点间距大约 3 厘米eps取 0.30.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 路由设计与请求参数校验接口一共三个不多方法路径作用POST/segment上传点云提交分割任务返回job_idGET/jobs/{job_id}查询任务状态和结果GET/health健康检查给网关用参数校验用 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, gt0.05, lt3.0, description邻域半径米) min_points: int Field(30, gt1, lt10000, description最小簇点数) voxel_size: float Field(0.05, gt0.0, lt1.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, detailf不支持的文件类型: {ext}) # 2) 参数走同一套校验逻辑非法值在这里被拦下 try: params SegmentParams(epseps, min_pointsmin_points, voxel_sizevoxel_size) except ValueError as e: raise HTTPException(422, detailstr(e)) job_id uuid.uuid4().hex[:12] dst fstorage/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因为不同客户端给的值差别太大靠不住。第二写盘用分块readawait 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_workers2) _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(statusdone, resultresult) except MemoryError: with _lock: _jobs[job_id].update(statusfailed, error点云规模超出内存上限) except Exception as e: with _lock: _jobs[job_id].update(statusfailed, errorf{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_FOUNDjob_id 不存在、PROCESS_FAILED计算异常message里带原始异常类型。有一点容易忽略run_pipeline里的异常必须被捕获并翻译成上面这套结构不能直接往上抛。否则 FastAPI 默认会返回 500 加一段堆栈调用方拿到一坨 HTML 完全不知道发生了什么。我在_run_job里统一包了 try/except把异常类型和消息塞进error字段GET /jobs/{id}再把它映射成对应的错误码。5. 避坑与排查调参调到怀疑人生的五个现场这个算法最大的问题是看起来能跑结果不对。它不会报错只会给你一个荒唐的树数量。下面五条是我在实际数据上反复撞过的墙按出现频率排。5.1 整片林子被聚成一棵树现象跑完返回tree_count: 1point_count等于全部点数。原因eps太大。欧氏距离聚类的连通性是有传递性的——A 连 B、B 连 CA 和 C 就算隔着老远也在同一个簇里。只要eps超过了相邻两棵树之间最近点的间距整个林分就会通过一条点桥串成一片。解决先算数据的平均最近邻距离把eps设在它的 23 倍。快速估算的方法是对下采样后的点求每个点的最近邻距离取中位数import numpy as np from scipy.spatial import cKDTree tree cKDTree(points) dists, _ tree.query(points, k2) # 第 0 列是自己距离 0 nn dists[:, 1] print(最近邻中位数:, np.median(nn), 95分位:, np.percentile(nn, 95))eps取中位数的 23 倍是个不错的起点然后在这个值附近做小范围扫描看树数量随eps变化的曲线找那段平台区。5.2 一棵大树被切成一堆碎块现象树数量明显偏多每棵树的point_count都很小可视化时看到一堆小色块挤在一起。原因eps太小或者下采样太狠。树冠外缘的点本来就稀疏如果voxel_size设得比树冠间隙还大下采样后冠层会变成一个个孤岛半径搜索根本连不起来。解决两个方向同时调。先把voxel_size降到eps的 1/10 以内看碎块是否消失如果还碎说明eps确实偏小。另外检查min_points它太小会把噪声当成树太大又会把真正的小树干掉一般设成下采样后单棵平均点数的 20%30%。5.3 下采样后树高系统性偏低现象算出来的树高比实测矮了 0.51.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_dz1.5, max_dxy2.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], rmax_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%重跑一遍看树数量和每棵树的高度变化幅度。变化大的话说明参数贴着临界值不可靠。校验项方法合格线经验值树数量偏差与样地真值比对偏差 15% 以内位置匹配率2 米阈值最近邻匹配85% 以上树高一致性与实测树高回归R² 0.85 以上抽稀稳定性80% 抽稀重跑树数量波动 10% 以内这套流程我后来做每个点云项目都强制走一遍哪怕只是改了一个默认参数。有次图省事跳过了抽稀校验上线之后换了个扫描角度同一片林子的树数量从 47 变成 63回头排查才发现是min_points卡在了临界值上。从那以后我每次改聚类参数都强制走一遍「抽稀 重跑 比对」这三步宁可多花二十分钟也不想到线上去猜。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/11/4 10:18:38 · 阅读 1.2k 激光点云单木分割：欧氏聚类与FastAPI API工程实践
简介面向点云处理、林业遥感与三维重建方向的学习者这份小型实践包给出基于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 一次请求的完整数据流从接口角度看一次分割请求经过这么几个阶段客户端POST /segment上传.las或.laz同时带上eps、min_points、voxel_size三个参数。服务端校验后缀、大小、参数范围把文件落到storage/uploads/生成一个job_id。任务进入执行队列小体量项目常见做法是BackgroundTasks或一个ThreadPoolExecutor体量大就上独立 worker。io_las.read_points()把文件读成(N, 3)的float64数组同时把intensity、classification这些属性暂存下来。preprocess按voxel_size做体素下采样把点数量降一个量级顺便剔除明显的离群点。cluster.euclidean_cluster()返回每个点的簇标签-1表示噪声。metrics按簇统计点数、最高点、中心 xy、冠幅用 xy 方向的外接圆或分位数直径估算。结果写storage/results/{job_id}.json同时可选地按簇导出成多个小 LAS。客户端拿job_id去GET /jobs/{job_id}轮询status从pending变成done后取result_url。这条链路里第 4 到第 7 步是纯计算第 2、3、8 步是 IO 和调度。分开之后最耗时的聚类就可以单独压测同样的点云改eps跑十次看耗时和簇数量怎么变不用起服务。2.3 依赖与配置项依赖不多都是点云和 Web 方向的常规选择依赖用途备注laspy读写 LAS/LAZ读 LAZ 需要额外的解压后端numpy数组运算全程主数据结构scipycKDTree半径搜索聚类性能的关键open3d可选可视化、下采样纯算法可以不用fastapiuvicornHTTP 服务自带 OpenAPI 文档pydantic参数校验请求体模型配置项集中在config.py一般来说这几项要暴露成环境变量因为不同数据源差别很大配置项默认值说明DEFAULT_EPS0.4聚类的邻域半径单位米DEFAULT_MIN_POINTS30一个簇最少多少个点小于它当噪声丢掉DEFAULT_VOXEL_SIZE0.05下采样体素边长单位米MAX_UPLOAD_MB512单文件上限MAX_POINTS20000000硬上限超了直接拒绝防止内存打爆WORKER_CONCURRENCY2同时跑几个分割任务最后两项是血泪经验换来的。点云文件不像图片一个 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(axis0) 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, kindstable) sorted_keys keys[order] boundaries np.flatnonzero(np.diff(sorted_keys)) 1 groups np.split(order, boundaries) centroids np.array([points[g].mean(axis0) for g in groups]) return centroids这段代码的逻辑是先把坐标平移到非负区间除以体素边长取整得到每个点所属的体素编号然后用三个大质数异或的方式把三维编号压成一维键接着排序、找相邻差异不为零的位置作为分组边界最后每组求质心。参数上voxel_size一般取eps的 1/5 到 1/10比如eps0.4时用0.05。如果你的数据是单木级别的背包激光雷达点密度本来就高这个参数直接决定了后面聚类跑 3 秒还是 3 分钟。提示下采样会改变点云的绝对数量min_points必须按降采样后的规模重新设。很多人把原始点云上标定的min_points200直接搬过来结果一棵树都留不下。下采样有一个副作用要注意质心会往点密集的方向偏树冠边缘变稀疏。如果后续要算树高最高点必须从原始点云里取不能从降采样结果里取否则树高会系统性偏低十几厘米。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([], dtypeint), 0 tree cKDTree(points) labels np.full(len(points), -1, dtypeint) 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], reps) 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而不是 listpopleft()是O(1)。labels数组一次开好用-1表示未访问避免额外维护一个 visited 集合。参数上eps是唯一真正影响分割结果的量。它的物理含义是同一棵树内部相邻点之间允许的最大间距所以它和点云密度强相关密度 1000 点/㎡ 的数据点间距大约 3 厘米eps取 0.30.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 路由设计与请求参数校验接口一共三个不多方法路径作用POST/segment上传点云提交分割任务返回job_idGET/jobs/{job_id}查询任务状态和结果GET/health健康检查给网关用参数校验用 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, gt0.05, lt3.0, description邻域半径米) min_points: int Field(30, gt1, lt10000, description最小簇点数) voxel_size: float Field(0.05, gt0.0, lt1.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, detailf不支持的文件类型: {ext}) # 2) 参数走同一套校验逻辑非法值在这里被拦下 try: params SegmentParams(epseps, min_pointsmin_points, voxel_sizevoxel_size) except ValueError as e: raise HTTPException(422, detailstr(e)) job_id uuid.uuid4().hex[:12] dst fstorage/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因为不同客户端给的值差别太大靠不住。第二写盘用分块readawait 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_workers2) _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(statusdone, resultresult) except MemoryError: with _lock: _jobs[job_id].update(statusfailed, error点云规模超出内存上限) except Exception as e: with _lock: _jobs[job_id].update(statusfailed, errorf{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_FOUNDjob_id 不存在、PROCESS_FAILED计算异常message里带原始异常类型。有一点容易忽略run_pipeline里的异常必须被捕获并翻译成上面这套结构不能直接往上抛。否则 FastAPI 默认会返回 500 加一段堆栈调用方拿到一坨 HTML 完全不知道发生了什么。我在_run_job里统一包了 try/except把异常类型和消息塞进error字段GET /jobs/{id}再把它映射成对应的错误码。5. 避坑与排查调参调到怀疑人生的五个现场这个算法最大的问题是看起来能跑结果不对。它不会报错只会给你一个荒唐的树数量。下面五条是我在实际数据上反复撞过的墙按出现频率排。5.1 整片林子被聚成一棵树现象跑完返回tree_count: 1point_count等于全部点数。原因eps太大。欧氏距离聚类的连通性是有传递性的——A 连 B、B 连 CA 和 C 就算隔着老远也在同一个簇里。只要eps超过了相邻两棵树之间最近点的间距整个林分就会通过一条点桥串成一片。解决先算数据的平均最近邻距离把eps设在它的 23 倍。快速估算的方法是对下采样后的点求每个点的最近邻距离取中位数import numpy as np from scipy.spatial import cKDTree tree cKDTree(points) dists, _ tree.query(points, k2) # 第 0 列是自己距离 0 nn dists[:, 1] print(最近邻中位数:, np.median(nn), 95分位:, np.percentile(nn, 95))eps取中位数的 23 倍是个不错的起点然后在这个值附近做小范围扫描看树数量随eps变化的曲线找那段平台区。5.2 一棵大树被切成一堆碎块现象树数量明显偏多每棵树的point_count都很小可视化时看到一堆小色块挤在一起。原因eps太小或者下采样太狠。树冠外缘的点本来就稀疏如果voxel_size设得比树冠间隙还大下采样后冠层会变成一个个孤岛半径搜索根本连不起来。解决两个方向同时调。先把voxel_size降到eps的 1/10 以内看碎块是否消失如果还碎说明eps确实偏小。另外检查min_points它太小会把噪声当成树太大又会把真正的小树干掉一般设成下采样后单棵平均点数的 20%30%。5.3 下采样后树高系统性偏低现象算出来的树高比实测矮了 0.51.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_dz1.5, max_dxy2.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], rmax_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%重跑一遍看树数量和每棵树的高度变化幅度。变化大的话说明参数贴着临界值不可靠。校验项方法合格线经验值树数量偏差与样地真值比对偏差 15% 以内位置匹配率2 米阈值最近邻匹配85% 以上树高一致性与实测树高回归R² 0.85 以上抽稀稳定性80% 抽稀重跑树数量波动 10% 以内这套流程我后来做每个点云项目都强制走一遍哪怕只是改了一个默认参数。有次图省事跳过了抽稀校验上线之后换了个扫描角度同一片林子的树数量从 47 变成 63回头排查才发现是min_points卡在了临界值上。从那以后我每次改聚类参数都强制走一遍「抽稀 重跑 比对」这三步宁可多花二十分钟也不想到线上去猜。希望帮到你。本文还有配套的精品资源点击获取