CN101515365A - 一种ct影像中的粘连血管型肺结节自动分割方法 - Google Patents

Abstract

本发明一种CT影像中的粘连血管型肺结节自动分割方法，包括以下步骤：(1)输入包含粘连血管型肺结节的CT影像的感兴趣体；(2)在上述感兴趣体上进行预处理，得到感兴趣体的前景区；(3)在上述前景区上提取基于关系矩阵的流向特征；(4)基于上述流向特征方向角度建立粘连血管型肺结节模型；(5)在上述模型中基于期望最大方法估计模型参数；(6)利用上述模型参数计算得到均值漂移带宽参数；(7)将所述均值漂移带宽参数带入均值漂移聚类算法进行自动分割。应用本发明方法可实现快速、准确地求得自适应带宽参数，同时将自适应带宽参数选择方法用于CT影像中的粘连血管型肺结节自动分割问题中，满足了实际应用中对速度、精度的双重需求。

Description

一种CT影像中的粘连血管型肺结节自动分割方法

技术领域

本发明涉及医学技术领域中的图像处理方法，具体的说是一种CT影像中的粘连血管型肺结节自动分割方法。

背景技术

肺癌是所有癌症中死亡率最高的癌症。在临床中，诊断肺癌的良恶性意义重大。肺结节是肺癌的影像表现形式。肺结节的生长速度是鉴别良恶性的指标。计算机辅助诊断(Computer aided detection，CAD)系统为定量地诊断肺结节良恶性提供了新的手段。它应能自动地帮助医生分割影像上已检测出的结节，目的是测量它的体积和计算一段时间内它的倍增率。倍增率指得是对比同一病人不同时间的两幅影像中结节的尺寸，如果其中结节体积有明显增大，则该结节为恶性结节。血管粘连型结节恶性概率最大，而这类结节的CT像素灰度值通常与血管的CT像素灰度值接近，仅考虑图像CT灰度值将无法准确分割。因此，精确地分割与血管粘连型肺结节是有意义的，并且具有很大难度。

图1为一幅肺CT影像中所示的二维黑色圈部分为粘连血管型结节的一层数据。图2(a)～2(d)所示的三维黑色圈部分为肺CT影像中的一个粘连血管型结节。

聚类分割方法中通常采用均值漂移(Mean-shift)进行结节分割，其主要流程如图3所示：首先输入合适的VOI(Volume of Interesting，感兴趣体)数据，然后自适应地确定带宽参数，再采用均值漂移(Mean-shift)进行结节的聚类分割后输出分割结果。然而，以往的Mean-shift算法需要采用循环迭代的方法自动并且自适应地确定带宽参数，其问题在于：带宽参数与迭代步长有关，步长如果太大，带宽参数就不准确，不准确的带宽参数会导致过分割和欠分割现象，步长如果太小，虽然能得到适当的带宽参数，但由于需要多次循环迭代会消耗大量的时间，从而很难在结节分割的速度、精度方面满足应用的要求。

发明内容

针对现有技术中采用的结节分割方法存在耗时大、速度及精度不能满足应用要求等问题，本发明要解决的技术问题是提供一种满足速度、精度的双重需求的CT影像中的粘连血管型肺结节自动分割方法。

为解决上述技术问题，本发明采用的技术方案是：

本发明一种CT影像中的粘连血管型肺结节自动分割方法，包括以下步骤：

(1)输入包含粘连血管型肺结节的CT影像的感兴趣区域；

(2)在上述感兴趣区域上进行预处理，得到感兴趣区域的前景区；

(3)在上述前景区上提取基于关系矩阵的流向特征；

(4)基于上述流向特征方向角度建立粘连血管型肺结节模型；

(5)在上述模型中基于期望最大方法估计模型参数；

(6)利用上述模型参数计算得到均值漂移带宽参数；

(7)将所述均值漂移带宽参数带入二维均值漂移聚类算法中进行自动分割。

所述预处理采用基于迭代阈值的区域生长方法。

所述在前景区上提取基于关系矩阵的流向特征包括：

(31)在前景区上定义像素j的关系矩阵；

(32)利用上述关系矩阵求出粘连血管型结节像素梯度的法向量即前景区上各个像素的流向特征，并对各个像素的流向特征进行一致化处理；

(33)根据一致化后的各个像素的流向特征画出流向特征向量方向角度直方图。

所述对各个像素的流向特征进行一致化处理包括以下步骤：

从指向同一个方向最多的向量中任取一个向量Arbv；

如果向量Arbv和其它向量v_i之间的夹角在(

π]范围内，把向量v_i转到相反的方向-v_i，即

