CN105551041A - 普适的血管分割方法及系统 - Google Patents

Abstract

本发明属于医疗图像领域，提供一种普适的血管分割方法及系统，所述方法包括如下步骤：对输入的血管造影原始图像数据进行多尺度滤波后获得增强后的图像；对增强后的直方图曲线拟合血管数据和拟合噪声数据得到增强后的图像的混合模型；对混合模型中的参数进行参数估计处理；通过最大后验概率估计法则将经过参数估计处理的混合模型中的血管从图像中分离出来。本发明提供的技术方案具有适应性普遍的优点。

Description

普适的血管分割方法及系统

技术领域

本发明属于医疗图像领域，尤其涉及一种普适的血管分割方法及系统。

背景技术

目前，心脑血管疾病已经严重威胁着人类的健康。随着CT血管造影、磁共振血管造影(MRA)等成像技术的发展，医学图像的后处理技术变得越来越重要,精确的血管分割是血管类疾病诊断与治疗的关键。

多尺度滤波、形变模型、统计模型及混合方法在血管分割中得到了广泛的研究和应用。多尺度滤波通过抑制背景噪声增强血管，但是血管并没有与背景区分开来。形变模型有效地整合了自底向上的信息和自上而下的先验知识，但分割质量主要依赖于模型参数。统计学模型主要是根据像素点的强度信息进行分类，但其模型的选择受限于图像的统计学特性。混合方法是将多种分割算法相结合，以降低血管误分割的概率。

到目前为止，基于统计学的分割算法已经引起了越来越多的关注，而模型选择则是基于统计学的分割算法的关键。针对脑部MRA图像，Wilson和Noble于1999年提出了一种基于统计学模型的血管分割方法。根据灰度值特性，他们用一个均匀分布拟合血管类，两个高斯分布拟合背景类。2006年，Hassouna对模型进行了改进，他的方法中用一个瑞利函数和两个高斯函数拟合背景噪声，同时用一个高斯函数拟合血管，并取得了更好的血管分割效果。2013年，将统计学模型和多模式邻域系统相结合，可以提取出更精细的血管目标。然而，不同模态的血管图像或者不同器官的血管成像具有不同的统计学特性，很难用同样的混合模型去精确拟合。

发明内容

本发明实施例的目的在于提供一种普适的血管分割方法，该方法能够普遍适应血管的精确拟合。

本发明提供一种普适的血管分割方法，所述方法包括如下步骤：

对输入的血管造影原始图像数据进行多尺度滤波后获得增强后的图像；

对增强后的直方图曲线拟合血管数据和拟合噪声数据得到增强后的图像的混合模型；

对混合模型中的参数进行参数估计处理；

通过最大后验概率估计法则将经过参数估计处理的混合模型中的血管从图像中分离出来。

可选的，所述对输入的血管造影原始图像数据进行多尺度滤波后获得增强后的图像具体，包括：

多尺度滤波算法通过尺度变换，利用Hessian矩阵的特征值计算局部血管出现的概率，所述血管增强函数如下：

