（三十六：2021.01.13）Pre-MICCAI 2019学习（一）《从植入粒子后的CT图像自动定位植入的粒子》

《Automatic localization of implanted seeds frompost-implant CT images》

• 讲在前面
• 摘要
• 论文内容
• • 1.介绍
  • 2. 方法和材料
  • • 2.1 算法描述
    • 2.2 基于强度的CT阈值分割
    • 2.3 几何区域分割
    • 2.4 几何区域连通
    • 2.5 确定种子位置和方向
  • 3.结果
  • • 3.1 最佳阈值确定
    • 3.2 使用体模进行算法测试
    • 3.3 使用实际的植入后CT图像进行算法测试
    • 3.4 算法可移植性测试
    • 3.5 植入后剂量
  • 4.讨论
  • • 4.1 算法的可移植性
    • 4.2 CT层厚上的选择建议
  • 总结

讲在前面

• 一.Pre-MICCAI是我为了进一步了解MICCAI某些相当有意义论文的发展历史而采取的一种追根溯源的方法，在这个系列中我会将这篇论文所提到的之前的论文进行一定的了解，便于更深层次的学习和了解整个算法要解决哪些之前的问题，同时借鉴前人的思想来丰富自己的知识库；
• 二.我设计了几种字体颜色用于更加醒目地表现关键的思想和主题：
  • 红色表示尚未理解透彻的一些概念
  • 蓝色表示对原来的理解做的一些修改或补充
  • 绿色表示此处需要参考的论文其他部分
  • 橙色表示本文的重要关键字
  • 我会用删除线将自己曾经不到位的理解进行删除 es》)
• 三、本文基于论文**《A Deep Regression Model for Seed Localizationin Prostate Brachytherapy》**

摘要

本文介绍了一种从一系列植入后CT图像中自动定位植入种子的方法。对接受前列腺近距离放射治疗的患者进行了植入后CT研究。 在每个CT切片中使用二进制阈值对明亮的区域进行分割，并收集这些区域的几何信息。 通过在每个切片中进行基于几何的过滤，将较大的区域（可能包含两个相连的种子）分成较小的区域。 使用基于几何学的连接搜索算法对每个区域切片进行纵向切片，分析沿纵向的区域连通性，从而将连通区域组合为一个对象。 将每个对象的加权质心作为这些位置。 该方法在含种子的前列腺幻象上进行了良好的测试，并通过患者的CT研究进行了统计分析。统计分析表明，该方法可以实现99％以上的检测率，并且定位准确，速度快。 对于在CT上定位植入的种子，它是可靠且方便的，可用于辅助前列腺近距离放射治疗的植入后剂量测定。

论文内容

1.介绍

经直肠超声（TRUS）引导的永久性近距离放射治疗植入物已被越来越多地用作前列腺癌的治疗选择，其中 I 125 I^{125} I125或者 P d 103 Pd^{103} Pd103放射性种子被永久植入前列腺中。 始终使用最佳的预计划来指导粒子植入过程，旨在将规定的剂量覆盖范围提供给前列腺和备用尿道和直肠。 但是，由于植入过程中的一些技术问题，例如前列腺运动和浮肿，针头放置和粒子下沉，实际的种子位置将与预先计划的位置不同。 因此，植入后剂量测定是验证实际剂量覆盖范围和评估植入物质量的重要步骤。

在最近的一份报告中，美国近距离放射治疗协会（ABS）建议对所有接受永久性前列腺近距离放射治疗的患者进行植入后剂量测定，原因如下：

• （a）记录前列腺和正常邻近组织在植入物的整个生命周期内的实际剂量;
• （b）通过使用补充的外部束照射或其他粒子植入来纠正由于前列腺剂量不足引起的可能的治疗失败;
• （c）帮助刚刚开始进行永久性粒子植入的医师评估和修改其植入技术，并帮助经验丰富的医师完善和完善该程序；
• （d）用于将来的结果分析，例如比较各个机构的治疗结果，并用作预期的多机构临床试验的质量保证工具（质控）。

通过定位粒子位置，已采用几种方法确定植入后剂量。 这些包括使用两到三张胶片，使用计算机断层摄影（CT）图像以及使用MRI图像。 目前，基于可用性，成本以及对前列腺以及粒子成像的能力，ABS推荐使用基于CT的剂量测定法。

