CN101133431A

CN101133431A - 能够减少物体运动造成的成像伪影的生物医学图像配准方法和计算机程序产品

Info

Publication number: CN101133431A
Application number: CNA200680003997XA
Authority: CN
Inventors: 托尼·维尔纳·沃姆维格; 海纳·法贝尔; 迪尔克·迈耶
Original assignee: Bracco Imaging SpA
Current assignee: Bracco Imaging SpA
Priority date: 2005-02-03
Filing date: 2006-01-31
Publication date: 2008-02-27
Anticipated expiration: 2026-01-31
Also published as: WO2006082192A1; CN101133431B; EP1844440B1; US20080292214A1; US8073290B2; JP2008538080A; ATE553456T1; EP1844440A1

Abstract

公开了一种用于减少由物体运动造成的成像伪影的生物医学图像的配准方法，在该方法中，在第一幅和第二幅待配准图像中自动判定并跟踪一定数量的特征或标志，以确定上述第一幅和第二幅图像之间的光流矢量。随后通过对第二幅图像施加逆光流进行配准。上述自动特征选择步骤是通过对每个像素确定其定义的邻域内的平均信号强度，随后将该平均信号强度与预先确定的阈值进行比较来实现的。如果所述邻域的平均信号强度值大于预先确定的阈值，那么对应的像素就被定义为特征，并添加到要跟踪的特征的列表中。

Description

能够减少物体运动造成的成像伪影的生物医学图像配准方法和计算机程序产品

技术领域

本发明涉及图像配准的方法，该方法由以下步骤组成：

a)提供同一物体的至少第一幅和第二幅数字或数字化图像，或者至少第一组和第二组该物体的截面图像，这些图像由像素或体素的二维或三维阵列形成；

b)在第一幅图像或第一组图像中，通过选取一定数量的像素或体素，将这些像素或体素设为标志或特征，来定义一定数量的标志，即所谓的特征，并生成要跟踪的所述特征的列表；

c)通过对选定为特征的每个像素或体素确定从第一幅到第二幅图像或从第一组到第二组图像的光流矢量，从第一幅到第二幅图像或从第一组到第二组图像，跟踪选定为特征的每个像素或体素的位置；

d)通过对第二幅或第二组图像的像素或体素施加逆光流，来对第一幅和第二幅图像或第一组和第二组图像进行配准；

背景技术

这种为时间序列的部分的图像的配准方法可以从例如文献US6,553,152中获悉。

文献US6,553,152提供了用于两幅图像的配准方法，该方法是通过由操作员对标志的目视识别以及对两幅图像上的对应像素进行标记来完成的。

这种方法使得标志的选择完全依赖操作者的技能，并且当图像中不存在显著的或意义明确的可识别结构时，该方法在实践中非常难于实现。考虑核磁共振成像(MRI)成像或超声图像或X射线图像之类的诊断图像可知，根据这些图像的结构特征进行对于标志的目视识别是极端困难的，并且可能造成重大的错误。

许多不同的算法试图识别轮廓，如Fischer，B.；Modersitzki，J.Curvature Based Registration with Applications to MRMammography；Springer-Verlag Berlin Heidelberg：ICCS 2002，LNCS 2331，2002；pp202-206所公开的算法，或者试图识别被测量过两次或多次或者甚至用多种医疗器械进行过研究的器官或部位的放射线图像内的标志，例如Shi J，Tomasi C.Good Features to Track.1994 IEEE Conference onComputer Vision and Pattern Recognition(CVPR’94)1994；593-600、Tommasini T，Fusiello A，Trucco E，Roberto V.Making Good FeaturesTrack Better Proc.IEEE Int.Conf.on Computer Vision and PatternRecognition(CVPR’98)1998；178-183、以及Ruiter，N.V.；Muller，T.O.；Stotzka，R.；Kaiser，W.A.Elastic registration of x-ray mammograms andthree-dimensional magnetic resonance imaging data.J Digit Imaging 2001，14，52-55中所公开的算法。

对在不同时刻拍摄的同一物体的平面图像进行配准以计算物体随时间的运动，这种配准称为二维配准。实现相同功能但在例如MRI或计算断层(CT)扫描仪的截面图像组的三维图像组上作用的算法，被称为三维配准。在这些图像组内，标志或物体的运动可能在任何方向上发生。

根据要比较的不同成像方式的数量，配准算法可以分为单模式配准算法和多模式配准算法。

比较两个MRI序列就是一种单模式配准。

通常，单模式配准算法设法识别两幅图像内的一定标志或复杂形状，大多数情况下处理的是二维图像。在第二步则比较两幅图像内对应标志的位置，并计算其运动矢量。这个信息用来将第二幅图像的像素移回到新的位置，以消除运动伪影。

“刚性配准”将二维或三维的图像/体积作为单个单元整体移动到一定方向。“非刚性配准”则仅仅将一定面积/体积的单个或多个像素/体素移往不同的方向。

考虑很柔软的乳房组织这一例子可知，需要非刚性配准来消除运动伪影。

由于乳房不具有任何像骨头一样的成形结构，这样特殊的一致性使得它的直接相邻的部分能够在不同的方向上移动。因此，显而易见的是，任何配准算法或任何其他现有的方法都必须在整个乳房组织各处识别和散布其标志。

有一组算法定义了必须寻找的一定数量的标志。它们将按其效价(valence)分类。对于乳房MR图像而言，由于MRI中脂肪/空气之间的较高的对比度，这些方法常常会在乳房的外边界处和胸腔结构处找到许多标志(Huwer，S.；Rahmel，J.；Wangenheim，A.v.Data-DrivenRegistration for Lacal Deformations.Pattern Recognition Letters 1996，17，951-957)。

使用压缩板或者至少通过俯卧压在双乳上，通常得到皮肤的所有部分的良好固定，并使其与压缩板直接接触(图1)。相反地，双乳的内部区域因心跳和呼吸而继续轻微移动。由于与外部皮肤/空气边界比较而言脂肪与软组织间的对比度低得多，因此这些算法中乳房中心内的有效标志数保持较少。

标志选取的一般问题存在于如下的事实中，即选取过少的标志可能导致不充分或不精确的配准。但是，选取额外的标志并不一定能保证精确的配准，这是因为观察人员或者任何程序可能不得不包括具有较低的重分配确定性(certainty of being reallocated)的标志，举例而言，可以由乳房中心的不同组织类型的低对比度导致该较低的重分配确定性。增加标志数总是会显著增大配准的计算复杂度。

一种特别适用于所谓的非刚性配准的特征选择和特征跟踪算法是所谓的Lucas&Kanade算法及其特定的金字塔式实现，该算法及其实现在Lucas BD，Kanade T.An Iterative Image Registration Techniquewith an Application to Stereo Vision，IJCAI81；674-679 1981和BouguetJY.Pyramidal Implementation of the Lucas Kanade Feature Tracker，TechReport as Part of Open Computer Vision Library Documentation，Intel Corporation enclosed hereinafter as appendix 1中进行了详细的披露。Lucas&Kanade技术及其金字塔式特征跟踪实现特别适合用于自动识别图像中的可靠标志，以及适用于在不同时刻拍摄的图像之间跟踪这些标志。选定特征的偏移或位移就确定为位移或所谓的光流矢量。这种技术得到了广泛的应用，例如，用于机器人学中的计算机视觉，并且成为一些基本课程中所教的数字图像处理和计算机视觉领域内的技术人员的常识的部分。