$v_0\left(s\right)=\left\{\begin{array}{c}0\\ (1-\exp (-\frac{{R_A}^2}{2{\mathrm{\α}}^2}\left)\right)\exp (-\frac{{R_B}^2}{2{\mathrm{\β}}^2})(1-\exp (-\frac{S^2}{2c^2}\left)\right)\end{array}\right.---\left(1\right)$

其中，R_A，R_B，S为三个测度函数，如公式(2)所示，且R_A用来区分片状和线状结构，R_B用来区分点状结构和线状结构，S用于区分背景像素，α,β和c分别可以为个测度的阈值用于控制；

$\left\{\begin{array}{c}{R}_A=\frac{\left|h_2\right|}{\left|h_3\right|}\\ R_B=\frac{\left|h_1\right|}{\sqrt{\left|h_2{h}_3\right|}}\\ S=\left|\right|H|{|}_F=\sqrt{\underset{j\≤D}{\Σ}{h}_j^2}\end{array}\right..$

可选的，所述对增强后的直方图曲线拟合血管数据和拟合噪声数据得到增强后的图像的混合模型具体，包括：

混合模型f(x)如式(3)所示：

$f\left(x\right)=\underset{l=1}{\overset{2}{\Σ}}{w}_{El}{f}_{El}\left(x\right)+w_G{f}_G\left(x\right)---\left(3\right)$

其中，f_El(x)(l＝1,2)和f_G(x)分别表示指数分布和高斯分布，如式(4)所示，w_El(l＝1,2)为指数分布所占的比例，w_G可以为高斯分布所占的比例，且满足 $\underset{l=1}{\overset{2}{\Σ}}{w}_{El}+w_G=1;$

$\left\{\begin{array}{c}{f}_{El}\left(x\right)={\mathrm{\λ}}_{El}\exp (-{\mathrm{\λ}}_{El}x),l\∈\[1,2\]\\ f_G\left(x\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_G}\exp \left(\frac{-{(x-{\mathrm{\μ}}_G)}^2}{2{{\mathrm{\σ}}^2}_G}\right)\end{array}\right.---\left(4\right).$

可选的，所述对混合模型中的参数进行参数估计处理具体，包括：

当原始血管造影图像的峰值点超过3个时，则用3个峰值点完成k均值的初始化完成参数估计处理，否则，采用如下规则完成k均值算法的初始化；

$\left\{\begin{array}{c}{\mathrm{\μ}}_{k\_1}^{init}=I_{\max }/1000\\ {\mathrm{\μ}}_{k\_2}^{init}=I_{\max }/100\\ {\mathrm{\μ}}_{k\_3}^{init}=I_{\max }/10\end{array}\right.---\left(5\right)$

其中，I_max是原始图像中像素最大值，和表示k均值算法的初始聚类中心，μ_{k_l},σ_{k_l},w_{k_l},(l＝1,2,3)是由k均值聚类算法估计得到的三个分布的模型参数，

$\left\{\begin{array}{c}{\mathrm{\λ}}_{E1}^{init}=1.0/{\mathrm{\μ}}_{k\_1},w_{E1}^{init}=w_{k\_1}\\ {\mathrm{\λ}}_{E2}^{init}=1.0/{\mathrm{\μ}}_{k\_2},=w_2^{init}=w_{k\_2}\\ {\mathrm{\μ}}_G^{init}={\mathrm{\μ}}_{k\_3},{\mathrm{\σ}}_G^{init}={\mathrm{\σ}}_{k\_3},w_G^{init}=w_{k\_3}\end{array}\right.---\left(6\right)$

由EM算法得到的迭代更新公式如公式(7)和(8)：

$\left\{\begin{array}{c}{\mathrm{\λ}}_{Ei}^{k+1}=\frac{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(Ei\right|x_j)}{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(Ei\right|x_j)x_j}\\ w_{Ei}^{k+1}=\frac{1}{N}\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(Ei\right|x_j)\end{array}\right.,(i=1,2)---\left(7\right)$

$\left\{\begin{array}{c}{\mathrm{\μ}}_G^{k+1}=\frac{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j)x_j}{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j)},w_G^{k+1}=\frac{1}{N}\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j)\\ {\left({\mathrm{\σ}}_G^{k+1}\right)}^2=\frac{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j){(x_j-{\mathrm{\μ}}_G^{k+1})}^2}{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j)}\end{array}\right.---\left(8\right)$

其中，x_j是第j个体素的强度值，N是体素的总数，后验概率f^k(Ei|x_j)和f^k(G|x_j)由贝叶斯公式计算如下：

$\left\{\begin{array}{c}{f}^k\left(Ei\right|x_j)=\frac{w_{Ei}^k{f}^k\left(x_j\right|Ei)}{\underset{i=1}{\overset{2}{\Σ}}{w}_{Ei}^k{f}^k\left(x_j\right|Ei)+w_G^k{f}^k\left(x_j\right|G)}\\ f^k\left(G\right|x_j)=\frac{w_G^k{f}^k\left(x_j\right|G)}{\underset{i=1}{\overset{2}{\Σ}}{w}_{Ei}^k{f}^k\left(x_j\right|Ei)+w_G^k{f}^k\left(x_j\right|G)}\end{array}\right.---\left(9\right).$

可选的，所述通过最大后验概率估计法则将经过参数估计处理的混合模型中的血管从图像中分离出来具体，包括：

由最大后验概率估计法则，根据公式(10)可以最终把血管从图像中分离出来：

w_Gf_G(x_j)＞w_E1f_E1(x_j)+w_E2f_E2(x_j)(10)。

另一方面，提供一种普适的血管分割系统，所述系统包括：

增强单元，用于对输入的血管造影原始图像数据进行多尺度滤波后获得增强后的图像；

拟合单元，用于对增强后的直方图曲线拟合血管数据和拟合噪声数据得到增强后的图像的混合模型；

参数估计单元，用于对混合模型中的参数进行参数估计处理；

分离单元，用于通过最大后验概率估计法则将经过参数估计处理的混合模型中的血管从图像中分离出来。

可选的，所述增强单元具体，用于多尺度滤波算法通过尺度变换，利用Hessian矩阵的特征值计算局部血管出现的概率，所述血管增强函数如下：

$v_0\left(s\right)=\left\{\begin{array}{c}0\\ (1-\exp (-\frac{{R_A}^2}{2{\mathrm{\α}}^2}\left)\right)\exp (-\frac{{R_B}^2}{2{\mathrm{\β}}^2})(1-\exp (-\frac{S^2}{2c^2}\left)\right)\end{array}\right.---\left(1\right)$

其中，R_A，R_B，S为三个测度函数，如公式(2)所示，且R_A用来区分片状和线状结构，R_B用来区分点状结构和线状结构，S用于区分背景像素，α,β和c分别可以为个测度的阈值用于控制；

$\left\{\begin{array}{c}{R}_A=\frac{\left|h_2\right|}{\left|h_3\right|}\\ R_B=\frac{\left|h_1\right|}{\sqrt{\left|h_2{h}_3\right|}}\\ S=\left|\right|H|{|}_F=\sqrt{\underset{j\≤D}{\Σ}{h}_j^2}\end{array}\right..$

可选的，所述拟合单元具体，用于

混合模型f(x)如式(3)所示：

$f\left(x\right)=\underset{l=1}{\overset{2}{\Σ}}{w}_{El}{f}_{El}\left(x\right)+w_G{f}_G\left(x\right)---\left(3\right)$

其中，f_El(x)(l＝1,2)和f_G(x)分别表示指数分布和高斯分布，如式(4)所示，w_El(l＝1,2)为指数分布所占的比例，w_G可以为高斯分布所占的比例，且满足 $\underset{l=1}{\overset{2}{\Σ}}{w}_{El}+w_G=1;$

$\left\{\begin{array}{c}{f}_{El}\left(x\right)={\mathrm{\λ}}_{El}\exp (-{\mathrm{\λ}}_{El}x),l\∈\[1,2\]\\ f_G\left(x\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_G}\exp \left(\frac{-{(x-{\mathrm{\μ}}_G)}^2}{2{{\mathrm{\σ}}^2}_G}\right)\end{array}\right.---\left(4\right).$

可选的，所述参数估计单元具体，用于

当原始血管造影图像的峰值点超过3个时，则用3个峰值点完成k均值的初始化完成参数估计处理，否则，采用如下规则完成k均值算法的初始化；

$\left\{\begin{array}{c}{\mathrm{\μ}}_{k\_1}^{init}=I_{\max }/1000\\ {\mathrm{\μ}}_{k\_2}^{init}=I_{\max }/100\\ {\mathrm{\μ}}_{k\_3}^{init}=I_{\max }/10\end{array}\right.---\left(5\right)$

其中，I_max是原始图像中像素最大值，和表示k均值算法的初始聚类中心，μ_{k_l},σ_{k_l},w_{k_l},(l＝1,2,3)是由k均值聚类算法估计得到的三个分布的模型参数，

$\left\{\begin{array}{c}{\mathrm{\λ}}_{E1}^{init}=1.0/{\mathrm{\μ}}_{k\_1},w_{E1}^{init}=w_{k\_1}\\ {\mathrm{\λ}}_{E2}^{init}=1.0/{\mathrm{\μ}}_{k\_2},=w_2^{init}=w_{k\_2}\\ {\mathrm{\μ}}_G^{init}={\mathrm{\μ}}_{k\_3},{\mathrm{\σ}}_G^{init}={\mathrm{\σ}}_{k\_3},w_G^{init}=w_{k\_3}\end{array}\right.---\left(6\right)$

由EM算法得到的迭代更新公式如公式(7)和(8)：

$\left\{\begin{array}{c}{\mathrm{\λ}}_{Ei}^{k+1}=\frac{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(Ei\right|x_j)}{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(Ei\right|x_j)x_j}\\ w_{Ei}^{k+1}=\frac{1}{N}\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(Ei\right|x_j)\end{array}\right.,(i=1,2)---\left(7\right)$

$\left\{\begin{array}{c}{\mathrm{\μ}}_G^{k+1}=\frac{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j)x_j}{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j)},w_G^{k+1}=\frac{1}{N}\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j)\\ {\left({\mathrm{\σ}}_G^{k+1}\right)}^2=\frac{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j){(x_j-{\mathrm{\μ}}_G^{k+1})}^2}{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j)}\end{array}\right.---\left(8\right)$

其中，x_j是第j个体素的强度值，N是体素的总数，后验概率f^k(Ei|x_j)和f^k(G|x_j)由贝叶斯公式计算如下：

$\left\{\begin{array}{c}{f}^k\left(Ei\right|x_j)=\frac{w_{Ei}^k{f}^k\left(x_j\right|Ei)}{\underset{i=1}{\overset{2}{\Σ}}{w}_{Ei}^k{f}^k\left(x_j\right|Ei)+w_G^k{f}^k\left(x_j\right|G)}\\ f^k\left(G\right|x_j)=\frac{w_G^k{f}^k\left(x_j\right|G)}{\underset{i=1}{\overset{2}{\Σ}}{w}_{Ei}^k{f}^k\left(x_j\right|Ei)+w_G^k{f}^k\left(x_j\right|G)}\end{array}\right.---\left(9\right).$

可选的，所述分离单元具体，用于

由最大后验概率估计法则，根据公式(10)可以最终把血管从图像中分离出来：

w_Gf_G(x_j)＞w_E1f_E1(x_j)+w_E2f_E2(x_j)(10)。

在本发明实施例中，本发明提供的技术方案将多尺度滤波和统计学方法相结合，提出了一种普适的血管分割方法。首先，不同模态的图像数据通过多尺度滤波进行血管增强，增强后的数据具有一个新的相近的统计学特性，针对滤波后的数据，利用一个高斯函数拟合血管类，同时根据直方图特性选用两个指数函数拟合背景噪声类，模型参数选用最大期望(EM)算法进行估算，最后，本发明利用最大后验概率分类将血管和背景标记出来。

附图说明

图1为本发明第一较佳实施方式提供的一种普适的血管分割方法的流程图；

图2为本发明第一较佳实施方式提供的增强后的直方图曲线示意图；

图3为本发明第一较佳实施方式的仿真数据示意图；

图4为本发明第一较佳实施方式的临床数据示意图；

图5为本发明第一较佳实施方式的血管增强前后图像直方图；

图6为本发明第一较佳实施方式的直方图拟合结果示意图；

图7为本发明第一较佳实施方式的仿真数据分割结果示意图；

图8为本发明第一较佳实施方式的临床数据分割后的MIP图像和三维视觉效果图；

图9为本发明第一较佳实施方式的效果对比示意图；

图10为本发明第二较佳实施方式提供的普适的血管分割系统结构图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明第一较佳实施方式提供一种普适的血管分割方法，该方法由智能设备执行，该智能设备包括但不限于：个人计算机、服务器或医疗设备等。该方法如图1所示，包括如下步骤：

步骤S101、对输入的血管造影原始图像数据进行多尺度滤波后获得增强后的图像；

由于不同组织器官或不同模态的血管造影图像的直方图曲线具有不同的统计学特性，很难找到一个固定的混合模型去拟合不同的直方图曲线。因此，本发明第一较佳实施方式利用血管增强算法进行图像预处理，使得增强后的直方图曲线(如图2所示)具有相对一致的统计学特性。多尺度滤波算法通过尺度变换，利用Hessian矩阵的特征值计算局部血管出现的概率，其血管增强函数如下：

$v_0\left(s\right)=\left\{\begin{array}{c}0\\ (1-\exp (-\frac{{R_A}^2}{2{\mathrm{\α}}^2}\left)\right)\exp (-\frac{{R_B}^2}{2{\mathrm{\β}}^2})(1-\exp (-\frac{S^2}{2c^2}\left)\right)\end{array}\right.---\left(1\right)$

其中，R_A，R_B，S为三个测度函数，如公式(2)所示，且R_A可以用来区分片状和线状结构，R_B可以用来区分点状结构和线状结构，S可以用于区分背景像素。α,β和c分别可以为个测度的阈值用于控制上述血管增强算法对R_A，R_B，S的敏感性，上述α,β和c的具体取值可以由使用者根据实际情况自行设定。

$\left\{\begin{array}{c}{R}_A=\frac{\left|h_2\right|}{\left|h_3\right|}\\ R_B=\frac{\left|h_1\right|}{\sqrt{\left|h_2{h}_3\right|}}\\ S=\left|\right|H|{|}_F=\sqrt{\underset{j\≤D}{\Σ}{h}_j^2}\end{array}\right.---\left(2\right)$

其中，D为图像的维度。

不同的血管造影图像经过公式(1)的血管增强函数处理后，管状目标灰度值得到增强，背景噪声得到抑制，从而使得增强后的图像的直方图特性达到了相对统一。

步骤S102、对增强后的直方图曲线拟合血管数据和拟合噪声数据得到增强后的图像的混合模型；

上述步骤中血管类主要分布在高亮度部分，且以某一个均值为中心存在一定的波动，因此，本发明第一较佳实施方式可以采用一个高斯分布拟合血管。背景噪声类主要分布在低灰度值部分，通过对图2进行分析及多次实验测试，本发明最终选用两个指数分布来拟合背景噪声。因此，本发明所提出的混合模型可以表示为式(3)：

$f\left(x\right)=\underset{l=1}{\overset{2}{\Σ}}{w}_{El}{f}_{El}\left(x\right)+w_G{f}_G\left(x\right)---\left(3\right)$

其中，f_El(x)(l＝1,2)和f_G(x)分别表示指数分布和高斯分布，如式(4)所示。λ_El为指数分布的率参数，μ_G和σ² _G分别为高斯分布的期望和方差。f(x)为增强后的图像整体的混合概率密度函数，即混合模型，w_El(l＝1,2)可以为指数分布所占的比例，w_G可以为高斯分布所占的比例，且满足

$\left\{\begin{array}{c}{f}_{El}\left(x\right)={\mathrm{\λ}}_{El}\exp (-{\mathrm{\λ}}_{El}x),l\∈\[1,2\]\\ f_G\left(x\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_G}\exp \left(\frac{-{(x-{\mathrm{\μ}}_G)}^2}{2{{\mathrm{\σ}}^2}_G}\right)\end{array}\right.---\left(4\right)$

步骤S103、对混合模型中的参数进行参数估计处理；

本步骤S103首先利用k均值聚类进行混合模型参数的初始估计：当原始血管造影图像的峰值点超过设定个数(例如3个)时，则用设定个数(例如3个)峰值点完成k均值的初始化完成参数估计处理，否则采用如下规则进行K均值初始化

利用如下规则进行k均值算法的初始化：

$\left\{\begin{array}{c}{\mathrm{\μ}}_{k\_1}^{init}=I_{\max }/1000\\ {\mathrm{\μ}}_{k\_2}^{init}=I_{\max }/100\\ {\mathrm{\μ}}_{k\_3}^{init}=I_{\max }/10\end{array}\right.---\left(5\right)$

其中，I_max是原始图像中像素最大值。和表示k均值算法的初始聚类中心。μ_{k_l},σ_{k_l},w_{k_l},(l＝1,2,3)是由k均值聚类算法估计得到的三个分布的模型参数。

为了提高参数的估计精度，本发明用EM算对k均值估计得到的参数进行优化。如式(6)所示，k均值得到的估计结果作为EM算法迭代的初始值。

$\left\{\begin{array}{c}{\mathrm{\λ}}_{E1}^{init}=1.0/{\mathrm{\μ}}_{k\_1},w_{E1}^{init}=w_{k\_1}\\ {\mathrm{\λ}}_{E2}^{init}=1.0/{\mathrm{\μ}}_{k\_2},=w_2^{init}=w_{k\_2}\\ {\mathrm{\μ}}_G^{init}={\mathrm{\μ}}_{k\_3},{\mathrm{\σ}}_G^{init}={\mathrm{\σ}}_{k\_3},w_G^{init}=w_{k\_3}\end{array}\right.---\left(6\right)$

由EM算法得到的迭代更新公式如公式(7)和(8)：

$\left\{\begin{array}{c}{\mathrm{\λ}}_{Ei}^{k+1}=\frac{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(Ei\right|x_j)}{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(Ei\right|x_j)x_j}\\ w_{Ei}^{k+1}=\frac{1}{N}\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(Ei\right|x_j)\end{array}\right.,(i=1,2)---\left(7\right)$

$\left\{\begin{array}{c}{\mathrm{\μ}}_G^{k+1}=\frac{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j)x_j}{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j)},w_G^{k+1}=\frac{1}{N}\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j)\\ {\left({\mathrm{\σ}}_G^{k+1}\right)}^2=\frac{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j){(x_j-{\mathrm{\μ}}_G^{k+1})}^2}{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j)}\end{array}\right.---\left(8\right)$

其中，x_j是第j个体素的强度值，N是体素的总数。后验概率f^k(Ei|x_j)和f^k(G|x_j)可以由贝叶斯公式计算如下：

$\left\{\begin{array}{c}{f}^k\left(Ei\right|x_j)=\frac{w_{Ei}^k{f}^k\left(x_j\right|Ei)}{\underset{i=1}{\overset{2}{\Σ}}{w}_{Ei}^k{f}^k\left(x_j\right|Ei)+w_G^k{f}^k\left(x_j\right|G)}\\ f^k\left(G\right|x_j)=\frac{w_G^k{f}^k\left(x_j\right|G)}{\underset{i=1}{\overset{2}{\Σ}}{w}_{Ei}^k{f}^k\left(x_j\right|Ei)+w_G^k{f}^k\left(x_j\right|G)}\end{array}\right.---\left(9\right).$

步骤S104、通过最大后验概率估计法则将经过参数估计处理的混合模型中的血管从图像中分离出来。

由最大后验概率估计法则，根据公式(10)可以最终把血管从图像中分离出来：

w_Gf_G(x_j)＞w_E1f_E1(x_j)+w_E2f_E2(x_j)(10)

本发明第一较佳实施方式提供的技术方案将多尺度滤波和统计学方法相结合，提出了一种普适的血管分割方法。首先，不同模态的图像数据通过多尺度滤波进行血管增强，增强后的数据具有一个新的相近的统计学特性。针对滤波后的数据，本发明用一个高斯函数拟合血管类，同时根据直方图特性选用两个指数函数拟合背景噪声类，模型参数选用最大期望(EM)算法进行估算。最后，本发明利用最大后验概率分类将血管和背景标记出来。

本发明第一较佳实施方式提出了一种普适的血管分割方法：首先利用多尺度滤波进行血管增强和背景抑制，具有不同统计学特性的血管造影图像经过血管增强后达到较为一致的统计学特性。其次，针对滤波增强后的数据，本发明第一较佳实施方式用一个高斯和两个指数组合的混合模型对其进行拟合，并获得较好的拟合效果。最后，由最大后验概率完成了血管和背景噪声的分割。大量的仿真数据和临床数据证明了本发明所提出的方法的有效性和鲁棒性。

图3为本发明进行测试的仿真数据，PhantomⅠ是孔径渐变的管状目标，PhantomⅡ是由多个管状目标组成的球状体。本发明分第一较佳实施方式别对管状目标叠加了数量不等的高斯噪声进行测试。其中，图3中，(a)可以为真实目标；(b)可以为真实目标的一个切片；(c)可以为加噪声的切片。

图4为本发明进行测试用到的三维临床数据，包括：脑部MRA数据，心脏CTA数据，胎盘MRA数据和腿部MRA数据。如图4所示，(a)可以为脑部MRA；(b)可以为心脏CTA；(c)胎盘可以为MRA；(d)腿部可以为MRA。

由图5可以看出，不同的血管图像具有不同的统计学特性。本发明中用血管增强算法进行图像预处理，经滤波增强后的图像直方图曲线达到了一致。(a)可以为PhantomⅠ的图像直方图曲线；(b)可以为PhantomⅡ图像直方图曲线；(c)可以为脑部MRA图像直方图曲线；(d)可以为心脏CTA图像直方图曲线；(e)可以为胎盘MRA图像直方图曲线；(f)可以为腿部MRA图像直方图曲线。

由图6中的直方图拟合结果可以看出，本发明第一较佳实施方式提出的一个高斯分布和两个指数分布组成的混合模型对增强后的数据具有很好的拟合效果。(a)可以为PhantomⅠ直方图拟合结果示意图；(b)可以为PhantomⅡ直方图拟合结果示意图；(c)可以为脑部MRA直方图拟合结果示意图；(d)可以为心脏CTA直方图拟合结果示意图；(e)可以为胎盘MRA直方图拟合结果示意图；(f)可以为腿部MRA直方图拟合结果示意图。

图7和图8给出了仿真数据和临床数据的分割结果，可见，本发明第一较佳实施方式的技术方案的基于统计学的混合模型和传统的统计学模型相比，应用范围非常广泛。图7中(a1)-(a2)PhantomⅠ分割后的切片和三维效果；图7中(b1)-(b2)PhantomⅡ分割后的切片和三维效果。图8中(a1)-(a2)脑部MRA分割后的MIP图像和三维效果；图8中(b1)-(b2)心脏CTA分割后的MIP图像和三维效果；图8中(c1)-(c2)胎盘MRA分割后的MIP图像和三维效果；图8中(d1)-(d2)腿部MRA分割后的MIP图像和三维效果。

如图9所示，与传统方法对比：(a1)-(b1)为采用现有技术中方法：对原始三维图像直方图进行拟合，用一个高斯分布拟合血管类，一个瑞利分布和两个高斯分布共同拟合背景类；(a2)-(b2)本发明第一较佳实施方式提供的方法：对血管增强后的图像进行拟合，用一个高斯分布拟合血管类，两个指数分布拟合背景类；两套来自不同成像仪器的数据的分割结果表明：本发明第一较佳实施方式中所提出的分割方法对具有不同统计学特性的数据具有很好的鲁棒性。

本发明第二较佳实施方式提供一种普适的血管分割系统，所述系统如图10所示，包括：

增强单元1001，用于对输入的血管造影原始图像数据进行多尺度滤波后获得增强后的图像；

拟合单元1002，用于对增强后的直方图曲线拟合血管数据和拟合噪声数据得到增强后的图像的混合模型；

参数估计单元1003，用于对混合模型中的参数进行参数估计处理；

分离单元1004，用于通过最大后验概率估计法则将经过参数估计处理的混合模型中的血管从图像中分离出来。

可选的，增强单元1001具体，用于多尺度滤波算法通过尺度变换，利用Hessian矩阵的特征值计算局部血管出现的概率，所述血管增强函数如下：

$v_0\left(s\right)=\left\{\begin{array}{c}0\\ (1-\exp (-\frac{{R_A}^2}{2{\mathrm{\α}}^2}\left)\right)\exp (-\frac{{R_B}^2}{2{\mathrm{\β}}^2})(1-\exp (-\frac{S^2}{2c^2}\left)\right)\end{array}\right.---\left(1\right)$

其中，R_A，R_B，S为三个测度函数，如公式(2)所示，且R_A用来区分片状和线状结构，R_B用来区分点状结构和线状结构，S用于区分背景像素，α,β和c分别可以为个测度的阈值用于控制；

$\left\{\begin{array}{c}{R}_A=\frac{\left|h_2\right|}{\left|h_3\right|}\\ R_B=\frac{\left|h_1\right|}{\sqrt{\left|h_2{h}_3\right|}}\\ S=\left|\right|H|{|}_F=\sqrt{\underset{j\≤D}{\Σ}{h}_j^2}\end{array}\right..$

可选的，拟合单元1002具体，用于

混合模型f(x)如式(3)所示：

$f\left(x\right)=\underset{l=1}{\overset{2}{\Σ}}{w}_{El}{f}_{El}\left(x\right)+w_G{f}_G\left(x\right)---\left(3\right)$

其中，f_El(x)(l＝1,2)和f_G(x)分别表示指数分布和高斯分布，如式(4)所示，w_El(l＝1,2)为指数分布所占的比例，w_G可以为高斯分布所占的比例，且满足 $\underset{l=1}{\overset{2}{\Σ}}{w}_{El}+w_G=1;$

$\left\{\begin{array}{c}{f}_{El}\left(x\right)={\mathrm{\λ}}_{El}\exp (-{\mathrm{\λ}}_{El}x),l\∈\[1,2\]\\ f_G\left(x\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_G}\exp \left(\frac{-{(x-{\mathrm{\μ}}_G)}^2}{2{{\mathrm{\σ}}^2}_G}\right)\end{array}\right.---\left(4\right).$

可选的，参数估计单元1003具体，用于

当原始血管造影图像的峰值点超过3个时，则用3个峰值点完成k均值的初始化完成参数估计处理，否则，采用如下规则完成k均值算法的初始化；

$\left\{\begin{array}{c}{\mathrm{\μ}}_{k\_1}^{init}=I_{\max }/1000\\ {\mathrm{\μ}}_{k\_2}^{init}=I_{\max }/100\\ {\mathrm{\μ}}_{k\_3}^{init}=I_{\max }/10\end{array}\right.---\left(5\right)$

其中，I_max是原始图像中像素最大值，和表示k均值算法的初始聚类中心，μ_{k_l},σ_{k_l},w_{k_l},(l＝1,2,3)是由k均值聚类算法估计得到的三个分布的模型参数，

$\left\{\begin{array}{c}{\mathrm{\λ}}_{E1}^{init}=1.0/{\mathrm{\μ}}_{k\_1},w_{E1}^{init}=w_{k\_1}\\ {\mathrm{\λ}}_{E2}^{init}=1.0/{\mathrm{\μ}}_{k\_2},=w_2^{init}=w_{k\_2}\\ {\mathrm{\μ}}_G^{init}={\mathrm{\μ}}_{k\_3},{\mathrm{\σ}}_G^{init}={\mathrm{\σ}}_{k\_3},w_G^{init}=w_{k\_3}\end{array}\right.---\left(6\right)$

由EM算法得到的迭代更新公式如公式(7)和(8)：

$\left\{\begin{array}{c}{\mathrm{\λ}}_{Ei}^{k+1}=\frac{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(Ei\right|x_j)}{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(Ei\right|x_j)x_j}\\ w_{Ei}^{k+1}=\frac{1}{N}\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(Ei\right|x_j)\end{array}\right.,(i=1,2)---\left(7\right)$

$\left\{\begin{array}{c}{\mathrm{\μ}}_G^{k+1}=\frac{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j)x_j}{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j)},w_G^{k+1}=\frac{1}{N}\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j)\\ {\left({\mathrm{\σ}}_G^{k+1}\right)}^2=\frac{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j){(x_j-{\mathrm{\μ}}_G^{k+1})}^2}{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G\right|x_j)}\end{array}\right.---\left(8\right)$

其中，x_j是第j个体素的强度值，N是体素的总数，后验概率f^k(Ei|x_j)和f^k(G|x_j)由贝叶斯公式计算如下：

$\left\{\begin{array}{c}{f}^k\left(Ei\right|x_j)=\frac{w_{Ei}^k{f}^k\left(x_j\right|Ei)}{\underset{i=1}{\overset{2}{\Σ}}{w}_{Ei}^k{f}^k\left(x_j\right|Ei)+w_G^k{f}^k\left(x_j\right|G)}\\ f^k\left(G\right|x_j)=\frac{w_G^k{f}^k\left(x_j\right|G)}{\underset{i=1}{\overset{2}{\Σ}}{w}_{Ei}^k{f}^k\left(x_j\right|Ei)+w_G^k{f}^k\left(x_j\right|G)}\end{array}\right.---\left(9\right).$

可选的，分离单元1004具体，用于

由最大后验概率估计法则，根据公式(10)可以最终把血管从图像中分离出来：

w_Gf_G(x_j)＞w_E1f_E1(x_j)+w_E2f_E2(x_j)(10)。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种普适的血管分割方法，其特征在于，所述方法包括如下步骤：

对输入的血管造影原始图像数据进行多尺度滤波后获得增强后的图像；

对增强后的直方图曲线拟合血管数据和拟合噪声数据得到增强后的图像的混合模型；

对混合模型中的参数进行参数估计处理；

通过最大后验概率估计法则将经过参数估计处理的混合模型中的血管从图像中分离出来。

2.根据权利要求1所述的方法，其特征在于，所述对输入的血管造影原始图像数据进行多尺度滤波后获得增强后的图像具体，包括：

多尺度滤波算法通过尺度变换，利用Hessian矩阵的特征值计算局部血管出现的概率，所述血管增强函数如下：

$v_0\left(s\right)=\left\{\begin{array}{c}0\\ \left(1-\exp \left(-\frac{{R_A}^2}{2{\mathrm{\α}}^2}\right)\right)\exp \left(-\frac{{R_B}^2}{2{\mathrm{\β}}^2}\right)\left(1-\exp \left(-\frac{S^2}{2c^2}\right)\right)\end{array}\right.---\left(1\right)$

其中，R_A，R_B，S为三个测度函数，如公式(2)所示，且R_A用来区分片状和线状结构，R_B用来区分点状结构和线状结构，S用于区分背景像素，α,β和c分别可以为个测度的阈值用于控制；

$\left\{\begin{array}{c}{R}_A=\frac{\left|h_2\right|}{\left|h_3\right|}\\ R_B=\frac{\left|h_1\right|}{\sqrt{\left|h_2{h}_3\right|}}\\ S=\left|\right|H|{|}_F=\sqrt{\underset{j\≤D}{\Σ}{h}_j^2}\end{array}\right..$

3.根据权利要求1所述的方法，其特征在于，所述对增强后的直方图曲线拟合血管数据和拟合噪声数据得到增强后的图像的混合模型具体，包括：

混合模型f(x)如式(3)所示：

$f\left(x\right)=\underset{l=1}{\overset{2}{\Σ}}{w}_{El}{f}_{El}\left(x\right)+w_G{f}_G\left(x\right)---\left(3\right)$

其中，f_El(x)(l＝1,2)和f_G(x)分别表示指数分布和高斯分布，如式(4)所示，w_El(l＝1,2)为指数分布所占的比例，w_G可以为高斯分布所占的比例，且满足 $\underset{l=1}{\overset{2}{\Σ}}{w}_{El}+w_G=1;$

$\left\{\begin{array}{c}{f}_{El}\left(x\right)={\mathrm{\λ}}_{El}\exp (-{\mathrm{\λ}}_{El}x),l\∈\[1,2\]\\ f_G\left(x\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_G}\exp \left(\frac{-{(x-{\mathrm{\μ}}_G)}^2}{2{{\mathrm{\σ}}^2}_G}\right)\end{array}\right.---\left(4\right).$

4.根据权利要求1所述的方法，其特征在于，所述对混合模型中的参数进行参数估计处理具体，包括：

当原始血管造影图像的峰值点超过3个时，则用前3个峰值点完成k均值的初始化完成参数估计处理，否则，采用如下规则完成k均值算法的初始化；

$\left\{\begin{array}{c}{\mathrm{\μ}}_{k\_1}^{init}=I_{\max }/1000\\ {\mathrm{\μ}}_{k\_2}^{init}=I_{\max }/100\\ {\mathrm{\μ}}_{k\_3}^{init}=I_{\max }/10\end{array}\right.---\left(5\right)$

其中，I_max是原始图像中像素最大值，和表示k均值算法的初始聚类中心，μ_{k_l},σ_{k_l},w_{k_l},(l＝1,2,3)是由k均值聚类算法估计得到的三个分布的模型参数，

$\left\{\begin{array}{c}{\mathrm{\λ}}_{E1}^{init}=1.0/{\mathrm{\μ}}_{k\_1},w_{E1}^{init}=w_{k\_1}\\ {\mathrm{\λ}}_{E2}^{init}=1.0/{\mathrm{\μ}}_{k\_2},w_{E2}^{init}=w_{k\_2}\\ {\mathrm{\μ}}_G^{init}={\mathrm{\μ}}_{k\_3},{\mathrm{\σ}}_G^{init}={\mathrm{\σ}}_{k\_3},w_G^{init}=w_{k\_3}\end{array}\right.---\left(6\right)$

由EM算法得到的迭代更新公式如公式(7)和(8)：

$\left\{\begin{array}{c}{\mathrm{\λ}}_{Ei}^{k+1}=\frac{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(Ei|x_j\right)}{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(Ei|x_j\right)x_j},\left(i=1,2\right)\\ w_{Ei}^{k+1}=\frac{1}{N}\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(Ei|x_j\right)\end{array}\right.---\left(7\right)$

$\left\{\begin{array}{c}{\mathrm{\μ}}_G^{k+1}=\frac{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G|x_j\right)x_j}{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G|x_j\right)},w_G^{k+1}=\frac{1}{N}\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G|x_j\right)\\ {\left({\mathrm{\σ}}_G^{k+1}\right)}^2=\frac{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G|x_j\right){\left(x_j-{\mathrm{\μ}}_G^{k+1}\right)}^2}{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G|x_j\right)}\end{array}\right.---\left(8\right)$

其中，x_j是第j个体素的强度值，N是体素的总数，后验概率f^k(Ei|x_j)和f^k(G|x_j)由贝叶斯公式计算如下：

$\left\{\begin{array}{c}{f}^k\left(Ei|x_j\right)=\frac{w_{Ei}^k{f}^k\left(x_j|Ei\right)}{\underset{i=1}{\overset{2}{\Σ}}{w}_{Ei}^k{f}^k\left(x_j|Ei\right)+w_G^k{f}^k\left(x_j|G\right)}\\ f^k\left(G|x_j\right)=\frac{w_G^k{f}^k\left(x_j|G\right)}{\underset{i=1}{\overset{2}{\Σ}}{w}_{Ei}^k{f}^k\left(x_j|Ei\right)+w_G^k{f}^k\left(x_j|G\right)}\end{array}\right.---\left(9\right).$

5.根据权利要求1所述的方法，其特征在于，所述通过最大后验概率估计法则将经过参数估计处理的混合模型中的血管从图像中分离出来具体，包括：

由最大后验概率估计法则，根据公式(10)可以最终把血管从图像中分离出来：

w_Gf_G(x_j)＞w_E1f_E1(x_j)+w_E2f_E2(x_j)(10)。

6.一种普适的血管分割系统，其特征在于，所述系统包括：

增强单元，用于对输入的血管造影原始图像数据进行多尺度滤波后获得增强后的图像；

拟合单元，用于对增强后的直方图曲线拟合血管数据和拟合噪声数据得到增强后的图像的混合模型；

参数估计单元，用于对混合模型中的参数进行参数估计处理；

分离单元，用于通过最大后验概率估计法则将经过参数估计处理的混合模型中的血管从图像中分离出来。

7.根据权利要求6所述的系统，其特征在于，所述增强单元具体，用于多尺度滤波算法通过尺度变换，利用Hessian矩阵的特征值计算局部血管出现的概率，所述血管增强函数如下：

$v_0\left(s\right)=\left\{\begin{array}{c}0\\ \left(1-\exp \left(-\frac{{R_A}^2}{2{\mathrm{\α}}^2}\right)\right)\exp \left(-\frac{{R_B}^2}{2{\mathrm{\β}}^2}\right)\left(1-\exp \left(-\frac{S^2}{2c^2}\right)\right)\end{array}\right.---\left(1\right)$

其中，R_A，R_B，S为三个测度函数，如公式(2)所示，且R_A用来区分片状和线状结构，R_B用来区分点状结构和线状结构，S用于区分背景像素，α,β和c分别可以为个测度的阈值用于控制；

$\left\{\begin{array}{c}{R}_A=\frac{\left|h_2\right|}{\left|h_3\right|}\\ R_B=\frac{\left|h_1\right|}{\sqrt{\left|h_2{h}_3\right|}}\\ S=\left|\right|H|{|}_F=\sqrt{\underset{j\≤D}{\Σ}{h}_j^2}\end{array}\right..$

8.根据权利要求6所述的系统，其特征在于，所述拟合单元具体，用于混合模型f(x)如式(3)所示：

$f\left(x\right)=\underset{l=1}{\overset{2}{\Σ}}{w}_{El}{f}_{El}\left(x\right)+w_G{f}_G\left(x\right)---\left(3\right)$

其中，f_El(x)(l＝1,2)和f_G(x)分别表示指数分布和高斯分布，如式(4)所示，w_El(l＝1,2)为指数分布所占的比例，w_G可以为高斯分布所占的比例，且满足 $\underset{l=1}{\overset{2}{\Σ}}{w}_{El}+w_G=1;$

$\left\{\begin{array}{c}{f}_{El}\left(x\right)={\mathrm{\λ}}_{El}\exp (-{\mathrm{\λ}}_{El}x),l\∈\[1,2\]\\ f_G\left(x\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_G}\exp \left(\frac{-{(x-{\mathrm{\μ}}_G)}^2}{2{{\mathrm{\σ}}^2}_G}\right)\end{array}\right.---\left(4\right).$

9.根据权利要求6所述的系统，其特征在于，所述参数估计单元具体，用于

当原始血管造影图像的峰值点超过3个时，则用前3个峰值点完成k均值的初始化完成参数估计处理，否则，采用如下规则完成k均值算法的初始化；

$\left\{\begin{array}{c}{\mathrm{\μ}}_{k\_1}^{init}=I_{\max }/1000\\ {\mathrm{\μ}}_{k\_2}^{init}=I_{\max }/100\\ {\mathrm{\μ}}_{k\_3}^{init}=I_{\max }/10\end{array}\right.---\left(5\right)$

其中，I_max是原始图像中像素最大值，和表示k均值算法的初始聚类中心，μ_{k_l},σ_{k_l},w_{k_l},(l＝1,2,3)是由k均值聚类算法估计得到的三个分布的模型参数，

$\left\{\begin{array}{c}{\mathrm{\λ}}_{E1}^{init}=1.0/{\mathrm{\μ}}_{k\_1},w_{E1}^{init}=w_{k\_1}\\ {\mathrm{\λ}}_{E2}^{init}=1.0/{\mathrm{\μ}}_{k\_2},w_{E2}^{init}=w_{k\_2}\\ {\mathrm{\μ}}_G^{init}={\mathrm{\μ}}_{k\_3},{\mathrm{\σ}}_G^{init}={\mathrm{\σ}}_{k\_3},w_G^{init}=w_{k\_3}\end{array}\right.---\left(6\right)$

由EM算法得到的迭代更新公式如公式(7)和(8)：

$\left\{\begin{array}{c}{\mathrm{\λ}}_{Ei}^{k+1}=\frac{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(Ei|x_j\right)}{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(Ei|x_j\right)x_j},\left(i=1,2\right)\\ w_{Ei}^{k+1}=\frac{1}{N}\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(Ei|x_j\right)\end{array}\right.---\left(7\right)$

$\left\{\begin{array}{c}{\mathrm{\μ}}_G^{k+1}=\frac{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G|x_j\right)x_j}{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G|x_j\right)},w_G^{k+1}=\frac{1}{N}\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G|x_j\right)\\ {\left({\mathrm{\σ}}_G^{k+1}\right)}^2=\frac{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G|x_j\right){\left(x_j-{\mathrm{\μ}}_G^{k+1}\right)}^2}{\underset{j=1}{\overset{N}{\Σ}}{f}^k\left(G|x_j\right)}\end{array}\right.---\left(8\right)$

其中，x_j是第j个体素的强度值，N是体素的总数，后验概率f^k(Ei|x_j)和f^k(G|x_j)由贝叶斯公式计算如下：

$\left\{\begin{array}{c}{f}^k\left(Ei|x_j\right)=\frac{w_{Ei}^k{f}^k\left(x_j|Ei\right)}{\underset{i=1}{\overset{2}{\Σ}}{w}_{Ei}^k{f}^k\left(x_j|Ei\right)+w_G^k{f}^k\left(x_j|G\right)}\\ f^k\left(G|x_j\right)=\frac{w_G^k{f}^k\left(x_j|G\right)}{\underset{i=1}{\overset{2}{\Σ}}{w}_{Ei}^k{f}^k\left(x_j|Ei\right)+w_G^k{f}^k\left(x_j|G\right)}\end{array}\right.---\left(9\right).$

10.根据权利要求6所述的系统，其特征在于，所述分离单元具体，用于

由最大后验概率估计法则，根据公式(10)可以最终把血管从图像中分离出来：

w_Gf_G(x_j)＞w_E1f_E1(x_j)+w_E2f_E2(x_j)(10)。