学习 《Local Parametric Surface Approximation With Automatic Order Selection From Position Data》

abstract:从位置数据获取解剖图对于导管与软组织相互作用的医学应用非常重要，为了在这些数据中改善自主导航，我们需要超出通常可用的非参数图的信息。我们提出了一种从位置数据中自动进行曲面顺序选择的局部曲面近似算法。传统的曲面拟合目标函数是从贝叶斯角度得出的。占用图中的后验概率作为权重合并到为曲面拟合选择的点上。我们新颖的交替最小化算法结合了使用贝叶斯信息准则的表面顺序选择。仿真证明了在存在噪声的情况下自动选择与潜在表面一致的表面顺序的能力。并且介绍了人工程序数据的结果。

I. INTRODUCTION

当与整合平滑位置数据的占用图进行对比时，人类导航员使用的解剖图显示出明显的插值（请参见图1）。在导管消融手术中，由于导管与组织之间的相互作用会影响导管的反应，因此，占用图的作用范围已超出传统的自主导航任务（例如，路径规划，避障）。在这里，我们提出了一种用于局部参数曲面逼近的鲁棒算法，可从嘈杂，不完整的数据中实时提供新信息

图1：左心房的解剖图数据，呈现给医师导航的表面标记为网格。占用图的级别集标记为V。此外，我们指出了边界体素X的局部选择。标记了左下肺静脉（LIPV），左上肺静脉（LSPV）和左心耳（LAA）。点云指示在肺静脉和LAA之间的脊上的一组已访问量的边界体素。作为参考，体素为1立方mm

本文主要有三点贡献。第一，我们描述了一种从非平面表面的占用图中选择点的新颖算法
    第二，我们扩大了用于曲面拟合的传统目标函数，以包括可从占用图中获得的后验概率
    第三，最重要的是，我们在交替最小化算法中纳入了自动曲面顺序选择功能，这对于区分潜在曲面的曲率与数据中的噪声（误差）至关重要。

在这里，我们仅使用治疗导管中的位置数据来近似解剖表面。处理流程如图2所示。位置数据通常具有亚毫米级的精度。这些数据受到强烈的心跳和呼吸运动的影响，我们使用未公开的算法对其进行抑制。由于可以使用本地化数据，因此我们的问题不包括同时进行本地化和映射（SLAM），并且构造占用图很简单。可以基于其他传感器数据（例如，尖端力或磁转矩）将后验概率分配给体素。从占用图来看，我们考虑一个密集的二元矩阵V，对已知内部空间进行阈值化。由此，我们确定了点的无组织集合X，我们将其解释为潜在表面非均匀采样的噪声，从这些点，我们确定一组控制点A，定义一个贝塞尔曲面。参数nu，nv设置表面顺序并自动确定。

关于Bezier曲线和点云的拟合的相关概括已发表了大量著作.从较高的层次上讲，我们将无梯度搜索方法与基于梯度的搜索方法（例如，高斯牛顿）（包括其他方法）区分开来，在大多数情况下，假定表面顺序是固定的。我们发现了Iglesias等人的工作。一个有趣的例外，因为它们包括用于表面顺序选择的贝叶斯信息准则（BIC）。我们的曲面拟合方法可以概括为使用距离最小化方法（PDM）对位置参数和控制点进行交替最小化。我们的方法是独特的，因为我们利用BIC在每次更新控制点时有条件地增加表面顺序。固定表面顺序后，其他方法已证明具有更快的收敛速度[2]，[16]，并且在必要时可以用于更有效地细化表面逼近。

本文的其余部分安排如下：在第二部分中，我们介绍了指导算法设计的分析和统计模型，在第三部分中，我们描述了用于点选择，模型拟合和模型选择的算法，在第四节中，我们使用分析模拟来量化性能，并论证关于人工程序数据的结果，最后的评论在第五章中给出。

第二章，models

这项工作的重点是将参数化曲面拟合到可用数据。为此，我们描述了两个模型：用于表面的分析模型和用于数据的随机模型。随后，我们将表面拟合描述为最大似然估计，并且将使用相同的随机模型自动确定表面顺序。    

A小节，参数面

