CN100550056C - 在图像序列中进行运动校正的方法及系统 - Google Patents

Abstract

一种用于在时间序列中的第一与第二图像之间进行运动补偿的方法包括：用减小过程处理第一与第二图像，用于提供相应的分辨率减小的第一与第二图像；由相应的分辨率减小的第一与第二图像导出相应的第一与第二特征映射，这些特征映射包括导出相应的分辨率减小的第一与第二图像中的图像数据的相应的拉普拉斯算子；通过根据配准算法处理第一与第二特征映射而导出位移场，该配准算法包括针对每个像素或体素求解局部高斯最小均方问题，以便导出形成位移场的相应矢量；以及利用位移场弯曲第二图像。

Description

在图像序列中进行运动校正的方法及系统

相关申请的交叉引用

优先权

这里特别参考以本申请的发明人CHRISTOPHE CHEFD’HOTEL的名义于2004年8月31日提交的申请号为No.60/605,759的美国临时专利申请，该美国临时专利申请的题为“用于乳腺MR运动校正的高斯加权最小均方配准算法(A GAUSSIAN WEIGHTED LEAST MEAN SQUAREREGISTRATION ALGORITHM FOR BREAST MR MOTION CORRECTION)”并且代理人案卷号为2004P14822US，而且本申请要求该美国临时专利申请的优先权，且其公开内容在此引入作为参考。

技术领域

本申请通常涉及成像中的运动校正，并且更特别地涉及可能属于校正对象运动的运动补偿，该对象运动如可能在(诸如可以通过医学成像过程得到的)图像时间序列中发生的那样。

背景技术

医学成像技术被用于许多医学程序，包括例如被用于对患者进行癌症或癌症前期状况的检测。重要的应用是对乳腺癌中的肿瘤或潜在肿瘤的检测。潜在肿瘤是难于检测的。在提供潜在地有帮助的信息的可得到的技术中，公知的是，例如，这样的涉及肿瘤的组织与邻近的、非肿瘤组织相比通常可能呈现对造影剂的更快速的吸入(灌注(wash-in))以及更快速地清除。诸如这些及其它的特性可能在某些诊断中是有帮助的，这些诊断涉及检测可疑组织及通过将对患者进行一程序(诸如灌注和/或清除造影剂)之前以及进行该程序之后的图像进行比较来识别组织特性。使用由成像技术(诸如磁共振成像(MRI))所形成的这种时序图像，可以在图像之间进行比较，以检测由所采集的MR体积的不同区域所呈现的性能差异。

用于执行该检测的技术有利地要求跟踪这样的体积的时间序列中的单个体素的强度。然而，困难在于，患者通常在连续的采集之间移动，并且由此引入所采集的图像之间的涉及运动的差异，由此不再能跟踪空间中的单个点，除非执行运动校正。如在此所用的那样，空间中的点不是指典型的用几何学所定义的没有尺寸的点，而是指由数字化程序产生的具有小尺寸元素的点，这些小尺寸元素组成数字化图像。

过去解决这个问题的现有技术方法已计算出两个图像之间的光流，在这两个图像中选择任意一个作为该序列图像中的参考。例如，所述两个图像可以通过计算拉普拉斯金字塔而从所采集的图像中获得。光流可以例如通过求解两个拉普拉斯图像之间的点对点差异的最小化问题来计算。

发明内容

本发明的目的是以有利的方式解决在例如乳腺潜在肿瘤的MR检测中的运动校正或补偿问题，这些潜在肿瘤被检测为具有快速吸入(灌注)造影剂以及快速清除造影剂的组织。

根据本发明的一方面，患者的第一与第二图像由成像过程(诸如MRI设备)依照时间顺序获得，其中第二及以后的图像可以包括因例如在取得第一图像之后患者运动而引起的与第一和更早的图像的差异。用减小过程(reduction process)处理这些图像，从而导致如与原始图像相比分辨率更低的相应的第一与第二图像。通过导出较低分辨率的第一与第二图像中的图像数据的相应的拉普拉斯算子，相应的第一与第二特征映射得自较低分辨率的第一与第二图像。基于给定的初始位移场，根据配准算法用配准程序处理这些特征映射。该算法包括针对每个体素求解局部高斯加权最小均方问题，以便导出形成用于模拟变形的密集位移场的相应矢量。利用该位移场来弯曲(warp)第二图像，以便获得运动校正过的第二图像。

给定的初始位移场的缺省条件可以例如是空集或零位移，或者该给定的初始位移场可以考虑患者运动的先验知识。

根据本发明的另一方面，用于在时间序列中的第一与第二图像之间进行运动补偿的方法包括：由该第一与第二图像导出相应的第一与第二特征映射；通过根据配准算法处理这些第一与第二特征映射而导出位移场，所述配准算法包括针对每个像素或体素求解局部高斯加权最小均方问题，以便导出形成所述位移场的相应矢量；以及利用所述位移场弯曲所述第二图像。

根据本发明的另一方面，导出位移场的步骤利用前面导出的位移场来导出所述位移场，而针对最初导出所述位移场，导出位移场的步骤利用给定的初始位移场。

给定的初始位移场的缺省条件是空集，或者例如，给定的初始位移场可以考虑患者运动的先验知识。

根据本发明的另一方面，重复导出位移场的步骤，其中在对应于比前次重复具有更高分辨率的第一与第二图像的第一与第二特征映射上执行每次重复。

根据本发明的另一方面，每次重复利用来自直接在前的步骤的前面导出的位移场。

根据本发明的另一方面，针对最初导出位移场，导出位移场的步骤利用给定的位移场。

根据本发明的另一方面，导出位移场的步骤包括将在前导出的位移场扩展到分辨率增加版本的分辨率水平。

根据本发明的另一方面，用于进行运动补偿的方法包括以下步骤：利用为了弯曲该第二图像而存在的、在最高分辨率导出的位移场来产生运动校正过的图像。

根据本发明的另一方面，导出相应的第一与第二特征映射的步骤包括导出相应的分辨率减小的第一与第二图像中的图像数据的相应拉普拉斯算子。

根据本发明的另一方面，用于在时间序列中的第一与第二图像之间进行运动补偿的方法包括：用减小过程处理第一与第二图像，用于提供相应的分辨率减小的第一与第二图像；由相应的分辨率减小的第一与第二图像导出相应的第一与第二特征映射，这些特征映射包括导出该相应的分辨率减小的第一与第二图像中的图像数据的相应的拉普拉斯算子；通过根据配准算法处理第一与第二特征映射而导出位移场，所述配准算法包括针对每个像素或体素求解局部高斯加权最小均方问题，以便导出形成该位移场的相应矢量；以及利用该位移场来弯曲该第二图像。

根据本发明的另一方面，用于在时间序列中的第一与第二图像之间进行运动补偿的方法包括：用减小过程处理第一与第二图像，用于提供相应的分辨率减小的第一与第二图像；由该相应的分辨率减小的第一与第二图像导出相应的第一与第二特征映射；通过根据配准算法处理第一与第二特征映射而导出位移场，所述配准算法包括针对每个像素或体素求解局部高斯加权最小均方问题，以便导出形成该位移场的相应矢量；以及利用该位移场弯曲该第二图像。

根据本发明的另一方面，用于在时间序列中的参考图像与浮动图像之间进行运动补偿的方法包括：由该参考图像导出分辨率逐渐减少的第一组修正图像；由该浮动图像导出分辨率逐渐减少的第二组修正图像；由该第一组修正图像导出第一组特征映射；由该第二组修正图像导出第二组特征映射；根据配准算法，由分别第一与第二组特征映射中的每个特征映射中的最低分辨率组成部分和给定的初始位移场导出第一位移场，所述配准算法包括针对每个像素或体素求解局部高斯加权最小均方问题，以便导出形成该位移场的相应矢量；根据所述配准算法，由分别第一与第二组特征映射中的每个特征映射中的次最低分辨率组成部分和在前一步骤中所获得的第一位移场导出第二位移场；针对分别第一与第二组特征映射中的每个特征映射中的分辨率连续较高的组成部分(如果存在)重复前述步骤，并且在每种情况下都使用在当前步骤之前的步骤中所获得的位移场，直至达到浮动图像的分辨率，而如果不存在分辨率较高的组成部分，则直接进行到下一步；以及利用最后获得的位移场弯曲该浮动图像。

根据本发明的另一方面，用于在时间序列中的参考图像与浮动图像之间进行运动补偿的方法包括：(a)用相应的第一与第二多个级联的分辨率减小过程处理该参考图像与浮动图像，用于提供相应多个分辨率连续减小的参考图像与浮动图像，在此对于最低分辨率水平称为水平0，而对于次较高分辨率水平称为水平1，对于第二次较高分辨率水平称为水平2，对于任何现有较高分辨率水平以此类推；(b)由分辨率水平为L0、L1、L2等的相应多个分辨率连续减小的参考图像与浮动图像导出相应多个参考特征映射与浮动特征映射；(c)通过根据配准算法处理与水平L0相对应的参考特征映射与浮动特征映射而导出第一位移场，所述配准算法包括基于给定的初始位移场针对每个像素或体素求解局部高斯加权最小均方问题，以便导出形成该第一位移场的相应矢量；(d)将该第一位移场扩展到与分辨率L1相适合的分辨率水平，以提供被扩展的第一位移场；(e)通过根据该配准算法基于被扩展的第一位移场来处理与分辨率水平L1相对应的参考特征映射与浮动特征映射，从而导出第二位移场，以便导出形成所述第二位移场的相应矢量；(f)将该第二位移场扩展到与分辨率L1相适合的分辨率水平，以提供被扩展的第二位移场；(g)通过根据该配准算法基于被扩展的第二位移场来处理与分辨率水平L2相对应的参考特征映射与浮动特征映射，从而导出第三位移场，以便导出形成所述第二位移场的相应矢量；以及(h)如果L2是参考图像与浮动图像的分辨率水平，那么通过利用第三位移场弯曲该浮动图像并结束；而如果L2不是参考图像与浮动图像的分辨率水平，则(i)将该第三位移场扩展到与次较高分辨率水平相适合的分辨率水平，以提供被扩展的第三位移场；以及(j)重复从步骤(g)开始的一系列步骤，适当修改分辨率水平，以基于前面最后获得的被扩展的位移场而获得对应于越来越高的分辨率水平的位移场，直到达到参考图像与浮动图像的分辨率水平，并因此利用所获得的最后的位移场弯曲该浮动图像以及结束。

根据本发明的另一方面，提供相应多个分辨率连续减小的参考图像与浮动图像的步骤包括导出相应多个分辨率连续减小的参考图像与浮动图像中的图像数据的相应的拉普拉斯算子。

根据本发明的另一方面，导出第一位移场的步骤包括输入给定的初始位移场。

根据本发明的另一方面，用于在时间序列中的第一与第二图像之间进行运动补偿的方法包括：由相应的分辨率减小的第一与第二图像导出相应的第一与第二特征映射；通过根据配准算法处理第一与第二特征映射而导出第一位移场，该配准算法包括针对每个像素或体素求解局部高斯加权最小均方问题，以便导出形成该位移场的相应矢量；扩展该第一位移场，以与该第一与第二图像的分辨率相符合；通过根据配准算法基于第一位移场来处理第一与第二图像而导出第二位移场；以及利用该第二位移场弯曲该第二图像，以便获得运动校正过的图像。

根据本发明的另一方面，用于在时间序列中的第一与第二图像之间进行运动补偿的方法包括：由第一图像导出具有逐渐减小的分辨率的第一系列图像；由第二图像导出具有逐渐减小的分辨率的第二系列图像；由第一系列图像导出第一系列特征映射；由第二系列图像导出第二系列特征映射；通过根据配准算法处理具有第二系列特征映射的相对应的特征映射的第一系列特征映射的特征映射(这些特征映射具有相同分辨率)而导出系列位移场，以提供相应的位移场，该配准算法包括针对每个像素或体素求解局部高斯加权最小均方问题，以便导出形成该位移场的相应矢量，以给定的位移场开始，该给定的位移场被用于导出分辨率水平最低的第一位移场，并且此后利用直接在前的位移场来导出对应于该系列中的次较高分辨率的下一个位移场；以及利用分辨率水平为第一与第二图像的分辨率水平的最后的位移场来弯曲该第二图像。

根据本发明的另一方面，用于在时间序列中的参考图像与浮动图像之间进行运动补偿的方法包括：(a)用相应的第一与第二多个级联的分辨率减小过程处理参考图像与浮动图像，用于提供相应多个分辨率连续减小的参考图像与浮动图像，在此对于最低分辨率水平称为水平0，而对于次较高分辨率水平称为水平1，对于第二次较高分辨率水平称为水平2，对于任何现有的较高分辨率水平以此类推；(b)由按照分辨率增加的顺序排列的、分辨率水平为L0、L1、L2的相应多个分辨率连续减小的参考图像与浮动图像导出相应多个参考特征映射与浮动特征映射，其中L2是该参考图像与浮动图像的分辨率水平；(c)通过根据配准算法处理与水平L0相对应的参考特征映射与浮动特征映射而导出第一位移场，该配准算法包括基于给定的初始位移场针对每个像素或体素求解局部高斯加权最小均方问题，以便导出形成该第一位移场的相应矢量；(d)将该第一位移场扩展到与分辨率L1相适合的分辨率水平，以提供被扩展的第一位移场；(e)通过根据该配准算法基于被扩展的第一位移场来处理与分辨率水平L1相对应的参考特征映射与浮动特征映射而导出第二位移场，以便导出形成该第二位移场的相应矢量；(f)将该第二位移场扩展到与分辨率L2相适合的分辨率水平，以提供被扩展的第二位移场；(g)根据该配准算法，基于被扩展的第二位移场，通过处理与分辨率水平L2相对应的参考特征映射与浮动特征映射而导出最后的位移场，以便导出形成最后的位移场的相应矢量；以及(h)通过利用该最后的位移场弯曲该浮动图像。

根据本发明的另一方面，用于执行图像运动补偿的系统包括：存储装置，用于存储程序及其它数据；以及与该存储装置通信的处理器，该处理器运行该程序来执行：用于在时间序列中的第一与第二图像之间进行运动补偿的方法，该方法包括：用减小过程处理该第一与第二图像，用于提供相应的分辨率减小的第一与第二图像；由该相应的分辨率减小的第一与第二图像导出相应的第一与第二特征映射；通过根据配准算法处理该第一与第二特征映射而导出位移场，该配准算法包括针对每个像素或体素求解局部高斯加权最小均方问题，以便导出形成该位移场的相应矢量；以及利用该位移场弯曲该第二图像。

根据本发明的另一方面，一种计算程序产品包括计算机可用介质，该计算机可用介质具有针对程序代码而被记录于其上的计算机程序逻辑，用于通过用于在时间序列中的第一与第二图像之间进行运动补偿的方法来执行图像运动补偿，该方法包括：用减小过程处理该第一与第二图像，用于提供相应的分辨率减小的第一与第二图像；由该相应的分辨率减小的第一与第二图像导出相应的第一与第二特征映射；通过根据配准算法处理该第一与第二特征映射而导出位移场，该配准算法包括针对每个像素或体素求解局部高斯加权最小均方问题，以便导出形成该位移场的相应矢量；以及利用该位移场弯曲该第二图像。

附图说明

本发明结合附图从详细说明中被更全面地理解，所述详细说明包括后面的示例性实施例，在所述附图中：

图1以流程图格式示出根据本发明原理的对图像时间序列的运动校正；以及

图2示意性地示出用于实现本发明的可编程的数字计算机的应用。

具体实施方式

根据本发明的实施例，患者的参考图像通过医学成像程序(诸如利用MR I设备)获得，而浮动图像通常在比参考图像稍晚的时刻取得。浮动图像可以包括由于例如在取得第一图像之后患者的移动而引起的与参考图像的差异。用减小过程处理这些图像，从而产生如与原始的相应图像相比较分辨率更低的相应的参考图像和浮动图像。相应的参考图像和浮动图像特征映射通过导出分辨率较低的参考图像和浮动图像中的图像数据的相应的拉普拉斯算子而得自分辨率较低的参考图像和浮动图像。基于给定的初始位移场，特征映射用配准程序根据配准算法来处理。该算法包括针对每个体素求解局部高斯加权最小均方问题，以便导出形成用于模拟变形的密集位移场的相应矢量。利用该位移场来使浮动图像弯曲，以便获得运动校正过的第二图像。如上所述，给定的初始位移场的可能的缺省条件是空集，或者例如，给定的初始位移场可以考虑患者运动的先验知识。

用于乳腺MR运动校正的高斯加权最小均方配准算法被用于在检测涉及肿瘤的组织中提供运动校正，而用于乳腺MR运动校正的非刚性配准算法被用于在检测涉及肿瘤的组织中提供运动校正。所述变形被模拟为密集位移场，其由通过求解一系列局部高斯加权最小均方问题而获得的矢量形成，如后面将描述的那样。

利用配准算法估计立体图像之间的密集位移场。目的是补偿磁共振(MR)图像时间序列中的运动伪影。通常，运动伪影主要是由于呼吸、心脏运动、及患者的运动而引起的。本发明的主要应用领域是乳腺MR研究中强度变化的研究。造影剂在乳腺组织中吸入的时间变化可被用于提供有助于检测损伤的信息。计算这种变化要求图像采集之间的乳腺组织的准确的空间对准。

在此认识到，引入使用快速滤波技术所实现的高斯加权提供至少两个优点，即：

引入高斯加权提供一种简单且计算有效的方式来确保位移场的平滑度，从而防止奇点并提供变形规则性的隐式模型(implicitmodel))；以及

引入高斯加权确保最小均方问题被适当地公式化。

另外，根据本发明原理的图像相似性量度依赖于图像拉普拉斯算子的计算，而不是使用图像强度。本发明人已发现，这种方法在处理因造影剂吸入而引起的强度变化方面是有效的。使用多分辨率策略来提高捕获范围、速度、以及该方法的稳健性。