但是，当考虑将Lucas & Kanade技术应用到数字的或数字化的诊断图像中，更特别的是应用到乳房(诊断)这样的实际用途时，所述的技术会识别太多特征，当然也包括位于兴趣区之外的太多的特征，并且在譬如包围乳房的空气区易受噪声的干扰。

从下面的例子中可以更详细地看到上述问题。

对比度增强(CE)核磁共振乳腺成像(MRM)是一种有益且重要的研究，对于侵略性肿瘤具有极高的灵敏度。它表现出与乳房软组织无关的较高的负预测值。近来的多中心研究表明，灵敏度和特异性随图像分析技术的不同而不同。利用1.5T扫描仪，分别计算出96％、93％或86％的灵敏度及30％、50％或70％的特异性。因此，MRM的主要缺点仍是较高的假阳性(false positive)诊断百分率。

通常，MR成像在1.0或1.5特斯拉成像仪上完成，并且使用了厂家的双乳线圈(在美国则经常使用单乳线圈设备)，而患者处于俯卧姿势。该治疗方案至少包含任意方向上的动态对比度增强T1加权序列，且采样时间在30s到180s之间。在快速浓注钆喷酸二葡甲胺盐或任何其他的顺磁性MR造影剂之前，至少进行一次测量，注射之后则要进行多次测量。

为了计算乳房组织的一个区域的实际造影剂摄取量，必须将施加造影剂之前的每个体素的信号强度从施加了造影剂之后的信号强度中减去。由于呼吸和心跳造成双乳的最小限度的运动以及MR图像的切片厚度带来的部分体积效应，不同时刻拍摄的同一图像位置的体素并不确切示出同一块乳房组织。因此，减影图像并不完全漆黑一片(除肿瘤外)。微小运动效应可以在乳房的外边界处最清楚地展现出来。由于脂肪组织信号强度较高，而周围空气的信号强度约为零，非常小的运动将脂肪置于空气的之前图像位置处。从脂肪的高信号强度中减去空气的零信号强度，不恰当地模仿了CM应用所引起的高信号增量。结果，减影图像在成像乳房的至少部分处示出表示造影剂大摄取量的粗白线(图2：乳房的白色外边界利用符号进行标记)。

考虑笛卡尔系统定义的三维或体积图像采集，除了沿x和y方向上的任意运动之外，人们总能发现沿z方向的某种运动(朝胸部方向)。这是由研究过程中患者的松弛造成的。

为了减少伪影，通常借助某种压缩方法(因厂家而异)将乳房固定在乳房线圈内。但是，尽管这样，最小限度的运动总会发生。因此，纵然对乳房加以固定，其图像在x、y和z方向总会显现轻微的运动伪影，它们在减影图像中看来是明亮的像素。如果没有运动伪影，那么减影图像除了肿瘤组织所在的区域外应该完全呈黑色。

在乳房外周围空气的噪声信号中可检测到一定数量的标志。除了必须对考虑的每个特征进行跟踪并带来不必要的计算负担这一事实外，由于噪声是个随机现象，噪声中发现的所有特征或标志在第二幅或随后的图像中将没有对应的标志。但是事实上，在第二幅或随后的图像或图像组中搜索与之前在噪声中发现的标志相对应的标志的算法，将通过某种方式在第二幅或第二组图像中重分配一定量的标志。这就会导致错误的图像配准，并对结果产生影响。

发明内容

本发明的目的在于提供一种配准图像尤其是三维结构的生物医学诊断图像的方法，它通过跟踪三维运动的能力，通过减少数字或数字化的图像中兴趣区之外的噪声所引起的伪影，以及通过与相邻组织类型的高或低信号强度无关地在整个乳房区各处扩散标志，来有效克服已知配准方法的缺陷。

依照本发明，图像配准方法由以下步骤组成：

a)提供通过对于组织或组织区或解剖区进行MRI、超声或X射线扫描而获取的第一幅数字或数字化的图像或第一组截面图像；所述图像由像素或体素的二维或三维阵列构成；提供通过在第二时刻对相同组织或组织区或解剖区进行MRI、超声或X射线扫描而获取的相同解剖区的至少第二幅数字或数字化的图像或第二组截面图像，在获取图像时可选的在所述组织或组织区或所述解剖区存在对比媒质；

b)在第一幅图像或第一组图像中通过选取一定数量的像素或体素，将该像素或体素设为标志或特征，来定义一定数量的标志，即所谓的特征，并生成要跟踪的所述特征的列表；

c)通过对每个选定为特征的像素或体素确定从第一幅到第二幅图像或从第一组到第二组图像的光流矢量，从第一幅到第二幅图像或从第一组到第二组图像跟踪选定为特征的每个像素或体素的位置；

e)在步骤d中配准之后，从第二幅图像的第二图像数据减去第一幅图像的图像数据；

f)显示通过从第二幅图像的第二图像数据减去第一幅图像的图像数据所得到的图像数据。

步骤b“特征选择”可以如下更详细地描述：

B1)在第一幅图像或第一组截面图像的每个像素或体素的周围定义像素或体素邻域，所述的像素或体素邻域包括有限数量的像素或体素；对于单幅图像，选用二维邻域；对于截面图像组，则选用三维邻域；

B2)对于每个像素或体素，确定所述像素或体素邻域的所有像素或体素的平均信号强度值；

B3)定义平均信号强度值阈值；

B4)比较步骤B2中确定的每个像素或体素邻域的平均信号强度，并将所述平均信号强度值与平均信号强度值阈值比较；

B5)如果步骤B4中所述邻域的平均信号强度值大于阈值，则将所述像素或体素定义为要跟踪的特征，并将其添加到要跟踪的特征的列表中。

所述平均信号强度阈值是可调的，可以发生变化。经验的或实验的标准可以用来确定这种阈值，例如，在一组测试或样本图像上执行本方法并比较取不同阈值时获得的结果，可以确定该阈值。

尽管按照本发明的方法，可以应用任何用于选择可被定义成要跟踪的图像特征的像素或体素的方法，使用根据著名的Lucas & Kanade方法的自动特征选择方法却能得到最好的结果，Lucas & Kanade方法在Lucas BD，Kanade T.An Iterative Image Registration Technique withan Application to Stereo Vision，IJCAI81；674-679 1981中进行了更为详细的描述。

利用位移或偏移矢量，即所谓的光流，来跟踪特征从第一幅图像到第二幅图像的偏移量，也是已知的方法。在这种情况下，本发明也可应用到利用了标志或特征的一些跟踪方法中。尤其是如果考虑了诊断成像期间受限乳房的微小移动和Lucas & Kanade特征选择方法的使用这种特殊情况，本发明优选地采用所谓的“Lucas & Kanade特征跟踪器的金字塔式实现”(见附录1)更好。这种实现为普通技术人员所熟悉，并在一些大学课程中被讲授。Lucas & Kanade特征跟踪器的这种金字塔式实现的精确而详细的描述在Bouguet J-Y.PyramidalImplementation of the Lucas Kanade Feature Tracker，Tech Report as Partof Open Computer Vision Library Documentation，Intel Corporation、以及Shi J，Tomasi C.Good Features to Track.1994 IEEE Conference onComputer Vision and Pattern Recognition(CVPR’94)1994；593-600、以及Tommasini T，Fusiello A，Trucco E，Roberto V.Making Good FeaturesTrack Better.Proc.IEEE Int.Conf.on Computer Vision and PatternRecognition(CVPR’98)1998；178-183中进行了公开。

