数字医学图像中解剖实体的灰度值模型和/或几何模型的构造方法

具体实施方式

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

对象表示

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

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

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

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

模型化系统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>)>>->ver>>g>‾>>>i>,>j>>>)>>T>sup>>S>>i>,>j>>>->1>sup>>>(>>g>>i>,>j>>>>(>p>)>>->ver>>g>‾>>>i>,>j>>>)>>>s>$

马氏距离越小意味着轮廓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个点即集合N₁，…，N_2ns+1(其中 $>sup>>\Sigma >>k>=>1>>>>>2>n>>s>>+>1>sup>>>N>k>>=>N>>s>)上划分N个表决。$

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

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

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

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

$>>R>=>\begin{array}{}>>>>N>1,1>>>>.>.>.>>>>N>>1>,>>>2>n>>s>>+>1>>>>>>>.>>>>>.>>>>>.>>>>N>>i>,>j>>>>>.>>>>>.>>>>>.>>>>>>N>>n>,>1>>>>>.>.>.>>>>N>>n>,>2>>n>s>>+>1>>>>>>>,>>s>其中$ >>>\Sigma >>j>=>1>>>2>>n>s>>+>1>>>>N>>i>,>j>>>=>N>\forall >i>=>1>,>.>.>.>,>n>.>>s>$\end{array}$

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

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

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

$>>C>=>\begin{array}{}>>>1>/sup>>N>1,1>msup>>>>.>.>.>>>1>/sup>>N>>1>,>>>2>n>>s>>+>1>>msup>>>>>>.>>>>>.>>>>>.>>>1>/sup>>N>>i>,>j>>msup>>>>.>>>>>.>>>>>.>>>>>1>/sup>>N>>n>,>1>>msup>>>>.>.>.>>>1>/sup>>N>>n>,>>>2>n>>s>>+>1>>msup>>>>>>,>>s>\end{array}$

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

$>>J>>(>>j>1>>,>>j>2>>,>\dots >,>>j>n>>)>>=>>\Sigma >>i>=>1>>n>>>C>>i>,>>j>i>>>>+>\alpha >>>\Sigma >>i>=>1>>n>>|>>\delta >>i>,>>j>i>>>>->>\delta >>i>->1>,>>j>>i>->1>>>>>|>,>>>s>$

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

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

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

在上一步骤中，每个标记i被赋予新位置j_i^*。序列(j₁^*，…，j_n^*)表示通过将其替换成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>-ver>>x>‾>>)>>T>>->>b>T>>>\Phi >T>>)>>>(>>(>x>-ver>>x>‾>>)>>->b\Phi >)>>.>>s>$

如果J(b)＝0，则J(b)最小化。J(b)的梯度J(b)可以写成：

J(b)＝-2Φ^TW²(x-x)+2Φ^TW²Φb.

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

J(b)＝0Φ^TΦb＝Φ^T(x-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的如下各点处采样该图像的灰度值函数：

$>>>(>>x>0>>+>>r>c>>cos>>>2>\pi >>>n>c>>>>(>i>->1>)>>,>>y>0>>+>>r>c>>sin>>>2>\pi >>>n>c>>>>(>i>->1>)>>)>>i>=>1>,>.>.>.>,>>n>c>>.>>s>$

采样被置于长度为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>)>>->ver>>g>‾>>>i>,>j>>>)>>T>sup>>S>>i>,>j>>>->1>sup>>>(>>g>>i>,>j>>>>(>p>)>>->ver>>g>‾>>>i>,>j>>>)>>>s>$

马氏距离越小意味着轮廓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>)>>=>>\Sigma >>j>=>1>>N>>>>(>>g>>i>,>j>>>>(>p>)>>->ver>>g>‾>>>i>,>j>>>)>>T>sup>>S>>i>,>j>>>->1>sup>>>(>>g>>i>,>j>>>>(>p>)>>->ver>>g>‾>>>i>,>j>>>)>>>s>$

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

$>>>p>0>>=>arg>>min>\rho >>>h>i>>>(>p>)>>.>>s>$

为了构造灰度级外观模型，必需从训练图像中评估特征轮廓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₁＝p_i+1-p_i，这反应了v_i的真实程度，即其概率分布。通过训练形状评估v₁，v₂，…，v_n的分布，假设在其平均值附近具有正态分布。v_i的平均值和协方差分别用v_i和S_vi表示。在图10中给出了这些向量分布的示例。这些向量分布可以被看作是由该标记集代表的形状的一阶(正切)集合分解。

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

$>>>f>i>>>(>>v>i>>)>>=>>>(>>v>i>>->ver>>v>‾>>i>>)>>T>sup>>S>>v>i>>>->1>sup>>>(>>v>i>>->ver>>v>‾>>i>>)>>.>>s>$

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

分割系统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被选择成使得：

