基于水平集函数的图像处理转让专利

申请号 : CN201280008480.5

文献号 : CN103403761B

文献日 :

基本信息:

PDF:

法律信息:

相似专利:

发明人 : 汉斯·弗里梅尔王纯亮厄尔扬·斯梅德比

申请人 : 诺华美亚有限责任公司

摘要 :

一种将表示图像的输入数字图像数据处理为输出数字图像数据的方法,包括以下步骤:提供(100)输入图像数据;提供并初始化(200)水平集函数以表示输出数字图像数据的近似值;提供(300)速度函数;确定(400)与水平集函数的图像点相关的演化方向趋势图,所述趋势图包括与所述图像点关联的一个或多个趋势方向,使得每个所述图像点都具有关联的趋势方向;利用速度函数和演化方向趋势图来更新(500,520)水平集函数,使得在速度函数的速度与关联的趋势方向不符的图像点上排除速度函数对水平集函数的更新的贡献;以及基于水平集函数提供(700)输出数字图像数据。选择性地排除速度函数的在图像点上的贡献使得水平集的演化具有连贯性,避免了传统方法中由于局部扭动引起的速度减慢。

权利要求 :

1.一种将表示图像的输入数字图像数据处理为输出数字图像数据的方法,包括以下步骤:提供(100)输入图像数据;

提供并初始化(200)水平集函数以表示输出数字图像数据的近似值;

提供(300)速度函数;

确定(400)与水平集函数的图像点相关的演化方向趋势图,所述趋势图包括与所述图像点关联的一个或多个趋势方向,使得每个所述图像点都具有关联的趋势方向;

利用速度函数和演化方向趋势图来更新(500,520)水平集函数,使得在速度函数的速度与关联的趋势方向不符的图像点上排除速度函数对水平集函数的更新的贡献;以及基于水平集函数提供(700)输出数字图像数据。

2.根据权利要求1所述的方法,其特征在于更新水平集函数的步骤被重复执行至少一次。

3.根据权利要求1所述的方法,其特征在于更新水平集函数的步骤被重复执行,直到演化条件满足。

4.根据权利要求3所述的方法,其特征在于演化条件基于在更新水平集函数的步骤中有多少个图像点被排除和/或基于更新水平集函数的步骤被执行的迭代次数。

5.根据权利要求3-4中任一项所述的方法,其特征在于满足演化条件包括预定数量或用户输入确定的数量的图像点在更新水平集函数的步骤中被排除,和/或更新水平集函数的步骤已被执行预定次数或用户输入确定的次数的迭代。

6.根据权利要求1-4中任一项所述的方法,其特征在于更新水平集函数的步骤还包括将速度函数的速度与关联的趋势方向不符的图像点从更新水平集函数的步骤中排除。

7.根据权利要求1-4中任一项所述的方法,其特征在于更新水平集函数的步骤还包括:如果将被排除的图像点位于被更新了的图像点的局部邻域内,优选速度与关联的趋势方向一致的图像点,则仍利用速度函数对这样的被排除的图像点的水平集函数进行更新。

8.根据权利要求1-4中任一项所述的方法,其特征在于演化方向趋势图是至少通过用户输入初始确定的或者是至少初始预定的。

9.根据权利要求1-4中任一项所述的方法,其特征在于演化方向趋势图的一个或者多个趋势方向是通过基于速度函数的数值的计算来确定的。

10.根据权利要求1-4中任一项所述的方法,其特征在于演化方向趋势图的一个或者多个趋势方向是通过利用由速度函数确定的速度的符号来计算的。

11.根据权利要求1-4中任一项所述的方法,其特征在于演化方向趋势图的一个或者多个趋势方向是基于关联的图像点中的和/或其局部邻域中的速度函数的数值来计算的。

12.根据权利要求1-4中任一项所述的方法,其特征在于演化方向趋势图包括与所有的或一组所述图像点关联的一个趋势方向。

13.根据权利要求1所述的方法,其特征在于所述方法还包括以下步骤:重复(600,620)确定演化方向趋势图的步骤和更新水平集函数的步骤至少一次。

14.根据权利要求13所述的方法,其特征在于在将确定演化方向趋势图的步骤和更新水平集函数的步骤进行重复的步骤中的重复之后的确定演化方向趋势图包括改变一个或多个趋势方向中的至少一个。

15.根据权利要求14所述的方法,其特征在于所述改变涉及转变为反方向。

16.根据权利要求14所述的方法,其特征在于,通过转变表示所述一个或多个趋势方向的值的符号,所述改变涉及转变为反方向。

17.根据权利要求13-16中任一项所述的方法,其特征在于将确定演化方向趋势图的步骤和更新水平集函数的步骤进行重复的步骤中的重复涉及将满足图像点收敛条件的图像点排除在更新之外。

18.根据权利要求17所述的方法,其特征在于,当在前一个或前几个更新期间水平集函数已改变小于优选预定或者用户输入确定的阈值时,图像点收敛条件被满足。

19.根据权利要求13-16中任一项所述的方法,其特征在于将确定演化方向趋势图的步骤和更新水平集函数的步骤进行重复的步骤还包括对于水平集函数的所有图像点或满足图像点收敛条件的图像点将水平集函数设定到基于最近的水平集函数值的值。

20.根据权利要求13-16中任一项所述的方法,其特征在于将确定演化方向趋势图的步骤和更新水平集函数的步骤进行重复的步骤还包括对于水平集函数的所有图像点或满足图像点收敛条件的图像点将水平集函数设定到当前和先前水平集函数值的平均值。

21.根据权利要求13-16中任一项所述的方法,其特征在于将确定演化方向趋势图的步骤和更新水平集函数的步骤进行重复的步骤被重复至少一次,或者直至满足水平集函数收敛条件。

22.根据权利要求21所述的方法,其特征在于满足水平集函数收敛条件基于水平集函数有多少个图像点满足图像点收敛条件,和/或在将确定演化方向趋势图的步骤和更新水平集函数的步骤进行重复的步骤中执行的迭代的次数。

23.根据权利要求21所述的方法,其特征在于满足水平集函数收敛条件基于全部图像点满足图像点收敛条件,和/或在将确定演化方向趋势图的步骤和更新水平集函数的步骤进行重复的步骤中执行的迭代的次数。

24.根据权利要求21所述的方法,其特征在于满足水平集函数收敛条件包括水平集函数的预定数量或用户输入确定的数量的图像点满足图像点收敛条件,和/或在将确定演化方向趋势图的步骤和更新水平集函数的步骤进行重复的步骤中已执行了预定次数或用户输入确定的次数的迭代。

25.根据权利要求1-4中任一项所述的方法,其特征在于将表示图像的输入数字图像数据处理为输出数字图像数据的方法是图像分割方法,其中,输出数字图像数据表示由输入数字图像数据表示的图像中的物体的边界。

26.根据权利要求25所述的方法,其特征在于水平集函数仅在图像点的子集中被更新。

27.根据权利要求25所述的方法,其特征在于水平集函数仅在包括零水平的窄带中被更新。

28.根据权利要求26所述的方法,其特征在于更新水平集函数的步骤还包括检测加入到子集中的新图像点以及关于新点更新演化方向趋势图并从水平集函数值被更新了的相邻图像点继承要与新图像点关联的一个或多个趋势方向并由此使得新图像点被加入到子集。

29.根据权利要求26-28中任一项所述的方法,其特征在于,当新图像点已被通过更新水平集函数的步骤中的更新加入到子集中时,演化条件被视为被满足。

30.根据权利要求26-28中任一项所述的方法,还包括提供(350)先验形状的步骤;并且其中,更新水平集函数的步骤还包括基于先验形状计算(450)变形的形状,并在更新水平集函数中利用变形的形状。

31.根据权利要求30所述的方法,其特征在于更新水平集函数的步骤还包括利用速度函数的速度与关联的趋势方向不符的图像点来计算变形的形状。

32.根据权利要求30所述的方法,其特征在于在更新水平集函数的步骤中,开始计算变形的形状的决策基于:与先前的迭代相比,图像点的速度相对于它们的关联的趋势方向是否已发生了任何变化。

33.根据权利要求1-4中任一项所述的方法,其特征在于将表示图像的输入数字图像数据处理为输出数字图像数据的方法是图像降噪方法,其中,输出数字图像数据表示由输入数字图像数据表示的图像的经过平滑和降噪后的图像。

34.一种计算机(1101),被配置为执行权利要求1-33中任一项的各个步骤。

35.一种医学成像系统(1102),其中集成了如权利要求32中所述的计算机(1101)。

36.根据权利要求35所述的医学成像系统(1102),其中所述医学成像系统(1102)是核磁共振成像系统或计算机断层摄影成像系统。

37.一种用于将表示图像的输入数字图像数据处理为输出数字图像数据的设备,其中所述设备被配置为:提供输入图像数据;

提供并初始化水平集函数以表示输出数字图像数据的近似值;

提供速度函数;

确定与水平集函数的图像点相关的演化方向趋势图,所述趋势图包括与所述图像点关联的一个或多个趋势方向,使得每个所述图像点都具有关联的趋势方向;

利用速度函数和演化方向趋势图来更新水平集函数,使得在速度函数的速度与关联的趋势方向不符的图像点上排除速度函数对水平集函数的更新的贡献;以及基于水平集函数提供输出数字图像数据。

说明书 :

基于水平集函数的图像处理

技术领域

[0001] 本发明总体涉及基于水平集函数的图像处理。

背景技术

[0002] 图像分割常用于医学图像分析过程中,如用于手术计划和血管狭窄的测量。虽然目前已经有许多图像分割算法被提出并发表于科学文献中,但是在图像分辨率有限或者图像噪声较大的情况下,大多数算法并不能给出满意的分割结果。在某些情况下,这些包含错误的分割结果可以通过人机交互的方式来改善,但是这样往往使操作流程变得复杂,并延长了图像分析的时间。上个世纪80年代提出的水平集算法一直以来被认为是一种可靠的图像分割方法,因为它可以自动地处理被分割物体的复杂拓扑变化,并且在高噪声的图像背景下保持分割边界的平滑从而实现亚像素的分割精度。关于这一算法的详细介绍参见Osher和Fedkiw的著作《Level set methods and dynamic implicit surfaces》(Springer,2003)。很多研究表明水平集算法可以有效提高图像分割精度并减少人机交互操作,相关研究可见参考文献:“User-guided 3D active contour segmentation of anatomical structures:significantly improved efficiency and reliability,Neuroimage31(3),1116–1128(2006)”和“Level set based cerebral vasculature segmentation and diameter quantification in CT angiography,Medical Image Analysis10(2),200–214(2006)”。即便如此多的文献证实水平集算法可以在实验室里提供准确稳定的结果,但是在临床工作中这一算法的应用并不很广,这主要归咎于这一算法的运行时间非常长。目前已经有几种加速方法被提出来以提高水平集算法的运算速度,如Sethian等提出的窄带法,这种方法将水平集函数的运算限制在一个几个像素宽的窄带上而不是整个图像,详见“AFast Level Set Method for Propagating Interfaces,Journal Of Computional Physics118,269--277(1994)”。Whitaker进一步提出了稀疏场的算法,将窄带缩小到只有一个像素的宽度,而在每一次迭代中通过距离变换来重建水平集函数,参见“A level-set approach to 3d reconstruction from range data,International Journal Of Computer Vision29,203--231(1998)”。但是在很多应用中,这些方法的运行速度还是不让人满意,图像分割还是需要很长时间。因此,进一步提高基于水平集方法图像分割的运行速度就显得尤为重要。并且这一需求不仅限于图像分割,快速的水平集算法对于其他的各种图像分析任务,例如图像降噪处理等,也具有重要意义。

发明内容