一旦对所有特征定义了光流矢量，即可通过对第二幅图像施加逆光流，以便将第二幅图像的每个特征及其周围像素移回到所述特征在第一幅图像中的位置，可以实现两幅图像的配准。

前面已经提到，上述方法可以用来配准两幅诊断图像或三维图像体积，这些图像核磁共振成像(MRI)或计算X射线断层成像(CT)中用相同或不同类型的扫描仪获取。

上述公开的配准方法还与用于MRI或超声图像之类的对比媒质增强诊断成像尤其是对比度增强核磁共振乳腺成像(MRM)的方法结合地提供。在该方法中，为了计算乳房组织的一个区域的实际造影剂摄取量，必须将造影剂施加之前的每个体素的信号强度从造影剂施加之后的信号强度中减去。因此，将某时刻拍摄的、成像组织中不存在造影剂的第一幅图像与第二或稍后时刻拍摄的、造影剂存在且灌注在成像组织中的第二幅图像配准，有助于在很大程度上成功地抑制伪影，如果不进行配准，这些伪影将因为患者的运动而出现在所述的减影图像数据中。这个简单的思想是，把所述的两幅图像相减，将使得与造影剂汇集的组织区相对应的像素或体素的强度水平很高。所有其他不存在造影剂的区看起来将是一些低强度像素或体素，且具体地呈暗灰色或黑色。

进一步的改进是从属权利要求的主题。

附图说明

从下文对本方法的更详细的描述以及所附绘图和图表可更清楚地看到本发明的特征及派生优点，其中：

图1示出了乳房MRI中使用的乳房固定方法的示意性例子。

图2示出减运算过程以及不执行配准步骤时出现的运动伪影：左上和右上部分是同一MR扫描仪位置的两幅乳房截面图。测量第一幅图像后约2.5分钟施加造影剂，施加造影剂之后测量第二幅图像。两幅图像都只是整个三维截面图像组的典型截面。

图3-6示意性地说明了由像素阵列形成的二维图像，其中特征选择标准由示例给出。

图7-10示意性地说明了特征跟踪是如何实现的。

图11通过示出对比媒质施加前后不同时刻拍摄的两幅乳房图像、执行所述图像的特征跟踪和配准之后对两幅图像相减得到的图像，并将这些图像与本图左侧的图2中的例子进行对比，说明了本发明的结果。图2和图11的所有图像均取自同一MR试验，并且没有进行任何除配准之外的进一步的处理。

具体实施方式

虽然本发明的一个相关方面是将Lucas&Kanade算法从二维推广到三维，下文的图例只关于二维图像以便理解起来更简单。所述方法到三维图像的推广由用数学公式表示的一般性描述中给出。

工作于三维从而工作于表示三维体积数据的截面图像组中的自动特征选择，实质上是对Lucas&Kanade的二维特征检测算法的推广，并在下面进行更详细的解释：

将很短时差内记录的同一患者的两个图像体积分别定义为I和J。

在第一图像体积I内，必须识别一定数量的已选取体素，它们将被认为两个体积间要跟踪的特征。

令I为原始/第一图像体积，J为下一体积，要求在其中寻找对应的特征。

第一步骤由在第一图像体积内计算最小特征值来构成。其如下实现：

将第一图像体积转换成三个梯度体积Ix，Iy，Iz(描述图像体积的笛卡儿坐标系统的x，y和z方向)，其中x，y和z为每个单个体素的坐标。

梯度体积是根据每个体素相对于其包含在所谓窗口内的相邻体素的强度梯度来计算的，所述窗口一般在x方向上的大小为(2x+1)、y方向上的大小为(2y+1)以及在z方向上的大小为(2z+1)，其中x，y，z＝1，2，3，4，…，n个体素。

考虑中心在目标体素的3×3×3体素阵列的邻域的大小，则所述三个梯度体积被定义为：

a)X方向上的梯度体积：

I_{Δx} (x, y, z) = \frac{I (x + 1, y, z) - I (x - 1, y, z)}{2}

b)Y方向上的梯度体积：

I_{Δy} (x, y, z) = \frac{I (x, y + 1, z) - I (x, y - 1, z)}{2}

c)Z方向上的梯度体积：

I_{Δz} (x, y, z) = \frac{I (x, y, z + 1) - I (x, y, z - 1)}{2}

作为下一步骤，对每个体素，使用如上表示的所有三个之前计算出的梯度体积，来构造梯度矩阵G。如下生成梯度矩阵：

G = [\begin{matrix} m_{200} & m_{110} & m_{101} \\ m_{110} & m_{020} & m_{011} \\ m_{101} & m_{011} & m_{002} \end{matrix}],

其中

m_{200} = I_{Δx}^{2} (x, y, z),

对于每个体素的梯度矩阵，如下计算该梯度矩阵的最小特征值λ_m：

定义：Graphics Gems IV，Paul S.Heckbert，1994，Academic Press，S193-98，

c＝m₂₀₀·m₀₂₀

{d = m}_{011}^{2}

e＝m₁₁₀·m₁₀₁ f＝m₁₀₁·m₁₀₁ p＝-m₂₀₀-m₀₂₀-m₀₀₂

q＝c+(m₂₀₀+m₀₂₀)m₀₀₂-d-e-f

r＝(e-c)m₀₀₂+dm₂₀₀-2(m₁₁₀·m₀₁₁·m₁₀₁)+f·m₀₂₀

a = q - \frac{p^{2}}{3}

b = \frac{{2 p}^{3}}{27} - \frac{pq}{3} + r

θ = \frac{1}{3} \cos^{- 1} (\frac{3 b}{an}),

其中

n = 2 \sqrt{\frac{- a}{3}}

注：由于G是个中心矩矩阵，故

梯度矩阵G的特征值计算为

λ_{1} = n \cdot \cos θ - \frac{p}{3}

λ_{2} = n [\frac{\cos θ + \sqrt{3} \sin θ}{2}] - \frac{p}{3} &DoubleRightArrow; λ_{m} = \min imum (λ_{1}, λ_{2}, λ_{3})

λ_{3} = n [\frac{\cos θ - \sqrt{3} \sin θ}{2}] - \frac{p}{3}

然后定义最小特征值λ_m的阈值。

之后，应用下述标准来考虑体素是否代表可跟踪特征：

1.如果λ_m大于阈值，则将该体素的实际位置保存在要跟踪特征的列表中。

2.如果λ_m小于所有最小特征值λ_m中最大值的某个百分比，则将其从要跟踪特征的列表中删除。

3.如果在实际特征周围的限定区域(它是可调的，如3×3×3)中存在另一特征，则将具有较小最小特征值λ_m的特征从要跟踪特征的列表中删除。

4.如果特征位置周围的三维邻域的平均信号强度值(例如1)小于或等于取值小得接近零(如10)的平均信号强度值阈值，则认为该特征位于乳房外的周围空气中，并将其从要跟踪特征的列表中去除。

5.最后，保留可限定最大数量的特征。如果要跟踪特征的列表中存在更多的特征，则从列表中删除那些具有较小最小特征值λ_m的特征。

对于二维情况，这些步骤在图3-6中进行了可视的解释。这种情况下，像素阵列用小方格阵列表示。对于与具有满足第一条标准的最小特征值的选定特征的位置对应的不同像素进行标记。在图3中，这些特征和位置通过赋予对应像素不同外观来突出。在图4中，5×5邻域是通过在这种5×5像素阵列中的内切圆来突出。

图1中给出了用于固定患者乳房的典型装置。通常，MR成像在1.0或1.5特斯拉成像仪上完成，并且使用了厂家的双乳线圈(在美国则经常使用单乳线圈设备)，而患者处于俯卧姿势。该治疗方案至少包含任意方向上的动态对比度增强T1加权序列，且采样时间在30s到180s之间(各个方案有所不同，此外还可能进行脂肪抑制)。在以0.1mmol/kg-0.2mmol/kg体重的剂量来快速浓注钆喷酸二葡甲胺盐或任何其他的顺磁性MR造影剂之前，至少进行一次测量，注射之后则要进行多次测量。

在造影剂注入肘前静脉并因而到达乳房组织之前的时刻t0，至少要获取第一MRI图像体积。而在紧接着造影剂注射后的时刻t1至少要获取另一图像体积。通常，获取的是时间延迟体积的序列(例如，在施加造影剂之后在某行对同一体积采集次数多达六次)。

为了简化当前图表，从截面图像的三维体积中取出的两幅图像具有一定的(切片)厚度。因此，讨论体素阵列而非像素更为准确。每个体素具有与对应信号强度相关的预定灰度级。在图3-10所示例子中，网格代表体素阵列。网格的每个单元就是形成图像的体素。

以图3为例，其中区分了作为标志的四种体素，它们按照对应梯度矩阵的最小特征值进行排序。用“102”表示的体素是具有最大特征值的体素，用“202”表示的体素是具有较大特征值的体素，用“302”表示的体素是具有较小特征值的体素，而用“402”表示的体素是具有最小特征值的体素。

在图4中，对于每个体素，定义了以选定为特征的像素为中心、直径为5像素/体素的欧氏距离的邻域。所述邻域的选取在图4中用圆圈表示出来。

在图4中，第一种情况用502’表示的较粗圆圈突显出来。该处有两个选定为特征的像素位于上述用像素的欧氏距离所定义的同一邻域内。这两个像素用202’和302’表示，按照上述定义，选定为特征的像素202’的梯度矩阵的特征值大于第二个像素302’。因此，该第二个像素302’就被去除。

图5给出了进一层的步骤，图3中选定为特征的像素302’不再出现。此处在图5中，像素202”设有表示其邻域的圆圈。此处也出现了与图4中相似的情况，即在像素202”的邻域内出现像素302”，该像素302”的梯度矩阵的特征值小于像素202”的梯度矩阵特征值。

因此，在这种情况下，现在将选定为特征的像素302”从要跟踪特征的列表中去除。

最后，在图6的像素阵列中，选定为特征的原始像素只有一些得到保留。现在在单个邻域内再也找不到一个以上的特征。因此，可跟踪特征在整个图像/体积上得到扩散。

根据本发明的进一层的步骤，像素的强度要进一步考察，并将那些强度小于预定值且选定为特征的像素也从对应于要跟踪特征的像素的列表中去除。近似为零的很低的阈值被用于表示周围空气中发现的特征。这个步骤没有在图中示出，但它对于普通技术人员是显而易见的。

一旦确定了被认为是二维或三维图像的特征的像素或体素，随后就得进行特征跟踪。根据本发明的优选实施例，这是通过使用所谓的“Lucas & Kanade特征跟踪器的金字塔式实现”来进行的。Jean YvesBouguet在其由Intel公司微处理器研究实验室出版的文章“PyramidalImplementation of the Lucas Kanade Feature Tracker，Description of theAlgorithm”中给出了详细的描述。

在上述出版物中公开了其二维例子的理论在下文中经过进一步改动从而适用于三维情况。下面的步骤必须对表示要跟踪特征的每个体素进行，各体素按照上述方法和算法进行了区分。

定义u为体积I中的一个点，v为体积J中的对应点。在第一步，必须由体积I和J构造两个金字塔，且I^L和J^L代表第L(0，…，m)级的体积，塔基(L＝0)表示原始体积。每个下一层体积的大小通过将直接相邻的一定数量的像素合并成一个均值而缩减。其余层的数目由体积的大小来决定。

下一步，在最高层初始化所谓的全局金字塔运动矢量g：

g^{L_{m}} = {[\begin{matrix} g_{x}^{L_{m}} & g_{y}^{L_{m}} & g_{z}^{L_{m}} \end{matrix}]}^{T} = {[\begin{matrix} 0 & 0 & 0 \end{matrix}]}^{T}

随后对所有的级L计算最终的运动矢量d^L。

这是按照如下步骤实现的：

a)在每个体积I^L内，点u^L的位置计算为：

u^{L} = {[p_{x}, p_{y}, p_{z}]}^{T} = (\begin{matrix} u \\ 2^{L} \end{matrix}),

式中p为实际的体积位置

b)在体积的金字塔表示的每级L，沿笛卡尔坐标x、y、z的每个方向计算梯度。

对于笛卡尔坐标x，按照下述等式计算体积梯度：

I_{Δx} (u_{x}, u_{y}, u_{z}) = \frac{I^{L} (u_{x} + 1, u_{y}, u_{z}) - I^{L} (u_{x} - 1, u_{y}, u_{z})}{2}

该等式相应地也适用于其他两个笛卡尔坐标y和z：

I_{Δy} (u_{x}, u_{y}, u_{z}) = \frac{I^{L} (u_{x}, u_{y} + 1, u_{z}) - I^{L} (u_{x}, u_{y} - 1, u_{z})}{2}

I_{Δz} (u_{x}, u_{y}, u_{z}) = \frac{I^{L} (u_{x}, u_{y}, u_{z} + 1) - I^{L} (u_{x}, u_{y}, u_{z} - 1)}{2}

c)随后用上述体积梯度按照下式计算梯度矩阵：

G = Σ_{x = p_{x} - ω}^{p_{x} + ω} Σ_{y = p_{y} - ω}^{p_{y} + ω} Σ_{z = p_{z} - ω}^{p_{z} + ω} [\begin{matrix} m_{200} & m_{110} & m_{101} \\ m_{110} & m_{020} & m_{011} \\ m_{101} & m_{011} & m_{002} \end{matrix}]

此处，值ω定义了影响体素的区域尺寸或邻域。

矩阵元素以类似于按照前文定义得到的矩阵元素的方式来获得，即：

m_{200} = I_{Δx}^{2} (u_{x}, u_{y}, u_{z});

m_{020} = I_{Δy}^{2} (u_{x}, u_{y}, u_{z});

m_{002} = I_{Δz}^{2} (u_{x}, u_{y}, u_{z}),

m₁₁₀＝I_Δx(u_x，u_y，u_z)·I_Δy(u_x，u_y，u_z)

m₁₀₁＝I_Δx(u_x，u_y，u_z)·I_Δz(u_x，u_y，u_z)

m₀₀₁＝I_Δy(u_x，u_y，u_z)·I_Δy(u_x，u_y，u_z)

d)对于每个级L，初始化迭代矢量

{\overset{&RightArrow;}{v}}^{0} = {[\begin{matrix} 0 & 0 & 0 \end{matrix}]}^{T}

e)其后对于k取从1到最大计数K或者直到

的最小偏移量，计算如下的体积差：

