专利名称：用于分割磁共振心脏图像中左心室的变化的方法
技术领域：
一般而言，本发明针对用于分割心脏图像的系统和方法。具体地说，本发明针对一种使用将视觉信息和解剖学约束结合的轮廓传播模型来分割心脏左心室的系统和方法。视觉信息包含基于梯度向量流的边界部分，和根据反映不同组织特性的全局明暗度属性来区分心脏轮廓/区域的区域部分。根据心脏轮廓的相对距离，解剖学约束将连接心脏轮廓的传播。
背景医学图像处理是一个不断增长的应用领域。磁共振成像(MRI)可用于提供心脏随时间变化的三维图像，该图像可用于诊断。该心脏图像可用通用计算机视觉技术进行处理。可是，这些技术不能充分利用其它领域现有技术的优点，如生理学，心脏病学等等。
在医学图像分析中心脏分割是一个大家都尽力研究的课题，因为该成果有很强的诊断功能。在医学结构的分割中，准确和精确是两个重要的需求，因此，开发了各种用于分割的边界驱动方法。这些方法基于边界图像的产生和连续结构的提取，该连续结构准确说明了边界信息。
例如，大家所熟知的“蛇”型模型是一种领先的框架，该框架是重要的边界驱动图像分割技术的基础(参看如，国际计算机视觉杂志1988年第一期321-332页，Kass等人的“SnakesActive Contour Models(蛇型活动轮廓模型)”)。简要地说，蛇型模型涉及能量最小化技术，该技术寻求基于曲线的目标函数的最低势能。该函数是边界图像驱动引力项和说明曲线的期望的内部特性的项之间的折衷。根据最初的蛇型模型，需恢复的结构涉及在二维或三维空间中的一组点，所述空间在满足某些内部特性的约束的同时，局部朝向期望的图像特性的变形。尽管该最初方法可得到显著的效果，但是该方法对于噪音或物理破坏的数据很敏感。
建议在最初的蛇型框架上使用可变形模板和参数化的蛇型模型如B型样条和活动外形来克服其在噪声数据上的较差性能。因此，在分割过程之前，使用一定数量的训练实例来建立外形模型。该模型涉及很少参数，可以表达相当大的一组全局和局部变形。在普通模型可以描述相当大部分的符合条件的分割解的情况下，使用这些方法可以提高分割性能。可是，这些模型对于初始条件相当敏感。
“近视”是常被用来描述蛇型模型对初始条件的依赖性的术语。该蛇型模型基于很局部的信息，因此初始条件必须和最佳解足够接近。建议使用气囊力(ballon force)使该模型脱离于初始条件。中心思想是引入一个不断扩展或收缩初始结构的恒力。很明显，在初始结构或者环绕着需恢复的区域，或者被要分割的区域完全包围的条件下，这一部分对最初的模型有好处。可是，对于一般医学图像分割应用而言，这一条件却不容易被满足。
为了克服这些约束，也考虑并广泛研究局部/全局信息的使用。演进的轮廓用来定义由两个区域组成的图像划分。内部区域是指要被恢复的区域，外部区域是指所述图像的其余部分(背景)。然后，使用全局同质局部属性来区分重要区域和背景。可以使用连续概率密度函数为这些特性建立模型，该函数被根据最新分割图自动更新。
基于边界的医学图像分割技术的演进形成了一组模块，所述模块使得初始结构(点集)向着期望的图像特性变形。基于这些考虑，很明显该分割结果将必然依赖于初始结构确定的参数化(位置、控制点数等等)。另外，用于重新对演进中的结构参数化的技术也将在分割过程中起重要作用。尽管使用各种技术来处理这些问题，但是这种依赖性对于图像分割方法来说是不正常的。
使用水平集表示作为用于演进的界面的一种可选的(针对拉格朗日)技术。这些表示对于计算机视觉中变化框架的实施来说是一个共同的选择。使用连续的更高维的零水平集函数表示该演进的轮廓。这种表示可以是固有的，性质的，免参数和拓扑的。使用用于演进界面的水平集方法导致了用于图像分割的边界驱动方法的扩展。
正如该领域内所熟知的，测地活动轮廓是指作为针对最初蛇型模型的几何可选方案被引入的最优化框架。这一模型的主要优点在于它隐含的参数化，所述隐含的参数化可导致拓扑变化(合并/分离)的自然处理(通过水平集)。
还有一种措施是使用水平集表示实施将边界驱动流和全局区域明暗度信息以及它们的实现加以组合。然后，分割程序变得相当独立，不受初始条件的约束。当可以使用全局区域信息处理任意初始条件时，可以通过水平集表示来处理拓扑变化。在结构非常复杂并且是由多个部分组成的医学图像处理中，上述方法非常重要。可是，当需要处理噪声和不完全数据时，所述方法仍然遭受稳健性的困扰。另外，这些方法不能充分利用如下的现有形状知识的优点，该知识可以从与要被分割的医学结构有关的生理学中获得。
为解决该约束，已经作出了一些努力。例如，现有形状知识已经被引入到测地活动轮廓模型。在另一种方法中，将一个形状影响项和边界及区域驱动的视觉信息组合起来，以便进一步增加基于水平集的方法针对噪声和不完全数据的稳健性。人们可以要求使用蛇型模型和可变形模板满足该目的。
许多医学应用都包含同时提取以约束方式(生理学)被置于图像平面中的多个结构。根据感兴趣结构的相对位置，这些高级(抽象)约束可被转换成低级分割模块。例如，根据大脑的一些物理属性，一种用于皮层分割的方法考虑到两个轮廓的约束(连接)传播。
如上所述，磁共振成像(MRI)提供了心脏随时间变化的三维成像法(imagery)。为帮助诊断疾病，内科医生们对识别心室，心内膜和心外膜很感兴趣。另外，测量心动周期不同阶段的心室血量，心室壁质量，心室壁运动和壁厚特性是一项富有挑战的任务。左心室特别重要，因为它将含氧血从心脏泵到整个身体的远距离组织。
有一些用于心脏分割的方法。例如，Argus是一种西门子通过MRease工作站商品化的心脏磁共振MR分析包，该工作站连接在磁共振“MAGNETOM”系统上。该系统可自动执行3D/4D数据集的分割(不同时间的体积切片)。该分割算法包含三个不同模块并且在如下参考文献中有所阐述，所述文献即在2001年加拿大Vancouver的IEEE计算机视觉国际会议，M.Jolly的“Combining Edge Region and ShapeInformation to Segment the Left Ventricle in Cardiac MRImages(将边缘区域和形状信息合成以分割在心脏磁共振图像中的左心室)”。自动定位模块能够最大辨别地近似确定出新图像中的心肌层位置。该系统通过使心脏正负例子分布的Kullbach距离最大化来获悉心脏的灰度级方面(被建立成马尔可夫链模型)。局部变形过程从一个近似轮廓开始，并使用Dijktra的最短路径算法对它进行变形。在近似轮廓周围的递增尺寸的搜索空间内多次迭代地施加该算法，因此加强了好的边缘而弱化了模糊边缘。然后使用图形切割算法以便最终选择作为心内膜的一部分的最好的边缘片。通过将样条曲线拟合到由Dijsktra的算法勾画出的边缘点，来获得心外膜。最后传播模块负责针对局部变形提供一个近似的起点。
发明简述本发明针对用于分割心脏图像的系统和方法，特别针对心脏的左心室。优选地，使用一种水平集方法来分割磁共振MR心脏图像，所述方法以两个心脏轮廓的连接传播为基础。一种优选的分割方法将视觉信息和解剖学约束结合起来。优选地，该视觉信息是通过(i)基于梯度向量流的边界部分和(ii)如下区域项表示的，该区域项根据反映不同组织属性的全局明暗度属性来区别心脏轮廓/区域。根据心脏轮廓的相对距离，使用解剖学约束来连接心脏轮廓的传播。优选地根据一种快速稳定的数值逼近方案，优选地是附加算子分离，使用水平集方法实现作为结果的运动方程。
本发明的一个方面，用于分割心脏图像的方法包含接收心脏图像，使用梯度向量流边界方法确定心脏图像中一个或多个轮廓的边界信息，和确定区域信息以便根据全局明暗度参数区分心脏图像中的区域。当应用解剖学约束以便根据轮廓的相对距离来连接轮廓传播时，应用一种合成边界和区域信息的传播模型，以便将每个轮廓传播到心脏图像中的目标边界。
本发明的另一方面，该目标边界包含心脏左心室的心内膜和心外膜。
本发明的另一方面，该传播模型的边界部分优选地以使用梯度向量流扩展测地活动轮廓模型为基础。
本发明的另一方面，该传播模型包含边界分割模型，区域模型和连接函数的加权的集成。
本发明的另一方面，解剖学约束包含如下的活动连接函数，该活动连接函数用于保存演进的轮廓的确定拓扑。
在以下对优选实施例的详细描述和附图中，将阐述本发明的这些和那些方面，特征和优点，或将变得很清晰。
附图简述
图1a和1b是心脏图像的例图，其中图1a是描述左心室的一部分的心脏图像例图，图1b说明了左心室的感兴趣的轮廓。
图2是根据本发明的一个实施例，用于分割心脏图像的系统的方框图。
图3a，b，c是说明一种根据本发明的边界分割方法的例图，其中图3a是描述左心室的心脏图像的例图，图3b是说明图像中检测到的边缘(如边缘检测器所检测的)的例图，和图3c是说明图像的梯度向量流表示的例图。
图4是说明基于概率密度的区域分割方法的例图。
图5是根据本发明的一个实施例，用于分割心脏图像的方法的流程图。
图6(a)-6(h)是说明根据本发明，在分割过程期间传播两个轮廓以集中心内膜和心外膜边界的例图。
优选实施例的详细描述一般而言，本发明针对用于磁共振MR成像法的心脏分割。目的在于恢复左心室，特别是如下心肌层，心肌层是在心外膜(即心包膜(心脏周围的囊)的最内层)和心内膜(排列在心脏空洞上的薄内膜)之间的区域。图1中示出了根据本发明的分割应用的前后关系。
图1a是心脏的二维图像例图，具体地说是显示左心室横截面的二维心脏切片。图1中的暗环描述了组成左心室壁的心肌层(肌肉)。根据本发明的一个优选实施例，作为一种三形态图像分割应用(i)心内膜RI，(ii)心肌层RO-RI，(iii)背景I-RO，分割方法基于如图1b所示的两个演进界面，心内膜(RI-内轮廓)或心肌层内壁(以I表示)和心外膜(RO-外轮廓)或心肌层外壁(以O表示)。我们回忆RO是心外膜区域，也包括心内膜区域。
尽管磁共振MR成像法的质量优于其它医学形式(CT，超声波等等)，但是和该考虑的应用有关的视觉信息可能会导致误解，物理损坏，有时会不完全。因此，假如处理时不考虑专用约束，它就可能导致非最佳效果。有以下几个原因。一是因为心内膜中乳突状肌肉的存在。如图1中清楚地显示，尽管该乳突状肌肉是心内膜的一部分，但它们涉及不同的视觉属性。结果，边界(因为很强的边缘)和区域(因为明暗度不一致)驱动的分割模块将不能将该肌肉恢复为心内膜的一部分。
右心室的明暗度特性和心肌层的明暗度特性不同，因此可以很容易地检测共同边界。可是，区别心肌层和非心脏部分(根据视觉信息)不总是可行的。和乳突状肌肉情况相反，一些心脏实体和心外膜之间没有任何边界。因此，边界以及全局明暗度驱动的模块可能不能提供近似的分割图。
根据本发明的一个优选实施例，为了解决上述局限性和物理损坏视觉信息，作以下假设(1)两个心脏轮廓，心内膜和心外膜都是平滑的；(2)心外膜(除了乳突状肌肉外)和心肌层的明暗度特性非常不同，沿着它们的边界存在很强的不连续性(高梯度)；和(3)如图1b所示，内轮廓(I)和外轮廓(O)之间法线方向(测地路径)上的距离相对恒定。
一般而言，根据本发明的系统和方法包含一种三形态分割框架，该框架确定两个轮廓，心内膜轮廓和心外膜轮廓的传播。在示例应用中，感兴趣区域(ROI)包含心内膜和心外膜轮廓之间的区域，即心肌层(如图1b中示出)。
根据本发明的一个优选实施例的分割过程是基于集成视觉信息和解剖学约束的几何流的。优选地，该视觉信息是通过(i)边界驱动的部分和(ii)基于全局区域的分割模块表达的。和传统的分割方法不同，根据本发明的分割过程优选地使用双向流，该双向流不受初始条件影响，并可从任一侧捕获对象边界。区域部分包含分组模块，该模块优选地区别心外膜，心内膜和背景之间的明暗度特性。通过EM算法全局地评估这些特性并自动更新它们。
优选地，使用解剖学约束来处理具有易误解的视觉信息的感兴趣区域的部分。特别地，基于心脏解剖学的先验知识，使用连接模型来约束心外膜和心内膜轮廓的相对位置。优选地，采用一种快速稳定的数值技术，附加算子分离，使用水平集表示来实现作为结果的运动方程。极好的试验结果验证了根据本发明的分割过程的优点和准确性。
需要理解的是这里描述的根据本发明的系统和方法可以以不同形式的硬件，软件，固件，专用计算机，或其组合来实现。优选地，本发明采用作为应用的如下软件来实现，所述软件包含如下程序指令，所述程序指令确实地被包括在一个或多个程序存储设备(如，软磁盘，RAM，CD ROM，ROM等等)上，并可由包含合适体系结构的任何设备或机器加以执行。另外还需了解，因为附图中描述的组成系统的模块和方法步骤优选地通过软件执行，所以根据本发明被编程的方式不同，系统部件之间的实际连接(或过程步骤的流程)会有所不同。根据这里的描述，本领域的普通技术人员将能够想到本发明的这些和类似的实现和配置。
图2是根据本发明的一个实施例，用于分割心脏图像的系统10的方框图。需要了解，尽管这里描述的优选实施例针对分割心脏左心室，但是本发明的原理可以与其它解剖部分有关的图像分割方法一起实现。
系统10处理心脏图像数据11，该数据优选地使用MRI方法捕获。使用MRI捕获的该视觉信息在MR序列中是高质量且很准确的。该系统10包含一种多形态分割框架，该框架基于两个轮廓，心内膜轮廓和心外膜轮廓的传播。由边界分割模块12和区域/明暗度分割模块13处理心脏图像数据11。
在一个优选实施例中，该边界驱动分割模块11实施如下的边界模型，该边界模型是众所周知的测地活动轮廓模型的扩展，该模型用于确定最好地考虑了希望的图像特性的最小长度测地曲线。传统测地活动轮廓的关键局限性在于它高度依赖于初始条件。根据本发明，为克服该局限性，该测地活动轮廓流与二维边界域，即GVF(梯度向量流)相集成。该二维边界域是指与如下最佳流相类似的结构，所述最佳流指向到达心脏边界所遵循的方向。如以下论证的，本发明的边界驱动分割模块12的作为结果的流保证了演进界面的规律性，并且能够从任何一侧都可到达心脏边界。
该区域/明暗度分割模块13产生全局/局部信息，以更使得边界驱动分割过程12的性能更稳健。优选地，使用区域项来提供如下分割图，所述分割图最好使用某些先验技术对图像特性进行分组(所述区域项例如在希腊Corfu的IEEE国际计算机视觉会议926-932页Paragios等人的“Geodesic Active Regions for SupervisedTexture Segmentation(IEEE用于受监督的纹理分割的测地活动区域)”中所描述的，在此引入以供参考)。优选地，根据观察的图像密度函数(柱状图)和期望值最小化(EM)算法来确定这样的特性。这一区域部分的最小化导致了几何流，该几何流可处理噪声/损坏数据，并且不受初始条件影响。
另外，因为由所述应用(如乳突状肌肉)的性质所提供的易误解的视觉信息，该系统10还包含解剖学模块15，该模块实施连接函数，以约束心内膜和心外膜界面的相对位置。这一约束的实施导致了如下双向流，该双向流以法线方向传播心脏轮廓，其中该传播将它们的距离保持在一些给定的限值内。
每个系统部分12，13和15都趋向于通过根据各自基础模型部分产生作用力来传播一些轮廓。换句话说，对于输入的心脏图像的每个相关像素而言，每个模型(边界，区域，解剖学的)将与由该模型确定的优先测量相关联，该优先测量将用于使得所述传播沿着法线方向朝着最佳解驱动。
该系统包含集成模块14，该模块集成数据驱动(视觉信息)和解剖学约束，并导致用于左心室分割的最终几何流。更具体地说，在一个优选实施例中，作为结果的几何流的解是通过在同一平台上集成模型(边界，区域，解剖学的)而获得的，在该平台上所述解包含这些模型的加权合成。
作为结果的解包含两个轮廓，一个对应内轮廓(心内膜)，一个对应外轮廓(心外膜)。图像再现单元16用来根据图像数据11再现估计的轮廓。
边界部分根据本发明的优选实施例，以下详细描述了边界驱动模块12和基础模型。
边界驱动分割模块12优选地实施作为测地活动轮廓模型的扩展的模型。该测地活动轮廓假定一些边界信息的存在。根据该模型，该考虑的应用的最佳分割图由两个界面(一个对应心内膜和一个对应心外膜)而加以确定，所述两个界面最小化以下的目标函数 其中RI(c)，RO(c)分别是心内膜和心外膜边界的平面参数，g()是正单调递减函数(如高斯函数)。
使用变分计算和用于使上述函数最小化的梯度下降法，可以获得以下运动方程(位置u)，以便使心内膜和心外膜向着最佳解传播 其中NI和NO分别是该演进的内界面和外界面的内向法线，K是曲率。
获得的运动方程(2)有简单的解释。它们将(两个)初始轮廓向着由曲率效应约束的心脏边界收缩。这种流可以通过使用水平集表示来加以实现。作为结果的流可以适应拓扑的变化，而不需要初始条件处在最佳解的邻区(它们可以远到满足一些具体约束)。由测地活动轮廓流的性质强加这些约束。该数据驱动项倾向于使得所述界面向着心脏边界收缩(或者假如初始条件保守的话则使之扩展)。因此，初始轮廓必须完全在需恢复的区域内，或者必须包围该区域。这一要求降低了传统测地活动轮廓模型的可用性。
将要理解的是，根据本发明的一个实施例，通过使用GVF(梯度向量流)，可以解决传统测地活动轮廓模型的这一局限性。在美国专利申请系列号No.＿＿(代理人文档号2001P09086US)，名为“GradientVector Flow Fast Geodesic Active Contour(梯度向量流快速测地活动轮廓)”中详细讨论了GVF，该申请同此案一起提出，在此共同指定并被引入以供参考。简要地说，GVF涉及用来从任意一侧捕获心脏边界并可处理凹面区域的二维区域的定义。为引入该流，首先定义基于连续边缘的信息空间。优选地，高斯边缘检测器(零平均值，σE方差)是基于以下梯度图像实现的g(p)=12πσEe-|▿(Gσ*I)(p)|22σE2,f(p)=1-g(p)--(3)]]>其中[Gσ*I]表示高斯内核同输入的卷积(平滑)。
GVF涉及二维边界驱动域[v(p)＝(a(p)，b(p))，p＝(x，y)]，通过最小化以下目标函数可恢复该二维边界驱动域E(v)=∫∫μ(ax2+ay2+bx2+by2)+|▿f|2|v▿f|2dxdy--(4)]]>其中μ是合成参数。大家可以试着解释该函数当缺少边界信息[|f|≈0]时，能量受控于该域的偏导数，导致该域和平滑图的扩散。另一方面，当观测到边界空间上的变化[|f|很大]时，控制能量的项值是第二位的，导致v＝f。在美国Puerto Rico的IEEE计算机视觉和模式识别会议66-71页Xu等人的“Gradient Vector FlowA newExternal Force for Snakes(梯度向量流蛇型模型的新外力”)上可以找到详细解释。
最初的GVF目标函数的一个局限在于并不直接使用边界信息(边界梯度影响流)。结果，因为边界信息的扩散，强边缘和弱边缘可以对该域具有相似的贡献(假定它们的梯度是可比的)。具有优势的是，根据本发明，为解决这种局限，将目标函数(4)作以下修改E(v)=∫∫μ(ax2+ay2+bx2+by2)+f|▿f|2|v-▿f|2dxdy]]>
(5)这种修改可以导致能够克服弱(噪声)边缘的扩散方程，并保证边界信息的合理扩散，其中强边界可以补偿弱边界产生的流。
使用变分计算可以将目标函数(5)最小化，得到以下扩散方程dvdt(p)=m▿2v(p)-f(p)(v(p)-▿f(x,y))|▿f(p)|2--(6)]]>方程(6)优选地用于估算梯度向量流。根据该偏微分方程(PDE)(6)，将首先在具有重要边界信息的图像位置中定义/估算该域，并将该域逐渐传播到远离该边界的区域。当该流收敛时，该域将包含关于最优地到达心脏边界所遵循的方向的信息。然后，人们可以将GVF域规一化[v^(p)=v(p)|v(p)|].]]>作为结果的二维流和光学流类似，可以用于定义一种新的边界驱动外力。由该部分提供的信息主要是有前后关系的，通常可以从任意一侧指向心脏边界。在图3中描述了GVF。具体地说，图3a是输入图像的例图，图3b说明了由边缘检测器检测到的相应的边界信息，图3c说明了相应的梯度向量流。
在上面提到的Xu等人的“Gradient Vector FlowA new ExternalForce for Snakes(梯度向量流蛇型模型的新外力)”中建议了基于蛇型的双向几何流的基础域。和最初的蛇型模型相比，他们的方法表现出很多优点，而与此同时又都有重要的局限性拉格朗日实现。可是，如果给定应用环境，侧这不是一个重要限制。事实上，根据本发明，左心室分割问题涉及凸起的单种子区域(心外膜，心内膜)的恢复，因此不要求对拓扑变化进行处理(这是拉格朗日方法的局限性)。
应当理解，根据本发明的一种分割方法方便地使得能够改变演进的界面的拓扑。优选地，使用多种子可以更好地初步估计心内膜和心外膜的全局/区域明暗度特性。因此，可以取得准确的分割结果。
另外，如以下进一步详细解释的，解剖学约束优选地基于心内膜和心外膜界面之间的距离，使用水平集表示是一种用于估计该距离的优选方法。另外，心脏轮廓的几何属性可以通过这种水平集表示被进行性质上的估计。
记住，该规一化的GVF涉及到达目标边界所遵循的方向。因此，如果给定该演进界面的最新位置，则到达心脏边界的最佳轨迹(从视图的环境来看)将沿着GVF的方向移动。
基于上述讨论，根据本发明的一个实施例的分割方法确定了一种用于心内膜和心外膜的、按照该规一化的GVF的(、基于测地活动轮廓的)前后传播，其中该传播发生在法线方向，并且基于在给定位置[u]的界面法线和GVF向量的内积 该运动方程式(7)沿着指向轮廓边界的方向，移动左心室的内界和外界(即心内膜，心外膜)的轮廓。该运动方程沿着基于边界位置的两个方向中的一个移动轮廓。更具体地说，对该运动方程(7)作以下解释。当GVF指向外且法线向量有相同或相反的方向时，该演进的界面最大地扩展。相反，当GVF指向内且单位法线向量有相同或相反方向时，该心脏轮廓最大地收缩。另外，当GVF与法线相切时，界面位置不被更新。
运动方程(7)不考虑视觉信息。GVF的规一化引起了前后边界驱动的信息。高梯度和强边缘的像素涉及与在这一信息不存在的地方的具有高梯度和强边缘的相类似的实体。在测地活动轮廓流的例子之后，我们可引入以下的图像驱动的边界信息 当边界信息不存在时，设计好的流表现出一种前后关系的传播，这一传播能够从任意一侧到达心脏边界。另一方面，因为假定已经恢复心脏边界是部分令人满意的(强边缘)，所以图像驱动边界信息的存在会降低GVF的共享。
内部约束是蛇型模型和其它边界驱动分割流的关键部分。这些约束可以保证演进界面的规律性和用于实现该相应流的许多方法的稳定性。
根据本发明的基于GVF的边界分割流不考虑这种类型的约束。事实上，人体心脏是一个物质实体，该物质实体包含时间和空间上的平滑结构。如果没有观察到不正常的条件，则这一先验知识就被投入到和由内部约束产生的模块相类似的低水平分割模块中。在测地活动轮廓的例子之后，弯曲部分优选地用来施加内部约束，并将先验的抽象知识转化成低水平分割模块(使用心脏形状的先验模型约束可以是引入这一知识的更合适的方式)。在以下方程中，基于轮廓是平滑的这一假设，施加约束 其中α是一个合成常数。项(1-α)KI(u)强加了一个约束，即轮廓是平滑的，项α(v^(u)·NI(u))]]>表示轮廓向着轮廓边界的移动。
概括来说，通过结合测地活动轮廓和梯度向量流，根据本发明的一个实施例的边界驱动心脏分割方法可以方便地处理任意初始条件，并能处理拓扑变化。
统计区域部分以下是对根据本发明优选实施例的区域/明暗度模块13和基础模型的详细讨论。一般而言，基于区域/明暗度的分割方法优选地依赖于空间定位特征和属性的一致性。根据本发明，左心室分割问题满足这一准则。该区域/明暗度模块13之后的中心思想是使用演进界面来定义如下图像划分，该图像划分相对于某个分组准则来说是最佳的。在医学成像中，观测到的明暗度取决于相应的被映射的组织的属性。
另外，对于左心室的MR序列而言，假定存在三种群体(i)血液(亮)，(ii)肌肉(灰)和(iii)充满空气的肺(深灰)。这些群体的特性随着空间和时间而变化，但是它们的明暗度属性可以被相当好地区别开来。因此，观察到的心外膜区域的分布(柱状图)可以认为是三部分的混合模型(假定是高斯型的)。假定pI表示心内膜密度函数，pO表示心肌层密度函数，pB表示心脏器官的其余部分(背景)的密度函数。那么，我们定义p(I)＝PIpI(I)+POpO(I)+PBpB(I)(10)其中PI，PO和PB是心内膜，心肌层和背景的假设的先验概率。可以使用期望值最大化原则来估算该基础模型的未知参数。图4是根据三个假设，说明视觉信息的统计学模型(柱状图)的例图。给定一个心脏图像，该图像的希望区域被通过计算该区域的柱状图来加以处理，以便获得密度。在一个实施例中，我们假定心内膜是亮的，肌肉是灰的，而背景是黑的(如图4a所示)。在图4b的示例柱状图中，曲线C1是在图像中观察到的总体数据，曲线C4代表用于心内膜的假设，曲线C3代表用于心肌层的假设，曲线C2代表用于背景的假设。
概率密度函数是基于区域的划分模块13的基础。后验分割概率的最大化可以被认为是最优化准则E(∂RI,∂RO)=]]> 其中r()是正单调递减函数，如高斯型的。在方程(11)中定义的模型被用于区别不同的部分。对于本领域的普通技术人员来说，该定义的模型的解释是很显然的。视觉信息必须最大地支持分割图。因此，和一些视觉信息相关联的一个给定的像素，当且仅当它被根据不同假设[pI，pO，pB]的概率法则分配到最有可能的标签时，才有助于目标函数的最低势能。这种分配最小化地给目标函数添加负载。
使用Stokes定理，变分计算和梯度下降方法，可以获得以下运动方程(该方程将能量方程(11)最小化)
(12)该流(公式12)朝向分割图(沿着法线方向)收缩或扩展演进界面，更好地说明了心外膜和心肌层的期望的密度属性。另外，它可以看作是自适应的气囊力，该气囊力根据观察到的图像特性使得演进界面收缩或膨胀。
实现该区域/明暗度模块13的好处是明显的因为该基于区域的流起到了数据驱动的自适应气囊力的作用，因此该建议的方法较少依赖初始条件。可是，由于易误解的视觉信息(由于该应用的性质，心外膜和背景之间基于明暗度的区分容易进行)，可能取得次于最佳结果的结果。
根据本发明，区域驱动流的结果(图2的模块13)和基于边界GVF的分割模块(图2的模块12)(经由图2的模块14)被集成起来，得到以下运动方程 其中β和ε是合成参数。
这些方程包含三项第一项是分组部分，该部分根据要恢复的区域的全局明暗度属性加以估算。该项说明了心肌层的同质性。第二项是内在部分，该部分保持了演进界面的规律性。最后一项是边界吸引部分，该部分从任意一侧向着心肌层边界传播心脏轮廓。
解剖学部分以下是根据本发明的优选实施例，对图2的解剖学模型15的详细讨论。
由于心脏组织的属性，视觉信息可能提供易误解的分割结果。纯粹基于图像/数据驱动技术的分割方法没有充分利用和感兴趣的医学结构有关的先验知识(如心脏生理学)。根据本发明的另一方面，解剖学的先验知识优选地被用作分割过程的一个部分。
在一个实施例中，使用解剖学驱动的低水平分割过程来约束心外膜和心内膜之间的相对位置。事实上，演进界面涉及作为物质实体的一部分的心脏部分，心脏结构。假定它们的位置和它们在时间和空间上的演进是连续的，并且是符合心脏的运动。
根据本发明的一个实施例，解剖学模型假定心外膜和心内膜之间的距离在空间上(切片)和时间上(心脏周期的实例)是变化的。通过确定该距离在时间和空间上的演进，该模块根据它们的最短(测地)路径连接该演进界面。为方便这一约束限制的引入，将忽视心内膜心外膜距离在时间上的演进。
考虑来自心内膜或心外膜轮廓上的一个像素(x，y)。假如d(RI，RO)是(x，y)和该演进界面之间的最小欧几里得距离，那么-((x，y)∈RI)→ d(RI(x，y)，RO)＝d((，y)，RO)，(14)((x，y)∈RO)→d(RIRO(x，y))＝d((x，y)，RI)(15)在美国Santa Barbara的IEEE计算机视觉和模式识别会议708-715页Zeng等人的“Volumetric Layer Segmentation UsingCoupled Surfaces Propagation(使用连接表面传播的体积层分割)”中，建议在皮层分割中使用一个被动连接项，取得了显著的结果。该方法后面的中心思想是引入一个连接函数，该连接函数用于在皮层轮廓的距离超过可接受的限值时，就停止它们的传播。可是，这是一种严格的限制，并不演进心脏轮廓，而且代表一种消极的行为。
根据本发明的一个优选实施例，使用一种活动连接模块，该模块满足用于心内膜(内轮廓)的以下条件(i)假如用于给定像素的从心外膜的距离在可接受的限制[m≤d(RI(u)，RO)≤M]内，那么就满足所述约束，不会采取行动；(ii)假如这一距离低于可接受的最小值[d(RI(u)，RO)＜m]，那么心内膜离心外膜非常近，并且心内膜必须收缩以保持距离约束；和(iii) 假如这一距离超出可接受的最大值[d(RI(u)，RO)＞M]，那么心内膜就远离心外膜，并且心内膜必须扩展以满足所述距离约束。
基于这一公式化，根据本发明的一个实施例的用于心内膜界面的连接函数优选地定义如下 对心内膜(外轮廓)使用相同的推理，可以获得类似的条件(i)假如[m≤d(RI，RO(u))≤M]，那么满足所述约束，并且不采取任何行动；(ii)假如[d(RI，RO(u))＜m]，那么心外膜离心内膜非常近，并且心外膜必须扩展以保持距离约束；和(iii)假如[d(RI，RO(u)＞M]，那么心外膜远离心内膜，并且心外膜必须收缩以满足所述约束。
基于这些条件，根据本发明的另一个实施例的连接函数针对心外膜界面(外部轮廓)被定义成如下 优选地实现(在图2的解剖学模块15中)连接函数(16)和(17)，以提供一个活动解剖学模型，该活动解剖学模型将心脏轮廓之间的距离保持在一些可接受的限值内。该模型提供了一种活动的传播力，该作用力沿着法线方向朝着满足一系列符合本发明的距离要求的位置局部演进心脏轮廓，如下 该连接函数(18)提供许多优点。例如，该连接函数可以对所希望的分割算法提供强的贡献或柔性贡献。另外，就恢复满足所述约束的演进界面的拓扑而言，该连接函数扮演了一个主动的角色。这和传统方法(Zeng等人)相反，传统方法定义了一种严格的约束，该约束减少了可接受的拓扑/解的集合。
完全模型应当理解，本发明将数据驱动(视觉信息)流和解剖学约束结合起来，以提供用于左心室分割的几何流。根据本发明的分割模型不涉及最优化准则，但是这样的模型充分利用了该应用框架的优点，并且使视觉信息的输出最大化。根据本发明的一个实施例，用于分割的一个完全模型优选地被定义如下 其中(是合成参数。这些几何流包含四种不同的作用力，这些作用力都沿着法线的方向作用(i)说明演进界面的内部属性的基于曲率的项；(ii)向着其真实边界传播心脏轮廓的边界驱动双向作用力；(iii)说明心脏区域同质性的明暗度驱动区域作用力；和(iv)和心内膜和心外膜的位置有关的解剖学约束。
实现问题拉格朗日方法是实现几何流的统一中公共方法。这些技术采用一小组(控制)点来表达演进界面(可使用关于两个连续点之间的距离的统一取样规则)。然后根据运动方程更新该控制点的位置。
在过去的十年期间，关于用于实现这些技术的数值方法取得了巨大的进步。因此，它们是许多有良好性能的基于蛇型模型的医学图像分割技术的基础。可是，上述这些技术在各种方面还有限制。
例如，改变所述演进界面的拓扑是拉格朗日技术的一个主要的局限性。已经提出了一些方法来克服这一局限性，为涉及曲线和表面的几何流提供良好的性能，但是上述方法却不能推广到更高维。另一个局限性是估计几何属性如该演进表面的法线和曲率。通常需重新对演进结构参数化。根据取样规则进行重新参数化。因此该界面特性的估计高度地依赖于这一规则。
水平集表示是在拉格朗日方法的候选项，(参看，例如计算物理学杂志79期，12-29页，Osher等人的“Fronts Propagating WithCurvature-Dependent SpeedAlgorithms based on theHamilton-Jacobi Formulation“以由曲率决定的速度进行的前向传播基于Hamilton-Jacobi公式的算法”，在此引入以供参考)。该方法后面的中心思想是将该演进界面RI(u；t)，RO(u；t)如下表示为更高维函数[φI，φO]的零水平集 [RI(u)，RO(u)]的这种表示是隐含的，无参数的和固有的。另外，它们是无拓扑的，因为零水平集的不同拓扑不意味着[φI，φO]的不同拓扑。
我们考虑一个几何流的最简单的形式∂∂t∂Rx(u,t)=Fx(K∂Rx(u,t))N(u,t)--(21)]]>这一几何流沿着根据取决于界面几何属性(即曲率)的某个标量函数F的法线方向演进初始轮廓。(嵌入的)函数φ的基本导数得出以下水平集流ddtφx(u,t)=Fx(K∂Rx(u,t))|▿φx(u,t)|--(22)]]>其中N=▿φ|▿φ|.]]>这一水平集流的解与嵌入的函数φ的选择无关。离演进界面的带符号的(signed) 欧几里得距离是φ的公共选择。
优选地，该水平集表示用于实现心脏分割流。唯一需要处理的问题涉及用于连续估算心外膜和心内膜的演进界面之间距离的机制。通过将带符号的欧几里得距离当作嵌入函数可以解决这一问题。这一选择反映了该演进界面之间的最短距离。根据本发明的一个实施例，将心脏分割流投影到该水平集空间优选地如下进行 可是连接项的定义涉及在该定义流的数值实现期间可能引起不稳定性的问题的间断函数。为解决该问题，我们考虑一个和一个初始函数有同样属性的平滑函数 需要注意，最终模型包含几个参数(α，β，ε，γ，δ，m，M)。当然，不同信息驱动的分割模块的使用可能提高作为结果的算法的性能。可是，将它们恰当地集成可能是一个问题，对我们的方法来说也是一样。
基于试验，已经确定该区域/明暗度信息是最可靠的视觉源(β＝0.5)。边界部分也有稳定的行为(这说明了明暗度信息的不连续性)(ε＝0.2)。另外，必须加强该演进界面的规律性(α＝0.25)。对于严格的限值来说，解剖学模块是柔性的(γ＝0.3)。根据考虑的切片的空间位置定义距离限值限制。中间心脏切片对应左心室的最大体积，因此该限值尽可能地宽(m＝4，M＝8，δ＝2)。
现在，我们处在开发自动方法来确定解剖学模块参数的过程中。在左心室的大多数部位种可以合适地提取心肌层。因此，和心肌层有类似的视觉属性的乳突状肌肉和结构相对于演进的心脏轮廓之间距离而言可以可外露层。基于稳健性统计学的方法可用于恢复这一距离的平均值和可变性，导致了解剖学约束在空间上和时间上的自动测量。
数值实现用于实现该建议框架的快速的稳健的数值方法优选地基于所述的熟知的附加算子分离(AOS)方案，如，例如I EEE图像处理学报1998年第7期398-410页Weickert等人的“″Efficient and ReliableScheme For Non-Linear Diffusion and Filtering(非线性扩散和过滤的有效可靠方案)”中描述的，该文在此引入以供参考。在计算机视觉比例空间原理国际会议(1999)34-35页Goldenber g等人的“Fast Geodesic Active Contours(快速测地活动轮廓)”中，已经成功地使用水平集方法在平面曲线演进中应用了该技术。
在计算机视觉中使用PDE(偏微分方程)的一个局限性是低效率。经典的数值逼近是不稳定，这导致了耗时的方法。这是由于需要用于保证稳定演进和收敛到该PDE的短时间步骤而导致的。在以上结合的Weichert等人的文章中介绍了一种克服该局限性的方法，该方法被有效地用于为各种PDE提供一种稳定的数值方法。
为了说明起见，通过假定以下形式的扩散方程来考虑一维情况tu＝div(g(|u|)u)(25)然后，该扩散方程可以作如下离散tu＝x(g(|u|)xu)(26)得到以下迭代方案
其中I是单位矩阵，τ是时间步骤。尽管该系统使用前次迭代的值显式地更新u值，但是当时间步骤受上界约束时并不稳定。我们可以考虑使用半隐方案 该方案具有稳定的行为，但是在计算上很开销很大。该AOS技术涉及该半隐方案的以下修改um+1＝[I-τ4(um)]-1um(29)附加算子分离的具有一些优势属性，因为它是(i)稳定的；(ii)满足离散的非线性扩散的所有准则；(iii)具有低的复杂度(和像素数目成线性)；和(iv)可以很容易地向更高维扩展。因此，优选地将AOS方案和这里描述的用于左心室的分割的水平集运动方程一起加以应用。
为进一步降低水平集传播所需的计算费用，该AOS方案可以有效地结合窄带方法(参看，如计算物理学杂志118期(1995)269-277页Adalsteinsson等人的“A Fast Level Set Method For PropagatingInterfaces(用于传播界面的快速水平集方法)，该文在此引入以供参考)。该方法的性质在于只在一个有限的带宽内执行该水平集的传播。根据该传播的轮廓的最新位置定义该带宽。因此显著减少了感兴趣区域，从而显著地降低了计算的复杂度。该方法需要频繁地重新初始化该水平集表示(距离函数)，可以使用快速行进算法有效地进行(参看，如剑桥大学出版社1996年J.Sethian的“Level Set Methods(水平集方法)，在此引入以供参考)。参考IEEE解剖学控制学报40期(1995)1528-1538页J.Tsitsiklis的“Efficient AlgorithmsFor Globally Optimal Trajectories(用于全局最佳轨迹的有效算法)”，提出了在自动控制区域的类似算法。
结论，试验结果概括地说，本发明提供了用于分割MR心脏图像中的左心室的几何流。我们的方法是基于在规律性，边界，区域和解剖驱动作用力的影响下的两条曲线(心脏轮廓)的传播。该框架使用了无参数曲线，并使用水平集表示来实现。几个MR心脏序列被用于验证我们的已经导致成功的试验结果的方法。
图5是根据本发明的一个实施例，用于分割心脏图像的方法的流程图。图6(a)到6(h)是说明传播两个轮廓，心内膜(内轮廓)和心外膜(外轮廓)以便将该轮廓收敛到左心室的心肌层的期望的内边界和外边界的例图。
参考图5，根据本发明的一个实施例的分割过程开始于接收图像数据(步骤20)。该图像数据包含在MRI采集期间一给定时间所捕获的几个二维切片。例如，图6(a)-(h)中的每一个都描述了在一给定时间，在三维体积的一给定深度上捕获的同一二维图像。为这些初始轮廓选择初始的任意条件，并再现这些初始轮廓(步骤21)。例如，图6(a)说明了初始(任意)的内轮廓和外轮廓，它们被(经由分割过程)向着图6(h)示出的最终解传播。优选地，初始的内轮廓在心内膜之内。
为捕获用于分割的视觉信息，例如使用方程(6)确定该图像的梯度向量流(步骤22)。然后，根据GVF确定边界项(步骤23)。如上所述，对于每个像素，GVF的收敛产生了一个(从该像素)指向朝着心脏边界的方向的向量(参看图3c)。同时(假定并行处理)，使用和该图像有关的柱状图数据确定概率密度(步骤24)。如上所述，每个概率密度都和感兴趣的区域和假设有关，这允许该图像的像素基于预定的准则分组。通过像素的分组确定区域项，以便描绘图像中的结构和区域。例如，正如以上参考图4所提到的，心肌层，心内膜和背景区域被区别开来，并从那里确定轮廓。
接下来，使用集成的传播模型处理边界和区域数据(视觉信息)，该传播模型应用连接函数(解剖学约束)来迭代地确定实际轮廓边界(步骤27)。例如，方程(19)的传播模型优选地用于确定左心室心肌层的内外壁(心内膜和心外膜)的实际边界。换句话说，从初始的任意轮廓开始，该轮廓以这样的方式传播，以至于在满足应用于完全传播模型的解剖学约束的同时，来考虑边界和区域信息。随着迭代过程的进行，在图像上可以再现轮廓的被更新位置，如图6a-h中所示。重复该过程(步骤26-27)，直到存在解的收敛(步骤28中的肯定结果)为止。对需处理的每个图像切片重复该分割过程(步骤21-28)。对于一个三维实现来说，通过逐片使用该二维方法(图5)，并且将该结果放在一起，可以恢复三维分割图。
概括地说，根据本发明的几何流模型和分割过程的各种新颖的特征包含(i)使用和初始条件无关的新颖的柔性外部边界项；(ii)将边界驱动的和基于区域的信息分割模块集成起来；(iii)使用应用的前后关系来定义如下传播约束，该传播限制可处理物理损坏和不完全数据；(iv)使用关于PDE的数值技术的最新发展来稳健稳健地且高收敛率地实现获得的运动方程。
该分割算法的扩展包括使用连接的表面传播的三维分割。我们期望在水平集框架中可以容易地进行这样的扩展。实现三维流的这些表示的使用在视觉和图像处理著作中已经是标准技术(现有技术)。
在本发明的其他实施例中，在水平集能量框架中使用该传播模型可以实现形状优先(参看，如哥本哈根2002年计算机视觉欧洲会议第二卷78-92页Paragios等人的“Shape Priors for Level SetRepresentations(水平集表示的形状优先)”，该文在此引入以供参考)。有几个问题和该目的有关。第一是形状表示的选择，第二是训练取样的记录，第三是优先项本身。
另外，在医学图像分析中图像记录是一个关键部分。像皮层，心脏心室，肾等等这样的结构是很感兴趣的。相对于一些现有模型或它们先前的状态(过去的成像法)记录这些结构是有很强诊断力的令人感兴趣的应用。为了这个目的，我们将乐于把记录和分割问题结合起来。分割医学结构，并同时对照先前状态将它们记录下来是一个需要考虑的具有挑战性的方向。
尽管在此参考附图描述了用以说明的实施例，应当了解本发明不限于这些精确的实施例，本领域的技术人员在不背离本发明范围的情况下可以进行各种改变和修改。
权利要求
1.一种用于分割心脏图像的方法，包含如下步骤接收心脏图像；使用梯度向量流边界方法，确定心脏图像中的一个或多个轮廓的边界信息；根据全局明暗度参数确定区域信息，以区别心脏图像中的区域；和在心脏图像中应用组合边界和区域信息的传播模型，以便将每个轮廓传播至目标边界，同时应用解剖学约束以便根据所述轮廓的相对距离连接所述轮廓的传播。
2.如权利要求1所述的方法，进一步包含为所述轮廓选择任意初始条件的步骤。
3.如权利要求1所述的方法，其中所述目标边界包含心脏左心室的心内膜和心外膜。
4.如权利要求1所述的方法，其中所述传播模型的边界部分基于使用梯度向量流的测地活动轮廓模型的扩展。
5.如权利要求1所述的方法，其中所述传播是使用可处理拓扑变化的水平集表示来实现的。
6.如权利要求5所述的方法，其中水平集传播是使用附加算子分离方法来实现的。
7.如权利要求1所述的方法，其中所述传播模型包含对边界分割模型，区域模型和连接函数的加权集成。
8.如权利要求1所述的方法，其中应用传播模型的步骤包含将约束应用到所述传播模型的边界部分，所述模型基于所述目标边界的物理性质。
9.如权利要求8所述的方法，其中所述边界部分的约束包含平滑度约束。
10.如权利要求1所述的方法，其中确定区域信息的步骤包含施加数据驱动的自适应气囊作用力，以收缩或膨胀所述演进轮廓，所述演进轮廓基于观察到的图像特性。
11.如权利要求1所述的方法，其中所述解剖学约束包含一个活动连接函数，所述活动连接函数保持所述演进轮廓的确定拓扑。
12.如权利要求1所述的方法，进一步包含如下步骤在分割过程期间在心脏图像上再现轮廓传播。
13.一种可由机器读取的程序存储设备，该程序存储设备具体包含可由所述机器执行的指令的程序，用以执行分割心脏图像的方法步骤，所述方法步骤包含接收心脏图像；使用梯度向量流边界方法，确定心脏图像中一个或多个轮廓的边界信息；根据全局明暗度参数确定区域信息，以区别心脏图像中的区域；和在心脏图像中应用组合边界和区域信息的传播模型，以便将每个轮廓传播至目标边界，同时应用解剖学约束以便根据所述轮廓的相对距离连接所述轮廓的传播。
14.如权利要求13的程序存储设备，进一步包含用于执行为所述轮廓选择任意初始条件的步骤的指令。
15.如权利要求13所述的程序存储设备，其中所述目标边界包含心脏左心室的心内膜和心外膜。
16.如权利要求13所述的程序存储设备，其中所述传播模型的边界部分基于使用梯度向量流的测地活动轮廓模型的扩展。
17.如权利要求13所述的程序存储设备，其中所述传播是使用可处理拓扑变化的水平集表示来实现的。
18.如权利要求17所述的程序存储设备，其中水平集传播是使用附加算子分离方法来实现的。
19.如权利要求13所述的程序存储设备，其中所述传播模型包含对边界分割模型，区域模型和连接函数的加权集成。
20.如权利要求13所述的程序存储设备，其中应用传播模型的指令包含将约束应用到所述传播模型的边界部分的指令，所述模型基于所述目标边界的物理性质。
21.如权利要求20所述的程序存储设备，其中所述边界部分的约束包含平滑度约束。
22.如权利要求13所述的程序存储设备，其中确定区域信息的指令包含施加数据驱动的自适应气囊作用力以收缩或膨胀所述演进轮廓的指令，所述演进轮廓基于观察到的图像特性。
23.如权利要求13所述的程序存储设备，其中所述解剖学约束包含一个活动连接函数，所述函数保持所述演进轮廓的确定拓扑。
24.如权利要求13所述的程序存储设备，进一步包含在分割过程期间在心脏图像上再现轮廓传播的指令。
全文摘要
一种系统和方法，用于分割心脏图像，具体而言，使用集成视觉信息和解剖学约束的轮廓传播模型来分割心脏的左心室。所述视觉信息包含基于梯度向量流的边界部分和根据反映不同组织属性的全局明暗度属性来区别心脏轮廓/区域的区域部分。所述解剖学约束根据心脏轮廓的相对距离来连接心脏轮廓的传播。所述传播模型包含对边界分割模型，区域模型和连接函数的加权集成。
文档编号A61B5/00GK1524247SQ02810121
公开日2004年8月25日 申请日期2002年5月17日 优先权日2001年5月17日
发明者N·帕拉伊奥斯, V·拉梅斯, N 帕拉伊奥斯, 匪 申请人:西门子共同研究公司