CN103608825A - 确定迭代影像重建中的像元映射的方法 - Google Patents

Abstract

一种使用像元方法迭代重建包含噪声的信号的方法从用以在迭代中更新影像的变量中确定像元映射。更新变量是基于一个最优评价函数的，并且在该迭代期间将该更新变量平滑化。在该迭代结束时还可以使用在该迭代期间确定的该像元映射任选地进一步将所更新的影像平滑化。

Description

确定迭代影像重建中的像元映射的方法

相关申请

本申请要求在2011年4月15日提交的美国临时申请号61/476,012的优先权，该申请的披露内容通过引用以其全部结合于此。

发明领域

本发明涉及一种用于信号重建和增强的方法，并且更确切地说，涉及一种自适应地减少输入对象中的噪声的方法。

发明背景

对包含在一个信号中的数据的最优提取需要消除多个缺陷例如噪声以及多个仪器限制。试图进行最优化提取和重建数据的主要领域是影像增强领域。即使在可以将多个仪器制造得实质上无噪声时，与有限的空间、频谱、或时间分辨率相关的多个仪器签名仍然存在。在此点上，需要影像重建来消除这些仪器签名。影像增强的多个应用和多个噪声源以及会负面地影响数据提取的多个其他因素覆盖一个很大的范围，包括：天文观察和行星探测，其中多个源可以是微弱的并且大气干扰引入噪声和失真；军事和安全监视，其中光可以是多个目标的较低且快速的移动，从而导致对比度差且模糊；医学成像，其经常清晰度不够；以及视频影像，其中传输和仪器限制以及对实时响应的需要会负面地影响影像清晰度和细节。

已经开发了数字影像处理来提供对由多种传感器收集的模糊且有噪声的数据的高质量、鲁棒性重建。该领域存在是因为创建产生未被测量噪声污染的清晰度一致的多个图像的多个成像仪器是不可能的。然而，从自现实世界的仪器获得的不理想的数据数学上重建底图像是可能的，这样使得可以将现存但隐藏在数据中的信息提取得具有较少模糊和噪声。许多此类方法利用如下过程：构建一个预测模型且将其与这些数据比较来评估该模型的数据拟合的准确度。

贯穿此书面描述，“数据”是指任何所测量的量值，通过影像重建过程从该所测量的量值估计出一个未知的“影像”。术语“影像”表示要么估计的解要么产生观察到的数据的真实的底图像。该讨论通常清楚地表明适用的上下文；在可能的歧义的情况下，“影像模型”用以表示估计的解。注意，数据和影像不需要是相似的，并且甚至可以具有不同的维度，例如，断层重建试图从投射的2D数据确定一个3D影像。影像的一个替代术语是“对象”，其传达了模型可以比影像更概括的思想。这两个术语此后用作同义词。

统计学家长期以来试图限制用以拟合数据的多个参数，以便提高该拟合的解释和预测准确度。标准的技术是：子集选择，其中这些参数中的一些被发现是不重要的并且被从该拟合中排除（例如，米勒（Miller）2002）；以及岭回归，其中这些参数的值通过向用于该拟合中的评价函数添加一个正则项而受到约束（例如，吉洪诺夫（Tikhonov），1963）。蒂施莱尼（Tibshirani）（1996）将这两种方法组合在一种技术中，该技术称为最小绝对缩减和选择算子（LASSO）。

限制一个拟合中的参数的数目的需要对欠定的或不充分地确定的问题而言是必要的，其中参数的数目大于数据点的数目或与其相当。如已经由最小二乘（高斯（Gauss），1809）和最大似然（费希尔（Fisher），1912，1922）的统计方法的创始人强调的，这些方法仅仅在渐近极限内有效，其中数据点的数目大大地超过拟合的参数的数目。远离此渐近极限，噪声被拟合为信号，并且该拟合失去其解释和预测能力。

佩特（Puetter）等人（2005）回顾了当前使用中的许多重建算法，包括迭代影像重建方法，其将多个影像模型迭代地拟合成数据。即使在其较早终止来避免过度拟合时，被设计来收敛到最大似然解的许多现有技术迭代方案缓慢地收敛。收敛可以通过使用评价函数相对于变量的二阶偏导数的海赛矩阵而更快速地达成（海赛（Hesse），1876）。不幸地，在处理在影像重建过程中经常遇到的大范围问题时，例如，在该影像包含具有显著发射的大量的像素的情况下，即使使用此方法，海赛矩阵也太大而不能被计算。在此类情况下，矩阵元素可以达到数万亿。这样大小的矩阵完全不能被现代计算机处理，或者甚至存储在存储器中。

限制欠定的或不充分地确定的问题中的参数的数目（如现今所谓的，稀疏性）的希望是基于更抽象的最小复杂性的概念（索洛莫诺夫（Solomonoff），1964；柯尔莫戈洛夫（Kolmogorov），1965；蔡廷（Chaitin），1966），其追溯到提倡假说简约性的奥卡姆的威廉的中世纪工作。简单来说，在其他事物相等的情况下，一个较简单的解释比一个较复杂的解释更好。

例如通过向评价函数添加一个l₀范数正则项来找到一个稀疏的解是一个N-P困难问题，其中计算工作量比许多参数中的任何多项式增加得更迅速。这已经导致用一个l₁范数来替换该l₀范数（陈、多纳霍（Donoho）和桑德斯（Saunders），1999），这样使得该拟合相对于这些参数的最优化变成一个可解的凸问题。坎德丝（Candés）、龙伯格（Romberg）和陶（2004）更进一步并且示出了如何在一个不连贯的条件下随机地减少该拟合所需的数据量，一种称为压缩感测的技术。多纳霍（2006）示出了在一个类似的不连贯的条件下在用于参数化的基函数当中最小l₁范数解也是最稀疏的解。

l₁范数方法的缺点有两部分。首先，许多重要的问题完全不满足不连贯条件并且不适用于l₁范数方法。其次，即使在其不满足该不连贯条件时，大范围问题也需要过度的计算工作量。因此，尽管这些问题是凸问题且原则上可解，其实际上不能应用于具有数百万个或更多参数的当今问题，这是一个同样使传统的统计方法苦恼的问题。多纳霍等人（2006）讨论了如何将随机性更有效地应用于这些大范围问题而不使用一个l₁范数。

像元方法是一种获得基于像素的数据的多个复杂性最小的解的有效技术，包括大范围问题，无需严格的稀疏性并且没有利用l₁范数的方法所需的不连贯条件（参见例如碧娜（Pina）和佩特，1993；佩特和亚希勒（Yahil），1999；佩特、格斯内尔（Gosnell）和亚希勒，2005；以及美国专利5,912,993、6,353,688、6,490,374、6,895,125、6,993,204、7,863,574、7,928,727、8,014,580、8,026,846、8,058,601、8,058,625、8,086,011、8,090,179、8,094,898、8,103,487，这些专利的披露内容通过引用结合于此）。因此，像元方法对于大范围欠定的或不充分地确定的反问题例如影像重建或频谱分析而言是有用的。最小复杂性是通过在每个像素位置处通过一个内核库当中的最宽的内核自适应地平滑化而达成的，如此使得用此以及所有较窄的内核的平滑化提供对考虑中的像素的数据足迹的充分拟合。指明在每个像素处使用哪个内核的映射叫做像元映射。

在其当前形式中，一个像元重建由三个步骤组成。首先，其重建一个没有任何像元约束的“伪影像”。其次，此伪影像用以确定像元映射。第三，通过由该像元映射引导的一个受约束的重建获得最终影像。可以将步骤二和三重复许多次，但实际上这通常是不必要的，条件是在第一步骤中获得了一个合理的伪影像。参见佩特和亚希勒（1999）、佩特等人（2005）以及美国专利号5,912,993、6,353,688、6,490,374、6,895,125、6,993,204、7,863,574、7,928,727、8,014,580、8,058,601、8,058,625、8,086,011、8,090,179、8,094,898以及8,103,487来获得对像元方法以及其应用的更完整的讨论。

图2图示了具有一个成像检测器210和一个像元重建单元220的一个通用成像系统200。重建是基于一种像元方法的，该方法使用一个像元映射P，其与一个像元重建算法230互动。该像元方法是指一种通过向每个对象点指派一个形状或容量作为像元平滑化的基础来将对象空间中的每个点（此后，一个“对象点”）平滑化的方法。该对象空间是定义影像重建结果的空间并且对应于使用成像系统200所成像的一个域。（应注意，“影像空间”是与“对象空间”同义的一个术语，并且此后可交换地使用这两个术语。）通过用成像检测器210测量的多个数据点来给出一个对应的数据空间。

