CN1781123A - 用于跟踪运动对象的整体形状的系统及方法 - Google Patents

Abstract

本发明公开了一种用于跟踪运动对象的整体形状的系统和方法。沿着整体形状的初始轮廓定义一个或多个控制点。当对象运动时，对一个或多个控制点中的每一个进行跟踪。使用多种技术来表示运动中的控制点的位置的不确定性。使用先验的形状模型，利用该不确定性来限制整体形状。在一种替代的实施例中，针对每个控制点建立多个外观模型，并结合每个模型所产生的运动向量以便跟踪对象的形状。可使用显示器和彩色向量来视觉跟踪对象形状的运动。

Description

用于跟踪运动对象的整体形状的系统及方法

相关申请的交叉参考

本申请要求2003年3月7日提交的序列号为60/452669的美国临时申请、2003年5月27日提交的序列号为60/473425的美国临时申请以及2003年10月3日提交的序列号为60/508367的美国临时申请的优先权，它们的内容通过引用而全部包含于本申请中。

技术领域

本发明提供了一种对运动对象中的形状进行跟踪的方法，更特别地，提供了一种在异方差噪声存在时对具有线性约束的形状进行跟踪的方法。

背景技术

对于绝大多数视觉跟踪应用来说，测量数据是不确定的，并且有时会被遗漏：获取的图像带有噪声并且失真，而掩蔽可导致感兴趣的对象的一部分不能被观察。不确定性可以是整体均一的；但在绝大多数现实世界的情况中，其本质上是异方差的，也就是说，既是各向异性的又是不均一的。一个很好的例子是超声心动图(超声心脏数据)。超声倾向于反射人造制品，例如镜面反射器，比如那些由膜制成的镜面反射器。由于单一“视图方向”，镜面结构的垂直表面产生强回声，但倾斜或“离轴”的表面可能产生弱回声，甚至根本不产生回声(声波“失落(drop out)”)。对于超声心动图来说，在组织表面与超声波束平行的情况下，在心脏区域可能发生这种失落。

由于其有效性、相对低的成本以及非侵入性，心脏超声图像被广泛用于心脏功能的评价。特别地，心室运动的分析是一种用于评价局部缺血以及梗塞的有效方法。对心内膜壁进行分段或检测是量化左心室的弹性及收缩性的首要步骤。一些现有方法的例子包括基于像素的分段/聚类方法(例如Color Kinesis)，光流变化，可变形模板及Markov随机过程/场以及活动轮廓/蛇形方法。一些方法被用于2维、3维或4维(3D+时间)空间。

然而，大多数现有的分段或检测方法并不尝试再现心内膜壁的精确的局部运动，而且在绝大多数情况下，沿着该壁的运动成分被忽略。这种简化处理也由仅仅沿着当前轮廓的法线进行搜寻的轮廓跟踪装置所采用。这对于局部壁异常的检测来说是不适用的，因为异常的左心室的局部运动可能离开轮廓的法线，更不必提及整体运动，例如平移或旋转(由于声谱技师手的运动或者患者的呼吸运动)也引起轮廓上的偏离法线的局部运动。因此需要跟踪心内膜壁的整体形状以及其局部运动，以检测局部壁运动的异常。这种信息可被用于局部缺血及梗塞的进一步诊断。

通常，可将协方差分配给反映潜在的异方差噪声的图像特征或者流估计。当数据被净化成具有低的总噪声水平时，异方差的特性可被忽略，并且整体不确定性可代替局部估计。然而，对于有很大噪声的输入来说，特别是对于那些具有空间变化的结构噪声来说，以局部协方差矩阵形式被编码的信息就成为确保对象或者潜在的图像结构的可靠及鲁棒的推断的关键所在。

通常习惯于在跟踪框架中施加模型限制。这样的例子包括简单模型，诸如斑点或参数化的椭圆；以及复杂模型，诸如识别模板。在大多数实际情况下，由于捕获主要形状变化的模式的数量是有限的并且通常比用于描述形状的特征成分的初始数量要少得多，所以子空间模型适于形状跟踪。此外，一种基于主成分分析(PCA)的特征形状子空间可捕获任何复杂的形状变化，这些变化在初始空间中即使利用非常简单的参数模型也是高度非线性的。

如果测量向量受到异方差噪声的影响，则到限制性子空间中的正交投影不仅是未对齐的，而且在信息丢失方面也会造成很大损害。仅仅在噪声既是各向同性的又是均一的特定情况下才能对齐。

然而，大多数现有的关于子空间受限制的跟踪的工作没有考虑测量中的异方差噪声。在“点分布模型”或者“活动形状模型”中，基于具有地标点一致性的训练形状，导出基于PCA的子空间形状模型。所产生的具有特征形状的子空间捕获训练数据集中最有意义的变化。在检测时间，使模型扰动以产生与在候选位置的测试图像匹配的合成图像。然而，在该过程中没有对测量噪声进行建模。

即使当异方差噪声特性是可利用的时候，它们通常在拟合子空间模型过程中也被忽视。例如，在一种已知的在测量时捕获完整的协方差矩阵的方法中，应用相当特别的阈值，以使测量平均值被限定于由模型协方差定义的超椭圆体限制。这种操作不依赖于测量噪声。

另一种已知的方法应用两步法在Kalman滤波框架中施加形状空间限制。这种形状空间是线性变换的仿射子空间或者特征子空间。然而，到形状空间中的投影是正交的，而没有考虑量的异方差噪声。因此，这种方法导致在投影期间信息的丢失。

另一种已知的方法是将Gaussian分布用于对感兴趣的对象对象的外观(在其情况下为脸)进行适应性建模，其中使用EM算法获悉该外观。如在本发明中，在协方差矩阵中捕获局部不确定性。不同之处在于，本发明专门研究了当各向异性不确定性存在时子空间模型的限制以及在投影上相交的临界选择。

另一种已知的方法隐含地在光流估计过程中使用子空间限制，其也利用流的不确定性。尽管在不同的应用中使用了不同的框架，本发明认识到“更可靠的流向量将在子空间投影过程中具有更大的影响”。

