CN101639949B - 基于时变插值算法的柔性体点分布模型的建立方法 - Google Patents

Abstract

一种基于插值算法的柔性体点分布模型的形状、形变及协方差矩阵计算方法，包括以下步骤：1)载入两个相邻时刻T_a和T_b的柔性体分布模型M_a和M_b，其中：M_a＝<<overscore>y</overscore>_a，Φ_a，λ_a>，M_b＝<<overscore>y</overscore>_b，Ф_b，λ_b>；2)通过所述柔性体分布模型M_a和M_b，分别计算平均形状<overscore>y</overscore>_a和<overscore>y</overscore>_b、特征矩阵Φa和Фb、以及特征值λ_a和λ_b，设定柔性体从T_a到T_b呈连续变化；2.1)中间模型的平均形状；2.2)中间模型的形变；2.3)中间模型的协方差矩阵；3)根据得到中间模型的平均形状、形变和协方差矩阵，得到t时刻的中间模型M_t。本发明能够形成实时分布的中间过渡模型、模型精度高、实用性良好。

Description

基于时变插值算法的柔性体点分布模型的建立方法

技术领域

本发明涉及柔性体分布模型，结合计算机图像处理技术、建模技术的柔性体分布模型的参数计算方法。

背景技术

柔性体相对于刚性体，在计算机中进行建模和分析比较复杂，而在现实生活中物体随时间柔性变化的情况是非常普遍的。对于许多生物医学组织，如人的心脏，双手，骨骼，肾脏、大脑、细胞等等，这些人们可以很好地描述其外部形状特征但却很难给出准确的刚性模型。点分布模型(PDM)是一种新近崛起的技术，它能很好的对研究对象进行描述，可以成功地应用于非刚性模型的分解与合成处理。近年来，柔性体的统计学建模及其与图像解释技术的结合催生了一些高级应用软件，这些软件不仅针对生物医学，还涉及目标跟踪、识别、以至影视等领域。

PDM是一种基于对象样本分析的统计模型。从Cootes等人开始，这种基于训练样本上标记点的形体建模思想一直被贯彻下来。最典型的想法是，使标记点自动对齐以确定图形的位置和形变类型。对齐图像后可以很容易的比较同组标记点在不同例图中的位置变化，只需要计算它们的坐标就可以了。因而，PDM模型对于形状的描述是通过轮廓点向量来完成的，它由中间形状和形状变化样式组成。通过对这种分布情况建模，我们可以得到与原始训练集相似的新形状。

在三维空间，如果训练集图像中的每个图形都由一组数量为n的标记点描述，那么它可以用一个3n维的矢量表示(如果是二维对象即为2n维矢量)。训练集中的标记工作非常重要，通常是手动完成的。另外，在研究中我们要收集数百名例柔性体的个例以得到比较好的统计结果。训练集形成的概率模型维度可能非常的大(在一个心脏模型中通常可达到上万维)，所以建立PDM的过程是非常枯燥又费时费力的。然而，在一些应用中，单独的一个PDM模型或是有限个PDM模型都不足以很好的描述对象变化。事实上，不论是连续PDM模型还是模型***方法都是为了满足对象的性质，即对象时间和空间上的变化。例如在心脏建模时，需要不同时间段的模型才能很好地描述心室形状在心脏跳动不同阶段的变化。但是，在构建整个心脏运动周期的连续模型时，我们只能由已知的模型序列得到特定几个时刻的模型，而其他任意时刻的模型就不得而知了，尽管在基于模型的图像分割领域PDM建模出现已有十来年，但对于模型插值及相关内容的深入研究还没有开展过。

发明内容

为了克服已有的柔性体点分布模型的不能形成实时分布的模型、模型精度差、实用性差的不足，本发明提供一种能够形成实时分布的中间过渡模型、模型精度高、实用性良好的基于时变插值算法的柔性体点分布模型的建立方法。

本发明解决其技术问题所采用的技术方案是：

一种基于时变插值算法的柔性体点分布模型的建立方法，包括以下步骤：

1)、载入两个相邻时刻T_a和T_b的柔性体分布模型M_a和M_b，其计算式为：

$M_a=<{\stackrel{\&OverBar;}{y}}_a,{\mathrm{\&Phi;}}_a,{\mathrm{\&lambda;}}_a>,$ $M_b=<{\stackrel{\&OverBar;}{y}}_b,{\mathrm{\&Phi;}}_b,{\mathrm{\&lambda;}}_b>;$

2)、通过所述柔性体分布模型M_a和M_b，分别计算平均形状和

