TPAMI 2024 | 多维敏感性感知密度估计

文摘   2024-11-13 19:00   辽宁  

点击下方PaperEveryday”,每天获得顶刊论文解读

点击加入论文投稿、写作、阅读分享交流群

论文信息

题目:Sensitivity-Aware Density Estimation in Multiple Dimensions

多维敏感性感知密度估计

作者:Aleix Boquet-Pujadas; Pol del Aguila Pla; Michael Unser

论文创新点

  1. 优化问题的新公式:我们构建了一个优化问题来估计概率密度,这在处理多维数据时尤其重要,尤其是在样本分布不均匀时。

  2. 探测器灵敏度的整合:将探测器灵敏度作为概率密度函数(pdf)整合到我们的模型中,这允许我们更准确地模拟和处理来自物理传感器的数据。

  3. 正则化Hessian的核范数:我们通过正则化样条的Hessian的核范数来促进稀疏性,这使得方法在空间上自适应,并且对正则化参数的选择更为稳定。

  4. 边界条件的显式考虑:我们的样条展开允许我们为多维域定制边界条件,这对于处理周期性数据和其他特定类型的边界条件至关重要。

摘要

我们构建了一个优化问题,用于估计在不均匀概率采样的多维问题中的的概率密度。该方法考虑了探测器灵敏度作为非均匀密度,并利用了网格上样条的计算速度和灵活的边界条件。我们选择通过对样条的Hessian的核范数进行正则化来促进稀疏性。结果表明,该方法在空间上是自适应的,并且对正则化参数的选择是稳定的,该参数扮演了带宽的角色。我们在标准密度上测试了我们的计算流程,并提供了软件。我们还提出了一种新的PET(正电子发射断层扫描)重采样方法作为我们框架的应用。

关键字

加权密度估计,Hessian-Schatten范数,重采样,成像,重分组,PET。

Ⅰ 引言

越来越多的成像方式进入了低光子计数的领域,因此需要进行统计考虑。超分辨率显微镜和最近发射断层扫描[1]、[2]、[3]、[4]、[5]就是这种情况。然而,在大多数密度估计(DE)方法中,维度问题早在二维或三维时就会出现。随着域的维度的增加,均匀采样空间变得更加繁琐。事件不仅变得更加稀少,而且检测概率的一致性也在传感器元件之间变差。例如,在某些断层成像方式中,每个传感器都有其自己检测光子的独立概率,这在扫描仪中可能会有相当大的变化[6]。这类考虑通常通过测量的逐点加权纳入DE程序[7]。然而,低灵敏度区域仍然存在问题,因为权重的使用在检测概率的反演时引入了稳定性问题[7]、[8]。边界条件(BCs)是DE方法经常忽视的另一个考虑因素,但它们对于描述像周期性正弦图这样的域很重要。
这些考虑因素对DE的应用很重要的其他领域包括流行病学,手机查询和自然现象的研究[9],[10],[11],[12]。它们也受到敏感性问题的困扰,即:诊断能力,手机接收和气象站的范围。它们也涉及BCS,例如,参数化地球上的坐标。
本文中,我们开发了一个DE方法,将这些考虑纳入一个统一的计算框架。
  • 我们重新构建了加权DE以避免不稳定的反演。我们通过将灵敏度作为一个概率密度函数(pdf)纳入进来,该函数表征了探测器系统。这对应于度量的选取。
  • 我们用多个维度的基本样条来表达密度。由于存在一个底层网格,许多操作可以表示为可分离的卷积。
  • 我们明确关注BCs。样条展开使我们能够根据所考虑的多维域定制BCs。我们通过将似然嵌入到一个近端优化框架中来处理带宽选择。
  • 我们的正则化基于底层样条的Hessian的核范数。这使我们能够控制分段线性样条的节点的稀疏性,沿着单一轴。在更高维度中,它倾向于局部仿射的样条。效果是不变的,可以旋转、平移和缩放。
与标准DE方法相比,我们的实验表明,我们的框架更好地处理了不均匀灵敏度,计算上独立于高维度中的样本数量,并且能够无缝适应BCs。它们还证明了我们方法对带宽参数的鲁棒性。
我们专注于正电子发射断层扫描(PET)作为说明性应用。我们的研究动机源于最初的观察,即新扫描仪将正弦图数据转换为来自潜在泊松过程的点云(图1)。我们展示了我们的DE方法是适合PET重分组和在最先进的扫描仪中重建的。