另一种已知的方法是将异方差回归用于拟合椭圆以及基本矩阵。在具有参数化模型的初始空间中实现拟合。在本发明中，避免了形状变化的参数化——其可以是非常复杂以及高度非线性。替代地，本发明通过例如PCA建立子空间线性随机模型，并且获得关于拟合数据的平均值和协方差的封闭形式的解。

基于M估计器或者RANSAC的鲁棒模型匹配也可被用于限制或者消除相对于模型为界外值的数据成分的影响。另外，在这些框架中没有对(在空间或时间上)局部变化的不确定性进行研究。

其它相关的方法包括数据归因，“填充”具有似是而非的值的丢失数据的实践。在该领域中的工作植根于在语音识别、医学图像分析以及社会科学等方面具有广泛应用的统计学。然而，数据归因问题的公式化通常假定0-1的可获得性，也就是说，数据成分要么丢失，要么可获得。需要提供一种统一的框架，以便将子空间模型限制与关于测量噪声的形状动态和异方差特性的信息以及关于形状动态的信息相融合。

发明内容

本发明提供一种用于跟踪运动对象的整体形状的系统和方法。沿着整体形状的初始轮廓定义一个或多个控制点。当对象运动时，对一个或多个控制点中的每一个进行跟踪。估计运动时控制点的位置的不确定性。表示不确定性的一种形式是协方差矩阵。当采用子空间形状限制模型时，利用非正交投影和/或信息融合来研究不确定性。显示每个作为结果发生的轮廓。

本发明还提供一种视觉跟踪对象形状的运动的系统。产生一个或多个第一彩色向量来表示控制点沿着形状的轮廓的收缩。产生一个或多个第二控制向量来表示控制点沿着形状的轮廓的扩张。第一及第二彩色向量被周期地显示，从而显示形状的运动。

本发明还提供一种用于跟踪运动对象的整体形状的方法。沿着整体形状定义一个或多个控制点。当对象运动时，对一个或多个控制点中的每一个进行跟踪。针对每个控制点建立多个外观模型。结合由每个模型所产生的运动向量以便跟踪形状。

附图说明

下面将参照附图对本发明的优选实施例进行更为详细的描述，其中相同的附图标记表示相同的元件。

图1是根据本发明的用于实施形状跟踪方法的示例系统的框图；

图2是一种心脏超声波心动描记图像，说明了声波失落区域和所估计的局部壁运动的不确定性；

图3是左心室的超声波心动描记图像，说明了具有控制点定位不确定性的心内膜轮廓，这些控制点具有强烈的各向异性及不均一性；

图4说明了具有不同α值的增量式PCA模型和强适应性PCA(SA-PCA，strongly-adapted-PCA)模型的实例；

图5a-5c说明了左心室的超声波心动描记图像，其中心内膜壁被初始化(a)，并利用IPCA(b)及SA-PCA(c)来跟踪；

图6a及6b说明了左心室的超声波心动描记图像，其中根据本发明对心内膜壁的运动进行跟踪；

图7a-7d说明了代表根据本发明的心尖四腔图中左心室心内膜壁的超声波心动描记图像的测试轮廓；

图8a-8d说明了代表根据本发明的胸骨旁短轴图中左心室心内膜壁的超声波心动描记图像的测试轮廓；

图9说明了根据本发明作为用于显示图像的方法的例子的超声波心动描记图像；以及

图10说明了根据本发明的一个实施例的基于多模型控制点的跟踪装置。

具体实施方式

本发明提供了一种在存在异方差噪声时对具有线性限制的形状进行跟踪的方法。将使用这种方法的一个例子是跟踪心肌壁的整体形状以及其局部运动，以检测心脏内局部壁运动的异常。该方法还可用于跟踪心脏的心内膜壁和心外膜壁。本领域技术人员应当理解，本发明可用于其它应用，在这些应用中形状跟踪是有用的，例如但不限于识别人类特征的运动、例如头部运动、面部特征、手部运动或其他身体运动。本发明还可用在2维、3维及4维(3D+时间)解剖结构、例如随着时间发展的心脏、肺或者肿瘤等的医学分析中。

为了描述本发明，将以对左心室心内膜壁的跟踪为例进行描述。图1说明了超声波心动描记系统的示例结构，该超声波心动描记系统使用根据本发明的用于对左心室心内膜壁的形状进行跟踪的方法。利用医学传感器102、例如超声换能器对患者进行检查。传感器102被用于获得与特定医学检查相一致的医学测量结果。例如，有心脏问题的患者可做超声心动图检查，以帮助诊断特定的心脏疾病。超声系统根据各种透视图提供心脏的二维、三维及四维(3D+时间)图像。

传感器102所获取的信息被传送给处理器104，其可以是工作站或者个人计算机。处理器104将传感器数据转换成图像并将其传送给显示器108。显示器108还可传送其它与该图像有关的图形信息或者信息表。根据本发明，还可为处理器104提供代表心内膜壁的初始轮廓的数据。该数据可由用户、例如医生或声谱技师手动地提供，或由处理器104自动提供。该轮廓包括一系列单独的点，这些点的运动由处理器104进行跟踪并在显示器108上说明。关于如何对这些单独的点进行跟踪的细节将在下面更详细地进行描述。

除了来自医学传感器102的数据外，处理器104还可接收其它数据输入。例如，处理器可接收来自与处理器104相关联的数据库106的数据。这种数据可包括代表心内膜壁的潜在轮廓形状的子空间模型。这些子空间模型可以是代表多个患者的左心室的图像，或者可以是由计算机基于统计信息生成的轮廓形状模型。处理器104利用已知的方法、例如Bayesian核匹配方法或者基于光流的方法对轮廓形状的各个点进行跟踪。通过使用多模板适应性匹配框架来校正跟踪过程中的误差积累。跟踪的不确定性在每个点处以协方差矩阵的形式来表示，协方差矩阵随后完全由利用非正交投影的子空间形状限制所利用。

