ITK 学到 Data Presentation 这一节终于遇到了图像之外的数据类型。前面一直跟 itk::Image 打交道像素规规矩矩排成网格处理方式非常统一。可真实项目里不是所有信息都在网格上配准时要标注的解剖点、激光扫描得到的散点、形状分析用的轮廓顶点都是离散坐标。它们需要一个自己的容器就是这篇笔记的主角 itk::PointSet。这篇内容把 PointSet 从创建、获取到访问的完整过程过一遍适合刚开始啃官方文档又不想只看概念的人代码都来自我本地实际跑通过的最小示例踩过的坑也一并写出来。1. 为什么稀疏离散点需要单独的结构1.1 图像用位置说话点集用坐标说话图像里每个像素的位置是由网格索引隐式确定的知道 (i,j,k) 就能定位数据连续紧凑处理起来很舒服。但现实中的点大多数时候没有这种隐形坐标格。比如配准任务里两幅影像上各标了几十个对应标记点这些标记点位置由物理坐标直接给出彼此之间没有固定的拓扑关系又比如统计分析里的关键点集可能每隔几个毫米采一个特征点数量和顺序都不固定。这类数据如果硬塞进图像里要么浪费大量空白像素要么丢失维度信息所以需要另一种更离散友好的结构。1.2 PointSet 的底层并不只是一个点的数组一开始我以为 PointSet 就是个存放 itk::Point 的 vector改改接口而已。真正去看实现才发现它把数据拆成两层一层是点坐标本身另一层是与点关联的属性数据PointData同时每一个点都有独立的非负整数标识符。如果把一堆点想象成通讯录标识符是联系人编号坐标是姓名PointData 是备注编号和内容分开存备注也可以为空。这个设计的意义在于很多算法只需要坐标、不需要属性另一些算法比如表面标签传播恰恰只关心属性而不关心坐标拆开之后两边都不受拖累。能跟普通 std::vector 拉开差距的原因其实也在这里把几个特征放一起对比会更明显。维度std::vectoritk::PointSet编号规则隐含在数组下标中必须连续每个点显式带非负整数标识符允许稀疏附加属性需要自己另建一个平行数组内置 PointData 容器与坐标分离存储访问模式按下标 O(1)默认 Map 容器按标识符 O(log n)稀疏支持差中间空位也得占用好只有存在的点才占条目算法协作需要自己封装与迭代器、变换、网格等 ITK 组件配套这张表给我最大的启发是PointSet 在 ITK 里不是点的数组的别名而是标识符到坐标、标识符到属性两组映射关系。理解这一点后面创建和访问时很多行为就说得通了。2. 模板参数与容器类型先把骨架搭对2.1 一行声明里的三个参数PointSet 最常见的声明是itk::PointSetdouble, 3但这两个模板参数的含义值得拆开看。第一个是 TPixelType指每个点关联的属性数据类型不是坐标类型第二个是 VDimension指空间维度。完整模板其实有三个参数templateclass TPixelType, unsigned int VDimension, class TMeshTraits class PointSet;第三个参数 TMeshTraits 是网格特征类一般不用手写默认值会帮你把内部容器、坐标表示类型都定好。坐标本身的类型是itk::Pointdouble, VDimension它由 MeshTraits 里的 CoordRepType 决定默认是 double。如果你想给每个点挂一个三维向量作为属性那第一个模板参数就写成itk::Vectordouble, 3比如itk::PointSetitk::Vectordouble, 3, 3。这个参数顺序很容易记混我自己的记忆方法是第一个参数描述点上挂了什么第二个参数描述空间是几维的。点上挂的东西和空间维度是两回事官方文档里经常看到PointSetdouble, 3就是因为最常见的需求是给每个点挂一个标量值。2.2 默认容器为什么是 Map 而不是 VectorPointSet 内部有两个容器PointsContainer 存坐标PointDataContainer 存属性。默认情况下这两个容器都是 MapContainer也就是像 std::map 一样按标识符排序存放。为什么不用连续数组因为 PointSet 要支持任意非负整数编号。比如你先后插入编号 0 和 100 两个点用 VectorContainer 就得给中间 99 个空位留内存用 MapContainer 则只需要两个条目编号 100 的点不会让内存膨胀到 101 个元素。换容器的代价主要落在访问模式上。Map 的随机访问、插入都是 O(log n)遍历时缓存不友好Vector 是连续内存顺序遍历快得多按编号随机访问 O(1)STL 互操作也方便。如果你的点编号从 0 开始严格连续递增而且主要做遍历和批量读取换成 VectorContainer 会有肉眼可见的性能提升如果编号稀疏、经常随机插入Map 更合适。这里给一句实操建议学习和原型阶段用默认 Map 足够不用折腾 MeshTraits。等真正做大规模点云处理时再根据访问模式考虑自定义容器类型到时候查 MeshTraits 头文件配置也不迟。3. 创建点集一个最小可运行示例讲清楚3.1 三步创建声明、放坐标、挂属性下面这个例子是我本地跑过的最小完整代码逻辑只有三步实例化 PointSet、设置点坐标、设置点数据。#include itkPointSet.h #include iostream int main() { using PointSetType itk::PointSetdouble, 3; auto pointSet PointSetType::New(); PointSetType::PointType p; p[0] 1.0; p[1] 2.0; p[2] 3.0; pointSet-SetPoint(0, p); PointSetType::PointType q; q[0] 4.0; q[1] 5.0; q[2] 6.0; pointSet-SetPoint(7, q); PointSetType::PixelType data 42.0; pointSet-SetPointData(0, data); std::cout NumberOfPoints pointSet-GetNumberOfPoints() std::endl; return 0; }跑完输出是NumberOfPoints 2而不是8。这一点非常关键尽管最大标识符是 7但 GetNumberOfPoints 返回的是内部容器里的元素个数只有 0 和 7 两个点存在所以是 2。Map 带来的稀疏感在这里体现得淋漓尽致。3.2 SetPoint 到底做了什么SetPoint(0, p) 并没有做任何空间预分配它只是把点作为 value、标识符作为 key 插到了 PointsContainer 里。如果再对同一个标识符调用 SetPoint旧坐标会被新坐标覆盖这是 Map 语义不是追加语义。SetPointData 同理往 PointDataContainer 里插一条记录。有一点要提醒点坐标和点数据的插入顺序没有强制要求。可以先 SetPointData 再 SetPoint也可以只挂坐标不挂属性。ITK 允许只有坐标没有属性或只有属性没有坐标的半空状态。这个设计初衷是灵活的但如果你不检查返回值就假设两边一定同步后面很容易踩坑第 5 章会专门展开。3.3 点数据和点坐标是两本账把 PointData 单独拿出来说是因为它和坐标的独立性是最容易被忽视的。PointDataContainer 有自己独立的 Size不跟随 PointsContainer 变化。PointSetType::PixelType value 99.0; pointSet-SetPointData(7, value); PointSetType::PixelType outValue 0.0; bool found pointSet-GetPointData(7, outValue); if (found) { std::cout PointData outValue std::endl; }这里 GetPointData 也返回 bool表示标识符 7 是否有对应的属性记录。如果你用一个还没有 SetPointData 的标识符去查拿到的是 falseoutValue 保持原值。记住程序里不要假设有坐标就有属性查询时一定要检查返回值。4. 获取与访问三种姿势和它们的边界4.1 GetPoint 返回 true/false而不是抛异常获取点是 PointSet 最常用的操作之一。ITK 的接口设计得很克制它不像 std::vector 的 at 那样抛异常也不会像 subscript 那样悄悄越界而是通过返回值告诉你这个编号到底存不存在。PointSetType::PointType outputPoint; bool ok pointSet-GetPoint(7, outputPoint); if (ok) { std::cout Point 7 outputPoint[0] , outputPoint[1] , outputPoint[2] std::endl; } else { std::cout Point 7 is not stored. std::endl; }这个接口对新手很友好因为你不用写 try-catch也不用担心程序崩溃。只要检查返回值就能安全处理不存在的编号。如果你用 GetPoint(100, p)它不会爆炸只会返回 false。不过这里藏着另一个隐患如果你绕过 GetPoint、直接对内部的容器调用 operator[]行为就完全不一样了我在第 5 章第一个坑里会详细说。4.2 遍历必须用迭代器不能默认按编号循环刚学 PointSet 时最容易犯的一个错就是把 GetNumberOfPoints 当成最大编号加一然后写一个标准 for 循环去访问// 错误示例编号稀疏时会漏点 for (unsigned int i 0; i pointSet-GetNumberOfPoints(); i) { PointSetType::PointType pt; if (pointSet-GetPoint(i, pt)) { /* ... */ } }前面已经说过编号 0 和 7 都存在时 GetNumberOfPoints 是 2这个循环只会访问 0 和 1编号 7 的点永远处理不到。正确做法是走迭代器让容器自己告诉你到底有哪些编号auto pointsContainer pointSet-GetPoints(); auto it pointsContainer-Begin(); auto end pointsContainer-End(); for (; it ! end; it) { auto id it-Index(); auto pt it-Value(); std::cout Point id pt[0] , pt[1] , pt[2] std::endl; }迭代器接口里Index()返回当前条目的标识符Value()返回点坐标。Map 容器遍历时按标识符升序Vector 容器按存储顺序。如果你在循环里还需要同时拿到属性可以按 id 再调一次 GetPointData。不过要注意坐标容器和属性容器可能大小不一致所以更稳妥的做法是先判断返回值。4.3 按访问模式选容器而不是看心情结合前两节访问模式决定了你该不该换容器。如果只是按标志符随机取几个点Map 完全够用如果是连续编号且要频繁按顺序遍历Vector 会快不少如果数据要交给外部库处理VectorContainer 还能通过 CastToSTLContainer 拿到内部 std::vector 的引用互操作非常方便。我自己在实际项目里的一般原则是点集规模在几千量级时容器选择对性能影响不大别过度设计规模到几十万、上百万时连续编号就果断上 Vector编号不规则但需要频繁随机访问时Map 的灵活性更值钱。没必要一开始就把代码写成最复杂的样子。5. 学习过程中实测踩过的几个坑5.1 坑一不要直接对容器做 operator[] 查点这是我最想写出来的一个坑。某次调试点云显示时我发现点数量莫名增加了很多而且新增的点全是 (0, 0, 0)。追了半天问题出在我查点时的写法// 危险写法 auto container pointSet-GetPoints(); PointSetType::PointType p container-operator[](100); // 100 不存在MapContainer 的 operator[] 和 std::map 行为一致如果 key 不存在它会先默认构造一个 Point也就是所有坐标都是 0然后插入容器最后返回这个新元素。你只想查一下编号 100 的点结果容器里真的多了一个编号 100 的点点数从 2 变成 3。后续遍历、保存、统计都会受影响。正确做法是永远用 PointSet 提供的 GetPoint或者在容器上用 GetElementIfIndexExists。GetPoint 内部走的就是安全查询不会修改容器。这个坑的教训是ITK 容器接口虽然暴露了 operator[]但那是给确定存在且想写入的场景准备的不是给查询准备的。5.2 坑二Point 是物理坐标不是像素索引PointSet 里存的是itk::Point它表示物理坐标系下的位置通常带单位比如毫米。而图像像素位置用的是 Index也就是第几行第几列第几层。两者是完全不同的坐标系。我一开始做配准时标注点明明在图像上却想直接把 Point 的值当成图像下标去取值结果自然全错。正确做法是先建立图像与物理坐标的关系。ITK 里图像本身知道自己的原点、间距和方向可以通过 TransformPhysicalPointToIndex 把物理坐标变成索引或者反过来用 TransformIndexToPhysicalPoint 把索引变成物理坐标。只要记住一句话就行Textures 里这叫 UV 坐标和世界坐标ITK 里这就叫 Index 和 Point别混用。itk::Index3 index; bool inside image-TransformPhysicalPointToIndex(physicalPoint, index); if (inside) { auto pixelValue image-GetPixel(index); }5.3 坑三修改点后忘记通知模板参数还容易看错还有一个隐蔽问题当你对已经创建好的 PointSet 做坐标修改时最好调用一次pointSet-Modified()。ITK 很多组件内部有管线缓存机制你把坐标改了但某些依赖这个 PointSet 的对象可能还在用旧缓存渲染出来要么没变化、要么数据错乱。我遇到的情况是可视化窗口里点完全不刷新查了代码才发现少了 Modified。另一个经典问题在模板参数上。itk::PointSetdouble, 3的第一个参数是 PointData 的类型不是坐标类型坐标类型藏在 PointType 里。如果你嫌 double 坐标占内存、想用 float不是把模板参数改成itk::PointSetfloat, 3这么简单——那样会让每个点的属性变成 float而坐标依然是 double。正确做法是通过 MeshTraits 里的 CoordRepType 去调整。这种模板参数表面上看对实际含义不同的坑只能靠多看源码和多写小例子来避免。写了这些其实最想说的还是那句PointSet 不是 vector 它是两本账加一堆标识符的结构。我的学习顺序是先跑通上面的 demo再故意制造几个坑比如稀疏编号、越界访问、物理坐标当索引亲眼看到输出不对之后再回去查文档印象非常深。下一步我打算把 PointSet 和 Mesh 的关系整理出来以及点和图像之间怎么互相转换那块在配准和形变场里太常用了。