我们为分析模型选择Bezier曲面。Bernstein多项式b是Bezier曲面定义的基础：

我们将所有基本函数的值收集为向量值函数b

参数n确定b的长度n + 1，当从上下文中清除时我们将忽略它

贝塞尔曲面只是伯恩斯坦多项式的加权组合.对于每个坐标，我们定义s:

在这里，我们使用分别表示方向u和v的多项式阶数:

该函数通过矩阵的一部分进行参数化，我们定义Bezier曲面，将（3）式的三个副本连接为映射s

在后续表达式中，将A重塑为二维矩阵会很方便,我们之所以做出这种区分，是因为从A可以清楚地看到，但是仅从那里是模棱两可的.
    重新声明式3和4:

式中⊗代表克罗内克积

B.随机数据：

令代表点的索引集.我们对数据建模。

此处表示随机误差项。此偏移量归因于许多因素，包括体积的量化，生物运动和模型误差。根据中心极限定理，我们认为误差是高斯分布的。为简单起见，我们将这些术语建模为独立的并按分布。权重代表编码我们对每个的置信度的占用图的后验概率。降低与相关的，增加了相应的模型方差。
    让统称模型参数。我们定义列向量，该向量包含与的索引元素关联的坐标。因此，数据的对数似然（由θ参数化）为

随后，我们寻求θ的最大似然估计。该表达式对于选择很有用，它们决定了A的大小，因此决定了θ中的参数总数。

第三章，算法

对于我们的问题，不能立即使用。必须根据已知体积的记录确定它们。在下文中，我们首先提出一种用于从二进制占用图中选择和的算法。然后，针对固定的拟合曲面参数。最后，我们考虑选择。

A. 选取点

令表示通过阈值占用图阈值获得的导管访问的体素。我们在本节中的目标是确定这些体素的子集，这些子集局部地近似单个表面。我们使用卷积和阈值识别内部边界体素。

选择内核矩阵有一定的灵活性。例如，要确定以至少包含个外部体素的立方邻域为中心的内部体素，我们选择

这里表示上限运算符，我们假设基于一个的索引.在应用中，我们发现的时候为我们的问题提供了敏感性和特异性之间的合理平衡。

表示的体素通常构成多个曲面。我们将此归因于矩阵V中贴错标签的体素。例如，导管可能无法访问某些内部体素。接下来，我们寻求识别与单个曲面相关联的局部体素子集（由表示）。

当两个点都由表示并且在同一邻域内时，我们认为它们与同一曲面关联。我们总结了下面的条件。

在最终等式中，我们使用表示指标矩阵，其中元素{i; j; k}是1.我们用表示所有元素都为1的矩阵。以的初始点（i; j; k）开始，我们扩展了与迭代卷积的点的集合。在每次迭代中，迭代次数近似于从新添加的点到初始点沿曲面的距离。算法1中总结了点选择过程

Algorithm 1 POINTSELECTION 。在本地确定与查询点有关的单个表面的近似点，输入二进制矩阵V表示内部体素，而δ表示查询点。我们假设内核的大小为W； J不超过V的大小。指示的点与R中的非零条目相关联。R的值近似等于表面到查询点的距离。我们用*表示3D卷积，用◦表示逐元素乘法。

    通常，我们假设算法1在本地应用，因此可以快速执行3D卷积。当V代表已知体积的截断集时，沿的周长的体素可能贴错标签。应该强制执行最大迭代次数，以使J * R永远不会指示周围体素。

在其余算法中，我们将不引用V或R。相反，我们使用R来确定点的索引集。每个元素与R的非零条目相关联。此外，我们考虑了一组权重索引。这些可用于强调更接近查询点的点，或表示占用地图中的随机先验。

 

B. Model 匹配
    仅（8）中的受参数A，u和v影响。发现这些参数的最大似然估计最小。

要使此表达式最小化具有挑战性，因为（3）涉及其参数的高阶乘积。但是，这些参数为块协调下降提供了自然分解

我们定义如下参数：

迭代更新参数导致交替最小化算法。

通过固定A，（12）可进行分解。这样，分解为个二维优化问题。对于每个，我们求解下面的公式：

