专利名称：基于插值算法的柔性体点分布模型的特征值和特征向量计算方法
技术领域：
本发明涉及柔性体点分布模型分析技术，结合图形学、矩阵理论、线性代数、医学图像分析、计算机图像处理技术和柔性体模型的参数计算方法。
背景技术：
柔性体相对于刚性体，在计算机中进行建模和分析比较复杂，而在现实生活中物体随时间柔性变化的情况是非常普遍的，如文献1王敏琴 韩国强 涂泳秋.基于形变模型的医学图像分割综述[J].医疗卫生装备，2009，30(2)-37-39；文献2“Munsell，B.C.；Dalal，P.；Song Wang，″Evaluating Shape Correspondence forStatistical Shape AnalysisA Benchmark Study，″IEEE Transactions on PatternAnalysis and Machine Intelligence，vol.30，no.11，pp.2023-2039，Nov.2008，即马塞尔，B.C；戴拉，P.；王松，”统计形状分析中的形状匹配估计一种校准研究”，IEEE模式分析与机器智能会刊，vol.30，no.11，2008(11)2023-2039；以及文献3龚永义 罗笑南 贾维嘉 黄厚生.基于改进的弹簧质子模型的医学图像配准.计算机学报，2008(7)1224-1233。对于许多生物医学组织，如人的心脏，双手，骨骼，肾脏、大脑、细胞等等，这些人们可以很好地描述其外部形状特征但却很难给出准确的刚性模型，如文献4陈昱 庄天戈 王合.直方图估计互信息在非刚性图像配准中的应用.计算机学报，2000(4)444-448；文献5简江涛 冯焕清 熊进.点分布模型约束的主动轮廓及其在脑MR图像分割中的应用.中国生物医学工程学报，2006，25(5)513-517；以及文献6周寿军，梁斌，陈武凡.心脏序列图像运动估计新方法基于广义模糊梯度矢量流场的形变曲线运动估计与跟踪.计算机学报.2003.26(11)。点分布模型(PDM)是一种新近崛起的技术，它能很好的对研究对象进行描述，可以成功地应用于非刚性模型的分解与合成处理。近年来，柔性体的统计学建模及其与图像解释技术的结合催生了一些高级应用软件，这些软件不仅针对生物医学，还涉及目标跟踪、识别、以至影视等领域，如文献7唐果 赵晓东 汪元美.三维医学图像分割与可视化研究.计算机学报，1998(3)204-209；文献8段侪杰 马竟锋 张艺宝 侯凯 包尚联.能量传导模型及在医学图像分割中的应用[J].软件学报，2009(5)-1106-1115；文献9侯立华.改进Snake模型在医学超声图像分割中的应用[J].计算机仿真，2008，25(8)-183-185，322；以及文献10闵小平 王博亮 戴培山 黄晓阳.基于主动形变模型的医学序列图像的分割[J].系统仿真学报，2007，19(22)-5331-5335。
PDM是一种基于对象样本分析的统计模型。从Cootes等人开始，这种基于训练样本上标记点的形体建模思想一直被贯彻下来，如文献11Timothy F.Cootes，Christopher J.Taylor，“Combining point distribution models with shape models basedon finite element analysis”，Image Vision Comput.13(5)403-409(1995)，即蒂姆F.古茨，克里斯托弗J.泰勒，“有限元分析中点分布模型与形状模型的结合”，计算机图像与视觉，1995，13(5)403-409。最典型的想法是，使标记点自动对齐以确定图形的位置和形变类型。对齐图像后可以很容易的比较同组标记点在不同例图中的位置变化，只需要计算它们的坐标就可以了。因而，PDM模型对于形状的描述是通过轮廓点向量来完成的，它由中间形状和形状变化样式组成。通过对这种分布情况建模，我们可以得到与原始训练集相似的新形状。
在三维空间，如果训练集图像中的每个图形都由一组数量为n的标记点描述，那么它可以用一个3n维的矢量表示(如果是二维对象即为2n维矢量)。训练集中的标记工作非常重要，通常是手动完成的。另外，在研究中我们要收集数百名例柔性体的个例以得到比较好的统计结果，如文献12Guoyan Zheng；XiaoDong；Rajamani，K.T.；Xuan Zhang；Styner，M.；Thoranaghatte，R.U.；Nolte，L.-P.；Ballester，M.A.G.，″Accurate and Robust Reconstruction of a Surface Model of theProximal Femur From Sparse-Point Data and a Dense-Point Distribution Model forSurgical Navigation，″IEEE Transactions on Biomedical Engineering，vol.54，no.12，pp.2109-2122，Dec.2007，即张国艳，董晓，阿雅曼尼，K，T；张璇，斯蒂那，M；索拉纳哈迪，R U.；诺蒂，L，P.；巴耶斯特尔，M.A.G.，“手术导航系统中基于稀疏点数据及密集点分布模型的股骨近端精确鲁棒重建”，IEEE生物医学工程会刊，vol.54，no.12，2007(12)2109-2122；文献13Koikkalainen，J.；Tolli，T.；Lauerma，K.；Antila，K.；Mattila，E.；Lilja，M.；Lotjonen，J.，″Methods of Artificial Enlargementof the Training Set for Statistical Shape Models，″IEEE Transactions on MedicalImaging，vol.27，no.11，pp.1643-1654，Nov.2008，即考克克莱纳，J.；托里，T.；劳尔玛，K.；按提拉，K.；马蒂拉，E.；莉亚，M.；罗提那，J.，“统计形状模型训练集人工扩大方法”，IEEE医学图像会刊，vol.27，2008(11)1643-1654。训练集形成的概率模型维度可能非常的大(在一个心脏模型中通常可达到上万维)，所以建立PDM的过程是非常枯燥又费时费力的。然而，在一些应用中，单独的一个PDM模型或是有限个PDM模型都不足以很好的描述对象变化。事实上，不论是连续PDM模型还是模型插入方法都是为了满足对象的性质，即对象时间和空间上的变化。例如在心脏建模时，需要不同时间段的模型才能很好地描述心室形状在心脏跳动不同阶段的变化。但是，在构建整个心脏运动周期的连续模型时，我们只能由已知的模型序列得到特定几个时刻的模型，而其他任意时刻的模型就不得而知了，尽管在基于模型的图像分割领域PDM建模出现已有十来年，但对于模型插值及相关内容的深入研究还没有开展过。
发明内容
为了克服已有的柔性体点分布模型的不能形成实时分布的模型、模型精度差、实用性差的不足，本发明提供一种能够形成实时分布的中间过渡模型、模型精度高、实用性良好的基于插值算法的柔性体点分布模型的特征值和特征向量计算方法。
本发明解决其技术问题所采用的技术方案是 一种基于插值算法的柔性体点分布模型的特征值和特征向量的计算方法，所述计算方法包括以下步骤 1)、载入两个相邻时刻Ta和Tb的柔性体分布模型Ma和Mb，其计算式为 Ma＝<ya，Φa，λa>，Mb＝<yb，Φb，λb>； 2)、通过所述柔性体分布模型Ma和Mb，分别计算平均形状ya和yb、特征矩阵Φa和Φb、以及特征值λa和λb，设定柔性体从Ta到Tb呈连续变化， 2.1)中间模型的特征值的计算 λt≈a2λa+b2λb，(25) 这里λa和λb分别为前后两个相邻PDM模型的特征值向量，即 λa＝(λa1，λa2，λa3，...)，(26) λb＝(λb1，λb2，λb3，...). (27) 2.2)、中间模型的特征向量的估计 插入中间模型的协方差矩阵近似为 St≈a2Sa+b2Sb (28) 其中，协方差矩阵Sa和Sb由特征矩阵Φa和Φb近似得到，Φa和Φb从插补前已知的点分布模型中得到；设定中间模型建立前先对训练集数据进行PCA分析以减少参数个数，并选取t种形变模式来描述整个对象的形变；由于Φ是协方差矩阵中最前面的t个特征矢量，则有 其中，Λa和Λb是由特征矢量λa和λb组合形成的对角矩阵；用奇异值分解的方法从式得到插入中间模型的特征矢量 由两个已知的相邻PDM模型及相关参数a、b、Λa、Λb、Φa和Φb确定正交矩阵Ut，所述正交矩阵Ut的前t行就是特征矩阵Φt。
进一步，所述步骤2)中还包括 2.3)、中间模型的协方差矩阵由下式推导得出 令 因为 显然有 根据协方差矩阵，矩阵St由Sa、Sb和Sab计算得到，通过如下公式 用奇异值分解法，可令 通常Φab≠Φba，因为矩阵Sab并非对称矩阵，Φab和Φba都包含t个矢量； 在插补时，点分布模型保留特征值和特征向量，中间模型的协方差矩阵通过下式计算 再进一步，所述步骤2)中还包括 2.3)、特征向量组的计算 通过式(31)来计算特征向量时，将(31)改写为 其中是一个r×N维的矩阵，其中r是主成分数目，N是样本维数，即标定点的数量与标记点维度的乘积； 如果Ui是矩阵G＝LΦa的主要特征向量，对应特征值为λi，则(Φaui)和λi分别为矩阵S＝ΦaL的特征向量和特征值。
由Gui＝λiui，，得 S(Φaui)＝ΦaL(Φaui)＝ΦaGui＝Φa(λiui)＝λi(Φaui)(36)。
本发明的技术构思为PDM方法是一种先验建模方法，它建立在目标对象的一组样本(即训练集)已知的基础上。这组样本集可以提供对象的形状及形变的统计学描述。在实际应用中，这组对象样本的形状描述通常是使用轮廓来进行的(如图片上的一组像素的坐标值)。进一步的，可以在轮廓上选取一些标记点来描述目标物形状，而这些标记点往往与目标物的某种特征相关联。
不同样本上的对应标记点在选取的时候遵循着“特征位置大致相同”的原则，但往往这种定位是不够精确的，存在一定偏差。这种偏差要归因于不同个体间的自然形变。我们期望相对于整个形体的尺度而言，这种偏差是十分微小的。PDM方法允许我们模仿这种微小的差异，并且指出哪些偏差的确是微小的，哪些相对而言更大一些。
由于训练集中的形体是从不同的图片集中获得的，而不同图片集的坐标系并不统一，因此有必要先将所有训练集形体进行坐标归一化，即将它们对齐。对齐的过程实际上就是对每一个样本模型寻找适当的相似性变换的过程，包括平移、缩放和旋转，以使得不同样本之间尽可能的接近，另一方面，选定的变换也应该令各个样本模型与由所有样本得到的平均模型尽可能接近。为了方便表述，假设我们的样本只有两个，归一化简化为对齐这两个样本。每个样本的形状都由一个标记点坐标组成的向量来表示 x＝(x1，y1，z1，...，xn，yn，zn)T， (1) 其中(xi，yi，zi)是某个标记点在三维图像中的空间坐标。如果是二维医学图像，则标记点用(xi，yi)表示。
图片上标记点的获取有多种方式，人工手动标记或采用一定得自动标记算法都是可以的。例如，Izard等人提出的一种MRI图像标记方法，使用统一的算法规则在不同人脑图片中进行组织结构标记，如文献14C.Izard，B.Jedynak，and C.Stark，Automatic Landmarking of Magnetic Resonance brain Images，Proc.of theSPIE International Symposium on Medical Imaging，12-17 February 2005，San Diego，California，USA，即C.伊扎德，B.扎第纳克，C.斯塔克，“大脑核磁共振图像的自动标记”，SPIE医学图像网络研讨会，2005年2月12-17日，美国，加利福尼亚，圣地亚哥。得到一组经过标记的训练样本后，我们需要将它们对齐到一个统一的坐标系下。可以使用广义对齐算法(GPA)来将训练样本对齐，并使得所有样本形状到平均模型的距离平方和最小，实际上也就是寻找变换Ti(平移，旋转，缩放)，使下式最小化 Dm＝∑|m-Ti(xi)|2， (2) 其中 且xi是训练集中的一个形状样本，它是一个3n维的向量表示(在三维空间中有n个标记点的情况下)。
对齐后的训练集组成了一个三维空间上的点云，可以看作是一个概率密度函数。为了降低运算成本和内存占用，我们采用主成分分析法(PCA)寻找点云中的主元，并以前列最大的少数分量建立模型。这些被选出的分量可以描述对象主要的形变。以心动周期内左心室的模型为例，我们经过主成分分析后的用60个特征向量已经能足够描述心室的形状及其变化了。最后，PDM模型可以表达为 x＝x+p1b1+…+ptbt＝x+Φc(4) 其中，x是对齐后的平均形状向量，Φ是一个3n×t的矩阵，Φ的每一列表示主成分方向上的单位向量，c是由形状参数组成的t维列向量。可见，经过PCA分析后形状的维度已经由3n下降到了t。
经过以上步骤我们得到了一个类似于点分布模型的统计学模型，该模型可应用于柔性体的物体建模，以及在新的图像中定位和分割新个例。在训练集的学习过程中我们可以得到形状参数b的变化范围，只要在该范围内变化，任意b的值都可以生成合理的新样本(仿真形体)。通常bi的变化为λj(矩阵Φ中对应第i大的特征值)或者
从样本集中训练出一个统计学形态模型后，就可以用PDM来解释新的图像个例了。为了将模型与图像匹配，通常采用一种主动重复迭代的方法。该方法使模型不断地变形来适应新图像中目标对象的轮廓。模型形状的变化受到一些统计数据的约束，也即只能在训练样本集中所得出的范围内进行变化。对于形状模型，我们还要求图像轮廓大致位于模型点的附近。可以表示为某点的梯度与轮廓上过该点的法线方向的变化。令统计模型的轮廓沿图像轮廓不断移动，同时计算两轮廓的马氏距离最小化就可以最终得到真实的边界位置。因此，要在特定实例中获得目标柔性对象的形状需要反复执行以下两步直到收敛(1)沿模型各点的法线寻找图像轮廓上相匹配的对应点(使得模型点与对应图像点的马氏距离最小)；(2)改变模型位置及形状参数使得模型各标记点更接近图像轮廓上找到的对应点。也就是说，模型的拟合过程实际上就是在图像上搜寻模型点的最相似位置及不断修正整体几何变换T及形状参数c的过程，使之最小化。数学上可以表示为 f＝|X-T(x+Pc)|2＝f(b，Xc，Yc，Zc，s，α，β，γ)(5) ＝|X-T(x+Pc；Xc，Yc，Zc，s，α，β，γ)|2 其中X是在当前迭代时得到的临时模型。
这是一个寻找最小值的问题，可以通过一些非线性优化的迭代方法求解。最终我们可以确定整体变换参数T，以及个体实例的形状参数b c＝PT(T-1(X)-x). (6) 近几年，人们已经就基于PDM的柔性体建模方法的许多相关问题进行了深入的研究，如形体对齐、自动标记、平均形体生成、形变建模、模态分析、图像分割等等问题，如文献15蒋晓悦，赵荣椿，一种改进的活动轮廓图像分割技术，中国图象图形学报A辑，2004，9(9)1019-1024；文献16马峰，唐泽圣，夏绍玮.多尺度几何活动曲线及MR图像边界提取.计算机学报.2000.23(8)；文献17Gilliam，A.D.；Epstein，F.H.；Acton，S.T.，″Cardiac Motion Recovery via ActiveTrajectory Field Models，″IEEE Transactions on Information Technology inBiomedicine，vol.13，no.2，pp.226-235，March 2009，即吉列姆，A.D.；爱泼斯坦，F.H.；埃克森，S.T.，“基于轨迹场模型的心动复原”，IEEE生物医学信息技术会刊，vol.13，2009(3)226-235。此外，在该思想广泛地应用于非刚体对象建模的基础上，并可进一步扩展为诸如活动形体模型(ASM)和活动表现建模(AAM)等以应用于更深入的图像分析。目前，对PDM模型最成功的表达是通过PCA方法得到的，如文献18Tobon-Gomez，C.；Butakoff，C.；Aguade，S.；Sukno，F.；Moragas，G.；Frangi，A.F.，″Automatic Construction of 3D-ASM Intensity Models by Simulating ImageAcquisitionApplication to Myocardial Gated SPECT Studies，″IEEE Transactions onMedical Imaging，vol.27，no.11，pp.1655-1667，Nov.2008，即托顿-哥麦茨，C.；巴塔考夫，C.；阿古爱德，S.；萨克诺，F.；莫拉盖斯，G.；弗兰基，A.F.，“基于模拟图像采集的三维ASM密度模型自动重建应用于心脏门闸SPECT图像研究”，IEEE医学图像会刊，vol.27，2008(11)1655-1667。PCA方法非常简单易懂，如文献19Engelen，S.，Hubert，M.，Vanden Branden，K.，A comparison of three procedures forrobust PCA in high dimensions，Workshop on Robustness for High-dimensional Data(RobHD2004)，2004，即英格伦，S.，赫伯特，M.，范登布兰登，K.，“高维度PCA鲁棒分析的三个过程对比”，高维数据鲁棒性学术讨论会，2004，首先把一个形体投影到经过学习PCA空间(由线性无关方向形成的正交空间)，然后从训练样本中提取少量与形变模式相关的参数(称为主要成分)，再用这些参数来控制形体的变形。显然，这种形变是沿着样本集所确定的最大形变方向进行的。由于只研究几个主要方向上的形变，无论是计算复杂度还是对存储空间的要求都大大降低了。从数学角度讲，首先是从训练集中估计得到平均形状y和形体协方差矩阵S 这样，一个合理的形状可以被表达为 y＝y+Φc (8) 式中向量c由对应于各个特征向量的形变参数组成。矩阵Φ则包含了一组由协方差矩阵S分解得到的特征向量，这组向量可线性组合成任意的新形状y。根据PCA理论，c中的每个参数都控制形体在某一方向上的形变，并且，各个方向都彼此独立，参数按照形体在对应方向上形变由大到小顺序排列。本发明考虑在一个PDM模型仅包含平均形状、特征值和特征向量而不包括原始训练集的数据的情况下进行插值运算。即 M＝<y，Φ，λ>(9) 如上所述，单个PDM模型可以由一组柔性对象的形体获得，例如，我们研究的医院患者的心脏图片(一些被观测的柔性体个例)。然而，有时我们要观测这些柔性体随时间的变化情况，这样可以得到一系列的形状信息，例如心脏在整个活动周期内的三维变形。而医院现有设备只能获得若干时刻的三维图像，即对应于某些Ti时刻的模型Mi。
对一个特定的动态形变物体，假设我们已经在离散的时间坐标(T1，T2，...)上建立起一组PDM模型(M1，M2，...)。我们可以把问题描述为给出已知的两个PDM模型，Ta时刻的模型Ma和Tb时刻的模型Mb(模型已知意味着我们知道了平均形状ya和yb、特征矩阵Φa和Φb、以及特征值λa和λb)，如何得到在t(Ta≤t≤Tb)时刻的中间模型Mt的相关参数。这些参数有平均形状、特征值和特征向量，即Mt＝<yt，Φt，λt>(图1)。这里我们假设在Ta到Tb的时间段中，一个形体从Yai连续变化到Ybi，二者都是n维的，如Yai＝(xai1，xai2，xai3，...，xain)T。这里还定义了两个常数 b＝1-a(11) a和b反映了中间模型相对于Ma和Mb的距离。
(1)中间模型的平均形状 根据以上的定义和假设，用线性插值(可满足大部分的实际应用)得到Yai和Ybi之间的中间模型形状为 yti＝ayai+bybi(12) a和b由式(10)和(11)定义，即a+b＝1。
插入的平均形状为 所以，我们得到以下结论 定理1线性插值时，插入的PDM中间模型平均形状是相邻模型平均形状的线性组合，其参数由式(13)确定。
(2)中间模型的形变 相似的，中间PDM模型的形变可以这样算出 Δyti＝yti-yt＝(ayai+bybi)-(aya+byb) ＝a(yai-ya)+b(ybi-yb) (14) ＝aΔyai+bΔybi 从上式可以得到以下结论 定理2线性插值时，插入的PDM中间模型的形状变化速度是相邻模型形变速度的线性组合，其参数由式(14)确定。
(3)中间模型的协方差矩阵 中间模型的协方差矩阵可以由下式推导得出 令 因为 显然有 最终，得到以下推论。
推论1线性插值时，中间模型的协方差矩阵可由以下表达式得到 这里Sab表示模型Ma和Mb的交变矩阵。我们知道协方差矩阵是通过对同一时间段内不同个体的采样数据进行统计分析得出的，而交变矩阵则描述了在不同采样时间得到的两个对象模型的形变过程。
中间模型的估算 (1)交变矩阵 为方便我们的讨论，假设模型是n维的，如 Δyai＝[da1i da2i ... dani]T (21) 则，交变矩阵为 通常情况下，样本是正态或平均分布的，交变矩阵中的各元素趋于零，即 所以，交变矩阵对角线元素远小于矩阵Sa和Sb中的对应元素，即 (2)特征值的估算 可以看出，交变矩阵Sab和Sba对于中间模型特征值的计算影响非常小，因此我们可以估计 λt≈a2λa+b2λb， (25) 这里λa和λb分别为前后两个相邻PDM模型的特征值向量，即 λa＝(λa1，λa2，λa3，...)， (26) λb＝(λb1，λb2，λb3，...).(27) 对λt中元素估计的精度主要取决于样本集的数量和分布情况。经过计算机模拟发现，当样本集为30个随机均匀分布时，误差率可维持在百分之五以下。
(3)特征向量的估计方法 当样本集很大时，插入中间模型的协方差矩阵可以近似为 St≈a2Sa+b2Sb(28) 所以，协方差矩阵Sa和Sb可以由特征矩阵Φa和Φb近似得到，而Φa和Φb是可以从插补前已知的点分布模型中得到。假设PDM模型建立前先对训练集数据进行PCA分析以减少参数个数，并选取t种形变模式来描述整个对象的形变。由于Φ是协方差矩阵中最前面的t个特征向量，我们知道 其中Λa和Λb是由特征向量λa和λb组合形成的对角矩阵。然后，我们可以用奇异值分解(SVD)的方法从式得到插入中间模型的特征向量 现在，已由两个已知的相邻PDM模型及相关参数a、b、Λa、Λb、Φa和Φb，确定了正交矩阵Ut，且Ut的前t行就是特征矩阵Φt。
(4)特征向量的精确计算方法 通过上面的公式，我们可以由两个已知PDM模型计算得到插入中间模型的各个参数(如平均形状、特征值、特征向量等)。然而，由于我们忽略了交变矩阵的作用，因此这种方法存在一定的误差，尤其是当样本集容量很小或其分布比较异常时这种误差就更明显。事实上我们有一种精确的方法，就是直接生成训练集中所有相邻两个PDM模型的交变矩阵。根据式(20)，矩阵St可由Sa、Sb和Sab计算得到，这意味着矩阵Sab需要在构建Ma和Mb时就得到。实际上，可以通过如下公式 用奇异值分解法，可令 通常Φab≠Φba，这是因为矩阵Sab并不是一个对称矩阵。Φab和Φba都包含t个矢量。
在插补时，由于点分布模型只保留了特征值和特征向量的信息，中间模型的协方差矩阵应通过下式计算 相似的，中间模型的特征向量可以由式(31)得到。这种方法由于不受训练样本分布的影响，因此是一种更加准确的解决方案。
(5)特征向量组的计算 当直接通过式(31)来计算特征向量时，我们会发现由于矩阵通常很大，计算过程对存储空间和运算效率都有很高的要求。例如，我们实验中有一个三维心脏模型，包含2848个标记点，通过主成分分析后选择了90种形变模式，那么该模型的一个协方差矩阵将占据584MB的内存空间，连最简单的加减操作可能也要用数分钟时间。因此，为了简化计算复杂度，我们将(31)改写为 其中是一个r×N维的矩阵，其中r是主成分数目，N是样本维数，即标定点的数量与标记点维度的乘积。下面我们来描述简化后的计算方法。
定理3如果Ui是矩阵G＝LΦa的主要特征向量，对应特征值为λi，则(Φaui)和λi分别为矩阵S＝ΦaL的特征向量和特征值。
其证明相当简单，由Gui＝λiui，，得 S(Φaui)＝ΦaL(Φaui)＝ΦaGui＝Φa(λiui)＝λi(Φaui) (36) 本发明的有益效果主要表现在能够形成实时分布的中间过渡模型、模型精度高、实用性良好。
图1是从时间序列中的临近模型插值到得的瞬时PDM模型的示意图。
图2是典型的心动周期分为六个阶段的示意图。
图3是心脏不同阶段特征的六个PDM模型的示意图。
图4是插值生成新的PDM模型以进行柔性体分析的示意图。
图5是通过插值生成的PDM心脏模型的时间平滑变化的示意图。
具体实施例方式 下面结合附图对本发明作进一步描述。
参照图1～图5，一种基于插值算法的柔性体点分布模型的特征值和特征向量的计算方法，其特征在于所述计算方法包括以下步骤 1)、载入两个相邻时刻Ta和Tb的柔性体分布模型Ma和Mb，其计算式为 Ma＝<ya，Φa，λa>，Mb＝<yb，Φb，λb>； 2)、通过所述柔性体分布模型Ma和Mb，分别计算平均形状ya和yb、特征矩阵Φa和Φb、以及特征值λa和λb，设定柔性体从Ta到Tb呈连续变化， 2.1)中间模型的特征值的计算 λt≈a2λa+b2λb，(25) 这里λa和λb分别为前后两个相邻PDM模型的特征值向量，即 λa＝(λa1，λa2，λa3，...)，(26) λb＝(λb1，λb2，λb3，...). (27) 2.2)、中间模型的特征向量的估计 插入中间模型的协方差矩阵近似为 St≈a2Sa+b2Sb (28) 其中，协方差矩阵Sa和Sb由特征矩阵Φa和Φb近似得到，Φa和Φb从插补前已知的点分布模型中得到；设定中间模型建立前先对训练集数据进行PCA分析以减少参数个数，并选取t种形变模式来描述整个对象的形变；由于Φ是协方差矩阵中最前面的t个特征矢量，则有 其中，Λa和Λb是由特征矢量λa和λb组合形成的对角矩阵；用奇异值分解的方法从式得到插入中间模型的特征矢量 由两个已知的相邻PDM模型及相关参数a、b、Λa、Λb、Φa和Φb确定正交矩阵Ut，所述正交矩阵Ut的前t行就是特征矩阵Φt； 2.3)、中间模型的协方差矩阵由下式推导得出 令 因为 显然有 根据协方差矩阵，矩阵St由Sa、Sb和Sab计算得到，通过如下公式 用奇异值分解法，可令 通常Φab≠Φba，因为矩阵Sab并非对称矩阵，Φab和Φba都包含t个矢量； 在插补时，点分布模型保留特征值和特征向量，中间模型的协方差矩阵通过下式计算 2.3)、特征向量组的计算 通过式(31)来计算特征向量时，将(31)改写为 其中是一个r×N维的矩阵，其中r是主成分数目，N是样本维数，即标定点的数量与标记点维度的乘积； 如果Ui是矩阵G＝LΦa的主要特征向量，对应特征值为λi，则(Φaui)和λi分别为矩阵S＝ΦaL的特征向量和特征值。
由Gui＝λiui，，得 S(Φaui)＝ΦaL(Φaui)＝ΦaGui＝Φa(λiui)＝λi(Φaui) (36)。
在我们一个关于心脏图像分析的研究项目中，可以根据心动周期定义一个规范化的时间坐标轴。例如，一个心跳周期可以分为六个生理阶段心房收缩、心室收缩、心室射血、心室舒张、心室快速填充及舒张后期休息(图2)。也有人进一步将心室射血的过程分为射血前期和射血后期。左心与右心的运动过程相似但不完全相同。主要的差别是，左心室与动脉的收缩压比右心的大三至四倍，而且左右心泵血过程的时间方面也有显著区别，如文献20杨柳，汪天富，林江莉等，旋转扫描超声心脏图像的中轴匹配与图像插值，生物医学工程学杂志，2004，21(1)28-33，41；文献21陈强，周则明，王平安，夏德深，带标记线左心室MR图像的自动分割，中国图象图形学报A辑，2004，9(6)666-673；文献22罗渝兰，郑昌琼等，基于Snake模型的B超心脏图像的腔体分割及算法研究，四川师范大学学报自然科学版，2002，25(1)66-69。由于整个心脏的空间形态在一个心动周期内变化很大，而单一的模型很难准确描述动态的心脏运动过程，所以我们需要建立一组模型来描述心动不同阶段的心脏形态。
在gated-SPECT的图像研究中，可参考心电图信号进行时间规范化。我们希望通过对心脏形状进行统计后可以得到全部的六个模型，并且这些模型能用来准确地描述心脏的特定状态，如心脏舒张末期和心脏收缩末期。图3在标准心电图曲线中标出了各个阶段模型出现的位置。图中我们已经把压力值幅度经过了修正，使得波峰处的压力值为正常人体的收缩压(120mm汞柱)，所以图中也表示一条典型的左心室压力曲线。
如果在正常的心动周期中有6个PDM模型，则式(6)中的Ti(代表某一相位所在的时刻)可以定义为
式中t是获得gated-SPECT图像的时间值，Hc是心率的的倒数，
表示不大于某小数的最大整数。
然而，在实际的例子中，常会产生不同数目的相位(如图4中有8个)，这是由于心脏图像的采样通常有固定的频率，如每100ms(图4)。那么就有必要通过插值得到采样时刻的PDM模型，以期得到更好的三维心脏形状的分割。
柔性体点分布模型插值算法，作者研究过程中分别采用计算机仿真及真实心脏数据进行了大量的实验。在仿真中，我们可以假设一组九维(或更高)的模型，yi＝[x1i，x2i，...，x9i]T。两个PDM模型M1和M2的数据集都是随机产生的。各组数据都含有100个样本，作为生成点分布模型的训练集，例如 为生成一个插入模型，我们从S1计算M1的特征值及特征向量，对S2也同样。假设插入模型的左右距离为a＝0.3，b＝0.7。为了进行对比，我们还计算它们的交变矩阵Sab和相应的特征向量组。实验结果显示，本中所述的方法能很好的计算出插入模型的各种参数。由式(31)计算出的特征值误差小于1.83％，由式(35)计算出的特征值误差小于0.05％。该仿真实验的C++程序已由所在实验室开发完成。
用真实的心脏数据进行了测试。通过灌注SPECT成像法得到的大量图像数据中，我们预先建立了左心室在不同时刻的五个PDM模型。每个模型由246组病例的真实图像的训练得到(图像主要由巴塞罗那医院提供)，包括一个平均形状的VTK表达，一组特征值矢量和特征向量。所有医学图像的处理过程和模型的插入运算也由都C++程序实现(研究组共开发了约600MB程序代码)。实验时，我们在每两个相邻的PDM之间均匀地插入四个中间模型，得到的实验结果令人非常满意。图5显示通过插值生成的PDM模型随时间能够平滑地变化左心室柔性模型。
点分布模型作为一种描述柔性对象的统计学模型，已被证明在很多应用上是非常有效的，尤其是在医学图像分析领域。为了分析柔性体的动态变化情况，本发明首创研究了点分布模型的插值算法。研究结果表明，可由PDM的平均形状、特征值和特征向量等参数生成任意时刻的中间模型。事实上，插入模型的平均形状及形变就是它相邻前后模型对应项的线性组合，即yt＝aya+byb，Δyti＝aΔyai+bΔybi。其特征值则可以采用式λt≈a2λa+b2λb计算得到。最后，为求出插入模型的协方差矩阵和特征向量，本发明公开了交变矩阵形式及其相应的计算方法。
权利要求
1、一种基于插值算法的柔性体点分布模型的特征值和特征向量的计算方法，其特征在于所述计算方法包括以下步骤
1)、载入两个相邻时刻Ta和Tb的柔性体分布模型Ma和Mb，其计算式为
Ma＝<ya，Φa，λa>，Mb＝<yb，Φb，λb>；
2)、所述柔性体分布模型Ma和Mb，包括计算平均形状ya和yb、特征矩阵Φa和Φb、以及特征值λa和λb，设定柔性体从Ta到Tb呈连续变化，则
2.1)中间模型的特征值的计算
λt≈a2λa+b2λb，(25)
这里λa和λb分别为前后两个相邻PDM模型的特征值向量，即
λa＝(λa1，λa2，λa3，...)，(26)
λb＝(λb1，λb2，λb3，...). (27)
2.2)、中间模型的特征向量的估算
插入中间模型的协方差矩阵近似为
St≈a2Sa+b2Sb (28)
其中，协方差矩阵Sa和Sb由特征矩阵Φa和Φb近似得到，Φa和Φb从插补前已知的点分布模型中得到；设定中间模型建立前先对训练集数据进行PCA分析以减少参数个数，并选取t种形变模式来描述整个对象的形变；由于Φ是协方差矩阵中最前面的t个特征矢量，则有
其中，Λa和Λb是由特征矢量λa和λb组合形成的对角矩阵；用奇异值分解的方法从式得到插入中间模型的特征矢量
由两个已知的相邻PDM模型及相关参数a、b、Λa、Λb、Φa和Φb确定正交矩阵Ut，所述正交矩阵Ut的前t行就是特征矩阵Φt。
2、如权利要求1所述的基于时变插值算法的柔性体点分布模型的特征值和特征向量的计算方法，其特征在于所述步骤2)中还包括
2.3)、中间模型的协方差矩阵由下式推导得出
令
因为
显然有
根据协方差矩阵，矩阵St由Sa、Sb和Sab计算得到，通过如下公式
用奇异值分解法，可令
通常Φab≠Φba，因为矩阵Sab并非对称矩阵，Φab和Φba都包含t个矢量；
在插值计算时，点分布模型保留特征值和特征向量，中间模型的协方差矩阵通过下式计算
3、如权利要求1或2所述的基于插值算法的柔性体点分布模型的特征值和特征向量的计算方法，其特征在于所述步骤2)中还包括
2.3)、特征向量组的计算
通过式(31)来计算特征向量时，将(31)改写为
其中是一个r×N维的矩阵，其中r是主成分数目，N是样本维数，即标定点的数量与标记点维度的乘积；
如果Ui是矩阵G＝LΦa的主要特征向量，对应特征值为λi，则(Φaui)和λi分别为矩阵S＝ΦaL的特征向量和特征值。
由Gui＝λiui，，得
S(Φaui)＝ΦaL(Φaui)＝ΦaGui＝Φa(λiui)＝λi(Φaui)(36)。
全文摘要
一种基于插值算法的柔性体点分布模型的特征值和特征向量计算方法，包括以下步骤1)载入两个相邻时刻Ta和Tb的柔性体分布模型Ma和Mb，其中Ma＝<ya，Φa，λa>，Mb＝<yb，Φb，λb>；2)通过所述柔性体分布模型Ma和Mb，分别计算平均形状ya和yb、特征矩阵Φa和Φb、以及特征值λa和λb，设定柔性体从Ta到Tb呈连续变化；2.1)中间模型的特征值的计算；2.2)中间模型的特征向量的估算。本发明能够形成实时分布的中间过渡模型、模型精度高、实用性良好。
文档编号G06T17/00GK101639948SQ20091010218
公开日2010年2月3日 申请日期2009年8月20日 优先权日2009年8月20日
发明者陈胜勇, 杜雅慧, 秋 管, 毛国红, 王万良 申请人:浙江工业大学