对3D数字医学图像中的解剖实体进行分割的方法

摘要

对于图像中的一定数目的地标当中的每一个定义该地标的初始位置。接下来，对包括所述地标的一定数目的候选位置的所述初始位置周围的邻域进行采样，并且把成本与每一个所述候选位置相关联。对表示所有候选位置的总体灰度成本和总体形状成本的加权和的成本函数进行优化。把已分割解剖实体定义为经过候选位置的所选组合的一条路径，其中针对所述组合优化了所述成本函数。在朝向最优已分割表面/体积的优化期间利用图遍历方法。

权利要求

1.一种用于对数字医学图像中的解剖实体进行分割的方法，其包括以下步骤：-对于所述图像中的一定数目的地标当中的每一个地标定义所述地标的初始位置；

-计算所有地标点的图数据结构，其中每一个节点(顶点)代表一个地标，并且每一条边(弧线)代表一对地标之间的连接；

-按照系统方式或随机方式计算经过所述图数据结构的路径；

-对所述初始位置周围的邻域进行采样，所述邻域包括所述地标的一定数目的候选位置；

-把成本与每一个所述候选位置相关联；

-对于所有候选位置优化表示总体灰度成本与总体形状成本的加权和的成本函数；

-通过每一个地标的所述候选位置的所选组合把已分割解剖实体计算为表面或体积，其中对于所述组合优化所述成本函数，访问地标的顺序由经过所述图数据结构的所述路径确定，所述路径是最小跨度树(MST)。

2.根据权利要求1的方法，其中，与候选位置相关联的所述成本代表所述候选位置处的采样灰度图样与所述地标处的预期灰度图样之间的相似度的度量，所述相似度由Mahalanobis距离表示，Mahalanobis距离被定义为所述采样灰度图样与利用逆协方差矩阵加权的均值灰度图样之间的距离。

3.根据权利要求2的方法，其中，与候选位置相关联的所述成本是通过对为所述医学图像计算的特征图像集合当中的每一个特征图像中的所述候选位置的灰度成本求和而获得的总灰度成本，所述灰度成本被表示为Mahalanobis距离，Mahalanobis距离被定义为所述候选位置处的灰度值轮廓与通过相应的协方差矩阵的逆矩阵加权的相应的均值轮廓之间的距离，所述均值轮廓和协方差矩阵是从与所述解剖实体相关联的灰度值模型获取的。

4.根据权利要求3的方法，其中，所述灰度值模型是通过以下步骤获得的：-在一定数目的地标点处对所述解剖实体的人工分割的表面或体积进行采样；

-通过分别连接多边形或多面体内的成组的相邻地标点而形成所述表面的多边形网格或所述体积的多面体网格；

-定义所有地标点的图数据结构，其中每一个节点(顶点)代表一个地标，每一条边(弧线)代表一对地标之间的连接；

-在每一个所述地标点的邻域内采样一定数目的点；

-对于每一个地标点，把在地标点周围的邻域内采样的点设置在轮廓内；

-计算至少一个特征图像；

-计算每一个地标点的以及每一个特征图像的均值轮廓；

-计算每一个地标点和每一个特征图像的轮廓的协方差矩阵；

-把所述图、均值轮廓和协方差矩阵标识为所述解剖实体的灰度值模型。

5.根据权利要求4的方法，其中，所述多边形是三角形。

6.根据权利要求4的方法，其中，所述多面体是四面体。

7.根据权利要求3的方法，其中，特征图像是预定义的尺度σ下的导数图像，或者是包括地标位置周围的宽度为α的窗口内的所述图像或导数图像的局部直方图的数学矩的图像表示。

8.根据权利要求1的方法，其中，所述总体灰度成本是一定数目的特征图像中的轮廓的所有单独灰度成本的组合，所述单独成本是Mahalanobis距离，Mahalanobis距离被定义为所述轮廓与利用所述协方差矩阵的逆矩阵加权的那个特征的均值轮廓之间的距离，所述均值轮廓和协方差矩阵是从所述解剖实体的灰度值模型获取的。

9.根据权利要求8的方法，其中，所述灰度值模型是通过以下步骤获得的：-在一定数目的地标点处对所述解剖实体的人工分割的表面或体积进行采样；

-通过分别连接多边形或多面体内的成组的相邻地标点而形成所述表面的多边形网格或所述体积的多面体网格；

-定义所有地标点的图数据结构，其中每一个节点(顶点)代表一个地标，每一条边(弧线)代表一对地标之间的连接；

-在每一个所述地标点的邻域内采样一定数目的点；

-对于每一个地标点把在地标点周围的邻域内采样的点设置在轮廓内；

-计算至少一个特征图像；

-计算每一个地标点的以及每一个特征图像的均值轮廓；

-计算每一个地标点和每一个特征图像的轮廓的协方差矩阵；

-把所述图、均值轮廓和协方差矩阵标识为所述解剖实体的灰度值模型。

10.根据权利要求1的方法，其中，所述总体形状成本是所有地标点的连接矢量的所有单独成本的组合，所述连接矢量连接相继地标，所述单独形状成本由Mahalanobis距离表示，Mahalanobis距离被定义为两个相继地标之间的连接矢量与利用所述协方差矩阵的逆矩阵加权的均值连接矢量的距离，所述均值连接矢量和协方差矩阵是从所述解剖实体的形状模型获取的。

11.根据权利要求10的方法，其中，所述形状模型是通过以下步骤获得的：-在一定数目的地标点处对所述解剖实体的人工分割的表面或体积进行采样；

-通过分别连接多边形或多面体内的成组的相邻地标点而形成所述表面的多边形网格或所述体积的多面体网格；

-定义所有地标点的图数据结构，其中每一个节点(顶点)代表一个地标，每一条边(弧线)代表一对地标之间的连接；

-计算成对的地标点之间的连接矢量；

-计算成对的地标点之间的均值连接矢量；

-计算所述图中的每一个连接的连接矢量的协方差矩阵；

-把所述图、均值连接矢量和协方差矩阵标识为所述解剖实体的几何模型。

12.根据权利要求1的方法，其中，通过动态编程优化所述成本函数。

13.根据权利要求1的方法，其中，通过选择具有最小总灰度值成本的那些候选位置而减少所述候选位置的数目，所述总灰度值成本是在所有特征图像上所述邻域内的所述候选位置的所有灰度值成本的总和，所述灰度值成本被表示为Mahalanobis距离，Mahalanobis距离被定义为所述候选位置处的灰度值轮廓与利用协方差矩阵的逆矩阵加权的相应的均值轮廓的距离，所述均值轮廓和协方差矩阵是从与所述解剖实体相关联的灰度值模型获取的。

14.根据权利要求1的方法，其中，所述邻域由采样点的矩形栅格构成。

15.根据权利要求1、4或11的方法，其中，所述灰度值模型和所述形状模型是从所述图像的多分辨率表示构造的并且被应用于所述多分辨率表示。

16.一种在数字医学图像中测量解剖实体的方法，其包括以下步骤：-利用根据权利要求1的方法计算对所述解剖实体的分割；

-基于对所述解剖实体的所述分割计算或选择特性点；

-基于所述特性点计算测量对象；

-基于所述测量对象导出对所述解剖实体的测量。

17.一种在数字医学图像中测量解剖实体的方法，其包括以下步骤：-利用根据权利要求1的方法计算对所述解剖实体的分割；

-把所述地标标识为测量点；

-基于所述测量点计算测量对象；

-基于所述测量对象导出对所述解剖实体的测量。

18.根据权利要求1的方法，其改进在于，通过广度优先遍历获得经过所述图数据结构的所述路径。

19.根据权利要求16的方法，其改进在于，通过广度优先遍历获得经过所述图数据结构的所述路径。

20.根据权利要求17的方法，其改进在于，通过广度优先遍历获得经过所述图数据结构的所述路径。

21.一种计算机程序产品，当运行在计算机上时，其适于实施权利要求1、16、17或18的方法。

22.一种计算机可读介质，其包括适于实施权利要求1、16、17或18的步骤的计算机可执行程序代码。

说明书

技术领域

本发明涉及一种对数字图像中的实体进行分割的方法，更具体来说涉及对数字医学图像中的解剖实体进行分割的方法。根据本发明的分割处理是基于对图像的几何模型及光度模型的使用。

本发明的方法特别可以用作在数字图像上执行几何测量的处理，更具体来说是在数字医学图像上执行几何测量的处理。

背景技术

在放射学实践中，几何测量被频繁地用于帮助诊断异常。为了执这些测量，必须把关键用户点放置在图像(例如在显示设备上显示的图像)中其相应的解剖地标位置处。诸如两点之间的距离或者线之间的角度之类的测量就是基于这些关键用户点的位置。此外还可以在整体上评估几何结构的常态或异常，包括对完整形状的分析。因此需要自动且客观地提取出嵌入在放射学图像中的定量信息。
所述频繁执行的测量的一个例子是计算胸廓RX图像(图6)中的心胸比值(CTR)。这一比值被定义为心脏在顶点水平的横直径与胸廓的内直径(ID)的比值，即CTR＝(MLD+MRD)/ID。
心脏的横直径由心脏左侧最大横直径(MLD)和心脏右侧最大横直径(MRD)构成。这一定义显然要求放射科医师沿着左、右肋廓的内边界进行搜索，以便定位计算所述内直径ID所需要的点对。该点对必须位于胸廓的最大内直径处。同样地，必须搜索左、右心影边界以便定位计算MLD与MRD的和所需要的点。更具体来说，这些点处于相对于脊柱中线的最远处。所述边界搜索处理要求放射科医师执行解剖结构分割并且定位已分割解剖结构上的各点(在本例中一共是四个点)。在CTR计算的情况下，所述分割步骤等于是描绘肺野。
3D应用的一个例子是对肺进行分割以及计算肺容积。
数字图像中的许多其他测量遵循类似的方法，其中涉及到解剖器官或实体的分割，在分割的几何结构上确定特征点和测量对象并且最后执行测量(图5)。
参照心胸指数计算的例子，为了自动定位所需各点，需要一种在胸腔射线照相图像上自动对肺野进行分割的方法。
取决于应用，可以通过几种方式解决所述分割问题。分割策略已经从早年间的计算机视觉的低级别策略发展到更近的基于模型的策略。
低级别方法依赖于局部图像算子把具有不同光度特性的像素分开并且把具有类似的局部光度特性的像素分组在一起。这两种类别的例子是边缘检测和区域生长。尽管这些低级别方法的性能差，但是它们在大多数商用图像分析工具中仍然非常受欢迎。其主要原因是这些低级别方法易于理解及实施。但是对于例如在医学图像中呈现并且以上述胸廓图像的内容为例的复杂图像数据来说，其效用就受到限制。
更为成功的方法结合了关于将要分割的形状以及关于图像中的对象的光度或灰度外观的先验知识。这些方法被称作基于模型的方法，其通常以模板匹配为基础。例如通过进行相关或者利用广义Hough变换技术来匹配模板。不幸的是，所述模板匹配在医学图像的情况中往往会失败。这是由于解剖对象可能会在形状和灰度外观方面表现出很大的变异性。
基于由Kass等人引入的主动外形(M.Kass、A.Witkin和D.Terzopoulos的“Snakes：active contour models(蛇：主动外形模型)”，Int.J.Computer Vision，1(4)：321-331，1988)以及水平集(J.A.Sethian的“Level set methods and fast marching methods(水平集方法和快速行进方法)”，Cambridge Univ.Press，Cambridge，U.K.，1999)的方法能够应对更大的形状变异性，但是仍然不适用于许多医学分割任务，这是因为只能结合关于待分割对象的很少的先验知识。
手工(handcrafted)参数模型克服了这一问题，但是其受限于单一应用。
鉴于上述缺点，很明显的是需要一种可以利用实例进行训练以便获得关于待分割对象的形状以及图像中的对象的灰度外观的知识的一般性分割方案。由Cootes和Taylor引入的主动形状模型(ASM)(T.F.Cootes、C.J.Taylor、D.Cooper、J.Graham的“Active Shape Models-theirtraining and applications(主动形状模型——其训练及应用)”，ComputerVision and Image Understanding，61(1)：38-59，1995)满足这种分割方案定义。形状模型由地标点的矢量的主分量给出。所述灰度外观模型描述了以每一个地标为中心的与对象外形(contour)垂直的轮廓(profile)的归一化一阶导数的统计。通过最小化一阶导数轮廓与轮廓的分布之间的Mahalanobis距离找到地标在新图像中的位置。这种算法开始于一个初始估计并且执行拟合程序，其是地标位移和形状模型拟合的交替。已经设想出类似的方法，这些方法全都采用一种三步程序。首先，所述方法都使用确保产生似真结果的形状模型。其次，所述方法使用灰度外观模型把对象放置在这样的位置处，在所述位置，边界周围或者对象内部的灰度图样与从训练实例所预期的类似。最后，所述算法通过最小化某一成本函数来拟合所述模型。
由主动形状模型所提出的方法仍然面对几方面的局限。
第一种局限是需要初始估计，这被称作初始定位问题。在某些情况下，在这种框架中需要针对适当的初始外形进行大量搜索。
第二种局限在于交替地使用形状模型和灰度外观模型。首先，利用所述灰度外观模型更新所述分割估计。其次，把所述形状模型拟合到所述更新后的外形。不幸的是，如果所述灰度外观模型错误地定位了地标，则会误导所述形状模型。
利用主动形状模型分割方案的一些实验表现出另一个问题。如果特定地标的灰度外观模型被用来找该地标的更好位置，则要求真实的地标位置处在所述算法所探索的区域内部。在某些情况下，所述算法可能在某一地标的区域内搜索另一个地标。这意味着在错误的区域内使用了错误的灰度外观模型，从而导致差的分割。
本发明的一个目的是提供一种对数字图像中的实体，更具体来说是数字医学图像中的解剖实体，进行分割的方法，该方法克服了现有技术的问题。