大多数基于CT的粒子定位方法会首先确定连续2D切片上的种子位置。 由于种子通常会暴露在多个相邻的切片上（取决于成像厚度），因此此过程往往会比患者体内的种子数量识别出更多的种子位置。 因此，使用迭代过程将初始位置的数量减少到种子的确切数量。 Brinkmann和Kline（1998）和Liet等人（2001）分别报告了自动种子定位的方法。 他们的方法基于数​​据聚类技术，该技术通过对CT图像中分离的斑点进行分组来确定种子位置。 由于这些方法仅取决于亮点的大小，因此检测结果将极大地依赖于CT研究的设置参数和阈值。随着CT切片厚度的增加，每个种子的像素数变得更加取决于粒子轴向平面的相对方向；另一个重要因素是，由于粒子丢失（例如，通过外部迁移或排出种子），在成像时前列腺中的粒子数量可能不等于植入的数量。 这两个因素会使使用这些方法识别单个粒子变得困难，尤其是当轴向分辨率降低时。

在本文中，我们描述了一种基于几何的识别方法，该方法可自动确定粒子中心的三维（3D）坐标并将其位置记录在植入种子后前列腺的CT研究中的患者解剖结构中。 这种方法并不是强制性找到预期的植入种子数量，因此，可以解释粒子丢失的影响。 分割的粒子位置可以被转移到治疗计划中以计算剂量分布。 在足够高的切片分辨率下（例如1.5毫米或1毫米切片厚度），该算法可以提取粒子的3D方向进行各向异性计算，从而比点源方法更精确地进行剂量测定。 已对包含粒子的前列腺体模以及来自患者的CT进行了测试。

2. 方法和材料

2.1 算法描述

图一：

图1显示了该算法的流程图，该算法以C++程序实现。 CT图像以DICOM 3.0格式导入到程序中。 经过二进制阈值处理后，程序会在每个CT切片中的明亮区域上收集相关的几何信息（例如大小，位置）。 初始的基于几何的过滤算法将大区域（可能包含连接的粒子）分割为单个较小的区域。 每个区域的描述（例如位置，大小，质心，长轴和短轴）都输出到多维矩阵中。 下面的基于几何的连接搜索算法逐一搜索从顶部切片到底部切片的每个区域，并根据几个标准确定连续切片上的连接区域是否属于一个粒子体积，每个最终连接体积均被视为一个粒子且加权质心被计算为种子位置。 如果粒子由两个或多个区域组成，则可以确定粒子方向。根据CT图像的像素间距和得到的粒子位置，将坐标从像素单位转换为厘米单位，以进行随后的剂量计算。 该算法包括三个关键步骤：图像分割，基于几何的初始过滤和基于几何的连接搜索。 我们将在本节的其余部分讨论这些步骤。

2.2 基于强度的CT阈值分割

本研究中使用的CT图像具有512×512像素，像素值介于0和4095之间，0为黑色，而4095为白色。 尽管像素值不在以0表示水而-1000表示空气的Hounsfield单位（HU）中，但是可以将它们视为与Hounsfield单位的线性偏移。

在我们的研究中，为了选择用于分割每个CT图像中放射性粒子的最佳值，我们首先对38个CT图像集进行了统计分析。 对于每种设置，首先选择一个感兴趣的区域（ROI）以覆盖前列腺区域（在我们的程序中，可以通过将鼠标拖到CT图像上来快速定义矩形ROI），然后使用低阈值（1500）进行分割 从背景组织中剔除出较高价值的像素。 然后根据像素连接性将这些像素聚类到一些分离的对象。 对于每个对象，记录其最大值和大小。 放射性粒子具有较高的最大值和较小的尺寸，而骨骼和钙化物具有较低的最大值和较大的尺寸（如图2所示）。 在通过CT图像观察确认骨骼和钙化的检测之后，记录每种情况下的钙化物数目，钙化物的最大像素值和骨骼的最大像素值。 本文描述的自动算法中使用的阈值大于从该统计研究中获得的骨骼和钙化物的最大值。（也就是说，仅仅统计38个CT的骨骼和钙化物来确定手动分割的阈值，在这个阈值之上过滤掉钙化物，怎么感觉这么草率呢？）

图2：

2.3 几何区域分割

上一步在每个CT切片中创建了许多分段区域。 这一步将设置一个截断尺寸，以区分单个种子与CT平面内的重叠种子对。 截断尺寸通过以下方法自动确定：