文章的其余部分如下:在第二节中回顾了DE的最新进展,我们在第三节中提出了我们的框架的公式,并在第四节中进行了优化。我们在第五节中对标准密度进行了框架测试。第六节包含了应用动机,特别是成像应用,以及对PET重分组框架的测试。

II. 密度估计的最新进展

除了基于直方图的估计(HE),核密度估计(KDE)是非参数方法中最为突出的[13]。KDE和HE因其简单性和强有力的理论保证而受欢迎。大多数编程库没有实现其他方法,这进一步证明了它们的普遍性[14]。我们认为DE面临的四个主要挑战是:维度、权重、边界和带宽。

A. 核密度估计(KDE)

  1. 维度:KDE在统计上是强一致的,并且在相对温和的假设下是渐近正态的。(随着d的增加,收敛速度会下降[13]。)KDE的计算成本很高:评估单个点的密度需要与样本数量一样多的操作,这很难并行化。目前仍在进行工作以加速KDE[15]、[16]、[17]。一种方法是截断核的求和[18]。(选择核并不随着维度的增加而扩展,然而。)另一种方法是将样本投影到网格上,这使我们可以使用基于卷积的方法[19]。许多加速方法的权衡通常是以牺牲准确性为代价的。
  2. 权重:对于有限集合的样本,加权KDE估计器可以写成
    其中是样本的权重,是带宽为的核函数。权重的目的是补偿灵敏度。它们通常被设置为,其中处的检测概率。这是不令人满意的,因为灵敏度只是被反转,并且只在整个核的支持范围内而不是在样本点上被考虑。
  3. 边界:边界条件(BCs)和紧凑域的存在引入了KDE的偏差。已经提出了缓解措施,将紧凑区间转换为实数线,将数据扩展到边界之外,或者使用定制的核[20]。
  4. 带宽:带宽选择在KDE中扮演重要角色,并且仍然是研究的主题[13]、[21]、[22]、[23]、[24]、[25]、[26]。一些(“经验规则”)方法给出了给定样本集的最优带宽,例如通过最小化均方误差在假设正态性等假设下[21]、[25];它们也可以考虑到权重。在一些(“插件”)方法中,核的宽度根据样本局部调整,以增加计算复杂性为代价[22]、[26]。

B. 直方图估计(HE)

尽管HE的速度随着维度的增加而扩展得很好,但由于维度的诅咒,超立方体的箱子很快就变空了。准确地说,带宽选择的挑战在HE中以箱大小的选择的形式出现。直方图对带宽特别敏感,因为箱子不仅是离散的,而且还是量化的。结果,HE方法的收敛速度很慢[13]。与KDE类似,权重通过反转检测概率并入HE。

C. 其他方法

另一种DE方法是将估计器表示为正交级数[27]。在那里,通过截断级数展开来调整带宽,从而得到更平滑的估计。在logspline方法中,将密度的对数展开,并将其表示为定义在实数线上的一系列自由结点的样条[28]、[29]。然后优化样条的系数以最大化似然[30]。通过选择结点的数量和位置来调整带宽,这是根据Akaike信息准则等标准进行优化的[31]、[32]。一种方法是在系数固定后启发式地删除结点。另一种方法是同时优化结点和参数,但得到的最小化问题严重非凸。一种较少采用的方法是采用平滑样条来规范问题[33]。
虽然logspline的结果在1D中是有希望的,但只有少数工作扩展到2D,在那里自由结点的启发式变得复杂,实现比KDE更加涉及[34]。此外,它们不考虑权重或BCs。

III. 提出的框架

我们现在制定了我们的敏感性感知DE问题,并在基于网格的B样条的基础上构建方法。我们通过近端可解的正则化来处理带宽选择。

A. 问题表述