本发明将通过图示实例进一步详细地进行解释。对本发明的图示实例来说，假设要被配准的图像对的强度是两个实值函数

和

的(规则网格上的)采样版本，其中那些符号具有其传统含义；因此，在本实例中，I₂是包括输入Ω和输出R的maplet的函数，其中Ω是R³的真子集，等等。注意的是，这里处理立体图像的情况；但随后的研发容易延伸到任意维数。

在下面，将I₁称为参考图像，而将I₂称为浮动图像。将配准问题限定为找到位移场

(这将以网格点来估计)，该位移场

将I₁与浮动图像的弯曲版本I₂ο(id+p)(其中id是标识映射，而^ο是合成算子)之间的相似性最大化。

图像强度不总是执行配准任务的最好特征。从这个角度来看，与由I₁和I₂计算出的辅助函数

和

一起工作有时是有用的，所述辅助函数代表其它图像特征(例如图像拉普拉斯算子，或者图像梯度的范数)。尤其是在对乳腺MR的应用中，在实践中计算图像拉普拉斯算子一直非常有效。

该算法是迭代的。该算法包括建立向“真正的”运动场收敛的位移场序列p⁰，p¹，...，p^k。最后的“运动校正过的”图像由I₂ο(id+p^k)给出。下面描述位移场的迭代构建，该位移场的迭代构建形成这个配准技术的核心。

迭代程序结合粗糙到精细的多分辨率策略来应用。该方法的优点是：

通过降低落入局部最小值的风险来提高捕获范围；以及

通过与更少的分辨率较低的数据一起工作来降低计算成本。

多分辨率策略可以通过如下例子进行描述：

通过取平均及子采样来针对参考图像及浮动图像建立两个多分辨率金字塔。这种类型的金字塔构建是图像处理的标准技术；详细内容例如参考由2.R.Gonzales和R.Woods编写的教科书Digital ImageProcessing(第二版，Prentice-Hall，2002年)。

令I₁ ^l和I₂ ^l是分辨率为1的所得到的图像。

计算被标记为J₁ ^l和J₂ ^l的对应特征图像。

令p^k，l为在分辨率1处第k次迭代之后获得的位移场。p^0，0表示分辨率最低时的初始位移场(通常是空白字段)。在分辨率1处，如下面将进一步解释的那样，运行配准算法固定数目的迭代。通过三线性内插法及缩放将位移场扩展到下一分辨率。将该结果用作初始位移场p^0，l+1。

执行前面的操作，直至达到原始分辨率。

在图1中以概括图表的方式图示说明了这些步骤。参照图1，分辨率被表示为“1”。通过示例性实施例以及出于图解说明的目的，在该图中示出了三个水平；更多或更少的扩展步骤被认为在本发明的范围内。在图1中，分别被示为I₁和I₂的、分辨率水平为l＝l_max＝2的参考图像及浮动图像经过第一减小过程，从而得到分辨率为l＝1的相应图像I₁ ¹和I₂ ¹，在又一减小过程之后，所得到的图像分别是分辨率水平为l＝0的I₁ ⁰和I₂ ⁰特征映射J₁ ⁰和J₂ ⁰分别得自图像I₁ ⁰和I₂ ⁰。特征映射J₁ ⁰和J₂ ⁰使用给定的初始位移场p^0，0经历如在此所描述的配准过程。该初始位移场可以是任意的且可以是零。配准过程得到位移场p^k，0，该位移场p^k，0与初始位移场相比被认为是改进的位移场。然后，位移场p^k，0经历与图像分辨率l＝1相对应的扩展或膨胀，从而导致产生被扩展的分辨率场p^0，1。

特征映射J₁ ¹和J₂ ¹分别得自图像I₁ ¹和I₂ ¹。特征映射J₁ ¹和J₂ ¹使用前面导出的被扩展的位移场p^0，1经历如在此所描述的配准过程，以便得到位移场p^k，1，该位移场p^k，1被扩展到p^0，2。对图像I₁ ²和I₂ ²执行与以前类似的步骤，图像I₁ ²和I₂ ²导致产生最后的变形场，该最后的变形场被用来执行对图像I2的最后弯曲，以产生如在这里所解释的由I₂ο(id+p)所代表的运动校正过的图像。

应注意到的是，仅仅必需将该算法的本征参数设定为如下：

金字塔水平的数目，

在每一水平迭代的数目，以及

控制位移场的平滑度的正则化参数σ：下面讨论该参数。

根据本发明的原理，在描述特征映射计算之后，接下来阐述根据本发明的高斯加权最小均方算法的更详细的描述。

首先将特别描述配准技术的核心概念。在给定的分辨率处执行下面的迭代程序。下面表示法中的分辨率水平1被省略。

令G_σ(x)为标准偏差σ的各向同性的三维高斯核：

$G_{\mathrm{\σ}}\left(x\right)=\frac{1}{\sqrt{{\left(2\mathrm{\π}\right)}^3{\mathrm{\σ}}^6}}{e}^{\frac{{\left|\right|x\left|\right|}^2}{{\mathrm{\σ}}^2}}$

每个体素由下标j来标记。其坐标由矢量x_j∈R³给出。将 $p_j^k\∈R^3$ 称为该算法第k次迭代之后在点x_j处重新获得的位移矢量。假设一组位移矢量p_j ^k，在第k次迭代时的浮动特征映射的校正版本被表示为J_2，k：

$\∀j,J_{2,k}\left(x_j\right)=J_2(x_j+p_j^k)$

在每个体素j处，提出用非线性最小均方问题的解s_j ^k来定义更新规则 $p_j^{k+1}=p_j^k+s_j^k:$

$\underset{s_j^k}{\min }\underset{i}{\mathrm{\Σ}}{G}_{\mathrm{\σ}}(x_j-x_i)\·{\left(z_i^k\left(s_j^k\right)\right)}^2$

其中， $z_i^k\left(s\right)=J_1\left(x_i\right)-J_{2,k}(x_i+s).$ 高斯加权允许控制被用来估计位移的邻域尺寸(的隐含性)。作为副作用，高斯加权也确保变形场的整个平滑度。