图2显示了心脏的典型的超声波心动描记图像。左心室心内膜壁的具有声波失落的那部分由实心椭圆208标出。局部壁运动的估计由虚线椭圆202及204所描绘。由于声波失落，心内膜壁并不总是位于图像中最强的边缘处。根据本发明，使用统一的框架来融合带有关于形状动态和测量噪声的异方差特性的信息以及关于形状动态的信息的子空间模型限制。子空间模型可以采取特定子空间分布的形式，例如Gaussian，或者简单的子空间限制，例如特征空间模型。

本发明跟踪代表心内膜壁的轮廓上的各个控制点。可利用Bayesian核匹配方法或者基于光流的方法来进行跟踪。在由作为共同发明人的Dorin Comaniciu在Annual Conf.of the German Societyfor Pattern Recognition(DAGM′02)(苏黎世，瑞士，438-445，2002年)上发表的题为“Bayesian Kernel Tracking”一文中描述了Bayesian核匹配方法的一个例子，通过引用而将该文的全部包含在内。在序列号为10/681,702的题为“Density Estimation-BasedInformation Fusion for Multiple Motion Computation”的共同未决的申请中描述了基于光流的跟踪各个点的方法一个例子，在此通过引用而将该申请的全部包含在内。根据本发明，使用多模板适应性匹配框架并利用心脏运动的周期特性来对跟踪过程中的误差积累进行校正。跟踪的不确定性在每个点处以协方差矩阵的形式来表示，该协方差矩阵随后完全由利用非正交投影的子空间形状限制所利用。

跟踪框架是对图像序列的两步迭代处理。在第一图像上画出带有控制点的初始轮廓(自动或手动)；然后，针对每一个按顺序的连续的图像，首先对每个控制点单独进行跟踪，同时还记录各向异性的不确定性。作为第二步，利用非正交投影将新的轮廓投影到可行的子空间中。基于训练轮廓获悉可行的子空间，并利用跟踪装置可获得的初始轮廓使可行的子空间与当前情况相适应，同时还考虑初始化的置信度(也就是说，手动初始化比完全自动的初始化具有更高的置信度)。

如上所述，在跟踪过程中采用多个模板。多个模板的使用导致当前情况的形状统计的表征更为准确。在Bayesian核匹配方法中，从已初始化的帧中提取第一模板。当后续模板与现有模板充分不同并且还充分提供用于定位目的的信息时，添加后续模板，这通过与其本身的核匹配来测量：提供更多信息的补丁是当与其本身匹配时具有较高置信度的补丁。

使用一个以上模板的决定基于以下观察结果，即心脏运动是周期性的，因此一个周期内的不同的外观模式都会在后面的周期中再出现。为了与多个模板相匹配，与每个模板进行匹配并选择具有最佳匹配的那个模板。或者，为了节省计算，可以仅仅与以前相匹配的模板相邻的模板进行匹配，从而再次利用周期运动。

在最佳位置的附近计算匹配位置的不确定性，同时使用核匹配方法以相同的方式计算相似图。接着通过例如与二维Gaussian分布拟合或者通过将Hessian矩阵的加权估计的反转，将相似表面用于估计协方差矩阵。

在基于流的方法的情况下，典型的光流实施仅仅利用相邻的帧，从而使对长序列的跟踪过程对误差积累以及偏移敏感。采用多个模板以与用于核匹配相同的方式进行流计算。只要流的不确定性高，就添加局部梯度与现有模板明显不同的新模板。

在得到每个控制点的新位置之后，下一步是通过捕获人类心脏的“合理”形状变化的统计模型来限制整体形状。使用活动形状模型或者基于PCA的形状模型，后者也被称为点分布模型。

由于心脏轮廓的不确定性既不是均一的，也就是说，由于信号失落，一些区域比另一些区域更糟；也不是各向同性的，例如，沿着边缘比沿着梯度方向有更糟的定位。图3显示了这种越过轮廓的各向异性的并且不均一的噪声的例子。如图所示，最初已沿着心内膜壁的轮廓识别出各个点、例如点302、304、306及308。对每一个点进行确定性测量，并环绕每一个点创建一个椭圆，该椭圆指示特定点位于正确位置的确定性水平；椭圆越大，关于特定点的位置的不确定性水平就越高。如图3所示，点302具有环绕它的相对小的椭圆，从而表明关于该点在轮廓上的位置的确定性程度高。点306和308具有环绕它们的相对大的椭圆，从而表明关于这些点的位置的不确定性程度高。

考虑到新轮廓点的一般不确定性模型，最佳投影将不再是正交的。实际上，子空间中的最佳解决方案是在形状模型子空间中相交分布上的最大可能形状。此外，如果在子空间中有分布模型，则没有理由忽视这种额外信息。下面提供关于怎样进行非正交投影的详细分析；以及在可获得子空间模型分布的情况下，怎么将这种模型的信息与不确定输入相融合。

本发明使用基于子空间模型的融合方法。假定同一n维变量x的两次有噪声的测量，每次测量由多维Gaussian分布p₁及p₂来表征，x的最大似然估计是具有与两质心的Mahalanobis距离的最小和的点：x*＝argmin d_m，其中

$d_m={(x-x_1)}^T{C}_1^{-1}(x-x_1)+{(x-x_2)}^T{C}_2^{-1}(x-x_2)---\left(1\right)$

对x求导，得到：

$x^{*}=C(C_1^{-1}x+C_2^{-1}{x}_2)$   (2)

$C={(C_1^{-1}+C_2^{-1})}^{-1}$

这也被称为x的最佳线性无偏估计(BLUE)。

假定Gaussian分布中的一个位于维度为p的子空间中，例如，C₂是奇的。在C₂＝UΛU^T的奇值分解的情况下，其中U＝[u₁，u₂，…，u_n]，u/s为标准正交并且Λ＝diag{λ₁，λ₂，…，λp，0，…，0}，将方程(1)中与x₂的Mahalanobis距离改写为规范形式：

$d_{m,2}={(x-x_2)}^T{C}_2^{-1}(x-x_2)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_i^{-1}{\left[U^T(x-x_2)\right]}^2---\left(3\right)$