[0003] 基于上述背景,本发明的目的是解决目前通行方法存在的问题,或者至少提供通行方法之外的另外一种解决方案。更具体的目的是阐述一种比传统方法更为快速的基于水平集的图像处理方法。
[0004] 为实现上述目的,本发明经过以下步骤将表示图像的输入数字图像数据处理为输出数字图像数据:提供输入图像数据;提供输入并初始化水平集函数以表示输出图像数据的近似值;提供速度函数;确定与水平集函数的图像点相关的演化方向趋势图,所述趋势图包括与所述图像点关联的一个或多个演化趋势方向,使得每个所述图像点都具有关联的趋势方向;利用速度函数和演化方向趋势图来更新水平集函数,使得在速度函数给出的速度与关联的演化趋势方向不符的图像点上排除速度函数对水平集函数的更新的贡献;基于输出数字图像数据。选择性的排除速度函数的在图像点上的贡献使得水平集的演化具有连贯性,避免了传统方法中由于局部扭动引起的速度减慢。
[0005] 上述演化趋势方向意指预测的速度函数的行进方向
[0006] 上述水平集函数的图像点意指有水平集函数定义的采样点,可以散布于整个图像,也可能局限于图像一部分,例如位于一个窄带上。
[0007] 演化方向趋势图中,可以是所有的图像点都被赋予同一个趋势方向,也可以是一组图像点被赋予同一个趋势方向,还可以是不同的图像点被赋予不同的趋势方向。本发明的图像处理步骤并不局限于以上描述的顺序。
[0008] 上述更新水平集函数的步骤可以被执行一次或者重复多次直到演化条件满足。
[0009] 演化条件往往根据应用的不同而不同,其设计思想是当水平集函数就某个应用的具体要求而言已经演化到满意的程度时,这种重复的更新应该停止。在本发明中,从计算水平集函数演化方向趋势图到重复更新水平集函数至满足演化条件被称为一个周期。在同一个周期内,演化方向趋势图是保持不变的。在演化条件满足后,也就是说在当前演化方向趋势图影响下的演化条件满足后,这一周期结束。
[0010] 上述演化条件可以是基于在更新水平集函数的步骤中有多少个图像点被排除和/或基于更新水平集函数的步骤被执行的迭代次数。满足演化条件包括预定数量或用户输入确定的数量的图像点在更新水平集函数的步骤中被排除,和/或更新水平集函数的步骤已被执行预定次数或用户输入确定的次数的迭代。上述重复更新水平集函数的步骤中还可以包括将速度函数的速度与关联的趋势方向不符的图像点从更新水平集函数的步骤中排除。这种方式避免了无用的重复计算。
[0011] 在某些基于本专利的应用中,上述重复更新水平集函数的步骤中还可以包括如果将被排除的图像点位于被更新了的图像点的局部邻域内,则仍利用速度函数对这样的被排除的图像点的水平集函数进行更新。在此临近区域指的是到所述图像点距离小于一个预设的或者用户给定值的区域。
[0012] 所述演化方向趋势图的初始值可以通过用户设定的方式或者程序推测的方式。
[0013] 所述演化方向趋势图的一个或者多个趋势方向可以通过计算速度函数的数值或者仅速度的正负符号来推测。所述演化方向趋势图的一个或者多个趋势方向可以通过计算在相应的临近区域内所有图像点的速度函数的数值或者仅速度的正负符号来推测。这里的临近区域与前述的临近区域的概念相同。但是这俩临近区域的作用不同因此其定义的区域大小可以不同。在此临近区域可以小到只涉及一个图像点,大到包含所有的图像点。在后一种情况下演化方向趋势图将只有一个单一的趋势方向。
[0014] 所述演化方向趋势图可以只包含一个趋势方向,也就是说所有图像点的预期趋势方向都一样,这样需要存储和处理的成本都较小,从而进一步提高算法的运算速度。
[0015] 所述方法中还可以包含以下额外的步骤:重复确定演化方向趋势图的步骤和更新水平集函数的步骤至少一次。这些重复须要在输出结果图像数据之前进行。在这种周期性重复过程中,一个新的周期开始于计算水平集函数演化方向趋势图的步骤。通过这种周期变换,水平集可以在前一个周期向一个趋势方向演化,而在下一个周期向另外一个趋势方向演化,例如分割边界的整体或者一部分可以在一个周期里扩大,而在下一个周期里缩小。
[0016] 在将确定演化方向趋势图的步骤和更新水平集函数的步骤进行重复的步骤中的重复之后的确定演化方向趋势图包括改变一个或多个趋势方向中的至少一个,这种改变可能涉及将趋势方向转变为其反方向,例如改变其值的正负符号。
[0017] 将确定演化方向趋势图的步骤和更新水平集函数的步骤进行重复的步骤中的重复涉及将满足图像点收敛条件的图像点排除在更新之外。演化方向趋势图的应用使得探测局部图像点的收敛成为可能,而在随后重复的演化方向趋势图计算和水平集函数更新的步骤中,这些的图像点就可以被跳过。这也对加快运行速度有很大帮助。
[0018] 图像点的收敛条件可以被定义为在前一个或前几个周期内,在相应图像点上水平集函数的数值变换小于一个预估或者用户设定的值。
[0019] 在某些基于本专利的应用中,将确定演化方向趋势图的步骤和更新水平集函数的步骤进行重复的步骤还包括对于水平集函数的所有图像点或满足图像点收敛条件的图像点将水平集函数设定到基于最近的水平集函数值的值,例如当前和先前水平集函数值的平均值。
[0020] 所述计算水平集函数演化方向趋势图和更新水平集函数的步骤可以运行一次或者重复多次直至满足一个水平集函数收敛的条件。
[0021] 满足水平集函数收敛的条件可以是指满足图像点收敛条件的图像点个数达到一定的数量(最好是全部)。水平集函数收敛的条件也可以是周期性将确定演化方向趋势图的步骤和更新水平集函数的步骤进行重复的次数。
[0022] 当满足图像点收敛条件的图像点个数达到一个预计或者用户设定的数量时,则可视为水平集函数收敛。或者是周期性将确定演化方向趋势图的步骤和更新水平集函数的步骤进行重复达到一个预计的或者用户设定的次数时,也可以视为水平集函数收敛。或者当以上两个条件都满足时,则视为水平集函数收敛。
[0023] 本发明所述将表示图像的输入数字图像数据处理为输出数字图像数据的过程可以是一个图像分割的过程。其中,输出数字图像数据表示由输入数字图像数据表示的图像中的物体的边界。水平集函数的更新可以局限于图像的一部分,例如在临近水平集的零水平的一个窄带上。
[0024] 在图像分割过程中,本发明所述更新水平集函数的步骤还可以包括检查新的需要更新的图像点并将其加入到窄带上,这些新加入图像点的趋势方向可以从临近的已经在窄带上的图像点继承而来。从而避免对整个演化方向趋势图进行更新,却可以使得新加入的图像点得到一个趋势方向以便进行下一步更新水平集函数的步骤。
[0025] 另外一种方法是当新的图像点被加入到窄带上时,前述演化条件也可以被视为满足,从而引发重新计算水平集函数演化方向趋势图的步骤。这样也可使新加入的图像点被赋予相应的趋势方向以便进行下一步更新水平集函数的步骤。
[0026] 在图像分割过程中,本发明所述的方法还可以包括输入一个先验形状的步骤,和在更新水平集函数的步骤还包括基于先验形状计算变形的形状,并在更新水平集函数中利用变形的形状。因为传统的基于先验图形的图像分割方法运行速度相对较慢,所以本发明所产生的加速效果对于这些需要使用先验图形的应用就更有意义。在更新水平集函数的步骤中还可以包括用速度函数给出的速度与函数趋势方向不符的图像点来计算先验形状的形变参数。这样做的原因是这些速度与函数趋势方向不符的图像点更接近收敛状态,也就是说更接近物体的边缘,因此比其他远离物体边缘的图像点对于估计先验形状的形变来说更有意义。而那些远离物体边缘的图像点,即速度函数给出的速度与函数趋势方向一致的图像点,则在计算先验形状的形变时被忽略掉了。
[0027] 在前述更新水平集函数的步骤中,对先验形状进行变形的操作可以由速度方向与趋势方向相对关系发生变化的图像点的数量来决定。如果相较前一次迭代,图像点的速度与趋势方向相对关系没有发生明显变化,那么对先验形状进行变形的操作就可以跳过。而明显变化指的是一个或者多个图像点的速度与趋势方向相对关系发生了变化。至于多少个图像点发生变化可以定义为明显变化可以需要根据具体的案例来分析。
[0028] 本发明所述将表示图像的输入数字图像数据处理为输出数字图像数据的过程可以是一个图像降噪的过程。其中,输出数字图像数据表示由输入数字图像数据表示的图像的经过平滑和降噪后的图像。
[0029] 本发明所述方法可以被整合到一个可以加载的计算机软件内,这个软件可以调用所述的各个步骤。本发明所述方法还可以被整合到任何计算机可读取的介质里,从而指导计算机完成所述的各个步骤。本发明所述方法还可以整合到任何可预先配置的计算设备内,由该计算设备完成所述的各个步骤。本发明所述方法还可以被整合到医疗影像设备中,比如核磁共振成像设备和CT成像设备中,由这些设备中的计算机完成所述的各个步骤。

