CN1924930B

CN1924930B - 分割数字医学图像中解剖实体的方法

Info

Publication number: CN1924930B
Application number: CN2006101266313A
Authority: CN
Inventors: P·德沃勒
Original assignee: Agfa HealthCare NV
Current assignee: Agfa NV
Priority date: 2005-08-30
Filing date: 2006-08-30
Publication date: 2011-10-19
Anticipated expiration: 2026-08-30
Also published as: CN1924929A; CN1924929B; EP1760658A1; EP1760658B1; CN1924930A; DE602005012261D1

Abstract

为图像内多个标记的每一个标记定义该标记的初始位置。接着，采样初始位置周围的邻域，该邻域包括该标记的多个候选位置，并将候选位置的每一个位置与成本相关联。优化成本函数，该成本函数表达了所有候选位置的总体灰度级成本和总体形状成本的加权总和。将分割解剖实体定义为穿过候选位置的选定组合的路径，该组合的成本函数被优化。

Description

分割数字医学图像中解剖实体的方法

技术领域

本发明涉及分割数字图像中实体的方法，具体地分割数字医学图像中的解剖学实体的方法.根据本发明的分割处理是基于利用图像的几何和光度模型。

本发明的方法还可以用于对数字图像，具体地对数字医学图像进行几何测量处理.

背景技术

在放射实践中，频繁地使用几何测量辅助对异常的诊断。为了进行这些测量，必需将关键用户点放置在图像中例如显示装置上显示的图像中的其相应解剖标记位置上。诸如两点之间距离或者直线之间夹角的测量是基于这些关键用户点的位置.另外，可以通过评估整体几何以确定正常或异常，这涉及对完整形状进行分析.因此，需要自动和客观地提取内嵌在放射图像内的定量信息.

这种频繁进行测量的示例为计算胸腔RX图像中的心胸比率(CTR)(图6)。这个比率定义为心脏在顶点水平处的横向直径与胸腔的内径(ID)的比例，即CTR＝(MLD+MRD)/ID。心脏的横向直径由心脏左侧的最大横向直径(MLD)和心脏右侧的最大横向直径(MRD)组成.显然，这种定义要求放射科医生沿左和右胸廓(ribcage boundary)的内图廓寻找，以定位计算内径ID所需的点对(point pair)。该点对必需位于胸腔的最大内径处。同样，必需找到左和右心脏阴影边缘，以定位计算MLD和MRD之和所需的点.更具体地，这些点位于相对于脊柱中线的最远端。寻找边缘的过程需要放射科医生进行解剖分割并在被分割的解剖上定位这些点(这个示例中总共有四个点)。该分割步骤，对于CTR计算的情形，相当于描绘肺野(lungfield)。

数字图像中的许多其它测量遵从相似的方法，涉及分割解剖器官或实体，其中在该解剖器官或实体上分割几何特征点和测量对象被确定并最终进行测量(图5)。

参考心胸比率计算的示例，为了自动定位所需的点，需要在胸部放射图像上自动地分割肺野的方法。

根据应用，可以通过许多种方式解决分割问题。分割策略从早些年的计算机视觉的低水平策略进展到最近的基于模型的策略.

低水平方法依靠局部图像算子，将具有不同光度特征的像素分离并将具有相似局部光度特征的像素分组。这两种的示例为边缘检测和区域生长。尽管这些低水平方法性能差，但这些方法在绝大多数商业图像分析工具中非常流行。主要原因为这些低水平方法易于理解和实施。然而对于复杂图像数据，例如医学图像中存在的以及如前所述的胸腔图像的内容为代表的复杂图像数据，这些方法的有效性是有限的。

更成功的方法包括有关待分割的形状以及有关图像中对象的光度或灰度级外观的先验知识.这些方法(称为基于模型的方法)经常是基于模板匹配。例如通过相关或使用广义霍夫变换技术来匹配模板.不幸的是，对于医学图像的情形，模板匹配可能不成功.这是因为解剖对象呈现的形状和灰度级外观具有很大的可变性。

基于由Kass等(M.Kass，A.Witkin，和D.Terzopoulos，Snakes：active contour models，Int.J.Computer Vision，1(4)：321-331，1988)提出的活性轮廓(contour)以及基于水平集合(J.A.Sethian，Levelset methods and fast matching methods，Cambridge Univ.Press，Cambridge，U.K.1999)的方法能应付更大的形状可变性，但仍然不适用于许多医学分割任务，因为这些方法包括很少关于待分割对象的先验知识.

手工参数模型克服了这个问题，但这些模型仅限于单个用途。

鉴于这些缺点，显然，需要一种通用分割方案，可以使用示例进行训练从而采集有关待分割对象的形状以及图像中该对象的灰度级外观的知识。由Cootes和Taylor(T.F.Cootes，C.J.Taylor，D.Cooper，J.Graham，Active Shape Models-theirt raining and applications，Computer Vision and Image Understanding，61(1)：38-59，1995)提出的主动形状模型(ASM)满足这种分割方案定义。由标记点向量的主分量给出该形状模型.灰度级外观模型描述了中心位于与对象轮廓垂直的每个标记处的轮廓的归一化一阶导数的统计。通过最小化该一阶导数轮廓与该轮廓分布之间的马氏(Mahalanobis)距离，可以在新图像中找到标记的位置。这种算法开始于初始化评估，并执行拟合程序，该拟合程序是标记位移和形状模型拟合的更迭。已经设计出许多相似的方法，所有这些方法都采用了三步程序.首先，这些方法都使用形状模型，保证产生看似可信的结果.其次，这些方法使用灰度级外观模型，以将对象放置在边界周围或对象内部的灰度级图形和从训练示例预计的灰度级图形相似的位置上。最终，该算法通过最小化一些成本函数而拟合该模型。

由主动形状模型提出的方法仍然面临诸多限制。

第一限制为需要初始评估，称之为初始定位问题。在一些情况下，需要在这个框架内进行广泛的搜寻以获得合适的初始轮廓。

第二限制为形状模型和灰度级外观模型的交替使用。首先，使用灰度级外观模型更新该分割评估。其次，将该形状模型拟合到被更新的轮廓。不幸的是，如果灰度级外观模型错误地定位标记，则该形状模型会受到误导.

使用主动形状模型分割方案的一些实验呈现另一个问题.如果特定标记的灰度级外观模型被用于寻找该标记的更佳位置，要求真实标记位置位于该算法所探索的区域内。在一些情况下，该算法会在另一个标记的区域内搜寻标记.这意味着在错误区域使用了错误的灰度级外观模型，导致不良分割。

本发明的目标是提供一种分割数字图像中实体的方法，具体地分割数字医学图像中的解剖学实体的方法，该方法解决了现有技术的问题。

发明内容

通过如权利要求1的方法实现上述方面。

在从属权利要求中陈述了本发明的优选实施例的具体特征.

本发明的另一个方面涉及如权利要求书中陈述的测量数字医学图像中解剖实体的方法。

在下文中，术语灰度值模型、灰度值外观模型、灰度级模型、灰度级外观模型以及光度模型被当作同义词使用。

同样，术语形状模型和几何模型被当作同义词使用。

术语特征图像和特征也被当作同义词使用.

本发明方法的实施例通常被实施成计算机程序产品的形式，当在计算机上运行时用于执行本发明的方法步骤。

该计算机程序产品通常被储存于计算机可读取载体介质例如DVD或CD-ROM中。备选地，该计算机程序产品采取电信号的形式，并可以通过电子通信与用户通信。

在EP A 1349098中公开了由双向链接外部图像模型和内部信息模型组成的，通过将测量对象和实体分组成计算机化的测量方案，而实现数字采集医学图像中测量自动化的方法。

在根据EP A1349098的测量对话期间，从计算机提取测量方案并激活该测量方案。在被激活的测量方案的指导下，随后对显示的图像进行测量。