,是对应于感兴趣现象的pdf。我们的目标是基于由物理传感器提供的观测来估计,这些传感器的灵敏度取决于它们在紧凑域内的位置。
灵敏度图模拟了检测系统注册现象的实例,该现象由描述的概率(图1(b))。我们假设这种现象和检测是独立的过程。灵敏度图可以实验测量或从某些底层检测模型中理论推导。我们的要求是中几乎处处可测量且为正。因此,我们将检测到的事件的密度写成
这个表达式等同于
对于任何。为了方便,这里我们使用作为目标密度的未归代理。
注意,选择等同于选择一个度量;例如,为我们在(5)中制定的向量指数族选择一个参考度量。它在归一化中的存在避免了大多数加权密度估计方法中的逐点除法,类似于(1)。这是因为不需要反转,并且它的效果是在整个基函数的支持范围内集成的,我们将在第三节B部分中介绍。为了强调这个观点,我们采用符号
它可以被解释为相对于与累积函数相关的勒贝格-斯蒂尔杰斯度量的期望,其中的导数,因此是相应的pdf,
我们的实验观测收集在有限集合中,card() = N < \infty独立同分布(i.i.d.)的的实现。我们的目标是从这些数据中恢复

B. 密度的样条参数化

我们通过系数参数化我们的估计器,通过将其表示为多维样条的指数族来实现
其中
对于。(当时,由域的BCs指定,见(38)的周期性示例)。符号分别代表逐元素除法和乘法。注意,为了简洁,我们省略了的依赖性。一个人决定一个来选择网格相对于域的细度。通过膨胀B样条基来各向异性地膨胀函数,从而膨胀函数本身。的作用是通过膨胀B样条基来实现的,因此,函数也是如此。这种膨胀保留了BC的类型。
张量样条基是从一维B样条以可分离的方式构建的
B样条的度数为n,支持大小为n+1,逼近阶数为n+1。几个度数的封闭形式表达式可以在[35]中找到。在下文中,我们为了简化符号省略了度数。
通过在M上移动基,指数内的B样条展开可以表示任何相同度数n和具有相同(均匀)结点的样条,根据[36]中的理论。例如,具有整数结点的分段线性样条由唯一表征,其中,并且来自三角函数。对于任意的,结点在处,得到的函数是样条的各向异性膨胀,这也是样条。
虽然标准logsplines定义在具有非均匀结点的整个实数线上,但表达式(5)定义在具有一般BCs的多维域上,但在均匀网格上。一个优点是B样条展开在数值上是高效的,因为它们在支持和准确性之间有良好的权衡[37]。许多相关操作可以写成这些展开下的卷积[35]。它们也非常适合处理适当BCs的有限域。通过在这些BCs下对系数进行卷积操作来纳入这些操作,这对于周期性域特别方便,因为它们可以使用FFT计算。这种效率的一个后果是,可以在不损失准确性的情况下快速评估密度。例如,可以通过卷积在均匀网格上精确计算评估。这对于成像中的迭代方法来说是理想的,因为可能需要重复评估。指数保证了(5)是非负的。在第三节E部分,我们进一步评论了我们的(正则化的)DE方法与标准logsplines和自由结点概念的关系。
为了说明在网格上使用样条的优势,让我们定义一个采样算子,它接受一个函数并在网格上评估它。这个网格是从原始的M=M_1放大到一个更细的尺度s。我们写作
因此,的值域是。采样算子允许我们将密度估计的评估写成卷积
在d维中,具有对应于X的BCs。垂直箭头下标在中指的是通过用(s-1)零扩展来上采样c,使得如果,否则。离散卷积核在(9)中是可分离的,具有
这是对应于评估n阶导数的B样条滤波器,其中n = (n_k),n_k 。这个符号对于Hessian基础的正则化很有用,因为导数的评估也可以写成卷积形式。

C. 似然作为数据保真度

旨在从测量中恢复,我们针对系数制定一个优化问题。数据保真度基于观察到的i.i.d.样本x X_N的似然,由下式给出
特别是,我们使用相应的对数似然
为了整洁,我们经常省略密度对c的依赖性,简单地写作而不是
对数似然具有导数
对于每个,k M。它还具有大小为的Hessian ,其元素为
对数似然的Hessian是负定的,因为
对于所有,因为,且最内层项不能为零,因为是线性独立的。因此,对数似然是严格凹的。因此,(12)在参数化的系数处有一个唯一的最大值。这里,最大似然估计的存在对于非退化情况(当有足够的数据且常数解不是更好时)是保证的。
在(3)中通过两种方式约束结果密度。第一种是通过使用B样条基,因为单个样本x对单个的影响贯穿整个基函数的支持。第二种是通过先验灵敏度,这根据检测概率隐式地加权密度。由于选择等同于选择一个度量,该方法保留了原始1D logsplines[30]的许多有趣属性。
注意,在优化中的作用仅通过(3)或(12)-(13)中的归一化发挥,因为如果观测到的未归一化,导数将独立于。梯度由局部数据点和局部(在基支持上)对归一化积分的贡献相对于灵敏度驱动。框架接受(无需更改)任何非负且可以计算(13)中积分的灵敏度函数;它不需要(乘法)逆(参见(1))。当灵敏度平滑或更好时,以与(5)相同的形式表达时,一些计算可以加速(或更准确地执行)。