附图说明

[0030] 以下所列附图可以帮助进一步说明本发明的具体步骤、目标和优点。但本发明并不局限于这些附图所示范畴。
[0031] 图1a显示传统算法中水平集函数演化过程中扭动的示意图
[0032] 图1b在相同的情况下使用本发明所述方法的示意图
[0033] 图2显示传统方法每次迭代的运算次数并与基于本发明的一种实施方式相对比[0034] 图3a-c基于本发明的一种实施方式的伪代码
[0035] 图4a基于本发明的一种实施方式的流程图
[0036] 图4b与图4a做对比的基于传统水平集方法的图像处理流程图,
[0037] 图5a-d使用先验形状的图像分割的示意图
[0038] 图6a基于传统水平集方法并使用先验形状的流程图
[0039] 图6b基于本发明使用先验形状的一种实施方式的流程图
[0040] 图7本发明作为计算机软件应用的示意图
[0041] 图8医学影像采集设备与计算机的关系
[0042] 不同流程图中各个操作步骤所标注的引用标号可能相同,这表示即便是在不同的实施方式中,这些步骤的特征相同。

具体实施方式

[0043] 虽然本发明所述方法理论上可以应用于各种可以使用水平集方法的图像处理任务中,但是它在图形识别和攫取等应用中,尤其是图像分割的任务中优势最为明显,如果与窄带法相结合就更佳完美。
[0044] 下面就基于本发明的一种实施方式进行详细地说明:
[0045] 在水平集函数演化的过程中,我们可以观察到在某些点上水平集传播地比临近点上要快一些,而在下一次迭代的时候,由于曲率力的作用,这些点又会向相反的方向运动以保持水平集(通常是零水平集)的平滑。这种反复出现的扭动现象大大减慢了水平集的传播速度。针对这一观察结果我们在本发明中提出一种改进的稀疏场水平集方法。这一方法使用一种周期性的单调的速度函数,从而使水平集函数的传播变得具有一致性。在水平集方法中,分割的边界由一个符号距离函数φ来隐含地表示。在内部力和外部力的联合作用下这一函数按照一个速度函数vt迭代演化,从而使水平集函数的零水平集按待分割物体的形状变形。在很多基于水平集的应用中,曲率力(所谓的内在因子或称内部力)被用于保持分割边界的平滑,而图像的信号强度,梯度和纹理等信息(所谓的外在因子或称外部力)则被用来将分割边界吸引到待分割物体的边缘上。公式l给出了一种水平集函数应用的示例。
[0046] (公式1) 这里vt=αC(x,t)+(1-α)D(x)
[0047] 在公式1中C(x,t)表示局部曲率,D(x)则表示一种基于信号强度的外部力。由于图像噪声和计算误差等原因,零水平集(分割边界)上的相邻的图像点不会按照同样速度传播。即便是他们开始排列在同一条直线上,很快他们之间就会形成锯齿样的分布,参见图1a(a-c)。那些传播较快的点(例如图1a中点2和4)在下一次迭代中就会在曲率力的作用下被“拽”回。这种局部扭动的现象在整个水平集演化中会不断地出现,大大影响了水平集的传播速度。在本实施方式中,我们用一个周期性单调的函数来取代传统的速度函数,参见公式2。
[0048] (公式2) 这里 在公式2中,trend是一个趋势函数,表示局部的分割边界是在扩展还是收缩。当一个像素或者体素被加入到窄带上时,它的趋势函数值被初始化为该点速度函数vt的正负符号,即sign(vt),或者取 ,在此Nr表示距离该点小于r的临近区域。在实际操作中,在第一次迭代以后,新加入点的趋势函数值可以无须计算,而是从引起加入点改变的临近点那里继承过来。虽然所有点的趋势函数值最好是一样的,并不是必须一样。当所有在窄带上的点的w(vt)的函数值都变为零(也就是说vt×trend≤0)时一个周期结束,所有点的趋势函数值都要同时改变符号。如图1b(d-f)所示,当使用这样周期性单调的速度函数时,那些传播较快的点(点2和4)不会像图1a那样由于传播慢的临近点的作用而向反方向移动(因为vt×trend<0),而是在原地等待其他的临近点赶上。这样一来,分割界面又变得平坦了(图1b(f)),点2和4的速度函数的值与趋势函数值变成同号,所有点又开始共同向前方传播。趋势函数的初始化对于结果的影响并不大。由于水平集能够自动处理拓扑变化,即便由于噪声等原因某些点的趋势方向被设成与临近点的趋势方向相反,水平集函数还是可以绕过这些点后重新融合为分割界面。而这些噪声点则变成水平集函数中一个个孤岛,在下一个周期中被淹没。通过这种方式,分割界面上扩张的部分会不断的扩张,而收缩的部分会不断地收缩,从而达到呈一致性地演化,直到它们碰到物体的边界。由于外在因子的作用,当分割边界跨越物体的边缘时,速度函数会改变符号,那么这些区域就不可能再扩展或者收缩了。而当分割边界伸展入一个锐利的夹角里时,曲率力也会使速度函数改变符号。这时所有点上的w(vt)都变为零。在下一周期中,趋势函数会改变符号,从而允许分割边界向相反方向传播,来纠正在上一周期中跨越物体边缘所造成的误差。就这样,趋势函数会周期性地改变符号,而分割边界则慢慢收敛于物体的边缘上。
[0049] 基于上述一致性传播算法,我们又进一步提出了一种懒惰窄带算法来跳过在物体边缘上无休止的冗余的计算。当一部分分割边界经过几次迭代的演化,跨越物体边缘时,在余下的迭代中可以只计算那些还在扩展或者收缩的部分。用公式1在每个像素或体素上求解水平集函数的计算量非常大。在稀疏场的水平集算法中,只有在零水平集经过的像素或体素上水平集函数的值会被更新,而其临近的点上(2-3个像素或者体素距离内)水平集函数值则是通过距离变换来更新。这样就使得运算量下降了一个数量级。但是在传统的稀疏场算法中,每次迭代中窄带上的所有点都要被更新,即便是分割边界的某些部分已经收敛于物体的边缘。而使用上述一致性传播方式时,在某一个周期内当分割边界的某些部分触及物体的边缘时,w(vt)就会变为零,并且如果其周围的点不发生改变的话,那么在余下的迭代中,w(vt)就一直为零,这些点便可以被标记为“休眠”状态,而不必再更新其水平集函数值。但是由于这些点只是处于一个暂时的“稳态”(w(vt)=0),须要有一种叫醒机制,以便当一个点的水平集函数值发生变化时唤醒周围1到2个像素或体素范围内的“休眠”点。并且最好是只唤醒那些和这些改变点具有相同趋势方向的点,因为这意味着临近点已经赶上了传播较快的点,否则则表示临近点更远离传播较快的点,w(vt)还是保持为零。在一个周期结束后,所有处于休眠的点又都被唤醒。当这些点都很接近物体的边缘时,它们会在有限的几次迭代后重新进入休眠状态。
[0050] 我们进一步提出了一种水平集函数收敛检测的方法,来停止对已经稳定于某些像素或体素上几个周期的零水平集区域进行更新,并且在所有的点都变得稳定时让整个分割过程自动地停止下来。在大多数稀疏场的水平集方法中,由于重构水平集函数(距离变换)和解偏微分方程过程中的计算误差,速度函数vt很少会真正地收敛于零,而是在零的周围进行小幅地跳跃(一般步长小于一个像素)。而当使用周期性单调速度函数时,分割边界虽然在多数情况下会像传统方法那样收敛,但在某些区域,尤其是在尖锐拐角处曲率力量相对较大的部分分割边界会以较大的步长来回跳跃。这是因为在两边关节节点上的一个小的跳跃使得这个弹性区域产生很大的形变。一种解决方案是,在进行了几个周期的一致性传播后,转换为传统算法。另外一种更为有效的方法是,当几个周期内某些点上水平集函数值的变化变得足够小时就永久性地停止对该区域进行更新。一个扭动计数器可以被用来记录在某个点上已经有多少个周期水平集函数φ甲的变化量小于一个给定值(比如1.0)。在每个周期末,如果在窄带上的点的水平集函数值的变化足够小,其扭动计数就会增加一。当一个点的扭动计数超过一个给定阈值(例如10),这个点就被标记为永久休眠,所有这个点上的计算将停止。而这个点上最终的水平集函数值则可以由计算最后两个周期末函数值的平均值来获得。
[0051] 由此可见本发明的方法可以通过在几个方面跳过不必要的计算并且加速收敛速度来对水平集方法进行加速。例如,一致性的传播方式可以使得零水平集更快地接近物体的边缘。这不仅是因为传播较快的点不再向相反方向运动以等待较慢的临近点,而且每次迭代的步长也可以变大,因为在新方法中Δt可以被设得更大。在所有的水平集算法中,为了保持水平集函数演化的稳定性,一般要求所有零水平集上的点满足Δt×vt小于一个像素或体素单位(否则零水平集的连续性就会被破坏)。为了选择一个统一的Δt来对所有点的水平集函数进行更新,一般算法会选择每次迭代中所有窄带点中最大的vt来计算Δt,以使Δt小于1/|vt-max|。在传统水平集方法中那些向相反方向扭动的点往往有较高的速率,Δt往往因为这些点而受到限制。在本发明的算法中,这些点的Δt被设为零,因此最大速度的绝对值,|vt-max|就会变小,从而使其他的点可以经过更少的迭代次数完成一个像素单位的移动。此外,与传统的窄带或者稀疏场水平集方法不同,上述的懒惰窄带法不仅可以避免远离零水平集的图像点上的计算,而且可以跳过那些已经收敛部分的分割边界上的运算,集中力量于那些活动的变化区域。
[0052] 图2比较了两种方法的每次迭代中具体参与计算的图像点的个数,这里2表示本发明介绍的方法,而1表示传统的稀疏场算法,两者使用了同样的实验数据和参数。分割的准确性对扭动计数器阈值的选择并不敏感。在实验中,在处理二维图像时我们将这一阈值设为10,而在处理三维数据时我们将之设为5。
[0053] 我们对二维和三维的版本在一台使用Intel Xeon处理器的Mac Pro的计算机上进行了实现。并设定这些程序为单线程。速度函数中vt的D(x)项被设定为一个基于阈值的函数:D(x)=ε-(I(x)-T),这里I(x)是输入图像的信号强度,T凡是阈值窗的中心而ε是阈值窗的宽度。在对二维合成图像和临床影像的测试中,基于本发明的算法较传统稀疏场算法要快5到15倍。分割一个512×512的图像一般需要0.1至0.2秒。而在三维数据的测试中我们测试了12个腹主动脉的CT血管造影的图像。每个病例都标记了一个种子区域。据我们的实验,本发明的方法较传统方法提速10至20倍。处理一个256×256×256的三维体数据大概需要5到12秒的时间。同时我们测量了两种方法窄带点之间的距离,从而比较新方法和传统方法分割结果的差别。由于在传统算法中分割边界倾向于在小血管中提早结束增长,我们在进行比较时,只比较了腹主动脉主干上的点。在上述12例病例中,两种算法的零水平集之间的距离在-1.80至1.46体素单位内变化(负值表示新方法的结果的边界在传统方法边界的外侧),平均距离为-0.12,标准方差为0.11。如果不考虑距离的正负号,即绝对距离,平均距离为0.12(标准方差为0.10,但其并不呈正态分布)。
[0054] 除了以上介绍的周期性单调水平集的实现方式以外,在一些步骤上的变种也值得一提。一种是将趋势函数提升到一个尺度空间,也就是说针对于每个像素的趋势方向trend(x,t)是由一个半径为r的区域中所有临近像素计算而来,而不仅仅是当前像素。显然,当r=1时,这一方法将返回和上述讨论的实现方式相同的结果。在另外一个极端的情况下,r可取无穷大,或者至少是能涵盖整个图像的所有像素,这样使在窄带上所有点的收缩或扩张都同步起来,也就是每个点的trend(x,t)都相同,因此就不必为每个像素分配额外的内存空间。这种变种使得基于GPU的运算更容易实现,因为在GPU上内存往往是有限的。同时它也使得懒惰窄带算法更简化,因为如果分割边界上收缩和扩张同时进行,那么在处理连接收缩和扩展区域的关节区域时就要格外小心,以避免水平集更新过程在这些节点上进入死循环而无法被标记为休眠状态。实验表明,这种同步的周期性单调水平集方法所给出的结果和上面讨论的实现方式是一致的。
[0055] 在上述懒惰窄带算法中,所有像素上trend(x,t)的改变是同步进行的。换句话说,所有点的周期变化是同步的。但是并不是所有的情况下都要要求已经触及物体边缘的像素来等待远处那些还在传播的像素。一种非同步的版本也是可行的。在这个变种里,休眠的像素可以在所有在s大小的临近区域内的所有点都进入休眠时就进入下一个周期,而不必考虑远处距离大于s还在传播的像素还处于上一个周期中。这种“局部周期性单调的水平集”的演化进程比上面讨论的“周期性单调的水平集”更接近传统水平集方法。
[0056] 另外一种变种是用一种阻尼因子的方法让分割边界在振荡了几个周期后慢慢地减速并最终停止下来,这种机制将取代上面提到的扭曲计数器的机制。一种实现这种阻尼式停止的方法是让Δt小于sp/|vt-max|,这里sp是周期p中的最大步长,并且sp随着周期数增加而呈线性或指数性减少,且so≤1.0。整个运算可以在几个周期后停止,具体周期的个数则可由用户来设定,和/或在sp减少到一个预设或者用户设定的阈值时终止运算。专业人士可以根据具体的应用相当容易的选择运算周期的个数或者sp的阈值。
[0057] 水平集方法的另一应用领域是图像降噪。在这一应用中的通行方法包括,全变分法(参见Osher,S.,and Fedkiw,R.P.“, Level set methods and dynamic implicit surfaces”,Springer,2003)和各向异性扩散滤波(参见G.Gerig,O.Kubler,R.Kikinis,and F.A.Jolesz“,Nonlinear anisotropic filtering of MRI data,”Medical Imaging,IEEE Transactions on11,221–232,2002)。在降噪的应用里,输入的图像被看做是初始化的水平集。与图像分割时不同,速度函数被用于减少所有水平的函数曲率,而不仅仅是零水平集。输入的图像同时还被用作外部力来防止图像在经过一段时间以后变得过度平滑。例如,在公式1中,如果设定C(x,t)为φt的局部曲率,D(x,t)则可以被设成φt-I,并且φ0=I,这里I是输入的图像。这样在水平集演化过程中曲率被减少的同时也保持输入和输出图像间的差别较小。本发明涉及的方法可以直接被用于降噪的应用中。与上述图像分割的实现方式不同的是,所有这些计算演化趋势和更新水平集函数的操作要在所有的图像点上进行,而不是仅限于一个窄带上。
[0058] 需要指出的是,公式1中的速度函数只是一个用来解释如何进行图像分割和图像降噪的示例,但是本发明并不仅限于这一系列的速度函数。其它的相关文献中记载的速度函数也都可以应用于本发明。更多的速度函数的例子可以在《Level set methods and dynamic implicit surfaces》中的第11、12章中找到。
[0059] 图3a-b中的伪代码对一种与上述方法相类似的实现方式进行了诠释。这些伪代码可以帮助相关专业人士深入而有效地了解本发明中所述的方法。图3a中显示的是伪代码的主程序,图3b-c中显示的是主程序调用的两个函数。在此,我们不对这些伪代码进行过多的解释,因为对于专业人士这些伪代码本身已经非常明了,而且前面已经对整体思路做了解析。
[0060] 图4a中的流程示意图展示了一种基于本发明的水平集方法图像处理实施方案。在步骤100,200和300中,输入图像数据,水平集函数φ和速度函数vt被输入到图像处理系统中。这些步骤可以与传统的水平集方法的最初几个步骤一致。输入图像数据一般是通过读取影像文件,或者通过用户人机交互等方式,但是在一些集成系统中也可能是通过一个图像采集装置采集而来。这里的水平集函数φ和速度函数vt一般是事先根据具体图像处理任务而设定或者估计的,但是它们中的一个或者两个函数可以是通过用户人机交互方式输入的,比如通过用户选择或其它的交互方式。在步骤200之后,或者至少是在更新水平集函数之前,水平集函数需要先被初始化,比如通过计算或者按照预设或用户给定数值等方式。
[0061] 在步骤400中,水平集函数的图像点的演化趋势方向被计算出来。这个趋势函数可以涵盖图像中所有的图像点,或者只是图像中有水平集函数定义的那部分图像点,比如在一个窄带上。前面提到的trend(x,t)就是趋势函数的一个例子。这个趋势函数一般包含了每个图像点的演化趋势方向,因此最好以演化方向趋势图的形式表示,而这个图上每个点的值就表示了相对应图像点的演化趋势方向。在某些实现方式中,这个演化方向趋势图可以由用户人机交互来初始化,或者通过估算来初始化。一般由用户人机交互来初始化更有优势,因为这样避免了一次初始化运算。这里的用户人机交互操作可以包括用户在待分割物体的内部或者外部指定一个种子点或者区域,从而保证演化趋势一致地向外扩张或者向内收缩。在图像分割的一些实例中,初始的演化趋势是可以推测的,比如当初始的分割边界被设定为图像的边缘时我们可以推断演化趋势一定是向内收缩的。演化趋势函数可以通过计算速度函数的值来推测。这种推测操作可以只用于初始化阶段,或者只用于初始化以后的演化方向趋势图的计算。在某些实现方式中,演化趋势函数可以由计算速度函数值的正负符号而得到,例如前面详细讨论的那个实现方式和伪代码展示的那个实现方式。每个图像点的演化趋势方向可以由这一个图像点或者其临近区域内的多个图像点的速度函数计算而来。在后一种情况下,这一图像点的趋势方向可以通过临近区域内所有的图像点的速度函数值的和来计算,如上面给出的伪代码中的实施方式所述。这里的方向可以只是所在点的当前速度函数值的正负符号。临近区域一般指的是距离当前点一定距离之内的所有图像点。这个区域可以只有一个图像点,也可以大到包含水平集函数上所有的图像点。在后一种情况下,趋势函数将只有同一个值。在某些实现方式中,演化方向趋势图可以是统一的,也就是所有图像点的趋势方向在同一周期内相同。而在另外一种实现方式中,趋势函数可以只是局部一致。在演化方向趋势图局部一致的情况下,这个值可以被存储在一个变量中而表示这个区域的演化趋势。演化方向趋势图可以用来判断当前的传播方向是否符合演化趋势,比如像上面讨论的实现方式中将演化方向表示为1或-1。
[0062] 在步骤500中,所有速度方向与趋势一致的图像点的水平集函数值得到更新,而那些所有速度方向与趋势不一致的图像点则被跳过。在那些更新的图像点上,更新过程可以与传统的水平集方法的更新过程一样使用速度函数。在某些实现方式中,这个更新水平集的操作可以是针对所有的点,也就是说包括那些已经收敛的点或者那些趋势方向与速度方向不符合的点。比如在某一实现方式中,在用速度函数和趋势函数更新水平集函数时,可以让那些速度方向和趋势方向不一致的点上的变化为零。上面的w(vt)函数就是一个这样的例子,这样做的好处是可以避免水平集传播过程中的扭动现象。但是,更有效的一种方式是直接跳过对收敛点和那些趋势方向与速度方向不符合的点的更新。这样可以进一步减少运算,让算法运行更快。判断速度方向和趋势方向是否一致的方法有很多,可以根据情况由专业人士进行选择。一种简单有效的方法是前面提到的用1或-1来表示每个点的可能的两个方向(比如扩张或收缩的趋势)。这个趋势方向可以直接与速度值相乘,其结果如果大于零则表示速度方向和趋势方向一致,水平集函数在这一图像点上得以更新,如果结果小于等于零,则水平集函数的更新被跳过。另外一种方法则如上述伪代码所示,趋势方向可以直接由正负号表示,直接与速度函数的符号进行对比。在前述的实现方式中,那些速度方向与趋势方向不一致的图像点被设成休眠状态。在某些实现方式中,在图4a所示步骤500中,如果某些已经被跳过的点恰好临近某个刚刚被更新过的图像点,那么这些点也会被更新。这在前面讨论的实现方式里被称为唤醒机制。这些被设置为暂时休眠的点就会被重新设成活动状态。
[0063] 还是在图4a中,步骤520表示检查是否演化条件已经满足。如果否,则上述步骤500中更新水平集函数的操作就会继续迭代。
[0064] 这里演化条件可以是一定数量的点在步骤500被设置为休眠和/或已经进行了一定次数的迭代。例如在伪代码表示的例子里,演化条件是所有的图像点已经处于“休眠”或者“停止”的状态。如果演化条件满足则意味着一个周期的结束,进入下一个步骤600。
[0065] 步骤600是一个检查水平集函数收敛条件的操作。如果这个条件满足,则意味着水平集函数已经整体收敛。如果不满足,则开始一个新的更新周期,重复步骤400来计算演化方向趋势图和步骤500来更新水平集函数。这与上面讨论的过程一样,只不过这次在新的周期里使用了新的演化方向趋势图。某些实现方式可能只包括一个更新周期。而在多个周期的实现方式中,后续的周期中可能只有那些还没有收敛的图像点会参与运算。步骤620体现了这种排除收敛点的操作。当一个图像点满足收敛条件时,该图像点被称为收敛点。在某些实现方式中,这种排除收敛点的操作并不是必须的。但是通常跳过这些收敛点可以使运算效率更高,因为这些收敛点上的运算一般是没有必要的。这里的收敛条件可以被设成在最后的几个周期中某个图像点的水平集函数值变化低于某个收敛阈值(一般很小)。而每个周期的结束是由上述演化条件来决定。这里的收敛阈值可以预先估算并由专业人士根据具体情况进行调节。同时,包括具体要观察周期的个数也可以由专业人士根据具体情况进行调节,这个数值一般在5至15之间。在前面讨论的实现方式中,我们使用了所谓的扭动计数器来检查某个图像点是否满足收敛条件。同时,在伪代码的例子中,当Δφ小于收敛阈值时,这一点的扭动计数器就加一,否则就被重置为零。当这个扭动计数器的值超过n时,也就是说在之前的n个周期中水平集函数在这一图像点上的变化小于收敛阈值,这一点上的水平集函数被视为收敛。而这一点上的水平集函数的值则可以由最近的几个周期末的函数值来计算,比如取最后两个周期末的平均值,如上述伪代码所示。
[0066] 如上所述,水平集函数按照图4中参照号900所包含的几个步骤进行迭代,直到水平集收敛条件得以满足。水平集收敛条件可以被设为一定数量的图像点(最好是全部)满足图像点收敛条件。在伪代码所示的实现方式里,主程序中的while-do循环对应的是水平集收敛条件;这个循环一直持续至所有的图像点都变得稳定,也就是说,所有图像点都满足图像点收敛条件。
[0067] 当水平集函数收敛后,当前的水平集函数就被用于步骤700中,来产生输出图像数据。具体如何产生输出图像数据根据具体应用的不同而不同。在传统的方法里这个步骤一般就是对收敛的水平集函数的简单处理。专业人士应能体会到,这里讨论的方法主要是输入和输出之间的步骤,如图4a所示,这里讨论的方法主要是参照号900所标注的几个步骤,而之前和后继的步骤则未作详细讨论。图4b里显示了传统的水平集方法的流程图,以与图4a进行对照。在图4b中,与图4a相对应的步骤包括:10和100,20和200,30和300,40和400,60和600,70和700。一般来讲,这里讨论的方法与传统方法中的输入和输出步骤可以是一样的,所不同的是核心的处理步骤,也就是图4b中参照号为90所包含的几个步骤,和图4a中参照号为900所包含的几个步骤。但如前面讨论的,由于本发明的方法是基于传统方法的改进,其中的一些步骤还是有其相似之处,比如用速度函数来对水平集函数进行更新的步骤(步骤60和600),虽然在两种方法中什么时候和在什么地方进行水平集的更新是不同的。可以说本发明是对传统算法的一种改进,并为其增加了一些特性。
[0068] 一般来说,对于那些用于图像分割的实现方式来说,输出的图像数据表示的是输入图像中待分割物体的边缘,像传统方法一样,这种情况下的水平集函数可以只在一个零水平集周围的窄带上进行更新。对于这些基于窄带水平集方法的实现方式,在更新水平集函数的步骤中,例如图4a中的步骤500,需要加入判断是否有新的点加入到窄带上来、对新加入的点的趋势方向进行推测和从临近的图像点继承趋势方向等几个步骤。所谓临近的图像点是指那些其水平集函数值刚刚发生变化的点,并且这一改变致使当前点被加入到窄带中来。在另外一些用到窄带法的实现方式中,在对水平集更新过程中(图4a步骤500),当一个新的图像点被加入到窄带中来时,可以认为演化条件得以满足,如图4a中步骤520所示。在这些变种中,这些新加入点的趋势方向可以经由重新计算演化方向趋势图的方法得以赋值。在接下来的步骤中水平集函数可以继续用新的趋势图来进行演化,重新回到图4a步骤
500。当然如果窄带外的图像点的趋势方向已经在加入窄带前就有定义,那么以上的这些操作就没有必要了。
[0069] 一般而言,在图像降噪的实现方式中,输出的图像数据代表的是对输入图像进行平滑和降噪处理后的影像。
[0070] 水平集方法的一个优势是,它可以利用曲率力的作用来防止分割边界泄露现象的发生。然而在某些情况下,两个信号强度相似的物体距离较近时,两者间相连接的边缘可能较大。如在图5a所示的例子里,一个矩形和一个圆形区域相重叠。在这种情况下,传统的水平集方法一般会给出错误的分割结果。Leventon等人在“Statistical shape influence in geodesic active contours,in Computer Vision and Pattern Recognition,2000.Proceedings.IEEE Conference on,vol.1,pp.316–323,2002”中提出了一种基于统计形状模型的水平集方法来解决这一问题。这里的统计形状模型,也被称为先验形状,是通过一组形状的平均形状(u)和他们的变异因子的主成分(a)来表示的。在每次水平集函数迭代中,这一统计形状模型被变形为最为符合当前物体的形状(u*|a,p;p是刚体变换参数,包含平移、旋转和防缩等),这里的形变参数可以用最大后验概率的方法来估计。这个变形的统计形状模型随后被用来限制水平集函数的演化,使后者不至于离当前物体的边缘(u*)太远。图5a-d给出了一个示例。在这个例子中,一个矩形和一个圆形重叠形成了图形3,但是我们的目标是希望能够应用先验知识分割出矩形。图5b显示了一个初始的水平集的种子区域
4和一个先验形状5。图5c显示了经过几次迭代后的运行结果。图5d显示了最终结果,这里的收敛的水平集函数4与先验形状5非常接近。
[0071] 人们还提出了很多基于先验形状的水平集方法的变种,比如Chen等人发表的“Using prior shapes in geometric active contours in a variational framework,International Journal of Computer Vision,vol.50,no.3,pp.315–328,2002”和Tsai等人发表的“A shape-based approach to the segmentation of medical imagery using level sets,IEEE Transactions on Medical Imaging,vol.22,no.2,pp.137-154,2003”。虽然这些方法使用了不同的先验形状、形状配准技术和含有先验形状项的速度函数,但是其总体工作流程是一样的。
[0072] 图6a给出了传统的基于先验形状的水平集方法的流程图。这一流程是在图4b的基础上演化而来的。以下我们只讨论两种流程不同的部分。首先在步骤35中输入先验形状。这个先验形状以及输入的图像数据和水平集函数随后被用在步骤45中,来计算符合当前物体的形变参数。这个变形后的形状被用在步骤50中来指导水平集函数的更新。
[0073] 基于先验形状的水平集函数的一个缺点是它的运行速度非常慢,因为在每次迭代中都有进行一次形状配准的操作。即便是使用窄带方法,只用窄带上的点来进行形状配准,在处理三维医学影像数据时,这种方法的运行速度还是很难让人满意。本发明所述方法则可以很好地解决这一问题,在本发明中所述的水平集方法中加入先验形状也并不复杂。因为本发明所述方法可以加速水平集函数的收敛,所以即便没有做任何针对先验形状方法的优化,直接套用本发明的方法来解决使用先验形状的水平集问题,也可以起到加速的作用。这包括以下一些步骤:在步骤35中提供一个先验形状,在步骤45中这个先验形状被变形为符合当前水平集函数和输入图像中目标物体的形状。这个变形后的先验形状随后被用在步骤50中来指导水平集函数的更新。
[0074] 图6b显示了将本发明应用于使用先验形状水平集方法的一种实现方式的流程图。这个流程图是从图4a修改而来的。为了方便,我们以下只讨论两者之间不同的部分。(这些不同之处与图6a中较一般水平集方法增加的步骤相类似)。在步骤350中输入先验形状,在步骤450中利用输入图像和水平集函数对先验形状进行形变。变形后的先验形状被用在步骤500中来指导水平集函数的更新。
[0075] 本发明所述使用先验形状的方法还可以引入新的优化手段。在传统的使用先验形状的方法中,所有在窄带上的点都被用作图像配准,也就是计算先验形状的形变参数,并且在每次迭代中都要进行重复的配准,例如图6a所示流程图中的循环部分。但是在本发明所述方法中,通过使用演化方向趋势图,我们可以检测出那些触及物体边缘后已经收敛的图像点(当他们的速度方向和趋势方向不一致时)。因此我们可以只利用这些速度方向和趋势方向不同的图像点来计算形变参数,而跳过其他的图像点。这样不仅通过跳过那些还在传播的图像点而减少了运算次数,而且去除了那些还在不断改变的图像点的干扰,增加了算法的稳定性。
[0076] 另外一种加速的手段是通过分析速度方向改变的图像点的个数来决定是否需要进行形变参数的更新,这里速度方向改变的图像点指的是相比上一次迭代的步骤500中速度方向与趋势方向的相对关系发生改变的图像点。例如要求必须有10%以上的图像点的速度方向发生改变时才对形变参数进行重新计算。可以预见,在开始阶段,形变参数的重新计算会较频繁,随着越来越多的图像点满足收敛条件,其频率会逐渐下降。结果表明,在接近结束时,先验形状与当前分割的结果已经非常吻合了。
[0077] 如图7所示,本发明可以被实现于某种计算机程序中,该计算机程序可以被调入到计算机内存中,而这段计算机程序在该计算机上运行时要包含能够控制和调用本发明所述步骤的相关指令。例如所述计算机程序可以作为可执行文件1003被存储于硬盘或是通其它方式存储1004,如存储于服务器上并通过网络1005来下载到某台计算机1006上。该计算机可以执行该程序,但也可以是作为一个文件传输的一个中转环节。该计算机程序还可以被存储于存储卡1001或存储盘片上1002,如CD或DVD。这里的存储卡1001和存储盘片1002都是来举例说明那些用来记录所述计算机程序的计算机存储媒介,而该计算机程序可以指导计算机运行本发明所述步骤。
[0078] 图8给出了一个计算机1101和一个医学影像采集设备1102的示意图。该影像采集设备可以是核磁共振扫描机或CT机等。本发明可以通过设置计算机1101来执行所述步骤来实现。这可能涉及上面讨论的计算机程序被调入到计算机内存中,而该计算机程序在该计算机上运行时要包含能够控制和调用本发明所述步骤的相关指令。该计算机程序可以是被记录在某种可以为计算机读取的媒介中。而该计算机1101则可以通整合等方式成为医学影像采集设备的一部分。
[0079] 所有上述对于本发明的解释和说明,包括附图,都是为了方便理解而做的举例,不应该被理解为本发明仅限于上述范畴和形式。本发明的实现方式并不限于这里讨论的几种方式。例如,可以用于本发明的速度函数就有很多不同的形势,不同的速度函数可以参考相关文献,几乎所有文献报道的速度函数都可以应用于本发明中。
[0080] 本发明由权利要求书中的权利要求来定义。上述关于实现方式的讨论是为了帮助专业人士更有效地实现本专利所述方法。在权利要求中用到的“包括”字眼并不意味着将其它成分或者步骤排除在外,所用“某个”等字眼并不排除所饰名词复数的可能。如果某些特征出现在不同的相对独立的权利要求项目中,这并不代表他们互相排斥,这些特征的组合也是可能的。权利要求中各种不同方法的表述顺序并不代表任何意义。某种方法内部不同步骤的顺序要根据具体语境来解释,并不一定表示这些步骤要按照表述的顺序执行。对于步骤执行顺序的唯一限制是步骤之间的依赖性。权利要求中出现的参照号标记只是为了方便理解,不应视为对相应权利要求的限制。