特征矩阵Φa和Φb、以及特征值λ_a和λ_b，设定柔性体从T_a到T_b呈连续变化，

2.1)中间模型的平均形状

用线性插值得到y_ai和y_bi之间的中间模型形状y_ti为：

y_ti＝ay_ai+by_bi                       (12)

a和b由式(10)和(11)定义：

$a=\frac{t-T_a}{T_b-T_a}---\left(10\right)$

b＝1-a                               (11)

a和b反映了中间模型相对于M_a和M_b的距离；

***的平均形状为：

${\stackrel{\&OverBar;}{y}}_t=\frac{1}{s}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{y}_{\mathrm{ti}}=\frac{1}{s}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}({\mathrm{ay}}_{\mathrm{ai}}+{\mathrm{by}}_{\mathrm{bi}})=\frac{a}{s}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{y}_{\mathrm{ai}}+\frac{b}{s}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{y}_{\mathrm{bi}}=a{\stackrel{\&OverBar;}{y}}_a+{b\stackrel{\&OverBar;}{y}}_b---\left(13\right);$

2.2)、中间模型的形变

中间模型的形变由下式算出：

${\mathrm{\&Delta;}}_{\mathrm{yti}}=y_{\mathrm{ti}}-{\stackrel{\&OverBar;}{y}}_t=({\mathrm{ay}}_{\mathrm{ai}}+{\mathrm{by}}_{\mathrm{bi}})-(a{\stackrel{\&OverBar;}{y}}_a+b{\stackrel{\&OverBar;}{y}}_b)$

$=a(y_{\mathrm{ai}}-{\stackrel{\&OverBar;}{y}}_a)+b(y_{\mathrm{bi}}-{\stackrel{\&OverBar;}{y}}_b)---\left(14\right)$

$=a{\mathrm{\&Delta;}}_{\mathrm{yai}}+b{\mathrm{\&Delta;}}_{\mathrm{ybi}}$

2.3)、中间模型的协方差矩阵

中间模型的协方差矩阵由下式推导得出：

$S_t=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}(y_{\mathrm{ti}}-{\stackrel{\&OverBar;}{y}}_t){(y_{\mathrm{ti}}-{\stackrel{\&OverBar;}{y}}_t)}^T=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yti}}{\mathrm{\&Delta;}}_{\mathrm{yti}}^T$

$=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}({\mathrm{a\&Delta;}}_{\mathrm{yai}}+{\mathrm{b\&Delta;}}_{\mathrm{ybi}}){({\mathrm{a\&Delta;}}_{\mathrm{yai}}+{\mathrm{b\&Delta;}}_{\mathrm{ybi}})}^T---\left(15\right)$

$=\frac{a^2}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yai}}{\mathrm{\&Delta;}}_{\mathrm{yai}}^T+\frac{b^2}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}^T+\frac{\mathrm{ab}}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yai}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}^T+\frac{\mathrm{ab}}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}{\mathrm{\&Delta;}}_{\mathrm{yai}}^T$

令

$S_{\mathrm{ab}}=\frac{\mathrm{ab}}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yai}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}^T---\left(16\right)$

因为

$S_a=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}(y_{\mathrm{ai}}-{\stackrel{\&OverBar;}{y}}_a){(y_{\mathrm{ai}}-{\stackrel{\&OverBar;}{y}}_a)}^T=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yai}}{\mathrm{\&Delta;}}_{\mathrm{yai}}^T---\left(17\right)$

$S_b=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}(y_{\mathrm{bi}}-{\stackrel{\&OverBar;}{y}}_b){(y_{\mathrm{bi}}-{\stackrel{\&OverBar;}{y}}_b)}^T=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}^T---\left(18\right)$

显然有

$S_{\mathrm{ab}}=S_{\mathrm{ba}}^T---\left(19\right);$

3)、根据得到中间模型的平均形状、形变和协方差矩阵，得到t时刻的中间模型Mt。

本发明的技术构思为：PDM方法是一种先验建模方法，它建立在目标对象的一组样本(即训练集)已知的基础上。这组样本集可以提供对象的形状及形变的统计学描述。在实际应用中，这组对象样本的形状描述通常是使用轮廓来进行的(如图片上的一组像素的坐标值)。进一步的，可以在轮廓上选取一些标记点来描述目标物形状，而这些标记点往往与目标物的某种特征相关联。