D. 基于近端的底层样条的正则化

为了进一步约束问题,我们在对数似然(12)中添加了一个不可微的正则化项R,结果是优化问题:
其中控制数据保真度(似然)和先验知识(正则化器)之间的平衡。我们直接在样条表述上强制所需的先验行为,作为R(),因为实验发现正则化对数密度通常是一种捕捉数据多模态性的良好方法[29]、[38]。

E. Hessian的核范数用于结点稀疏性

我们选择的正则化基于对数密度的Hessian(H)的Schatten p-范数,其中
是矩阵A的Schatten p-范数,其向量包含奇异值,p
我们为一般p开发了理论和软件。然而,在实践中,我们选择p = 1的核范数。它也被称为迹范数,对于正方形半正定矩阵,。由于核范数是秩函数的凸包,p = 1在这里作为促进具有一个主要奇异值或主曲率的低秩Hessians的凸替代。核范数也由于它是奇异值向量的范数而对等距变换不变。
我们的选择鼓励了在信息缺失时样条的仿射性,这转化为密度的指数行为。这种行为同时解决了标准logsplines的两个挑战。
  1. 尾部行为:正则化鼓励了在密度尾部的线性行为,即使在几乎没有样本的情况下也是如此。这种属性已知可以减少边界处的方差[39]。在传统的logspline拟合中,必须通过混合不同度数的样条来特别强制执行[31]。
  2. 带宽自适应性:类似于KDE中的自适应性,结点放置和结点删除是典型的logspline拟合中的重要步骤,涉及优化某些信息标准[31]、[32]。对于d = 1,我们的方法在B样条线性时自动停用/消除最不相关的结点——诱导结点稀疏性。对于更高的d和任意度数,它促进了Hessian的稀疏性和数据稀缺时对数的线性(图2)。反过来,这种稀疏性减少了网格粗糙度的影响(图3)。见附录A,在线可用,有关正则化效果的更理论描述及其与自由结点的关系。
实际上,我们通过在网格上评估Hessian和混合Sp-ℓ1范数来近似R
其中。在第四节中,我们将展示如何有效地评估(20)的最后一项,以及它的共轭和近端算子。

IV. 提出框架的优化

函数J从(− log(L_N))继承了严格凸性,因为R是凸的。因此,我们考虑使用加速近端梯度算法来解决(16)中的问题[40]、[41]。对数似然的梯度是从(13)计算的。正则化项R的近端算子也可以有效计算,如第四节C部分详细说明。
尽管负对数似然(− log(L_N))是两次连续可微的,但它只是局部梯度Lipschitz和局部强凸的。实际上,除非底层域是紧凑的,否则其二阶导数的范数不能被上界或(正)下界限制,因为我们是在一个指数族上工作。即使对于紧凑域,最坏情况的界限也会在其他地方太松。其结果是,设置近端梯度算法的全局下降步长大小是低效的。它需要在每一步适应(见图4中局部Lipschitz常数的界限)。

A. 自适应Lipschitz步长

为了利用加速算法的收敛速率,我们基于Lipschitz界限设计了一种自适应策略。这种方法的吸引力基于观察到Hessian(14)中的张量
可以写成外积,其中
是一个值,结果是梯度(13)的重要成分。剩余项
源自(14)中的,通常可能有多个非零奇异值,但总是带状的、对称的和正定的。
我们通过用Hessian来界定梯度g的最佳Lipschitz常数Lip{g}来制定自适应策略
最右边的项通过利用外积来计算,导致
其中下标“F”和“2”分别表示Frobenius和谱范数。为了界定,我们使用Gershgorin圆定理,因为我们通常知道D通常是对角占优的。这产生了一个界限定义为
注意,最多有个非零元素,并且由于D的对称性,只需要计算个。零元素的数量很大是由于B样条基的小支持。更一般的界限是
其中是大小为的m周围的掩模。我们在(26)是较松的界限的罕见情况下使用(27)。为了设置近端梯度下降的步长,我们使用界限BLip(c)的倒数