$>sup>>\int >>x>min>>>x>max>sup>>>p>x>>>(>x>)>>dx>=>>f>c>>>s>$

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

$>>>1>>>\sigma >x>>>2>\pi >>>sup>>\int ver>>x>‾>>>x>max>sup>>exp>>(>->>>>(>x>-ver>>x>‾>>)>>2>sup>>>2>\sigma >>x>2sup>>>)>>dx>=>>1>2>>>f>c>>.>>s>$

下界则为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{array}{}>>>>h>1>>>(>>p>1,1>>)>>>>.>.>.>>>>h>1>>>(>>p>>1>,>k>>>)>>>>.>.>.>>>>h>1>>>(>>p>>1>,>m>>>)>>>>>>.>>>>>.>>>>>.>>>>>.>>>>>.>>>>>.>>>>>.>>>>>.>>>>>.>>>>>>h>i>>>(>>p>>i>,>1>>>)>>>>.>.>.>>>>h>i>>>(>>p>>i>,>k>>>)>>>>.>.>.>>>>h>i>>>(>>p>>i>,>m>>>)>>>>>>.>>>>>.>>>>>.>>>>>.>>>>>.>>>>>.>>>>>.>>>>>.>>>>>.>>>>>>h>n>>>(>>p>>n>,>1>>>)>>>>.>.>.>>>>h>n>>>(>>p>>n>,>k>>>)>>>>.>.>.>>>>h>n>>>(>>p>>n>,>m>>>)>>>>>>.>>s>\end{array}$

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

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

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

确定分割对象的轮廓被简化为通过选择每行的一个元素在矩阵C中彻底地寻找路径。将第i行中选定元素的指数写成k_i，则曲率变为点序列(p_1，k1，p_2，k2，…，p_n，kn)。最佳路径(k₁^*，k₂^*，…，k_n^*)是指使成本函数J(k₁，…，k_n)最小的路径：

$>>>(sup>>k>1>*sup>>,sup>>k>2>*sup>>,>.>.>.>,sup>>k>n>*sup>>)>>=>arg>>min>>>k>1>>,>.>.>.>,>>k>n>>>>J>>(>>k>1>>,>.>.>.>,>>k>n>>)>>.>>s>$

上面提出的模型允许进行大量的成本测量。考虑两个连续标记p_i，ki和p_i+1，ki+1，根据灰度级外观模型的成本分量以及形状模型的成本分量为：

-与标记p_i，ki及p_i+1，ki+1相对应的灰度级成本h_i(p_i，ki)和h_i+1(p_i+1，ki+1)；

-验证从p_i，ki到p_i+1，ki+1连接的真实性的形状成本f_i(p_i+1，ki+1-p_i，ki)。

通过对总体灰度级成本和总体形状成本的加权求和，将这两种成本类型都结合到成本函数J(k₁，…，k_n)内。加权因子为γ₁的总体灰度级成本是标记的单个灰度级成本之和h_i(p_i，ki)，i＝1，…，n。加权因子为γ₂的总体形状成本是各形状成本之和f_i(p_i+1，ki+1-p_i，ki)，i＝1，…，n。因此成本函数J(k₁，…，k_n)变为：

$>>J>>(>>k>1>>,>.>.>.>,>>k>n>>)>>=>>\gamma >1>>>\Sigma >>i>=>1>>n>>>h>i>>>(>>p>>i>,>>k>i>>>>)>>+>>\gamma >2>>>\Sigma >>i>=>1>>>n>->1>>>>f>i>>>(>>p>>i>+>1>,>>k>>1>+>i>>>>>->>p>>i>,>>k>1>>>>)>>+>>\gamma >2>>>f>n>>>(>>p>>1>,>>k>1>>>>->>p>>n>,>>k>n>>>>)>>.>>s>$

使J(k₁，…，k_n)最小的最佳路径(k₁^*，k₂^*，…，k_n^*)称为最小成本路径。使用例如由G.Behiels等在Evaluation of image feature and searchstrategies for segmentation of bone structures using activeshape models，Medical Image Analysis，6(1)：47-62，2002中提出的现有技术中已知的动态编程技术，可以计算该最佳路径。灰度级成本的典型加权因子为γ₁＝1，形状成本的典型加权因子为γ₂＝0.25。

图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)收集线性或圆形轮廓在标记p₁的灰度值分布，即计算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和协方差S_vi。这个步骤产生统计形状知识。

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

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

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

$>>CTR>=>>>MLD>+>MRD>>ID>>>s>$

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 automatically mapping ofgeometric objects in digital medical images”中公开方法的搜寻栅格定位点映射的概念可以和本发明结合使用。这种改进的方法于是对于分割整形外科放射照片中多骨的身体部分是特别让人感兴趣的，因为其特征在于在图像中具有相对固定的形状，这些部分可以被定位到图像中非常明显的多个标记。

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

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

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