在阈值化之后，每个分割的区域具有大小，即，大于阈值的像素数。 由于大多数区域仅包含一个种子，因此它们的大小相似。几个区域的大小可能几乎翻倍，因为它们包含两个相连的种子区域。 将整个CT图像集中所有分割区域的大小收集为数据数组，并进行统计分析。 将平均值加此数据数组标准差的三倍设置为截断大小。 大于此大小的区域被视为两个相连的粒子区域，因此从阵列中删除。 重新计算了新阵列的截断大小，并且再次消除了大于该大小的区域。 重复此过程，直到没有任何区域的大小大于当前的截断尺寸为止，该截断尺寸因此被设置为最终截断尺寸。

对于比截断尺寸大的区域，基于几何的分割算法（图3）被分成了两部分。 首先，通过在区域边界中搜索最大距离来找到区域的长轴。 然后计算该区域的质心，并将质心投影到长轴上，将长轴分为两部分。 计算这两个部分的各个质心，并将它们投影到长轴上作为起点和终点。 沿着长轴，从起点到终点搜索，找到了垂直距离最短的位置。 该位置用作分割点，从该分割点创建了两个子区域。

图3：

2.4 几何区域连通

在此步骤中，如果满足多个规则，则将可能属于一个粒子体积的多个连续切片中的粒子区域连接起来。 该算法从头到尾连续搜索每个区域，并将每个区域分配为

起始

粒子区域或

后续

粒子区域。 第一个CT切片中的所有区域都被视为

起始

粒子区域。 对于每个后续切片中的每个区域，该算法首先检查它是否连接到上一个切片中的任何粒子区域。 计算了从该粒子区域到先前切片中所有粒子区域的横向距离。 然后选择距离最小的种子区域。如果这两个粒子区域具有连接的体素，则将应用

最大粒子体积

标准。

最大粒子体积

标准意味着，如果添加此粒子区域使粒子量的总大小超过

最大粒子体积

，则该粒子区域将被分配为新的

起始

种子区域； 否则，它是

后续

粒子区域。 根据所植入粒子的数量 N N N 自适应地确定

最大粒子体积

，具体实现方法如下：

• 首先，确定

  平均粒子体积

  ：为所有种子区域的总大小除以 N N N ；
• 然后，两倍的

  平均粒子体积

  作为初始的

  最大粒子体积

  ,总共的待侦测粒子数目也会被确定。 因为

  平均粒子体积

  足够大，大到可以连接物理上分离的粒子，因此检测到的粒子总数 N 1 N1 N1小于预期的粒子数 N N N。这些检测到的粒子的体积按降序计算和排序。 设 M ＝ N − N 1 M ＝ N-N1 M＝N−N1。 选择第 M M M 个粒子体积作为最终的

  最大粒子体积

  ，然后将再次执行连接算法。

如果

后续

粒子区域位于

起始

粒子区域的第三（或更多）切片上，则

直线规则

也将适用。 这意味着，如果该区域的质心在粒子区域的直线上，它将连接到粒子上； 否则，它将被分配为新的

起始

种子区域。

2.5 确定种子位置和方向

最后一步是确定每个粒子的3D位置和方向（如果仅在一个切片上出现了，则无法确定方向）。 对于显示在两个或更多CT切片上的粒子，使用面积加权方法确定3D坐标，方法如下：

x = ∑ i = 1 N x i A i ∑ i = 1 N A i , y = ∑ i = 1 N y i A i ∑ i = 1 N A i , z = ∑ i = 1 N z i A i ∑ i = 1 N A i x=\frac{\sum^N_{i=1}x_iA_i}{\sum^N_{i=1}A_i}, y=\frac{\sum^N_{i=1}y_iA_i}{\sum^N_{i=1}A_i}, z=\frac{\sum^N_{i=1}z_iA_i}{\sum^N_{i=1}A_i} x=∑i=1N​Ai​∑i=1N​xi​Ai​​,y=∑i=1N​Ai​∑i=1N​yi​Ai​​,z=∑i=1N​Ai​∑i=1N​zi​Ai​​

• N N N是所有包含粒子的切片数目；
• ( x i , y i , z i ) (x_i,y_i,z_i) (xi​,yi​,zi​)是第 i i i个切片上粒子的坐标；
• A i A_i Ai​是第 i i i个切片上的粒子的面积；
• ( x , y , z ) (x,y,z) (x,y,z)是结果：面积加权坐标。

粒子的方向通过计算 ( x , y , z ) (x,y,z) (x,y,z)得到的 α \alpha α和 β \beta β来表示：

α = x N − x 1 ( N − 1 ) Δ , β = y N − y 1 ( N − 1 ) Δ \alpha=\frac{x_N-x_1}{(N-1)\Delta}, \beta=\frac{y_N-y_1}{(N-1)\Delta} α=(N−1)ΔxN​−x1​​,β=(N−1)ΔyN​−y1​​