不同样本上的对应标记点在选取的时候遵循着“特征位置大致相同”的原则，但往往这种定位是不够精确的，存在一定偏差。这种偏差要归因于不同个体间的自然形变。我们期望相对于整个形体的尺度而言，这种偏差是十分微小的。PDM方法允许我们模仿这种微小的差异，并且指出哪些偏差的确是微小的，哪些相对而言更大一些。

由于训练集中的形体是从不同的图片集中获得的，而不同图片集的坐标系并不统一，因此有必要先将所有训练集形体进行坐标归一化，即将它们对齐。对齐的过程实际上就是对每一个样本模型寻找适当的相似性变换的过程，包括平移、缩放和旋转，以使得不同样本之间尽可能的接近，另一方面，选定的变换也应该令各个样本模型与由所有样本得到的平均模型尽可能接近。为了方便表述，假设我们的样本只有两个，归一化简化为对齐这两个样本。每个样本的形状都由一个标记点坐标组成的向量来表示：

x＝(x₁，y₁，z₁，...，x_n，y_n，z_n)^T，                (1)

其中(x_i，y_i，z_i)是某个标记点在三维图像中的空间坐标。如果是二维医学图像，则标记点用(x_i，y_i)表示。

图片上标记点的获取有多种方式，人工手动标记或采用一定得自动标记算法都是可以的。例如，Izard等人提出的一种MRI图像标记方法，使用统一的算法规则在不同人脑图片中进行组织结构标记。得到一组经过标记的训练样本后，我们需要将它们对齐到一个统一的坐标系下。可以使用广义对齐算法(GPA)来将训练样本对齐，并使得所有样本形状到平均模型的距离平方和最小，实际上也就是寻找变换T_i(平移，旋转，缩放)，使下式最小化

D_m＝∑|m-T_i(x_i)|²，                (2)

其中

$m=\frac{1}{n}\mathrm{\&Sigma;}{T}_i\left(x_i\right),\left|m\right|=1,---\left(3\right)$

且x_i是训练集中的一个形状样本，它是一个3n维的向量表示(在三维空间中有n个标记点的情况下)。

对齐后的训练集组成了一个三维空间上的点云，可以看作是一个概率密度函数。为了降低运算成本和内存占用，我们采用主成分分析法(PCA)寻找点云中的主元，并以前列最大的少数分量建立模型。这些被选出的分量可以描述对象主要的形变。以心动周期内左心室的模型为例，我们经过主成分分析后的用60个特征向量已经能足够描述心室的形状及其变化了。最后，PDM模型可以表达为：

$x=\stackrel{\&OverBar;}{x}+p_1{b}_1+...+p_t{b}_t=\stackrel{\&OverBar;}{x}+\mathrm{\&Phi;c}---\left(4\right)$

其中，

是对齐后的平均形状向量，Φ是一个3n×t的矩阵，Φ的每一列表示主成分方向上的单位向量，c是由形状参数组成的t维列向量。可见，经过PCA分析后形状的维度已经由3n下降到了t。

经过以上步骤我们得到了一个类似于点分布模型的统计学模型，该模型可应用于柔性体的物体建模，以及在新的图像中定位和分割新个例。在训练集的学习过程中我们可以得到形状参数b的变化范围，只要在该范围内变化，任意b的值都可以生成合理的新样本(仿真形体)。通常b_i的变化为λ_j(矩阵Φ中对应第i大的特征值)或者

从样本集中训练出一个统计学形态模型后，就可以用PDM来解释新的图像个例了。为了将模型与图像匹配，通常采用一种主动重复迭代的方法。该方法使模型不断地变形来适应新图像中目标对象的轮廓。模型形状的变化受到一些统计数据的约束，也即只能在训练样本集中所得出的范围内进行变化。对于形状模型，我们还要求图像轮廓大致位于模型点的附近。可以表示为某点的梯度与轮廓上过该点的法线方向的变化。令统计模型的轮廓沿图像轮廓不断移动，同时计算两轮廓的马氏距离最小化就可以最终得到真实的边界位置。因此，要在特定实例中获得目标柔性对象的形状需要反复执行以下两步直到收敛：(1)沿模型各点的法线寻找图像轮廓上相匹配的对应点(使得模型点与对应图像点的马氏距离最小)；(2)改变模型位置及形状参数使得模型各标记点更接近图像轮廓上找到的对应点。也就是说，模型的拟合过程实际上就是在图像上搜寻模型点的最相似位置及不断修正整体几何变换T及形状参数c的过程，使之最小化。数学上可以表示为：

$f=|{X-T(\stackrel{\&OverBar;}{x}+\mathrm{Pc})|}^2=f(b,X_c,Y_c,Z_c,s,\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&gamma;})$ (5)

$={|X-T(\stackrel{\&OverBar;}{x}+\mathrm{Pc};X_c,Y_c,Z_c,s,\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&gamma;})|}^2$

其中X是在当前迭代时得到的临时模型。

这是一个寻找最小值的问题，可以通过一些非线性优化的迭代方法求解。最终我们可以确定整体变换参数T，以及个体实例的形状参数b：

$c=P^T(T^{-1}\left(X\right)-\stackrel{\&OverBar;}{x}).---\left(6\right)$

近几年，人们已经就基于PDM的柔性体建模方法的许多相关问题进行了深入的研究，如形体对齐、自动标记、平均形体生成、形变建模、模态分析、图像分割等等问题。此外，在该思想广泛地应用于非刚体对象建模的基础上，并可进一步扩展为诸如活动形体模型(ASM)和活动表现建模(AAM)等以应用于更深入的图像分析。目前，对PDM模型最成功的表达是通过PCA方法得到的。PCA方法非常简单易懂，首先把一个形体投影到经过学***均形状

和形体协方差矩阵S：

$\stackrel{\&OverBar;}{y}=\frac{1}{s}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{y}_i,$ $S=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}(y_i-\stackrel{\&OverBar;}{y}){(y_i-\stackrel{\&OverBar;}{y})}^T---\left(7\right)$

这样，一个合理的形状可以被表达为

$y=\stackrel{\&OverBar;}{y}+\mathrm{\&Phi;c}---\left(8\right)$

式中向量c由对应于各个特征向量的形变参数组成。矩阵Φ则包含了一组由协方差矩阵S分解得到的特征向量，这组向量可线性组合成任意的新形状y。根据PCA理论，c中的每个参数都控制形体在某一方向上的形变，并且，各个方向都彼此独立，参数按照形体在对应方向上形变由大到小顺序排列。本文考虑在一个PDM模型仅包含平均形状、特征值和特征向量而不包括原始训练集的数据的情况下进行插值运算。即：

$M=<\stackrel{\&OverBar;}{y},\mathrm{\&Phi;},\mathrm{\&lambda;}>---\left(9\right)$

如上所述，单个PDM模型可以由一组柔性对象的形体获得，例如，我们研究的医院患者的心脏图片(一些被观测的柔性体个例)。然而，有时我们要观测这些柔性体随时间的变化情况，这样可以得到一系列的形状信息，例如心脏在整个活动周期内的三维变形。而医院现有设备只能获得若干时刻的三维图像，即对应于某些Ti时刻的模型Mi。

对一个特定的动态形变物体，假设我们已经在离散的时间坐标(T1，T2，...)上建立起一组PDM模型(M₁，M₂，...)。我们可以把问题描述为：给出已知的两个PDM模型，T_a时刻的模型M_a和T_b时刻的模型M_b(模型已知意味着我们知道了平均形状

和

特征矩阵Φa和Φb、以及特征值λ_a和λ_b)，如何得到在t(T_a≤t≤T_b)时刻的中间模型M_t的相关参数。这些参数有平均形状、特征值和特征向量，即(图1)。这里我们假设在T_a到T_b的时间段中，一个形体从Y_ai连续变化到Y_bi，二者都是n维的，如Y_ai＝(x_ai1，x_ai2，x_ai3，...，x_ain)^T。这里还定义了两个常数：

$a=\frac{t-T_a}{T_b-T_a}---\left(10\right)$

b＝1-a                (11)

a和b反映了中间模型相对于M_a和M_b的距离。

(1)中间模型的平均形状

根据以上的定义和假设，用线性插值(可满足大部分的实际应用)得到Y_ai和Y_bi之间的中间模型形状为：

y_ti＝ay_ai+by_bi                        (12)

a和b由式(10)和(11)定义，即a+b＝1。

***的平均形状为：

${\stackrel{\&OverBar;}{y}}_t=\frac{1}{s}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{y}_{\mathrm{ti}}=\frac{1}{s}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}({\mathrm{ay}}_{\mathrm{ai}}+{\mathrm{by}}_{\mathrm{bi}})=\frac{a}{s}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{y}_{\mathrm{ai}}+\frac{b}{s}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{y}_{\mathrm{bi}}=a{\stackrel{\&OverBar;}{y}}_a+{b\stackrel{\&OverBar;}{y}}_b---\left(13\right)$

所以，我们得到以下结论：