该像元方法提供对来自数据空间中的一个所测量的数据集d的对象空间中的一个影像对象I的高质量重建。作为一种空间上自适应的重建方法，该像元方法将一个数据激发的平滑化操作应用于每个对象点。在这样做时，该像元方法在向每个对象点指派一个像元内核函数时使用最小复杂性的原则，这是平滑化操作的基础。在像元重建单元220内，像元映射P定义将这些像元内核函数中的哪一个指派给这些对象点中的每一个。

在成像系统200中，成像检测器210检测所测量的数据集d并且将其传递给像元重建单元220。像元重建单元220使用多个特别地适配的像元重建算法230来将所获得的数据集d重建成一个影像对象I。在这样做时，像元重建算法230使用一个系统矩阵H来描述成像系统200的多个性质，并且通过调整数据模型估计一个迭代地改进的影像对象，该数据模型是影像对象I的基础。使用熟知的渲染技术将影像对象I例如显示在一个显示器240上。

针对每个对象点，像元映射P提供一个基于最小复杂性方法确定的像元内核函数。此像元内核函数用于一个应用于对象空间中的像元平滑化操作中。

该像元方法还可以允许以超分辨率形式的压缩感测，从而利用影像的非负性以及最小复杂性来重建一个影像，该影像具有比用来获得数据的那些像素更精细的像素。因为非负性和最小复杂性的附加条件（例如，佩特等人，2005），这并没有因为奈奎斯特（Nyquist）（1928）和香农（Shanon）（1949）而违反取样定理。在数据中被截短的超过衍射极限的空间频率可以在该影像中类似地重建。

图3图示了像元方法的一个示例性流程。像元平滑化顺序地应用于一个标准的重建算法。

使用一个标准的重建算法将输入影像拟合成一个所测量的数据集d（步骤300）。根据像元内核算子K的上文讨论的使用，该影像的所产生的估计叫做一个伪影像。使用该伪影像以及所测量的数据集d确定像元映射P（步骤310）。该伪影像还是像元平滑化操作（步骤320）的最初对象。在该像元平滑化操作（步骤320）期间，该伪影像的每个对象点通过一个像元内核函数平滑化。（在现有的像元方法的一些变体中，像元映射还可以在每个迭代中通过从更新的影像计算来更新。）

迭代的影像重建方法将多个影像模型迭代地拟合成所测量的数据，并且因此最小化噪声对最终影像的影响。一个重建算法的结果是根据该算法的多个规则拟合成所测量的数据集d的一个近似的影像。

在像元方法中，一个近似的影像可以用作像元平滑化的一个输入对象，以用于像元重建以及用于确定像元映射。

像元方法包括在对象空间中的每个点处对多个最广的可能像元内核函数的搜寻，这些像元内核函数一起支持一个对象向所测量的数据集d的一个充分拟合。具体地，该像元映射向每个对象点指派一个特定的像元内核函数。

首先计算一个伪影像，接着可以从该伪影像确定像元映射的步骤可能存在多个弊端。此过程需要更多的计算并且承受着将多个伪像引入该伪影像中的风险，这可以使该像元映射以及因此最终的影像重建产生偏差。此外，在从对象空间向数据空间的变换是非局部的情况下（巴特纳格尔（Bhatnagar）和康威尔（Cornwell），2004），像元映射的确定效果较差。例如，在干涉法以及磁共振成像过程中，数据是影像（加上噪声）的多个傅立叶变换，其中每个傅立叶波（该影像的基函数）在整个影像之上传播。另一个示例是在断层成像中的部分非局部的变换。这些数据是一个3D影像（加上噪声）的多个2D投影，一个横向于投影方向是局部的但沿着投影方向是非局部的变换。

鉴于前述内容，需要一种改进的方法来确定该像元方法内的像元映射。

发明概述

根据此中的发明，通过在迭代期间计算以及更新像元映射避免三阶段像元方法的上述多个困难。改进的方法从一个变量确定该像元映射，该变量用以在该迭代中更新影像，即，一个“更新变量”，并且在该迭代期间将此更新变量平滑化。在该迭代结束时还经常使用在该迭代期间确定的该像元映射进一步将更新的影像平滑化。相比之下，现有的像元方法从已经更新之后的影像确定像元映射并且继续用彼像元映射将该影像平滑化。

根据本发明，更新变量取决于重建方法，但典型地是一个评价函数的梯度或者一个乘法更新因子（例如，佩特等人，2005）。此更新变量通过多个像元内核平滑化，并且在每个影像位置处选择的内核是最宽的内核，其中该更新变量中因为平滑化的改变的平方与该更新变量的方差之间的比率小于该内核以及所有较窄的内核的一个预定的阈值。相同的内核接着通常还用以在继续进行到下一个迭代之前在彼位置处将更新的影像平滑化。一个进一步的细化允许“多个插入的内核”。

在本发明的一个方面中，一种从自物理过程获得的数据集重建对象模型的方法，其中该数据集包含噪声，该方法包括：接收在一个数据空间中定义的该数据集；在一个对象空间中构建一个对象模型，其中，该对象模型包括多个对象点；设计该对象模型从对象空间向数据空间的一个变换，从而产生一个数据模型，其中，该变换对应于获得该数据集的物理过程；选择用于确定一个数据模型向该数据集的一个拟合的一个评价函数；基于该评价函数确定对象空间中的该对象模型的一个更新变量；将该更新变量平滑化以通过以下方法确定一个平滑化的更新变量：使该更新变量与多个像元内核中的每一个卷积；以及为该输入对象的每个对象点选择一个像元内核，该像元内核具有符合一个预定的最小标准的一个最大大小；通过在对应于所选的像元内核的该输入对象的每个对象点处指派多个索引来产生一个像元映射；以及基于该像元映射内的这些索引产生一个包括一个实质上降噪的对象模型的输出。在一个实施例中，使用一种共轭梯度法来确定该评价函数并且该更新变量是一个负梯度。

在本发明的另一个方面中，一种非瞬态计算机可读介质，具有附带的软件，该软件包括用于从自一个物理过程获得的一个数据集重建一个对象模型的多个指令，其中，该数据集包含噪声，并且这些指令包括：接收在一个数据空间中定义的该数据集；在一个对象空间中构建一个对象模型，其中，该对象模型包括多个对象点；设计该对象模型从对象空间向数据空间的一个变换，从而产生一个数据模型，其中，该变换对应于获得该数据集的物理过程；选择用于确定一个数据模型向该数据集的一个拟合的一个评价函数；基于该评价函数确定对象空间中的该对象模型的一个更新变量；将该更新变量平滑化以通过以下方法确定一个平滑化的更新变量：使该更新变量与多个像元内核中的每一个卷积；以及为该输入对象的每个对象点选择一个像元内核，该像元内核具有符合一个预定的最小标准的一个最大大小；以及通过在对应于所选的像元内核的该输入对象的每个对象点处指派多个索引来产生一个像元映射。在一个示例性实施例中，使用一种共轭梯度法来确定该评价函数并且该更新变量是一个负梯度。

附图简要说明

图1是可以实践本发明的一般计算环境的方块图。

图2是基于像元方法的一个现有技术成像系统的概述图。

图3是图示将现有技术方法用于像元重建的影像重建的一个示例的流程图。

图4是图示使用根据本发明的像元方法的影像重建的一个示例的流程图。

图5是图示使用像元方法的一个替代影像重建过程的一个示例的流程图。

图6示出了一个无噪声的、归一化的、质子密度（PD）SBD“真”的脑影像。

图7示出了图6的影像向数据空间的傅立叶变换的绝对值的对数（基数10）。

图8示出了具有随机（白）高斯噪声的图6的影像，该随机（白）高斯噪声添加到该影像，具有最亮的影像体素的3%的标准差。

图9示出了图8的影像向数据空间的傅立叶变换的绝对值的对数（基数10）。

图10示出了与图9相同的数据，其中高频率被截短，对应于限制到可用频率的三分之一的频率覆盖范围。图10中的数据包括向图11至图14中所示的影像的多个重建的输入。

图11示出了图10的数据的一个简单的逆傅立叶变换，揭示了两种类型的伪像：关联的（非白）噪声以及在多个锐边处的振影。

图12示出了图10的数据的一个像元重建的结果，其中像元因子为0.3。

图13示出了图10的数据的一个像元重建的结果，其中像元因子为0.5。

图14示出了图10的数据的多个像元重建的结果，其中像元因子为1.0。

发明详述

在描述本发明的多个方面之前，提供对可以实现本发明的一个合适的计算系统环境100（图1）的一个简要描述可能是有用的。计算系统环境100仅仅是一个合适的计算环境的一个示例，并且不旨在建议关于本发明的使用或功能性的范围的任何限制。也不应将计算环境100解释为具有与在示例性运行环境100中图示的多个部件中的任一者或组合有关的任何从属性或要求。

本发明是用许多其他一般目的或特殊目的计算系统环境或配置可运行的。可能适合于与本发明一起使用的多个熟知的计算系统、环境、和/或配置的多个示例包括（但不限于）个人计算机、服务器计算机、手持式或膝上型装置、多处理器系统、基于微处理器的系统、机顶盒、可编程的消费电子器件、网络PC、小型计算机、大型计算机、电话系统、包括上述系统或装置中的任一者的分布式计算环境，等等。

可以在多个计算机可执行的指令例如由一个计算机执行的多个程序模块的一般上下文中描述本发明。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、分量、数据结构等。本领域技术人员可以将此中的描述和/或多个图实现为多个计算机可执行的指令，这些指令可以在下文讨论的任何形式的计算机可读介质上体现。

本发明还可以在多个分布式计算环境中实践，其中通过多个远程处理装置执行多个任务，这些远程处理装置通过一个通信网络链接。在一个分布式计算环境中，多个程序模块可以定位在本地和远程计算机存储介质两者中，这些计算机存储介质包括存储器存储装置。

参见图1，一种用于实现本发明的示例性系统包括呈一个计算机110形式的一个通用计算装置。计算机110的多个部件可以包括（但不限于）一个处理单元120、一个系统存储器130、以及一个系统总线121，该系统总线将包括该系统存储器的不同系统部件耦合至处理单元120。系统总线121可以是若干类型的总线结构中的任何一种，包括存储器总线或存储器控制器、外围总线、以及使用多种总线架构中的任一种的局部总线。通过举例，并非限制，此类架构包括工业标准架构（ISA）总线、微通道架构（MCA）总线、增强型ISA（EISA）总线、视频电子标准协会（VESA）内部总线、以及外围部件互连（PCI）总线，也称为夹层（Mezzanine）总线。

计算机110典型地包括多种计算机可读介质。计算机可读介质可以是可以由计算机110访问的任何可用的介质并且包括易失性和非易失性介质、可移除和不可移除介质两者。通过举例，并非限制，计算机可读介质可以包括计算机存储介质和通信介质。计算机存储介质包括易失性和非易失性、可移除和不可移除介质两者，这些介质以用于存储信息例如计算机可读指令、数据结构、程序模块或其他数据的任何方法或技术实现。计算机存储介质包括（但不限于）RAM、ROM、EEPROM、闪存或其他存储器技术、CD-ROM、数字多用磁盘（DVD）或其他光盘存储器、磁盒、磁带、磁盘存储器或其他磁性存储装置、或可用以存储期望的信息并且可以由计算机110访问的任何其他介质。通信介质典型地体现计算机可读指令、数据结构、程序模块或一个已调数据信号例如载波或其他传送机构中的其他数据，并且包括任何信号输送介质。术语“已调制数据信号”意味着具有其特征集中的一个或多个或被改变的方式为将该信号中的信息编码的一个信号。通过举例，并非限制，通信介质包括有线介质，例如有线网络或直接有线连接，以及无线介质，例如声学、RF、红外线以及其他无线介质。上述中的任一者的组合也应该包括在计算机可读介质的范围内。

系统存储器130包括呈易失性和/或非易失性存储器形式的计算机存储介质，例如只读存储器（ROM）131以及随机访问存储器（RAM）132。基本输入/输出系统133（BIOS）典型地存储在ROM131中，该基本输入/输出系统包含多个基本例程，这些基本例程帮助例如在启动期间在计算机110内的多个元素之间传递信息。RAM132典型地包含处理单元120立即可访问和/或目前由该处理单元对其运行的数据和/或程序模块。通过举例，并非限制，图1图示了一个操作系统134、多个应用程序135、多个其他程序模块136、以及程序数据137。

计算机110还可以包括其他可移除/不可移除易失性/非易失性计算机存储介质。仅通过举例，图1图示了从不可移除、非易失性磁性介质读取或向该介质写入的一个硬磁盘驱动器141、从一个可移除、非易失性磁盘152读取或向该磁盘写入的一个磁盘驱动器151、以及从一个可移除、非易失性光盘156例如一个CD ROM或其他光学介质读取或向该光盘写入的一个光盘驱动器155。可用于示例性运行环境中的其他可移除/不可移除、易失性/非易失性计算机存储介质包括（但不限于）盒式磁带、闪存卡、数字多用磁盘、数字视频磁带、固态RAM、固态ROM，等等。硬磁盘驱动器141典型地通过一个不可移除的存储器接口例如接口140连接至系统总线121，并且磁盘驱动器151和光盘驱动器155典型地通过一个可移除的存储器接口例如接口150连接至系统总线121。

上文讨论并且在图1中图示的这些驱动器及其相关联的计算机存储介质为计算机110提供对计算机可读指令、数据结构、程序模块以及其他数据的存储。在图1中，例如，硬磁盘驱动器141图示为存储操作系统144、多个应用程序145、多个其他程序模块146以及程序数据147。注意，这些部件可以与操作系统134、多个应用程序135、多个其他程序模块136以及程序数据137要么相同要么不同。操作系统144、多个应用程序145、多个其他程序模块146、以及程序数据147在这里被给予不同的数字来说明在最低限度其为不同的副本。

使用者可以通过多个输入装置向计算机110输入命令和信息，这些输入装置例如键盘162、麦克风163（其还可以表示通过电话提供的输入）、以及指点装置161，例如鼠标、轨迹球或触摸板。其他输入装置（未示出）可以包括操纵杆、游戏键盘、圆盘式卫星天线、扫描仪，等等。这些和其他输入装置经常通过耦合至该系统总线的使用者输入接口160连接至处理单元120，但是可以通过其他接口和总线结构例如并行端口、游戏杆端口或通用串行总线（USB）连接。监视器191或其他类型的显示装置也经由一个接口例如视频接口190连接至系统总线121。除该监视器之外，多个计算机还可以包括其他外围输出装置，例如扬声器197和打印机196，这些装置可以通过一个输出外围接口195连接。

计算机110可以在一个联网环境中使用至一个或多个远程计算机例如一个远程计算机180的逻辑连接来运行。远程计算机180可以是个人计算机、手持式装置、服务器、路由器、网络PC、对等装置或其他公用网络节点，并且典型地包括上文相对于计算机110描述的多个元素中的许多或全部。图1中所描绘的多个逻辑连接包括一个局域网（LAN）171以及一个广域网（WAN）173，但是还可以包括其他网络。此类联网环境在办公室、企业范围的计算机网络、企业内部网以及互联网中是司空见惯的。

当用于一个LAN联网环境中时，计算机110通过一个网络接口或适配器170连接至LAN171。当用于一个WAN联网环境中时，计算机110典型地包括一个调制解调器172或其他装置来在WAN173例如互联网之上建立通信。可以是内部的或外部的调制解调器172可以经由使用者输入接口160或其他适当的机构连接至系统总线121。在一个联网的环境中，相对于计算机110或其多个部分描绘的多个程序模块可以存储在远程存储器存储装置中。通过举例，并非限制，图1将多个远程应用程序185图示为驻存在远程计算机180上。应了解，所示的多个网络连接是示例性的并且可以使用在多个计算机之间建立一个通信链路的其他装置。

在此申请内，“像元”用以指示一个术语、一种方法、一个对象等是指该像元方法，即，是指多个对象点特定的形状在将如在美国专利5,912,993、6,895,125以及其他相关专利中所描述的一个影像对象平滑化时的使用。例如，多个所指派的形状由多个像元内核函数定义，并且像元映射P存储关于这些像元内核函数中的哪一个/哪些被指派给每个对象点的信息。

数据的拟合是一个噪声逆问题的解。这些数据被模型化成

d=H(l)+ε，    (1)

其中d是包含所测量的数据的一个阵列，I是包含多个拟合参数的一个阵列，H是从参数空间到数据空间的变换函数，并且ε是包括向数据的噪声贡献的一个阵列。出于图示的目的，在Eq.(1)中使用适合于影像重建的一个符号，即，d为像素数据并且I为重建的影像。然而，此示例不旨在是限制性的，因为此中的发明具有向可以使用该像元方法解的任何噪声逆问题的广泛的应用。

针对许多应用，变换H是线性的，并且Eq.(1)可以写成

d=HI+ε，    (2)

其中d和I被当作向量，并且H被当作一个矩阵，称为系统矩阵。以下讨论局限于Eq.(2)的线性情况，但是其可以通过线性化推广到非线性问题，即，考虑H(I)是近于线性的一个有限的影像域。

注意，事实上，d和I可以是多维的并且甚至不需要具有相同的维度（例如，在断层成像中）。但是多维的像素总是可以安排成长向量，并且人们可以接着将H当作一个矩阵。

数据模型

m=HI=d-ε    (3)

是数据的无噪声的信号部分。其还是以下数据的期望值：

m=E(d)，    (4)

因为在没有通用性损失的情况下可以为以下所有数据点将噪声的期望值同样地设置成零：

E(ε)≡0。    (5)

噪声的协方差矩阵一般地是非零的并且可以是位置（像素位置）和/或数据模型值的一个函数

V(ε)=V(m)。    (6)

正常地，多个数据点是独立的，所以V是一个对角矩阵，但是还可以存在关联的数据，在此情况下，V具有多个非零的、非对角的元素。

给定数据d，系统矩阵H以及噪声ε的统计模型，影像重建是针对影像I解Eq.(2)——或者，更一般地，Eq.(1)——的逆问题。给定以下数据模型，影像重建是通过最优化数据的一个评价函数而求解的：

L≡L(d|m)=L(d|HI)。    (7)

评价函数的选择取决于噪声统计。评价函数通常被看作是对数似然函数（LLF），其与最大似然（例如，佩特等人，2005）相关。

针对高斯噪声，LLF减少至x²

$L={\mathrm{\χ}}^2=\underset{i}{\mathrm{\Σ}}{\mathrm{\σ}}_i^{-2}{[d_i-m_i]}^2=\underset{i}{\mathrm{\Σ}}{\mathrm{\σ}}_i^{-2}{[d_i-\underset{\mathrm{\α}}{\mathrm{\Σ}}{H}_{\mathrm{i\α}}{I}_{\mathrm{\α}}]}^2$ （高斯LLF），(8)

其中i的求和超过数据像素，a的求和超过影像像素，并且σi是像素i中的噪声的高斯标准差。高斯LLF具有在影像I中是二次的优点，所以其相对于I的梯度在I中是线性的。

针对泊松噪声，LLF是一个更复杂的非线性函数

$L=2\underset{i}{\mathrm{\Σ}}[m_i-d_i\ln \left(m_i\right)]=2\underset{i}{\mathrm{\Σ}}[\underset{\mathrm{\α}}{\mathrm{\Σ}}{H}_{\mathrm{i\α}}{I}_{\mathrm{\α}}-d_i\ln \left(\underset{\mathrm{\α}}{\mathrm{\Σ}}{H}_{\mathrm{i\α}}{I}_{\mathrm{\α}}\right)]$ （泊松LLF），(9)

其可以通过非线性最优化（例如，普雷斯等人，2007）或通过使用一个替代函数（兰格、亨特和杨，1999）求解。代替地，迈厄尔（Mighell）（1999）提议使用泊松LLF的一个二次近似，其即使在以下低计数的极限内也产生一个公正的数据模型：

$\begin{array}{c}{\mathrm{\χ}}_{\mathrm{\γ}}^2=\underset{i}{\mathrm{\Σ}}{[d_i+1]}^{-1}{[d_i+\mathrm{Min}(d_i,1)-m_i]}^2\\ =\underset{i}{\mathrm{\Σ}}{[d_i+1]}^{-1}{[d_i+\mathrm{Min}(d_i,1)-\underset{\mathrm{\α}}{\mathrm{\Σ}}{H}_{\mathrm{i\α}}{I}_{\mathrm{\α}}]}^2\end{array}$ （迈厄尔x_y ²）。(10)

存在最优化评价函数的许多方法，但是很少适合于大范围的问题，因为系统矩阵H以及其转置矩阵H^T变得太大而既不能计算也不能存储——即使使用稀疏矩阵计数——并且仅可以应用为算子。这基本上留下了两种方法：期望最大化（EM）法（理查德森（Richardson），1972；露西（Lucy），1974；邓普斯特（Dempster）、莱尔德（Laird）和鲁宾（Rubin），1977；谢普（Shepp）和瓦迪（Vardi）；1982）或者共轭梯度（CG）法（福克斯（FOX）、赫斯基（Huskey）以及威尔金森（Wilkinson），1948；海斯顿斯（Hestenes）和斯蒂菲尔（Stiefel），1952；普雷斯等人，2007）。有序子集期望最小化（OSEM）是EM的一个加速变体，其仅使用每个迭代中的数据的一个子集（哈德逊（Hudson）和拉金（Larkin），1994）。相同的过程同样可以应用于CG。

CG法通过发现评价函数的连续最小值下降梯度而进行。更精确地，其使用共轭梯度方向，这些共轭梯度方向是梯度的线性组合（例如，普雷斯等人，2007）。一个有用的概念是“负梯度（negradient）”的概念，其定义为相对于影像的评价函数的梯度的负数的一半。针对高斯评价函数，Eq.(8)，负梯度是

$G_{\mathrm{\α}}=-\frac{\∂L}{{2\∂I}_{\mathrm{\α}}}=\underset{i}{\mathrm{\Σ}}{H}_{\mathrm{i\α}}{\mathrm{\σ}}_i^{-1}[d_i-\underset{\mathrm{\β}}{\mathrm{\Σ}}{H}_{\mathrm{i\β}}{I}_{\mathrm{\β}}]$ （高斯负梯度），(11)

在向量矩阵符号中，Eq.(11)可以写为

G=H^TV^-1r，    (12)

其中

r=d-m=d–HI    (13)

是残差，并且噪声的方差是

V=Diag(σ²)（高斯噪声方差），(14)

其中Diag(σ²)指示对角矩阵，其对角元素为标准差的平方σ²。

针对泊松噪声，迈厄尔（1999）

评价函数的负梯度是

$G_{\mathrm{\α}}=-\frac{\∂L}{{2\∂I}_{\mathrm{\α}}}=\underset{i}{\mathrm{\Σ}}{H}_{\mathrm{i\α}}{[d_i+1]}^{-1}[d_i+\mathrm{Min}(d_i,1)-\underset{\mathrm{\β}}{\mathrm{\Σ}}{H}_{\mathrm{i\β}}{I}_{\mathrm{\β}}]$ （迈厄尔负梯度）。(15)

迈厄尔泊松负梯度在影像I中也是线性的，并且可以写成Eq.(12)的紧凑的向量-矩阵符号，其中残差为

r=d+Min(d,l)-m=d+Min(d,l)-HI（迈厄尔残差），(16)

并且噪声方差为

V=Diag(d+l)（迈厄尔噪声方差），(17)

CG法可以通过添加一个预处理算子（例如，戈卢布（Golub）和范龙（Van Loan），1996；萨阿德（Saad），2003）而极大地加速，从而以以下方程替代评价函数的梯度

G=P^-1H^TV^-1r（具有预处理算子的负梯度），(18)

其中线性方程

Px=b    (19)

可以容易地求解。所产生的方法称为预处理共轭梯度（PCG）法。

多个约束还经常可以通过向凸集（比蒙德（Biemond）、拉根迪克（Lagendijk）和莫塞里奥（Mersereau），1990）上投影的方法处理。例如，可以在每个迭代之后通过仅仅将影像的负分量截短成零并且好像没有截短发生一样继续来强加非负性，I≥0（佩特等人，2005）。还有用的是，在每个迭代中，将在彼迭代中具有零影像值的多个像素处的任何负共轭梯度分量设置为零。这防止这些像素参与最小化，从而只看到其结果负影像值在最小化结束时截短成零。

在此中的发明中，从更新变量且并非从已更新之后的影像确定像元映射，并且根据此像元映射平滑化该更新变量。这排除了计算一个最初伪影像的需要，因为像元映射被更新并且更新变量在每个迭代中被平滑化。