发明内容

上述各方面是通过如权利要求1所述的方法实现的。
在从属权利要求中阐述了本发明的优选实施例的具体特征。
本发明的另一方面涉及一种如在所附权利要求书中阐述的对数字医学图像中的解剖实体进行测量的方法。
下面所使用的术语灰度值模型、灰度值外观模型、灰度模型、灰度外观模型和光度模型都是同义词。
同样地，术语形状模型和几何模型被用作同义词。
术语特征图像和特征也被用作同义词。
本发明的方法的实施例通常以计算机程序产品的形式实现，当在计算机上运行时，所述计算机程序产品适于实施本发明的方法步骤。
所述计算机程序产品通常被存储在诸如DVD或CD-ROM的计算机可读载体介质中。或者，所述计算机程序产品可以具有电信号的形式，并且可以通过电子通信被传送给用户。
在EP A 1349098中公开了一种在数字采集的医学图像中自动进行测量的方法，该方法是通过把测量对象和实体分组到由双向链接的外部图形模型和内部信息模型的计算机化测量方案中而实现的。
在根据EP A 1349098的测量期内，从计算机获取测量方案并且将其激活。随后在所激活的测量方案的指导下在所显示的图像上执行测量。
在这种计算机化的方法中，在所述数字图像中映射多个几何对象，其他测量对象以及最终的测量实体(比如距离和角度)以所述几何对象为基础。基本几何对象通常是关键用户点，其定义作为几何测量的基础的其他几何对象。但是上面提到的专利申请没有公开如何能够在不需要用户定位及操纵关键测量点的情况下自动实施所述映射。根据本发明的方法可以被用来实施这一映射。
从如在已公布的欧洲专利申请EP 1760658和EP 1760659中所描述的二维情况出发概述所述方法，在所述欧洲专利申请中分割得到一个或更多外形，所述外形是开口的或者闭合的，并且围绕感兴趣的解剖对象。
根据本发明，提供一种分割方法，所述分割方法包括构造医学图像中的诸如器官或结构之类的解剖实体的模型，可以从实例训练所述模型，从而获得关于灰度值外观的知识以及关于待分割的解剖实体的形状的知识。
灰度值外观模型包括灰度值的概率分布以及在分布于解剖轮廓线(outline)上的一组地标处的灰度值的多个多阶导数。
形状模型包括沿着一个解剖实体的轮廓线的相继地标之间的每一个连接矢量的概率分布。
引入了基于离散点的对象表示来描述医学图像中的解剖轮廓线。
公开了用于生成模型的两种策略和相关联的系统以及应用所述模型来分割实际图像。这些策略虽然分解成具有相同目标的构建块，但是其构造并且采用不同的光度和几何模型知识。
首先，从训练图像构造新的灰度外观模型，编码地标位置处的光度知识。这一步骤利用了在每一个地标周围采样的邻域内的强度相关。
其次，构造形状模型，编码了各地标之间的几何知识。这一步骤利用了已分割对象的地标之间的空间相关。
在分割步骤中，联合使用光度和几何知识来分割新图像上的一个或更多解剖外形。
所得到的分割可以用来在图像上执行几何测量的过程中导出所期望的测量点的位置。
可以按照直接或间接的方式来定义测量点。
在直接方式下，所述地标的位置被用作各个测量点。
在间接方式下，通过一条曲线把从所述分割处理得到的地标互连，并且把该曲线上的特性点定义为测量点。
灰度外观模型
每一种建模系统中的第一步都是把所述地标的几何位置约束到位于与当前外形垂直的线性路径上，或者位于以所述地标为中心的稀疏栅格上。在图像中的地标点处的光度图样由在特征图像中采样的轮廓矢量基于在一定阶数下提取的图像导数捕获。所述轮廓采样可以沿着所述图像中的线性路径进行，或者可以在围绕所述地标点的圆形路径上进行。
特征图像
在对于图像中的所有图像点计算光度图样时，获得了特征图像，其尺寸与原始图像的尺寸相等。通常采用邻域运算来计算所述特征图像。对应于特定点的所述特征图像中的特征的值可以设置在特征矢量中。特征图像的计算通常包括两个步骤：(a)使用一组卷积内核的线性滤波步骤；以及(b)涉及到与局部图像统计量的计算相组合的非线性点运算的子步骤的后处理步骤。在步骤(b)中，可以单独执行所述子步骤中的一个。
不同的滤波器和不同类型的后处理导致不同类型的特征图像。举例来说，Laws构造了25个二维卷积内核作为类似诸如水平、边缘、点、波以及波纹之类的理想特征轮廓的矢量对的外积。所述后处理步骤包括非线性加窗，其中在更大的窗口中对已滤波值的绝对值求平均，以便获得图像纹理能量度量。Unser使用这样的滤波器组，该滤波器组被构造为对应于公知的变换(比如离散正弦变换(DST)、离散余弦变换(SCT)、离散偶正弦变换(DEST)、离散实偶傅立叶变换(DREFT)、离散实奇傅立叶变换(DROFT))的矢量的外积，并且用于对通道方差的计算进行后处理。最后，利用被调谐到不同空间频率和不同取向的组合的对称或反对称Gabor滤波器，构造Gabor滤波器组，使得所有输出都具有相同的频带宽度。所述非线性点运算步骤在这种情况下可以包括阈值处理。
可以利用基于多分辨率小波滤波器、秩值(rank-value)滤波器或自适应滤波器的特征图像来测量其他图像特性。
基于局部无序图像(LOI)的特征图像
医学图像中的解剖地标周围的图像结构表现出典型的灰度图样。在本发明中，该结构在后面概述的N个数学特征中被捕获。较早前在B.VanGinneken等人的“Active shape model segmentation with optimal features(利用最优特征的主动形状模型分割)”(IEEE Trans.on MedicalImaging，21(8)：924-933，2002年)中介绍了使用特征图像来建立灰度外观模型的概念，其是基于由J.J.Koenderink和A.J.Vandoorn在“Thestructure of locally orderless image(局部无序图像的结构)”(Int.J.ofComputer Vision，31(2)：159-168，1999年)中提出的所谓的局部无序图像(LOI)。
Taylor展开可以用来通过K次多项式来近似位置x0处的感兴趣点(即地标)周围的灰度函数f。每一项的系数由x0处的导数f(i)给出：
$f\left(x\right)\approx \underset{i=0}{\overset{K}{\Sigma }}\frac{1}{i!}{f}^{\left(n\right)}\left(x_0\right){(x-x_0)}^i$
该函数近似中的所有信息都被包含在所述导数f(i)中。类似地，可以通过包含原始图像的导数(L，Lx，Ly，Lxx，Lyy，Lxy，...)的图像滤波器组来描述所述图像。还可以在一定数目的模糊或内尺度σ下计算每一个导数图像。最终通过把所述后处理步骤计算为每一个导数图像的局部图像统计量而获得所述特征图像。所述局部图像统计量包括每一个位置x0周围的宽度为α的窗口内的局部直方图的一定数目的矩。在一个具体实施例中，计算一阶和二阶直方图矩、二次及以下的所有导数(L，Lx，Ly，Lxx，Lyy，Lxy)以及5个内尺度(σ＝0.5、1、2、4、8个像素)，从而得到一共N＝2×6×5＝60个特征图像。用来计算所述直方图的窗口尺度是根据所述内尺度，即α＝2σ。
直接使用所述特征矢量中的行、列(即分别是x、y)内的导数的局限是关于变换和旋转(即欧氏变换)缺乏不变性。除非把所述图像反旋转到标准取向，否则就只能在具有相同取向的实例上训练这些算子。笛卡尔微分不变式(CDI)独立于所选择的笛卡尔坐标系描述了图像的微分结构；因此可以省略所述反旋转步骤。在L.M.J.Florack等人的“Scaleand the differential structure of images(图像的尺度和微分结构)”(Imageand Vision Computing，Vol.10(6)：376-388，1992年)一文中描述了从高斯微分算子构造CDI。使用二次及以下导数的CDI如下：$(L,L_{\mathrm{xx}}+L_{\mathrm{yy}},L_x^2+L_y^2,L_x^2{L}_{\mathrm{xx}}+2L_{\mathrm{xy}}{L}_x{L}_y+L_y^2{L}_{\mathrm{yy}},L_{\mathrm{xx}}^2+2L_{\mathrm{xy}}^2+L_{\mathrm{yy}}^2).$同样地，在每一个像素位置处，通过在一定数目的内尺度σ下使用这些算子并且计算范围为α的窗口内的局部无序直方图的前两阶矩来构造各特征图像，其中α与所述内尺度相关联，例如α＝2σ。这些算子独立于图像几何结构的欧氏变换，但是仍然与尺度参数σ的选择有关。
轮廓和特征相似度测量
正如下面所解释的那样，所述灰度外观模型建立在平均轮廓和方差-协方差矩阵的基础上，并且利用了所述轮廓的点之间的强度相关。一轮廓处的强度由如上概述的强度特征矢量表征。
为了把一个轮廓与另一个轮廓做比较，必须知道所述轮廓中的每一个点与任意另一个点如何相关，并且还必须测量所述轮廓的特定点的值与任意其他点的均值的差异程度。这种统计关系由协方差度量所捕获，其中所述协方差度量总是在两个点(或者这些点处的特征)之间计算。如果在一点的特征值与其自身之间测量协方差，则获得方差。当所述轮廓在所述外形的每一侧包括k个点时，对于总共2k+1个点存在(2k+1)k个关系。所述轮廓的一对点之间的协方差表明该对中的第一值相对于第二值如何变化。如果所述协方差大并且为正，则表明当第一点的特征值增大时，第二点的特征值也按照成比例的方式增大。如果所述协方差大并且为负，则表明当第一点的特征值增大时，第二点的特征值按照成比例的方式减小。如果所述协方差为零，则意味着在所述特征值之间没有系统性关联，也就说它们彼此独立。可以把所有的方差-协方差对设置在一个协方差矩阵S中，其关于主对角线对称。因此，所述协方差矩阵捕获了所述轮廓的特征值对之间的所有结构相关性。在实践中，所述协方差矩阵是从学习样本建立的，在本例中所述学习样本是一定数目的地标处的轮廓的集合以及一组训练图像中的一定数目的强度特征。
所述协方差涉及到如何在更高维度的空间(即给定轮廓中的所述2k+1个点的空间)内计算距离，以便计算当前轮廓与模型轮廓之间的距离。当所述轮廓仅仅包括值为x的单一点时，可以仅仅通过把当前轮廓点与模型轮廓点的值相减并且取绝对值，即$d=|x-\bar{x}|,$来比较所述各轮廓与模型点值的相似度。当所述轮廓包括两个点时，可以把2维矢量x与模型2维矢量之间的距离计算为平面内的欧氏距离，即：
$d=\sqrt{{(x-\bar{x})}^{\prime }(x-\bar{x})}=\sqrt{{|x_1-{\bar{x}}_1|}^2+{|x_2-{\bar{x}}_2|}^2}.$
欧氏距离在所有方向上的权重相等。但是当x1和x2分别具有不相等的标准偏差σ1和σ2时，等距离线的外形是主轴与坐标轴平行的椭圆。可以用所述标准偏差的倒数对所述差进行加权，从而得到加权欧氏距离：
$d_w=\sqrt{{(x-\bar{x})}^{\prime }W(x-\bar{x})}=\sqrt{{\left|\frac{x_1-{\bar{x}}_1}{{\sigma }_1}\right|}^2+{\left|\frac{x_2-{\bar{x}}_2}{{\sigma }_2}\right|}^2},$
矩阵W包含变量的方差的倒数，即：
$W=\left[\begin{array}{c}1/{\sigma }_1^2\\ 1/{\sigma }_2^2\end{array}\right].$
如在线性或圆形轮廓的情况中，当彼此接近的像素属于具有平滑变化的强度水平的相同解剖区域时，所述像素的强度水平的相关性将非常高。当像素对被视为属于不同组织时，相关性可能成反比例。在任一种情况下，成像噪声以及不重要的解剖变化都会引入像素自身的强度水平的非零方差。
因此，如果变量也是相关的并且具有不同的方差，则必须插入所述协方差矩阵的倒数，从而得到所谓的Mahalanobis距离：
$d_m=\sqrt{{(x-\bar{x})}^{\prime }{S}^{-1}(x-\bar{x}).}$
所述Mahalanobis距离对多维数据点x与其均值之间的距离进行加权，从而使得处于相同的多元正态强度外形上的观察将具有相同的距离dm。这些Mahalanobis等距离外形在二维空间内具有椭圆的形式，并且在三维空间内具有椭球的形式。显而易见的是，所述轮廓中的点的数目2k+1通常远远多于3个，在这种情况下，所述等距离轨迹在多维空间内具有超椭球的形式。
所述灰度外观模型是从训练图像集获得的，并且其包括对应于每一个地标的平均轮廓以及对应于每一个所述特征的协方差矩阵。因此，所存储的模型的总数与地标的数目和特征的数目成比例。
形状模型
在形状建模步骤中，把沿着所述外形的地标之间的位置相关性纳入考虑。由于旨在分割具有相对确定的形状的解剖结构，因此关于所述形状的知识必须优选地编码确定性分量连同有关的解剖变化，并且还优选地忽略无关的形状细节。
概述了两个实施例以考虑位置相关性。第一个实施例基于涉及到同时使用所有地标位置的全局方法；第二个实施例采用一种局部方法，其对相继地标的空间布置进行建模。
PCA分析及拟合
在所述第一个实施例中，从一组对准的训练图像中的人工放置的地标构造地标坐标位置的方差-协方差矩阵。所述方差-协方差矩阵中的每一个对角项代表所述地标的坐标值的方差。非对角项代表地标坐标值对的协方差，也就是说，其表示两个地标的坐标相对于彼此如何改变。所述变化可以关于其均值系统地处于相等的“正方向”上从而产生高的正协方差，或者可以关于其均值系统地处于相等的“负方向”上从而产生高的负协方差，或者关于其均值随机地处于相反的方向上从而产生零协方差。在现有技术中称为主分量分析(PCA)的特征向量/特征值分析将检测出所述地标关于所述平均地标位置的任何全局的系统性变化模式。PCA旨在通过仅仅保持具有最大变化的坐标系(coordinate set)分量来降低数据集的维度；因此其在统计上将原始坐标位置去相关。通过投影到这些高度变化的子空间，可以通过保留重要形状特征的维度更少的系统来近似有关的统计量。数目足够多的所述主模式被保留并且被用作所述形状模型，忽略可以归因于噪声变化的其他模式。
地标连接矢量模型
所述第二个实施例旨在在局部水平对相继地标的布置进行建模。但是这一限制并不意味着不捕获全局形状。在给出起始点坐标和连接矢量的情况下，通过从所述起始点开始迭代地把所述连接矢量加到当前点上，可以完全重建形状外形线。例如通过链码在像素水平的形状描述中利用了这一属性。所述连接矢量是局部切线矢量对曲线的离散近似。在考虑相继地标之间的方向矢量(即一阶几何导数)时，其允许描述对象的形状而不管该对象的平移。通过考虑例如由两个相继连接矢量所成的角度表示的相继地标之间的曲率(即二阶几何微分)，可以描述所述对象的形状而不管所述形状的平移和旋转(所谓的刚体变换)。在这种情况下，当所述起始位置和一阶几何导数(即起始连接矢量)已知时，通过迭代地添加相继切线矢量的差矢量可以完全重建所述曲线。基于像素水平的链码的曲率通常与噪声高度相关，这是因为该尺度处于过细水平。这一缺陷在本公开内容中得到了解决，这是因为所述地标的间隔足够大，从而允许计算独立于噪声但是仍然精确的一阶或二阶导数信息。由于所述切线或曲率矢量的量值与尺度有关，因此例如通过把所述训练集合的每一条曲线的尺寸归一化到“一”还允许描述曲线达到欧氏相似变换。一个对准的训练图像集中的人工放置的地标之间的位移矢量(即一阶导数或切线)或者所述位移矢量的矢量差(即二阶导数或曲率)的均值和协方差统计量构成了本实施例的形状模型。
分割系统
还可以沿着分别由具有相关目的的构建块组成的两条轨迹来组织分割系统。
第一个构建块旨在把实例化模型地标位置的可允许位置限制到目标图像中的一定数目的最可能位置。
如果地标由特征图像中的线性轮廓表示，则可以采用对穿过实例化地标的垂直线上的一组候选位置进行排序的最佳点票选方法。
对于圆形轮廓的情况，可以在每一个实例化地标处考虑矩形搜索栅格，并且根据所有特征图像的圆形轮廓到所述模型圆形轮廓的组合Mahalanobis距离对栅格交叉位置排序。
第二个和第三个构建块联合构成优化步骤。可以利用与所保留的每一个地标的可允许位置相关联的成本度量构造成本矩阵。所述成本包括与灰度外观中的一致性相关联的项。可以利用加权PCA外形拟合在全局水平验证与所述形状模型的一致性。或者，可以通过合并与所述地标的候选位置相关联的形状成本来验证所述形状一致性。可以使用动态编程算法来构造穿过所述成本矩阵的最小成本路径，从而得到作为具有最低总体成本的路径的最终分割。
最后，所述最终分割上的特性点可以间接地定义关键测量点的几何位置，并且可以在这些关键测量点的基础上生成所有依赖性的测量对象和测量量。
或者，用在所述分割中的地标点可以用来直接定义所述关键测量点，并且可以为每一个测量点提供搜索栅格以便定位该测量点的可允许位置。
通过下面的描述和附图，本发明的其他具体实施例和相关联的优点将变得显而易见。