• Δ \Delta Δ是CT的厚度。

3.结果

3.1 最佳阈值确定

在我们研究的38个前列腺患者CT中，在19个病例上发现了41处钙化物，这些钙化物的最大HU值在1232到2338之间，骨骼的最大HU值在1870到2322之间。这表明，2338可能是最好的过滤掉钙化物和骨骼的阈值。但是考虑到算法向其他数据的迁移，我们将阈值设置为2400。

3.2 使用体模进行算法测试

在植入了56个假粒子的前列腺体模上测试了这种自动粒子定位方法。 体模的尺寸为5.0×4.5×4.0 cm3（长×宽×高），直径为1cm的尿道嵌入模型中。 假粒子的直径为0.8毫米，长度为4.5毫米。模型的ACT切片如图4（a）所示。 这项研究评估了该算法在种子定位中的准确性及其分离位置紧密的种子的能力。更近似的模仿实际植入物，在不同角度引入某些针头，以验证算法能够以随机方向准确定位种子。 使一些粒子彼此接触以形成多个粒子的交叉。 在使用1.5毫米厚度的连续切片获取的CT数据集上测试了自动方法。 正交膜定位技术是由一名观察员执行的。 观察者获悉图像中包含的种子总数，将种子坐标与自动化方法的结果进行比较，最大，最小和平均差异分别为2.72 mm，0.17 mm和1.63 mm。 据报道，约11％的种子位置超过2 mm。 识别出的种子的3D坐标被投影到CT扫描之前获取的AP和侧面侦察视图上。 图4（b）和（c）显示了CT检测结果与实际X射线侦察图像之间的匹配：

图四：

3.3 使用实际的植入后CT图像进行算法测试

使用来自两名患者的影像学研究，使用CT和正交膜技术通过自动方法确定了粒子坐标。 对植入了86125粒子的患者1进行了1毫米和3毫米厚度的连续CT图像扫描。 图5显示了该患者CT检测到的粒子位置与X射线胶片之间的比较。 植入了69103 P d Pd Pd粒子的患者2仅在3毫米厚的CT图像中进行了扫描。 在CT扫描的同一天为每位患者获取AP和横向正交胶片。 将所得的粒子坐标与相应3D位置之间的最大，最小和平均差异的自动方法的结果进行比较（表1）。 对于CT厚度为1 mm的患者1，所有种子均被识别，并且85个种子具有3D方向（其中78％记录在三个相邻的切片上）。剩余的一颗种子仅记录在一个切片上，因此无法确定其方向。

图5：(a)是将CT检测的结果和X光进行覆盖的AP显示；(b)是侧面显示

表一：

3.4 算法可移植性测试

为了测试该算法在不同临床情况下的可移植性，使用四种不同类型的种子（ 103 P d ^{103}Pd 103Pd TheraSeed， 103 P d ^{103}Pd 103Pd PdGold， 125 I ^{125}I 125I 6711和 125 I ^{125} I 125I IoGold），三种不同的CT切片厚度（1 mm，3 mm和5 mm）和两种不同的FOV（ 16厘米和36厘米）。 所有测试用例均获得满意的结果。 该算法针对不同情况自动确定不同的参数。 例如，TheraSeed，PdGold和IoGold种子的截断尺寸通常在16 cm FOV下为27–32像素，对于3 mm切片厚度在36 cm FOV下为7-11像素，对于16 cm FOV为17–18像素。 切片厚度为5 mm，而对于3 mm切片厚度，碘6711种子的截断尺寸在16 cm FOV下为45–58像素，在36 cm FOV下为9–13像素。 表2列出了20个测试案例的种子类型，CT切片厚度，FOV，算法确定的参数以及植入和检测到的粒子数量。（由此也可以看出一个问题，对于不同性质的数据，算法内部的参数都需要人工调整，这里的问题在于，论文没有进行横向的对比，没有在控制住其他变量只改变病例的情况下做测试。）

表二：

使用Pentium 4 2.4 GHz计算机的平均处理速度为1秒。 与手动操作相比，这节省了大量时间。

3.5 植入后剂量