在实践中，前面问题的线性化版本通过最小化下面的基准来求解：

$E\left(s_j^k\right)=\underset{i}{\mathrm{\Σ}}{G}_{\mathrm{\σ}}(x_j-x_i)\·{(z_i^k\left(0\right)+\▿z_i^k{\left(0\right)}^T\·s_j^k)}^2$

其中，

被定义为 $\▿z_i^k\left(s\right)=-{\▿J}_{2,k}(x_i+s).$

通过计算E(s_j ^k)的第一变型，最优性的必要条件

$\∀h,\mathrm{dE}\left(s_j^k\right)\·h=0$

产生s_j ^k的递归、闭合形式解(closed-form solution)：

$s_j^k=-{(\underset{t}{\mathrm{\Σ}}{G}_{\mathrm{\σ}}(x_j-x_i)\▿z_i^k\left(0\right)\▿z_i^k{\left(0\right)}^T)}^{-1}(\underset{i}{\mathrm{\Σ}}{G}_{\mathrm{\σ}}(x_j-x_i)\▿z_i^k\left(0\right)z_i^k\left(0\right)).$

注意到的是，矩阵

是对称的且可以由6个系数来表征。矢量

由3个系数定义。

一旦计算出了这些系数，可以看出，估计位移更新减少到执行通过高斯核对矩阵和矢量系数进行卷积，并且减少到求解每个网格点处的3×3对称线性问题。

卷积使用如在R.Deriche的文献“Fast algorithms for low-level vision”(IEEE Transactions on Pattern Analysis andMachine Intelligence，1(12)，1990年，第78-88页)所提出的二阶递归滤波器来进行近似。这些滤波器减少了近似高斯平滑及其导数所要求的计算工作量。利用这种方法，与滤波器大小无关地对每个输出点进行固定数目的乘法和加法来执行这些操作。

递归、闭合形式解被用于3×3线性问题。如果 $\underset{i}{\mathrm{\Σ}}{G}_{\mathrm{\σ}}(x_j-x_i)\▿z_i^k\left(0\right)\▿z_i^k{\left(0\right)}^T$ 由下面的对称矩阵给出：

$A=\left(\begin{array}{ccc}{a}_1& a_2& a_3\\ a_2& a_4& a_5\\ a_3& a_5& a_6\end{array}\right)$

那么其逆矩阵为

$A^{-1}=\frac{1}{a_6{a}_1{a}_4-a_1{a_5}^2-a_6{a_2}^2+2a_2{a}_3{a}_5-a_4{a_3}^2}\left(\begin{array}{ccc}{c}_1& c_2& c_3\\ c_2& c_4& c_5\\ c_3& c_5& c_6\end{array}\right)$

其中