针对CG法，更新变量被看作是负梯度，Eq.(12)。（以类似方式处理一个预处理的负梯度，Eq.(18)。）负梯度通过一个像元内核库中的多个内核平滑化，并且在每个影像位置处选择的内核是最宽的内核，其中负梯度中因为平滑化的改变的平方与负梯度的方差之间的比率小于该内核以及所有较窄的内核的一个预定的阈值。

确切地说，在以下对象空间之上选择及定义内核函数库

K^(j)(τ)j=1,...,J.   (20)

内核函数的第零矩和第一矩设置为

${\mathrm{\μ}}_0^{\left(j\right)}=\∫K^{\left(j\right)}\left(\mathrm{\τ}\right)\mathrm{d\τ}\≡1,---\left(21\right)$

${\mathrm{\μ}}_1^{\left(j\right)}=\∫{\mathrm{\τK}}^{\left(j\right)}\left(\mathrm{\τ}\right)\mathrm{d\τ}\≡0,---\left(22\right)$

并且第二矩以升序排列

${\mathrm{\&mu;}}_2^{\left(j\right)}=\&Integral;{\mathrm{\&tau;}}^2{K}^{\left(j\right)}\left(\mathrm{\&tau;}\right)\mathrm{d\&tau;},{\mathrm{\&mu;}}_2^{\left(1\right)}<{\mathrm{\&mu;}}_2^{\left(2\right)}<...{\mathrm{\&mu;}}_2^{\left(J\right)},---\left(23\right)$

其中积分实际上通过像素的求和近似。第二矩通常被设计来形成一个上升的几何级数以用于多尺度影像结构的最优特征化。注意，维度在Eq.(21)至Eq.(23)中是重要的；因此，先前使用的1D影像像素索引α由向量索引τ和x替代。

接着通过使负梯度依次与内核中的每一个卷积并且减去以下原始负梯度来获得一系列平滑化的负梯度改变ΔG^(j)

$\mathrm{\&Delta;}{G}^{\left(j\right)}\left(x\right)=\&Integral;K^{\left(j\right)}\left(\mathrm{\&tau;}\right)G(x+\mathrm{\&tau;})\mathrm{d\&tau;}-G\left(x\right).---\left(24\right)$

接着，像元映射是索引的影像，如此使得在每个像素处，所选的索引是最大的索引，其中平滑化负梯度改变的平方小于或等于负梯度的p²倍

M(x)=Max_j{ΔG^(j)(x)²≤p²V(G(x))}，(₂₅)

其中p是下文讨论的一个像元因子。通过在每个影像像素处向负梯度添加对应于以下其像元索引的平滑化负梯度改变来获得平滑负梯度

$\tilde{G}\left(x\right)=G\left(x\right)+\mathrm{\&Delta;}{G}^{\left(M\left(r\right)\right)}\left(x\right).---\left(26\right)$

从Eq.(25)将非常明显的是，使用者可调整的像元因子p控制像元平滑化的强度。较大的p可以允许较大的j来在任何像素x处满足Eq.(25)，由此增加M(x)值。结果是负梯度平滑化得更强烈，Eq.(26)。相比之下，较小的p可以将j限制到一个较小的值，由此减小M(x)并且较少地平滑化。结局是使用者可以调整p来设置噪声与平滑化之间的折衷。较大的p更好地抑制噪声，代价是分辨率有些损失，而较小的p更好地保护分辨率，代价是噪声水平较高。

图4图示了确定一个影像重建过程中的像元映射的改进方法的一个示例性流程图。在步骤400中，给定Eq.(7)的数据模型，影像重建以最优化数据d的一个评价函数开始。在示例性实施例中，使用一种共轭梯度（CG）法。在步骤402中，CG法用以确定负梯度，其将用作更新变量（UV）。在步骤404中，使负梯度与来自像元内核库406的多个内核中的每一个卷积来将该UV平滑化。在每个影像位置处选择的内核是最宽的内核，其中更新变量中的因为平滑化的改变的平方与更新变量的方差之间的比率小于该内核以及所有较窄的内核的一个预定的阈值。使用多个所选的内核的多个索引的影像形成像元映射410。相同的内核通常还用以在继续进行到下一个迭代之前在彼相同位置处将更新的影像408平滑化。按照现有的像元方法，迭代地改进的影像是输出影像I的基础，该输出影像可以显示在一个显示器上和/存储在一个存储器介质中。

如先前所提及的，影像重建的示例应用不旨在是限制性的，并且改进的方法适用于使用像元方法的任何信号重建过程。将内核函数推广成非局部函数是可能的，这些非局部函数执行更复杂的平滑化。例如，在包含许多相同窗户的一个公寓建筑的一个影像中，使用一个单一的、非邻接内核来将该影像中的所有窗户一起平滑化是可能的。或者，数据可以不完全不具有几何意义，例如，股票的价格，在此情况下，以股票的某个索引方案替代像素化。针对此类系统用以下一个更一般的关系替代Eq.(24)的卷积积分

$G^{\left(j\right)}\left(x\right)=\&Integral;K^{\left(j\right)}(x,y)G(x,y)\mathrm{dy}---\left(27\right)$

维持类似Eq.(21)的以下一个归一化条件通常是有益的：

${\mathrm{\&mu;}}_0^{\left(j\right)}=\&Integral;K^{\left(j\right)}(x,y)\mathrm{dy}\&equiv;1---\left(28\right)$

但是Eq.(22)至Eq.(23)不再有意义，并且必须由问题特定的条件替代。至于涉及像元方法，要点是内核函数维持优先次序，这样使得Eq.(25)至Eq.(26)仍然有意义。

一个进一步的细化是用“分式指数”和“中间内核”通过在试验负梯度之间插入来形成一个连续的像元映射。这将负梯度进一步平滑化和/或允许使用较少的内核函数。中间像元索引定义为M(x)+δ(x)，其中M(x)是由Eq.(25)指定的整数像素索引，并且δ(x)是将ΔG(x)²/V(G(x))线性地插入p²的增量

δ(x)=[p²V(G(x))-ΔG^(M)(x)²]/[ΔG^(M+1)(v,)²-ΔG^(M)(x)²]0≤δ(x)<1,（29）

并且Eq.(26)由以下方程替代

$\tilde{G}\left(x\right)=G\left(x\right)+[1-\mathrm{\&delta;}\left(x\right)]{\mathrm{\&Delta;G}}^{\left(M\left(x\right)\right)}+\mathrm{\&delta;}\left(x\right){\mathrm{\&Delta;G}}^{(M\left(x\right)+1)},---\left(30\right)$

其对应于以下插入的内核的使用

$\tilde{K}\left(x\right)=[1-\mathrm{\&delta;}\left(x\right)]K^{\left(M\left(x\right)\right)}+\mathrm{\&delta;}\left(x\right)K^{(M\left(x\right)+1)}.---\left(31\right)$

迭代CG法照常继续进行（例如，普雷斯等人，2007），其中平滑化的负梯度替代原始的、未平滑化的负梯度。

使用Eq.(31)的插入的内核将连续的像元映射M(x)+δ(x)与其自身平滑化也是有用的。此背后的原理是像元映射应该在影像中的每个位置处以适合于在彼位置处的像元索引的像元内核的宽度的规模平滑化。因此，将该像元映射与其自身平滑化应该不将显著的额外平滑化引入到该影像中，而是可以将伪像平滑掉，这些伪像可能在计算像元映射的过程中已经被引入。

更新的影像通常在迭代之后由Eq.(31)的多个插入的内核进一步平滑化。然而，额外平滑化可能导致解随迭代的某些变动，因为将负梯度正确地平滑化的多个内核倾向于将整个影像过度平滑化。缓和该变动的一种简单的方式是用在当前以及先前的迭代中确定的内核的平均值来将该影像平滑化。（针对第一迭代，“先前”内核是一个差量函数。）

图5图示了确定使用多个插入的内核的一个影像重建过程中的像元映射的改进方法的一个示例性流程图。在步骤500中，给定Eq.(7)的数据模型，影像重建以最优化数据d的一个评价函数开始。在示例性实施例中，使用一种共轭梯度（CG）法。在步骤502中，CG法用以确定负梯度，其将用作更新变量（UV）。在步骤504中，使负梯度与来自像元内核库506的多个内核中的每一个卷积来将该UV平滑化。在每个影像位置处选择的内核是最宽的内核，其中更新变量中的改变的平方与更新变量的方差之间的比率小于该内核以及所有较窄的内核的一个预定的阈值。使用多个所选的内核的多个索引的影像形成像元映射510。在与来自该像元内核库的多个内核的卷积已经完成之后，可以在步骤508中产生多个插入的内核来将用以产生该像元映射的多个UV进一步平滑化。相同的插入的内核可以用以在继续进行到下一个迭代之前在彼位置处将更新的影像512平滑化。按照现有的像元方法，迭代地改进的影像是输出影像I的基础，该输出影像可以显示在一个显示器上和/存储在一个存储器介质中。