定理1线性插值时，***的PDM中间模型平均形状是相邻模型平均形状的线性组合，其参数由式(13)确定。

(2)中间模型的形变

相似的，中间PDM模型的形变可以这样算出：

${\mathrm{\&Delta;}}_{\mathrm{yti}}=y_{\mathrm{ti}}-{\stackrel{\&OverBar;}{y}}_t=({\mathrm{ay}}_{\mathrm{ai}}+{\mathrm{by}}_{\mathrm{bi}})-(a{\stackrel{\&OverBar;}{y}}_a+b{\stackrel{\&OverBar;}{y}}_b)$

$=a(y_{\mathrm{ai}}-{\stackrel{\&OverBar;}{y}}_a)+b(y_{\mathrm{bi}}-{\stackrel{\&OverBar;}{y}}_b)---\left(14\right)$

$=a{\mathrm{\&Delta;}}_{\mathrm{yai}}+b{\mathrm{\&Delta;}}_{\mathrm{ybi}}$

从上式可以得到以下结论：

定理2线性插值时，***的PDM中间模型的形状变化速度是相邻模型形变速度的线性组合，其参数由式(14)确定。

(3)中间模型的协方差矩阵

中间模型的协方差矩阵可以由下式推导得出

$S_t=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}(y_{\mathrm{ti}}-{\stackrel{\&OverBar;}{y}}_t){(y_{\mathrm{ti}}-{\stackrel{\&OverBar;}{y}}_t)}^T=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yti}}{\mathrm{\&Delta;}}_{\mathrm{yti}}^T$

$=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}({\mathrm{a\&Delta;}}_{\mathrm{yai}}+{\mathrm{b\&Delta;}}_{\mathrm{ybi}}){({\mathrm{a\&Delta;}}_{\mathrm{yai}}+{\mathrm{b\&Delta;}}_{\mathrm{ybi}})}^T---\left(15\right)$

$=\frac{a^2}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yai}}{\mathrm{\&Delta;}}_{\mathrm{yai}}^T+\frac{b^2}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}^T+\frac{\mathrm{ab}}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yai}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}^T+\frac{\mathrm{ab}}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}{\mathrm{\&Delta;}}_{\mathrm{yai}}^T$

令

$S_{\mathrm{ab}}=\frac{\mathrm{ab}}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yai}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}^T---\left(16\right)$

因为

$S_a=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}(y_{\mathrm{ai}}-{\stackrel{\&OverBar;}{y}}_a){(y_{\mathrm{ai}}-{\stackrel{\&OverBar;}{y}}_a)}^T=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yai}}{\mathrm{\&Delta;}}_{\mathrm{yai}}^T---\left(17\right)$

$S_b=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}(y_{\mathrm{bi}}-{\stackrel{\&OverBar;}{y}}_b){(y_{\mathrm{bi}}-{\stackrel{\&OverBar;}{y}}_b)}^T=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}^T---\left(18\right)$

显然有

$S_{\mathrm{ab}}=S_{\mathrm{ba}}^T---\left(19\right)$

最终，得到以下推论。

推论1线性插值时，中间模型的协方差矩阵可由以下表达式得到

$S_t=a^2{S}_a+b^2{S}_b+\mathrm{ab}{S}_{\mathrm{ab}}+{\mathrm{abS}}_{\mathrm{ab}}^T---\left(20\right)$

这里S_ab表示模型M_a和M_b的交变矩阵。我们知道协方差矩阵是通过对同一时间段内不同个体的采样数据进行统计分析得出的，而交变矩阵则描述了在不同采样时间得到的两个对象模型的形变过程。

本发明的有益效果主要表现在：能够形成实时分布的中间过渡模型、模型精度高、实用性良好。

附图说明

图1是从时间序列中的临近模型插值到得的瞬时PDM模型的示意图。

图2是典型的心动周期分为六个阶段的示意图。

图3是心脏不同阶段特征的六个PDM模型的示意图。

图4是插值生成新的PDM模型以进行柔性体分析的示意图。

图5是通过插值生成的PDM心脏模型的时间平滑变化的示意图。

具体实施方式

下面结合附图对本发明作进一步描述。

参照图1～图5，一种基于时变插值算法的柔性体点分布模型的建立方法，包括以下步骤：

1)、载入两个相邻时刻T_a和T_b的柔性体分布模型M_a和M_b，其计算式为：

$M_a=<{\stackrel{\&OverBar;}{y}}_a,{\mathrm{\&Phi;}}_a,{\mathrm{\&lambda;}}_a>,$ $M_b=<{\stackrel{\&OverBar;}{y}}_b,{\mathrm{\&Phi;}}_b,{\mathrm{\&lambda;}}_b>;$

