带Classification的LAS跑通PointNet++点云语义分割

简介面向需要将 PointNet 或 PointNet 落地到自有 LAS 点云的研究者与工程人员这份资料聚焦带 Classification 属性的 LAS 点云如何组织训练与预测基于 yanx27 的 PyTorch 实现整理可在两种网络结构间按需选择适合已有一定深度学习基础、希望跑通自定义点云分类或分割流程的读者。压缩包共 19 个文件、约 26KB主体为 15 个 Python 脚本覆盖网络模型、数据加载、LAS 数据处理及训练测试入口另附一个已分类 LAS 文件用于训练与验证、一个 LAS 文件用于预测以及一份说明代码结构的 docx 文档便于对照目录快速定位模块。已有 896 人学习。读者可据此理解从原始 LAS 到数据集构建、再到模型训练与推理的完整链路并借助现成脚本减少环境配置与数据接口改造中的反复试错快速形成可复现的自定义点云实验基线。1. 从一份带 Classification 的 LAS 说起PointNet 语义分割落地要跨过哪几道坎手上拿到一份 LAS每条记录里都带着 Classification 字段看起来标签已经现成直接丢进 PointNet / PointNet 就能跑点云语义分割。真动手就会发现中间隔着好几道坎LAS 里的分类码是测绘体系的编码地面是 2、建筑是 6、植被是 3/4/5不是从 0 开始的连续类别一整幅点云动辄几百万点显存根本塞不下点云没有网格结构分块、采样、归一化任何一步训练和推理不一致结果就会翻车。这篇笔记按 pytorch 点云语义分割的实际流程走一遍从 LAS 读 Classification、做类别映射和分块采样、改 PointNet 的分割头、训练评估到最后把预测结果写回 LAS。适合手里有带分类码的机载或车载点云、想拿 PointNet 系列做一个能跑通的 baseline 的从业者。2. 拆开 LAS 的 Classification用 laspy 把点云读成可训练的数组2.1 LAS 的坐标、尺度与 Classification 字段语义LAS 文件分成 header 和 point records 两块。header 里最该盯的是版本号、point format、scale/offset 和点数point records 里才是每个点的坐标和属性。坐标在文件里存的是整数读出来要按x X_record * scale_x offset_x还原成浮点。很多人第一次读 LAS 觉得坐标“怎么是几百万”就是因为只看了整数部分没管 offset。laspy 的las.x已经替你做了这层换算直接拿来用就行但 header 里的 scales 和 offsets 建议打印出来看一眼机载数据的 scale 常见是 0.001offset 一般是大数丢了它就等于把整幅点云平移到了别处。Classification 在 LAS 1.0 到 1.3 里就是一个 uint8 字段取值遵循 ASPRS 的分类码约定LAS 1.4 的 point format 6 以上把分类标志位synthetic、key-point、withheld、overlap拆成了独立字节Classification 本身仍然是 uint8。这意味着你拿到的分类码大概率是测绘软件写进去的语义码而不是神经网络要的 0 起始连续标签。下表是机载点云里最常出现的几个码做映射时基本绕不开LAS 分类码常见含义处理建议0未分类一般设为 ignore别当负样本硬训1未指定视数据质量决定 ignore 或并入某类2地面保留通常占比最高3 / 4 / 5低/中/高植被按任务合并成一类或拆成三类6建筑保留9水体保留或 ignore看任务是否需要7 / 17 / 18噪声、桥梁等建议 ignore噪声点会污染梯度提示分类码 0 和 1 在多数开源数据集里都是“没标”的意思把它们当正常类别训模型会学到一个语义混杂的类mIoU 看着不掉但可视化会很难看。2.2 先统计再动手读取与类别分布不要一上来就写转换脚本。第一件事是把文件读进来看点数和类别直方图这一步能提前暴露 80% 的坑。import numpy as np import laspy # 读入单幅 LASlaspy 新版把子字段直接挂在对象上 las laspy.read(tile_01.las) hdr las.header print(version:, hdr.version, point_format:, hdr.point_format.id) print(point_count:, hdr.point_count) print(scales:, hdr.scales, offsets:, hdr.offsets) # x/y/z 已经按 scale/offset 还原为浮点直接堆成 N×3 xyz np.stack([las.x, las.y, las.z], axis1).astype(np.float32) cls np.asarray(las.classification).astype(np.uint8) print(xyz range:, xyz.min(0), xyz.max(0)) uniq, cnt np.unique(cls, return_countsTrue) for u, c in zip(uniq, cnt): print(fclass {u:3d}: {c:10d} {c / len(cls) * 100:6.2f}%)逻辑很直白先把 header 里的版本、点格式、缩放参数打出来再用las.classification拿到分类数组做直方图。参数上唯一要注意的是 dtypelas.classification读出来是 uint8转 numpy 时别再降精度xyz必须转 float32float64 会让后面的 Dataset 内存翻倍。跑完这段你会得到两个关键信息类别码都有哪些、地面类占比多少。地面占到 70% 以上是常态这直接决定了后面损失函数和采样策略怎么设。2.3 类别码映射用 LUT 把测绘编码压成 0..N-1网络输出的是 0 到 N-1 的连续类别LAS 给的是稀疏的分类码中间必须有一张映射表。我一般写成字典再展开成 256 长度的查找表比np.vectorize快一个数量级而且反查也方便。import numpy as np # 左边是 LAS 分类码右边是网络类别 id没列出的码统一落到 ignore CLASS_MAP {2: 0, 6: 1, 3: 2, 4: 2, 5: 2, 9: 3} IGNORE 255 lut np.full(256, IGNORE, dtypenp.uint8) for src, dst in CLASS_MAP.items(): lut[src] dst label lut[cls] # 映射后的标签取值 0..3 或 255 print(np.unique(label, return_countsTrue)) # 反查表把网络预测的类别 id 还原成 LAS 分类码 inv_lut np.zeros(256, dtypenp.uint8) for src, dst in CLASS_MAP.items(): inv_lut[dst] src这段的核心是lut和inv_lut两张表。lut用分类码直接索引256 项覆盖 uint8 全空间没配的码默认 255天然就是ignore_index不需要额外写 if。inv_lut是写回 LAS 时的后悔药——网络输出 1你得知道它对应 LAS 的 6 而不是原样写 1。参数上要留意一个细节多个源码映射到同一个目标类3/4/5 都映射到 2时反查只能还原成其中一个码写回前想清楚这是否可接受如果下游要区分植被高度就得在映射阶段就拆开。2.4 落盘npy 存坐标、标签和原始索引转换脚本跑一次要几分钟到几十分钟别每次都重读 LAS。我习惯把一幅图拆成三个数组存坐标、标签、原始点序号。第三个最容易被忽略但写回预测结果时全靠它。import numpy as np import laspy las laspy.read(tile_01.las) xyz np.stack([las.x, las.y, las.z], axis1).astype(np.float32) cls np.asarray(las.classification).astype(np.uint8) lut np.full(256, 255, dtypenp.uint8) for src, dst in {2: 0, 6: 1, 3: 2, 4: 2, 5: 2, 9: 3}.items(): lut[src] dst np.save(tile_01_xyz.npy, xyz) np.save(tile_01_label.npy, lut[cls]) np.save(tile_01_idx.npy, np.arange(len(xyz), dtypenp.int64))idx就是 0 到 N-1 的自然序列看着没用但分块采样会打乱点序推理阶段要靠它把每个点的预测投回原始位置。坐标用 float32、标签用 uint8、索引用 int64是内存和精度之间比较划算的组合。文件大到内存放不下时换成np.memmap或者按 tile 分文件存不要试图一次性读完整测区。3. 把整幅点云切成训练样本分块、采样与标签对齐3.1 为什么必须分块显存、感受野与边界效应PointNet 的 set abstraction 层做的是球查询加最远点采样计算量随点数呈超线性增长。一幅机载点云几百万点直接喂进去不是显存爆就是跑一晚上一个 epoch。分块是刚需但块大小不是随便定的块太小单块里可能只有地面模型学不到建筑和植被的上下文块太大显存吃不住而且最远点采样在超大点集上会漏掉小目标。工程上常用的起点是水平块 2 到 5 米、每块采 4096 个点机载数据可以放宽到 8192车载或室内数据往 1 米、2048 点收。另一个容易被忽视的是边界效应。滑窗切块时落在块边缘的建筑会被切断模型在边缘学到的形状是残缺的。训练阶段我一般用不重叠切块stride 等于 block推理阶段才用 50% 重叠加投票两者分开处理别用同一套参数硬套。3.2 滑窗分块与点数补齐的实现分块函数要解决三件事按平面网格取点、点数不足时补齐、点数过多时降采样。补齐用有放回采样比补零好补零会让网络学到一堆位于原点的假点。import numpy as np def make_blocks(xyz, label, idx, block2.0, stride2.0, num_point4096, min_pts512, seed0): rng np.random.default_rng(seed) xmin, ymin xyz[:, 0].min(), xyz[:, 1].min() xmax, ymax xyz[:, 0].max(), xyz[:, 1].max() xs np.arange(xmin, xmax, stride) ys np.arange(ymin, ymax, stride) blocks [] for bx in xs: for by in ys: m ((xyz[:, 0] bx) (xyz[:, 0] bx block) (xyz[:, 1] by) (xyz[:, 1] by block)) n int(m.sum()) if n min_pts: # 点太少放弃这一块 continue sel np.nonzero(m)[0] if n num_point: sel rng.choice(sel, num_point, replaceFalse) # 降采样 else: extra rng.choice(sel, num_point - n, replaceTrue) # 有放回补齐 sel np.concatenate([sel, extra]) blocks.append((xyz[sel], label[sel], idx[sel])) return blocks参数含义逐个说清楚block是水平块边长单位跟坐标系一致投影坐标下就是米stride是滑窗步长训练时设成和block相等做不重叠切块推理时设成block/2换取重叠num_point决定单样本点数直接对应显存min_pts是块的最小点数阈值低于它就跳过避免大量只有几十个点的垃圾块进入训练集。seed固定住能保证每次转换结果一致调参排查时很省事。注意这里只按 x/y 切块z 方向不切。机载数据 z 跨度大按三维切会产生大量空块但如果你的数据是高层建筑立面加地面块内高差可能超过 100 米归一化时要留意。3.3 归一化绝对坐标是训练里的隐形炸弹投影坐标的数值动辄几十万直接送进网络第一层 MLP 的激活值就炸了。归一化必须做而且训练和推理必须用同一个函数这是血泪经验。import numpy as np def normalize(xyz): # 先平移到块质心再缩放到单位球 center xyz.mean(axis0, keepdimsTrue) p xyz - center scale np.max(np.linalg.norm(p, axis1)) return (p / max(scale, 1e-6)).astype(np.float32)这段做的是“去质心 球归一化”。先用块的质心做平移消掉绝对坐标再除以块内点到质心的最大距离把所有点压进单位球。为什么不逐轴 min-max因为那会改变点云的几何比例一个狭长的建筑块会被拉成方块法向和曲率全变了。为什么不除以固定常数比如 50 米因为不同块的尺度差异很大固定常数在稀疏块上会把点压成一团。有个边界情况必须提如果块里只有地面点近似共面scale可能很小除完会放大噪声点的扰动。max(scale, 1e-6)只是防除零真正的解法是在分块阶段就把纯地面块按比例丢弃或者给scale设一个下限。3.4 训练/验证划分按区域切不要按块随机切按块随机划分是最常见的错误之一。同一幅点云相邻块之间有大量重叠点随机切会让验证集里的点出现在训练集中指标虚高得离谱上线就现原形。正确做法是按 tile 或按空间区域划分拿几幅独立的 LAS 做验证或者在同一幅图里按坐标范围切出验证区。类别平衡同样要处理。地面占 70% 时模型把所有点预测成地面就能拿到很高的 OA但 IoU 惨不忍睹。两个手段配合用一是在分块阶段对含稀有类的块做上采样建筑块重复采样两次二是损失函数加类别权重权重按1 / log(1.02 freq)这类公式算别用 1/freq容易把稀有类权重拉到几百倍导致训练不稳。4. 模型侧改造把 PointNet/PointNet 的类别数对齐到你的标注体系4.1 分割头与 num_classes、ignore_index 的对应关系PointNet 系列做分割时分类分支和分割分支共享特征提取分割头的输出通道数必须等于你的类别数。原始实现里num_classes通常是 13 或 16对应公开数据集换成自己的数据必须改改完还要检查两处最后一层 conv 的 out_channels以及损失函数里的ignore_index。import torch import torch.nn as nn NUM_CLASSES 4 # 对应 CLASS_MAP 里的 0..3 IGNORE 255 # 分割头输入是 (B, C, N)输出 (B, num_classes, N) seg_head nn.Sequential( nn.Conv1d(128, 128, 1), nn.BatchNorm1d(128), nn.ReLU(inplaceTrue), nn.Dropout(0.5), nn.Conv1d(128, NUM_CLASSES, 1), ) criterion nn.CrossEntropyLoss(ignore_indexIGNORE)NUM_CLASSES必须和CLASS_MAP的值域上限加一严格一致写成 4 就意味着映射后的标签只能是 0/1/2/3/255。ignore_index255让被标为 ignore 的点不参与梯度这是处理未分类点的标准做法。Dropout(0.5)放在分割头前是原始实现的配置小数据集上可以保留数据量大时降到 0.3 更稳。提示PyTorch 的CrossEntropyLoss默认期望输入是(B, C, d1, d2, ...)点云模型输出通常是(B, N, C)调用前记得logits.transpose(1, 2)这个顺序错了不会报错但训练结果就是玄学。4.2 Dataset 与 DataLoader把归一化锁在一个函数里Dataset 只做两件事取块、归一化。所有跟坐标变换相关的逻辑都收在这里训练和推理共用同一个类杜绝“训练用 A 方法、推理用 B 方法”的低级错误。import numpy as np import torch from torch.utils.data import Dataset def normalize(xyz): center xyz.mean(axis0, keepdimsTrue) p xyz - center scale np.max(np.linalg.norm(p, axis1)) return (p / max(scale, 1e-6)).astype(np.float32) class BlockDataset(Dataset): def __init__(self, blocks, augmentFalse): self.blocks blocks self.augment augment def __len__(self): return len(self.blocks) def __getitem__(self, i): xyz, label, idx self.blocks[i] if self.augment: # 绕 z 轴随机旋转 小幅抖动是最便宜的增广 theta np.random.uniform(0, 2 * np.pi) c, s np.cos(theta), np.sin(theta) rot np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]], dtypenp.float32) xyz xyz rot.T xyz normalize(xyz) return (torch.from_numpy(xyz), torch.from_numpy(label.astype(np.int64)), torch.from_numpy(idx.astype(np.int64)))坐标、标签、索引三件套一起返回idx在训练时用不上但推理阶段必须靠它把预测投回原始点序。增广放在归一化之前做旋转是因为旋转不改变尺度顺序反过来结果一样但读起来更顺。label转 int64 是CrossEntropyLoss的硬性要求忘了转会在第一个 batch 就报 dtype 错误。DataLoader 的num_workers设成 CPU 核数的 1/4 到 1/2设太大反而因为进程切换变慢。4.3 训练循环lr、batch、epoch 的实际取值超参没有万能值但有一组能跑通的起点Adam初始 lr 1e-3weight decay 1e-4batch 8 到 16epoch 80 到 120学习率用多项式衰减 power 0.9。import torch from torch.optim.lr_scheduler import PolynomialLR device torch.device(cuda if torch.cuda.is_available() else cpu) model PointNet2Seg(num_classes4).to(device) # 换成你实际的模型类 opt torch.optim.Adam(model.parameters(), lr1e-3, weight_decay1e-4) crit torch.nn.CrossEntropyLoss(ignore_index255) sched PolynomialLR(opt, total_iters100, power0.9) for epoch in range(100): model.train() for xyz, label, _ in train_loader: xyz, label xyz.to(device), label.to(device) logits model(xyz.transpose(1, 2)) # (B, N, C) - (B, C, N) loss crit(logits, label) opt.zero_grad() loss.backward() torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0) opt.step() sched.step()transpose(1, 2)那行是整段最容易翻车的地方前面提过。clip_grad_norm_设 1.0 是防 PointNet 在稀疏块上梯度爆炸的保险去掉也能跑但偶尔会碰到 loss 变 NaN。batch 大小主要受num_point约束4096 点时 batch 16 大概吃 8 到 12 GB 显存显存不够先降 batch 再降点数别反过来点数降太狠会破坏局部结构。4.4 评估混淆矩阵算 mIoU 与 OA语义分割不能只看 loss。每个 epoch 结束在验证集上累计混淆矩阵最后一次性算 IoU 和 OA比逐 batch 平均靠谱。import numpy as np def update_cm(cm, pred, gt, num_classes, ignore255): m gt ! ignore idx gt[m] * num_classes pred[m] cm np.bincount(idx, minlengthnum_classes ** 2).reshape(num_classes, num_classes) return cm def compute_metrics(cm): tp np.diag(cm).astype(np.float64) fp cm.sum(axis0) - tp fn cm.sum(axis1) - tp iou tp / np.maximum(tp fp fn, 1e-9) oa tp.sum() / np.maximum(cm.sum(), 1e-9) return iou, iou.mean(), oaupdate_cm把预测和真值压成一个整数索引再bincount比双重循环快几个数量级768 类以内都不会有性能问题。compute_metrics返回逐类 IoU、mIoU 和整体准确率。看结果时别只盯 mIoU逐类 IoU 里出现 0 的类要单独查——通常是这一类在验证集里样本太少或者映射表里根本没配。5. 避坑清单LAS 标签训练里最容易翻车的五件事5.1 现象loss 卡在某个值不降预测全是同一类原因通常不是模型不行而是标签被大面积 ignore 掉了。CLASS_MAP里只写了几个码剩下的全落进 255如果数据来自不同测绘软件分类码体系可能完全不一样实际有效标签只剩 5%。另一种情况是稀有类占比太低梯度被地面类淹没。解决办法转换前先跑 2.2 的直方图把类别分布打出来训练第一个 epoch 后统计 batch 里非 ignore 点的比例低于 60% 就回去改映射表。稀有类用分块上采样加类别权重双管齐下权重上限压到 20 倍以内。5.2 现象写回 LAS 后点云错位或者整片变成一个颜色原因是分块采样打乱了点序写回时按预测数组的顺序往原始文件里塞。分块时有放回补齐还会产生重复点一个原始点可能对应多条预测。解决办法分块时把idx一起带出来推理时按idx把概率累加回原始数组最后对每个原始点取 argmax。这一步没有捷径idx是唯一的后悔药转换阶段就必须存下来。5.3 现象验证集 mIoU 高得离谱换一块数据就崩原因是训练集和验证集存在空间泄漏。随机按块划分时相邻块共享大量边界点模型等于在背答案。解决办法按 tile 或按坐标范围划分验证区保证验证区的点从未在训练中出现。划分完抽查一下两边的 x/y 范围有没有交叠交叠了就是切错了。5.4 现象训练正常推理时同一条航带的数据预测全乱原因是归一化方式不一致。训练时在块内做去质心加球归一化推理时图省事用了全局 min-max同一栋建筑在两种变换下的尺度差了好几倍。解决办法把normalize抽成独立模块Dataset 和推理脚本 import 同一个函数禁止各写一份。推理时依然按块处理每块单独归一化不要试图对整幅图归一化后再切。5.5 现象预测结果里出现了根本没定义的类别号原因是映射表和反查表方向搞反了或者写回时直接把网络输出的 0..N-1 当成 LAS 分类码写进去了。LAS 里 0 是“未分类”、1 是“未指定”网络类别 1 写进去就变成了未指定。解决办法写回前强制走一遍inv_lut写完再读回来做一次np.unique检查确保所有值都在CLASS_MAP的键集合里。这一步我一般写成断言不在集合里直接抛异常别让它静默通过。6. 从预测回到 LAS重叠推理、投票与写回自查训练跑通只是半程真正交付的是带预测分类码的 LAS 文件。推理阶段和训练阶段有两个关键差异块要重叠、结果要投票。用 50% 重叠滑窗每个点会被预测 2 到 4 次把这些预测的概率累加起来再取 argmax比单次预测稳得多尤其在块的边界处。import numpy as np import torch def predict_tile(model, xyz, block2.0, stride1.0, num_point4096, num_classes4, devicecuda): model.eval() prob_sum np.zeros((len(xyz), num_classes), dtypenp.float32) hit np.zeros(len(xyz), dtypenp.int32) xs np.arange(xyz[:, 0].min(), xyz[:, 0].max(), stride) ys np.arange(xyz[:, 1].min(), xyz[:, 1].max(), stride) for bx in xs: for by in ys: m ((xyz[:, 0] bx) (xyz[:, 0] bx block) (xyz[:, 1] by) (xyz[:, 1] by block)) if m.sum() 64: continue sel np.nonzero(m)[0] if len(sel) num_point: sel np.random.choice(sel, num_point, replaceFalse) pts xyz[sel] center pts.mean(axis0, keepdimsTrue) p pts - center p p / max(np.max(np.linalg.norm(p, axis1)), 1e-6) t torch.from_numpy(p.astype(np.float32)).unsqueeze(0).to(device) with torch.no_grad(): logits model(t.transpose(1, 2)) prob torch.softmax(logits, dim1)[0].cpu().numpy().T # (n, C) # 只把参与预测的点累加回原始位置 prob_sum[np.nonzero(m)[0][:len(sel)] if len(sel) m.sum() else sel] prob hit[np.nonzero(m)[0][:len(sel)] if len(sel) m.sum() else sel] 1 hit[hit 0] 1 return (prob_sum / hit[:, None]).argmax(axis1).astype(np.uint8)上面这段里有一处需要按实际降采样逻辑对齐的写法降采样后sel和块内全部点不再一一对应累加时要严格按sel索引回原始数组不能按m的顺序。stride取block/2是重叠推理的关键重叠率越高边界越稳代价是推理时间线性增长。hit记录每个点被预测了几次没被任何块覆盖的点通常是研究区边缘会保持 0除之前要兜一下底。拿到逐点类别后写回这一步反而简单但两个检查不能省import numpy as np import laspy inv_lut np.zeros(256, dtypenp.uint8) for src, dst in {2: 0, 6: 1, 3: 2, 4: 2, 5: 2, 9: 3}.items(): inv_lut[dst] src pred predict_tile(model, xyz) # 网络类别 0..3 las_code inv_lut[pred] # 还原成 LAS 分类码 assert set(np.unique(las_code)).issubset({2, 6, 3, 4, 5, 9}), 出现未定义分类码 out laspy.read(tile_01.las) assert len(out.points) len(las_code), 点数必须和原文件一致 out.classification las_code.astype(np.uint8) out.write(tile_01_pred.las)断言是给自己看的。第一处检查反查后有没有冒出没定义过的码第二处检查点数有没有对上——predict_tile如果中途丢过点这里会立刻报出来而不是等到打开查看器才发现少了一片。写回时保持原 header 不动scale/offset 沿用原文件别新建 header否则坐标会漂。一个能省很多时间的自查习惯挑一块既有地面又有建筑的区域把真值和预测并排渲染成两幅图对比重点看建筑边缘和植被交界。mIoU 涨两个点可能只是地面拟合更好了边缘才是真正决定可用性的地方。从那以后我每次做 LAS 标签训练都先把类别直方图打出来把CLASS_MAP和inv_lut一起存进实验目录跑完第一个 epoch 就抽一块可视化核对宁可 early 阶段多花二十分钟也不要在写回阶段对着一整片单色点云找问题。希望帮到你。本文还有配套的精品资源点击获取