B. Hessian-Schatten范数的滤波器

为了计算(20),我们利用了在网格上的样条的便利性。在所有网格点上评估对数密度估计器的Hessian涉及个卷积。Hessian在网格点m M的第n个分量可以表示为
滤波器如第三节B部分所述。我们澄清Hessian的第n个条目n = (n_1, n_2) 对应于的偏导数。因此,一次卷积就足以在所有网格点上计算第n个Hessian元素。它们是类似于中心有限差分核,如[1, -2, 1]的滤波器的张量组合。一维B样条滤波器的递归表达式可以在[42]中找到。

C. Hessian-Schatten范数的近端

我们通过解决相应的最小化问题来计算(20),(29)的近端算子
我们使用基于对偶公式最大化的梯度投影算法来迭代求解,具体方法见[43]。实施它,我们需要三个定制的组成部分:算子在Sp-ℓ1范数内的共轭,对偶空间Sq-ℓ∞范数单位球的投影,以及合适的步长。
  1. 共轭:我们通过执行共轭定义来派生Hessian在网格上的评估算子的共轭。它是在共轭BCs下与滤波器共轭的卷积之和,
    其中,上标*表示共轭,下标表示在维度上进行卷积。原则上,共轭的脉冲响应将是原始滤波器沿适当轴的反转脉冲响应,但实际上它就是原始滤波器的脉冲响应,这是由于对称性考虑。
  2. 投影:将矩阵投影到Schatten q-范数的单位球上,需要将矩阵分解为奇异值并将其投影到q-范数的单位球上[44, Proposition 1 and Equation (34)]。然后使用原始奇异向量重新组合矩阵。对于小d,存在封闭形式的解,这导致了快速实现。
  3. 步长:为了最大化(30)中的对偶函数,我们需要选择梯度上升算法的步长。我们根据对偶梯度的Lipschitz常数的上界来设置它。通过将[44, Proposition 2]中的论证推广到任意维度和样条滤波器,我们可以展示是这样的界限。关键在于与的卷积被限制在1以内。
我们提前停止了基于对偶的算法的迭代,因为即使近端算子的应用近似粗糙,外层近端梯度算法的收敛也会发生[45]。

D. 算法和实现

有了梯度(13)、近端算子(30)和自适应Lipschitz界限(28),我们设置了一个带有重启方案的加速近端梯度下降算法,概述在算法1中。我们还推导出卷积表达式来加速积分的归一化和密度评估的计算(附录B-C,在线可用)。

在我们的算法1实现中,我们选择了相对收敛标准
容忍度为。加速被选择为标准
以及
重启方案被实现为
算法以负系数初始化。
产生的库将在线提供。1它完全基于NumPy和SciPy;它还通过它们的CuPy对应物提供GPU支持。

V. 实验

我们在标准理论分布和PET正弦图上测试了我们的框架。我们通过将得到的估计器与真实密度进行比较来评估每种方法的准确性,根据
这涉及到在一些精细网格上的均方误差(MSE)。在这里,我们将我们的框架称为RDS(正则化密度样条)。
为了模拟来自观测密度的样本x ν,我们首先从感兴趣的底层密度中取样本p 。然后我们根据概率ξ/ max(ξ)对它们进行稀释。更准确地说,如果ξ(p)/ max(ξ) > u,其中u U(0, 1)来自区间(0,1)上的均匀分布,我们就保留样本p。这被称为拒绝采样。

A. 标准分布

我们考虑了两种不同的分布ν作为标准分布:UGL和GG。UGL是由均匀分布、高斯分布和拉普拉斯分布(图5,最左)的和组成的复合分布。我们通过考虑每个维度上的独立拉普拉斯分布将其扩展到nD。GG是两个高斯的和。(在线1上找到所有参数的确切列表。)分布的域被视为周期性的。我们测试了维度d = 2和d = 3。灵敏度图被选择为ξ(x) = 1和ξ(x) = sin²(x₁/T + φ) + ϵ,其中φ是相位,T是周期,ϵ = 10⁻³(图5,中间)。我们称这些灵敏度图为sU和sS。