当λ_i趋向于0时，d_m，2趋于无穷大，除非U₀ ^TX＝O，其中U₀＝[U_p+1，U_p+2，…，U_n]。(这里我们假定在不丧失一般性的情况下，子空间通过初始空间的原点。由于x₂位于子空间中，所以U₀ ^TX₂＝0。)

因此，仅仅需要在方程(3)中保留那些与非0λ_I对应的项：

$d_{m,2}=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{\mathrm{\λ}}_i^{-1}{\left[U_p^T(x-x_2)\right]}^2\≡{(x-x_2)}^T{C}_2^{+}(x-x_2),---\left(4\right)$

其中C⁺为C的伪逆，并且U_p＝[u₁，u₂，…，u_p]。

此外，由于U₀ ^Tx＝0，所以x可以用另一种形式来表示以便反映这种限制，对于1×p向量y来说：

x＝UU^Tx＝U([U_p|U₀]^Tx)＝U[y|0]＝U_py    (5)

方程(1)现在可以采取下列尽管奇异性但却一般的形式：

$d_m={(U_{p}y-x_1)}^T{C}_1^{-1}(U_{p}y-x_1)+{(U_{p}y-x_2)}^T{C}_2^{+}(U_{p}y-x_2)---\left(6\right)$

对y求导得到

$y^{*}=C_{y^{*}}{U}_p^T(C_1^{-1}{x}_1+C_2^{+}{x}_2)---\left(7\right)$

$C_{y^{*}}={\left[U_p^T(C_1^{-1}+C_2^{+})U_p\right]}^{-1}---\left(8\right)$

$x^{*}=U_p{y}^{*}=C_{x^{*}}(C_1^{-1}{x}_1+C_2^{+}{x}_2)---\left(9\right)$

$C_{x^{*}}=U_p{C}_{y^{*}}{U}_p^T---\left(10\right)$

方程(7)表示信息空间中的融合，随后转换到子空间中。方程(9)为返回到初始空间中的坐标变换。可以看出，C_x*和C_y*是与x*及y*相应的协方差矩阵。

注意，这种解决方案不是通过用伪逆来代替常规的逆的方程(2)的简单概括，这不会将x*限制于子空间中。

作为替代方案，可以将方程(7)和(8)写成

$y^{*}={(U_p^T{C}_1^{-1}{U}_p+{\mathrm{\Λ}}_p^{-1})}^{-1}(U_p^T{C}_1^{-1}{x}_1+{\mathrm{\Λ}}_p^{-1}{y}_2)---\left(11\right)$

这里y₂是x2在由U_p所跨越的子空间中的变换坐标，并且Λ_p＝diag{λ₁，λ₂，...，λ_p}。

方程(11)还可被显示为两种分布的子空间中的BLUE融合，一种分布是N(y₂，Λ_p)，另一种分布为N(x₁，C₁)在子空间中的“交集”(不是投影！)N((U^T _pC^-1 ₁U_p)^-1U^T _pC^-1 ₁x₁，(U^T _pC^-1 ₁U_p)^-1)。

将从大的训练样品库获悉的统计形状模型用于引导特定心脏的轮廓可能有时是有问题的。根据本发明，训练样品被用于获得当前心脏的形状模型，而不是一般心脏的形状模型。因此，存在强烈的动机以使一般模型适应已知的当前情况。根据本发明，对患者心脏的心内膜壁的初始模型(手动或者通过自动检测)进行检查，以调整现有的PCA模型。

在确定心内膜壁的实际轮廓并对其运动进行跟踪时，考虑两种方法：一种方法中初始轮廓是确定性的，另一种方法中初始轮廓是不确定的(这可以是初始轮廓来自自动检测算法的情况，其也提供不确定性)。

当初始轮廓被假定是一个点时(确定性地必然)，使用强适应性PCA(SA-PCA)模型来跟踪该点的运动。假定老的PCA模型(排除当前情况)和针对当前情况初始化了的轮廓共同表示当前情况的变化，但相对能量(即典型功率)分别为α和(1-α)，其中0＜α＜1。换句话说，假定当前情况的形状变化的一部分将在通用模型中被表示，而其余部分在初始轮廓的方向上进行捕获。

PCA模型由分别为X_m、Λ以及U的其平均值、特征值矩阵以及特征向量矩阵表示。如果存储初始全协方差矩阵C(这将是初始维度不是可怕地高的情况)，所调节的平均值和协方差矩阵简单地为两贡献源的加权和：

x_m，new＝αx_m+(1-α)x                   (12)

C_new＝α(C+(x_m-x_m，new)(x_m-x_m，new)^T)+(1-α)(x-x_m，new)(x_m-x_m，new)^T

＝αC+α(1-α)(x-x_m)(x-x_m)^T

                                        (13)

可对C_new进行特征分析以获得新的子空间模型。

在不存储C但在子空间中只有{X_m，Λ，U}是可获取的情况下，通过直接的代数操作，可得到如下的调整后的特征分析结果{X_m，new，Λ_new，U_new}：

初始轮廓x具有诸如x_s＝U^Tx_d的子空间分量，其中x_d＝x-x_m，并且具有诸如x_r＝(x-x_m)-Ux_s的剩余向量。使x_ru为范数为1的x_r的标准化版本(或者如果x_r为零范数，则为零)。

代表所结合的能量的调整后的特征向量矩阵将具有如下形式：

U_new＝[U，x_ru]R                         (14)

R和Λ_new将是下列特征分析问题的解：

$(\mathrm{\α}\left[\begin{array}{cc}\mathrm{\Λ}& 0\\ 0^T& 0\end{array}\right]+\mathrm{\α}(1-\mathrm{\α})\left[\begin{array}{cc}{x}_s{x}_s^T& e_r{x}_s\\ e_r{x}_s^T& e_r^2\end{array}\right])R=R{\mathrm{\Λ}}_{\mathrm{new}}---\left(15\right)$

其中e_r＝X_ru ^T(x-x_m)为剩余能量。