附图说明

图1是示出了根据本发明的第一种建模系统(建模系统I)的流程图；
图2是示出了根据本发明的第一种分割系统(分割系统I)的流程图；
图3是示出了根据本发明的第二种建模系统(建模系统II)的流程图；
图4是示出了根据本发明的第二种分割系统(分割系统II)的流程图；
图5是示出了基于根据本发明所分割的解剖对象上的特性点的测量计算的流程图；
图6是用于从胸腔射线照片测量心胸比的测量方案；
图7(a)示出了胸腔射线照片中左肺的人工描绘，在图7(b)中所述左肺被表示为一组地标在本例中n＝14；
图8示出了通过一组n个地标p1＝(x1，y1)，...，pn＝(xn，yn)来离散地表示肺野。所述图像尺寸被归一化到0...1之间。图中的连续曲线是通过所述离散地标内插得到的；
图9示出了作为分别对一条线和围绕所述地标的半径为rc的一个圆上的一定数目的点处的灰度值的采样的线性和圆形灰度轮廓；
图10示出了包括14个地标的左肺形状上的一些连接矢量的实验估计的分布：(a)v1，(b)v2，(c)v3；
图11示出了胸腔射线照片的一阶矩特征图像，对应于二次及以下导数(L，Lx，Ly，Lxx，Lyy，Lxy)以及等于0.5、2和8个像素的尺度σ；
图12示出了胸腔射线照片的二阶矩特征图像，对应于二次及以下导数(L，Lx，Ly，Lxx，Lyy，Lxy)以及等于0.5、2和8个像素的尺度σ。
图13示出了本发明应用于胸腔射线照片上的肺野分割的实例：(a)左、右肺野的人工分割；(b)左、右肺野的自动分割；
图14示出了(a)利用人工肺野分割的心胸比计算，其得到CTRman＝0.47；(b)利用根据本发明的自动导出的肺野分割的心胸比计算，其得到CTRautomatic＝0.45。

具体实施方式