2)、通过所述柔性体分布模型M_a和M_b，分别计算平均形状

和

特征矩阵Φa和Φb、以及特征值λ_a和λ_b，设定柔性体从T_a到T_b呈连续变化，

2.1)中间模型的平均形状

用线性插值得到y_ai和y_bi之间的中间模型形状y_ti为：

y_ti＝ay_ai+by_bi                        (12)

a和b由式(10)和(11)定义：

$a=\frac{t-T_a}{T_b-T_a}---\left(10\right)$

b＝1-a                                (11)

a和b反映了中间模型相对于M_a和M_b的距离；

***的平均形状为：

${\stackrel{\&OverBar;}{y}}_t=\frac{1}{s}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{y}_{\mathrm{ti}}=\frac{1}{s}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}({\mathrm{ay}}_{\mathrm{ai}}+{\mathrm{by}}_{\mathrm{bi}})=\frac{a}{s}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{y}_{\mathrm{ai}}+\frac{b}{s}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{y}_{\mathrm{bi}}=a{\stackrel{\&OverBar;}{y}}_a+{b\stackrel{\&OverBar;}{y}}_b---\left(13\right);$

2.2)、中间模型的形变

中间模型的形变由下式算出：

${\mathrm{\&Delta;}}_{\mathrm{yti}}=y_{\mathrm{ti}}-{\stackrel{\&OverBar;}{y}}_t=({\mathrm{ay}}_{\mathrm{ai}}+{\mathrm{by}}_{\mathrm{bi}})-(a{\stackrel{\&OverBar;}{y}}_a+b{\stackrel{\&OverBar;}{y}}_b)$

$=a(y_{\mathrm{ai}}-{\stackrel{\&OverBar;}{y}}_a)+b(y_{\mathrm{bi}}-{\stackrel{\&OverBar;}{y}}_b)---\left(14\right)$

$=a{\mathrm{\&Delta;}}_{\mathrm{yai}}+b{\mathrm{\&Delta;}}_{\mathrm{ybi}}$

2.3)、中间模型的协方差矩阵

中间模型的协方差矩阵由下式推导得出：

$S_t=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}(y_{\mathrm{ti}}-{\stackrel{\&OverBar;}{y}}_t){(y_{\mathrm{ti}}-{\stackrel{\&OverBar;}{y}}_t)}^T=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yti}}{\mathrm{\&Delta;}}_{\mathrm{yti}}^T$

$=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}({\mathrm{a\&Delta;}}_{\mathrm{yai}}+{\mathrm{b\&Delta;}}_{\mathrm{ybi}}){({\mathrm{a\&Delta;}}_{\mathrm{yai}}+{\mathrm{b\&Delta;}}_{\mathrm{ybi}})}^T---\left(15\right)$

$=\frac{a^2}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yai}}{\mathrm{\&Delta;}}_{\mathrm{yai}}^T+\frac{b^2}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}^T+\frac{\mathrm{ab}}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yai}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}^T+\frac{\mathrm{ab}}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}{\mathrm{\&Delta;}}_{\mathrm{yai}}^T$

令

$S_{\mathrm{ab}}=\frac{\mathrm{ab}}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yai}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}^T---\left(16\right)$

因为

$S_a=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}(y_{\mathrm{ai}}-{\stackrel{\&OverBar;}{y}}_a){(y_{\mathrm{ai}}-{\stackrel{\&OverBar;}{y}}_a)}^T=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yai}}{\mathrm{\&Delta;}}_{\mathrm{yai}}^T---\left(17\right)$

$S_b=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}(y_{\mathrm{bi}}-{\stackrel{\&OverBar;}{y}}_b){(y_{\mathrm{bi}}-{\stackrel{\&OverBar;}{y}}_b)}^T=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}^T---\left(18\right)$

显然有

$S_{\mathrm{ab}}=S_{\mathrm{ba}}^T---\left(19\right);$

3)、根据得到中间模型的平均形状、形变和协方差矩阵，得到t时刻的中间模型Mt。

在我们一个关于心脏图像分析的研究项目中，可以根据心动周期定义一个规范化的时间坐标轴。例如，一个心跳周期可以分为六个生理阶段：心房收缩、心室收缩、心室射血、心室舒张、心室快速填充及舒张后期休息(图2)。也有人进一步将心室射血的过程分为射血前期和射血后期。左心与右心的运动过程相似但不完全相同。主要的差别是，左心室与动脉的收缩压比右心的大三至四倍，而且左右心泵血过程的时间方面也有显著区别。由于整个心脏的空间形态在一个心动周期内变化很大，而单一的模型很难准确描述动态的心脏运动过程，所以我们需要建立一组模型来描述心动不同阶段的心脏形态。