上述公式比IPCA更为常用，具有新数据和老数据之间的可调能量比。如果将α设定为模型中的点比点的总数的分数，这些公式将变得与IPCA等同。通常，由于训练集中的轮廓数通常很大，这将是与1非常接近的数。在α被设置在较小的值(即0.5)的情况下，根据当前情况强烈地调整PCA模型，由此得名。图4显示了具有不同α值的IPCA和SA-PCA的简单的二维图表。

点404代表当前情况，每个x 402表示与特定模型相对应的训练点。椭圆406显示了初始模型分布。椭圆408说明了增量式PCA模型，其对应于强适应性PCA模型，其中α＝0.99。椭圆410显示了强适应性PCA模型，其中α＝0.5。椭圆412显示了强适应性PCA模型，其中α＝0.1。每个椭圆描绘了相应分布的90％等同可能的轮廓。

实际上，与一般训练集中的那些轮廓相比，当前心脏的轮廓更可能与同一心脏的初始轮廓相似，尤其是在当前心脏具有在训练集中不能很好地表示的不规则的形状的情况下。在我们的系统中，将α设定为0.5。这允许来自初始轮廓的较强的影响：模型能量的50％来自初始轮廓(如果与用99例训练的模型一起应用IPCA，这将减少到仅仅1％)。

图5a-c显示了IPCA和SA-PCA的比较。该胸骨旁短轴图具有不规则的形状(带有凹陷部分)，而训练集中圆形形状占绝对多数。考虑初始轮廓(图5a)但具有非常小的权重(＜0.01％)的增量式PCA模型不能捕获当前形状的凹陷特性；而且将轮廓限制为典型的圆形形状(图5b)。事实上，这种结果与使用不具备增量步骤的老PCA模型得到的结果相同。具有α＝0.5的适应性PCA模型更好地匹配真实的边界(图5c)。

利用融合方法的非正交投影与SA-PCA模型调整之间微妙的但重要的相互作用如下：与模型平均值和协方差的融合对于滤出子空间中但距离模型分布太远的轮廓来说是必要的。然而，这种(与正交投影相比)更强的限制将不可避免地改变低的可能性，或者心脏轮廓的界外值。SA-PCA模型通过使用在给定的初始轮廓中提供的额外信息而朝当前情况的方向强烈移动模型来调和这种两难问题。

卡尔曼滤波器融合来自由动态处理所限定的预测的信息和来自有噪声的测量的信息。当用于形状跟踪时，需要另外的全局限制来稳定可行范围内的整体形状。本发明利用统一的融合框架将子空间模型限制合并到卡尔曼滤波器中。

对于卡尔曼滤波器来说，测量更新方程具有如下形式：

x_k+1|k+1＝x_k+1|k+K(z_k+1-Hx_k+1|k)             (16)

其中

K＝P_k+1|kH^T(HP_k+1|kH^T+R)^-1               (17)

P_k+1|k+1＝(I-KH)P_k+1|k                        (18)

P_k+1|k＝SP_k|kS^T+Q                          (19)

这里P为估计误差协方差，而x_i|j为给定在时间j的状态时在时间i的状态估计。测量模型为z_k＝H x_k+r_k，其中r_k表示协方差为R的测量噪声。系统/过程模型为x_k+1＝S x_k+q_k，其中q_k表示协方差为Q的系统噪声。

使用上述方法，通过下式给出具有子空间限制和异方差噪声的卡尔曼滤波器的更新方程：

$x_{k+1|k+1}=P_{k+1|k+1}({(P_{k+1|k+1}+Q)}^{+}{x}_{k|k}+R^{-1}z+C_2^{+}{x}_2)---\left(20\right)$

$p_{k+1|k+1}=U_p{\left[U_p^T({(P_{k|k}+Q)}^{+}+R^{-1}+C_2^{+})U_p\right]}^{-1}{U}_p^T---\left(21\right)$

其中假定系统噪声协方差Q被包含在子空间中。观察这种结合信息空间中所有可获取的知识的解决方案的对称性。这些方程提供了系统动态、子空间限制和噪声信息的统一的融合。它们代表影响跟踪系统的各种不确定性的全部表示。

根据本发明的一种替代的实施例，使用一种跟踪技术，其依赖于基于控制点的对象表示以及用于结合多帧的模型信息的鲁棒的融合。在不同时间维持并获取基于一组控制点的对象表示。鲁棒地结合由控制点建议的估计运动来确定对象的下一位置。通过随着时间保持几个模型来获得对对象的视觉跟踪。该结果是表征对象外观的概率密度函数的非参数表示。通过独立地从每个模型通过光流获得运动估计及其不确定性来进行跟踪。利用鲁棒的融合技术、例如基于可变带宽密度的融合(Variable-Bandwidth Density-based Fusion(VBDF))程序来计算对每个控制点的最终估计。VBDF计算位移密度函数的最有意义模式的位置，同时考虑它们的不确定性。VBDF程序管理多个数据源以及位移估计中的界外值。通过大的残留误差的估计不确定性，自然地对掩蔽进行处理。模型被分成多个区域，针对这多个区域对流独立地进行计算。残留校准误差被用于计算估计的协方差矩阵的比例，因而减小了不可靠的位移的影响。

在本发明中，通过基于多重模态控制点的方法来处理由于掩蔽和外观变化而引起的对象跟踪挑战。保持二维外观模型的几个代表并不将其限制为单峰分布，并且VDBF机制鲁棒地结合多个估计以确定每个控制点的最主要的运动。为了在跟踪过程中对变化进行建模，保持对象外观随着时间的多个样本。

图10说明了根据本发明的示例性的基于多模型控制点的跟踪装置。顶部的三个帧1002、1004及1006说明了模型集中的当前样本，每个当前样本与一组重叠的控制点相关联。基于控制点的方法比整体表示更鲁棒，并且因此对照明变化以及姿势不太敏感。本发明的另一个优点在于，可通过分析匹配可能性在控制点水平对部分掩蔽进行处理。

独立地对每个控制点进行处理；根据所有的模型模板，在当前图像中对其位置和协方差矩阵进行估计。例如，控制点中的一个由灰色矩形1014来示出，其相对于每一模型的位置及不确定性在I_new 1008中示出。如帧1010中所示，VBDF鲁棒融合程序被用于确定具有相关的不确定性的最主要的运动(模式)。注意由于掩蔽或者外观变化而引起的每一控制点的估计位置的方差。