$\left\{\begin{array}{c}{c}_1=a_6{a}_4-{a_5}^2\\ c_2=a_3{a}_5-a_6{a}_2\\ c_3=a_2{a}_5-a_3{a}_4\\ c_4=a_1{a}_6-{a_3}^2\\ c_5=a_3{a}_2-a_1{a}_5\\ c_6=a_1{a}_4-{a_2}^2\end{array}\right.$

在前面所有的计算中，当在网格点之外需要数值时，使用三线性内插法。

特征映射计算的描述如下。在乳腺MR中，造影剂吸入可能明显地或局部地改变了所观察的组织的强度值。不依赖于针对配准的强度，取函数J₁和J₂的值作为原始数据的拉普拉斯算子。这个想法主要集中在图像中的边缘信息(表征边缘的拉普拉斯算子的过零)上，这很少受造影剂吸入的影响。该方法被证明在实践中非常有效。现在使用笛卡尔坐标系统用三个一组的整数(i₁，i₂，i₃)来标识体素，并且下面的有限差分法被用于拉普拉斯算子：

$J_2(i_1,i_2,i_3)=\frac{I_2(i_1,i_2,i_3+1)}{{\mathrm{\λ}}^2}+\frac{I_2(i_1,i_2,i_3-1)}{{\mathrm{\λ}}^2}+I_2(i_1,i_2+1,i_3)+I_2(i_1,i_2-1,i_3)$

$+I_2(i_1+1,i_2,i_3)+I_2(i_1-1,i_2,i_3)-(4+\frac{2}{{\mathrm{\λ}}^2})I_2(i_1,i_2,i_3)$

注意的是，被用来估计拉普拉斯算子的中心差分法考虑体素的各向异性。系数λ对应于沿第三轴的体素间距(切片之间的距离)与沿前两个轴的像素间距(这里被假设为是相同的)之间的比率。

同样的等式用于由I₁计算J₁。

如将会是显而易见的那样，本发明意图使用和应用编程数字计算机来实现。图2以基本示意的形式示出了被耦合用于与输入装置、输出装置、以及用于存储程序及其它数据的存储装置进行双向数据通信的数字处理器。输入装置用广泛的术语被指定为用于提供适当的一个或多个图像的装置，这些图像用于根据本发明进行处理。例如，输入可以来自于成像装置，诸如来自被并入CATSCAN、x射线机器、MRI或其它装置中的装置，或该输入是所存储的图像，或该输入通过直接连接、调制过的红外射束、无线电、陆上线路、传真、或卫星(例如，通过万维网或因特网)与另一计算机或装置进行通信，或该输入是任何其它适当的这种数据源。输出装置可以包括使用任何适当设备(诸如阴极射线显象管、等离子显示器、液晶显示器等等)的计算机型显示装置，或该输出装置可以或可以不包括用于再现图像的装置且可以为了方便包括图2中的存储装置或部分该存储装置，该存储装置用于存储图像以进一步进行处理或进行查看或评估，或该输出装置可以利用包括如上面提到的与输入装置有关的连接或耦合。处理器运行根据本发明所建立的用于实现本发明的步骤的程序。这种被编程的计算机易于经通信介质(诸如陆上线路、无线电、因特网等)连接，用于采集和传输图像数据。

本发明可以容易地(至少部分)用软件存储装置实现或被封装在形成为软件产品的软件存储装置中。这可以是计算机程序产品的形式，该计算机程序产品包括计算机可用介质，该计算机可用介质具有针对程序代码而被记录于其上的计算机程序逻辑，用于利用本发明的方法执行图像运动补偿。

虽然本发明已经通过使用图示说明的示例性实施例的例子进行了解释，该实施例涉及MRI检测人体乳房的潜在肿瘤中的图像时间序列的运动补偿，但是本发明通常也可被用于解决要求其它场中的空间对准(诸如，但不限于，PET-CT配准的例子)的问题。

应理解，通过示例性实施例的说明不意图来限制，而在不脱离其范围由随后的权利要求所限定的本发明的精神的情况下，也可以进行在此没有明确描述的各种变化及置换。

Claims (43)

1.一种用于在时间序列中的第一与第二图像之间进行运动补偿的方法，所述方法包括：

由所述第一与第二图像导出相应的第一与第二特征映射；

通过根据配准算法处理所述第一与第二特征映射而导出位移场，所述配准算法包括针对每个像素或体素求解局部高斯加权最小均方问题，以便导出形成所述位移场的相应矢量；以及

利用所述位移场弯曲所述第二图像。

2.根据权利要求1所述的用于进行运动补偿的方法，其中，所述导出位移场的步骤利用前面导出的位移场来导出所述位移场。

3.根据权利要求2所述的用于进行运动补偿的方法，其中，所述导出位移场的步骤针对最初导出所述位移场而利用给定的初始位移场。

4.根据权利要求3所述的用于进行运动补偿的方法，其中：所述给定的初始位移场的缺省条件是空集。

5.根据权利要求3所述的用于进行运动补偿的方法，其中：所述给定的初始位移场考虑患者运动的先验知识。

6.根据权利要求2所述的用于进行运动补偿的方法，包括重复导出位移场的步骤，其中在对应于比前次重复具有更高分辨率的所述第一与第二图像的第一与第二特征映射上执行每次重复。

7.根据权利要求2所述的用于进行运动补偿的方法，其中，每次重复利用由来自直接在前步骤的之前导出的位移场。

8.根据权利要求7所述的用于进行运动补偿的方法，其中，所述导出位移场的步骤针对最初导出所述位移场而利用给定的位移场。

9.根据权利要求7所述的用于进行运动补偿的方法，其中：

所述导出位移场的步骤包括将所述前面导出的位移场扩展到所述分辨率增加的版本的分辨率水平。

10.根据权利要求7所述的用于进行运动补偿的方法，包括利用为了弯曲所述第二图像而存在的、在最高分辨率导出的位移场来产生运动校正过的图像。

11.根据权利要求7所述的用于进行运动补偿的方法，其中：

所述导出相应的第一与第二特征映射的步骤包括导出所述相应的分辨率减小的第一与第二图像中的图像数据的相应拉普拉斯算子。

12.一种用于在时间序列中的第一与第二图像之间进行运动补偿的方法，所述方法包括：

用减小过程处理所述第一与第二图像，用于提供相应的分辨率减小的第一与第二图像；

由所述相应的分辨率减小的第一与第二图像导出相应的第一与第二特征映射，所述特征映射包括导出所述相应的分辨率减小的第一与第二图像中的图像数据的相应的拉普拉斯算子；

通过根据配准算法处理所述第一与第二特征映射而导出位移场，所述配准算法包括对每个像素或体素求解局部高斯最小均方问题，以便导出形成该位移场的相应矢量；以及

利用所述位移场弯曲该第二图像。

13.一种用于在时间序列中的第一与第二图像之间进行运动补偿的方法，所述方法包括：

用减小过程处理所述第一与第二图像，用于提供相应的分辨率减小的第一与第二图像；

由所述相应的分辨率减小的第一与第二图像导出相应的第一与第二特征映射；

通过根据配准算法处理所述第一与第二特征映射而导出位移场，所述配准算法包括对每个像素或体素求解局部高斯加权最小均方问题，以便导出形成所述位移场的相应矢量；以及

利用所述位移场弯曲所述第二图像。

14.根据权利要求13所述的用于进行运动补偿的方法，其中：

所述导出相应的第一与第二特征映射的步骤包括导出所述相应的分辨率减小的第一与第二图像中的图像数据的相应的拉普拉斯算子。

15.根据权利要求13所述的用于进行运动补偿的方法，其中：

所述导出位移场的步骤利用初始位移场。

16.根据权利要求15所述的用于进行运动补偿的方法，其中：

所述初始位移场的缺省条件是空集。

17.根据权利要求15所述的用于进行运动补偿的方法，其中：所述初始位移场考虑患者运动的先验知识。

18.根据权利要求13所述的用于进行运动补偿的方法，其中：

所述根据配准算法处理所述第一与第二特征映射的步骤包括将所述位移场扩展到与所述第二图像的分辨率水平相适合的分辨率水平的步骤。

19.一种用于在时间序列中的参考图像与浮动图像之间进行运动补偿的方法，所述方法包括：

由所述参考图像导出分辨率逐渐减少的第一组修正图像；

由所述浮动图像导出分辨率逐渐减少的第二组修正图像；

由所述第一组修正图像导出第一组特征映射；

由所述第二组修正图像导出第二组特征映射；

根据配准算法，从分别第一与第二组特征映射中的每个特征映射中的最低分辨率组成部分以及给定的初始位移场导出第一位移场，所述配准算法包括针对每个像素或体素求解局部高斯最小均方问题，以便导出形成所述位移场的相应矢量；

根据所述配准算法，从分别第一与第二组特征映射中的每个特征映射中的次最低分辨率组成部分以及在前一步骤中所获得的所述第一位移场导出第二位移场；

如果第一与第二组特征映射中的每个特征映射中存在分辨率连续较高的组成部分，则针对这些组成部分重复前述步骤，并且在每种情况下都使用在当前步骤之前的步骤中所获得的位移场，直至达到所述浮动图像的分辨率，如果不存在分辨率较高的组成部分，则直接进行到下一步；以及

利用最后获得的位移场弯曲所述浮动图像。

20.根据权利要求19所述的用于进行运动补偿的方法，其中：

所述导出特征映射的步骤包括导出相应图像中的图像数据的拉普拉斯算子。

21.根据权利要求19所述的用于进行运动补偿的方法，其中：

所述导出位移场的步骤包括以下步骤：如果需要扩展，则将所述位移场扩展到与所述次较高分辨率水平相适合的分辨率水平。

22.根据权利要求19所述的用于进行运动补偿的方法，其中：

所述给定的初始位移场的缺省条件是空集。

23.根据权利要求19所述的用于进行运动补偿的方法，其中：

所述给定的初始位移场考虑患者运动的先验知识。

24.一种用于在时间序列中的参考图像与浮动图像之间进行运动补偿的方法，所述方法包括：

(a)用相应的第一与第二多个级联的分辨率减小过程处理所述参考图像与浮动图像，用于提供相应多个分辨率连续减小的参考图像与浮动图像，在此对于最低分辨率水平称为水平0，而对于次较高分辨率水平称为水平1，对第二次较高分辨率水平称为水平2，对任何现有的较高分辨率水平以此类推；

(b)由分辨率水平为L0、L1、L2等的所述相应多个分辨率连续减小的参考图像与浮动图像导出相应的多个参考特征映射与浮动特征映射；

(c)通过根据配准算法处理与水平L0相对应的参考特征映射与浮动特征映射而导出第一位移场，所述配准算法包括基于给定的初始位移场针对每个像素或体素求解局部高斯加权最小均方问题，以便导出形成所述第一位移场的相应矢量；

(d)将所述第一位移场扩展到与分辨率L1相适合的分辨率水平，以提供被扩展的第一位移场；

(e)通过根据所述配准算法基于所述被扩展的第一位移场来处理与分辨率水平L1相对应的参考特征映射与浮动特征映射而导出第二位移场，以便导出形成所述第二位移场的相应矢量；

(f)将所述第二位移场扩展到与分辨率L1的位移场相适合的分辨率水平，以提供被扩展的第二位移场；

(g)通过根据所述配准算法基于所述被扩展的第二位移场来处理与分辨率水平L2相对应的参考特征映射与浮动特征映射而导出第三位移场，以便导出形成所述第二位移场的相应矢量；以及

(h)如果L2是所述参考图像与浮动图像的分辨率水平，那么通过利用所述第三位移场来弯曲所述浮动图像并结束；而如果L2不是所述参考图像与浮动图像的分辨率水平，则

(i)将所述第三位移场扩展到与次较高分辨率水平相适合的分辨率水平，以提供被扩展的第三位移场；以及

(j)重复从步骤(g)开始的一系列步骤，适当修改分辨率水平，以基于前面最后获得的被扩展的位移场而获得与越来越高的分辨率水平相对应的位移场，直到达到所述参考图像与浮动图像的分辨率水平，并因此利用所获得的最后位移场弯曲所述浮动图像并结束。

25.根据权利要求24所述的用于进行运动补偿的方法，其中：所述提供相应多个分辨率连续减小的参考图像与浮动图像的步骤包括导出所述相应多个分辨率连续减小的参考图像与浮动图像中的图像数据的相应的拉普拉斯算子。

26.根据权利要求24所述的用于进行运动补偿的方法，其中：所述导出第一位移场的步骤包括输入给定的初始位移场。

27.一种用于在时间序列中的第一与第二图像之间进行运动补偿的方法，所述方法包括：

由相应的分辨率减小的第一与第二图像导出相应的第一与第二特征映射；

通过根据配准算法处理所述第一与第二特征映射而导出第一位移场，所述配准算法包括针对每个像素或体素求解局部高斯加权最小均方问题，以便导出形成所述位移场的相应矢量；

扩展所述第一位移场，以与所述第一与第二图像的分辨率相符合；

通过根据所述配准算法基于所述第一位移场来处理所述第一与第二图像而导出第二位移场；以及

利用所述第二位移场弯曲所述第二图像，以便获得运动校正过的图像。

28.根据权利要求27所述的用于进行运动补偿的方法，其中：所述导出第一位移场的步骤包括输入给定的初始位移场的步骤。

29.根据权利要求28所述的用于进行运动补偿的方法，其中：所述给定的初始位移场的缺省条件是空集。

30.根据权利要求28所述的用于进行运动补偿的方法，其中：所述给定的初始位移场考虑患者运动的先验知识。

31.根据权利要求27所述的用于进行运动补偿的方法，其中：所述导出相应的第一与第二特征映射的步骤包括导出所述相应的分辨率减小的第一与第二图像中的图像数据的相应的拉普拉斯算子。

32.一种用于在时间序列中的第一与第二图像之间进行运动补偿的方法，所述方法包括：

从所述第一图像导出具有逐渐减小的分辨率的第一系列图像；

从所述第二图像导出具有逐渐减小的分辨率的第二系列图像；

从所述第一系列图像导出第一系列特征映射；

从所述第二系列图像导出第二系列特征映射；

通过根据用于提供相应位移场的配准算法来处理具有所述第二系列特征映射的相对应的特征映射的所述第一系列特征映射的特征映射而导出系列位移场，所述第二系列特征映射与所述第一系列特征映射属于相同分辨率，所述配准算法包括针对每个像素或体素求解局部高斯加权最小均方问题，以便导出形成所述位移场的相应矢量，以给定的位移场开始，该给定的位移场被用来导出分辨率水平最低的第一位移场，并且此后利用直接在前的位移场来导出与所述系列中的次较高分辨率相对应的下一位移场；以及

利用在所述第一与第二图像的分辨率水平导出的最后的位移场来弯曲所述第二图像。

33.根据权利要求32所述的用于进行运动补偿的方法，包括：

除了所述给定场以及所述最后场之外，扩展所述每个分辨率场，以与次较高分辨率相符合。

34.一种用于在时间序列中的参考图像与浮动图像之间进行运动补偿的方法，所述方法包括：

(a)用相应的第一与第二多个级联的分辨率减小过程处理所述参考图像与浮动图像，用于提供相应多个分辨率连续减小的参考图像与浮动图像，在此对于最低分辨率水平称为水平0，而对于次较高分辨率水平称为水平1，对于第二次较高分辨率水平称为水平2，对任何现有的较高分辨率水平以此类推；

(b)由按照分辨率增加的顺序排列的、分辨率水平为L0、L1、L2的所述相应多个分辨率连续减小的参考图像与浮动图像导出相应多个参考特征映射与浮动特征映射，其中L2是所述参考图像与浮动图像的分辨率水平；

(c)通过根据配准算法处理与水平L0相对应的参考特征映射与浮动特征映射而导出第一位移场，所述配准算法包括基于给定的初始位移场对每个像素或体素求解局部高斯加权最小均方问题，以便导出形成所述第一位移场的相应矢量；

(d)将所述第一位移场扩展到与分辨率L1相适合的分辨率水平，以提供被扩展的第一位移场；

(e)通过根据所述配准算法基于所述被扩展的第一位移场来处理与分辨率水平L1相对应的参考特征映射与浮动特征映射而导出第二位移场，以便导出形成所述第二位移场的相应矢量；

(f)将所述第二位移场扩展到与分辨率L2相适合的分辨率水平，以提供被扩展的第二位移场；

(g)通过根据所述配准算法基于所述被扩展的第二位移场处理与分辨率水平L2相对应的参考特征映射与浮动特征映射来导出最后的位移场，以便导出形成所述最后的位移场的相应矢量；以及

(h)通过利用所述最后的位移场弯曲所述浮动图像。

35.根据权利要求34所述的用于进行运动补偿的方法，其中：所述给定的初始位移场的缺省条件是空集。

36.根据权利要求34所述的用于进行运动补偿的方法，其中：所述给定的初始位移场考虑患者运动的先验知识。

37.根据权利要求34所述的用于进行运动补偿的方法，其中：

所述导出特征映射的步骤包括导出所述相应的分辨率减小的第一与第二图像中的图像数据的相应的拉普拉斯算子。

38.根据权利要求34所述的用于进行运动补偿的方法，其中：

导出第一位移场的步骤(c)包括以下步骤：

用下标j标记每个体素，该体素的坐标由矢量x_j∈R³给出；

将 $p_j^k\∈R^3$ 定义为所述算法进行第k次迭代之后在x_j点处重新获得的位移矢量；

给定一组位移矢量p_j ^k，由J_2，k表示在迭代k时所述浮动特征映射的校正过的版本：

$\∀j,J_{2,k}\left(x_j\right)=J_2(x_j+p_j^k);$

在每个体素j用s_j ^k将更新规则 $p_j^{k+1}=p_j^k+s_j^k$ 定义为非线性最小均方问题的解：

$\underset{s_j^k}{\min }\underset{i}{\mathrm{\Σ}}{G}_{\mathrm{\σ}}(x_j-x_i)\·{\left(z_i^k\left(s_j^k\right)\right)}^2$

其中 $z_i^k\left(s\right)=J_1\left(x_i\right)-J_{2.k}(x_i+s)$ 和G_σ(x)是标准偏差σ的各向同性三维高斯核：

$G_{\mathrm{\σ}}\left(x\right)=\frac{1}{\sqrt{{\left(2\mathrm{\π}\right)}^3{\mathrm{\σ}}^6}}{e}^{\frac{{\left|s\right|}^2}{{\mathrm{\σ}}^2}};$

通过最小化下列标准来求解所述非线性最小均方问题的线性化版本：

$E\left(s_j^k\right)=\underset{i}{\mathrm{\Σ}}{G}_{\mathrm{\σ}}(x_j-x_i)\·{(z_i^k\left(0\right)+\▿z_i^k{\left(0\right)}^T\·s_j^k)}^2$

其中，

被定义为 $\▿z_i^k\left(s\right)=-\▿J_{2,k}(x_i+s);$

计算E(s_j ^k)的第一变型，最优性的必要条件

$\∀h,\mathrm{dE}\left(s_j^k\right)\·h=0$

以产生s_j ^k的递归、闭合形式解：

$s_j^k=-{(\underset{i}{\mathrm{\Σ}}{G}_{\mathrm{\σ}}(x_j-x_i)\▿z_i^k\left(0\right)\▿z_i^k{\left(0\right)}^T)}^{-1}(\underset{i}{\mathrm{\Σ}}{G}_{\mathrm{\σ}}(x_j-x_i)\▿z_i^k\left(0\right)z_i^k\left(0\right)).$

39.一种用于执行图像运动补偿的系统，其包括：

存储装置，用于存储程序及其它数据；以及

与所述存储装置进行通信的处理器，所述处理器运行所述程序来执行：

用于在时间序列中的第一与第二图像之间进行运动补偿的方法，所述方法包括：

用减小过程处理所述第一与第二图像，用于提供相应的分辨率减小的第一与第二图像；

由所述相应的分辨率减小的第一与第二图像导出相应的第一与第二特征映射；

通过根据配准算法处理所述第一与第二特征映射而导出位移场，所述配准算法包括对每个像素或体素求解局部高斯加权最小均方问题，以便导出形成所述位移场的相应矢量；以及

利用所述位移场弯曲所述第二图像。

40.一种用于在时间序列中的第一与第二图像之间进行运动补偿的方法，所述方法包括：

通过根据一过程处理所述第一与第二图像而导出位移场，所述过程包括配准算法，所述配准算法包括针对每个像素或体素求解局部高斯加权最小均方问题，以便导出形成所述位移场的相应矢量；以及

利用所述位移场弯曲所述第二图像。

41.根据权利要求40所述的用于进行运动补偿的方法，其中，所述根据一过程处理所述第一与第二图像的步骤包括分别由所述第一与第二图像导出相应的第一与第二特征映射以及将所述配准算法用于所述特征映射的步骤。

42.根据权利要求40所述的用于进行运动补偿的方法，其中，所述导出位移场的步骤利用前面导出的位移场来导出所述位移场。

43.根据权利要求40所述的用于进行运动补偿的方法，其中，所述导出位移场的步骤针对最初导出所述位移场而利用给定的初始位移场。

Images