必须仔细地计算对确定该像元映射至关重要的负梯度的方差，因为即使在数据是统计学上独立的情况下，Eq.(12)的系统矩阵转置矩阵H^T也在多个影像像素当中产生相关性。在一些情况下，例如，当H^T是一个傅立叶变换时，可以解析地计算负梯度的方差。在其他情况下，一种计算大范围问题的方差的方便的蒙特卡洛（Monte Carlo）方法是将数据分解成许多不重叠的数据子集并且针对每个子集从Eq.(12)计算负梯度G，从而用噪声ε的一个随机实现来替代残差r。

G^(s)=H^TV^-1ε^(s)（子集梯度），(32)

其中s上标表明子集。

由于在每个像素中噪声的期望值是零——Eq.(5)，因此每个子集梯度的期望值也是零

E(G^(s))=0.    (33)

在大量子集的极限内，负梯度的方差可以因此近似为多个子集梯度的平方的总和

$V\left(G\right)\&ap;\underset{k}{\mathrm{\&Sigma;}}{\left(G^{\left(s\right)}\right)}^2.\left(34\right)$

方差的蒙特卡洛估计的相对误差是

Δ(V(G((/V(G)～S^-1/2，(35)

其中S是子集的数目。

注意，针对高斯噪声，负梯度的方差仅取决于噪声的标准差σ，并且不取决于数据或影像。因此，其不随迭代改变，并且仅需要在影像重建开始时计算一次。

该负梯度的方差并不严格地适用于迈厄尔

的泊松负梯度，因为Eq.(15)中的[d_i+1]^-1项以一种复杂的、非线性的方式取决于噪声数据。然而，实践已经表明，这些项中的方差没有Eq.(16)的残差的方差重要，残差的该方差可以从以下数据近似

V(r)≈d（泊松残差的方差）。(36)

因此，当Eq.(36)中的残差的近似方差不取决于数据时，该方差也不取决于影像。因此，该方差也不随迭代改变，并且也可能在影像重建开始时计算一次。

一种类似的像元方法还可以应用于EM重建。然而，其中的问题是乘法更新因子的方差（理查德森，1972；露西，1974；谢普和瓦迪，1982）随迭代改变并且需要在每个迭代中重新计算。除EM法的熟知的较慢分辨率恢复之外，这显著地增加计算工作量，该EM法针对相同程度的分辨率恢复需要显著地更多的迭代。综合起来，这些考虑提出理由赞成使用CG法而非EM法。

以下是使用像元方法将发明方法应用于一个输入信号的重建的多个示例。这些示例仅仅是说明性的并且不旨在是限制性的。

示例

示例1：孔径合成

孔径合成是一种类型的干涉，其混合来自一个望远镜阵列的多个信号来产生具有与具有整个阵列的大小的大小的仪器相同的角分辨率的多个影像。在每个分离与取向处，干涉仪的波瓣图产生一个输出，该输出是所观察的对象的亮度的空间分布的傅立叶变换的一个分量。从这些测量值产生源的影像（或“映射”）。孔径合成只有在输入信号的振幅和相位两者被每个望远镜测量的情况下是可能的。对于射频，这通过电子器件是可能的，而对于可见光，电磁场不能在软件中直接地测量和关联，而是必须通过灵敏的光学器件传播并且光学地干涉。

为了产生一个高质量影像，需要不同望远镜之间的许多不同的分离（如从无线电发射源所见的任意两个望远镜之间的投影的分离叫做一个基线）——需要尽可能多的不同的基线以便获得一个高质量影像。n个望远镜构成的一个阵列的基线数目（nb）通过nb=(n²-n)/2给出。例如，在天文无线电望远镜中，特大天线阵（Very Large Array；VLA）具有27个望远镜，从而立刻给出351个独立的基线，而现在正在构建中的阿塔卡玛大型毫米波/亚毫米波天线阵（Atacama Large Millimeter/sub-millimeter Array；ALMA）在完成时，将具有66个望远镜，从而给出2145个独立的基线。大多数孔径合成干涉仪使用地球的自转来增加包括在一次观察中的基线取向的数目。在不同的时间取得数据提供具有不同的望远镜分离和角度的测量值，无需购买额外的望远镜或者手动地移动这些望远镜，因为地球的自转将这些望远镜移动到新的基线上。通过允许多个个别的望远镜移动到不同的配置中提供进一步的灵活性，这给予该阵列一个强大的、可变的“缩放”。

孔径合成的其他应用包括干涉合成孔径雷达（IfSAR或InSAR）、合成孔径雷达（SAR）和逆合成孔径雷达（ISAR）、合成孔径声纳、波束成形、以及合成孔径磁场定位法。

起初认为将以本质上每个基线长度和取向的测量值理解为某个最大值是必要的：这种完全取样的傅立叶变换形式上包含完全等于来自一个常规的望远镜的影像的信息，该常规的望远镜具有一个等于最大基线的孔径直径，因此具有名称“孔径合成”。迅速地发现的是，在许多情况下，可以在多个非线性影像重建方法的帮助下用一组相对稀疏并且不规则的基线产生有用的影像。

像元方法是一种如此强大的非线性影像重建方法，但是其使用受傅立叶变换的非局部本性阻碍（巴特纳格尔和康威尔，2004）。在此中的发明中的确定像元映射的新方法可以克服此障碍。

示例2：磁共振成像

磁共振成像（MRI）是一种医学成像技术，用于放射学中来显现详细的内部结构。MRI使用核磁共振的性质来将身体内部的原子的核成像。一种MRI机器使用一个强大的磁场来校准身体中的某些原子核的磁化强度，以及多个射频场来系统地更改此磁化强度的校准。这引起这些核产生一个可由扫描仪检测的旋转磁场，并且记录此信息来构建身体的扫描区域的一个影像。磁场梯度引起不同位置处的核以不同的速度旋转。通过使用不同方向上的梯度，可以在任何任意的取向上获得多个2D影像或多个3D体积。MRI用以将身体的每个部分成像，并且特别地有用于具有许多氢原子核和小密度差的组织，例如大脑、肌肉、结缔组织以及大多数肿瘤。

身体最主要地由水分子构成。每个水分子具有两个氢原子核或质子。当一个人在扫描仪的强大的磁场内部时，许多质子的平均磁矩变得与磁场的方向对齐。暂时地接通一个射频发射器，从而产生一个变化的电磁场。此电磁场恰好具有将被吸收并且转动磁场中的质子的自旋的正确频率，称为共振频率。在关闭电磁场之后，质子的自旋返回到热力学平衡并且绝大部分磁化变得与静态磁场重新对齐。在此弛豫期间，产生一个射频信号，其可以用接收器线圈测量。不同组织中的质子以不同的弛豫速率返回到其平衡状态。不同的组织变量（包括自旋密度、T1和T2弛豫时间以及流动和频谱移位）可以用来构建影像。如在功能MRI（fMRI）以及扩散MRI中，通过改变扫描仪上的多个设置，此作用用来产生不同类型的身体组织之间的或其他性质之间的反差。

关于3D空间中的信号的起源的信息可以通过在扫描期间施加额外的磁场而得知。通过使电流穿过梯度线圈而产生的这些场使得磁场强度取决于在磁体内的位置而变化。因为这使得释放的无线电信号的频率还以一种可预测的方式取决于其起源，所以可以从该信号数学上地恢复质子在身体中的分布，典型地通过使用逆傅立叶变换。

MRI的一个主要限制是扫描的持续时间，这导致患者不适增加并且患者周转量降低。MRI扫描可以用两种方式中的一种或者其组合加速。每个射频的扫描持续时间可以缩短，和/或所扫描的射频的数目可以减少。前者导致噪声增加，而后者导致频率覆盖不完全，从而需要类似于孔径合成的技术。非线性像元影像重建可以帮助减少噪声并且补偿缺失的傅立叶数据，如下文所示。

不幸地，对于分析活体所获得的MRI数据，不存在“地面实况”或者黄金标准。因此，在麦吉尔大学构建了一个模拟大脑数据库（SBD）（在万维网上以bic.mni.mcgill.ca/brainweb可公开地获得的数据库；关、埃文斯（Evans）和派克（Pike），1999）。该SBD包含一组由一个MRI模拟器产生的实际的MRI数据量。这些数据可以由神经成像团体用来在真相是已知的背景下估计不同影像分析方法的性能。

当前，该SBD包含基于以下两个解剖模型的模拟的大脑MRI数据：正常的和多发性硬化症（MS）。针对这两者，已经使用三个序列（T1加权的、T2加权的和质子密度（PD）加权的）以及多种切片厚度、噪声水平以及强度不均匀性水平来模拟完整的三维数据量。这些数据可用于在三个正交视图（横向的、矢状的、以及冠状的）中查看，并且用于下载。关于该SBD的创建的进一步细节可通过公共资源获得。

图6示出了一个无噪声的、归一化的、质子密度（PD）SBD“真”的脑影像，并且图7示出了其向数据空间的变换的绝对值的对数（基数10）。在两个图中，上部画面示出了一个冠状的切片（左画面）以及一个矢状的切片（右画面）。下部画面是一个横向的切片。从对象空间向数据空间的变换由多个2D傅立叶变换组成，一次一个横向切片，并且显示器循环地横向地移位，以便显示零频率位于中心处。图8示出了具有随机（白）高斯噪声的图6的影像，该随机（白）高斯噪声添加到该影像，具有最亮的影像体素的3%的标准差。图9示出了图8的向数据空间的傅立叶变换。噪声被认为支配高频率的数据（数据显示的横向外围）。图10示出了与图9相同的数据，其中高频率被截短，对应于限制到可用频率的三分之一的频率覆盖范围。图10中的数据包括向影像的多个重建的输入。

图11示出了数据的一个简单的逆傅立叶变换，其示出了两个类型的伪像。首先，噪声在影像上到处清晰地可见；因为数据的频率截止，其不再是未校正的白噪声。其次，因为首先由吉布斯（Gibbs）（1898，1899）描述的现象，紧跟着多个锐边可见振影。数据（未图示）的一个非负最小二乘（NNLS）拟合本质上与直接逆傅立叶变换相同。这些影像之间的区别仅仅在体外背景下显出，该体外背景无足轻重并且总之，非常小。更显著的影像改进在图12至图14中可见，图12至图14示出了像元重建的结果，其中像元因子分别为0.3、0.5、以及1.0。所有重建显著地减少噪声以及吉布斯振影，但是并不减少相同的量。如上文在Eq.(25)至Eq.(26)的上下文中所解释，噪声与平滑化之间存在折衷，使用者可以用像元因子p控制该折衷。较大的p更好地抑制噪声，代价是分辨率有些损失，而较小的p更好地保护分辨率，代价是噪声水平较高。此折衷在图12至图14中是明显的，图12至图14随着p增加渐进地示出了较强的平滑化，其中产生分辨率的一些损失。

示例3：计算机断层成像

借助于可以检查一个患者或试验对象的哪些内部结构，而无需在该过程中对该患者进行外科手术或者伤害该试验对象，计算机断层成像（CT）为药物以及试验技术提供一种诊断和测量方法。在此情况下，存在从不同角度记录的待检查的对象的许多投影，从这些投影计算该对象的一个3D描述是可能的。

通过将多个观察到的投影（数据）转换成一个影像而产生多个断层影像。例如，在X射线CT成像中，X射线束指向一个对象并且这些射线束因为该对象内的变化的结构而衰减了不同量。在该对象的另一侧上，用多个检测器测量这些衰减的射线束。此类投影在围绕该对象的许多不同角度上产生。这些测量不仅是有噪声的，而且相对噪声水平取决于该衰减量。通过密致材料例如骨头以及尤其金属的投影比通过肉体、水或其他较不密致材料的投影具有较低的信噪比。应对许多检测到的光子中的大的并且空间上变化的波动经常需要一种统计平滑化技术来改进该影像。

在多个统计影像重建方法中，问题提出为找到最好地拟合根据该（可能地非线性的）物理模型和一个统计模型的多个测量值的多个影像。恰当的统计建模可以导致较低的噪声影像，由此允许到达患者的X射线剂量减少。

示例4：发射断层成像

发射断层成像是一种核医学成像技术，其使用类似于用于CT中的那些断层技术产生身体中的多个功能过程的一个3D影像。区别是一个发射伽玛射线的或发射正电子的放射性同位素（叫做放射性核素）被注入患者的血流中。发射伽玛射线的放射性核素发射一个单一的光子，并且该成像方法称为单光子发射计算机断层成像（SPECT，或有时SPET）。相比之下，发射的正电子消灭身体中的电子来形成两个在相反方向上移动的光子，这两个光子是同时检测到的；该成像方法称为正电子发射断层成像（PET）。

在大多数时候，仅对其放射性性质重要的标记放射性同位素附加到一个特定的配体来产生一个放射性配体，该放射性配体对其与某些类型的组织的化学粘合性质是重要的。此密切结合允许配体和放射性同位素的组合（放射性药物）被运送到身体中的一个重要的地方并且束缚于此，其接着（因为该同位素的直接或间接的伽玛射线发射）允许配体浓度被成像。

核扫描连同CT或MRI扫描正越来越多地被了解，其中组合给出了解剖信息和新陈代谢信息两者（即，结构是什么，并且其正在生化地做什么）。因为核成像与解剖成像结合是最有用的，所以现代的核扫描仪现在是可与以下仪器一起使用的：集成的高端多排检测器CT扫描仪或更近地，MRI。因为这两种扫描可以按即时序列甚至同时地执行，所以在同一会话期间，伴随着病人不改变在两种类型的扫描之间的位置，两组影像被更精确地标示，这样使得核影像上的异常区域可以与CT或MRI影像上的解剖更完美地关联。这在示出具有较高解剖变异的移动的器官或结构的详细视图方面是非常有用的，该较高解剖变异在大脑外部更常见。

一种与计算机断层成像（CT）的重建非常像的技术通常用以产生3D影像，但是在核成像过程中收集的数据集具有比CT少得多的光子，因此重建技术是更困难的。像元处理可以接着用以提高影像质量和/或减少注入患者体内的剂量。

示例5：频谱分析

一种频谱分析器测量一个输入信号对在该仪器的完整频率范围内的频率的大小。主要使用是测量已知信号和未知信号的频谱的功率。一个频谱分析器测量的输入信号是电的，但是其他信号例如声压波和可见光波的频谱组分可以通过一个适当的换能器的使用而考虑。通过分析这些电信号的频谱，可以观察一个信号的主频率、功率、失真、谐波、带宽和其他频谱分量，这些在时域波形中并不是容易地可检测的。这些参数在多个电子装置例如无线发射器的特征化过程中是有用的。

频谱分析器类型是通过用以获得一个信号的频谱的方法而指定的。存在扫描调谐式以及基于快速傅立叶变换（FFT）的频谱分析器：一种扫描调谐式频谱分析器使用一个超外差接收器来将输入信号频谱的一部分下转换成一个带通滤波器的中心频率（使用一个受电压控制的振荡器和一个混频器）。使用一个超外差架构，该受电压控制的振荡器扫过一系列频率，从而允许考虑该仪器的完整频率范围。一种FFT频谱分析器计算该快速傅立叶变换，由此将输入信号变换成其频率谱的多个分量。一些频谱分析器（例如实时频谱分析器）使用一种混合的技术，其中首先使用多种超外差技术将输入信号下转换成一个较低的频率并且接着使用多种FFT技术分析该输入信号。

这是一个1D示例，其中频谱（“影像”）通过一个非局部的傅立叶变换而与输入数据相关。该像元方法可以用与所引用的其他傅立叶示例类似方式应用于此类数据。

参考文献

1.巴特纳格尔·S、康威尔·TJ，2004，天文与天体物理学报，第426卷：747-754页

2.比蒙德·J、拉根迪克·RL和莫塞里奥·RM，1990，电子和电气工程师协会会报，第78卷：856-883页

3.坎德丝·EJ、龙伯格·J和陶·T，2004，电子和电气工程师协会导航信息论，第52卷：489-509页

4.蔡廷·GJ，1966，J.Assoc.Comput.Mach，第13卷：547-569页

5.陈·S、多纳霍·D和桑德斯·MA，1999，SIAM J.Sci.Comp.，第20卷：33-61页

6.邓普斯特·AP、莱尔德·NM和鲁宾·DB，1977，皇家统计学协会期刊B系列，第39卷：1-38页

7.多纳霍·DL，2006，商业纯应用数学，第59卷：797-829页

8.多纳霍·DL、曾（Tsaig）·Y、德罗里·I、斯达克·JL，2006，“阶梯式正交匹配追踪未确定的线性方程的稀疏解（Sparse Solution ofUnderdetermined Linear Equations by Stagewise Orthogonal MatchingPursuit）”，斯坦福技术报告

9.费希尔·RA，1912，信使数学（Messenger Math），第41卷：155-160页

10.费希尔·RA，1922，Phil.Trans.Roy.Soc.A，第222卷：309-368页

11.福克斯·L、赫斯基·HD、威尔金森·JH，1948，Quart.J.Mech.Appl.Math，第1卷：149-173页

12.高斯·CF，1809，Theoria Motus Corporum Coelestium in sectionibusconicis solem ambientium，英译本，1963，天体运动论（纽约：多佛）

13.吉布斯·JW，1898，自然，第59卷：200页

14.吉布斯·JW，1899，自然，第59卷：606页

15.戈卢布·GH、范龙·CF，1996，矩阵计算，第三版（巴尔的摩：约翰〃霍普金斯大学出版社）

16.海赛·LO，1876，Vorlesungen iiber Analytische Geometrie desRaumes，Insbesondere iiber Oberflachen Zweiter Ordnung，（莱比锡：托伊布纳）

17.海斯顿斯·MR、斯蒂菲尔·E，1952，J.Res.Nat.Bur.Standards，第49卷：409-436页

18.哈德逊·HM、拉金·RS，1994，IEEE Trans.Med.Imag.，第13卷：601-609页

19.柯尔莫戈洛夫·AN，1965，Problems Inf.Transm.，第1卷：4-7页

20.关·RKS、埃文斯·AC和派克·GB，1999，IEEE Trans.Med.Imag.，第18卷：1085-1097页

21.兰格·K、亨特·DR和杨·I，2000，J.Comp.Graph.Statist.，第9卷：1-59页

22.露西·LB，1974，Astron.J.，第79卷：745-754页

23.迈厄尔·KJ，1999，Astrophys.J.，第518卷：380-393页

24.米勒·A，2002，回归子集选择，第二版（博卡拉顿：查普曼和霍尔/CRC）

25.奈奎斯特·H，1928，Trans.AIEE，第47卷：617-644页

26.碧娜·RK、佩特·RC，1993，Publ.Astron.Soc.Pac.，第105卷，630-637页

27.普雷斯·WH、图科斯基（Teukolsky）·SA、维特琳（Vetterling）·WT和弗兰诺雷（Flannery）·BP，2007，数值分析方法库，第三版，（剑桥：剑桥大学出版社）

28.佩特·RC、格斯内尔·TR和亚希勒·A，2005，Annu.Rev.Astron.Astrophys.，第43卷，139-194页

29.佩特·RC、亚希勒·A，1999，Astron.Soc.Pac.Conf.Ser.，第172卷：307-316页

30.理查德森·W，1972，J.Opt.Soc.Am.，第62卷：55-59页

31.萨阿德·Y，2003，稀疏线性系统的迭代方法，第二版，（费城：SIAM）

32.香农·CE，1949，Proc.Institu.Radio Engineers，第37卷：10-21页

33.谢普·LA、瓦迪·Y，1982，IEEE Trans.Med.Imaging，第1卷：113-122页

34.索洛莫诺夫·RJ，1964，Inf.Control，第7卷：1-22页

35.蒂施莱尼·R，1996，J.R.Statist.Soc.B，第58卷：267-288页

36.吉洪诺夫·AN，1963，苏联数学，第4卷：1035-1038页

Claims (21)

1.一种从自物理过程获得的数据集重建对象模型的方法，其中，该数据集包含噪声，该方法包括：

接收在一个数据空间中定义的该数据集；

在一个对象空间中构建一个对象模型，其中，该对象模型包括多个对象点；

设计该对象模型从对象空间向数据空间的一个变换，从而产生一个数据模型，其中，该变换对应于获得该数据集的物理过程；

选择用于确定一个数据模型向该数据集的一个拟合的一个评价函数；

基于该评价函数确定对象空间中的该对象模型的一个更新变量；

将该更新变量平滑化以通过以下方法确定一个平滑化的更新变量：

使该更新变量与多个像元内核中的每一个卷积；以及

为该输入对象的每个对象点选择一个像元内核，该像元内核具有符合一个预定的最小标准的一个最大大小；

通过在对应于所选的像元内核的每个对象点处指派多个索引来产生一个像元映射；以及

基于该像元映射内的这些索引产生一个包括一个实质上降噪的对象模型的输出。

2.如权利要求1所述的方法，进一步包括根据该像元映射通过使每个对象点处的对象模型与该所选的像元内核卷积来将该对象模型平滑化。

3.如权利要求1所述的方法，其中，使用一种共轭梯度法来确定该评价函数。

4.如权利要求3所述的方法，其中，该更新变量是一个负梯度。

5.如权利要求3所述的方法，进一步包括添加一个预处理算子。

6.如权利要求1所述的方法，进一步包括：

迭代地将该更新变量平滑化并且基于该平滑化的更新变量更新该像元映射以及该对象模型，直到将该评价函数最优化为止。

7.如权利要求1所述的方法，进一步包括在对该更新变量进行卷积的步骤之后，减去该更新变量来确定一个更新变量改变。

8.如权利要求7所述的方法，其中，该预定的最小标准是基于该更新变量改变的平方与该更新变量的一个方差之间的比率的。

9.如权利要求1所述的方法，其中，该数据集包括结合一个分析过程产生的干涉数据，该分析过程选自由以下各项组成的群组：孔径合成、干涉合成孔径雷达（IfSAR或InSAR）、合成孔径雷达（SAR）和逆合成孔径雷达（ISAR）、合成孔径声纳、波束成形、以及合成孔径磁场定位法。

10.如权利要求1所述的方法，其中，该数据集包括磁共振成像数据。

11.如权利要求1所述的方法，其中，该数据集包括计算机或发射断层成像数据。

12.如权利要求1所述的方法，其中，该数据集包括频谱分析数据。

13.如权利要求1所述的方法，其中，该数据集包括来自不同仪器的多个输入信号，并且该降噪的对象模型包括对应于这些输入信号的一个混合物的一个信号输出。

14.一种非瞬态计算机可读介质，其上具有附带的软件，该软件包括用于从自物理过程获得的数据集重建对象模型的多个指令，其中，该数据集包含噪声，该方法包括：

接收在一个数据空间中定义的该数据集；

在一个对象空间中构建一个对象模型，其中，该对象模型包括多个对象点；

设计该对象模型从对象空间向数据空间的一个变换，从而产生一个数据模型，其中，该变换对应于获得该数据集的物理过程；

选择用于确定一个数据模型向该数据集的一个拟合的一个评价函数；

基于该评价函数确定对象空间中的该对象模型的一个更新变量；

将该更新变量平滑化以通过以下方法确定一个平滑化的更新变量：

使该更新变量与多个像元内核中的每一个卷积；以及

为该输入对象的每个对象点选择一个像元内核，该像元内核具有符合一个预定的最小标准的一个最大大小；以及

通过在对应于所选的像元内核的该输入对象的每个对象点处指派多个索引来产生一个像元映射。

15.如权利要求14所述的非瞬态计算机可读介质，进一步包括根据该像元映射通过使每个对象点处的对象模型与该所选的像元内核卷积来将该对象模型平滑化。

16.如权利要求14所述的非瞬态计算机可读介质，其中，使用一种共轭梯度法来将该评价函数最优化。

17.如权利要求16所述的非瞬态计算机可读介质，其中，该更新变量是一个负梯度。

18.如权利要求17所述的非瞬态计算机可读介质，进一步包括添加一个预处理算子。

19.如权利要求14所述的非瞬态计算机可读介质，进一步包括：

迭代地将该更新变量平滑化并且基于该迭代地平滑化的更新变量更新该像元映射以及该对象模型，直到将该评价函数最优化为止。

20.如权利要求19所述的非瞬态计算机可读介质，进一步包括在对该更新变量进行卷积的步骤之后，减去该更新变量来确定一个更新变量改变。

21.如权利要求20所述的非瞬态计算机可读介质，其中，该预定的最小标准是基于该更新变量改变的平方与该更新变量的一个方差之间的比率的。

Images