我们选择了网格M = ∏_{k=1}^d {0, 1, ..., N_k - 1},并根据周期性域设置了基
注意“缺少”基在Nk,因为它是0处的基。我们考虑的样本数量范围从10²到10⁶。

B. 初步分析

  1. 步长自适应性:我们验证了在优化问题的每次迭代中自适应Lipschitz常数的重要性。为此,我们比较了两种方案的算法收敛性:一个在每次迭代中更新Lipschitz“常数”,我们称之为BLip(ctemp,k);另一个在c1重新初始化后在整个迭代中保持不变,我们称之为BLip(c1)。在这两种情况下,界限都是根据(24)计算的。在图4中,可以看到BLip(ctemp,k)的演变导致了快速收敛。相反,当使用非变化的BLip(c1)时,我们观察到有时优化会立即对某些初始化发散,有时在一些下降后会出现强烈的振荡。
  2. 度数:在我们的实验中,我们观察到从0度样条移动到1度样条时改进最大。3度样条与它们的高计算成本相比并没有增加太多准确性(图6)。这与我们之前使用样条的经验一致。因此,我们将n = 1设置为以后使用。
  3. 正则化效果:我们首先检查线性样条沿轴诱导的稀疏性。我们使用d = 2的拉普拉斯分布在这个测试中。我们观察到,最小化核范数导致了一个稀疏的Hessian,在这里它作为沿第一轴的结点稀疏性的替代(图2)。
在图3中,我们观察了RDS估计对正则化参数的敏感性。对于给定的样本数量,最优λ并不依赖于测试范围内的网格大小(图3,垂直线)。样本数量以可预测的方式影响最优λ:将样本数量乘以10的因子需要将最优λ乘以2的因子(图3,垂直线)。超过一定的网格大小时(相对于样本数量),MSE对λ的依赖性在最小错误处达到平台(图3)。这种相变行为让人想起稀疏优化,并在选择最优λ时留有余地。我们还观察到,随着网格规模或样本数量的增加,MSE的减少回报递减。

C. 标准分布中的误差和速度

我们将RDS与KDE和HE进行了比较。HE的带宽根据Freedman-Diaconis和Sturges估计量中的最大值进行调整,而KDE的带宽则遵循Scott的规则。KDE和HE的实现取自SciPy。我们为KDE的核选择了高斯分布(径向基函数)。

1) 定量评估

我们评估了样本数量函数的MSE。RDS在所有情况下都优于HE和KDE,并且在整个样本范围内(图7,UGL,d ∈ {2, 3},sU和sS)都是如此。对于成像方式相关的样本数量(≥ 10^5),改进达到了一个数量级(图7,最左)。RDS对sS灵敏度的鲁棒性仅在样本范围中损失了大约0.1 dB,与sU相比(图7,中间)。KDE也相当鲁棒(损失了0.6 dB)。KDE中的不规则权重通过过度平滑得到了补偿。这是以估计器的定性价值为代价的,因为模式混合在一起。相比之下,HE变得更糟——其MSE大约下降了4 dB,但仍然保持清晰。

在3D中,与RDS相比的相对改进甚至更大,对于10^6个样本达到了15 dB的最大值(图7,最右)。从2D到3D,由于维度的诅咒,HE与KDE相比变得更糟,因为超立方体的箱子很快就变空了。
为了进一步探索该方法对λ的鲁棒性,我们比较了两个版本的RDS。一个在每个样本数量中调整λ,根据五个值(对应五个数量级)中的最好MSE。另一个则将λ固定在10^2个样本找到的最优值。它们的表现相似(图7)。

2) 定性评估

为了更好地支持KDE的高斯核,我们首先在5·10^3 个样本下评估GG(图8)。我们观察到RDS很好地适应了不同底层方差区域,即使在不均匀的灵敏度下也是如此。它还很好地补偿了数据稀缺的不那么敏感区域。相反,KDE通过过度平滑大部分领域来补偿样本的缺乏,而HE则无法纠正灵敏度。在UGL分布的估计中也可以观察到类似的结果(图9和10,分别为5·10^3和5·10^5个样本)。

3) 速度