控制点在当前帧中的位置还通过全局参数运动模型来限制。利用每个控制点位置的置信度来估计相似变换模型及其参数。因此，可靠的控制点对全局运动估计的贡献更大。如果参考外观的残余误差相对低，则将当前帧添加到模型集中。这样选择阈值，使得在对象具有明显的掩蔽的情况下不添加图像。模型中模板的数目可以是可变的或固定的。如果模板的数目是固定的，则可提供抛弃某些模板的方案(例如，抛弃最老的模板)。

一种VBDF估计器基于具有适应性核带宽的非参数密度估计。该VBDF估计器被限定为密度函数的最有意义模式的位置。模式计算是基于多标度最佳化框架中的可变带宽平均移动技术。

使x_i∈R^d，i＝1...n为可得到的d维估计，每个估计具有由协方差矩阵C_i所给定的相关的不确定性。以多标度方式迭代地确定密度函数的最有意义模式。带宽矩阵H_i＝C_i+α²I与每一个点x_i有关，其中I为单位矩阵，参数α决定了分析的标度。在位置x处的采样点密度估计量由下式限定：

$\hat{f}\left(x\right)=\frac{1}{n{\left(2\mathrm{\π}\right)}^{d/2}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\exp (-\frac{1}{2}{D}^2(x,x_i,H_i))---\left(22\right)$

其中D表示x与x_i之间的Mahalanobis距离

$D^2(x,x_i,H_i)={(x-x_i)}^T{H}_i^{-1}(x-x_i)---\left(23\right)$

在位置x处的可变带宽平均移动向量由下式给出：

$m\left(x\right)=H_h\left(x\right)\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i\left(x\right)H_i^{-1}{x}_i-x---\left(24\right)$

其中H_h表示用依赖于数据的权重ω_i(x)加权的带宽矩阵的调和平均值

$H_h\left(x\right)={\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i\left(x\right)H_i^{-1}\right)}^{-1}---\left(25\right)$

在当前位置x计算出的依赖于数据的权重具有下列表达式：

${\mathrm{\ω}}_i\left(x\right)=\frac{\frac{1}{{\left|H_i\right|}^{1/2}}\exp (-\frac{1}{2}{D}^2(x,x_i,H_i))}{{\mathrm{\Σ}}_{i=1}^n\frac{1}{{\left|H_i\right|}^{1/2}}\exp (-\frac{1}{2}{D}^2(x,x_i,H_i))}---\left(26\right)$

并且注意，它们满足 ${\mathrm{\Σ}}_{i=1}^n{\mathrm{\ω}}_i\left(x\right)=1.$

可以看出，对应于点x+m(x)的密度总是大于或等于对应于点x的密度。因此，利用平均移动向量对该位置进行重复更新，得到爬山过程，其集中于具有基本密度的固定点。

VBDF估计器通过以多个标度重复应用适应性平均移动过程而找到最重要的模式。其通过根据点x_i的分布选择大的参数α而从大标度开始。在这种情况下，密度表面为单峰的，因而所确定的模式将对应于整体上密度最大的区域。在减小参数α的值并且从以前一标度所确定的模式开始平均移动迭代时，重复该过程。对最后一个步骤来说，与每一个点相关的带宽矩阵与协方差矩阵相等，即H_i＝C_i。

VBDF估计器是一种用于信息融合的强有力工具，其具有处理多源模型的能力。由于在局部邻域中的点可能表现出多种运动，所以这一点对运动估计来说是重要的。最有意义模式对应于最相关的运动。

考虑具有n个模型M₀，M₁，...M_n。针对每个图像，维持c个成分，它们的位置由x_ij表示，i＝1...n，j＝1...c。当可得到新的图像时，对每一成分及每一模型的位置和不确定性进行估计。这一步骤可利用多种技术、例如基于图像相关性、空间梯度或时空能量的调整的技术来完成。利用VBDF技术，结果为每一成分及其不确定性C_ij的运动估计x_ij。因此x_ij表示成分j相对于模型i的位置估计。协方差矩阵的标度还可根据匹配的残余误差来估计。当相应的控制点被掩蔽，因而在控制点水平对掩蔽进行处理时，这将增大协方差矩阵的尺寸。

VBDF鲁棒融合技术被用于确定成分j在当前帧中的最相关的位置x_j。通过标度的模式跟踪得出：

${\hat{x}}_j=C\left({\hat{x}}_j\right)\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i\left({\hat{x}}_j\right){\hat{C}}_{\mathrm{ij}}^{-1}{\hat{x}}_{\mathrm{ij}}$

$C\left({\hat{x}}_j\right)={\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i\left({\hat{x}}_j\right){\hat{C}}_{\mathrm{ij}}^{-1}\right)}^{-1}$   (27)

其中权重ω_i如在公式(26)中那样定义：

${\mathrm{\ω}}_i\left({\hat{x}}_j\right)=\frac{\frac{1}{{\left|C_{\mathrm{ij}}\right|}^{1/2}}\exp (-\frac{1}{2}{D}^2({\hat{x}}_j,{\hat{x}}_{\mathrm{ij}},C_{\mathrm{ij}}))}{{\mathrm{\Σ}}_{i=1}^n\frac{1}{{\left|C_{\mathrm{ij}}\right|}^{1/2}}\exp (-\frac{1}{2}{D}^2({\hat{x}}_i,{\hat{x}}_{\mathrm{ij}},C_{\mathrm{ij}}))}---\left(28\right)$

并且注意，它们满足 ${\mathrm{\Σ}}_{i=1}^n{\mathrm{\ω}}_i\left(x\right)=1.$

可以看出，对应于点x+m(x)的密度总是高于或等于对应于点x的密度。因此，利用平均移动向量对当前位置进行重复更新产生爬山过程，其集中于具有基本密度的固定点。