下面将关于具体应用，即医学图像中的肺野分割，详细解释本发明。
对象表示
在下面描述的本发明的方法的各具体实施例中，图像中的解剖对象被数学地表示为位于围绕该对象的外形上的固定数目的离散的已标记点，即p1＝(x1，y1)，...，pn＝(xn，yn)。
所述外形从p1行进到pn又回到p1。因此，所述对象可以被离散形状矢量x＝(x1，y1，...，xn，yn)T捕获。坐标系被选择成使得所述图像区域内的所有点都位于[0，1]x[0，1]域内(图7)。
此外，位于由所述外形围绕的区域内的特性点也可以被添加到所述对象表示当中，从而例如允许测量位于所述对象内部的实体。
下面描述的分割方案需要一定数目的训练图像，在所述训练图像中所述形状矢量x被人工确定。一旦在所述数据集上训练了所述算法之后，其就能够在新图像中产生所述形状矢量。
可以通过从所述点进行曲线内插来重建解剖对象的连续外形。一阶内插用线段链接各相继点，从而得到多边形近似。更高阶的内插可以被用来获得更平滑的轮廓线。不管顺序如何，所述连续表示可以被用来导出所述曲线的空间属性，比如沿着所述曲线的切线方向或者所述曲线上的特定点处的法线方向。
在三维中，添加了第三个轴z，并且所述地标点现在分布在围绕所述对象的一个表面上，即p1＝(x，y1，z1)，...，pn＝(xn，yn，zn)。因此，在这种基于点的表示中，可以通过离散形状矢量x＝(x1，y1，z1，...，xn，yn，zn)T来捕获所述对象。
通过链接边形中的成组的表面点获得一阶连续表示，其中每一个多边形定义多边形网格的一个内插表面贴片，连接成使得每一条边被最多两个多边形共享。一条边连接两个顶点，并且多边形是边的闭合序列。当所述点被分组成三元组时，每一个三元组恰好定义一个平面表面贴片的方程，在这种情况下所述多边形化处理变成表面三角剖分。这种处理与现有技术中的方法不同，所述现有技术方法提取出等表面(iso-surface)并且对所述表面进行局部三角剖分(这些方法的实例包括在W.Lorensen、H.E.Cline的“Marching Cubes：a high resolution 3Dsurface construction algorithm(行进立方体：高分辨率3D表面构造算法)”(Computer Graphics，Vol.21，No.4，1987年7月，p.163-169)中介绍的行进立方体算法)以及被称作NOISE、WISE和SAGE算法的后来的细化。不同之处在于，这里所采用的3D对象表示不假设所述点具有相等的灰度值，因此假设所述对象形状或器官轮廓线不同于等灰度值表面，从而所述分割不能以等表面提取为基础。相反，将对全体地标进行几何及辐射建模，并且把所述模型应用在新颖的基于模型的3D分割方案中。
此外还可以把位于由所述表面围绕的体积内的特性地标点添加到所述对象表示中，从而(a)使得能够由器官内部的主导解剖特征影响及操纵分割结果，或者(b)例如允许测量位于所述对象内部的实体。现在可以构造多面体网格，其把所述表面上以及所述表面内部的点链接成多个多面体体积单元的排列。例如在3D肺部扫描中，这些附加的点可以包括血管树和支气管树的特性点(比如分叉点)。具体来说，所述多面体可以仅仅包括四个点，每一组形成一个四面体，并且所述多面体化现在变成四面体化。代表所述体积的四面体结构的外壳是三角剖分的表面。
建模系统I(图1)
灰度外观模型
所述灰度外观模型(也称作光度模型)捕获所谓的地标(也称作地标点)的邻域内的空间强度分布。已经与灰度轮廓的使用一同设想了所述灰度外观模型的空间组成部分。从其构造可以明显看出，也可以采用其他方案。
例如被用来精确地描述肺野轮廓线的地标数目n可以低至14个，例如地标1表示每一个肺野的最高点，地标7被定位在左心影或右心影与其对应的半膈的横截面处，地标9表示肋膈角。其他地标可以等距离地定位在这些主要地标之间。很明显，可以使用数目更多的地标。
-线性灰度轮廓(图9)
在每一个地标处，垂直于所述外形采样描述围绕该地标的典型图像结构的灰度外观模型。在所述地标的每一侧，利用特定步长采样k个像素，这给出长度为2k+1的轮廓矢量。这种采样方案的优点在于，其能够对所述地标处的线性强度梯度进行建模。所述轮廓可以被归一化，从而使得所述轮廓中的元素的绝对值之和为1。采样的方向可以被计算为通过将所述地标与其前一个或下一个相邻地标相连而得到的线段上的法线之间的平均法线方向，或者可以被计算为连接所述前一个与下一个地标的线段上的法线。
在三维中，所述法线方向的计算更为复杂，因为需要找到表示为分布点集的表面的法线。在三维中，所述法线是关于平面定义的。该平面可以利用欧氏距离作为接近度标准，通过一组至少三个所选的最近相邻点拟合一个平面ax+by+cz+d＝0来获得。所得到的平面的法线具有分量(ai，bi，ci)并且被用作方向矢量，沿着所述方向矢量在特定地标(xi，yi，zi)处垂直于所述对象表面对所述体积进行采样，以便获得3D体素水平轮廓。
当把所述体积内的地标点添加到所述模型当中时，可以例如沿着所述地标点与其周围的位于所述表面内部或者位于所述表面上的地标点的连接矢量中的一些，在所述体积内对轮廓进行采样。
-灰度特征图像
地标周围的图像结构所表现出的被捕获在一定数目的N个数学特征中，正如前面所解释的那样。通过计算每一个位置x0周围的宽度为α的窗口内的导数图像的局部直方图的一定数目的矩而获得所述特征图像。在一个优选实施例中，计算一阶和二阶直方图矩、二阶及以下的所有导数(L，Lx，Ly，Lxx，Lyy，Lxy)以及5个内尺度(σ＝0.5、1、2、4、8个像素)，从而得到一共N＝2×6×5＝60个特征图像。用来计算所述直方图的窗口尺度是根据所述内尺度，即α＝2σ。所述特征图像集合可以被视为地标的邻域内的灰度函数的灰度分解。
在图11(三个尺度下的原始图像和5个导数图像的一阶直方图矩)和图12(三个尺度下的原始图像和5个导数图像的二阶直方图矩)中给出了胸廓图像中的所计算的特征图像的一个实例。胸廓图像通常出现在所述图像中的右上位置处；关于其他射线照相检查的解剖图像可能会表现出非标准化的位置和旋转，在这种情况下，特征图像的构造可以采用笛卡尔微分不变式。
在三维中，所述导数是3D算子。可以通过包含原始图像的部分导数(L，Lx，Ly，Lz，Lxx，Lyy，Lzz，Lxy，Lyz，Lxz，...)的图像来局部地描述所述体积数据。与2D情况类似，可以在一定数目的模糊尺度σ下计算所述导数，并且最终利用所述局部直方图的矩将所述导数用来计算一定数目的特征图像。在一个优选实施例中，计算一阶和二阶直方图矩、二阶及以下的所有导数(L，Lx，Ly，Lz，Lxx，Lyy，Lzz，Lxy，Lyz，Lxz)以及5个内尺度(σ＝0.5、1、2、4、8个像素)，从而得到一共N＝2×10×5＝100个特征图像。用来计算所述直方图的窗口尺度同样是根据所述内尺度，即α＝2σ。
-灰度外观模型
一个训练集合的所有图像被用来建立对应于每一个特征和每一个地标的统计模型。通过把在特征j(例如对于2D图像一共是60个)中的地标i(例如对于胸廓肺野的2D的情况一共是14个)处采样的归一化轮廓表示为矢量gi，j，可以从所述训练图像估计gi，j的概率分布，从而得到均值和协方差Si，j。对于2k+1的轮廓矢量长度，该协方差矩阵的尺寸为(2k+1)x(2k+1)。特定地标i的灰度外观模型包括所述均值轮廓(j＝1...N)和协方差矩阵Si，j(j＝1...N)。整个肺野的灰度外观模型包括所有均值轮廓(i＝1...n，j＝1...N)和协方差矩阵Si，j(i＝1...n，j＝1...N)的集合。
在三维中，在所述N个特征图像当中的每一个特征图像内，获得所述n个3D地标当中的每一个3D地标的体素水平采样(也称作“轮廓”)。所述地标的3D统计模型包括所有均值轮廓(i＝1...n，j＝1...N)和协方差矩阵Si，j(i＝1...n，j＝1...N)的集合。
-灰度成本
在对应于特定地标i的一个新特征图像j当中的一点p处采样的新轮廓gi，j之间的Mahalanobis距离被如下计算：
$h_{i,j}\left(p\right)={(g_{i,j}\left(p\right)-{\bar{g}}_{i,j})}^T{S}_{i,j}^{-1}(g_{i,j}\left(p\right)-{\bar{g}}_{i,j})$
更小的Mahalanobis距离意味着所述轮廓gi，j(p)源自具有均值和协方差Si，j的高斯分布的概率更大。因此，所述Mahalanobis距离可以被用作灰度成本度量，其被表示为hi，j(p)。该成本度量是一个给定图像中的位置p的函数。根据特征j，最有可能是地标i的真实位置的位置是使得hi，j(p)最小化的点p。因此可以为每一个特征j定义灰度成本。
在三维中，同样地把特定3D点p处的成本函数估计为在对应于特定3D地标i的新3D特征图像j中的3D点p处采样的新3D体素水平轮廓gi，j与所述均值之间的Mahalanobis距离hi，j(p)。
形状模型
所述曲线轮廓线，在肺野的情况下：每一个所述肺野有对应的一条，由n个地标点描述(图8)。在一组s个训练图像中人工地确定这些地标点，从而产生一个点序列(x1，y1)...(xn，yn)。这些坐标元组随后被放置在代表所述曲线的矢量x＝(x1，y1，...，xn，yn)T中。接下来对所述训练图像的形状矢量x应用主分量分析(PCA)。所述PCA把曲线投影到较低维空间中，涵盖了所述肺野的最重要的几何变化模式。可以用b∈□t来近似每一条曲线x∈□2n，其中t□2n：
$x\approx \bar{x}+\mathrm{\Phi b}$
其中，x是均值形状，并且Φ∈□2n×t是对应于t个最大特征值的x的协方差矩阵的特征向量。每一个特征值决定该形状模式的形状的方差。由矢量b表示的所述曲线近似构成所述形状模型，并且将所述曲线x拟合到该形状模型中。所述特征形状可以被视为由所述地标集合所代表的形状的零阶(位置)几何分解。所述形状模型的作用是约束所述对象在由所述训练集合施加的变化界限之间的形变，例如关于所述均值形状的三个标准偏差。
在三维中，通过n个3D点(地标)来描述所述表面轮廓线(即对应于每一个解剖对象有一条轮廓线)。所述轮廓线例如可以代表所述3D肺野。可以在一组s个训练图像中人工确定这些地标，或者利用自举技术从所述训练图像的先前子集s’获得这些地标，从而得到一个3D点序列(x1，y1，z1)...(xn，yn，zn)。这些坐标三元组随后被放置在代表所述表面的矢量x＝(x1，y1，z1，...，xn，yn，zn)T中。如前所述，由所述表面所描绘的所述体积内部的附加特性点可以补充所述形状模型。
接下来对每一个3D训练图像中覆盖所述(多个)对象的表面的3D地标集合应用主分量分析。所述PCA分析把所述表面投影到较低维空间中，涵盖了所述(多个)解剖对象(例如所述3D肺野)的最重要的几何变化模式。可以用b∈□t来近似每一个表面x∈□3n，其中t□3n：
$x\approx \bar{x}+\mathrm{\Phi b}$
其中，x是均值形状，并且Φ∈□3n×t是对应于t个最大特征值的x的协方差矩阵的特征向量。每一个特征值决定该形状模式的形状的方差。由矢量b表示的该表面近似构成所述形状模型，并且将所述表面x拟合到该形状模型中。与2D情况类似，所述特征形状可以被视为由所述地标集合所代表的表面形状的零阶(位置)几何分解。所述表面形状模型的作用是约束所述对象在由所述训练集合施加的变化界限之间的形变，例如关于所述均值形状的三个标准偏差。
分割系统I(图2)
所述分割算法最初将把代表肺野轮廓线的曲线放置在所述图像中，并且通过迭代地更新该曲线的组成地标的位置一直到优化标准达到最优值，而找到分割解决方案。每一次迭代包括以下步骤。
步骤1、最佳点票选
假定特定地标的当前位置，将在垂直于外形的线上搜索该地标的更好位置。在所述地标的当前位置的每一侧估计ns个位置。2ns+1个点当中的一个将被票选为最佳的下一点。所述一组N个特征当中的每一个特征为所述2ns+1个点之一贡献一票。根据所述Mahalanobis距离标准，选择所述2ns+1个点当中的具有最低成本的那点。对于所述2ns+1个点当中的每一个，包括在特征图像j中被采样的2k+1个点的轮廓被放置在填有适当的均值和协方差Si，j的hi，j(p)的方程中。具有最小Mahalanobis距离的点被该特征票选。每一个特征选择一个点，从而在2ns+1个点上分N票，从而得到集合，其中${\Sigma }_{k=1}^{2n_s+1}{N}_k=N.$显而易见的是，所述特征的数目N必须显著大于可选点的数目2ns+1，否则所能分配给可选点集合当中的点的票数就会过少。
特定地标i的所考虑的2ns+1个点所构成的集合当中的每一个点所接收到的票数Ni，j，可以视为该特定点是所搜索的地标的置信率。
在三维中，取通过三维特征票选最多的点，而获得2ns+1个3D点当中的选择。
步骤2、最小成本路径优化
根据第一步，每一个地标可以分别地在2ns+1个可能的位置上移位。所述置信率Ni，j可以如下放置在矩阵中：
$R=\left[\begin{array}{ccc}{N}_{\mathrm{1,1}}& ...& N_{1,{2n}_s+1}\\ .& & .\\ .& N_{i,j}& .\\ .& & .\\ N_{n,1}& ...& N_{n,{2n}_s+1}\end{array}\right]$其中$\underset{j=1}{\overset{2n_s+1}{\Sigma }}{N}_{i,j}=N$$\forall i=1,...,n.$
对所述外形进行更新涉及到在所述矩阵C的每一行当中选择一个元素；形成经过所述矩阵的路径使得成本标准被优化(即最小化)的所有元素的串连就是最小成本路径(MCP)搜索。
一种符合逻辑的选择将是取每一行当中的具有最高置信度的元素；这种方法的缺点在于，其没有考虑到离群值(outlier)的存在。因此，该方法在所述路径中可能会包括突然改变。这种缺点可以通过利用排除所述突然改变的平滑度约束来克服。
为此目的，首先如下构造成本矩阵：
$C=\left[\begin{array}{ccc}1/N_{\mathrm{1,1}}^m& ...& 1/N_{\mathrm{1,2}{n}_s+1}^m\\ .& & .\\ .& 1/N_{i,j}^m& .\\ .& & .\\ 1/N_{n,1}^m& ...& 1/N_{n,2n_s+1}^m\end{array}\right],$
其中，分母中的幂m是正数，其会影响所计算的路径(更大的m对点产生更小的成本贡献，因此把所述路径朝向该点前进)。其次，如下构造成本函数：
$J(j_1,j_2,...,j_n)=\underset{i=1}{\overset{n}{\Sigma }}{C}_{i,j_i}+\alpha \underset{i=1}{\overset{n}{\Sigma }}|{\delta }_{i,j_i}-{\delta }_{i-1,j_{i-1}}|,$
其中，是点i沿着所述形状的法线朝向其目标点ji的位移，α是权重因数。该等式中的第二项是所述平滑度约束，其对这样的目标点配置进行惩罚，对于所述目标点配置，一个目标点的位移相对于前一个目标点的位移沿着所述外形突然改变。例如利用现有技术中称为动态编程的技术对所述成本函数J进行优化，从而得到把J最小化为最小成本路径(MCP)的目标点(j1*，...，jn*)，这种技术例如在D.H.Ballard、C.M.Brown的“Computer Vision(计算机视觉)”(Englewood Cliffs，Prentice HallInc.，1982年，pp 137-143)中做了描述。
作为所述矩阵数据结构的替换方案，还可以使用分层图。每一层代表一个地标处的垂直采样的图像点的阵列。利用从最后一层回到第一层的一组弧扩增所述图。所述动态编程技术提取边界作为经过所述图的最优闭合路径。
从下文中可以明显看出，在三维中求解所述最小成本路径的问题更为复杂，并且无法从2D情况直接导出。首先，在所述2D外形检测方法中构造的矩阵数据结构无法被扩展用于3D情况，这是因为没有明显的方式来进行再循环(在2D情况下这是通过在所述矩阵中的最后一行之后拷贝第一行而实现的)。此外，所述2D矩阵方法实现了访问所有地标(每一个对应于一行)的系统方式；一种直接的方式是通过串连相继地标而得到的顺序的次序。在3D中，每一个地标具有跨越在所述3D表面上的一定数目的最近相邻地标，为了相等地对待这些最近相邻地标，必须并行地对其进行访问。其次，用于最小化的变量数目要多很多，这是因为3D图像中的对象的表面所包括的点通常大大多于2D图像中的对象的外形所包括的点。因此，所述最小化算法可能变得更容易固定在局部最小值处，这是因为无法穷尽估计所有的位置组合。概述用来应对这一问题的仿真退火优化策略。
下面概述适用于3D表面优化的两种优化策略。
最小成本表面优化
1、加权无向图数据结构G通过从所述三角剖分的表面开始连接所述表面上的所有地标点来创建，并且有可能利用朝向所述表面内部的点的连接来扩增，从而形成总体的多面体网格。每一个节点(顶点)g包括点的1D阵列以及所述点的相关联的成本(例如反梯度)。朝向所述表面的内部和外部垂直于所述表面对所述点阵列中的点进行采样。每一个节点被链接到通过从该节点出发遵循边(弧)a而获得的所有其相邻节点。为了避免循环，不再访问先前已经访问过的任何节点。可以利用以下对所述图的边进行加权：(a)在通过所述边链接的点对之间所计算的欧氏距离；(b)一地标及其相邻地标的几何属性(比如局部曲率)。所得到的图是稀疏的，这是因为并不是在所有非相邻地标点之间都存在边。此外还可以使用备选的表示，比如邻近矩阵表示，其中一对点之间的现有邻近关系的权重被作为条目放置在所述矩阵中；所有非相邻点都接收不存在的权重。所述图的另一种备选表示是邻近列表，其更适用于所述图稀疏的情况，即从表面三角剖分创建图的情况。在这里，对于每一个地标保持一个链接列表，存储所述列表的后续单元中的所有相邻点；可以提供这些列表之间的交叉链接以便在算法上易于删除节点。
每一个都包含1D点阵列的节点的集合可以被视为围绕当前表面的点的外壳。必须对于每一个节点计算由所述阵列中的所计算的新位置构成的新表面。这一处理在下面的步骤中概述，且可以视为可变形表面分割。
2、从所述图中的节点g开始，按照系统的方式访问所述图中的所有其他节点。在关于比如R.Sedgewick的“Algorithms in C++(C++算法)”(Addison-Wesley，1992年))的算法的标准教科书中描述的图算法可以被应用于这一目的。
a、最小跨度树(MST)
利用现有技术中称为最小跨度树的方法解决这一步骤。加权图的MST是连接所有节点的边的集合，其中所述边的互连权重总和至少与任何其他连接所有节点的边的集合的权重总和一样小。
这一问题的解决方案基于图论的属性，即假设把图中的节点任意划分成两个集合，所述MST包含把一个集合中的节点连接到另一个集合中的节点的最短边。
从任一节点(顶点)开始建立所述MST，并且所述MST总是取“最靠近”已取节点的下一个节点(得到欧氏MST，但是极值性可以是根据其他标准)。换句话说，在把已经存在于所述树中的节点连接到尚未处于所述树中的节点的所有边当中遵循具有最低权重的边，并且把该边所导向的节点添加到所述树中。上面提到的最小跨度树的属性确保所添加的每一条边是所述最小跨度树的一部分。所得到的MST不需要是唯一的。
找到所述MST的实际方法在现有技术中已知为Kruskal方法或Prim算法。
b、广度优先遍历
系统地访问所述图中的节点的另一种方式在现有技术中已知为从特定节点开始的广度优先遍历。广度优先搜索开始于覆盖所述表面的给定点节点周围的区域，并且只有在已经看了靠近的一切时才在所述表面上向更远处移动。这种图遍历可以被视为从所述起始点传播开去，以便逐渐包含所述表面的所有部分。对于闭合表面(即没有边界的表面)，传播前沿的各部分将在所述表面的背面彼此相撞。
如果所有的权重在所述图G中都是1，则所述广度优先遍历树将导致最短路径跨度树。
一般来说，MST树将不同于所述广度优先树。
虽然每一个节点与其相邻节点之间的距离会改变，因此必须重新计算所述MST或广度优先树，但是这些树在下一步骤的动态编程算法的所有迭代期间可以保持恒定。应当注意到，虽然初始的图通常将不是平面的，但是所得到的访问MST和广度优先树则是平面的。
按照这种方式，所分割表面的所有节点及其互补内地标按照系统的方式被链接在一起。
3、现在将利用每一个采样点中的梯度以及一节点的位移与相邻节点的位移的差，来对所述MST或广度优先遍历树执行所述动态编程(DP)方法，即：
$J(j_1,j_2,...,j_n)=\underset{i=1}{\overset{n}{\Sigma }}{C}_{i,j_i}+\alpha \underset{i=1}{\overset{n}{\Sigma }}|{\delta }_{i,j_i}-{\delta }_{i-1,j_{i-1}}|$
其中，n是所述树中的节点(点)数目，是点i沿着所述表面的法线朝向其目标点ji的位移，是所述树中的亲代节点点的位移，并且α是权重因数。节点i和i-1代表所述表面上的相邻点；ji和ji-1是分别与节点i和i-1相关联的点阵列中的索引。该等式中的第二项是平滑度约束，其对这样的目标点配置进行惩罚，对于所述目标点配置，在所述表面上一个目标点的位移相对于前一个目标点的位移突然改变。例如利用现有技术中称为动态编程的技术对所述成本函数J进行优化，从而得到把J最小化为最小成本路径(MCP)的目标点(j1*，...，jn*)。在把内部点添加到所述模型当中时，对应于这些节点的位移可以沿着与周围节点的连接矢量之一。
从为所述地标获得的位置开始，现在可以从所述图中的另一个节点g开始迭代步骤2和3并且由此生长新的MST以及计算新的地标位置，直到满足停止标准，所述停止标准例如是两次迭代之间的所有节点的绝对位移和不再显著改变。
仿真退火优化
在现有技术中已经由S.Kirckpatrick、C.D.Gelatt、M.P.Vecchi于“Optimization by simulated annealing(通过仿真退火进行优化)”(Science，Vol.220，No.4598，1983年5月，p.671-680)一文中引入了仿真退火。下面概述关于根据本发明的实现方式的细节。在这种优化策略中，不需要从所述图构造树以按照系统的方式访问所述表面上的点。对于当前表面上的地标位置的每一个集合，可以如下从所述图G计算与之相关联的能量：
$E=\underset{i=1}{\overset{n}{\Sigma }}(C_{i,j_i}+\alpha \underset{n_i\in N_i}{\Sigma }|{\delta }_{i,j_i}-{\delta }_{n_i,j_{n_i}}\left|\right).$
所述能量是必须被最小化的目标函数。其由与节点i相关联的一项以及与节点i的邻域节点集合Ni相关联的一项构成。例如当所述点具有较高的局部梯度时，所述节点项分配较低的能量。例如当每一个相邻节点的位移相对于被施加到当前节点的位移较小时，所述邻域项分配较低的能量，这意味着表面形状发生形变从而产生局部突然改变时对所述表面形状进行惩罚。所述仿真退火优化包括在最优表面的计算中应用Metropolis算法：当新迭代的能量E′低于前一次迭代的能量E时，所述能量总是取对应于较低能量状态的选项(下坡步骤)，而当能量E′高于能量E时有时采取上坡步骤。能量配置发生这种改变的概率采用自热力学系统，并且遵循所谓的Bolzmann概率分布p＝exp[-(E′-E)/kT]。控制参数T(对温度的模拟)和退火程序(annealing schedule)告知如何把所述能量E从高值降低到低值，例如在多少次随机配置改变之后采取T中的每一个向下的步骤以及该步骤有多大。与前面概述的基于树的方法相反，并不按照系统的方式访问所述表面上的点；相反，在所述表面上随机取所述点。SA算法的步骤如下：
1、把T设置到一个大值。
2、(配置步骤)把所述表面初始化为点集(x1，y1，z1)...(xn，yn，zn)并且计算其能量E。
3、(重新排列步骤)通过沿着所述地标垂直于所述表面的方向生成所述地标的新位置(x1′，y1′，z1′)...(xn′，yn′，zn′)而令所述表面发生形变。所述搜索算法可以选择每次把一个随机选择的地标移位特定步长，或者同时移位所有地标。通过追踪每一点发生移动的指向(sense)(表面向内或表面向外)并且在其处于较低能量的方向上时继续，可以获得高效的移动集合。
4、(目标函数步骤)计算所得到的表面的能量E′。
如果E′＜E则接受所述改变，并且去往6。
5、(退火程序)计算p＝exp[-(E′-E)/kT]。
如果p＞P，则接受所述改变，也就是说，虽然新的表面涉及相对于先前表面的能量增大，但是所述概率高到足以接受上坡步骤。P是0到1之间的随机数。
6、(退火程序)略微减小T。如果T＞Tmin则去往3，否则停止。
随着这一步骤逐渐降低温度，所述上坡移动变得逐渐较不可能，直到对于低的T值这种移动变为不可能为止，这意味着所述解决方案开始朝向最终分割迁入，只允许小的细化。
由于所述基于树的方法遵循经过所述表面上的点的系统路径，并且由于不存在唯一的树，因此所述表面分割的解决方案可能不是最优的。所述仿真退火方法没有预先定义的路径，使得能够有更大的机会进入全局最优，但是代价是需要更长的计算时间。
步骤3、无加权及加权PCA外形拟合
在前一步骤中，为每一个地标i赋予一个新的位置ji*。序列(j1*，...，jn*)代表一个外形，通过用该序列替换x并且通过求解$\mathrm{\Phi b}=x-\bar{x}$导出b可以将其拟合到形状模型$x\approx \bar{x}+\mathrm{\Phi b}$中。把x拟合到所述形状模型中意味着找到一个b，从而使得在高置信度地标中以小的误差近似x。
无加权PCA外形拟合
由于Φ∈□2n×t，因此等式多于未知数，从而所述系统是超定的并且没有确切解。因此，对于p的某些适当选择必须使得最小化。不同的范数给出不同的最优解。在p＝2的具体情况下，最小平方问题更易处理，并且在现有技术中通过基于正规方程和QR因数分解的方法得到解决。把实际点序列与拟合点序列之间的差异表示为拟合误差$e=x-(\bar{x}+\mathrm{\Phi b}),$该问题是寻求最小化eTe的b。替换e之后，该优化问题变为下式：
$\underset{b}{\min }J\left(b\right)=\underset{b}{\min }{e}^{T}e=\underset{b}{\min }({(x-\bar{x})}^T-b^T{\Phi }^T)((x-\bar{x})-\mathrm{b\Phi }).$
如果▽J(b)＝0，则J(b)被最小化。J(b)的梯度▽J(b)可以被写成下式：
$▿J\left(b\right)=-2{\Phi }^T{W}^2(x-\bar{x})+2{\Phi }^T{W}^2\mathrm{\Phi b}.$
因此，所述优化问题得到下面的等式：
$▿J\left(b\right)=0\iff {\Phi }^T\mathrm{\Phi b}={\Phi }^T(x-\bar{x}),$
其被称作正规方程。对于b求解这些方程得到下式：
$b={\left({\Phi }^T\Phi \right)}^{-1}{\Phi }^T(x-\bar{x})$
该等式把图像空间的数据矢量x投影到所述PCA模型空间中，并且返回近似系数b。可以通过估计所述形状模型$x=\bar{x}+\mathrm{\Phi b}$而把所述近似结果b投影回图像空间中。
加权PCA外形拟合
在把外形拟合到所述形状模型中时，可以把相应的票数用作权重。作为结果，所述PCA模型往往对具有高置信度的点给出高优先级。具有低置信率的点对于所述模型所生成的外形的影响将会变小。所述优化问题现在包括最小化加权二次误差(We)T(We)。通过像所述无加权情况那样使用利用正规方程的方法，该问题得到如下表示的针对b的解：
$b={\left({\Phi }^T{W}^2\Phi \right)}^{-1}{\Phi }^T{W}^2(x-\bar{x}).$
对该等式进行估计以便把x投影到所述PCA空间中，并且可以利用所述形状模型$x=\bar{x}+\mathrm{\Phi b}$投影回去。
迭代步骤1到3，直到满足停止标准为止，例如当所述地标的位置在相继迭代之间不发生显著改变。
所述无加权和加权PCA拟合都产生作为所述图像中的一个地标点集合x的最终分割结果。可以通过这些点内插连续的样条曲线以便解析地表示所述分割。
建模系统II(图3)
灰度外观模型
所述灰度外观模型捕获地标邻域内的空间强度分布。在本实施例中，所述空间模型从在地标的邻域内采样的圆形轮廓出发。被用来精确地描述肺野轮廓线的地标数目n可以低至14个，例如地标1表示每一个肺野的最高点，地标7位于左心影或右心影与其对应的半膈的横截面处，地标9表示肋膈角。其他地标等距离地定位在这些主要地标之间。很明显，在根据本发明的计算方案中也可以使用数目更多的地标。
-圆形轮廓
用于选择地标的邻域内的点的一种备选方式是采样以该地标为中心并且具有至少一个半径的圆上的点。在图9中给出了圆形采样的一个例子。如果所述地标位于(x0，y0)处，则在半径rc处在如下点对所述图像的灰度值函数进行采样：
$(x_0+r_c\cos \frac{2\pi }{n_c}(i-1),y_0+r_c\sin \frac{2\pi }{n_c}(i-1))i=1,...,n_c.$
所述样本被放入长度为nc的轮廓矢量中。对于nc的适当选择有12、8、6、4和3(对应于30、45、60、90和120度角间距)。可以同时考虑在一组半径下的多个圆形子采样。被无量纲地表示为图像尺寸的分数的适当的半径rc有0.05和0.25。
这种方案相比于所述线性轮廓的优点在于，其捕获邻域周围的2D结构。与线性轮廓相比，圆形轮廓更好地建模诸如肋膈角点或者左心影与左膈的交点之类的特定解剖地标，这是因为它们具有肺野外形的不连续的一阶导数。如果线性轮廓的长度过长，则其也可能会不明确地表示这些特定解剖地标附近的沿着所述肺野外形的地标，因为此时所述轮廓可能会穿越并置的肺野边界。
所述圆形轮廓可以被归一化，从而使得所述轮廓内的元素的绝对值之和为1。子采样的程度范围可以是从没有子采样(从而在图像的原始分辨率下产生重采样)到程度很高的子采样，例如在圆形轮廓的情况下仅主轴上的4个点。
在三维中，采用一种用来对体积内的地标的邻域进行采样的新颖方案。一种适当的选择是在以坐标为(x0，y0，z0)的地标为中心并且具有半径rc的球面上进行采样。必须实现对所述球面上的(x，y，z)坐标的各向同性的3D采样，即所述球面上的点密度必须是各处相等的。在采用球极坐标系时，其中表示在xy平面内相对于x轴的角度，θ表示相对于z轴的角度，rc是相对于以所述地标为中心的原点的半径(图15)，则到笛卡儿坐标的映射如下：