RDS比KDE的另一个优势是减少了计算时间。RDS在2D中优化大约需要1秒,评估需要10^-3 秒。在我们的测试中,这与样本数量无关(图11)。KDE的优化成本与其评估成本相比是微不足道的。一次评估从10^-2 秒到10^2  秒不等,这取决于我们考虑的样本范围。这在只有10^ 4个样本时就达到了RDS的优化时间。HE的计算时间也随着样本数量的增加而增加,但总体上保持较低。因此,RDS在存在许多样本或需要重复评估时特别有利。例如,我们预计在PET这样的成像方式中将节省大量的计算资源,因为它们通常处理超过10^6个样本。

在3D中的结果与2D一致。RDS仍然与样本数量无关。RDS在3D中的优化比2D慢了一个数量级,而评估只慢了1.2倍。KDE评估比其2D对应物慢了一个数量级。

VI. 成像应用

超分辨率显微镜和最近的单光子发射计算机断层扫描(SPECT)和PET都使用低光子计数。电子激发和核衰变的随机性导致光子发射被描述为具有空间变化强度的不均匀泊松过程。一旦最终检测到N个光子,我们认为问题归结为估计pdf 的源点,该源点生成点云,其中。我们将用PET示例来说明这些想法。

A. PET背景

PET根据与样本中放射性分布成比例的成对光子发射来重建图像。这些成对光子发射的线称为响应线(LOR),构成测量[47]。扫描仪对LORs的灵敏度通过对参考扫描进行逐点除法来校正。事实证明,在一些最先进的扫描仪中,探测器非常小[3]、[6]、[48]、[49]。因此,大多数LOR累积了一对光子(约0.95%)或根本没有(约99%)而不是几百对(图1(a))[3]、[6]、[48]、[49]。因此,我们认为典型的PET重分组中的插值可以更好地从DE的角度来考虑。这些扫描仪的不均匀灵敏度的校正是另一个挑战,因为它可能跨越一个数量级(图1(b))。
测量集合称为断层扫描正弦图。在那里,我们将LORs的角坐标和“距离”坐标的对视为样本x = (θ, s)设置在X = [0, π) × [−ρ, ρ]上,其中ρ是扫描仪的视野半径——没有LOR可以在其外存在。这也意味着对于所有的θ ∈ [0, π),。角BC是周期性的,距离BC是恒定的。每个样本对应于连接在大约相同的时间但在不同探测器处检测到的一对反平行光子的LOR。

B. PET实验

我们在实验小动物PET扫描仪的背景下测试了RDS与KDE和HE。该扫描仪的探测器非常小[3]、[6]、[48]。我们还评估了从正弦图重建的图像的质量。KDE和HE的带宽选择如第五节C部分所述。
典型的PET重分组的函数插值方法(与DE相对)被省略了,因为它不适用于这种稀疏数据区域。DE更合适,因为与探测器的大小相比,样本数量很小(图1(a))。
我们首先评估了扫描仪的灵敏度ξ。在正弦图域X中的测量表明,检测LOR的机会是非常不均匀的(图1(b))。我们记录了最高和最低灵敏度之间高达一个数量级的差异。这对结果正弦图有相当大的影响(图1(c)-(d))。
然后我们假设放射性物质(放射性示踪剂)的底层浓度是四个高斯的和(数据未显示,见[50, 图3-4])。原因有两个:一是与KDE进行比较,二是便于解释正弦图。由于PET是一种断层成像方式,因此在结果正弦图上保留了高斯性。PET采集模拟了相当于10^2到10^6个样本的扫描时间。结果MSE的行为几乎与图7中的相同。在PET相关的10^6个样本的标记处,RDS的性能比KDE和HE高出一个数量级以上。我们的定性观察结果也与之前对标准分布的分析相同。
RDS对正弦图估计的改进转化为在图像域上使用滤波反投影(FBP)重建时更好的MSE。这是在Derenzo幻影中的情况,这是PET系统最常见的基准(图12和1(d))。RDS成功减少了FBP引起的条纹伪影,而没有受到KDE的过度平滑的影响。由于FBP要求在网格上执行评估,RDS的卷积属性也提供了性能优势。

  1. 数据和幻影:我们接下来在更实际的条件下测试了算法。