通过以多个标度迭代地应用适应性平均移动过程，VBDF估计器找到最重要的模式。其通过根据点x_i的分布选择大的参数α而从大的标度开始。在这种情况下，密度表面是单峰的，因而所确定的模式将对应于整体上密度最大的区域。在减小参数α的值并且从以前一标度所确定的模式开始平均移动迭代时，重复该过程。对最后一个步骤来说，与每一个点相关的带宽矩阵与协方差矩阵相等，即H_i＝C_i。

VBDF估计器是一种用于信息融合的强有力的工具，其具有处理多源模型的能力。由于在局部邻域中的点可能表现出多种运动，所以这一点对运动估计来说是重要的。最有意义模式对应于最相关的运动。

现在将描述根据本发明的跟踪多个控制点模型的例子。有n个模型M₀，M₁，...M_n。针对每个图像，保持c个控制点，其中用x_ij，i＝1...n，j＝1...c来表示它们的位置。当可得到新的图像时，对每一控制点及每一模型的位置和不确定性进行估计。这一步骤可利用多种技术、例如基于图像相关性、空间梯度或时空能量的调整的技术来完成。使用VBDF技术，结果为每一控制点及其不确定性C_ij的运动估计x_ij。因此x_ij表示成分j相对于模型i的位置估计。协方差矩阵的标度还可根据匹配的残余误差来估计。当相应的控制点被掩蔽，因而在控制点水平对掩蔽进行处理时，这将增大协方差矩阵的尺寸。

VBDF鲁棒融合技术被用于确定成分j在当前帧中的最相关的位置x_j。通过标度的模式跟踪得出：

${\hat{x}}_j=C\left({\hat{x}}_j\right)\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i\left({\hat{x}}_j\right){\hat{C}}_{\mathrm{ij}}^{-1}{\hat{x}}_{\mathrm{ij}}$

$C\left({\hat{x}}_j\right)={\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i\left({\hat{x}}_j\right){\hat{C}}_{\mathrm{ij}}^{-1}\right)}^{-1}$    (29)

其中权重ω_i如在公式(28)中那样定以：

${\mathrm{\ω}}_i\left({\hat{x}}_j\right)=\frac{\frac{1}{{\left|C_{\mathrm{ij}}\right|}^{1/2}}\exp (-\frac{1}{2}{D}^2({\hat{x}}_j,{\hat{x}}_{\mathrm{ij}},C_{\mathrm{ij}}))}{{\mathrm{\Σ}}_{i=1}^n\frac{1}{{\left|C_{\mathrm{ij}}\right|}^{1/2}}\exp (-\frac{1}{2}{D}^2({\hat{x}}_i,{\hat{x}}_{\mathrm{ij}},C_{\mathrm{ij}}))}---\left(30\right)$

在每个控制点的计算之后，利用由估计的协方差矩阵所给定的权重执行加权的矩形拟合。图像补丁通过由4个参数限定的相似变换T相关联。动态控制点位置x的相似变换由下列方程来表征：

$T\left(x\right)=\left(\begin{array}{cc}a& -b\\ b& a\end{array}\right)x+\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)---\left(31\right)$

其中t_x、t_y为平移参数，并且a、b用参数表示2D旋转和比例缩放。

最小化的标准为参考位置x⁰ _j与估计位置x_j(当前帧中第j个控制点位置)之间的Mahalanobis距离之和。

${({\hat{x}}_j-T\left(x_j^0\right))}^{T}C{\left({\hat{x}}_j\right)}^{-1}({\hat{x}}_j-T\left(x_j^0\right))---\left(32\right)$

通过标准加权的最小平方来完成最小化。注意，由于协方差矩阵被用于每一控制点，所以具有高不确定性的点的影响被减小。

在矩形与所跟踪的控制点拟合之后，对矩形内的动态成分待选者均一地进行再采样。假定每个控制点相对于矩形的相对位置不会改变很多。如果再采样位置与通过某一控制点的光流计算出的跟踪位置之间的距离大于可容许的阈值，则跟踪位置被认为是界外值并用再采样点来代替。如果有足够的控制点具有低的残余误差，则将当前图像添加到模型集中。将模型与当前帧之间的中值残余误差与预定阈值Th相比较。

假设一组模型M₀，M₁，...，M_n，其中成分j在帧i中具有位置x_ij，我们的对象跟踪算法可通过下列步骤来概括：

1.给定一幅新的图像I_f，通过始自在前一帧中所估计的位置_j ^(f-1)的鲁棒光流[8]来计算_ij ^(f)；

2.对于j＝1...c，利用VBDF估计器(Subsection 3.2)估计成分j的位置_j ^(f)，得出(6)；

3.利用通过最小化(9)所计算的变换来限制成分的位置；

4.如果其中值残余误差小于T_h，则将新的外观添加到模型集中。

所建议的多模板框架可直接用于形状跟踪的环境。如果所跟踪的点表示通过样条函数来建模的形状的控制点，则多位置估计的鲁棒融合的使用增加了形状的位置估计的可靠性。当形状空间受限于所获悉的子空间限制时，其还导致较小的校正。如果可获取轮廓，用于跟踪的模型模板则可从基于形状之间的距离的模型集中在线地进行选取。

图6a及6b示出了左心室的超声波心动描记图像，其中根据本发明对心内膜壁的运动进行跟踪。每个图602、610的帧1说明了壁的初始轮廓，如图像中的点所示。后续帧604、606、608、612、614、616说明了如何随着时间对壁的运动进行跟踪。根据本发明对每一个点进行测量。

实例

在本实例中将在超声波心动描记图像中手动跟踪的心脏轮廓用作训练集。心尖四腔(开放轮廓)与胸骨旁短轴(封闭轮廓)图都如图7a-7d以及图8a-8d所示的那样进行测试。还针对每个轮廓对地标点作注释。图7b和图8b显示了画在一起的训练轮廓集。在每个轮廓是以有序地标点的坐标作为其分量的向量的情况下(开放轮廓具有34个分量，封闭轮廓具有36个分量)，应用特征分析方法。然后对列为训练向量的矩阵执行PCA。并且将特征值用于形成作为子空间中的模型协方差矩阵的诊断矩阵。调节测量协方差矩阵来测试不同的情形。