z＝rc cosθ与所述2D情况一样，可以按照均匀的方式在[0..2π]之间直接对(经度)进行采样，即其中i＝1，...，n1。但是对于θ(纬度)的等间距，在所述z轴上的投影在靠近极点处的更密集，因此只有在[-1..1]之间对cosθ进行均匀采样而不是在之间对θ进行均匀采样的情况下，才能使得z上的采样是均匀的。因此，$\theta =\arccos \left(\frac{j}{n_2}\right)$(其中j＝-n2，...，n2)是对于θ的采样。$t_1=\frac{2\pi }{n_1}$和$t_2=\frac{1}{n_2}$是分别对应于经度和纬度的周期。
可以同时考虑多个同心球面来对所述邻域进行采样。球面采样的优点在于，所有样本与所述地标位置的欧氏距离(定义为d2((x1，y1，z1)，(x2，y2，z2))＝((x2-x1)2+(y2-y1)2+(z2-z1)2)1/2)都是相等的。或者可以使用计算代价没有那么高的城市街区(city-block)距离(定义为d1((x1，y1，z1)，(x2，y2，z2))＝|x2-x1|+|y2-y1|+|z2-z1|)，从而在以所述地标为中心并且顶点在轴上的钻石形状上产生采样；或者还可以使用L∞范数或者棋盘距离(定义为d∞((x1，y1，z1)，(x2，y2，z2))＝max(|x2-x1|，|y2-y1|，|z2-z1|))，从而产生侧面平行于轴的作为等距离表面的立方体。后一种方案的优点在于，当把与所述中心的距离选择为体素尺寸的倍数时，其可以避免在各向同性的体积内进行内插。
-灰度特征图像
如在发明内容中所概述的那样，在N个数学特征中捕获围绕地标的图像结构。通过计算每一个位置x0周围的宽度为α的窗口内的导数图像的局部直方图的一定数目的矩而获得所述特征图像。在一个优选实施例中，计算一阶和二阶直方图矩、二次及以下的所有导数(L，Lx，Ly，Lxx，Lyy，Lxy)以及5个内尺度(σ＝0.5、1、2、4、8个像素)，从而得到一共N＝2×6×5＝60个特征图像。用来计算所述直方图的窗口尺度是根据所述内尺度，即α＝2σ。所述特征图像的集合可以被视为地标的邻域内的灰度函数的灰度分解。
-灰度成本
在特定地标i的新特征图像j当中采样的新轮廓gi，j之间的Mahalanobis距离被如下计算：
$h_{i,j}\left(p\right)={(g_{i,j}\left(p\right)-{\bar{g}}_{i,j})}^T{S}_{i,j}^{-1}(g_{i,j}\left(p\right)-{\bar{g}}_{i,j})$
更小的Mahalanobis距离意味着所述轮廓gi，j(p)源自具有均值和协方差Si，j的高斯分布的概率更大。因此，所述Mahalanobis距离可以被用作灰度成本度量，其被表示为hi，j(p)。该成本度量是一个给定图像中的位置p的函数。根据特征j，最有可能是地标i的真实位置的位置是使得hi，j(p)最小化的点p。因此可以为每一个特征j定义该灰度成本。
通过将所有特征j的所有灰度成本组合到下述和中，获得总灰度成本：
$h_i\left(p\right)={\underset{j=1}{\overset{N}{\Sigma }}(g_{i,j}\left(p\right)-{\bar{g}}_{i,j})}^T{S}_{i,j}^{-1}(g_{i,j}\left(p\right)-{\bar{g}}_{i,j})$
该和反映了p处的灰度图样与地标i处的预期灰度图样之间的相似度。根据灰度外观，最有可能与地标i重合的位置就是下述优化问题的解：
$p_o=\arg \underset{p}{\min }{h}_i\left(p\right).$
为了构造所述灰度外观模型，必须从训练图像估计特征轮廓gi，j(对于所有地标i＝1，...，n，对于所有特征图像j＝1，...，N)的分布
在三维中，同样地把特定3D点p处的成本函数估计为在特定3D地标i的新3D特征图像j中的该3D点p处球形采样的新3D体素水平轮廓gi，j之间的Mahalanobis距离hi，j(p)。同样地通过在一个总和hi(p)中组合所有特征j的所有灰度成本而获得一个总灰度成本，其反映了在一个球形采样邻域内在p处的灰度图样与地标i处的预期灰度图样之间的相似度。为了搜索所述3D地标i的最有可能的位置，必须解决优化问题。
为了在3D情况下构造所述灰度外观模型，必须从训练体积图像估计球形特征图像采样gi，j(对于所有地标i＝1，...，n并且对于所有特征图像j＝1，...，N)的分布
形状模型
灰度成本验证灰度图样或特征图样，而形状成本则被定义成验证两个相继地标之间的连接。
曲线轮廓线，每一个肺野一条，由n个点(地标)描述。这些地标在一组s个训练图像中被人工确定，从而产生一个点序列(x1，y1)...(xn，yn)＝(p1，...，pn)。这些坐标元组随后被放置在矢量x＝(x1，y1，...，xn，yn)T中，表示所述形状。考虑相继地标的一对所估计的位置(pi，pi+1)。为矢量vi＝pi+1-pi分配一个形状成本，其反映了vi关于其概率分布的似真性。从训练形状估计v1，v2，...，vn的分布，假设v1，v2，...，vn具有围绕其均值的正态分布。vi的均值和协方差分别被表示为和。在图10中给出了这些矢量分布的一个例子。所述矢量分布可以被视为由所述地标集合表示的形状的一阶(切线)几何分解。
验证两个相继地标pi、pi+1之间的连接矢量vi(完全由其取向及其长度表征的平面内的矢量)的一种新颖的形状成本是vi与其均值之间的Mahalanobis距离：
$f_i\left(v_i\right)={(v_i-{\bar{v}}_i)}^T{S}_{v_i}^{-1}(v_i-{\bar{v}}_i).$
因为其Mahalanobis距离较大而不太可能发生的连接将获得高的形状成本。另一方面，将为与预期值相等的连接分配零成本。
在三维中，每一个地标i是通过边连接到与所述表面上或者所述表面内部的相邻地标节点相对应的一定数目的顶点vi，j(j＝1..Vi)的基于点的表示中的节点或顶点。每一条边是由其在3D空间中的取向和长度表征的一个3D矢量，并且每一条边的3D分布是从训练图像集合导出的。在实践中，这意味着覆盖所述表面的多边形的每一条边是从训练实例获知的。具体来说，当所述多边形是三角形时，对所有三角形的三条边的连接矢量的三元组进行建模。当所述体积内部的点也被包括并且在所述表面上和所述体积内部的全体点之间定义多面体单元时，从所述训练集合获知所述多面体分解的边的分布。具体来说，当所述多面体是四面体(其由四个点和六条连接边构成，从而形成四个三角形面)时，对所述六条边当中的每一条的空间分布进行建模。所述协方差矩阵有效地对所述连接矢量之间的空间相关性进行建模，这是因为其捕获可能存在于不同解剖地标的位置之间的相关变化。或者，利用所述三角剖分表面表示或者所述四面体体积表示，还可以对由三角形形成的平面贴片的法线的分布进行建模。
分割系统II(图4)
在训练阶段采集关于图像中的对象的灰度外观以及关于待分割对象的形状的知识。该知识随后被用来根据以下步骤在新的图像中分割所述对象。
步骤1、矩形搜索栅格构造
构造每一个地标i的矩形搜索栅格以便约束该地标的可能位置。每一个矩形栅格由几何范围和栅格间距这样的参数表征。
所述栅格间距rg应当是根据所述灰度外观模型的圆形轮廓的半径rc。一种适当的关系是rg＝fgrc，其中rc是半径，fg是一个固定分数。所述栅格间距的一个典型分数是fg＝0.5。
所述栅格范围和栅格位置是特定于地标的，并且由xmin、xmax、ymin和ymax表示。真实的(未知)地标位置(x*，y*)被认为位于下边界与上边界之间：
xmin＜x*＜xmax。
ymin＜y*＜ymax
只有在所述栅格跨越整个图像区域(即xmin＝ymin＝0并且xmax＝ymax＝1)的情况下才能确保这些条件。一种更为有效的方法是把所述搜索约束到所述图像的相关部分。假设所述地标的x坐标的概率分布px(x)是以为均值并且以σx为标准偏差(从所述训练形状估计)的正态分布。联合限定x方向上的栅格范围的栅格边界xmin和xmax，被选择成使得下式成立：
${\int }_{x_{\min }}^{x_{\max }}{p}_x\left(x\right)\mathrm{dx}=f_c$
其中fc是代表条件xmin＜x*＜xmax为真的概率的某一分数。该分数fc的一个适当的值是0.995。如果要求所述范围围绕所述均值对称，即$x_{\max }-\bar{x}=\bar{x}-x_{\min },$则所述上边界xmax是满足下式的值：
$\frac{1}{{\sigma }_x\sqrt{2\pi }}{\int }_{\bar{x}}^{x_{\max }}\exp (-\frac{{(x-\bar{x})}^2}{2{\sigma }_x^2})\mathrm{dx}=\frac{1}{2}{f}_c.$
于是所述下边界就是$x_{\min }=2\bar{x}-x_{\max }.$可以类似地获得y坐标的边界。
在三维中，每一个地标i的搜索栅格是与坐标轴对准的平行六面体。所述栅格间距rg应当是根据所述灰度外观模型的圆形轮廓的半径rc。一种适当的关系是rg＝fgrc，其中rc是所述半径，fg是一个固定分数。对于所述栅格间距，一个典型分数是fg＝0.5。
所述栅格范围和栅格位置是特定于地标的，并且由xmin、xmax、ymin、ymax、zmin和zmax表示。真实的(未知)地标位置(x*，y*，z*)被认为位于下边界与上边界之间：
xmin＜x*＜xmax
ymin＜y*＜ymax
zmin＜z*＜zmax。
只有在所述栅格跨越整个图像体积(即xmin＝ymin＝zmin＝0并且xmax＝ymax＝zmax＝1)的情况下才能确保这些条件。一种更为有效的方法是把所述搜索约束到所述图像的相关部分。与所述2D方法类似，所述栅格范围是从所述训练图像集合导出的，以便覆盖3D体积内的地标位置的积分概率分布函数的分数fc。该分数fc的一个适当的值是0.995。
与把地标的解限制到位于与当前分割轮廓线垂直的1D上的所述线性轮廓采样不同，所述搜索栅格把解空间扩大到位于围绕当前地标的3D邻域内。
步骤2、灰度成本矩阵构造
所述灰度外观模型被用来找到每一个地标的适当位置。根据所述灰度外观模型，如下选择地标i的m个最佳位置：
-首先，根据前一步骤定义覆盖每一个地标i的预期位置周围的图像区域的相关部分的矩形栅格。
-其次，对于每一个栅格点，计算总灰度成本hi(p)。选择出对应于m个最低总灰度成本的点pi，1，pi，2，...，pi，m。
-第三，把灰度成本矩阵C构造成使得每一行i包含地标i的m个最可能位置的成本：
$C=\left[\begin{array}{ccccc}{h}_1\left(p_{\mathrm{1,1}}\right)& ···& h_1\left(p_{1,k}\right)& ···& h_1\left(p_{1,m}\right)\\ ·& & ·& & ·\\ ·& & ·& & ·\\ ·& & ·& & ·\\ h_i\left(p_{i,1}\right)& ···& h_i\left(p_{i,k}\right)& ···& h_i\left(p_{i,m}\right)\\ ·& & ·& & ·\\ ·& & ·& & ·\\ ·& & ·& & ·\\ h_n\left(p_{n,1}\right)& ···& h_n\left(p_{n,k}\right)& ···& h_n\left(p_{n,m}\right)\end{array}\right].$
每个地标的最佳点的典型数目m的范围是从10个到40个。
由于每个地标的所述m个最佳点是独立于相邻地标的m个最佳点的集合选择的，因此可能会出现这些点当中的一个或更多个更靠近非相邻地标的情况。通过在下一步骤中考虑所述形状成本，将很有可能在最终的外形分割中忽略这些点。
在三维中，灰度成本矩阵的等价物是灰度成本图，其中每一个节点具有代表所述地标的m个最低总灰度成本点的数目的重数(multiplicity)。与所述2D情况一样，所述灰度外观模型被用来找到每一个地标的适当位置。根据所述灰度外观模型，如下选择地标i的m个最佳位置：
-首先，根据前一步骤定义覆盖每一个地标i的预期位置周围的图像区域的相关部分的3D栅格。
-其次，对于每一个3D栅格点，利用所述球形采样的邻域计算总灰度成本hi(p)。选择出对应于所述m个最低总灰度成本的点pi，1，pi，2，...，pi，m。
-第三，把灰度成本图G构造成包含地标i的m个最可能3D位置的成本，其中每一个节点i具有一个重数m。从所述表面三角剖分的节点开始并且遵循相继链接所述表面上的所有最近相邻节点的边，构造所述图。在每一个节点中存储前一步骤的m个最低总灰度成本点。每个3D地标的最佳点的典型数目m的范围是从10个到40个。
由于唯一的标准是灰度成本，因此这一步骤可能会导致一个节点在其m个最可能点的集合当中包含在3D空间中更靠近节点i的各非相邻节点的点，特别在所述表面上密集采样所述地标的情况下尤其如此。通过在下一步骤中考虑所述形状成本，有可能在最终的表面分割中忽略这些点，也就是说它们将代表不太可能的连接。
步骤3、最小成本路径构造
确定分割所述对象的外形被简化成通过每行选择一个元素找到一条从上到下穿过所述矩阵C的路径。把行i中的所选元素的索引表示为ki，则所述曲线变为点序列最优路径(k1*，k2*，...，kn*)是使得成本函数J(k1，...，kn)最小化的路径：
$(k_1^{*},k_2^{*},...,k_n^{*})=\arg \underset{k_1,...,k_n}{\min }J(k_1,...,k_n)$
上面引入的模型允许多种成本度量。考虑两个相继地标和根据所述灰度外观模型的成本分量和根据所述形状模型的成本分量如下：
-对应于所述地标和的灰度成本和
-形状成本其验证从到的连接的似真性。可以借助于总体灰度成本与总体形状成本的加权和把这两种成本类型都合并到所述成本函数J(k1，...，kn)中。具有权重因数γ1的总体灰度成本是地标个体灰度成本(i＝1，...，n)的和。具有权重因数γ2的总体形状成本是形状成本(i＝1，...，n)的和。因此，所述成本函数J(k1，...，kn)变成下式：
$J(k_1,...,k_n)={\gamma }_1\underset{i=1}{\overset{n}{\Sigma }}{h}_i\left(p_{i,k_i}\right)+{\gamma }_2\underset{i=1}{\overset{n-1}{\Sigma }}{f}_i(p_{i+1,k_{i+1}}-p_{i,k_i})+{\gamma }_2{f}_n(p_{1,k_1}-p_{n,k_n}).$
使得J(k1，...，kn)最小化的最优路径(k1*，k2*，...，kn*)被称作最小成本路径。利用现有技术中已知的动态编程技术计算所述最优路径，比如在G.Behiels等人的“Evaluation of image features and search strategies for segmentationofbone structures using active shape models(利用主动形状模型估计用于分割骨骼结构的图像特征和搜索策略)”(Medical Image Analysis，6(1)：47-62，2002年)中介绍的动态编程技术。所述灰度成本的典型权重因数是γ1＝1，所述形状成本的典型权重因数是γ2＝0.25。
图13示出了根据前面概述的系统对三个胸廓图像的自动分割，并且将其与由经验丰富的放射科医师执行的人工分割进行比较。在每一个肺野上示出了14个自动放置的地标以及经过这些地标内插的连续样条。
在三维中，用于确定最优分割表面的最小成本表面构造可以同样使用灰度成本与形状成本的组合。与利用线性采样的先前模型中一样，经过所述表面图的点的路径可以具有对最小跨度树(MST)或深度优先树的遍历的形式。
确定分割所述对象的表面被简化成对于所述表面上的每一个地标节点找到所述m个点的其中之一。把节点i中的所选元素的索引表示为ki，则所述表面变为所遍历的树的点序列最优表面(k1*，k2*，...，Kn*)是使得成本函数J(k1，...，kn)最小化的点集合：
$(k_1^{*},k_2^{*},...,k_n^{*})=\arg \underset{k_1,...,k_n}{\min }J(k_1,...,k_n)$
上面引入的模型允许多种成本度量。考虑通过一条弧线相连的两个相邻地标和根据所述灰度外观模型的成本分量和根据所述形状模型的成本分量如下：
-对应于通过一条弧线相连的所述表面上(或者附加地，所述表面内部)的相邻地标和的灰度成本和
-形状成本其验证3D空间中的从到的连接矢量的似真性，其中所述矢量完全由其取向和长度表征。
可以借助于总体灰度成本与总体形状成本的加权和把这两种成本类型合并到成本函数J(k1，...，kn)中。具有权重因数γ1的总体灰度成本是地标个体灰度成本(i＝1，...，n)的和。具有权重因数γ2的总体形状成本是形状成本(i＝1，...，n)的和。因此，所述成本函数J(k1，...，kn)变成下式：
$J(k_1,...,k_n)={\gamma }_1\underset{i=1}{\overset{n}{\Sigma }}{h}_i\left(p_{i,k_i}\right)+{\gamma }_2\underset{i=1}{\overset{n-1}{\Sigma }}{f}_i(p_{i+1,k_{i+1}}-p_{i,k_i})+{\gamma }_2{f}_n(p_{1,k_1}-p_{n,k_n}).$
所述序列p1，...，pn代表随着所述树被遍历而被访问的节点的链。使得J(k1，...，kn)最小化的最优表面(k1*，k2*，...，kn*)被称作最小成本表面。
计算加速细化
包括步骤1-3的所述分割方法可以用来分割一个或更多解剖实体的轮廓线的仅仅一部分。为此目的，仅仅对于连续地标的子集构造所述矩形搜索栅格。所述灰度成本矩阵将包括等于所保留的地标数目的行数。所述最小成本路径构造将使得仅仅包括所保留的地标的灰度项和形状项的成本函数最小化。所得到的分割仅仅追踪由所保留的地标覆盖的轮廓线部分。当人们仅仅对解剖实体轮廓线的一部分感兴趣时或者当只需要定位某些特定地标以获得特定测量时，部分分割的明显优点在于计算速度。
用于加速分割的计算的另一种算法细化是使用从粗到细的多分辨率方法。在每一个地标的多分辨率表示上构造所述灰度外观模型和形状模型，引入了空间尺度作为新的维度。为了实现这一目的而实施两种备选方案，这两种备选方案都在相继更细的分辨率水平上应用所述分割概念。
在一种备选方案中，在分割期间，在原始图像中，减少轮廓中的点数、所述搜索栅格中的点数以及沿着所述外形的地标点数。随后通过相继地增大在前一个解的附近的搜索栅格的分辨率(因此在每次迭代中可以把所述轮廓和搜索栅格中所考虑的总点数保持恒定)并且应用分辨率更精细的灰度和形状模型，来提高地标的位置精度。
在另一种备选方案中，可以利用一定数目的水平的高斯和拉普拉斯金字塔及其相关联的导数特征图像而实现所述多分辨率方法。地标的初始位置在粗糙的水平上确定，并且利用前一水平的中间位置作为下一水平的起始位置在最精细的水平上被跟踪到其最终位置。由于粗糙分辨率图像包含较少细节，因此所述搜索空间很可能将包含较少的假最小值，从而使得能够在给定水平下进行更快的优化，并且允许更快地朝向其最终位置定位地标。
训练所述分割模型
从图像存储库收集一定数目的胸廓图像并且在图形用户界面上顺序地呈现。执行以下步骤以便建立所述分割模型：
-经验丰富的用户通过(利用鼠标左键点击)沿着描绘所述肺野的外形指示出一定数目的点，而在每一个所显示的图像中人工分割所述肺野。这些点不需要是沿着所述肺野轮廓线等间距的；唯一的要求是所述点集合共同以足够高的程度近似所述肺野。为了评估所述解剖拟合，经过到目前为止确定的人工定位的点连续内插样条曲线，直到通过鼠标右键点击从最后一点到第一点闭合所述曲线。可以通过把单独一点拖曳到一个新位置来进行人工分割调整。再次关于解剖正确性评估所得到的轮廓线。
-接下来，将把一定数目的地标自动放置在所述人工分割的肺野外形上。为了使得所述训练集合当中的所有图像上的相同的地标映射在彼此之上，用户可以定位少数几个容易辨别的地标，并且通过把点等距地放置在所述肺野外形上而获得其他地标。在肺野轮廓线的情况下，一定数目的容易辨别的地标是最高肺野点、表示肋膈角的点以及心影与半膈的结点，小计是三个点。接下来选择地标总数，一种适当的选择范围例如是从14个到40个。在14个地标的情况下，点p1、p7和p9代表三个固定的地标，五个点被均匀地分布在p1与p7之间，一个点分布在p7与p9之间，另外五个点分布在p9与p1之间。在左、右肺野上分别执行这一步骤。
-接下来，询问关于训练阶段的参数：(a)用于训练(和分割)的图像尺寸，一个典型值是256或512；(b)将由所述主分量分析(PCA)解释的形状方差的分数，一个典型值是0.9999；(c)线性或圆形轮廓中的点数，圆形轮廓的典型值是12、8、6、4或3；(d)作为所述图像维度的分数的圆形轮廓的半径，典型值是0.04、0.03、0.02。
-对所述灰度外观模型进行训练：(a)分解地标周围的灰度函数，即，计算N个特征图像，例如前面概述的灰度值导数的局部直方图矩；以及(b)收集地标pi处的线性或圆形轮廓的灰度值分布，即计算(i＝1...n，j＝1...N)和协方差矩阵Si，j(i＝1...n，j＝1...N)(即对于特征图像j内的每一个地标i)。这一步骤生成统计灰度外观知识。
-对所述形状模型进行训练：(a)计算包含在所述地标中的零阶(位置)信息的几何分解，即，计算均值形状和逐列地设置在矩阵b中的t个主变化模式(特征形状)；(b)计算包含在地标的连接序列中的矢量分布(一阶切线信息)，即，计算vi的均值和协方差这一步骤生成统计形状知识。
-存储所述统计灰度外观和形状知识，以便例如用于根据上面给出的基于模型的分割子系统分割新图像中的肺野。
把2D图像中的基于模型的分割应用于心胸比(CTR)的计算
自动肺野分割的一种应用是心胸比(CTR)的计算。所述CTR(图6)被定义为心脏的横直径与胸廓的内直径(ID)的比值：
$\mathrm{CTR}=\frac{\mathrm{MLD}+\mathrm{MRD}}{\mathrm{ID}}$
其中，MLD是心脏左侧的最大横直径，MRD是心脏右侧的最大横直径。这一指数是一个重要的临床参数，其对于成年人在39％到50％之间变化，并且平均值是大约45％。高于50％的心胸指数被视为异常。可能的原因有心力衰竭、心包积液以及左心室肥厚或右心室肥厚。有可能利用本发明中所公开的自动肺野分割自动计算所述心胸比。
参照图8，通过选择右肺分割的地标p1与p9之间的拟合外形段上的最左笛卡尔点而获得定义MRD的特性点；通过选择左肺分割的地标p1与p7之间的拟合外形段上的最右笛卡尔点而获得定义MLD的特性点。通过把这些特性点的列坐标相减并且取绝对值而获得MLD+MRD这一和。类似地，通过以下步骤获得所述ID：选择右肺分割的地标p1与p7之间的拟合外形段上的最左笛卡尔点和左肺分割的地标p1与p9之间的拟合外形段上的最右笛卡尔点、把这些特性点的列坐标相减并且取绝对值。
图14示出了对已分割肺野上的特性点的计算，(a)利用所述人工肺野分割的心胸比计算得到CTRman＝0.47；(b)利用根据本发明的自动导出的肺野分割的心胸比计算得到CTRautomatic＝0.45。
把所述基于模型的分割和测量系统应用于其他身体部分
所述肺野分割中的每一个地标的搜索栅格的空间范围，是在训练胸廓图像集合中的所有类似地标的位置的基础上导出的。使用搜索栅格来约束给定地标的候选位置的概念可以被扩展用于其他身体部分，所述其他身体部分在图像中可能具有比胸廓图像中的肺野更宽的位置、旋转和尺寸变化以及确实占据整个图像区域(这与胸廓图像中的肺野不同)。
为了允许所述身体部分的这种位置、旋转和尺寸变化，可以与本发明相结合地应用搜索栅格的锚点映射的概念，其中利用了在标题为“Method for automatically mapping of geometric obj ects in digital medicalimages(用于自动映射数字医学图像中的几何对象的方法)”的欧洲专利申请04076454中公开的方法。所述经过细化的方法于是变得对于分割整形外科射线照片中的骨骼身体部分特别有意义，这是因为所述骨骼身体部分在图像中由相对恒定的形状表征并且它们可以被锚定到所述图像中的明确指明的骨骼地标。
举例来说，在骨盆、臀部或全腿检查中，股骨轮廓线可以被锚定到公知的地标，比如大转子的尖端、膝关节中心以及股骨头中心。这些锚点被用来建立模型锚点与在实际图像中选择的相关联的锚点之间的仿射变换。由于搜索栅格包括设置在模型图像中的分割模型外形上的地标点周围的矩形格子内的候选点的集合，因此通过对所述模型点的坐标应用所述变换把每一个构成候选点依次映射到所述实际图像中。按照这种方式，在给定地标点的最有可能的位置周围重建该点的搜索栅格。随后，通过优化经过所选择的候选位置组合(每一个搜索栅格上一个位置)的路径的组合灰度值与形状成本，按照上面所公开的方式继续所述优化处理。具有最优成本的路径就是对所述身体部分的最终分割。其他骨骼的实例有手部和上肢骨骼、足部和下肢骨骼、骨盆以及脊椎骨和结构。
适用于这种类型的基于地标的分割的其他身体部分还有MR图像的3D CT中的软组织器官。一方面，在这里可以采取一种逐切片的方法，其确定每一个切片上的解剖结构的已分割外形，并且把所有切片的结果组合成3D分割表面。另一方面，可以采用一种在3D体积上扩展当前概念的全3D方法，其对限制在3D搜索栅格内的3D地标建立并且应用模型。肾脏、肺、心脏、肝脏、胃、脾脏、前列腺以及大脑结构是本发明的方法所适用的器官的实例。
在所提到的所有应用实例中，测量点可以基于所述分割上的地标点的位置，或者是基于几个地标的组合(比如一段长形骨骼的皮质的每一侧上的两个并置点的中点，一对这样的中点确定所述骨骼的解剖轴)，随后可以把所述地标点馈送到比如在EP A 1349098中公开的测量系统中。
把3D图像中的基于模型的分割应用于肺总容积(TLV)的计算
用于分割2D医学图像中的对象的方法可以如前所述地被扩展用于分割嵌入在诸如CT或MR的3D体积图像数据中的对象。所得到的分割可以用于导出定量测量以便帮助诊断。这些测量包括诸如器官容积或者病理组织体积之类的尺寸测量。定量测量可以被用来把单个患者与其先前的测量进行比较以便跟踪疾病发展或治疗进展，或者把一个患者与其他患者或参考值进行比较。经常从体积图像中导出的测量量是对象或器官的容积V，其由已分割表面S所包围的体素表示，并且被如下计算：
$V=\underset{v\in V}{\Sigma }\mathrm{\Delta x\Delta y\Delta z},$
其中，Δx和Δy代表切片内体素维度，Δz代表切片之间的间距。
举例来说，可以在3D肺野分割的基础上计算肺总容积(TLV)，并且将其与先前的TLV进行比较以便评估改变。与根据年龄和性别的TLV参考值的比较可以被用来评估所述肺容积是否表明肺气肿。显而易见，在CT肺野分割的基础上计算的TLV可以被分成左、右肺容积，这种区分无法通过呼吸肺功能测试获得。
根据本发明的用来在医学图像中构造解剖对象的模型并且分割这些解剖对象的方法可以被应用于其他器官。所述器官的实例有2D射线照片上的臀部或膝关节轮廓线，以及CT或MR图像上的骨骼、大脑、胸廓和腹部结构。后续的测量可以基于包含在所述(多个)对象的已分割表面内的几何信息，或者可以更直接地基于位于所得到的已分割表面上或其内部的定义地标的位置。