在gated-SPECT的图像研究中，可参考心电图信号进行时间规范化^[19]。我们希望通过对心脏形状进行统计后可以得到全部的六个模型，并且这些模型能用来准确地描述心脏的特定状态，如心脏舒张末期和心脏收缩末期。图3在标准心电图曲线中标出了各个阶段模型出现的位置。图中我们已经把压力值幅度经过了修正，使得波峰处的压力值为正常人体的收缩压(120mm汞柱)，所以图中也表示一条典型的左心室压力曲线。

如果在正常的心动周期中有6个PDM模型，则式(6)中的Ti(代表某一相位所在的时刻)可以定义为：

式中t是获得gated-SPECT图像的时间值，H_c是心率的的倒数，

表示不大于某小数的最大整数。

然而，在实际的例子中，常会产生不同数目的相位(如图4中有8个)，这是由于心脏图像的采样通常有固定的频率，如每100ms(图4)。那么就有必要通过插值得到采样时刻的PDM模型，以期得到更好的三维心脏形状的分割。

柔性体点分布模型插值算法，作者研究过程中分别采用计算机仿真及真实心脏数据进行了大量的实验。在仿真中，我们可以假设一组九维(或更高)的模型，yi＝[x_1i，x_2i，...，x_9i]^T。两个PDM模型M1和M2的数据集都是随机产生的。各组数据都含有100个样本，作为生成点分布模型的训练集，例如

${\stackrel{\&OverBar;}{y}}_1=\frac{1}{100}\underset{i=1}{\overset{100}{\mathrm{\&Sigma;}}}{y}_{1i},$ $S_1=\frac{1}{99}\underset{i=1}{\overset{100}{\mathrm{\&Sigma;}}}{(y_{1i}-{\stackrel{\&OverBar;}{y}}_1)(y_{1i}-{\stackrel{\&OverBar;}{y}}_1)}^T$

为生成一个***模型，我们从S₁计算M1的特征值及特征向量，对S₂也同样。假设***模型的左右距离为a＝0.3，b＝0.7。为了进行对比，我们还计算它们的交变矩阵S_ab和相应的特征向量组。实验结果显示，本中所述的方法能很好的计算出***模型的各种参数。由式(31)计算出的特征值误差小于1.83％，由式(35)计算出的特征值误差小于0.05％。该仿真实验的C++程序已由所在实验室开发完成。

用真实的心脏数据进行了测试。通过灌注SPECT成像法得到的大量图像数据中，我们预先建立了左心室在不同时刻的五个PDM模型。每个模型由246组病例的真实图像的训练得到(图像主要由巴塞罗那医院提供)，包括一个平均形状的VTK表达，一组特征值矢量和特征向量。所有医学图像的处理过程和模型的***运算也由都C++程序实现(研究组共开发了约600MB程序代码)。实验时，我们在每两个相邻的PDM之间均匀地***四个中间模型，得到的实验结果令人非常满意。图5显示通过插值生成的PDM模型随时间能够平滑地变化左心室柔性模型。

点分布模型作为一种描述柔性对象的统计学模型，已被证明在很多应用上是非常有效的，尤其是在医学图像分析领域。为了分析柔性体的动态变化情况，本文首创研究了点分布模型的插值算法。研究结果表明，可由PDM的平均形状、特征值和特征向量等参数生成任意时刻的中间模型。事实上，***模型的平均形状及形变就是它相邻前后模型对应项的线性组合，即

Δ_yti＝aΔ_yai+bΔ_ybi。其特征值则可以采用式λ_t≈a²λ_a+b²λ_b计算得到。最后，为求出***模型的协方差矩阵和特征向量，本文提出了交变矩阵形式并给出了相应的计算方法。论文还分别采用计算机仿真和真实心脏数据进行了大量的实验验证。

Claims (1)

1.一种基于时变插值算法的左心室柔性体点分布模型建立方法，其特征在于：所述建立方法包括以下步骤：

1)、载入两个相邻时刻T_a和T_b的左心室柔性体分布模型M_a和M_b，其计算式为：