$v_i=\left\{\begin{array}{cc}-v_i& -1\&le;\frac{v_i\&CenterDot;\mathrm{Arbv}}{\left|v_i\right|*\left|\mathrm{Arbv}\right|}<0\\ v_i& \mathrm{otherwise}\end{array}\right..$

基于流向特征方向角度建立粘连血管型肺结节模型为：

$f\left(\mathrm{\&theta;}\right)=w_1*\frac{1}{\sqrt{2*\mathrm{\&pi;}}}\exp (-\frac{{(\mathrm{\&theta;}-u)}^2}{2*{\mathrm{\&sigma;}}^2})+\frac{w_2}{\mathrm{\&pi;}}$

其中：w₁和w₂分别为正态分布的权重和均匀分布的权重，且有w₁+w₂＝1，μ和σ分别为正态分布的均值和方差。

所述利用模型参数计算得到均值漂移带宽参包括以下步骤：

(61)带宽选择：假设X＝Xi服从正态分布N(μ，∑)，采用正态核函数进行均值漂移，当带宽H＝∑，则均值漂移向量m(x)的模取最大值；

(62)依据带宽选择定理计算均值漂移带宽参数；

$\left\{\begin{array}{c}{h}_x=\mathrm{SXY}/4;\\ h_y=\mathrm{SXY}/4;\\ h_r=\mathrm{SR}/4\end{array}\right.$

其中：h_x、h_y以及h_r分别为空间位置特征x，y和灰度特征r方向上的带宽参数； $\mathrm{SXY}=\mathrm{Round}\left(\sqrt{N_{\mathrm{wu}}*w_2}\right),$ 为结节在x，y方向上的跨度，N_wu为二维平面上前景区的总像素个数，粘连血管型结节模型中结节部分的权重w₂，Round为取整运算；SR＝CT_max-CT_min，为结节在r方向的跨度，CT_max，CT_min分别为N_wu中像素的最大和最小CT值。

本发明一种CT影像中的粘连血管型肺结节自动分割方法，还可以通过以下步骤获得：

(1)输入包含粘连血管型肺结节的CT影像的感兴趣体；

(2)对上述感兴趣体进行预处理，得到感兴趣体的前景区，将该前景区利用最大密度投影方法投影到二维平面上；

(3)在上述二维平面上的前景区投影中提取基于关系矩阵的流向特征；

(4)基于上述流向特征方向角度建立粘连血管型肺结节模型；

(5)在上述模型中基于期望最大方法估计模型参数；

(6)利用上述模型参数计算得到均值漂移带宽参数；

(7)将所述均值漂移带宽参数带入三维均值漂移聚类算法进行自动分割。

所述预处理采用基于迭代阈值的区域生长方法。

所述在二维平面上的前景区投影中提取基于关系矩阵的流向特征包括：

(31)在二维平面上的前景区投影中定义像素j的关系矩阵；

(32)利用上述关系矩阵求出粘连血管型结节像素梯度的法向量即为二维平面上前景区投影中各个像素的流向特征，并对各个像素的流向特征进行一致化处理；

(33)根据一致化后的各个像素的流向特征画出流向特征向量方向角度直方图。

所述对各个像素的流向特征进行一致化处理包括以下步骤：

从指向同一个方向最多的向量中任取一个向量Arbv；

如果向量Arbv和其它向量v_i之间的夹角在(

π]范围内，把向量v_i转到相反的方向-v_i，即

$v_i=\left\{\begin{array}{cc}-v_i& -1\&le;\frac{v_i\&CenterDot;\mathrm{Arbv}}{\left|v_i\right|*\left|\mathrm{Arbv}\right|}<0\\ v_i& \mathrm{otherwise}\end{array}\right.$

所述基于流向特征方向角度建立粘连血管型肺结节模型为：

$f\left(\mathrm{\&theta;}\right)=w_1*\frac{1}{\sqrt{2*\mathrm{\&pi;}}}\exp (-\frac{{(\mathrm{\&theta;}-u)}^2}{2*{\mathrm{\&sigma;}}^2})+\frac{w_2}{\mathrm{\&pi;}}$

其中：w₁和w₂分别为正态分布的权重和均匀分布的权重，且有w₁+w₂＝1，μ和σ分别为正态分布的均值和方差。

所述计算均值漂移带宽参数包括以下步骤：

(61)带宽选择：假设X＝Xi服从正态分布N(μ，∑)，采用正态核函数进行均值漂移，当带宽H＝∑，则均值漂移向量m(x)的模取最大值；

(62)依据带宽选择定理计算均值漂移带宽参数；

$\left\{\begin{array}{c}{h}_x=\mathrm{SXY}/4;\\ h_y=\mathrm{SXY}/4;\\ h_z=\mathrm{SZ}/4;\\ h_r=\mathrm{SR}/4\end{array}\right.$

其中：h_x、h_y、h_z以及h_r分别为空间位置特征x，y，z和灰度特征r方向上的带宽参数； $\mathrm{SXY}=\mathrm{Round}\left(\sqrt{N_{\mathrm{wu}}*w_2}\right),$ 为结节在x，y方向上的跨度，N_wu为二维平面上前景区的总像素个数，粘连血管型结节模型中结节部分的权重w₂，Round为取整运算；SZ＝Round(SXY*xmm/zmm)为结节在z方向上的跨度，xmm为CT影像在x、y方向的分辨率，zmm为CT影像在z方向的分辨率；SR＝CT_max-CT_min，为结节在r方向的跨度，CT_max，CT_min分别为N_wu中像素的最大和最小CT值。

本发明具有以下有益效果及优点：

1.应用本发明方法对肺CT影像中的粘连血管型肺结节自动分割，可以快速自适应地求得Mean-shift(均值漂移)带宽参数，同时将自适应带宽参数选择方法用于CT影像中的粘连血管型肺结节自动分割问题中，对于已检测出来的肺结节可实现快速、准确地分割，分割速度达200层左右的影像数据需要时间少于30秒，很好地满足了实际应用中对速度、精度的双重需求。

附图说明

图1为一层肺CT影像中的粘连血管型结节图像；

图2(a)～(d)分别为三维肺CT影像中的粘连血管型结节在不同扫描层的图像；

图3为现有技术中均值漂移分割方法流程图；

图4为本发明自适应求出均值漂移带宽参数流程图；

图5为图2经过区域生长、迭代阈值和MIP处理后的图像；

图6为图5中各个像素的流向特征图；

图7为图6一致化后的各个像素的流向特征；

图8为图6一致化后的各个像素的流向特征方向角度直方图；

图9(a)为血管粘连型结节人造数据MIP图像；

图9(b)为图9(a)中血管粘连型结节人造数据流向特征方向角度直方图；

图10(a)～(d)分别为采用区域生长方法分割图2(a)～2(d)所示结节在不同层的结果图片(一)～(四)；

图11(a)～(d)分别为采用本发明方法分割图2(a)～2(d)所示结节在不同层的结果图片(一)～(四)；

图12(a)～(d)分别为采用本发明方法分割的血管粘连型结节的3D可视化的效果图(一)～(四)。

具体实施方式

实施例1

本实施例以在二维CT影像中粘连血管型肺结节的自动分割为例。

本发明一种CT影像中的粘连血管型肺结节自动分割方法，可以快速自适应地求得Mean-shift(均值漂移)带宽参数，并用于CT影像中的粘连血管型肺结节自动分割问题中，从而同时满足了血管粘连型肺结节分割算法在速度和精度上的要求。如图4所示，本发明方法的具体步骤如下：

(1)输入包含粘连血管型肺结节的CT影像的感兴趣区域；

(2)在上述感兴趣区域上进行预处理，得到感兴趣区域的前景区；

(3)在上述前景区上提取基于关系矩阵的流向特征；

(4)基于上述流向特征方向角度建立粘连血管型肺结节模型；

(5)在上述模型中基于期望最大方法估计模型参数；

(6)利用上述模型参数计算得到均值漂移带宽参数；

(7)将所述均值漂移带宽参数带入二维均值漂移聚类算法进行自动分割。

上述步骤(1)中，二维感兴趣区域的确定包括以下步骤：

(11)确定肺CT影像结节候选点集S：用基于Hessian矩阵的圆点滤波器卷积肺CT影像，滤波后图像灰度值大于τ的像素为结节候选点，τ为滤波后图像像素灰度值的方差，在候选点集S中用8-领域连通区域生长算法得到结节候选点的中心点c；S中有多个中心点；

(12)选择结节候选点集S中的感兴趣区域：以中心点c作为感兴趣区域的中心，感兴趣区域的尺寸从初始的3*3逐步扩大，直到感兴趣区域中候选点的个数Nh占总区域像素个数Nz的比率R_c ^k小于给定阈值R，R是0～0.5之间小数且与Nh成正比例关系；

上述步骤(2)对感兴趣区域进行预处理是基于区域生长和迭代阈值方法得到感兴趣区域中(Region of Interesting，ROI)的前景区即粘连血管型结节；

在二维前景区上提取基于关系矩阵的流向特征采用以下方法实现：

(31)在前景区上定义像素j的关系矩阵：

本实施例定义像素j的关系矩阵为GG^T，其中 $G=\frac{1}{\sqrt{n}}{\mathrm{\&zeta;}}_j,$ 以像素j为中心的3*3窗中的二维高斯梯度向量ζ_j＝[g₁...g_n]，g_i为像素i的梯度向量，n＝3*3。

(32)利用上述关系矩阵求出粘连血管型结节像素梯度的法向量即为前景区上各个像素的流向特征，并对各个像素的流向特征进行一致化处理：

设λ₁，λ₂为GG^T的特征值且λ₁＜λ₂，与λ₁对应的特征向量是e₁，则像素i的流向特征向量为：

$v_i=\sqrt{{\mathrm{\&lambda;}}_2}*e_1---\left(1\right)$

向量v_i的方向表示像素i上的流向特征方向；

前景区上各个像素的流向特征如图6，流向特征方向角度为[0，2*π]。对各个像素的流向特征进行一致化处理是为了保证血管上的像素的流向特征方向一致，需要把一些向量的方向转换到其相反的方向，从指向同一个方向最多的向量中任取一个向量Arbv，如果它和其它向量v_i之间的夹角在(

π]范围内，把向量v_i转到相反的方向-v_i，即

$v_i=\left\{\begin{array}{cc}-v_i& -1\&le;\frac{v_i\&CenterDot;\mathrm{Avbv}}{\left|v_i\right|*\left|\mathrm{Arbv}\right|}<0\\ v_i& \mathrm{otherwise}\end{array}\right.---\left(2\right)$

如图7所示，为对图6一致化后的各个像素的流向特征

(33)根据一致化后的各个像素的流向特征画出流向特征向量方向角度直方图：

如果直方图不连续，把部分角度加2π后使得直方图连续，并使所有角度在[0，π]范围内。一致化后的各个像素的流向特征方向角度直方图如图8所示。

上述步骤(4)基于上述流向特征方向角度建立粘连血管型肺结节模型为：

$f\left(\mathrm{\&theta;}\right)=w_1*\frac{1}{\sqrt{2*\mathrm{\&pi;}}}\exp (-\frac{{(\mathrm{\&theta;}-u)}^2}{2*{\mathrm{\&sigma;}}^2})+\frac{w_2}{\mathrm{\&pi;}}---\left(3\right)$

其中：w₁和w₂分别为正态分布的权重和均匀分布的权重，且有w₁+w₂＝1，μ和σ分别为正态分布的均值和方差；θ_i∈Θ，i＝1...N_wu，N_wu为集合中总个数，Θ为前景区中流向特征的角度集合。

为简化研究问题，在建立上述模型时，血管通常被看作是符合高斯分布的管状物，而结节通常被看作是符合高斯分布的球状物，血管里充满着向同一方向流动的血液，这就使得血管上像素梯度的垂直方向基本一致，服从正态分布，而结节上像素梯度的垂直方向指向四面八方，服从均匀分布。

上述步骤(5)在上述模型中基于期望最大(Expectation Maximum，EM)方法估计模型参数通过以下步骤实现：

(51)用K均值(K-Means)聚类算法对粘连血管型结节各个像素流向特征方向角度集合聚类(设初始聚类数C＝2)，得到的两类中类个数较小的均值和方差为EM算法中正态分布的初始参数，两类中个数与总像素个数的比值为权重的初始参数；

(52)对于粘连血管型结节各个像素流向特征方向角度集合数据集X的对数似然函数L(Θ|X)，找到使L(Θ|X)最大的Θ：

$L\left(\mathrm{\&Theta;}\right|X)=\log \underset{i=1}{\overset{N_{\mathrm{wu}}}{\mathrm{\&Pi;}}}p\left(x_i\right|\mathrm{\&Theta;})$

$=\underset{i=1}{\overset{N_{\mathrm{wu}}}{\mathrm{\&Sigma;}}}\log (w_1*\frac{1}{\sqrt{2*\mathrm{\&pi;}}}\exp (-\frac{{(x_i-u_1)}^2}{2*{{\mathrm{\&sigma;}}_1}^2})+\frac{w_2}{\mathrm{\&pi;}})---\left(4\right)$

其中：x_i∈X，i＝1...N_wu，N_wu为集合中总个数，Θ为参数，w₁为正态分布的权重，w₂为均匀分布的权重，μ_１，σ₁分别为正态分布的均值和方差；

(53)通过以下公式计算后验期望概率

$p(l=1|x_i,{\mathrm{\&Theta;}}^{\left(s\right)})=\frac{w_1^{\left(s\right)}*\frac{1}{\sqrt{2*\mathrm{\&pi;}}}\exp (-\frac{{(x_i-u_1^{\left(s\right)})}^2}{2*{\mathrm{\&sigma;}}_1^{\left(s\right)2}})}{w_1^{\left(s\right)}*\frac{1}{\sqrt{2*\mathrm{\&pi;}}}\exp (-\frac{{(x_i-u_1^{\left(s\right)})}^2}{2*{\mathrm{\&sigma;}}_1^{\left(s\right)2}})+w_2^{\left(s\right)}*\frac{1}{\mathrm{\&pi;}}}---\left(5\right)$

$p(l=2|x_i,{\mathrm{\&Theta;}}^{\left(s\right)})=\frac{w_2^{\left(s\right)}*\frac{1}{\mathrm{\&pi;}}}{w_1^{\left(s\right)}*\frac{1}{\sqrt{2*\mathrm{\&pi;}}}\exp (-\frac{{(x_i-u_1^{\left(s\right)})}^2}{2*{\mathrm{\&sigma;}}_1^{\left(s\right)2}})+w_2^{\left(s\right)}*\frac{1}{\mathrm{\&pi;}}}---\left(6\right)$

其中：p(l＝1|x_i，Θ^(s))为样本x_i在取得Θ^(s)参数时，样本x_i属于第一类即样本x_i属于正态分布的概率；Θ^(s)是参数Θ在第s次迭代的估计值；x_i∈X，i＝1...N_wu，N_wu为集合中总个数，Θ为参数，w₁ ^(s)，w₂ ^(s)分别为第s次迭代权重参数，μ₁ ^(s)，σ₁ ^(s)分别为第s次迭代正态分布的均值和方差；p(l＝2|x_i，Θ^(s))为样本x_i在取得Θ^(s)参数时，样本x_i属于第二类即样本x_i属于均匀分布的概率；

(54)通过以下公式对各参数进行估计

$w_1^{(s+1)}=\frac{1}{N_{\mathrm{wu}}}\underset{i=1}{\overset{N_{\mathrm{wu}}}{\mathrm{\&Sigma;}}}p(l=1|x_i,{\mathrm{\&Theta;}}^{\left(s\right)})---\left(7\right)$

$w_2^{(s+1)}=\frac{1}{N_{\mathrm{wu}}}\underset{i=1}{\overset{N_{\mathrm{wu}}}{\mathrm{\&Sigma;}}}p(l=2|x_i,{\mathrm{\&Theta;}}^{\left(s\right)})---\left(8\right)$

${\mathrm{\&mu;}}_1^{(s+1)}=\frac{\underset{i=1}{\overset{N_{\mathrm{wu}}}{\mathrm{\&Sigma;}}}{x}_i*p(l=1|x_i,{\mathrm{\&Theta;}}^{\left(s\right)})}{N_{\mathrm{wu}}*w_1^{(s+1)}}---\left(9\right)$

${\mathrm{\&sigma;}}_1^{(s+1)}=\mathrm{Sqrt}\left[\frac{\underset{i=1}{\overset{N_{\mathrm{wu}}}{\mathrm{\&Sigma;}}}p(l=1|x_i,{\mathrm{\&Theta;}}^{\left(s\right)})*(x_i-{\mathrm{\&mu;}}_1^{(s+1)})*(x_i-{\mathrm{\&mu;}}_1^{(s+1)})}{N_{\mathrm{wu}}*w_1^{(s+1)}}\right]---\left(10\right)$

其中，w₁ ^(s+1)，w₂ ^(s+1)分别为第s+1次迭代权重参数，μ₁ ^(s+1)，σ₁ ^(s+1)分别为第s+1次迭代正态分布的均值和方差。

(55)当本次迭代和前次迭代对数似然估计值相同时，所得到的参数即为所求，即满足以下公式：

L^(s+1)(Θ|X)-L^(s)(Θ|X)＜ε

其中：ε是一个接近于0的数，0.001＜ε＜0.01，L^(s+1)为第s+1次迭代的对数似然函数值，x_i∈X，i＝1...N_wu，N_wu为集合中总个数，X为前景区中流向特征的角度集合；Θ为参数，包含正态分布的权重w₁，均匀分布的权重w₂，正态分布的均值和方差μ₁，σ₁。

(56)根据估计出的均匀分布的权重w₂，确定结节像素个数占整个粘连血管型结节像素个数的比例。

所述步骤(6)计算均值漂移带宽参数包括以下步骤：

(61)带宽选择定理：假设样本服从正态分布N(μ，∑)，采用正态核函数进行均值漂移，当带宽H＝∑，则均值漂移向量m(x)的模取最大值；

(62)依据带宽选择定理计算均值漂移带宽参数：

$\left\{\begin{array}{ccc}{h}_x& =& \mathrm{SXY}/4;\\ h_y& =& \mathrm{SXY}/4;\\ h_r& =& \mathrm{SR}/4;\end{array}\right.$

其中：h_x、h_y以及h_r分别为空间位置特征x，y和灰度特征r方向上的带宽参数； $\mathrm{SXY}=\mathrm{Round}\left(\sqrt{N_{\mathrm{wu}}*w_2}\right),$ 为结节在x，y方向上的跨度，N_wu为二维平面上前景区的总像素个数，粘连血管型结节模型中结节部分的权重w₂，Round为取整运算；SR＝CT_max-CT_min，为结节在r方向的跨度，CT_max，CT_min分别为N_wu中像素的最大和最小CT值。

综上所述，本发明方法利用关系矩阵求得粘连血管型结节像素梯度的法向量方向即为粘连血管型结节像素流向特征，假设粘连血管型结节流向特征向量的方向服从正态分布和均匀分布的混合概率分布，基于期望最大(Estimate Maximum，EM)拟合出正态分布均值、方差及权重。均匀分布的权重可以被看作是粘连血管型结节中的结节部分的像素个数占粘连血管型结节总像素个数的比重，再根据均匀分布的权重求得Mean-shift核函数的带宽参数，代入二维Mean-shift聚类算法分割出精确的肺结节边缘。二维Mean-shift的聚类过程如下：

对于ROI(Region ofInteresting，感兴趣区域)的前景区的二维粘连血管型结节数据集在空间位置x，y和灰度r这三维上，X_i∈R³，i＝1，...，N，则样本点X处的多变量高斯核概率密度估计函数为：N为前景区的像素个数，H为带宽参数矩阵，D表示在四维特征空间上样本X和X_i间的距离。

$\hat{f}\left(X\right)=\frac{1}{{N\left(2\mathrm{\&pi;}\right)}^2{\left|H\right|}^{1/2}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\exp (-{\frac{1}{2}D}^2(X,X_i,H));---\left(31\right)$

D²(X，X_i，H)＝(X-X_i)^T*H^-1(X-X_i)；                    (32)

H＝diag(h_x，h_y，h_r)                                  (33)

特征空间中，密度的极大值点为聚类中心点。也就是，密度函数梯度等于零的点，

$\&dtri;\hat{f}\left(X\right)=\frac{1}{{N\left(2\mathrm{\&pi;}\right)}^2}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{H^{-1}(X_i-X)}{{\left|H\right|}^{-1/2}}\exp (-{\frac{1}{2}D}^2(X,X_i,H))---\left(34\right),$

$m\left(X\right)=\frac{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{X}_i\exp \left({-\frac{1}{2}D}^2(X,X_i,H)\right)}{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\exp \left({-\frac{1}{2}D}^2(X,X_i,H)\right)}$ 是样本X漂移到其聚类中心过程中的偏移向量。

Mean-shift向量总是指向密度增长最快的方向。Mean-shift过程对于每个样本点X通过迭代执行下面a，b两步的操作，保证了收敛到聚类中心：

a.计算偏移向量m(X)。

b.平移窗口，窗口中心位置移到X^t+1，X^t+1＝X^t+m(X^t)。

结束条件当X^t+1＝X^t，表示该样本聚类过程收敛，停止迭代。X^t是样本X∈R³在第t次迭代中的坐标和CT灰度值，X^t+1是样本X在第t+1次迭代中的坐标和CT灰度值。

对于每个数据点，运行上述过程，收敛到相同聚类中心的点的集合被认为属于同一个聚类簇。

实施例2

本实施例以在三维CT影像中粘连血管型肺结节的自动分割为例。

本发明方法的具体步骤如下：

(1)输入包含粘连血管型肺结节的CT影像的感兴趣体；

(2)对上述感兴趣体进行预处理，得到感兴趣体的前景区，将该前景区利用最大密度投影方法投影到二维平面上；

(3)在上述二维平面上的前景区投影中提取基于关系矩阵的流向特征；

(4)基于上述流向特征方向角度建立粘连血管型肺结节模型；

(5)在上述模型中基于期望最大方法估计模型参数；

(6)利用上述模型参数计算得到均值漂移带宽参数；

(7)将所述均值漂移带宽参数带入三维均值漂移聚类算法进行自动分割。

上述步骤(1)中，三维感兴趣体的确定包括以下步骤：

(11)确定肺CT影像结节候选点集S：用基于Hessian矩阵的圆点滤波器卷积肺CT影像，滤波后图像灰度值大于τ的为结节候选点，τ为滤波后图像像素灰度值的方差，在候选点集S中用26领域连通区域生长算法得到结节候选点的中心点c；S中有多个中心点。

(12)选择结节候选点集S中的感兴趣体：感兴趣体以中心点c作为感兴趣体的中心，感兴趣体的尺寸从初始的3*3*3逐步扩大，直到感兴趣体中候选点的个数Nh占总感兴趣体像素个数Nz的比率R_c ^k小于给定阈值R，R是0～0.5之间小数且与Nh成正比例关系。

步骤(2)是采用以下方法进行的：

(21)基于区域生长和迭代阈值方法得到感兴趣体中(Volume ofInteresting，VOI)的前景区即粘连血管型结节；

(22)将感兴趣体的前景区利用最大密度方法(Maximum IntensityProject，MIP)投影到二维平面上，效果如图5所示。

上述步骤(3)在上述二维平面上的前景区投影中提取基于关系矩阵的流向特征采用以下方法实现：

(31)在二维平面的前景区上定义像素j的关系矩阵：

本实施例定义像素j的关系矩阵为GG^T，其中 $G=\frac{1}{\sqrt{n}}{\mathrm{\&zeta;}}_j,$ 以像素j为中心的3*3窗中的二维高斯梯度向量ζ_j＝[g₁...g_n]，g_i为像素i的梯度向量，n＝3*3。

(32)利用上述关系矩阵求出粘连血管型结节像素梯度的法向量即为二维平面的前景区上各个像素的流向特征，并对各个像素的流向特征进行一致化处理：

设λ₁，λ₂为GG^T的特征值且λ₁＜λ₂，与λ₁对应的特征向量是e₁，则像素i的流向特征向量为：

$v_i=\sqrt{{\mathrm{\&lambda;}}_2}*e_1---\left(1\right)$

向量v_i的方向表示像素i上的流向特征方向；

二维平面的前景区上各个像素的流向特征如图6，流向特征方向角度为[0，2*π]。对各个像素的流向特征进行一致化处理是为了保证血管上的像素的流向特征一致，需要把一些向量的方向转换到其相反的方向，从指向同一个方向最多的向量中任取一个向量Arbv，如果它和其它向量v_i之间的夹角在(

π]范围内，把向量v_i转到相反的方向-v_i，即

$v_i=\left\{\begin{array}{cc}-v_i& -1\&le;\frac{v_i\&CenterDot;\mathrm{Arbv}}{\left|v_i\right|*\left|\mathrm{Arbv}\right|}<0\\ v_i& \mathrm{otherwise}\end{array}\right.---\left(2\right)$

如图7所示，为对图6一致化后的各个像素的流向特征。

(33)根据一致化后的各个像素的流向特征画出流向特征向量方向角度直方图：

如果直方图不连续，把部分角度加2π后使得直方图连续，并使所有角度在[0，π]范围内。一致化后的各个像素的流向特征方向角度直方图如图8所示。

本实施例中，步骤(4)及步骤(5)与实施例1的相应步骤一致，此处不再重述。

所述步骤(6)计算均值漂移带宽参数包括以下步骤：

(61)带宽选择定理：假设样本服从正态分布N(μ，∑)，采用正态核函数进行均值漂移，当带宽H＝∑，则均值漂移向量m(x)的模取最大值；

(62)依据带宽选择定理计算均值漂移带宽参数：

$\left\{\begin{array}{c}{h}_x=\mathrm{SXY}/4;\\ h_y=\mathrm{SXY}/4;\\ h_z=\mathrm{SZ}/4;\\ h_r=\mathrm{SR}/4\end{array}\right.$

其中：h_x、h_y、h_z以及h_r分别为空间位置特征x，y，z和灰度特征r方向上的带宽参数； $\mathrm{SXY}=\mathrm{Round}\left(\sqrt{N_{\mathrm{wu}}*w_2}\right),$ 为结节在x，y方向上的跨度，N_wu为二维平面上前景区的总像素个数，粘连血管型结节模型中结节部分的权重w₂，Round为取整运算；SZ＝Round(SXY*xmm/zmm)为结节在z方向上的跨度，xmm为CT影像在x、y方向的分辨率，zmm为CT影像在z方向的分辨率；SR＝CT_max-CT_min，为结节在r方向的跨度，CT_max，CT_min分别为N_wu中像素的最大和最小CT值。

综上所述，本发明方法利用关系矩阵求得粘连血管型结节像素梯度的法向量方向即为粘连血管型结节像素流向特征，假设粘连血管型结节流向特征向量的方向服从正态分布和均匀分布的混合概率分布，基于期望最大(Estimate Maximum，EM)拟合出正态分布均值、方差及权重。均匀分布的权重可以被看作是粘连血管型结节中的结节部分的像素个数占粘连血管型结节总像素个数的比重，再根据均匀分布的权重求得Mean-shift核函数的带宽参数，代入三维Mean-shift聚类算法分割出精确的肺结节边缘。三维Mean-shift的聚类过程如下：

对于VOI前景区(三维的粘连血管型结节)数据集在空间位置x，y，z和灰度r这四维上，X_i∈R⁴，i＝1，...，N，则样本点X处的多变量高斯核概率密度估计函数为：N为前景区的像素个数，H为带宽参数矩阵，D表示在四维特征空间上样本X和X_i间的距离。

$\hat{f}\left(X\right)=\frac{1}{N{\left(2\mathrm{\&pi;}\right)}^2{\left|H\right|}^{1/2}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\exp (-\frac{1}{2}{D}^2(X,X_i,H));---\left(31\right)$

D²(X，X_i，H)＝(X-X_i)^T*H^-1(X-X_i)；                        (32)

H＝diag(h_x，h_y，h_z，h_r)                                  (33)。

特征空间中，密度的极大值点为聚类中心点。也就是，密度函数梯度等于零的点，

$\&dtri;\hat{f}\left(X\right)=\frac{1}{N{\left(2\mathrm{\&pi;}\right)}^2}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{H^{-1}(X_i-X)}{{\left|H\right|}^{-1/2}}\exp (-\frac{1}{2}{D}^2(X,X_i,H))---\left(34\right)$

$m\left(X\right)=\frac{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{X}_i\exp (-\frac{1}{2}{D}^2(X,X_i,H))}{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\exp (-\frac{1}{2}{D}^2(X,X_i,H))}$ 是样本X漂移到其聚类中心过程中的偏移向量。

Mean-shift向量总是指向密度增长最快的方向。Mean-shift过程对于每个样本点X通过迭代的执行下面a，b两步操作，保证了收敛到聚类中心：

c.计算偏移向量m(X)。

d.平移窗口，窗口中心位置移到X^t+1，X^t+1＝X^t+m(X^t)。

结束条件当X^t+1＝X^t，表示该样本聚类过程收敛，停止迭代。X^t是样本X∈R⁴在第t次迭代中的坐标和CT灰度值，X^t+1是样本X在第t+1次迭代中的坐标和CT灰度值。

对于每个数据点，运行上述过程，收敛到相同聚类中心的点的集合被认为属于同一个聚类簇。

通过人造数据和肺结节真实影像数据的实验可以表明本发明方法能够快速而有效地解决粘连血管型结节的分割问题，所述人造数据设计如下：

一个半径为6(单位为像素，以下同)，尺度为1的三维服从高斯分布的球和一个半径为2，高度为21，尺度为1的三维服从高斯分布的圆柱相叠加模拟血管粘连型结节。最大密度投影后如图9(a)所示，流向特征方向角度直方图如图9(b)所示。可以看出人造血管粘连型结节的流向特征方向角度直方图服从正态分布和均匀分布的叠加。图9(a)中的人造数据和19套真实CT影像数据都能得到正确的肺结节分割结果。其中，19套真实CT影像数据分割结果是和医生的金标准(医生手工画出的结节轮廓线)相比较后，得出结论的。

图10(a)～10(d)分别显示采用区域生长方法分割图2(a)～2(d)所示结节在不同层的结果，图11(a)～11(d)分别显示采用本发明方法分割图2(a)～2(d)所示结节在不同层的结果。将图10(a)～10(d)和图11(a)～11(d)比较后，可以清楚的看出，本发明方法可以把结节和血管精确地分割开，而传统的区域生长方法则无法把血管和结节精确分割。

图12(a)～(d)中的四幅图为采用本发明方法分割的血管粘连型结节的三维可视化的效果图。图中深色圆球状物体即为分割出的结节。

应用本发明方法分割CT影像中的粘连血管型肺结节速度快，每分割200左右层的CT影像数据需要时间少于30秒。满足速度的要求体现在不用重复迭代地计算自适应的带宽参数，而是直接求得自适应带宽参数，这样在求得带宽参数的时间上就能缩短，从而分割算法的时间也就相应缩短了。

Claims (12)

1.一种CT影像中的粘连血管型肺结节自动分割方法，其特征在于包括以下步骤：

(1)输入包含粘连血管型肺结节的CT影像的感兴趣区域；

(2)在上述感兴趣区域上进行预处理，得到感兴趣区域的前景区；

(3)在上述前景区上提取基于关系矩阵的流向特征；

(4)基于上述流向特征方向角度建立粘连血管型肺结节模型；

(5)在上述模型中基于期望最大方法估计模型参数；

(6)利用上述模型参数计算得到均值漂移带宽参数；

(7)将所述均值漂移带宽参数带入二维均值漂移聚类算法进行自动分割。

2.按权利要求1所述的CT影像中的粘连血管型肺结节自动分割方法，其特征在于：所述预处理采用基于迭代阈值的区域生长方法。

3.按权利要求1所述的CT影像中的粘连血管型肺结节自动分割方法，其特征在于所述在前景区上提取基于关系矩阵的流向特征包括：

(31)在前景区上定义像素j的关系矩阵；

(32)利用上述关系矩阵求出粘连血管型结像素梯度的法向量即前景区上各个像素的流向特征，并对各个像素的流向特征进行一致化处理；

(33)根据一致化后的各个像素的流向特征画出流向特征向量方向角度直方图。

4.按权利要求3所述的CT影像中的粘连血管型肺结节自动分割方法，其特征在于所述对各个像素的流向特征进行一致化处理包括以下步骤：

从指向同一个方向最多的向量中任取一个向量Arbv；

如果向量Arbv和其它向量v_i之间的夹角在(

π]范围内，把向量v_i转到相反的方向-v_i，即

$v_i=\left\{\begin{array}{cc}-v_i& -1\&le;\frac{v_i\&CenterDot;\mathrm{Arbv}}{\left|v_i\right|*\left|\mathrm{Arbv}\right|}<0\\ v_i& \mathrm{otherwise}\end{array}\right..$

5.按权利要求1所述的CT影像中的粘连血管型肺结节自动分割方法，其特征在于所述基于流向特征方向角度建立粘连血管型肺结节模型为：

$f\left(\mathrm{\&theta;}\right)=w_1*\frac{1}{\sqrt{2*\mathrm{\&pi;}}}\exp (-\frac{{(\mathrm{\&theta;}-u)}^2}{2*{\mathrm{\&sigma;}}^2})+\frac{w_2}{\mathrm{\&pi;}}$

其中：w₁和w₂分别为正态分布的权重和均匀分布的权重，且有w₁+w₂＝1，μ和σ分别为正态分布的均值和方差。

6.按权利要求1所述的CT影像中的粘连血管型肺结节自动分割方法，其特征在于所述利用模型参数计算得到均值漂移带宽参包括以下步骤：

(61)带宽选择：假设X＝Xi服从正态分布N(μ，∑)，采用正态核函数进行均值漂移，当带宽H＝∑，则均值漂移向量m(x)的模取最大值；

(62)依据带宽选择定理计算均值漂移带宽参数；

$\left\{\begin{array}{c}{h}_x=\mathrm{SXY}/4;\\ h_y=\mathrm{SXY}/4;\\ h_r=\mathrm{SR}/4\end{array}\right.$

其中：h_x、h_y以及h_r分别为空间位置特征x，y和灰度特征r方向上的带宽参数； $\mathrm{SXY}=\mathrm{Round}\left(\sqrt{N_{\mathrm{wu}}*w_2}\right),$ 为结节在x，y方向上的跨度，N_wu为二维平面上前景区的总像素个数，粘连血管型结节模型中结节部分的权重w₂，Round为取整运算；SR＝CT_max-CT_min，为结节在r方向的跨度，CT_max，CT_min分别为N_wu中像素的最大和最小CT值。

7.一种CT影像中的粘连血管型肺结节自动分割方法，其特征在于包括以下步骤：

(1)输入包含粘连血管型肺结节的CT影像的感兴趣体；

(2)对上述感兴趣体进行预处理，得到感兴趣体的前景区，将该前景区利用最大密度投影方法投影到二维平面上；

(3)在上述二维平面上的前景区投影中提取基于关系矩阵的流向特征；

(4)基于上述流向特征方向角度建立粘连血管型肺结节模型；

(5)在上述模型中基于期望最大方法估计模型参数；

(6)利用上述模型参数计算得到均值漂移带宽参数；

(7)将所述均值漂移带宽参数带入三维均值漂移聚类算法进行自动分割。

8.按权利要求7所述的CT影像中的粘连血管型肺结节自动分割方法，其特征在于：所述预处理采用基于迭代阈值的区域生长方法。

9.按权利要求7所述的CT影像中的粘连血管型肺结节自动分割方法，其特征在于所述在二维平面上的前景区投影中提取基于关系矩阵的流向特征包括：

(31)在二维平面上的前景区投影中定义像素j的关系矩阵；

(32)利用上述关系矩阵求出粘连血管型结节像素梯度的法向量即为二维平面上前景区投影中各个像素的流向特征，并对各个像素的流向特征进行一致化处理；

(33)根据一致化后的各个像素的流向特征画出流向特征向量方向角度直方图。

10.按权利要求9所述的CT影像中的粘连血管型肺结节自动分割方法，其特征在于：所述对各个像素的流向特征进行一致化处理包括以下步骤：

从指向同一个方向最多的向量中任取一个向量Arbv；

如果向量Arbv和其它向量v_i之间的夹角在(

π]范围内，把向量v_i转到相反的方向-v_i，即

$v_i=\left\{\begin{array}{cc}-v_i& -1\&le;\frac{v_i\&CenterDot;\mathrm{Arbv}}{\left|v_i\right|*\left|\mathrm{Arbv}\right|}<0\\ v_i& \mathrm{otherwise}\end{array}\right.$

11.按权利要求7所述的CT影像中的粘连血管型肺结节自动分割方法，其特征在于：所述基于流向特征方向角度建立粘连血管型肺结节模型为：

$f\left(\mathrm{\&theta;}\right)=w_1*\frac{1}{\sqrt{2*\mathrm{\&pi;}}}\exp (-\frac{{(\mathrm{\&theta;}-u)}^2}{2*{\mathrm{\&sigma;}}^2})+\frac{w_2}{\mathrm{\&pi;}}$

其中：w₁和w₂分别为正态分布的权重和均匀分布的权重，且有w₁+w₂＝1，μ和σ分别为正态分布的均值和方差。

12.按权利要求7所述的CT影像中的粘连血管型肺结节自动分割方法，其特征在于：所述计算均值漂移带宽参数包括以下步骤：

(61)带宽选择：假设X＝Xi服从正态分布N(μ，∑)，采用正态核函数进行均值漂移，当带宽H＝∑，则均值漂移向量m(x)的模取最大值；

(62)依据带宽选择定理计算均值漂移带宽参数；

$\left\{\begin{array}{c}{h}_x=\mathrm{SXY}/4;\\ h_y=\mathrm{SXY}/4;\\ h_z=\mathrm{SZ}/4;\\ h_r=\mathrm{SR}/4\end{array}\right.$

其中：h_x、h_y、h_r以及h_r分别为空间位置特征x，y，z和灰度特征r方向上的带宽参数； $\mathrm{SXY}=\mathrm{Round}\left(\sqrt{N_{\mathrm{wu}}*w_2}\right),$ 为结节在x，y方向上的跨度，N_wu为二维平面上前景区的总像素个数，粘连血管型结节模型中结节部分的权重w₂，Round为取整运算；SZ＝Round(SXY*xmm/zmm)为结节在z方向上的跨度，xmm为CT影像在x、y方向的分辨率，zmm为CT影像在z方向的分辨率；SR＝CT_max-CT_min，为结节在r方向的跨度，CT_max，CT_min分别为N_wu中像素的最大和最小CT值。

Images