图7a和8a示出了测试轮廓，其中实线曲线704、802为真值(ground truth)，而虚曲线702、804为有噪声的测量。期望在特征轮廓子空间中找到与当前有噪声的轮廓最接近的轮廓，其中更多不确定的点需要更多的调节。图7c和图8c显示了当将正交投影应用于图7a和图8a的轮廓时的结果(其中各向同性的协方差由围绕地标点的小圆712、812表示)。图7d和图8d显示了结果。可以看出，该结果更接近于真值。由于训练数据与测试数据相当不同，因而不是完美的拟合。因此，测试轮廓中一些小的形状变形可能在子空间中是不能实现的。

根据本发明，可利用显像工具对轮廓进行观察以帮助诊断。为了使医生更容易地从超声波心动描记图像中诊断出心脏疾病，通过彩色向量实时地显示LV轮廓。向量的长度表示LV壁轮廓上的点的位移大小，而向量的方向表示LV壁轮廓上的点的运动方向。出于示例的目的，橙色被用于描绘收缩运动，深蓝色被用于描绘扩张运动。图9显示了一种表示结果。LV轮廓上的点的运动被跟踪。通过Gaussian在时域和空间域中对运动进行平滑。为了使医生更容易看到，向量的长度被放大三倍。

通过这种显像方法，医生可以很容易地看到心内膜的每段的具有大小和方向的运动。结合全局运动补偿和这种显像方法可使医生很容易看到LV的所有段中的收缩方向和大小。局部缺血区域或其它异常区域可直接并且容易地用人眼识别出来。

已经描述了在医学检查过程中确定特征灵敏性的方法的实施例，需要说明的是，根据上述教导，本领域的普通技术人员可进行修改和变化。因此应当理解，可在由所附的权利要求书所限定的本发明的精神和范围内公开的本发明的特定实施例中作出改变。因此，已经描述了具有专利法所要求的细节及特征的本发明，该特征在所附的权利要求书中阐明，其要求并期望受到专利特许证的保护。

Claims (20)

1.一种用于跟踪运动对象的整体形状的方法，该方法包括如下步骤：

沿着所述整体形状定义一个或多个控制点；

当所述对象运动时对所述一个或多个控制点中的每一个进行跟踪；

表示运动中的控制点的位置的不确定性；以及

使用先验的形状模型，利用所述不确定性来限制所述整体形状。

2.根据权利要求1的方法，其特征在于，对所述一个或多个控制点中的每一个进行跟踪的步骤采用Bayesian核匹配方法。

3.根据权利要求1的方法，其特征在于，对所述一个或多个控制点中的每一个进行跟踪的步骤采用基于光流的方法。

4.根据权利要求1的方法，还包括如下步骤：

在对所述一个或多个控制点进行跟踪的过程中，通过使用多模板适应性匹配框架来校正误差积累。

5.根据权利要求1的方法，其特征在于，所述先验的形状模型为子空间形状限制模型。

6.根据权利要求1的方法，其特征在于，强适应性主成分分析(SA-PCA)模型被用于限制形状。

7.根据权利要求1的方法，其特征在于，利用所述不确定性的步骤还包括使用协方差矩阵以便采用利用非正交投影或者信息融合的子空间形状限制模型。

8.一种用于跟踪运动对象的整体形状的系统，包括：

用于沿着所述整体形状的初始轮廓定义一个或多个控制点的装置；

用于在所述对象运动时对所述一个或多个控制点中的每一个进行跟踪的装置；

用于表示运动中的控制点的位置的不确定性的装置；以及

用于在使用先验的形状模型时利用所述不确定性来限制整体形状的装置。

9.根据权利要求8的系统，其特征在于，跟踪装置采用Bayesian核匹配方法。

10.根据权利要求8的系统，其特征在于，跟踪装置采用基于光流的方法。

11.根据权利要求8的系统，还包括：

用于在对所述一个或多个控制点进行跟踪的过程中通过使用多模板适应性匹配框架来校正误差积累的装置。

12.根据权利要求8的系统，其特征在于，将强适应性主成分分析(SA-PCA)模型用于限制轮廓的形状。

13.根据权利要求8的系统，其特征在于，利用所述不确定性的步骤还包括使用协方差矩阵以便采用利用非正交投影或者信息融合的子空间形状限制模型。

14.一种用于视觉跟踪对象形状的运动的方法，包括如下步骤：

产生一个或多个第一彩色向量来表示控制点沿着形状的轮廓的收缩；

产生一个或多个第二彩色向量来表示控制点沿着形状的轮廓的扩张；

周期性地显示第一和第二彩色向量，从而显示形状的运动。

15.一种用于视觉跟踪对象形状的运动的系统，包括：

用于产生一个或多个第一彩色向量来表示控制点沿着形状的轮廓的收缩的装置；

用于产生一个或多个第二彩色向量来表示控制点沿着形状的轮廓的扩张的装置；

用于周期性地显示第一和第二彩色向量从而显示形状的运动的装置。

16.一种用于跟踪运动对象的整体形状的方法，该方法包括如下步骤：

沿着所述整体形状定义一个或多个控制点；

当所述对象运动时对所述一个或多个控制点中的每一个进行跟踪；

针对每一个控制点建立多个外观模型；以及

结合由每一个模型所产生的运动向量以便跟踪形状。

17.根据权利要求16的方法，还包括如下步骤：

利用鲁棒的信息融合来计算每个控制点的位置的估计值。

18.根据权利要求16的方法，其特征在于，基于一个或多个控制点的位置的不确定性来建立所述多个外观模型。

19.根据权利要求18的方法，其特征在于，所述鲁棒的信息融合是基于可变带宽密度的融合(VBDF)。

20.根据权利要求19的方法，其特征在于，所述鲁棒的信息融合包括下列步骤：

输入由平均向量和协方差矩阵所表示的多个测量结果；

通过计算可变带宽密度估计值，将所述测量结果融合在一起；以及

检测可变带宽密度估计值的最有意义的节点。

Images