$M_a=<{\stackrel{\&OverBar;}{y}}_a,{\mathrm{\&Phi;}}_a,{\mathrm{\&lambda;}}_a>,$ $M_b=<{\stackrel{\&OverBar;}{y}}_b,{\mathrm{\&Phi;}}_b,{\mathrm{\&lambda;}}_b>;$

2)、通过所述柔性体分布模型M_a和M_b，分别计算平均形状

和特征矩阵Φa和Φb、以及特征值λ_a和λ_b，设定柔性体从T_a到T_b呈连续变化，

2.1)中间模型的平均形状

用线性插值得到y_ai和y_bi之间的中间模型形状y_ti为：

y_ti＝ay_ai+by_bi    (12)

a和b由式(10)和(11)定义：

$a=\frac{t-T_a}{T_b-T_a}---\left(10\right)$

b＝1-a    (11)

a和b反映了中间模型相对于M_a和M_b的距离；

***的平均形状为：

${\stackrel{\&OverBar;}{y}}_t=\frac{1}{s}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{y}_{\mathrm{ti}}=\frac{1}{s}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}({\mathrm{ay}}_{\mathrm{ai}}+{\mathrm{by}}_{\mathrm{bi}})=\frac{a}{s}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{y}_{\mathrm{ai}}+\frac{b}{s}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{y}_{\mathrm{bi}}=a{\stackrel{\&OverBar;}{y}}_a+b{\stackrel{\&OverBar;}{y}}_b---\left(13\right);$

2.2)、中间模型的形变

中间模型的形变由下式算出：

${\mathrm{\&Delta;}}_{\mathrm{yti}}=y_{\mathrm{ti}}-{\stackrel{\&OverBar;}{y}}_t=({\mathrm{ay}}_{\mathrm{ai}}+{\mathrm{by}}_{\mathrm{bi}})-(a{\stackrel{\&OverBar;}{y}}_a+b{\stackrel{\&OverBar;}{y}}_b)$

$=a(y_{\mathrm{ai}}-{\stackrel{\&OverBar;}{y}}_a)+b(y_{\mathrm{bi}}-{\stackrel{\&OverBar;}{y}}_b)---\left(14\right)$

$=a{\mathrm{\&Delta;}}_{\mathrm{yai}}+b{\mathrm{\&Delta;}}_{\mathrm{ybi}}$

2.3)、中间模型的协方差矩阵

中间模型的协方差矩阵由下式推导得出：

$S_t=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}(y_{\mathrm{ti}}-{\stackrel{\&OverBar;}{y}}_t){(y_{\mathrm{ti}}-{\stackrel{\&OverBar;}{y}}_t)}^T=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yti}}{\mathrm{\&Delta;}}_{\mathrm{yti}}^T$

$=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{(a{\mathrm{\&Delta;}}_{\mathrm{yai}}+b{\mathrm{\&Delta;}}_{\mathrm{ybi}})(a{\mathrm{\&Delta;}}_{\mathrm{yai}}+b{\mathrm{\&Delta;}}_{\mathrm{ybi}})}^T---\left(15\right)$

$=\frac{a^2}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yai}}{\mathrm{\&Delta;}}_{\mathrm{yai}}^T+\frac{b^2}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}^T+\frac{\mathrm{ab}}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yai}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}^T+\frac{\mathrm{ab}}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}{\mathrm{\&Delta;}}_{\mathrm{yai}}^T$

令

$S_{\mathrm{ab}}=\frac{\mathrm{ab}}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yai}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}^T---\left(16\right)$

因为

$S_a=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}(y_{\mathrm{ai}}-{\stackrel{\&OverBar;}{y}}_a){(y_{\mathrm{ai}}-{\stackrel{\&OverBar;}{y}}_a)}^T=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{yai}}{\mathrm{\&Delta;}}_{\mathrm{yai}}^T---\left(17\right)$

$S_b=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}(y_{\mathrm{bi}}-{\stackrel{\&OverBar;}{y}}_b){(y_{\mathrm{bi}}-{\stackrel{\&OverBar;}{y}}_b)}^T=\frac{1}{s-1}\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}{\mathrm{\&Delta;}}_{\mathrm{ybi}}^T---\left(18\right)$

显然有

$S_{\mathrm{ab}}=S_{\mathrm{ba}}^T---\left(19\right);$

3)、根据得到中间模型的平均形状、形变和协方差矩阵，得到t时刻的中间模型Mt，所述中间模型Mt的特征值采用式λ_t≈a²λ_a+b²λ_b计算得到，中间模型Mt的特征向量Φ_t则包含了一组由协方差矩阵S_t分解得到的特征向量。

Images