我们将这种自动粒子定位算法应用于近距离放疗后的28个病例。 这些包括18个 103 P d ^{103}Pd 103Pd 植入物和101个 25 I ^{25}I 25I植入物。 所有病例都以3毫米的切片厚度扫描。 根据每个CT图像序列的前列腺边界的手动轮廓线，计算每个患者的植入后体积和剂量体积直方图（DVH）。表3列出了研究结果。 在该表中，植入前的体积是从超声研究中获得的。 经验丰富的观察员对每个病例​​进行目视检查，以确认粒子的正确性。 最后一行给出统计总结，包括植入粒子总数（2634）和已鉴定种子（2610），平均鉴定率99.1±1.0％（平均值±标准差），平均体积增加1.25±0.26（平均值±标准差） ）和平均D90 / Rx剂量1.08±0.18（平均值±标准偏差）。 在此，D90是前列腺体积的90％所接受的剂量，Rx剂量是前列腺的规定剂量。 D90 / Rx剂量之比是反映植入后剂量测定中植入质量的重要因素，这是美国医学物理学会（AAPM）第64号任务组的建议。 125 I ^{125 }I 125I的D90 / Rx剂量（1.23）高于 103 P d ^{103}Pd 103Pd（1.00）。 这可能是由于 103 P d ^{103}Pd 103Pd病例的体积增加较高（1.33比1.10）和 103 P d ^{103}Pd 103Pd来源的剂量衰减更快。

一名经验丰富的观察员根据人工粒子定位对2例（ I − 125 I-125 I−125病例编号为10， P d − 103 Pd-103 Pd−103病例编号为18）进行了植入后剂量测定.D90的值分别为201和94。 这两种情况下的粒子自动定位结果均在1％的范围内。这表明，自动粒子定位方法与植入后剂量计算的专家手动定位相当。 粒子定位差异（与正交X射线胶片定位技术相比，平均差异为3 mm）没有转化为明显的剂量学差异。

表三：

4.讨论

4.1 算法的可移植性

该算法应该相对容易在其他诊所中实施，因为它与粒子的类型和CT扫描参数无关。

在本研究中，通过实验确定了用于将粒子区域分割出前列腺组织，钙化和耻骨弓骨的阈值。 在本研究中，发现2400的像素值足够高，可以消除那些属于所有测试用例的钙化和骨骼像素。 搜索属于非粒子区域的最大像素值的方法是该程序的一项独立功能，因此可以很容易地在其他诊所中实施以找到其最佳阈值。

该算法中使用的两个关键参数是来自CT，而不是来自外部输入。 分割双粒子区域的截断面积被计算为平均值加上单粒子区域大小分布的标准偏差的三倍。 测试结果表明，这是一个合理的选择，并使算法与粒子类型和FOV无关。 确定体素数量的

最大粒子体积

是基于知识的自适应过程，因为最初植入的粒子数是前提条件。 由于此输入可能不是CT上存在的确切种子数，因此不会强制算法匹配此输入数，因此本地化粒子的最终数可能会有小的差异。 该参数的自动确定使该算法与粒子类型，FOV和CT切片厚度无关。

4.2 CT层厚上的选择建议

在该算法中，较小的CT切片厚度（1毫米）将提高粒子定位的准确性，因为将加权的质心位置计算为这些质心。 表1体现了使用1mm CT切片厚度比使用3mmCT导致更高的精度。此外，在较小的CT切片厚度中，更多的种子将具有取向信息，因为粒子将出现在两个或更多个切片上。 研究发现，使用1毫米CT切片时，可以3D方向计算98％以上的粒子，而使用3毫米CT切片时可以减少到67％的种子，但是在常规情况下使用1毫米CT可能会被禁止（由于成本和时间的临床实践）。

另一方面，使用1 mm和3 mm CT（平均差0.2 mm）的结果之间的一致性表明，当不需要粒子方向信息时，使用3 mm CT是一个很好的折衷方案。 布林克曼还研究了5毫米CT切片，发现在1毫米和5毫米CT中匹配坐标之间的距离大于2.0毫米。 我们的研究表明3 mm CT是常规临床实践中的合理选择。

总结

本文从植入后的CT图像中描述了一种自动粒子定位算法，该算法与先前发布的方法之间的主要区别在于该算法依赖于3D粒子分布上的几何和强度信息，因此在处理粒子重叠和计算方向方面可能更准确。 对于那些仍处于选定阈值连接状态的粒子区域，该算法使用粒子大小标准来拆分连接的粒子区域，而不是像以前的方法那样将组合对象的质心坐标分配给每个连接的粒子区域。

该算法已使用20个病例进行了测试，其中有四种不同种子类型、三种不同CT切片厚度和两种不同FOV的病例。结果对于所有情况都是令人满意的，因此表明该算法与这些变量无关。 在Pentium III 600 MHz处理器计算机上，典型CT序列的处理速度为4 s，在Pentium 4 2.4GHz处理器计算机上为1s。