在这种计算机化方法中，大量几何对象被映射到数字图像中，其它测量对象和最终测量实体(例如距离和角度)都是基于这些几何对象。基础的几何对象典型地为关键用户点，这些关键用户点定义了几何测量所基于的其它几何对象。

然而上述专利申请没有公开如何可以自动地实现该映射而无需用户定位和操纵关键测量点.根据本发明的方法可以用于实现该映射.

根据本发明，提供了一种分割方法，包括构造诸如医学图像中器官或结构的解剖实体的模型，可以通过示例训练这些模型，从而获得有关灰度值外观的知识以及有关待分割解剖实体形状的知识.

该灰度值外观模型包括灰度值的概率分布以及灰度值在分布于解剖剖面上一组标记处的多个多尺度导数。

该形状模型包括沿解剖实体轮廓的连续标记之间的每个连接向量的概率分布。

提出基于离散点的对象表示，用于描述医学图像中的解剖轮廓。

公开了用于产生模型的两种策略以及相关系统以及将这些模型应用于分割真实图像。这些策略，尽管分解成具有相同目标的多个结构单元，但构造和采用了不同的光度和几何模型知识。

首先，通过训练图像而构造新的灰度级外观模型，编码在标记位置的光度知识。这个步骤利用了每个标记周围被采样邻城内的强度关联.

其次，构造形状模型，编码标记之间的几何知识。这个步骤利用被分割对象的标记之间的空间关联。

在进一步的分割步骤中，光度和几何知识被共同地用于分割新图像上的一个或多个解剖轮廓。

结果的分割可以在对图像执行几何测量的过程中被用于推导预期测量点的位置。

可以以直接或者间接的方式定义测量点.

在直接方式中，标记的位置被用作独立的测量点.

在间接方式中，通过分割过程得到的标记通过曲线被互连，且该曲线上的特征点被定义成测量点.

灰度级外观模型

每个模型化系统中的第一步骤将标记的几何位置限制成位于线性路径上，该线性路径垂直于位于标记中心上的现有轮廓或稀疏栅格。图像中标记点处的光度图形被轮廓向量捕获，其中在基于多个尺度下提取的图像导数的特征图像中采样该轮廓向量。该轮廓采样可以沿图像内的线性路径，或者围绕标记点的圆形路径。

特征图像

当计算图像内所有图像点的光度图形时，获得特征图像，该特征图像的大小等于原始图像。相邻运算通常被用于计算这些特征图像.特定点的特征图像内特征的数值可被排列成特征向量。特征图像的计算通常包括两个步骤：(a)使用一组卷积核进行线性过滤的步骤，和(b)涉及结合计算局部图像统计的非线性点运算的多个子步骤的后期处理步骤。在步骤(b)中，可以单独地执行这些子步骤之一。

不同滤波器和不同类型的后期处理导致不同类型的特征图像。例如，Laws构造25个二维卷积核，作为重组理想特征轮廓例如Leve1、Edge、Spot、Wave和Ripple的向量对的外积。后期处理步骤包括非线性开窗口(windowing)，其中将绝对过滤数值在更大窗口内求平均以获得图像纹理能量测量。Unser使用的滤波器组构造成对应于公知变换的向量的外积，该公知变换是例如离散正弦(DST)、离散余弦(DCT)、离散偶正弦(DEST)、离散实偶傅立叶(DREFT)、离散实奇傅立叶(DROFT)变换，并用于后期处理通道方差的计算。最终，使用被调谐成不同空间频率和不同取向的组合的对称或反对称Gabor滤波器构造Gabor滤波器组，使得所有输出具有相等的频率带宽。在这种情况下，该非线性点运算步骤可包括设定阈值.

使用基于多分辨率小波滤波器、Rank值过滤器或自适应滤波器的特征图像，可以测量其它图像特性。

特征图像基于局部无序图像(LOI)。医学图像中解剖标记周围的图像结构呈现典型的灰度级图形。在本发明中，在如后文中勾画出的N个数学特征中捕获这个结构.使用特征图像建造灰度级外观模型的概念是稍早由B.Van Ginneken等，Active shape model segmentation withoptimal features，IEEETrans.On Medical Imaging，21(8)：924-933，2002提出的，该模型是基于由J.J.Koenderink和A.J.Vandoorn，Thestructure of locally orderless images，Int.J.of ComputerVision，31(2)：159-168，1999提出的所谓局部无序图像(LOI)。

可以利用泰勒展开式，使用K阶多项式近似在感兴趣点即标记附近的位置x₀的灰度值函数f。各项的系数由在x₀的导数f⁽ⁱ⁾给出：

f (x) \approx Σ_{i = 0}^{K} \frac{1}{i!} f^{(n)} (x_{0}) {(x - x_{0})}^{i}

这个函数近似中的所有信息都包括在导数f⁽ⁱ⁾内.类似地，可由包括原始图像的导数的图像滤波器组(L，L_x，L_y，L_xx，L_yy，L_xy，...)描述该图像.可以使用许多模糊或内尺度σ进一步计算每个导数图像(derivative image)。通过计算后期处理步骤，最后获得这些特征图像，作为每个导数图像的局部图像统计。该局部图像统计包括每个位置x₀附近宽度为α的窗口内的许多时刻(moment)的局部直方图.在具体实施例中，计算第一和第二直方图时刻、一直到二阶的所有导数(L，L_x，L_y，L_xx，L_yy，L_xy)以及5个内尺度(σ＝0.5、1、2、4、8像素)，形成总共N＝2×6×5＝60个特征图像。计算直方图的窗口尺度是根据内尺度，即，α＝2σ.

在特征向量中直接在行和列中使用这些导数(即，x和y)的缺点为平移和旋转缺乏不变性(即欧几里得变换)。除非这些图像抗旋转成标准取向，否则这些算子只能在具有相同取向的示例上进行训练.笛卡儿微分不变式(CDI)描述了与所选择的笛卡儿坐标系无关的图像微分结构，因此可以省略抗旋转步骤。由高斯微分算子构造CDI在L.M.J.Florack等，Scale and the differential structure of images，Image and Vision Computing，Vol.10(6)：376-388，1992中被描述。下面为使用一直到二阶导数的CDI：

(L, L_{xx} + L_{yy}, L_{x}^{2} + L_{y}^{2}, L_{x}^{2} L_{xx} + 2 L_{xy} L_{x} L_{y} + L_{y}^{2} L_{yy}, L_{xx}^{2} + 2 L_{xy}^{2} + L_{yy}^{2}) .

同样地，在每个像素位置，在多个内尺度σ使用这些算子并计算范围为α的窗口内局部无序直方图的头两个时刻，从而构造特征图像，其中α和内尺度相关联，例如α＝2σ。这些算子与图像几何的欧几里得变换无关，但仍然与所选择的尺度参数σ有关。

轮廓和特征相似性测量

灰度级外观模型是基于平均轮廓和方差-协方差矩阵而建立，并利用轮廓这些点之间的强度关联，这将在下面解释.轮廓处的强度由前面所勾画出的强度特征向量表征.

为了将一个轮廓和另一个轮廓比较，需要知道轮廓中每个点是如何与任何其它点关联的，并测量该轮廓的特定点的数值相对于任一其它点偏离平均值的程度。通过协方差测量捕捉这个统计关系，总是在两个点(或这些点处的特征)之间计算该协方差。如果测量某一个点和它自身的特征值之间的协方差，则可以获得方差.对于总的点数为2k+1时，若轮廓在轮廓的任意一侧包括k个点，则存在(2k+1)k个关系.来自轮廓的点对之间协方差告知点对的第一个数值如何相对于第二对发生变化.当协方差大且为正值时，表示当第一点的特征值增大时，第二点成比例地增大。当协方差大且为负值时，表示当第一点的特征值增大时，第二点成比例地减小。当协方差为零时，意味着这些特征值之间不存在系统耦合，即，这些特征值是相互独立的.所有方差-协方差对可以排列成协方差矩阵S，该矩阵沿主对角线对称.协方差矩阵因此捕捉了轮廓中特征值对之间的所有结构关联。实践中，协方差矩阵是由学习采样构建成的，这种情形下为许多标记处的轮廓汇集以及训练图像集中的许多强度特征.

协方差涉及如何在高维空间，即，给定轮廓中2k+1个点的空间内计算距离，从而计算当前轮廓和模型轮廓之间的距离.当轮廓只包括值为x的单个点时，通过将当前轮廓点和模型轮廓点的数值相减并取绝对值，即d＝|x-x|，就可以将这些轮廓和模型点数值x进行相似性比较。当轮廓包括两个点时，二元向量x和模型二元向量x之间的距离可以计算成平面内的欧几里得距离，即：

d = \sqrt{(x - \overset{&OverBar;}{x})' (x - \overset{&OverBar;}{x})} = \sqrt{{| x_{1} - {\overset{&OverBar;}{x}}_{1} |}^{2} + {| x_{2} - {\overset{&OverBar;}{x}}_{2} |}^{2}} .

在所有方向上衡量欧几里得距离都是等同的.然而，当变量x1和x2具有不等的标准差σ₁及σ₂时，等距离线的轮廓则为主轴平行于坐标轴的椭圆。可以使用标准差的倒数衡量这些差异，得到加权的欧几里得距离：

矩阵W包括变量的方差的导数，即

W = [\begin{matrix} 1 / σ_{1}^{2} \\ 1 / σ_{2}^{2} \end{matrix}] .

在线性或圆形轮廓中相互靠近的像素的强度水平，当其属于具有平滑变化强度水平的相同解剖区域时，将是强关联的。当像素对被认为属于不同组织时，关联将成反比。在任一情形中，成像噪声和不重要的解剖变化会引入像素自身强度水平的非零方差.

因此，如果方差也是关联的并具有不同方差，则必需插入协方差矩阵的逆矩阵，得到所谓的马氏距离：

d_m＝√(x-x)′S^-1(x-x).

马氏距离衡量多维数据点x与其平均值x的距离，使得在相同的多变量正常密度轮廓上的观察将具有相同的距离d_m.这些马氏等距离轮廓在二维空间内呈椭圆的形式，在三维空间中呈椭圆体。显然，轮廓2k+1中点的数目通常远大于3，这种情况下等距离轨迹在多维空间中呈超椭圆体的形式。

灰度级外观模型是从训练图像集中获得的，包括平均轮廓中的每个标记以及每个特征的协方差矩阵。所储存模型的总数因此正比于标记的数目和特征的数目。

形状模型

在形状模型化步骤中，考虑了沿轮廓的标记之间的位置关联。由于对具有相对确定性形状的解剖结构的分割是针对有关该形状必须优选编码该确定性成分的知识，以及相关的解剖变形，这样进一步优选地忽略了不相关的形状细节。

描述了两种实施例来考虑位置关联。第一个实施例是基于涉及同时使用所有标记位置的整体法；第二个实施例采用局部法，模型化连续标记的空间部署。

PCA分析和拟合

在第一实施例中，由手动放置于排列训练图像集中的标记构成标记坐标位置的方差-协方差矩阵。该方差-协方差矩阵中的每个对角元代表这些标记的坐标值的方差.非对角元代表标记坐标值对的协方差，即，表示两个标记的坐标是如何彼此相对变化的。该变化可以系统地相对其每个平均值沿相同的“正方向”而形成大的正协方差，或者可以系统地相对其每个平均值沿相同的“负方向”而形成大的负协方差，或者相对于其每个平均值随机地沿相对的方向而产生零协方差。现有技术主成分分析(PCA)中已知的本征向量/本征值分析将检测这些标记相对于平均标记位置的任何全局系统变化模式.PCA旨在通过只保留具有最大变化的坐标分量而减小数据集合的维数，因此PCA对原始坐标位置进行统计学地去关联(decorrelate)。通过投影到这些高度变化的子空间，可以由保持重要形状特征的维数更少的系统近似这些相关的统计量。足够数量的这些主要模式被保留并用于形状模型，忽略了由于噪声变化引起的其它模式。

标记连接向量模型

第二实施例旨在模型化在局部水平布置连续标记。然而，这种约束并不意味着不捕捉全局形状.给出起始点坐标以及这些连接向量，从该起始点开始，通过迭代地将连接向量加到当前的点，就可以完整地重构形状轮廓。在像素水平的形状描述中，由例如链式码利用这种性能。该连接向量为该曲线的局部切向量的离散近似.当考虑到连续标记之间的方向向量时，即一阶几何导数，可以描述对象的形状而与对象的平移无关。通过考虑例如由两个连续连接向量夹角表达的连续标记间的曲率，即二阶几何微分，则可以描述对象的形状而与该形状的平移和旋转无关(所谓的刚性体变换).这种情况下，当起始点和一阶几何导数(即，起始连接向量)已知时，通过迭代地累加连续切向量的差向量，可以完整地重构该曲线。基于像素水平的链式码的曲率通常对噪声是高度相关的，这是因为其尺度处于极精细的水平。这种缺点在本发明公开中得到解决，这是因为标记间隔足够远，使得可以实现与噪声无关并且精确地计算一阶或二阶导数信息.由于切向量或曲率向量的大小与尺度有关，将训练集的每个曲线的大小规一化成例如单位元素(unity)，可进一步将曲线描绘成相当于欧几里得相似变换。排列训练图像集中手动放置的标记之间的平移向量(即，一阶导数或切向量)或平移向量的向量差(即，二阶导数或曲率)的平均和协方差统计构成了本实施例的形状模型。

分割系统

分割系统也可以根据每个包括具有相关目标的构建单元的两种路径进行组织。

上述模型可以用于分割系统.

分割系统也可以根据分别包括具有相关目标的构建单元的两种路径进行组织。

第一构建单元旨在将示例模型标记位置的容许位置限制在目标图像中的多个最可能位置。

对于使用特征图像中线性轮廓代表标记的情形，可以采用最佳点表决方法，该方法在穿过示例标记的垂直线上对候选位置集进行分级。

对于圆形轮廓的情形，可以在每个示例标记和栅格交点位置考虑矩形搜寻栅格(search grid)，并根据穿过所有特征图像的圆形轮廓到模型圆形轮廓的组合马氏距离对栅格交叉位置进行分级。

第二和第三构建单元共同组成优化步骤。使用与每个标记的保留容许位置相关联的成本测量，可以构造成本矩阵.该成本包括与灰度级外观中一致性相关联的项。可以使用加权PCA轮廓拟合在全局水平校验与形状模型的一致性。备选地，可以通过结合与标记候选位置相关联的形状成本，以验证该形状一致性。动态编程算法可以用于构造穿过该成本矩阵的最小成本路径，进而得到路径具有最低总成本的最终分割。

最后，最终分割上的特征点可间接地定义关键测量点的几何位置，且可以基于这些关键测量点产生所有相关测量对象和测量数量.

备选地，分割中使用的标记点可用于直接定义测量点，且可提供每个测量点的搜寻栅格以定位测量点的允许位置。

通过下述描述和附图，本发明另外的具体实施例和相关优点将变得显而易见。

附图说明

图1为说明根据本发明的第一模型化系统(模型化系统I)的流程图。

图2为说明根据本发明的第一分割系统(分割系统I)的流程图。

图3为说明根据本发明的第二模型化系统(模型化系统II)的流程图。

图4为说明根据本发明的第二分割系统(分割系统II)的流程图。

图5为说明基于根据本发明被分割的解剖对象上的特征点测量计算的流程图。

图6为用于根据胸部放射照片测量心胸比率的测量方案.

图7(a)说明了胸部放射照片中左肺的手动描绘，在图7(b)中使用标记集表达该左肺，本示例中n＝14。

图8说明了使用n个标记集p₁＝(x₁，y₁)，…，p_n＝(x_n，y_n)对肺野的离散表达。该图像尺寸被规一化在0...1之间.通过该离散标记插入连续曲线。

图9说明了分别在标记附近的直线上和半径为r_c圆上的多个点的灰度值采样的线性以及圆形灰度级轮廓。

图10示出了包括14个标记的左肺形状上一些连接向量的实验评估分布：(a)v₁，(b)v₂，(c)v₃。

图11示出了胸部放射照片的第一时刻特征图像，对应于一直到二阶的导数(L，L_x，L_y，L_xx，L_yy，L_xy)以及等于0.5、2和8像素的尺度σ.

图12示出了胸部放射照片的第二时刻特征图像，对应于一直到二阶的导数(L，L_x，L_y，L_xx，L_yy，L_xy)以及等于0.5、2和8像素的尺度σ。

图13示出了将本发明应用于胸部放射照片上肺野的分割的示例：(a)手动分割左肺野和右肺野，(b)自动分割左肺野和右肺野。

图14说明了心胸比率的计算：(a)使用手动肺野分割得到CTR_man＝0.47，(b)使用根据本发明的自动导出肺野分割得到CTR_automatic＝0.45。

具体实施方式

将参考具体应用，即医学图像中肺野的分割，详细地解释本发明。

对象表示

在下述本发明方法的具体实施例中，图像中的解剖对象被数学地表示成位于封闭该对象轮廓上的固定数目的离散的标号点，即p₁＝(x₁，y₁)，…，p_n＝(x_n，y_n)。

轮廓从p₁前进到p_n并返回到_p1。因此，可以通过离散形状向量x＝(x₁，y₁，…x_n，y_n)^T捕获该对象。选择坐标系统，使得该图像区域内的所有点都落在域[0，1]×[0，1]内(图7)。

下述分割方案需要许多训练图像，其中手动确定形状向量x。一旦在该数据集上训练该算法，则该算法可以在新图像内产生该形状向量。可以由来自点

的曲线插值重构解剖对象的连续轮廓。一阶插值使用线段连接连续的点，形成多边形近似。高阶插值可以用于获得更平滑的轮廓。与阶数无关地，该连续表示可以用于导出曲线的空间性能，例如沿曲线的切线方向，或者曲线上特定点处的法线方向。

模型化系统I(图1)

灰度级外观模型

灰度级外观模型(也称为光度模型)捕捉在所谓标记(也称为标记点)附近范围内的空间强度分布。通过使用灰度级轮廓已经实现了该灰度级外观模型的空间分量。从其构造可知，可以采用其它方案。

例如被用于精确地描述肺野轮廓的标记的数目n可以仅为14，例如标记1表示每个肺野的最高点，标记7置于左心或右心阴影及其相应半横隔膜的截面处，标记9表示肋膈角。其它标记可以等距离地置于这些主要标记之间.显然，可以使用更多数目的标记。

线性灰度级轮廓(图9)

在每个标记，描述标记周围的典型图像结构的灰度级外观模型垂直于轮廓地被采样。在标记的任何一侧，使用特定步进尺寸采样k个像素，这给出了长度为2k+1的轮廓向量。这种采样方案的优点为能够模型化在该标记处的线性强度梯度。这些轮廓可以被归一化，使得轮廓中元素的绝对值的总和为1.采样的方向可以被计算为将标记与其前一个标记或后一个标记连接而得到线段上的法线之间的平均法线方向，或者计算为连接该前一个标记和后一个标记的线段上的法线.

灰度级特征图像

标记周围呈现的图像结构被捕捉为N个数学特征，这里所解释的数目更高。通过在每个位置x₀周围宽度为α的窗口内计算导数图像在多个时刻的局部直方图，由此获得这些特征图像.在优选实施例第一和第二直方图时刻中，计算一直到二阶的所有导数(L，L_x，L_y，L_xx，L_yy，L_xy)和5个内尺度(σ＝0.5、1、2、4、8像素)，得到总数为N＝2×6×5＝60的特征图像。计算直方图的该窗口尺度是根据内尺度，即α＝2σ.特征图像集可以被看作是灰度值函数在标记附近的灰度值分解。

胸腔图像内的计算特征图像的示例在图11(原始图像以及5个导数图像在三个尺度的第一直方图时刻)和图12(原始图像以及5个导数图像在三个尺度的第二直方图时刻)中示出。胸腔图像通常在图像中显示成处于竖直位置，关于其它放射检查的解剖对象可能呈现非标准的位置以及旋转，这种情况下可以采用笛卡儿微分不变式构造特征图像。

灰度级外观模型

训练集的所有图像被用于建立每个特征和每个标记的统计模型。将在特征j(例如2D图像中其总数为60)中在标记i(例如胸腔肺野的2D情形中其总数为14)采样的归一化轮廓示为向量g_i，j，从训练图像评估g_i，j的概率分布，得到平均值g_i，j和协方差S_i，j。对于长度为2k+1的轮廓向量，该协方差矩阵的大小为(2k+1)×(2k+1).特定标记i的灰度级外观模型包括平均轮廓g_i，j，j＝1...N和协方差矩阵S_i，j，j＝1...N。总肺野的灰度级外观模型包括所有平均轮廓g_i，j，i＝1...n，j＝1...N和协方差矩阵S_i，j，i＝1...n，j＝1...N的集合.

灰度级成本

在新特征图像j中点p处为特定标记i采样的新轮廓g_i，j之间的马氏距离计算成：

h_{i, j} (p) = {(g_{i, j} (p) - {\overset{&OverBar;}{g}}_{i, j})}^{T} S_{i, j}^{- 1} (g_{i, j} (p) - {\overset{&OverBar;}{g}}_{i, j})

马氏距离越小意味着轮廓g_i，j(p)源于平均值为g_i，j且协方差为S_i，j的高斯分布的可能性越大.因此，马氏距离可用作灰度级成本测量，用h_i，j(p)表示。该成本测量为给定图像中位置p的函数。最可能根据特征j的标记i的真实位置的位置为使h_i，j(p)最小的点p。这样就可以为每个特征j定义灰度级成本。

形状模型

由n个标记点描述曲线轮廓，在肺野的情形中每个肺野具有一个曲线轮廓(图8)。在s个训练图像集合中手动地确定这些标记点，得到一系列点(x₁，y₁)，…，(x_n，y_n).这些坐标元组被依次置于表示该曲线的向量x＝(x₁，y₁…，x_n，y_n)^T中。接着，对这些训练图像的形状向量x进行主成分分析(PCA)。PCA将曲线投影到覆盖了肺野的几何变形的最重要模式的更低维空间内。可以使用

t<<2n近似每个曲线

x≈x+Φb

其中x为平均形状，为与第t最大本征值相对应的协方差矩阵x的本征向量。每个本征值决定该形状模式的形状的方差。由向量b表示的该曲线近似组成了该形状模型，并将曲线x拟合到该形状模型.这些本征形状可以被看作是由标记集表示的形状的零阶(位置的)几何分解.形状模型的作用是将对象的形变约束在由训练集设置的变化界限之间，例如相对于平均形状的三个标准偏差。

分割系统I(图2)

分割算法最初将把代表肺野轮廓的曲线放置于图像内，并通过迭代地更新曲线组成标记的位置直至优化标准达到最佳值为止而找到分割解决办法。每个迭代包括下述步骤。

步骤1.最佳点表决

假设具体标记的当前位置，将在垂直于轮廓的线上寻找该标记的更佳定位。在该标记的当前位置的各侧上评估了n_s个位置.这些2n_s+1个点之一将被表决成最佳下一个点。N个特征集合的每个特征对该2n_s+1个点之一贡献一个表决。该2n_s+1点中被选择的点根据马氏距离标准具有最小的成本。对于每个该2n_s+1点，由特征图像j中被采样的2k+1个点组成的轮廓被置于h_i，j(p)方程中，且该方程中填充了恰当的平均g_i，j和协方差S_i，j。由该特征来表决具有最小马氏距离的点。每个特征选择一个点，导致在2n_s+1个点即集合

(其中

Σ_{k = 1}^{2 n_{s} + 1} N_{k} = N

)上划分N个表决.显然，特征的数目N必需远大于可选择的点的数目2n_s+1，否则可能分配到可选择点集的点的表决太少.

特定标记i所考虑的2n_s+1个点的集合中的每个点接收到的表决数N_i，j，可以被看作是该具体点正是标记所搜寻的点的置信等级.

步骤2.最小成本路径优化

根据第一步骤，每个标记可以分离地置换到2n_s+1个可能的位置.可将置信等级N_i，j置于如下矩阵中：

R = [\begin{matrix} N_{1,1} & \cdot \cdot \cdot & N_{1, {2 n}_{s} + 1} \\ \cdot & \cdot \\ \cdot & N_{i, j} & \cdot \\ \cdot & \cdot \\ N_{n, 1} & \cdot \cdot \cdot & N_{n, {2 n}_{s} + 1} \end{matrix}],

其中

\begin{matrix} Σ_{j = 1}^{{2 n}_{s} + 1} N_{i, j} = N & &ForAll; i = 1, \cdot \cdot \cdot, n . \end{matrix}

更新轮廓涉及选择矩阵C各行内的一个元素；最小成本路径(MCP)搜寻是指，将形成贯穿该矩阵的路径的所有元素连接，使得成本标准最优化(即最小化)。

一种合乎逻辑的选择是取各行中置信最高的元素，这种途径的缺点为它不考虑存在界外值(outlier)的情形.因此可能包括路径内的陡变.通过施加排除该陡变的平滑约束，这个缺点得到克服.

为此目的，首先如下所述地构造成本矩阵：

C = [\begin{matrix} 1 / N_{1,1}^{m} & \cdot \cdot \cdot & 1 / N_{1, {2 n}_{s} + 1}^{m} \\ \cdot & \cdot \\ \cdot & 1 / N_{i, j}^{m} & \cdot \\ \cdot & \cdot \\ 1 / N_{n, 1}^{m} & \cdot \cdot \cdot & 1 / N_{n, {2 n}_{s} + 1}^{m} \end{matrix}],

其中，分母中的幂m为影响计算路径的正数(m越大则对点的成本贡献越小，因此指导路径朝向该点)。其次，首先如下所述地构造成本函数：

J (j_{1}, j_{2}, \cdot \cdot \cdot, j_{n}) = Σ_{i = 1}^{n} C_{i, j_{i}} + α Σ_{i = 1}^{n} | δ i, j_{i} - δ_{i - 1, j_{i - 1}} |,

其中为点i朝其目标点j_i沿该形状法线的位移，α为加权因子.该方程中的第二项为平滑约束，该约束处罚如下情形的目标点配置：一个目标点相对于(w.r.t.)前一目标点的位移沿该轮廓陡变的位移。使用例如D.H.Ballard，C.M.Brown，ComputerVision，Englewood Cliffs，PrenticeHallInc.1982，pp137-143中所描述的现有技术中已知的动态编程技术，优化成本函数J，得到使J最小化的目标点

作为最低成本路径(MCP)。

作为矩阵数据结构的备选，还可以使用分层曲线图。各层代表在标记的垂直采样图像点阵列。使用一组弧从最后一层返回到第一层，由此增加该曲线图.动态编程技术提取边界作为穿过该曲线图的最佳闭合路径。

步骤3.未加权和加权的PCA轮廓拟合

在上一步骤中，每个标记i被赋予新位置

序列

表示通过将其替换成x并通过求解Φb＝x-x导出b而可以被拟合到形状模型x≈x+Φb的轮廓。在该形状模型中拟合x是指找到b值，使得x+Φb在高置信标记中可以小误差地近似x.

未加权的PCA轮廓拟合

由于

因此方程数目多于未知量，该方程组为超定方程组，没有精确解.因此，‖Φb-(x-x)‖_p必需被最小化，以恰当地选择p。不同的规范得到不同的最佳解。具体地p＝2的情形，即最小二乘方问题，更容易处理，在现有技术中采用基于标准方程和QR因式分解的方法可以求解该方程。将实际点序列和拟合点序列之间的差异用拟合误差e＝x-(x+Φb)表示，则问题变为寻找最小化e^Te的b。替换e，该优化问题变为：

\min_{b} J (b) = \min_{b} e^{T} e = \min_{b} ({(x - \overset{&OverBar;}{x})}^{T} - b^{T} Φ^{T}) ((x - \overset{&OverBar;}{x}) - bΦ) .

如果

则J(b)最小化.J(b)的梯度可以写成：

&dtri; J (b) = - 2 Φ^{T} W^{2} (x - \overset{&OverBar;}{x}) + 2 Φ^{T} W^{2} Φb .

因此，由优化方程得到已知的标准方程：

&dtri; J (b) = 0 &DoubleLeftRightArrow; Φ^{T} Φb = Φ^{T} (x - \overset{&OverBar;}{x}),

其中求解上述方程得到b：

b＝(Φ^TΦ)^-1Φ^T(x-x)

该方程将图像空间的数据向量x投影到PCA模型空间并返回近似系数b。可以通过评估形状模型x＝x+Φb而将近似结果b投影回到图像空间。

加权PCA轮廓拟合

当将轮廓拟合到形状模型时，相应数目的表决

可以用作加权.因此，PCA模型将趋于赋予高置信点高的优先级。置信等级低的点对该模型产生的轮廓的影响减小。优化问题现在则包括最小化加权二次误差(We)^T(We)。使用和未加权情形中采用标准方程的相似方法，该问题的解b可表达成：

b＝(Φ^TW²Φ)^-1Φ^TW²(x-x)

计算该方程以将x投影在PCA空间内，且使用形状模型x＝x+Φb可以往回投影。

迭代步骤1至3，直到满足停止判据，例如在连续迭代之间标记位置基本上不再改变时。

未加权和加权PCA拟合都得到图像内标记点集x的最终分割结果.可在这些点之间插入连续样条曲线以解析地表示该分割.

模型化系统II(图3)

灰度级外观模型

灰度级外观模型捕捉标记邻城内的空间强度分布.在本实施例中，空间模型偏离在标记邻城内采样的圆形轮廓.用于精确地描述肺野轮廓的标记的数目n可以仅为14，例如标记1表示每个肺野的最高点，标记7置于左心或右心阴影及其相应半横隔膜的截面处，标记9表示肋膈角。其它标记可以等距离地置于这些主要标记之间。显然，在根据本发明的计算方案中还可以使用更多数目的标记。

圆形轮廓

在标记邻域中选择点的一个备选方式为，在以标记为中心并具有至少一个半径的圆上对点进行采样。图9中给出了圆形采样的示例.如果标记位于(x₀，y₀)，则在半径为r_c的如下各点处采样该图像的灰度值函数：

\begin{matrix} (x_{0} + r_{c} \cos \frac{2 π}{n_{c}} (i - 1), y_{0} + r_{c} \sin \frac{2 π}{n_{c}} (i - 1)) & i = 1, . . ., n_{c} . \end{matrix}

采样被置于长度为n_c的轮廓向量中.n_c的适当选择为12、8、6、4和3(对应于30、45、60、90和120度角间隔)。可以同时考虑在一系列半径处的多个圆形子采样.表达成与图像尺寸之比的无量纲量分数的适当的半径r_c为0.05和0.025.

和线性轮廓相比，该方案的优点为在邻域周围捕捉2D结构.和线性轮廓相比，圆形轮廓可以更好地模型化诸如肋膈角落点或左心阴影与左横隔膜交点的特殊解剖标记，因为这些标记具有肺野轮廓的不连续一阶导数。在这些特殊解剖标记附近区域中沿肺野轮廓的标记，如果其长度太长，则还可能被线性轮廓模糊地表示，因为该轮廓可能跨过并列的肺野边界。

圆形轮廓可以被归一化，使得轮廓中元素的绝对值的总和为1。子采样的角度可以从无子采样(得到原始分辨率的该图像的再次采样)变化到对于圆形轮廓的情形，在主轴上例如仅仅4点的子采样.

灰度级特征图像

标记周围的图像结构被捕捉为N个数学特征，正如本发明的发明内容中所述。通过在每个位置x₀周围宽度为α的窗口内计算导数图像在多个时刻的局部直方图，由此获得这些特征图像.在优选实施例的第一和第二直方图时刻中，计算一直到二阶的所有导数(L，L_x，L_y，L_xx，L_yy，L_xy)和5个内尺度(σ＝0.5、1、2、4、8像素)，得到总数为N＝2×6×5＝60的特征图像。计算直方图的该窗口尺度是根据内尺度，即α＝2σ.特征图像集可以被看作是灰度值函数在标记附近的灰度值分解.

灰度级成本

在新特征图像j中为特定标记i采样的新轮廓g_i，j之间的马氏距离计算成：

h_{i, j} (p) = {(g_{i, j} (p) - {\overset{&OverBar;}{g}}_{i, j})}^{T} S_{i, j}^{- 1} (g_{i, j} (p) - {\overset{&OverBar;}{g}}_{i, j})

马氏距离越小意味着轮廓g_i，j(p)源于平均值为g_i，j且协方差为S_i，j的高斯分布的可能性越大.因此，马氏距离可用作灰度级成本测量，用h_i，j(p)表示。该成本测量为给定图像中位置p的函数.根据特征j最可能为标记i的真实位置的位置为使h_i，j(p)最小的点p。这样就可以为每个特征j定义该灰度级成本。

总体灰度级成本可通过组合所有特征j的所有灰度级成本的总和而得到，即：

h_{i} (p) = Σ_{j = 1}^{N} {(g_{i, j} (p) - {\overset{&OverBar;}{g}}_{i, j})}^{T} S_{i, j}^{- 1} (g_{i, j} (p) - {\overset{&OverBar;}{g}}_{i, j})

该总和反应了位置p处的灰度级图形和标记i处的预计灰度级图形之间的相似性。根据灰度级外观，最可能与标记i一致的位置为下述优化问题的解：

p_{0} = \arg \min_{p} h_{i} (p) .

为了构造灰度级外观模型，必需从训练图像中评估特征轮廓g_i，j(对于所有标记i＝1，…，n，对于所有特征图像j＝1，…，N)的分布(g_i，j，S_i，j)。

形状模型

尽管灰度级成本验证了灰度级图形或特征的图形，但仍定义了形状成本以验证两个连续标记之间的连接。

由n个标记点描述曲线轮廓，每个肺野具有一个曲线轮廓。在s个训练图像集合中手动地确定这些标记点，得到一系列点(x₁，y₁)…(x_n，y_n)＝(p₁，...，p_n).这些坐标元组被依次置于表示该形状的向量x＝(x₁，y₁…，x_n，y_n)^T中.考虑连续标记的评估位置对(p_i，p_i+1)。将形状成本赋予向量v_i＝p_i+1-p_i，这反应了v_i的真实程度，即其概率分布。通过训练形状评估v₁，v₂，…，v_n的分布，假设在其平均值附近具有正态分布.v_i的平均值和协方差分别用v_i和表示.在图10中给出了这些向量分布的示例.这些向量分布可以被看作是由该标记集代表的形状的一阶(正切)集合分解。

验证两个连续标记p_i，p_i+1之间的连接向量v_i(在平面内可完全由其取向和长度表征的向量)的一种新颖形状成本为v_i与其平均值v_i之间的马氏距离：

f_{i} (v_{i}) = {(v_{i} - {\overset{&OverBar;}{v}}_{i})}^{T} S_{v_{i}}^{- 1} (v_{i} - {\overset{&OverBar;}{v}}_{i}) .

由于具有大的马氏距离而不可能发生的连接将具有高的形状成本。另一方面，零成本将被赋予等于预期值的连接.

分割系统II(图4)

在训练阶段期间采集有关图像中对象的灰度级外观的知识以及有关待分割对象的形状的知识.这些知识随后被用于根据下述步骤在新图像中分割对象。

步骤1.矩形搜寻栅格构造

构造每个标记i的矩形搜寻栅格，以约束该标记的可能位置.每个矩形栅格的特征在于几何范围以及栅格间隔的参数。

栅格间距r_g应该与灰度级外观模型的圆形轮廓的半径r_c相对应.适当的关系为r_g＝f_gr_c，其中r_c为半径，f_g为固定分数.栅格间距的典型分数为f_g＝0.5。

栅格范围和栅格位置是针对每个标记的，由x_min、x_max、y_min和y_max表示。真实(未知)标记位置(x^*，y^*)被认为是位于下界和上界之间：

x_min<x<x_max

y_min<y<y_max

只有当该栅格跨越整个图像区域时才能满足这些条件(即，x_min＝y_min＝0以及x_max＝y_max＝1).更高效的方法是将搜寻约束在图像的相关部分。假设该标记的x坐标的概率分布p_x(x)是平均值为x且标准偏差为σ_x(从训练形状评估)的正态分布。在x方向上共同定义栅格范围的栅格边界x_min和x_max被选择成使得：

{&Integral;}_{x_{mim}}^{x_{\max}} p_{x} (x) dx = f_{c}

其中f_c是代表条件x_min<x^*<x_max为真的概率的分数.分数f_c的恰当值为0.995。如果施加的条件为该范围围绕平均值x是对称的，即x_max-x＝x-x_min，则上界x_max的值满足：

\frac{1}{σ_{x} \sqrt{2 π}} {&Integral;}_{\overset{&OverBar;}{x}}^{x_{\max}} \exp (- \frac{{(x - \overset{&OverBar;}{x})}^{2}}{2 σ_{x}^{2}}) dx = \frac{1}{2} f_{c} .

下界则为x_min＝2x-x_max.类似地获得y坐标的边界。

步骤2.灰度级成本矩阵构造

灰度级外观模型被用于寻找每个标记的恰当位置.如下所述地根据灰度级外观模型选择标记i的m个最佳位置.

首先，根据前述步骤定义覆盖图像区域相关部分的矩形栅格，该图像区域相关部分在每个标记i的预计位置附近。

其次，对于每个栅格点，计算总体灰度级成本h_i(p)。选择和m个最低总体灰度级成本相对应的点p_i，1，p_i，2，…，p_i，m。

第三，构造灰度级成本矩阵C，使得每行i包括标记i的m个最可能位置的成本：

C [\begin{matrix} h_{1} (p_{1,1}) & \cdot \cdot \cdot & h_{1} (p_{1, k}) & \cdot \cdot \cdot & h_{1} (p_{1, m}) \\ \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot \\ h_{i} (p_{i . 1}) & \cdot \cdot \cdot & h_{i} (p_{i, k}) & \cdot \cdot \cdot & h_{i} (p_{i, m}) \\ \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot \\ h_{n} (p_{n, 1}) & \cdot \cdot \cdot & h_{n} (p_{n, k}) & \cdot \cdot \cdot & h_{n} (p_{n, m}) \end{matrix}] .

每个标记典型的最佳点数目m为10至40.

由于每个标记的m个最佳点的选择与邻近标记的m个最佳点集合相互独立，因此会出现如下情况，即一个或多个这些点更靠近非邻近的标记。通过在下一步骤中考虑形状成本，这些点可能在最后的轮廓分割中被忽略。

步骤3.最小成本路径构造

确定分割对象的轮廓被简化为通过选择每行的一个元素在矩阵C中彻底地寻找路径。将第i行中选定元素的指数写成k_i，则曲率变为点序列

最佳路径是指使成本函数J(k₁，…，k_n)最小的路径：

(k_{1}^{*}, k_{2}^{*}, . . ., k_{n}^{*}) = \arg \min_{k_{1}, . . ., k_{n}} J (k_{1}, . . ., k_{n}) .

上面提出的模型允许进行大量的成本测量。考虑两个连续标记和

根据灰度级外观模型的成本分量以及形状模型的成本分量为：

-与标记及相对应的灰度级成本

和

-验证从到连接的真实性的形状成本

通过对总体灰度级成本和总体形状成本的加权求和，将这两种成本类型都结合到成本函数J(k₁，…，k_n)内。加权因子为γ₁的总体灰度级成本是标记的单个灰度级成本之和

h_{i} (p_{i, k_{i}}), i = 1, \cdot \cdot \cdot, n

加权因子为γ₂的总体形状成本是各形状成本之和

f_{i} (p_{i + 1, k_{i + 1}} - p_{i, k_{i}}), i = 1, \cdot \cdot \cdot, n .

因此成本函数J(k₁，…，k_n)变为：

J (k_{1}, . . ., k_{n}) = γ_{1} Σ_{i = 1}^{n} h_{i} (p_{i, k_{i}}) + γ_{2} Σ_{i = 1}^{n - 1} f_{i} (p_{i + 1, k_{i + 1}} - p_{i, k_{i}}) + γ_{2} f_{n} (p_{1, k_{1}} - p_{n, k_{n}}) .

使J(k₁，…，k_n)最小的最佳路径称为最小成本路径.使用例如由G.Behiels等在 Evaluation of image fea tureands earchstrategies for segmentation ofbone structures usingactiveshape

图13示出了根据前述系统的三个胸腔图像的自动分割，并将其与由经验的放射科医生的手动分割进行比较。示出了每个肺野上的14个自动放置的标记，同时还示出了穿过这些标记的连续样条曲线.

计算加速改进

包括步骤1至3的分割方法可以被应用于分割一个或多个解剖实体的轮廓的仅仅一部分。为此目的，只为这些连续标记的子集构造矩形搜寻栅格。该灰度级成本矩阵将包括数目与残留标记的数目相等的行。该最小成本路径构造将最小化成本函数，该成本函数仅包括残留标记的灰度值项和形状项。最终分割仅跟踪被残留标记覆盖的轮廓部分.当一个标记仅在部分解剖实体轮廓中感兴趣时，或者当需要定位仅部分特殊标记以获得具体测量时，部分分割的明显优点为计算速度.

加速分割计算的另一个算法改进是使用由粗到细的多分辨率方法。基于每个标记的多分辨率表示而构造灰度级外观模型和形状模型，引入空间尺寸作为一个新的维数。为此目的采用了两种备选，这两种备选都在连续更精细分辨率水平上采用该分割概念.

在一个备选中，在分割期间，减少原始图像中轮廓内点的数目、搜寻栅格内点的数目以及沿轮廓的标记点的数目。通过在前一分辨率的附近范围内连续地提高搜寻栅格的分辨率(轮廓和搜寻栅格中所考虑的点的总数因此可以在每次迭代时保持不变)并采用更精细分辨率灰度值和形状模型，于是可以提高这些标记的位置精度.

在另一个备选中，可以使用多个水平的高斯和拉普拉斯金字塔及其相关导数特征图像，实施该多分辨率方法.在粗糙水平上确定这些标记的初始位置，并使用前一水平的中间位置作为下一水平的开始位置，由此跟踪到最精细水平的最后位置。由于粗糙分辨率图像包括更少的细节，搜寻空间将可能包括更少的错误最小值，这使得可以在给定水平下实现更快的优化并可朝其最后位置更快地定位这些标记。

训练分割模型

从图像储存器收集许多胸腔图像，并将这些图像依次表示在图形用户接口上。执行下述步骤以建立分割模型：

—有经验的用户通过描绘肺野的轮廓(使用鼠标左击)标示出多个点而手动地分割每个被显示图像中的肺野.这些点无需沿肺野轮廓等间距地排布，唯一的要求就是这些点的集合整体上足够高程度地近似该肺野。为了评估解剖拟合，将样条曲线连续地插入穿过到目前为止已经确定的手动定位点，直到通过鼠标右击使该曲线从最后一个点闭合到第一个点。通过将单个点朝新位置拖曳，可以进行手动分割调整。再次评估结果轮廓的解剖准确性.

—接着，许多标记将被自动放置到手动分割的肺野轮廓上.为了实现训练集的所有图像上的相同标记彼此映射，用户定位几个容易识别的标记，通过等间距地放置于肺野轮廓上而获得其它标记.对于肺野轮廓的情形，多个容易识别的标记为肺野最顶点、指示肋膈角的点以及心脏阴影和半横隔膜交叉点，总共三个。接着选择全部标记，适当的选择范围为例如14至40.对于14个标记的情形，点p₁、p₇和p₉代表三个固定的标记，5个点平均地分布在p₁和p₇之间，1个点位于p₇和p₉之间，另外5个点位于p₉和p₁之间。分别对左肺野和右肺野执行这个步骤。

—接着，询问有关训练阶段的参数：(a)用于训练(和分割)的图像尺寸，典型值为256或512；(b)将由主成分分析(PCA)解释的形状变化的分数，典型值为0.9999；(c)线性或圆形轮廓中点的数目，圆形轮廓中的典型值为12、8、6、4或3；(d)作为图像尺寸分数的圆形轮廓半径，典型值为0.04、0.03和0.02。

—训练灰度级外观模型：(a)在这些标记附近分解该灰度值函数，即计算N个特征图像，例如如前所述的灰度值导数局部直方图时刻，和(b)收集线性或圆形轮廓在标记pi的灰度值分布，即计算g_i，j，i＝1…n，j＝1…N以及协方差矩阵S_i，j，i＝1…n，j＝1…N(即，对于特征图像j内的每个标记i)。这个步骤产生统计灰度级外观知识。

—训练形状模型：(a)计算包括在这些标记中的零阶(位置)信息的几何分解，即计算平均形状x以及在矩阵b中沿列方向排列的t个主要变形模式(本征形状)；(b)计算标记连接序列中包括的向量分布(一阶切线信息)，即计算v_i的平均值v_i和协方差这个步骤产生统计形状知识。

—储存统计灰度级外观和形状知识，例如根据上面给出的基于模型的分割子系统来分割新图像中的肺野.

将基于模型的2D图像分割应用于心胸比率(CTR)的计算

自动肺野分割的一个用途为计算心胸比率(CTR)。CTR(图6)定义为心脏的横向直径与胸腔的内径(ID)的比例：

CTR = \frac{MLD + MRD}{ID}

MLD为心脏左侧的最大横向直径，MRD为心脏右侧的最大横向直径。这个比率是重要的临床参数，成人的变化范围为39％至50％，平均值约为45％。高于50％的心胸比率被认为是异常的.可能的原因为心力衰竭、心包积液以及左或右心室肥厚.使用本发明中公开的自动肺野分割，可以自动地计算心胸比率。

参考图8，通过选择右肺分割的标记p₁和p₉之间的拟合轮廓片断上最靠左的笛卡儿点，从而获得定义MRD的特征点；通过选择左肺分割的标记p₁和p₇之间的拟合轮廓片断上最靠右的笛卡儿点，从而获得定义MLD的特征点。将这些特征点的列坐标相减并取绝对值，得到和MLD+MRD。类似地，通过下述步骤获得ID：选择右肺分割的标记p₁和p₇之间的拟合轮廓片断上最靠左的笛卡儿点以及左肺分割的标记p₁和p₉之间的拟合轮廓片断上最靠右的笛卡儿点；将这些特征点的列坐标相减；并取绝对值。

图14示出了分割肺野的特征点的计算，以及心胸比率的计算(a)使用手动肺野分割得到CTR_man＝0.47，(b)使用根据本发明的自动导出肺野分割得到CTR_automatic＝0.45。

将基于模型的分割和测量系统应用于其它身体部分

肺野分割中每个标记的搜寻栅格的空间范围是基于训练胸腔图像集中所有相似标记的位置而导出的.使用搜寻栅格用于约束给定标记的候选位置的概念可以扩展到其它身体部分，这些身体部分在图像中可能具有比胸腔图像中肺野更大的位置、旋转和尺寸变化，且和胸腔图像中肺野相反地占据了整个图像区域。

为了考虑到身体部分的这种位置、旋转和尺寸变化，使用欧洲专利申请第04076454号标题为“Method for automaticallymapping ofgeometric objects indigital medical images”中公开方法的搜寻栅格定位点映射的概念可以和本发明结合使用.这种改进的方法于是对于分割整形外科放射照片中多骨的身体部分是特别让人感兴趣的，因为其特征在于在图像中具有相对固定的形状，这些部分可以被定位到图像中非常明显的多个标记.

例如，在骨盆、髋部或整条腿检查中，股骨轮廓可以被定位到熟知的标记，例如稍大的转节(trochanter)的尖端、膝盖中心以及股骨头部中心。这些定位点被用于建立模型定位点和实际图像中选定的相关定位点之间的仿射变换。因为搜寻栅格包括收集排列在模型图像中分割模型轮廓上标记点附近矩形格子的候选点，通过将该变换应用于模型点的坐标，每个组成的候选点依次被映射到实际图像中.按照这个方式，在该点的最可能位置附近重构特定标记点的搜寻栅格。和前述方式相同地进行优化过程，即通过优化穿过候选位置选定组合的路径的组合灰度值和形状成本，每个搜寻栅格进行一次。具有最佳成本的路径为该身体部分的最终分割.其它骨头的示例为手和上肢骨头、脚和下肢骨头、骨盆和脊柱椎骨及结构。

满足这种类型的基于标记分割的其它身体部分为MR图像的3D CT中的软组织器官.另一方面，这里可以采用基于逐个切片的方法，这种方法确定每个切片上解剖的分割轮廓，并将所有切片的结果组合成3D分割表面。另一方面，可以采用将本发明的概念扩展到3D体积的完全3D方法，该方法建立并应用模型限制在3D搜寻栅格内的3D标记。肾、肺、心、肝、胃、脾、前列腺和大脑结构为可以应用本发明方法的器官示例。

在所有提及的应用示例中，测量点可以是基于分割标记的位置或者基于多个标记的组合(例如伸长骨头的皮层任意一侧上两个并列点的中点，一对该中点确定该骨头的解剖轴)，且可以将这些标记点随后输入到诸如EP A 1349098中公开的测量系统内.

Claims

1.一种分割数字医学图像中解剖实体的方法，包括下述步骤：

为所述图像内多个标记的每一个标记定义所述标记的初始位置；

采样所述初始位置周围的邻域，所述邻域包括所述标记的多个候选位置；

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

优化成本函数，该成本函数表达了所有候选位置的总体灰度级成本和总体形状成本的加权总和，其中，总体灰度级成本是标记的单个灰度级成本之和，总体形状成本是各形状成本之和；

将分割解剖实体定义为穿过所述候选位置的选定组合的路径，该组合的所述成本函数被优化，

其中和候选位置相关联的所述成本是通过对为所述医学图像计算的特征图像集的每个特征图像中所述候选位置的灰度级成本求和而获得的总体灰度级成本，

所述灰度级成本表达为马氏距离，所述马氏距离定义为由相应协方差矩阵的逆矩阵加权的所述候选位置处的灰度值轮廓与相应的平均轮廓之间的距离，所述平均轮廓和协方差矩阵是从与所述解剖实体相关联的灰度值模型中提取的，

其中通过下述步骤获得所述灰度值模型：

在多个标记点处采样所述解剖实体的手动分割剖面；

在每个所述标记点的邻域中采样多个点；

对于每个标记点，将采样点排列在轮廓标记点周围的邻域内；

计算至少一个特征图像；

计算每个标记点和每个特征图像的平均轮廓；

计算每个标记点和每个特征图像的轮廓的协方差矩阵，

其中特征图像为预定尺度σ处的导数图像，或者为包括标记位置周围宽度为α的窗口内所述图像或导数图像的局部直方图的数学时刻的图像表示。

2.根据权利要求1的方法，其中所述总体灰度级成本为多个特征图像中轮廓的所有单个灰度级成本的组合，所述单个成本指马氏距离，所述马氏距离定义为由协方差矩阵的逆矩阵加权的所述轮廓与该特征的平均轮廓之间的距离，所述平均轮廓和协方差矩阵是从与所述解剖实体的灰度值模型中提取的。

3.根据权利要求2的方法，其中通过下述步骤获得所述灰度值模型：

在多个标记点处采样所述解剖实体的手动分割剖面；

在每个所述标记点的邻域中采样多个点；

计算至少一个特征图像；

计算每个标记点和每个特征图像的平均轮廓；

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

确定所述平均轮廓和协方差矩阵为所述解剖实体的灰度值模型。

4.根据权利要求1的方法，其中所述总体形状成本为所有单个连接向量成本的组合，所述连接向量连接连续标记，对于所有标记，所述单个形状成本由此表达为马氏距离，所述马氏距离定义为由协方差矩阵的逆矩阵加权的两个连续标记之间连接向量与平均连接向量之间的距离，所述平均连接向量和协方差矩阵是从所述解剖实体的形状模型中提取的，

其中通过下述步骤获得所述形状模型：

在多个标记点处采样所述解剖实体的手动分割剖面；

计算连接向量或者连续标记点之间的连接向量差；

计算连续标记点对的平均连接向量或平均连接向量差；

计算连接向量或连接向量差的协方差矩阵；

确定所述平均连接向量和协方差矩阵为所述解剖实体的几何模型。

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

6.根据权利要求1的方法，其中通过选择具有最小总体灰度值成本的位置而减小所述多个候选位置，所述总体灰度值成本为所述邻域中候选位置的所有灰度值成本的所有特征图像的总和，所述灰度值成本被表达为马氏距离，所述马氏距离定义为由协方差矩阵的逆矩阵加权的候选位置处的灰度值轮廓与相应的平均轮廓之间的距离，所述平均轮廓和协方差矩阵是从与所述解剖实体相关联的灰度值模型中提取的。

7.根据权利要求1的方法，其中所述邻域包括矩形栅格采样点。

8.根据权利要求1的方法，其中所述邻域包括圆形轮廓。

9.根据权利要求1的方法，其中所述灰度值模型和所述几何模型是从图像的多分辨率表示构造而成的，并被应用于所述图像的多分辨率表示。

10.一种测量数字医学图像中解剖实体的方法，包括步骤：

使用根据权利要求1的方法计算所述解剖实体的分割；

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

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

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

11.一种测量数字医学图像中解剖实体的方法，包括步骤：

使用根据权利要求1的方法计算所述解剖实体的分割；

确定所述标记为测量点；

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

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

12.一种分割数字医学图像中解剖实体的方法，包括下述步骤：

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

将分割解剖实体定义为穿过所述候选位置的选定组合的路径，该组合的所述成本函数被优化，其中

所述总体形状成本为所有单个连接向量成本的组合，所述连接向量连接连续标记，对于所有标记，所述单个形状成本由此表达为马氏距离，所述马氏距离定义为由协方差矩阵的逆矩阵加权的两个连续标记之间连接向量与平均连接向量之间的距离，所述平均连接向量和协方差矩阵是从所述解剖实体的形状模型中提取的，

其中通过下述步骤获得所述形状模型：

在多个标记点处采样所述解剖实体的手动分割剖面；

计算连接向量或者连续标记点之间的连接向量差；

计算连续标记点对的平均连接向量或平均连接向量差；

计算连接向量或连接向量差的协方差矩阵；