相对于u和v的（16）的梯度和Hessian可通过解析获得（有关详细信息，请参阅附录A）。它们对于u和v而言不是恒定的，并且不能保证Hessian是正定的。但是，可以使用牛顿法和Armijo回溯法快速确定局部最小值。

相反，具有封闭形式的解决方案。为方便起见，我们将水平连接为矩阵，并定义。

使用组成对角矩阵Λ。然后我们将（12）重新声明为矩阵方程式:

用矩阵范数替换向量范数：通过这种形式，的解决方案显而易见。

通常的缩放比例很差，并且如果不进行正则化，则最佳A将包含远离元素的控制点。正则化先前已解决（例如Tikhonov [7]和整流罩[10]，[16]）。使用Tikhonov正则化，该解决方案仍然可以以封闭形式使用。不失一般性，我们假设以原点为中心。该假设可能需要在拟合之前转换初始点云，在拟合之后将A转换为:

   在这里，我们用表示所有值都为1的向量。我们定义正则化函数为

A最小化（19）被发现求解方程的线性系统

到目前为止，我们假设常数。它们确定A和B的大小，但对u或v的大小没有影响。这样，的每次迭代都提供了更改的机会。接下来，我们考虑如何选择这些参数

C. Model 选择:

我们可以通过选择大的来强制（8）任意小。但是，增加这些参数会导致在解剖学上不准确的高阶表面。我们提出把模型选择的问题转换为选择。为此，我们使用贝叶斯信息准则（BIC）[17]

让q索引候选模型。在我们的方案中，每个包括与唯一表面顺序相关的最佳参数。随着的增加，BIC渐近地逼近:

表示模型q的最大似然估计。标量表示与模型q相关的模型参数总数。

 

表达式（21）使用（8）求解，需要σ平方。给定最大似然估计A，u，v（如第III-B节所述），则σ2的最大似然估计为：

将（23）插入（8），我们发现最大化（21）的模型选择等效于选择q最大化统计：

q索引包含σq2的候选模型θq，并根据（22）确定。算法2描述了估算多对的A并选择产生最大t的结果的过程

Algorithm 2 MDLSELECT：
    在增加模型阶数的同时估计A和统计量t。
    返回具有最大统计量的模型。
    这里COMPUTEA指的是（20）的解决方案，而COMPUTESIGMA2指的是（23）。
    虽然参数u，v以θ返回，但它们没有改变。

我们承认，我们对BIC的应用并非完全基于统计理论。对于我们的问题大小，BIC的渐近逼近并不准确。此外，该标准要求对参数进行最大似然估计。这将需要将迭代算法收敛到全局最小值。为了减少计算量，我们随后在迭代改进θ的同时应用（24）。如果没有收敛，θq不会代表最大似然估计。尽管如此，在应用中我们发现（24）在数据拟合和模型顺序之间提供了合理的平衡。

算法（3）中给出了表面拟合过程。我们假设u，v已初始化，例如，使用SVD将X投影到2D子空间上。在MDLSELECT中，每次更新A时，我们都有可能根据BIC增加表面顺序。对u，v和σ2的更改对于停止条件都是有用的（此处未解决）
    
    Algorithm 3 FITSURFACE：
    使表面适合点的集合，以迭代方式更新u，v和A。UPDATEPOINT引用（15），它需要一个迭代求解器。

第四章，结果。

我们的方法的主要优点是快速，自动的表面顺序选择，对过度拟合和不足拟合均具有强大的效果。为了演示性能，我们考虑根据噪声数据进行插值。我们将名义参考面定义为点的无限集合。然后，我们定义两个唯一的子集，分别用于训练和测试目的的。为了方便起见，我们将使用引用这些集合的基数。水平连接的元素，我们定义矩阵。我们不假定直接可用。取而代之的是，我们假设观测数据受到与（7）一致的加性噪声Xtr的影响

从可用的XTR，我们拟合一个表面，估计，我们用它来确定。这些数量在图3中显示为示例表面。

图3说明：
    表面近似值，训练和测试数据的可视化。参考曲面是缩放和旋转的Rosenbrock函数。通过对参考表面进行采样并添加噪声来生成训练点.我们同时覆盖参考表面和近似表面.