{δI}_{k} (u_{x}, u_{y}, u_{z}) = I^{L} (u_{x}, u_{y}, u_{z}) - J^{L} (u_{x} + g_{x}^{L} + v_{x}^{k - 1}, u_{y} + g_{y}^{L} + v_{y}^{k - 1}, u_{z} + g_{z}^{L} + v_{z}^{k - 1})

根据下式计算失配矢量：

{\overset{&RightArrow;}{b}}_{k} = Σ_{x = u_{x} - ω}^{u_{x} + ω} Σ_{y = u_{y} - ω}^{u_{y} + ω} Σ_{z = u_{z} - ω}^{u_{z} + ω} δ I_{k} (u_{x}, u_{y}, u_{z}) \cdot (\begin{matrix} I_{x} (u_{x}, u_{y}, u_{z}) \\ I_{y} (u_{x}, u_{y}, u_{z}) \\ I_{z} (u_{x}, u_{y}, u_{z}) \end{matrix})

之后如下定义光流矢量

{\overset{&RightArrow;}{η}}^{k} = G^{- 1} {\overset{&RightArrow;}{b}}_{k}

借助上面定义的光流矢量可将迭代运动矢量计算为：

{\overset{&RightArrow;}{v}}^{k} = {\overset{&RightArrow;}{v}}^{k - 1} + {\overset{&RightArrow;}{η}}^{k}

并且使其等于级L的最终光流：

d^{L} = {\overset{&RightArrow;}{v}}^{k}

上述步骤对每个级重复进行，直到最后一级L＝0。

这样对于下一级L-1，全局光流由上述计算的值来计算，为：

g^{L - 1} = [\begin{matrix} g_{x}^{L - 1} & g_{y}^{L - 1} & g_{z}^{L - 1} \end{matrix}] = 2 (g^{L} + d^{L})

f)结果，最终的光流矢量为

d＝g⁰+d⁰

g)现在，与体积I中的点U对应的体积J中的点v的坐标可计算为：

v＝u+d

图7-10试图用图形示意性地说明上述跟踪方法的效果，当然，这只是真实情况的简化表示。

像在图3-6中的例子一样，在图7-10中不同时刻T0和T1得到的两幅图像I和J被表示成二维像素阵列。

对应于要跟踪特征并按上面公开的方法选取的像素20在图像I的二维像素阵列中指出。成像物体(此例中为乳房)不发生运动时，图像J中相同位置的像素20”将对应于像素20。像素20”的某像素邻域用灰色阴影方框表示，对应于中心在像素20”的7×7像素阵列，并用数字10表示。

图像I的特征20的新位置用20’指示。这意味着像素20的位置处的乳房组织随着时间推移移动到像素20’的位置。当然，图中给出的所述像素20’仅仅出于解释跟踪的实现方式的目的。实际情况中该像素是未知的，否则就不必进行跟踪了。

前面用数学公式描述的、跟踪对应于选定特征的像素的一般方法在下文中通过图中的实例加以描述，并且为简便起见限于二维情况。不过，记住跟踪的一般数学描述的普通技术人员将能理解和体会三维情况下跟踪发生的方式。

跟踪开始时，假定图像I中的像素20位置不变，即具有和图像J中的像素20”相同的位置。该像素20”在下文中就称为跟踪像素，而选定为特征的图像I的像素20就称为特征像素。应用本算法时，首先计算对应于跟踪像素20”的邻域、在本例情况下即中心在图像J中的像素20”处的所述7×7像素阵列(该值可调)的某区域的逆梯度矩阵。这个邻域对应于所谓的跟踪区域。通过应用已改编从而适用于三维图像体积并将稍后详加解释的Lucas & Kanade的特征跟踪算法的金字塔式实现，就可找到跟踪特征在图像J中的新位置。

图像J中跟踪特征在20’的新位置决定了光流矢量。这在第一个迭代步骤中定义和计算。比较图7和8可见，跟踪区域以使其中心更靠近相同组织区的新位置20’的方式发生了位移。

通过重复计算表示图像I中特征的像素20与图像J中跟踪区域中心的新识别位置之间的失配矢量，进一步迭代地进行跟踪。这就导致新的光流矢量定义以及所述像素邻域或跟踪区域的新位移，如图9所不。

确定图像I中特征像素20和跟踪区域的中心之间的失配矢量的步骤一直执行，直到失配矢量达到某个可调节的最小值。

这种情形如图10所示。参见图7-9可知，像素20’的位置保持不变以便突出这一事实，即移动跟踪区域，使其最终以新正确位置20’为中心，该新正确位置20’和图像I中的像素20表示相同的乳房部位。

上述例子仅仅用了由一个像素代表的一个特征来实现。上述步骤必须针对图像中的每个选定特征进行。

跟踪所有的特征到其新位置后，整个图像/体积中、不过仅限于代表乳房组织的区域内的不同位置处存在一定量的光流矢量。所述算法相对于其他配准算法的最大进步在于，它将光流矢量遍布图像/体积地展开。

其后，通过对图像/体积J的每个特征及其周围邻域施加逆光流矢量，来实现变形(morphing)。这样得到新的虚拟图像/图像体积J’。

图11展示了所述方法处理实际诊断图像的效果。左侧第一幅图像反映的是未施加造影剂在时间码t0时进行的乳房试验的典型切片。左侧第二幅图像是在施加造影剂后的时间码t1时且在相同的切片位置处获取的切片。左侧第三幅图像显示了从时间码t1时间码t0时的每个像素的信号强度的逐像素相减所得到的结果，即所谓的减影图像。在时间间隔t1-t0内，如果不存在乳房的相对移动，所得图像在没有造影剂的地方将包含较宽部分的信号强度为零且呈黑色的像素。时间间隔t1-t0内的运动效果是，由于图像/体积I和J的相同位置的像素不对应于相同的成像组织区域，而产生伪影。在本例中，减影图像显示有个乳房包含对比度增强肿瘤，如小箭头所示。大多数类型的恶性和良性乳房肿瘤均表现出由血管增生造成的附加造影剂摄取量。图中双乳顶部的另外白色结构属于运动伪影。双乳中央的其余白色结构是由正常乳房软组织的微少造影剂摄取量和运动结合所造成的。

图11的右侧展示了应用所述配准方法得到的相应图像。时间码t0时的图像没有被本算法改变，因而仍保持原样。但是经过三维特征跟踪后，计算出了新的虚拟图像体积。右侧时间码t1’时的第二幅图像展示了实行配准后相同切片位置处的对应虚拟图像。右侧最后一幅图像显示了时间码t1和时间码t0时的每个像素的信号强度的逐像素相减所得到的结果，即新的(虚拟的)减影图像。仍然有个乳房示出了对比度增强肿瘤，如小箭头所示。乳房内代表正常乳房软组织的区域仍然反映出适度的造影剂摄取量，并且由于运动伪影的减少而不那么亮。由于运动伪影补偿，双乳顶部较宽部分的白色结构消失了。

图11的例子是一种简化的情况，这是因为仅仅考虑了两个图像的序列，且它们分别在施加造影剂前和后获取。

实际上，可以在造影剂注射前后每隔固定时间间隔来采集一系列截面图像(例如，超过50幅)，并可在每对连续图像之间计算减影图像。

除了在一个系列内采集一幅以上的截面图像外，也可一个系列接一个系列地采集两个系列以上(如5个甚至更多)的图像。这就得到一些减影图像/体积。当然，本公开算法可以应用于所有的系列组合，而不管这些系列是在造影剂施加前或施加后获取的，目的是减小其间发生的任何运动。

在下面的附录中，给出了Lucas & Kanade跟踪算法的金字塔式实现的详细理论。

依照本发明的上述方法可以通过被编码成能由计算机运行的程序的算法来实现。

该程序可以保存在例如磁性、光磁或光学基片的便携式或静态存储设备上，如一张或多张软盘、一张或多张磁光盘、一张或多张CD-ROM或DVD-ROM，或者也可以保存在移动或固定硬盘上。另外，该程序也可以保存在譬如笔式驱动或类似设备之类的固态电子设备上。

编码程序的伪代码例子如下：

//build gradient volumes in X，Y and Z direction；

buildGradientVolumes()

{

for all positions in volume do

{

X-GradientVolume@position＝

GreyValue(originalVolume@position(x-1，y，z))-

GreyValue(originalVolume@position(x+1，y，z))；

Y-GradientVolume@position＝

GreyValue(originalVolume@position(x，y-1，z))-

GreyValue(originalVolume@position(x，y+1，z))；

Z-GradientVolume@position＝

GreyValue(originalVolume@position(x，y，z-1))-

GreyValue(originalVolume@position(x，y，z+1))；

}

//////////////////////////////////////////////////////////////////////////////////////////////////

//////

//selection of features

//////

for all Voxel in First Volume do

{

//build matrix with X-，Y-and Z-gradients；

gradientMatrix＝buildGradientMatrix(position(x，y，z)，

i_SearchRegionSize)；

calcMnimumEigenvalue of gradientMatrix；

if(MnimumEigenvalue＞i_GlobalMnimumQuality)

{

addToFeatureList(x，y，z，MnEigenValue)；

}

for all Features in FeatureList from highest to lowest eigenvalue to

{

if(Feature-＞MnEigenValue＜i_percentvalue of

OverallMaximumEigenvalue){

deleteFeatureFromList()；

}

if(distance to another Feature＜i_MnDistance){

deleteFeatureFromList()；

}

if(Number of Features＞i_MaximumNumberOfFeatures)

{

deleteFromList(FeaturesWithSmallerMnEigenValue)；

}

//////

//tracking of selected features

//////

for each Feature in FeatureList do

{

//build matrix from X-，Y-and Z-gradients of the original Volume

at feature position

GradientMatrix＝buildGradientMatrix(featurePosition，

i_TrackingRegionSize)；

invertedGradientMatrix＝invert(GradientMatrix)；

do

{

if(FeatureMovementVector＞i_MaximumAllowedDisplacement)

{

featureTracked＝false；

break；

}

for all Offsets in i_TrackingRegionSize do

{

diff＝

GreyValue(OriginalVolume@OriginalFeaturePosition+Offset)

-

GreyValue(MovedVolume@TrackedFeaturePosition+Offset)；

vectorX＝vectorX+diff*

GreyValue(X-GradientVolume@actualFeaturePosition)；

vectorY＝vectorY+diff*

GreyValue(Y-GradientVolume@actualFeaturePosition)；

vectorZ＝vectorZ+diff*

GreyValue(Z-GradientVolume@actualFeaturePosition)；

}

mismatchVector＝mathVector(vectorX，vectorY，vectorZ)；

currentMovementVector＝invertedGradientMatrix*

mismatchVector；

FeatureMovementVector＝FeatureMovementVector+

currentMovementVector；

if(currentMovementVector＜i_MnimumDisplacement){

featureTracked＝true；

break；

}

if(NumberOfIterations＞i_MaximumAllowedIterations){

featureTracked＝false；

break；

}

NumberOfIterations＝NumberOfIterations+1；

}

While(true)

}

//////

//calculation of gradient matrix

function buildGradientMatrix(position，regionSize)

{

for all Offsets in regionSize do

{

Ix＝GreyValue(X-GradientVolume@positon+Offset)；

Iy＝GreyValue(Y-GradientVolume@positon+Offset)；

Iz＝GreyValue(Z-GradientVolume@positon+Offset)；

Ixx＝Ixx+Ix*Ix；

Iyy＝Iyy+Iy*Iy；

Izz＝Izz+Iz*Iz；

Ixy＝Ixy+Ix*Iy；

Ixz＝Ixz+Ix*Iz；

Iyz＝Iyz+Iy*Iz；

}

return gradientMatrix(Ixx Ixy Ixz

Ixy Iyy Iyz

Ixz Iyz Izz)；

}

Claims

1.用于减少物体运动造成的成像伪影的生物医学图像配准方法，其包括以下步骤：

a)提供同一物体的至少第一幅和第二幅数字或数字化图像，或者至少第一组和第二组该物体的截面图像，所述图像由像素或体素的二维或三维阵列形成；

b)通过选取一定数量的被设为标志或特征的像素或体素，在第一幅图像或第一组图像中定义一定数量的标志，即所谓的特征，并生成要跟踪的所述特征的列表；

c)通过对被选定为特征的每个像素或体素确定从第一幅到第二幅图像或从第一组到第二组图像的光流矢量，从第一幅到第二幅图像或从第一组到第二组图像跟踪选定为特征的每个像素或体素的位置；

d)通过对第二幅或第二组图像的像素或体素施加逆光流，来对第一幅和第二幅图像、或对第一组和第二组图像进行配准；

其特征在于

执行自动特征选择步骤，该步骤包括：

B1)在第一幅图像或第一组截面图像的每个像素或体素周围定义像素或体素邻域，所述的像素或体素邻域包括有限数量的像素或体素；对于单幅图像，选用二维邻域；对于截面图像组，则选用三维邻域；

B3)定义平均信号强度值阈值；

B4)比较在步骤B2中为每个像素或体素邻域确定的平均信号强度，并将所述平均信号强度值与平均信号强度值阈值进行比较；

B5)如果步骤B4中所述邻域的平均信号强度值大于阈值，则将所述像素或体素定义为要跟踪特征，并将其添加到要跟踪特征列表中。

2.根据权利要求1所述的方法，其特征在于，借助实验，通过在一组测试或样本图像上执行所述方法，并比较取不同阈值时获得的结果，来确定所述平均信号强度阈值。

3.根据权利要求1或2所述的方法，其特征在于，通过所谓的Lucas & Kanade自动特征跟踪方法或算法，来自动地跟踪被选作在第二幅图像中要跟踪的特征或标志的图像的像素或体素。

4.根据权利要求3所述的方法，其特征在于，将发明出来仅作用于二维图像的Lucas & Kanade算法的原始实现推广到作用于三维图像体积而非二维图像，以适用于生物医学截面图像体积，其包含如下步骤：

a)提供由三维体素阵列构成的第一三维图像体积I；

b)对于每个目标体素定义体素邻域，该邻域包围所述目标体素，并具有中心在所述目标体素的体素阵列，该体素阵列具有可调节的大小，例如3×3×3；

c)定义笛卡尔坐标系统，以及相对于所述笛卡尔坐标系统的每个轴，计算所谓的梯度体积，由第一图像体积的每个体素相对于其邻域体素的强度梯度来形成所述梯度体积，且定义为：

a)X方向的梯度体积：

I_{Δx} (x, y, z) = \frac{I (x + 1, y, z) - I (x + 1, y, z)}{2}

b)Y方向的梯度体积：

I_{Δy} (x, y, z) = \frac{I (x, y + 1, z) - I (x, y - 1, z)}{2}

c)Z方向的梯度体积：

I_{Δz} (x, y, z) = \frac{I (x, y, z + 1) - I (x, y, z - 1)}{2}

d)对于每个目标体素，计算所谓的梯度矩阵，所述梯度矩阵定义为：

G = [\begin{matrix} m_{200} & m_{110} & m_{101} \\ m_{110} & m_{020} & m_{011} \\ m_{101} & m_{011} & m_{002} \end{matrix}],

其中

\begin{matrix} m_{200} = I_{Δx}^{2} (x, y, z) \\ m_{020} = I_{Δy}^{2} (x, y, z), m_{002} = I_{Δz}^{2} (x, y, z), \\ m_{110} = I_{Δx} (x, y, z) \cdot I_{Δy} (x, y, z), \\ m_{101} = I_{Δx} (x, y, z) \cdot I_{Δz} (x, y, z), \\ m_{011} = I_{Δy} (x, y, z) \cdot I_{Δz} (x, y, z) \end{matrix}

e)对于三维图像体积的每个目标体素的每个梯度矩阵，应用下式来计算最小特征值λ_m：

定义：

c＝m₂₀₀·m₀₂₀

d = m_{011}^{2}

q＝c+(m₂₀₀+m₀₂₀)m₀₀₂-d-e-f

a = q - \frac{p^{2}}{3}

b = \frac{2 p^{3}}{27} - \frac{pq}{3} + r

θ = \frac{1}{3} \cos^{- 1} (\frac{3 b}{an}),

其中

n = 2 \sqrt{\frac{- a}{3}}

注意：由于G是中心矩矩阵，所以

| 3 \frac{b}{an} | \leq 1^a \leq 0

如下，计算所述梯度矩阵的最小特征值：

λ_{1} = n \cdot \cos θ - \frac{p}{3}

λ_{2} = n [\frac{\cos θ + \sqrt{3} \sin θ}{2}] - \frac{p}{3}

λ_m＝minimum(λ₁，λ₂，λ₃)

λ_{3} = n [\frac{\cos θ - \sqrt{3} \sin θ}{2}] - \frac{p}{3}

f)定义最小特征值λ_m的阈值；

g)将对应梯度矩阵的最小特征值λ_m满足一定标准的每个体素选择作为图像I中要跟踪的特征，这些标准如下：

i)λ_m大于所述阈值；

ii)λ_m不小于所有最小特征值λ_m中的最大值的某个百分比；

iii)如果在保持为选定特征的目标体素周围的已定义体素邻域内存在另一个被选作特征的体素，则仅将其梯度矩阵具有较大的最小特征值的一个体素选为特征，而将另一个体素从可跟踪特征列表中删除；

iv)选作特征的体素周围的三维邻域的平均信号强度值大于可调节平均信号强度值阈值；

5.根据权利要求4所述的方法，其特征在于，选作特征的体素的邻域的平均信号强度的阈值被设为大约10。

6.根据上述一项或多项权利要求所述的方法，其特征在于，对于第一图像体积I中被选作特征的每个体素，相对于所述特征在第二较晚时刻得到的相同物体的另一图像J的体素阵列中的位置，对所述特征进行跟踪，且通过所谓的Lucas & Kanade特征跟踪器的金字塔式实现来执行所述跟踪。

7.根据权利要求6所述的方法，其特征在于，所谓的Lucas &Kanade特征跟踪器的金字塔式实现包含如下步骤：

将点u定义为对应于图像I的三维体素阵列的图像体积I中的点，将v定义为对应于图像J的三维体素阵列的图像体积J中的对应点，这些点具有图像I中被选作特征的体素的坐标；

-由图像体积I和J构造两个理想金字塔，且I^L和J^L表示第L(0，…，m)级体积，塔基(L＝0)表示原始体积；

-每个下一层体积的大小通过将在直接邻域中的一定数量的像素合并成一个均值而缩减；

-定义所谓的全局金字塔运动矢量g，并在最高层L_m上对其进行初始化：

g^{L_{m}} = {[\begin{matrix} g_{x}^{L_{m}} & g_{y}^{L_{m}} & g_{z}^{L_{m}} \end{matrix}]}^{T} = {[\begin{matrix} 0 & 0 & 0 \end{matrix}]}^{T}

i)将点u的位置定义为u^L，并通过以下等式来计算所述位置，其中L具有金字塔的实际层或级的值：

u^{L} = {[p_{x}, p_{y}, p_{z}]}^{T} = (\begin{matrix} u \\ 2^{L} \end{matrix}),

式中p为实际的体积位置

ii)根据体积梯度的以下等式，计算笛卡尔坐标x、y和z的每个方向上的体积梯度：

I_{Δx} (u_{x}, u_{y}, u_{z}) = \frac{I^{L} (u_{x} + 1, u_{y}, u_{z}) - I^{L} (u_{x} - 1, u_{y}, u_{z})}{2}

I_{Δy} (u_{x}, u_{y}, u_{z}) = \frac{I^{L} (u_{x}, u_{y} + 1, u_{z}) - I^{L} (u_{x}, u_{y} - 1, u_{z})}{2}

I_{Δz} (u_{x}, u_{y}, u_{z}) = \frac{I^{L} (u_{x}, u_{y}, u_{z} + 1) - I^{L} (u_{x}, u_{y}, u_{z} - 1)}{2}

iii)使用梯度体积，按照下式来计算梯度矩阵：

G = Σ_{x = p_{x} - ω}^{p_{x} + ω} Σ_{y = p_{y} - ω}^{p_{y} + ω} Σ_{z = p_{z} - ω}^{p_{z} + ω} [\begin{matrix} m_{200} & m_{110} & m_{101} \\ m_{110} & m_{020} & m_{011} \\ m_{101} & m_{011} & m_{002} \end{matrix}]

其中

m_{200} = I_{Δx}^{2} (u_{x}, u_{y}, u_{z});

m_{020} = I_{Δy}^{2} (u_{x}, u_{y}, u_{z});

m_{002} = I_{Δz}^{2} (u_{x}, u_{y}, u_{z}),

m₁₁₀＝I_Δx(u_x，u_y，u_z)·I_Δy(u_x，u_y，u_z)

m₁₀₁＝I_Δx(u_x，u_y，u_z)·I_Δz(u_x，u_y，u_z)

m₀₀₁＝I_Δy(u_x，u_y，u_z)·I_Δy(u_x，u_y，u_z)

此处，值ω定义了影响表示跟踪特征的目标体素的区域尺寸或体素邻域。

iv)对于每个级L，初始化如下定义的迭代矢量

{\overset{&RightArrow;}{v}}^{0} = {[\begin{matrix} 0 & 0 & 0 \end{matrix}]}^{T}

v)对于k取从1到最大计数K或者直到的最小偏移，计算如下的体积差：

δ I_{k} (u_{x}, u_{y}, u_{z}) = I^{L} (u_{x}, u_{y}, u_{z}) - J^{L} (u_{x} + g_{x}^{L} + v_{x}^{k - 1}, u_{y} + g_{y}^{L} + v_{y}^{k - 1}, u_{z} + g_{z}^{L} + v_{z}^{k - 1})

vi)根据下式计算失配矢量：

{\overset{&RightArrow;}{b}}_{k} = Σ_{x = u_{x} - ω}^{u_{x} + ω} Σ_{y = u_{y} - ω}^{u_{y} + ω} Σ_{z = u_{z} - ω}^{u_{z} + ω} δ I_{k} (u_{x}, u_{y}, u_{z}) \cdot (\begin{matrix} I_{x} (u_{x}, u_{y}, u_{z}) \\ I_{y} (u_{x}, u_{y}, u_{z}) \\ I_{z} (u_{x}, u_{y}, u_{z}) \end{matrix})

vii)按照下式确定光流矢量

\overset{&RightArrow; k}{η} = G^{- 1} {\overset{&RightArrow;}{b}}_{k}

式中，G^-1为在步骤iii)中确定的梯度矩阵G的逆矩阵

viii)使用上面定义的光流矢量如下计算迭代运动矢量：

\overset{&RightArrow; k}{v} = \overset{&RightArrow; k - 1}{v} + \overset{&RightArrow; k}{η}

并且使其等于级L的最终光流：

d^{L} = \overset{&RightArrow; k}{v}

ix)对各级重复步骤i)-viii)，直到最后级L＝0，得到下式定义的最终的光流矢量：

d＝g⁰+d⁰

x)通过下式，确定与代表体积I中要跟踪特征的点u对应的体积J中的点v的坐标：

v＝u+d

xi)对于对应于图像I中选定特征的每个体素，重复执行上述方法。

8.根据权利要求7所述的方法，其特征在于，两幅图像I和J的配准通过如下方式实现：对图像J中识别为与表示某特征的体素相对应的每个点施加逆光流矢量，其中该特征与图像I中对应于某选定特征的体素的点U相对应；或者对与图像I中识别为选定特征的体素对应的每个点u施加光流矢量。

9.根据前面的一项或多项权利要求所述的方法，其特征在于，它与用于对比媒质增强诊断成像，尤其是对比媒质增强MRI或超声成像的方法结合提供，该方法包含下列步骤：

a)提供由MRI或超声获取的同一物体的至少第一幅和第二幅数字或数字化图像，或者至少第一组和第二组该物体的截面图像，所述图像由像素或体素的二维或三维阵列形成，其中对于相同组织或组织区或者相同的解剖区的扫描，是在第二时刻或者任何随后时刻、且在所述组织或组织区或者所述解剖区中存在对比媒质时进行的；

b)在第一幅图像或第一组图像中，通过选取一定数量的像素或体素，将该像素或体素设置作为标记或特征，来定义一定数量的标志，即所谓的特征，并生成要跟踪的所述特征的列表；

c)通过对选定为特征的每个像素或体素确定从第一幅到第二幅图像或从第一组到第二组图像的光流矢量，来从第一幅到第二幅图像或从第一组到第二组图像跟踪选定为特征的每个像素或体素的位置；

d)通过对第二幅或第二组图像的像素或体素施加逆光流，来对第一幅和第二幅图像或第一组和第二组图像进行配准。

10.根据权利要求1-9所述的方法，其特征在于，在所述特征选择步骤之后以及执行所述特征跟踪步骤之前，进行进一步的自动特征选择步骤，该步骤包括：

B2)对于每个像素或体素，确定所述像素或体素的邻域的所有像素或体素的平均信号强度值；

B3)定义平均信号强度值阈值；

B4)比较在步骤B2中确定的每个像素或体素邻域的平均信号强度值，并将所述平均信号强度值与所述平均信号强度值阈值进行比较；

11.用于配准两个数字或数字化图像/图像体积的计算机程序，在该计算机程序中对根据前面的一项或多项权利要求的所述方法进行编码，该计算机程序的特征在于以下伪代码：

//build gradient volumes in X，Y and Z direction；

buildGradientVolumes()

{

for all positions in volume do

{

X-GradientVolume@position＝

GreyValue(originalVolume@position(x-1，y，z))-

GreyValue(originalVolume@position(x+1，y，z))；

Y-GradientVolume@position＝

GreyValue(originalVolume@position(x，y-1，z))-

GreyValue(originalVolume@position(x，y+1，z))；

Z-GradientVolume@position＝

GreyValue(originalVolume@position(x，y，z-1))-

GreyValue(originalVolume@position(x，y，z+1))；

}

//////

// selection of features

//////

for all Voxel in First Volume do

{

//build matrix with X-，Y-and Z-gradients；

gradientMatrix ＝ buildGradientMatrix(position(x，y，z)，

i_SearchRegionSize)；

calcMnimumEigenvalue of gradientMatrix；

if(MnimumEigenvalue＞i_GlobalMnimumQuality)

{

addToFeatureList(x，y，z，MnEigenValue)；

}

for all Features in FeatureList from highest to lowest eigenvalue to

{

if (Feature->MnEigenValue ＜ i_percentvalue of

OverallMaximumEigenvalue){

deleteFeatureFromList()；

}

if(distance to another Feature＜i_MnDistance){

deleteFeatureFromList()；

}

if(Number of Features＞i_MaximumNumberOfFeatures)

{

deleteFromList(FeaturesWithSmallerMnEigenValue)；

}

//////

//tracking of selected features

//////

for each Feature in FeatureList do

{

//build matrix from X-，Y-and Z-gradients of the original Volume

at feature position

GradientMatrix ＝ buildGradientMatrix(featurePosition，

i_TrackingRegionSize)；

invertedGradientMatrix＝invert(GradientMatrix)；

do

{

if(FeatureMovementVector＞i_MaximumAllowedDisplacement)

{

featureTracked＝false；

break；

}

for all Offsets in i_TrackingRegionSize do

{

diff＝

GreyValue(OriginalVolume@OriginalFeaturePosition+Offset)

-

GreyValue(MovedVolume@TrackedFeaturePosition+Offset)；

vectorX ＝ vectorX + diff *

GreyValue(x-GradientVolume@actualFeaturePosition)；

vectorY ＝ vectorY + diff *

GreyValue(Y-GradientVolume@actualFeaturePosition)；

vectorZ ＝ vectorZ + diff *

GreyValue(Z-GradientVolume@actualFeaturePosition)；

}

mismatchVector＝mathVector(vectorX，vectorY，vectorZ)；

currentMovementVector ＝ invertedGradientMatrix *

mismatchVector；

FeatureMovementVector ＝ FeatureMovementVector +

currentMovementVector；

if(currentMovementVector＜i_MnimumDisplacement){

featureTracked＝true；

break；

}

if(NumberOflterations＞i_MaximumAllowedIterations){

featureTracked＝false；

break；

}

NumberOflterations＝NumberOflterations+1；

}

While(true)

}

//////

//calculation of gradient matrix

function buildGradientMatrix(position，regionSize)

{

for all Offsets in regionSize do

{

Ix＝GreyValue(X-GradientVolume@positon+Offset)；

Iy＝GreyValue(Y-GradientVolume@positon+Offset)；

Iz＝GreyValue(Z-GradientVolume@positon+Offset)；

Ixx＝Ixx+Ix*Ix；

Iyy＝Iyy+Iy*Iy；

Izz＝Izz+Iz*Iz；

Ixy＝Ixy+Ix*Iy；

Ixz＝Ixz+Ix*Iz；

Iyz＝Iyz+Iy*Iz；

}

return gradientMatrix(Ixx Ixy Ixz

Ixy Iyy Iyz

Ixz Iyz Izz)；

}