我们首先使用了[51]提出的超高分辨率PET幻影(Brain Phantom)。我们在[49]中的PET扫描仪中模拟了其发射和采集。(有关PET模拟的更多详细信息和评论,请参阅在线附录D。)我们使用上述描述的三种DE方法重新采样了正弦图。然后使用FBP(根据NEMA标准)从得到的正弦图中重建PET图像。我们还使用总变分正则化重建了PET图像。重新采样的正弦图及其相应的重建如图13、14和15所示。

我们还测试了我们的方法在Digimouse[52]以及在广泛使用的PET-MR扫描仪Siemens Biograph mMR 3T上获得的Amyloid真实数据集上[53]、[54](有关更多详细信息,请参阅在线附录D)。它们的正弦图和重建也图13、14和15中。
我们对所有三个实验(Brain Phantom、Digimouse、Amyloid)的观察结果相似。正弦图中的低灵敏度区域对HE来说是有问题的,部分原因是由于灵敏度反转引入的稳定性问题。在这些情况下,KDE过度平滑了图像,试图补偿低灵敏度。这忽视了潜在的尖锐区域,这些区域可能存在于底层正弦图中。当低灵敏度区域(或正弦图值低的区域)具有不同的尺度或锐度时,这个问题就更加复杂了。例如,参见Amyloid数据集中的“间隙”区域网格或Digimouse中的区域。相比之下,RDS对底层正弦图和灵敏度的不均匀性表现得相对较好。这在断层成像的情况下是显而易见的,即使X射线变换(导致正弦图)是平滑的[55]、[56]。此外,RDS还更好地保持了原始对比度。这在PET中尤其重要,因为对比度描绘了放射性示踪剂吸收不同的区域。在模拟和真实实验中都可以观察到类似的效果。在使用总变分正则化重建(与FBP相比)时,所有三种方法都相对于背景获得了更好的对比度。

VII. 结论和讨论

正则化密度样条(RDS)补偿了不均匀灵敏度,并且随着维度的增加而扩展。该方法对正则化参数的选择是鲁棒的。这是密度估计(DE)的一个好属性,因为它扮演了带宽选择的角色。正则化的作用也得到了很好的理论支持:它在不变的框架下诱导了Hessian的稀疏性。优化和评估时间与样本数量无关,评估可以完全并行化。RDS的实现可通过支持GPU的库获得。
所有这些特性使RDS成为多种应用的良好候选者,特别是成像应用。我们验证了具有小探测器的现代PET扫描仪可以通过DE来处理,RDS是一个特别好的选择。这种方法也适用于低分辨率扫描仪,如果考虑了量化[57]。随着纳米技术推动探测器变得更小,我们预计其他成像方式将进入加权统计采样的领域。这是因为定位变得越来越准确,但信号功率可能保持不变——受到剂量或曝光的限制。
未来的研究方向将包括进一步调查正则化项与样条稀疏性之间的联系。Box样条在这方面有很大的潜力,因为它们是分段线性的,但可能需要增加计算和实现成本。其他工作将集中在根据RDS估计的不确定性定制数据项以进行图像重建。

声明

本文内容为论文学习收获分享,受限于知识能力,本文对原文的理解可能存在偏差,最终内容以原论文为准。本文信息旨在传播和学术交流,其内容由作者负责,不代表本号观点。文中作品文字、图片等如涉及内容、版权和其他问题,请及时与我们联系,我们将在第一时间回复并处理。

#论  文  推  广#

 让你的论文工作被更多人看到 


你是否有这样的苦恼:自己辛苦的论文工作,几乎没有任何的引用。为什么会这样?主要是自己的工作没有被更多的人了解。


计算机书童为各位推广自己的论文搭建一个平台,让更多的人了解自己的工作,同时促使不同背景、不同方向的学者和学术灵感相互碰撞,迸发出更多的可能性。 计算机书童 鼓励高校实验室或个人,在我们的平台上分享自己论文的介绍、解读等。


稿件基本要求:

• 文章确系个人论文的解读,未曾在公众号平台标记原创发表, 

• 稿件建议以 markdown 格式撰写,文中配图要求图片清晰,无版权问题


投稿通道:

• 添加小编微信协商投稿事宜,备注:姓名-投稿

△长按添加 PaperEveryday 小编


PaperEveryday
为大家分享计算机和机器人领域顶级期刊
 最新文章