参数，或更确切地说是A，也定义了一个表面。我们将此表面称为由坐标索引的点的集合。为了进行测试，我们对每个元素及其最近邻居之间的距离感兴趣。我们确定最近的邻居求解（15）。然后，我们定义水平的矩阵。

为了量化算法性能，我们考虑下面两个指标：

注意，当wi = 1时（26）等于（23）。

训练数据的两个属性对算法性能有重大影响。特别是和Y的方差。为了证明这一点，我们考虑Rosenbrock函数的缩放版本。我们对参考函数进行采样，并对所有采样点应用一致的随机生成的旋转。然后，我们将旋转的采样点划分为，随机生成元素，我们使用（25）生成。根据，我们使用算法3估算，该算法还确定。从通过对每个求解（15）来估计。使用，我们根据（27）确定。对多个对重复此过程100次。平均值如表1所示。

    在表I中，iter列表示AM迭代的平均次数（算法3中的t）。我们将最大迭代次数限制为10，这足以确定适当的曲面顺序。表I中的size列表示A中控制点的数量 。我们强调这不是为了避免依赖于。有两个趋势可以观察。首先，size随着nTR而增加。其次，尺寸随着σy^ 2的减小而增加。这两种趋势都是取决于BIC。类似地，我们发现较大的模型尺寸与较低的相关联，因为BIC可以避免过度拟合。最后一栏描述了曲面拟合的平均运行时间（不包括点选择）。

为了进一步展示自动顺序选择的优势，我们对比了固定顺序模型的拟合误差。除了Rosenbrock测试面之外，我们现在还将平面视为第二个潜在测试面。直观地，我们期望低阶模型在Rosenbrock测试表面上表现不佳（拟合不足），而高阶模型在平面上表现不佳（过度拟合）。如图4所示：此外，我们发现我们的方法在两种情况下均表现良好，可以从嘈杂的测量中自动为潜表面选择合理的顺序。

图4说明：对比测试误差与噪声的关系。对于图4a和图4b，潜表面分别表示一个平面和Rosenbrock函数。图例指示固定顺序模型中假定的控制点数。

为了演示完整的算法性能，我们再次回顾了图1中描述的人工程序数据。在应用未公开的算法[5]以减少生物运动[4]之后，根据位置数据构建占用图。从占用地图中，我们提取出一个密集的二元矩阵V，该矩阵表示内部体积为带有边长为1mm的立方体素。我们确定将算法1应用于边长为15mm的立方体积的。此集合在图1中表示为点云。使用算法3，我们将曲面拟合到。结果如图5所示。在这种情况下，，我们使用。对于拟合的曲面，。这种情况表示解剖结构的区域具有较高的曲率，但是产生的误差远低于数据的量化。为了将结果与呈现给医师导航器的表面进行对比，我们将显示限制为仅表示与元素之一最接近的点的网格顶点。到非参数显示表面（最近顶点）的中值距离为3.38毫米，而到我们的参数表面的中值距离为0.24毫米（正交距离）。

图5说明：从人类程序数据近似的局部表面。呈现给人类导航员的非参数表面被标记为网格。我们的参数曲面逼近和中间点云分别标记为。的标记颜色表示距近似表面的距离（法向），以mm为单位。

第五章. 推论

我们已经展示了一种用于表面逼近的迭代算法，该算法可以快速收敛并且对于过度拟合具有鲁棒性。我们认为这种方法非常适合与消融导管的自主导航相关的独特挑战。

我们的曲面拟合算法仅解决局部搜索问题。换句话说，我们假设初始参数选择接近全局最小值。我们通过将X投影到二维子空间（SVD）来选择。该策略的有效性将取决于潜在表面的曲率。例如，当点云近似于圆柱体时，这是不合理的。限制表面近似程度是一种缓解策略。较大的曲面可能需要全局搜索或将多个曲面近似拼接在一起，这为进一步研究提供了机会。

附录
可分离问题的梯度和海森矩阵。
（6）的雅可比式可以表示如下：

这里我们用表示（2）关于u的导数。随后，我们将使用表示二阶导数。他们的解析推导直接来自（1）。则（16）的梯度为:

为了方便起见，我们定义辅助标量值:

则（16）的Hessian为: