CN104637044A - 钙化斑块及其声影的超声图像提取系统 - Google Patents

Abstract

本发明提供了一种钙化斑块及其声影的超声图像提取系统，所述系统包括：血管内超声图像获取模块，用于获取血管内超声图像；感兴趣区域确定模块，用于根据所述血管内超声图像确定包含血管组织图像的感兴趣区域；第一图像获取模块，用于以所述感兴趣区域的中心像素为直角坐标系的坐标原点，将所述感兴趣区域变换到极坐标系下，获得第一图像；第一区域获取模块，用于根据所述第一图像判断包括钙化斑块和声影的图像的区域作为第一区域；提取模块，用于根据所述第一区域提取钙化斑块及其声影的图像。本发明提供的系统实现了钙化斑块图像及钙化斑块的声影图像的自动提取，极大地提高了提取钙化斑块及其声影的图像的效率。

Description

钙化斑块及其声影的超声图像提取系统

技术领域

本发明涉及生物信息技术领域，特别是涉及一种钙化斑块及其声影的超声图像提取系统。

背景技术

心脑血管疾病已成为人类健康的头号杀手。动脉粥样硬化及其并发症引起心脑血管疾病的常见原因。钙化斑块的大小与形状是临床诊断中衡量动脉粥样硬化程度的重要指标。

为了方便获知钙化斑块的位置信息，常通过超声检查进行拍摄图像。传统的提取钙化斑块图像的方法是通过人工观察，由富有经验的医生在拍摄的血管内超声(Intravascular ultrasound,IVUS)图像中判断出钙化斑块的位置和大小。由于为同一个病人拍摄的血管内超声图像经常包含上千幅图像，因而由医生人工观察的方法不但耗时，而且重复性差，容易受到医生经验和主观因素的影响。

发明内容

基于此，有必要由人工从血管内超声图像中判断出钙化斑块的位置效率低的问题，提供一种钙化斑块及其声影的超声图像提取系统。

一种钙化斑块及其声影的超声图像提取系统，所述系统包括：

血管内超声图像获取模块，用于获取血管内超声图像；

感兴趣区域确定模块，用于根据所述血管内超声图像确定包含血管组织图像的感兴趣区域；

第一图像获取模块，用于以所述感兴趣区域的中心像素为直角坐标系的坐标原点，将所述感兴趣区域变换到极坐标系下，获得第一图像；

第一区域获取模块，用于根据所述第一图像判断包括钙化斑块和声影的图像的区域作为第一区域；

提取模块，用于根据所述第一区域提取钙化斑块及其声影的图像。

上述钙化斑块及其声影的超声图像提取系统，在血管内超声图像中确定包含血管组织图像的感兴趣区域In1后，将感兴趣区域In1变换到极坐标系下，获得第一图像I。再将第一图像I中判断为包括钙化斑块和声影的图像的第一区域R^mrf，从第一区域R^mrf中提出钙化斑块图像和钙化斑块的声影图像。实现了钙化斑块图像及钙化斑块的声影图像的自动提取，极大地提高了提取钙化斑块及其声影的图像的效率。

附图说明

图1为一个实施例中钙化斑块及其声影的图像提取方法的流程示意图；

图2为一个实施例中获取第二图像的步骤的流程示意图；

图3为一个实施例中根据所述第一图像I判断包括钙化斑块和声影的图像的区域作为第一区域的步骤的流程示意图；

图4为一个实施例中提取第一图像中每列的特征的步骤的流程示意图；

图5为一个实施例中根据所述第一图像I中每列的特征将所述第一图像I中的各列分为含有钙化斑块和声影的图像的列以及不含有钙化斑块和声影的图像的列两类的步骤的流程示意图；

图6为一个实施例中根据第一区域提取钙化斑块及其声影的图像的步骤的流程示意图；

图7为一个实施例中根据过滤后的第三区域提取声影图像所在的区域的步骤的流程示意图；

图8为一个实施例中钙化斑块及其声影的超声图像提取系统的结构框图；

图9为另一个实施例中钙化斑块及其声影的超声图像提取系统的结构框图；

图10为图9中一个实施例的聚类模块的结构框图；

图11为图10中一个实施例的参数求解模块的结构框图；

图12为图9中一个实施例的特征提取模块的结构框图；

图13为图9中一个实施例的分类模块的结构框图；

图14为图13中一个实施例的置信计算模块的结构框图；

图15为图9中一个实施例的提取模块的结构框图；

图16为图15中一个实施例的提取执行模块的结构框图；

图17为一个实施例中的血管内超声图像；

图18为一个实施例中在血管内超声图像中确定感兴趣区域的示意图；

图19为一个实施例中对感兴趣区域中的像素进行聚类获得聚类图像；

图20为一个实施例中在第一图像中确定最大值线的示意图；

图21为一个实施例中以聚类图像的中心像素为直角坐标系的坐标原点，将聚类图像变换到极坐标系，获得的第二图像；

图22为一个实施例中马尔可夫随机场的概率图模型的示意图；

图23为一个实施例中在第三图像中确定的属于第二区域的像素设为1，不属于第二区域的像素的灰度值设为0的示意图；

图24为一个实施例中第一图像每列对应的置信的示意图；

图25为一个实施例中将第一图像中的各列分为含有钙化斑块和声影图像的列和不含有钙化斑块和声影图像的列两类后，第一图像每列对应的置信的示意图；

图26为一个实施例中一个第三区域的示意图；

图27为一个实施例中在第一图像中提取的钙化斑块及其声影的图像的示意图；

图28为一个实施例中将第一图像转化到直角坐标系下后，提取的钙化斑块及其声影的图像的边界。

具体实施方式

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

如图1所示，在一个实施例中提供了一种钙化斑块及其声影的图像提取方法，包括：

步骤102，获取血管内超声图像。

血管内超声是将无创性的超声技术和有创性的心导管技术相结合，对心血管病变进行探测的一种方法。通过心导管将微型化的超声探头插入心血管腔内进行探测，再经电子成像系统可形成血管内超声图像，可表示心血管断面的形态和血流图形，如图17所示。

步骤104，根据所述血管内超声图像确定包含血管组织图像的感兴趣区域In1。

在血管内超声图像中，感兴趣区域In1为包括超声图像中血管内组织的所有信息的大致为圆形的区域。血管内超声图像中感兴趣区域In1之外的区域不包含有效信息，剔除这部分无效区域，可避免干扰。在血管内超声图像中，该感兴趣区域In1为一圆形区域，如图18所示，圆形1801内的区域就是感兴趣区域In1。

步骤106，以所述感兴趣区域In1的中心像素为直角坐标系的坐标原点，将所述感兴趣区域In1变换到极坐标系下，获得第一图像I。

以感兴趣区域In1的中心像素为直角坐标原点，将圆形的感兴趣区域In1从直角坐标系变换到极坐标系，得到矩形的第一图像I，如图20所示。令H表示第一图像I的行数，W表示第一图像I的列数。

步骤108，根据所述第一图像I判断包括钙化斑块和声影的图像的区域作为第一区域R^mrf。

第一图像I是矩形图像，可通过扫描第一图像I的每一列，通过先验知识判断第一图像I每一列中是否含有钙化斑块图像和钙化斑块的声影图像，从而初步判断出包括钙化斑块和声影的图像的第一区域其中表示第i个第一区域，N_mrf表示第一区域在第一图像I中的个数，N_mrf可以等于1。需要说明的是，这里所说的图像的一列含有钙化斑块图像和钙化斑块的声影图像，是指该列包含钙化斑块和声影的一部分图像。

由于声影是因钙化斑块阻挡超声信号，从而在血管内超声图像上形成的在钙化斑块后端留下的阴影，因此在感兴趣区域In1中，包含钙化斑块与其对应的声影的区域大致为一扇形区域。而在极坐标系下的第一图像I中，该扇形区域对应矩形的第一区域R^mrf，根据矩形的第一区域R^mrf判断是否包含钙化斑块及其声影的图像，与根据扇形区域判断相比，可以降低计算复杂度。

步骤110，根据第一区域R^mrf提取钙化斑块及其声影的图像。

确定了包含带声影的钙化斑块图像的第一区域R^mrf后，可使用边缘检测从第一区域R^mrf中提出钙化斑块图像和钙化斑块的声影图像。

上述钙化斑块及其声影的图像提取方法，在血管内超声图像中确定包含血管组织图像的感兴趣区域In1后，将感兴趣区域In1变换到极坐标系下，获得第一图像I。再将第一图像I中判断为包括钙化斑块和声影的图像的第一区域R^mrf，从第一区域R^mrf中提出钙化斑块图像和钙化斑块的声影图像。实现了钙化斑块图像及钙化斑块的声影图像的自动提取，极大地提高了提取钙化斑块及其声影的图像的效率。

如图2所示，在一个实施例中，该钙化斑块及其声影的图像提取方法还包括获取第二图像I_R的步骤，包括：

步骤202，对所述感兴趣区域In1的像素进行聚类。

对感兴趣区域In1的像素聚类的数学本质是为了数据降维。因为在原始数据中，灰度值的灰度范围是0-255，换成2进制码就是8位00000000–11111111，而聚类之后得到了K类，这个K通常很小，然后换成2进制码后就是<8位，实际上相当于对原始灰度数据进行了压缩，提取了关键特征。

具体地，步骤202包括以下步骤：

（1）、使用瑞利混合分布模型描述感兴趣区域In1的每个像素的概率分布。

具体地，感兴趣区域In1中的每个像素的概率分布如公式1所示：

则所述感兴趣区域的每个像素的概率分布为：

$p\left(y_i\right)=\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_{j}p\left(y_i\right|{\mathrm{\&theta;}}_j)$                                 ；公式1

$p\left(y_i\right|{\mathrm{\&theta;}}_j)=\frac{y_i-a_j}{{\mathrm{\&sigma;}}_j^2}\exp (-\frac{{(y_i-a_j)}^2}{{2\mathrm{\&sigma;}}_j^2})$

其中，Y＝{y₁,...,y_N}表示感兴趣区域In1的像素集，y_i表示感兴趣区域的第i个像素的灰度值，θ_j＝{a_j,σ_j}是瑞利分布的参数，具体地，a_j表示第j个瑞利分布在横轴上的平移量，σ_j表示瑞利分布的众数（mode，众数是指一组数据中出现次数最多的数值）；π_j表示瑞利混合分布中每个分量的权重，K表示瑞利混合分布中混合分量的个数，p(y_i)表示第i个像素的混合概率，p(y_i|θ_j)表示第i个像素属于第j类的概率。

对于公式1，瑞利混合分布的似然函数为：

$L\left(\mathrm{\&Theta;}\right)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log \left(\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_{j}p\left(y_i\right|{\mathrm{\&theta;}}_j)\right)$                公式2

其中，N表示感兴趣区域中的像素总数，Θ表示参数集；

然后，定义一个函数表示感兴趣区域In1中的第i个像素对瑞利混合分布中的第j个分量的权重：

${\mathrm{\&xi;}}_j\left(y_i\right)=\frac{{\stackrel{\&OverBar;}{y}}_i-c_j}{b_j^2}\exp (-\frac{{({\stackrel{\&OverBar;}{y}}_i-c_j)}^2}{{2b}_j^2})$                      公式3

其中c_j和b_j表示计算参数，具体地，c_j表示所满足的瑞利分布的横轴偏移量，b_j表示所满足的瑞利分布的众数，表示第i个像素的8邻域的均值，即可以写成

${\stackrel{\&OverBar;}{y}}_i=\frac{1}{9}\underset{m=-1}{\overset{1}{\mathrm{\&Sigma;}}}\underset{n=-1}{\overset{1}{\mathrm{\&Sigma;}}}y(u+m,v+n)$                               公式4

其中，坐标为(u,v)的像素为中心点，调整步长m、n可取遍该像素领域的各个点，y(u+m,v+n)表示血管内超声图像中坐标为(u+m,v+n)的点的灰度值。

举例说明，假定有一个3×3的窗口，(u,v)就是窗口在血管内超声图像上的中心坐标，假设该坐标为(100,100)，即u=100,v=100。而m、n就是从(-l,l)的连续变化量，表示对窗口中的数据进行遍历。在上面假设中l=1，那么m和n都是从-1到+1。这样，公式4就是求(99,99)、(99,100)、(99,101)、(100,99)、(100,100)、(100,101)、(101,99)、(101,100)和(101,101)这9个点的灰度值的和的平均值。

对于第i个像素的邻域N_i，定义一个对瑞利混合分布的第j个分量为对象的权重函数：

                                公式5

其中M是邻域N_i中像素的个数，α是一个控制变量用来控制公式5中值的大小。这里的邻域N_i是8-邻域。

最后定义一个新的先验概率π_ij，表示第i个像素的邻域N_i的权重函数在第j个瑞利分布中所占有权重。

                              公式6

（2）、使用最大期望算法（Expectation-Maximization算法，简称EM算法）求解所述瑞利混合模型的参数。具体地，包括以下步骤：

A、为了估计出公式1中混合模型的参数θ＝{a_j,σ_j},j＝1,...,K，需要最大化公式2中的似然函数，即

${\mathrm{\&Theta;}}^{*}=\arg \underset{\mathrm{\&Theta;}}{\max }L\left(\mathrm{\&Theta;}\right)$                       公式7

则使用EM算法得到的目标函数为：

$Q(\mathrm{\&Theta;},{\mathrm{\&Theta;}}^{\left(t\right)})=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\log \left({\mathrm{\&pi;}}_{j}P\left({\mathrm{\&theta;}}_j\right|y_i;{\mathrm{\&Theta;}}^{\left(t\right)})\right)+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\log \left(P\left(y_i\right|{\mathrm{\&theta;}}_j;\mathrm{\&Theta;})\right)P\left({\mathrm{\&theta;}}_j\right|y_i;{\mathrm{\&Theta;}}^{\left(t\right)});$      公式8

其中，θ_j＝{a_j,σ_j,c_j,b_j,α}是参数向量，具体地，a_j表示第j个瑞利分布在横轴上的平移量，σ_j表示瑞利分布的众数，c_j和b_j是计算参数，具体地，c_j表示所满足的瑞利分布的横轴偏移量，b_j表示所满足的瑞利分布的众数，α表示控制变量；π_j是先验概率，P(θ_j|y_i;Θ^(t))是后验概率，P(y_i|θ_j;Θ)是类条件概率密度，Θ^(t)表示第t次迭代中已知的参数集，Θ表示第t次迭代中未知的参数集。

B、初始化参数集Θ为Θ⁽⁰⁾。

初始化参数集Θ，记作Θ⁽⁰⁾，给定分类个数K＝5和参数α＝10^-7，使用K均值算法（一种聚类算法，算法步骤为：1、输入数据集合和类别数K(用户指定)；2、随机分配类别中心点的位置；3、将每个点放在离他最近的类别中心点所在的集合；4、移动类别中心点到它所在集合的中心；5、转到第3步，直到收敛。）计算出感兴趣区域In1每一类的均值

令a_j表示上一步的K均值分类后，第j类的所有像素中的最小灰度值。

从下式公式9中得到 ${\mathrm{\&sigma;}}^{\left(0\right)}=[{\mathrm{\&sigma;}}_1^{\left(0\right)},{\mathrm{\&sigma;}}_2^{\left(0\right)},...,{\mathrm{\&sigma;}}_K^{\left(0\right)}],$

${\mathrm{\&sigma;}}^{\left(0\right)}=(m^{\left(0\right)}-a^{\left(0\right)})\sqrt{\frac{2}{\mathrm{\&pi;}}}$                                       公式9

令c⁽⁰⁾＝a⁽⁰⁾，b⁽⁰⁾＝σ⁽⁰⁾。

C、根据初始化的参数集Θ⁽⁰⁾计算参数向量，并使用最速下降算法更新所述参数向量，直到最大期望算法收敛，获得最终参数集Θ^*。

参数集Θ的计算使用到了最速下降法(steepest-descent method)：

${\mathrm{\&Theta;}}^{(t+1)}={\mathrm{\&Theta;}}^{\left(t\right)}-\mathrm{\&eta;}\frac{\&PartialD;Q(\mathrm{\&Theta;},{\mathrm{\&Theta;}}^{\left(t\right)})}{\&PartialD;\mathrm{\&Theta;}}$                         公式10

通过公式1和公式6，由贝叶斯定理可以得到像素i在模型中的后验概率为公式11，表示第i个像素属于第j类的概率：

$P^{\left(t\right)}\left({\mathrm{\&theta;}}_j\right|y_i)=\frac{{\mathrm{\&pi;}}_{\mathrm{ij}}^{\left(t\right)}{P}^{\left(t\right)}\left(y_i\right|{\mathrm{\&theta;}}_j)}{\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_{\mathrm{ik}}^{\left(t\right)}{P}^{\left(t\right)}\left(y_i\right|{\mathrm{\&theta;}}_k)}$                           公式11

计算参数向量为：

$\frac{\&PartialD;Q}{\&PartialD;\mathrm{\&Theta;}}={\left[\begin{array}{ccccc}\frac{\&PartialD;Q}{\&PartialD;a}& \frac{\&PartialD;Q}{\&PartialD;\mathrm{\&sigma;}}& \frac{\&PartialD;Q}{\&PartialD;c}& \frac{\&PartialD;Q}{\&PartialD;b}& \frac{\&PartialD;Q}{\&PartialD;\mathrm{\&alpha;}}\end{array}\right]}^T$                       公式12

且为

$\frac{\&PartialD;Q}{\&PartialD;a_j}=-\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}P\left({\mathrm{\&theta;}}_j\right|y_i)(\frac{y_i-a_j}{{\mathrm{\&sigma;}}_j^2}-\frac{1}{y_i-a_j})$                       公式13

$\frac{\&PartialD;Q}{\&PartialD;{\mathrm{\&sigma;}}_j}=-\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}P\left({\mathrm{\&theta;}}_j\right|y_i)(-\frac{2}{{\mathrm{\&sigma;}}_j}+\frac{{(y_i-a_j)}^2}{{\mathrm{\&sigma;}}_j^3})$                    公式14

     公式15

       公式16

             公式17

然后可由公式10更新参数向量。

当公式10中的参数向量不再变换时，EM算法收敛，记此时计算出的参数集为Θ^*；否则，令Θ^(t)＝Θ^(t+1)，然后通过公式1和公式6，由贝叶斯定理可以得到像素i在模型中的后验概率，继续进行计算。

（3）、使用最大后验概率准则对所述感兴趣区域In1的像素进行聚类。

得到最终的参数集Θ^*后，使用最大后验概率准则去对感兴趣区域In1中的每个像素进行聚类；聚类个数等于瑞利混合分布的分量的个数κ。把第i个像素归为第j类，如果

          公式18

其中，表示正整数集，K表示聚类个数。

步骤204，将聚类结果中属于同一聚类簇的像素的灰度值置为相同的值，且属于不同聚类簇的像素的灰度值不同，获得聚类图像。

具体地，将感兴趣区域In1聚类为K类，每一类称为一个聚类簇。将同一聚类簇的像素的灰度值置为相同的值，不同聚类簇的像素的灰度值置为不同的值，且不同聚类簇的像素的灰度值所置的值可与聚类簇中像素的灰度值均值正相关。进一步地，若聚类结果中像素属于第k类，可将该像素的灰度值设置为k。由于最后获得的聚类图像中就只有K种灰度值，达到了降维的效果，可以降低计算复杂度。获得的聚类图像如图19所示。

步骤206，以所述聚类图像的中心像素为直角坐标系的坐标原点，将所述聚类图像变换到极坐标系，获得第二图像I_R。

将把感兴趣区域In1聚类降维后获得的大致为圆形的聚类图像变换到极坐标系，获得的矩形的第二图像I_R，如图21所示，便于后续计算。

本实施例中对感兴趣区域In1进行聚类降维，可降低计算复杂度。

在一个实施例中，该方法还包括步骤：确定所述第一图像I每列的像素中灰度值最大的像素，获得最大值线L^mvl。

具体地，令C＝{1,2,...,W}表示第一图像I的列编号的集合，对于任意i∈C，把I的第i列的具有最大灰度值的像素（比如在图像中最白的像素）的行数记作于是对于第一图像I的每一列，可以得到记这里L^mvl即为最大值线。如附图20中线段2201所示。第一图像I中行坐标从左到右递增，列坐标从上到下递增。

如图3所示，在一个实施例中，步骤108包括：

步骤302，提取所述第一图像I中每列的特征。

第一图像I本身包括大量的像素，为了方便计算，需要对第一图像I中每列提取特征，降低计算复杂度。

如图4所示，在一个实施例中，步骤302具体包括：

步骤402，计算所述第二图像I_R中每个聚类簇对应的感兴趣区域In1中的像素的灰度值均值，找到灰度值均值最小的预设个数的聚类簇，将所述找到的聚类簇中的像素构成的区域作为第二区域R^D。

由于第二图像I_R是有K类，计算出每类的灰度值均值，找到灰度值均值最小的两类，将这两类中所有像素的位置定义为第二区域R^D。当预设个数为2时，能够达到最好的效果，识别正确率显著提高。

步骤404，生成与所述第二图像I_R尺寸相同的第三图像I_D，将所述第三图像I_D中属于第二区域R^D的像素的灰度值和不属于第二区域R^D的像素的灰度值分别设置为不同的值。

可定义第三图像I_D，其大小与第二图像I_R相同，第三图像I_D中属于第二区域R^D的像素的灰度值设为1（如图23中白色区域2301），不属于第二区域R^D的像素的灰度值设为0。

步骤406，在所述第三图像I_D中确定第二区域R^D的上边界。

可从上向下扫描第三图像I_D的第i列(i＝1,2,...,W)，找到属于第二区域R^D的第一个像素（灰度值为1的第一个像素），把该像素的行坐标记作于是，我们可以在第三图像I_D中定义第二区域R^D的上边界，记为

步骤408，计算所述第三图像I_D每列的底部与所述第二区域R^D的上边界L^ubrd的距离、所述第二区域R^D的上边界与所述最大值线L^mvl的距离以及所述第一图像I中每列对应的所述第二区域的上边界与所述最大值线L^mvl之间的像素的平均灰度值。

计算三个自定义的特征参数F^h,F^d,F^v，其中，表示第三图像I_D每列的底部与所述第二区域R^D的上边界L^ubrd的距离。

$F_i^h=H-p_i^{\mathrm{ubrd}},i=\mathrm{1,2},...,W$                  公式20

表示第二区域R^D的上边界L^ubrd与最大值线L^mvl的距离，

$F_i^d=p_i^{\mathrm{ubrd}}-p_i^{\mathrm{mvl}},i=\mathrm{1,2},...,W$                    公式21

表示第一图像I中每列对应的第二区域R^D的上边界L^ubrd与最大值线L^mvl之间的像素的平均灰度值。

$F_i^v=\frac{1}{H-p_i^{\mathrm{ubrd}}+1}\underset{j=p_i^{\mathrm{ubrd}}}{\overset{p_i^{\mathrm{mvl}}}{\mathrm{\&Sigma;}}}I(p_i^{\mathrm{ubrd}},i),i=\mathrm{1,2},...,W$              公式22

步骤410，根据所述第三图像I_D每列的底部与所述第二区域R^D的上边界L^ubrd的距离、所述第二区域R^D的上边界与所述最大值线L^mvl的距离以及所述第一图像I中每列对应的所述第二区域的上边界与所述最大值线L^mvl之间的像素的平均灰度值确定所述第一图像每列的特征。

利用上面三个参数可以定义F＝{F₁,F₂,...,F_W}为第一图像I每列的特征：

$F=\frac{h_1{F}^h+h_2{F}^d+h_3{F}^v}{255}$                  公式23

F_i＝{F₁,F₂,...,F_W}，其中h₁,h₂,h₃是权重,优选的,权重h₁,h₂,h₃的值分别为h₁＝5,h₂＝-0.5,h₃＝-1。当h₁＝5,h₂＝-0.5,h₃＝-1时，第一图像I每列的特征能够很好地反映该列是否含有钙化斑块图像和声影的图像，性能最优。

步骤304，根据所述第一图像I中每列的特征将所述第一图像I中的各列分为含有钙化斑块和声影的图像的列以及不含有钙化斑块和声影的图像的列两类。

如图5所示，在一个实施例中，步骤304具体包括：

步骤502，根据所述第一图像I每列的特征，使用置信传播算法（BeliefPropagation算法）在马尔可夫随机场上计算出第一图像I每列含有钙化斑块和声影的图像的置信。具体地，步骤502包括以下四个步骤：

A1、定义马尔可夫随机场的位置集S和状态集L分别为

S＝{1,2,...,W}                       公式19

L＝{-1,+1}

其中，S＝{1,2,...,W}表示该列的位置；L＝{-1,+1}表示该列的状态，若第一图像I的一列的状态为“+1”，表示该列含有钙化斑块及其声影的图像；如果第一图像I的一列的状态为“-1”，则表示该列不含有钙化斑块及其声影的图像。

马尔可夫随机场可由概率图模型表示，如图22所示；图中节点χ₁,...,χ_W是观测变量，表示图像I每列的特征，节点z₁,...,z_W是隐藏变量，表示图像I每列的状态。

A2、对于马尔可夫随机场上的每个隐藏变量z_i，令所有的初始状态的概率满足均匀分布，即隐藏变量z_i的边缘概率为：

P(z_i＝-1)＝P(z_i＝1)＝0.5                           公式24

且隐藏变量z_i的置信初始化为：

$b_i^{\left(0\right)}(z_i=-1)=b_i^{\left(0\right)}(z_i=1)=0.5$               公式25

局部信息φ(z_i,χ_i)由公式22和公式23决定：

φ_i(z_i,χ_i)＝F_i                                     公式26

F_i是第一图像第i列的特征，χ_i表示观测变量χ₁,...,χ_W。

相容函数ψ(z_i,z_j)是一个2×2的矩阵

${\mathrm{\&psi;}}_{i,j}(z_i,z_j)=\left[\begin{array}{cc}0.8& 0.2\\ 0.2& 0.8\end{array}\right]$                                  公式27

当节点z_i和z_j是互为邻域，那么从z_i传递到z_j的信息初始化为：

m_i,j(z_j)＝1                                      公式28

A3、在第t次迭代中，计算从z_i传递到z_j的信息

$m_{i,j}^{\left(t\right)}\left(z_j\right)=\underset{z_i}{\mathrm{\&Sigma;}}{\mathrm{\&phi;}}_i(z_i,{\mathrm{\&chi;}}_i){\mathrm{\&psi;}}_{i,j}(z_i,z_j){\mathrm{\&Pi;}}_{k\&Element;N\left(i\right)\backslash j}{m}_{k,i}^{\left(t\right)}\left(z_i\right)$               公式29

并计算节点z_i的置信

$b_i^{\left(t\right)}\left(z_i\right)=k{\mathrm{\&phi;}}_i\left(z_i\right)\underset{j\&Element;N\left(i\right)}{\mathrm{\&Pi;}}{m}_{j,i}^{\left(t\right)}\left(z_i\right)$                            公式30

其中N(i)是节点z_i的邻域。

A4、当在t+1次迭代中，满足公式31时，迭代算法收敛。

$\frac{1}{W}|\underset{i=1}{\overset{W}{\mathrm{\&Sigma;}}}\left(b^{(t+1)}\left(z_i\right)\right)-\underset{i=1}{\overset{W}{\mathrm{\&Sigma;}}}\left(b^{\left(t\right)}\left(z_i\right)\right)|<\mathrm{\&epsiv;}$                        公式31

其中ε为预设定值。

则每个隐藏变量的置信，即每列含有钙化斑块和声影的图像的置信为：

b^*(z₁)＝b^(t+1)(z₁),b^*(z₂)＝b^(t+1)(z₂),...,b^*(z_W)＝b^(t+1)(z_W)        公式32

其中，W表示第一图像的列数。

最终第一图像I每列含有钙化斑块和声影的图像的置信如图24所示，图24中的线2401表示每列含有钙化斑块和声影的图像的置信，其中横轴表示第一图像I的列编号，纵轴表示置信，线2402表示置信等于0.5。

步骤504，根据所述第一图像I每列的置信将所述第一图像中的各列分为含有钙化斑块和声影的图像的列和不含有钙化斑块和声影的图像的列两类。

对第一图像I的列进行分类使用的是最大后验准则。具体地，对于第一图像I的第i列的状态z_i，如果满足：

b^*(z_i＝-1)＞b^*(z_i＝+1)               公式33

那么将第i列分到类“-1”。如果满足

b^*(z_i＝-1)≤b^*(z_i＝+1)                 公式34

那么将第i列分到类“+1”。其中，若该列属于类“+1”表示该列含有钙化斑块和声影的图像，若该列不属于类“-1”则表示该列不含有钙化斑块和声影的图像。

步骤306，将所述第一图像I中连续的含有钙化斑块和声影的图像的列构成的区域判断为包括钙化斑块和声影的图像的第一区域R^mrf。

具体地，可对第一图像I从左向右每列依次扫描，把连续都是类“+1”的列当作同一个区域，记第一区域其中表示第i个第一区域，N_mrf表示这种区域在第一图像I中的个数，N_mrf可以等于1，结果见图25，其中2501表示第一图像I各列的置信，线2502表示置信为0.5的界限。对于任意一个第一区域其最左边列称作左列边界，其最右边的列称作右列边界，其中置信为1的列属于类“+1”，置信为0的列属于类“-1”；类“+1”表示含有钙化斑块和声影的图像，类“-1”表示不含有钙化斑块和声影的图像。

如图6所示，在一个实施例中，步骤110具体包括：

步骤602，根据所述第一区域R^mrf和所述第一图像I的最大值线L^mvl确定第三区域R^cs。

第一区域R^mrf是初步判断包括钙化斑块和声影的图像的区域，但这种判断并不精确，为了精确判断出钙化斑块和声影图像所在位置，需要先根据所述第一区域R^mrf和所述第一图像I的最大值线L^mvl定义第三区域R^cs，再根据第三区域R^cs精确判断钙化斑块及其声影的图像所在的位置。

在第一图像I中定义最大值线的上升沿和下降沿。最大值线的上升沿是指最大值线上的一个连续线段；最大值线上的所有上升沿记作其中N_re表示上升沿的个数，表示最大值线上的第i个上升沿，且可以表示成

$L_i^{\mathrm{re}}=\{P_l^{{\mathrm{re}}_i},p_r^{{\mathrm{re}}_i}\}$                  公式35

其中，分别表示第i个上升沿的最左边点和最右边点在第一图像I中的列坐标；对于i＝1,2,...,N_re，上升沿有以下性质：

$p_{p_l^{{\mathrm{re}}_i}}^{\mathrm{mvl}}\&GreaterEqual;p_{p_l^{{\mathrm{re}}_i}+1}^{\mathrm{mvl}}\&GreaterEqual;...\&GreaterEqual;p_{p_r^{{\mathrm{re}}_i}}^{\mathrm{mvl}}$                   公式36

最大值线的下降沿可表示成 $L^{\mathrm{fe}}=\{L_1^{\mathrm{fe}},L_2^{\mathrm{fe}},...,L_{N_{\mathrm{fe}}}^{\mathrm{fe}}\},L_i^{\mathrm{fe}}=\{p_l^{{\mathrm{fe}}_i},p_r^{{\mathrm{fe}}_i}\},$ 其中N_fe表示下降沿的个数，分别表示第i个下降沿的最左边点和最右边点在第一图像I中的列坐标。对于i＝1,2,...,N_fe，下降沿有以下性质：

$p_{p_l^{{\mathrm{fe}}_i}}^{\mathrm{mvl}}\&le;p_{p_l^{{\mathrm{fe}}_i}+1}^{\mathrm{mvl}}\&le;...\&le;p_{p_r^{{\mathrm{fe}}_i}}^{\mathrm{mvl}}$                             公式37

然后筛选上升沿和下降沿。对于某个可能包含钙化斑块图像和钙化斑块的声影图像的区域其左列边界必定和一个上升沿很接近，其右列边界必定和一个下降沿很接近；对于记其左列边界和右列边界为和我们在上升沿集L^re中寻找一个接近上升沿满足条件：

            公式38

我们在下降沿集L^re中寻找一个接近下降沿满足条件

             公式39

对于不满足公式38的上升沿和不满足公式39的下降沿，将其从上升沿集L^re和下降沿集L^fe中去除掉。

通过第一区域以及筛选后的上升沿集L^re和下降沿集L^fe，定义第三区域为钙化斑块及其声影的图像所在区域，其中表示第一图像I中第i个第三区域，N_cs＝N_mrf表示这种区域的个数。任意一个第三区域的边界是由3条直线和1条曲线组成，如图26所示，线2601包围的区域2602即一个第三区域记最左边的像素和最右边的像素所在的列为和分别与的左边界、右边界所在的位置是一样的，那么的左边界是一条竖直直线，其列坐标是右边界是一条竖直直线，其列坐标是下边界是一条水平直线，其横坐标是H；上边界是最大值线中的一段，上边界的最左边点是最大值线的第个点，下边界的最右边点是最大值线的第个点。筛选后的上升沿集中，分别表示筛选后的上升沿集中第i个上升沿的最左边点和最右边点在第一图像I中的列坐标，筛选后的下降沿集中，分别表示第i个下降沿的最左边点和最右边点在第一图像I中的列坐标。则对于一个第三区域其左边界为筛选后的上升沿中的右边界为筛选后的下降沿中的

步骤604，过滤掉不符合预设约束条件的第三区域

过滤第三区域需要考虑以下5个约束条件：

约束条件1：所述第三区域的上边界对应的像素中，灰度值大于第一阈值的像素数量与所述第三区域的上边界对应的所有像素的数量的比大于第二阈值。优选地，第一阈值T₁＝200，第二阈值T₂＝0.3。不满足该约束条件1的第三区域要被去掉。

约束条件2：所述第三区域的上边界的最左端和最右端的像素的列坐标的差与行坐标的差的比的绝对值小于第三阈值。

具体地，对于一个第三区域其上边界的最左端和最右端的点的坐标是和给定第三阈值T₃＝1，有

$\left|\frac{d_r^i-d_l^i}{L_{d_r^i}^{\mathrm{mvl}}-L_{d_l^i}^{\mathrm{mvl}}}\right|<T_3$                            公式40

不满足该约束条件2的第三区域要被去掉。

约束条件3：所述第三区域的左边界和右边界的之间的距离小于第四阈值。

具体地，对于一个第三区域其左边界和右边界的之间距离需要小于第四阈值T₄＝300，即不满足该约束条件3的第三区域要被去掉。

约束条件4：所述第三区域中像素的灰度值均值为m_m，所述第三区域左邻的不包含钙化斑块及其声影的图像的区域中像素的灰度均值为m_l，所述第三区域右邻的不包含钙化斑块及其声影的图像的区域中像素的灰度均值为m_r，T₅为第五阈值，则 $\begin{array}{c}{m}_l>m_m,m_r>m_m\\ \frac{1}{2}(m_l+m_r)-m_m>T_5\end{array}.$

具体地，对于一个第三区域考虑与之间的区域，记作与之间的区域，记作的左边界是竖直直线，其列坐标是右边界是一条竖直直线，其列坐标是下边界是一条水平直线，其横坐标是H；上边界是最大值线中的一段，上边界的最左边点是最大值线的第个点，下边界的最右边点是最大值线的第个点；的左边界是竖直直线，其列坐标是右边界是一条竖直直线，其列坐标是下边界是一条水平直线，其横坐标是H；上边界是最大值线中的一段，上边界的最左边点是最大值线的第个点，下边界的最右边点是最大值线的第个点；计算中所有像素的像素均值分别是m_m，m_l，m_r。给定第五阈值T₅＝20，上述3个均值m_m，m_l，m_r必须满足：

m_l＞m_m,m_r＞m_m

                                             公式41

$\frac{1}{2}(m_l+m_r)-m_m>T_5$

不满足该约束条件4的第三区域要被去掉。

约束条件5：对于一个第三区域若以下三个条件B1、B2和B3只要有一个条件满足，该第三区域就保留。

条件B1，所述第三区域对应的上升沿和下降沿的所有像素中最大的行坐标与所述第三区域的上边界中的像素的灰度值均值的差大于第六阈值。

具体地，对于对应的上升沿和下降沿的所有点中，最大的行坐标是令m_ub表示上边界所有的像素的均值，则

$\max (p_{p_l^{{\mathrm{re}}_i}}^{\mathrm{mvl}},p_{p_r^{{\mathrm{fe}}_i}}^{\mathrm{mvl}})-m_{\mathrm{uv}}>T_6$                             公式42

优选地，第六阈值T₆＝0；

条件B2，用N_ab表示列坐标在和之间的，行坐标在1到之间的像素的个数，表示所述第三区域对应的上升沿最左边点的列坐标，表示所述第三区域对应的下降沿最右边点的列坐标，表示所述第三区域对应的上升沿和下降沿的所有像素中最大的行坐标，T₇表示第七阈值，则

$N_{\mathrm{ab}}/(p_r^{{\mathrm{fe}}_i}-p_l^{{\mathrm{re}}_i})>T_7$                                    公式43

其中优选地，第七阈值T₇＝0.3。

条件B3，用N_ub表示列坐标在和之间的，行坐标在到H之间的像素的个数，表示所述第三区域对应的上升沿最左边点的列坐标，表示所述第三区域对应的下降沿最右边点的列坐标，表示所述第三区域对应的上升沿和下降沿的所有像素中最大的行坐标，T₈表示第八阈值，则

$N_{\mathrm{ub}}/(p_r^{{\mathrm{fe}}_i}-p_l^{{\mathrm{re}}_i})<T_8$                     公式44

其中优选地，第八阈值T₈＝5。

如果式公式42满足，则该第三区域被保留下来；如果公式42不满足，而公式42满足，则第三区域被保留下来；如果公式42和公式43都不满足，而公式44满足，则该第三区域被保留下来；如果公式42、公式43和公式44都不满足，则该第三区域被丢弃。

步骤606，根据过滤后的第三区域R^cs提取钙化斑块及其声影的图像。

具体地，在一个实施例中，步骤606包括：

步骤702，根据过滤后的第三区域R^cs提取声影图像所在的区域。

由于最大值线是穿过所有的带声影的钙化斑块，而最大值线是任意第三区域的上边界，任意第三区域都含有一部分的钙化斑块图像和全部的声影图像。第三区域的数量等于钙化斑块的数量。钙化斑块和声影的图像之间存在边缘，可使用图搜索算法（比如深度优先遍历算法、广度优先遍历算法和双向广度优先遍历算法等）来提取钙化斑块和声影之间的边界，该边界也是声影的上边界。声影的下边界是一条水平直线，其行坐标是H；声影的左右边界各是一条竖直直线，其列坐标和的左右边界的列坐标相同。提取声影的结果如图27所示，区域2701就是声影图像所在区域。

步骤704，根据过滤后的第三区域R^cs提取钙化斑块所在的区域。

对于第i个第三区域其左边界和右边界的列坐标分别为和在和之间的每一列中（假设是第j列），可以找到第k_j行，该行的要求是，从第k_j行到第之间的灰度值全部是1。于是(k_j,j)就是第j列处的钙化斑块的上边界坐标。然后我们可以找到区域钙化斑块的左边界坐标是由组成，右边界坐标是由组成。提取钙化斑块的结果如图27所示，其中区域2702就是钙化斑块图像所在区域。

步骤706，提取钙化斑块及其声影的图像的边界。

如图27，在第一图像I中奖钙化斑块图像所在区域和声影图像所在区域分别用第一颜色和第二颜色标识出，并将该第一图像I从极坐标系变换到直角坐标系。然后使用Canny（坎尼）边缘算子分别提取第一颜色区域和第二颜色区域的边界。得到结果如图13所示，第一曲线2801包围的是钙化斑块图像，第二曲线2802包围的是声影图像。

上述钙化斑块及其声影的图像提取方法，使用图像处理技术自动从血管内超声图像中提取钙化斑块图像及钙化斑块的声影图像，无需人工干预，可自动判断出钙化斑块的位置，自动化程度高，提高了提取钙化斑块图像及钙化斑块的声影图像的效率。

如图8所示，在一个实施例中，提供了一种钙化斑块及其声影的超声图像提取系统，包括血管内超声图像获取模块8020、感兴趣区域确定模块8040、第一图像获取模块8060、第一区域获取模块8080和提取模块8100。

血管内超声图像获取模块8020，用于获取血管内超声图像。

血管内超声是将无创性的超声技术和有创性的心导管技术相结合，对心血管病变进行探测的一种方法。通过心导管将微型化的超声探头插入心血管腔内进行探测，再经电子成像系统可形成血管内超声图像，可表示心血管断面的形态和血流图形，如图17所示。

感兴趣区域确定模块8040用于根据所述血管内超声图像确定包含血管组织图像的感兴趣区域In1。

在血管内超声图像中，感兴趣区域In1为包括超声图像中血管内组织的所有信息的大致为圆形的区域。血管内超声图像中感兴趣区域In1之外的区域不包含有效信息，感兴趣区域确定模块8040用于剔除这部分无效区域，可避免干扰。在血管内超声图像中，该感兴趣区域In1为一圆形区域，如图18所示，圆形1801内的区域就是感兴趣区域In1。

第一图像获取模块8060用于以所述感兴趣区域In1的中心像素为直角坐标系的坐标原点，将所述感兴趣区域In1变换到极坐标系下，获得第一图像I。

第一图像获取模块8060用于以感兴趣区域In1的中心像素为直角坐标原点，将圆形的感兴趣区域In1从直角坐标系变换到极坐标系，得到矩形的第一图像I，如图20所示。令H表示第一图像I的行数，W表示第一图像I的列数。

第一区域获取模块8080用于根据所述第一图像I判断包括钙化斑块和声影的图像的区域作为第一区域R^mrf。

第一图像I是矩形图像，第一区域获取模块8080可用于通过扫描第一图像I的每一列，通过先验知识判断第一图像I每一列中是否含有钙化斑块图像和钙化斑块的声影图像，从而初步判断出包括钙化斑块和声影的图像的第一区域其中表示第i个第一区域，N_mrf表示第一区域在第一图像I中的个数，N_mrf可以等于1。需要说明的是，这里所说的图像的一列含有钙化斑块图像和钙化斑块的声影图像，是指该列包含钙化斑块和声影的一部分图像。

由于声影是因钙化斑块阻挡超声信号，从而在血管内超声图像上形成的在钙化斑块后端留下的阴影，因此在感兴趣区域In1中，包含钙化斑块与其对应的声影的区域大致为一扇形区域。而在极坐标系下的第一图像I中，该扇形区域对应矩形的第一区域R^mrf，根据矩形的第一区域R^mrf判断是否包含钙化斑块及其声影的图像，与根据扇形区域判断相比，可以降低计算复杂度。

提取模块8100用于根据所述第一区域R^mrf提取钙化斑块及其声影的图像。

确定了包含带声影的钙化斑块图像的第一区域R^mrf后，提取模块8100可用于使用边缘检测从第一区域R^mrf中提出钙化斑块图像和钙化斑块的声影图像。

上述钙化斑块及其声影的超声图像提取系统，在血管内超声图像中确定包含血管组织图像的感兴趣区域In1后，将感兴趣区域In1变换到极坐标系下，获得第一图像I。再将第一图像I中判断为包括钙化斑块和声影的图像的第一区域R^mrf，从第一区域R^mrf中提出钙化斑块图像和钙化斑块的声影图像。实现了钙化斑块图像及钙化斑块的声影图像的自动提取，极大地提高了提取钙化斑块及其声影的图像的效率。

如图9所示，在一个实施例中，该所述系统还包括聚类模块8010、聚类图像获取模块8030和第二图像获取模块8050。

聚类模块8010，用于对所述感兴趣区域In1的像素进行聚类。

具体地，如图10所示，在一个实施例中，所述聚类模块包括概率分布描述模块8012、参数求解模块8014和聚类执行模块8016。

概率分布描述模块8012用于使用瑞利混合分布模型描述感兴趣区域In1的每个像素的概率分布。

具体地，感兴趣区域In1中的每个像素的概率分布如公式1所示：

则所述感兴趣区域的每个像素的概率分布为：

$p\left(y_i\right)=\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_{j}p\left(y_i\right|{\mathrm{\&theta;}}_j)$                              ；公式1

$p\left(y_i\right|{\mathrm{\&theta;}}_j)=\frac{y_i-a_j}{{\mathrm{\&sigma;}}_j^2}\exp (-\frac{{(y_i-a_j)}^2}{{2\mathrm{\&sigma;}}_j^2})$

其中，Y＝{y₁,...,y_N}表示感兴趣区域In1的像素集，y_i表示感兴趣区域的第i个像素的灰度值，θ_j＝{a_j,σ_j}是瑞利分布的参数，具体地，a_j表示第j个瑞利分布在横轴上的平移量，σ_j表示瑞利分布的众数（mode，众数是指一组数据中出现次数最多的数值）；π_j表示瑞利混合分布中每个分量的权重，K表示瑞利混合分布中混合分量的个数，p(y_i)表示第i个像素的混合概率，p(y_i|θ_j)表示第i个像素属于第j类的概率。

对于公式1，瑞利混合分布的似然函数为：

$L\left(\mathrm{\&Theta;}\right)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log \left(\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_{j}p\left(y_i\right|{\mathrm{\&theta;}}_j)\right)$                               公式2

其中，N表示感兴趣区域中的像素总数，Θ表示参数集；

然后，聚类模块包括概率分布描述模块8012用于定义一个函数表示感兴趣区域In1中的第i个像素对瑞利混合分布中的第j个分量的权重：

${\mathrm{\&xi;}}_j\left(y_i\right)=\frac{{\stackrel{\&OverBar;}{y}}_i-c_j}{b_j^2}\exp (-\frac{{({\stackrel{\&OverBar;}{y}}_i-c_j)}^2}{{2b}_j^2})$                      公式3

其中c_j和b_j表示参数，具体地，c_j表示所满足的瑞利分布的横轴偏移量，b_j表示所满足的瑞利分布的众数，表示第i个像素的8邻域的均值，即可以写成

${\stackrel{\&OverBar;}{y}}_i=\frac{1}{9}\underset{m=-1}{\overset{1}{\mathrm{\&Sigma;}}}\underset{n=-1}{\overset{1}{\mathrm{\&Sigma;}}}y(u+m,v+n)$                             公式4

其中，坐标为(u,v)的像素为中心点，调整步长m、n可取遍该像素领域的各个点，y(u+m,v+n)表示血管内超声图像中坐标为(u+m,v+n)的点的灰度值。

举例说明，假定有一个3×3的窗口，(u,v)就是窗口在血管内超声图像上的中心坐标，假设该坐标为(100,100)，即u=100,v=100。而m、n就是从(-l,l)的连续变化量，表示对窗口中的数据进行遍历。在上面假设中l=1，那么m和n都是从-1到+1。这样，公式4就是求(99,99)、(99,100)、(99,101)、(100,99)、(100,100)、(100,101)、(101,99)、(101,100)和(101,101)这9个点的灰度值的和的平均值。

对于第i个像素的邻域N_i，定义一个对瑞利混合分布的第j个分量为对象的权重函数：

                公式5

其中M是邻域N_i中像素的个数，α是一个控制变量用来控制公式5中值的大小。这里的邻域N_i是8-邻域。

最后聚类模块包括概率分布描述模块8012用于定义一个新的先验概率π_ij，表示第i个像素的邻域N_i的权重函数在第j个瑞利分布中所占有权重。

                公式6

参数求解模块8014用于使用最大期望算法求解所述瑞利混合模型的参数。

具体地，如图11所示，所述参数求解模块8014包括似然函数最大化模块8014a、参数集初始化模块8014b和参数集计算模块8014c。

似然函数最大化模块8014a用于为了估计出公式1中混合模型的参数θ＝{a_j,σ_j},j＝1,...,K，最大化公式2中的似然函数，即

${\mathrm{\&Theta;}}^{*}=\arg \underset{\mathrm{\&Theta;}}{\max }L\left(\mathrm{\&Theta;}\right)$                           公式7

则使用EM算法得到的目标函数为：

$Q(\mathrm{\&Theta;},{\mathrm{\&Theta;}}^{\left(t\right)})=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\log \left({\mathrm{\&pi;}}_{j}P\left({\mathrm{\&theta;}}_j\right|y_i;{\mathrm{\&Theta;}}^{\left(t\right)})\right)+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\log \left(P\left(y_i\right|{\mathrm{\&theta;}}_j;\mathrm{\&Theta;})\right)P\left({\mathrm{\&theta;}}_j\right|y_i;{\mathrm{\&Theta;}}^{\left(t\right)});$        公式8

其中，θ_j＝{a_j,σ_j,c_j,b_j,α}是参数向量，具体地，a_j表示第j个瑞利分布在横轴上的平移量，σ_j表示瑞利分布的众数，c_j和b_j是计算参数，具体地，c_j表示所满足的瑞利分布的横轴偏移量，b_j表示所满足的瑞利分布的众数，α表示控制变量；π_j是先验概率，P(θ_j|y_i;Θ^(t))是后验概率，P(y_i|θ_j;Θ)是类条件概率密度，Θ^(t)表示第t次迭代中已知的参数集，Θ表示第t次迭代中未知的参数集。

参数集初始化模块8014b用于初始化参数集Θ为Θ⁽⁰⁾。

参数集初始化模块8014b用于初始化参数集Θ，记作Θ⁽⁰⁾。具体地，参数集初始化模块8014b用于给定分类个数K＝5和参数α＝10^-7，使用K均值算法（一种聚类算法，算法步骤为：1、输入数据集合和类别数K(用户指定)；2、随机分配类别中心点的位置；3、将每个点放在离他最近的类别中心点所在的集合；4、移动类别中心点到它所在集合的中心；5、转到第3步，直到收敛。）计算出感兴趣区域In1每一类的均值 $m^{\left(0\right)}=[{\mathrm{\&mu;}}_1^{\left(0\right)},{\mathrm{\&mu;}}_2^{\left(0\right)},...,{\mathrm{\&mu;}}_K^{\left(0\right)}].$

令a_j表示上一步的K均值分类后，第j类的所有像素中的最小灰度值。

从下式公式9中得到 ${\mathrm{\&sigma;}}^{\left(0\right)}=[{\mathrm{\&sigma;}}_1^{\left(0\right)},{\mathrm{\&sigma;}}_2^{\left(0\right)},...,{\mathrm{\&sigma;}}_K^{\left(0\right)}],$

${\mathrm{\&sigma;}}^{\left(0\right)}=(m^{\left(0\right)}-a^{\left(0\right)})\sqrt{\frac{2}{\mathrm{\&pi;}}}$                   公式9

令c⁽⁰⁾＝a⁽⁰⁾，b⁽⁰⁾＝σ⁽⁰⁾。

参数集计算模块8014c用于根据初始化的参数集Θ⁽⁰⁾计算参数向量，并使用最速下降算法更新所述参数向量，直到最大期望算法收敛，获得最终参数集Θ^*。

参数求解模块8014还用于根据初始化的参数集Θ⁽⁰⁾计算参数向量，并使用最速下降算法更新所述参数向量，直到最大期望算法收敛，获得最终参数集Θ^*。

参数集Θ的计算使用到了最速下降法(steepest-descent method)：

${\mathrm{\&Theta;}}^{(t+1)}={\mathrm{\&Theta;}}^{\left(t\right)}-\mathrm{\&eta;}\frac{\&PartialD;Q(\mathrm{\&Theta;},{\mathrm{\&Theta;}}^{\left(t\right)})}{\&PartialD;\mathrm{\&Theta;}}$                       公式10

参数集计算模块8014c用于通过公式1和公式6，由贝叶斯定理得到像素i在模型中的后验概率为公式11，表示第i个像素属于第j类的概率：

$P^{\left(t\right)}\left({\mathrm{\&theta;}}_j\right|y_i)=\frac{{\mathrm{\&pi;}}_{\mathrm{ij}}^{\left(t\right)}{P}^{\left(t\right)}\left(y_i\right|{\mathrm{\&theta;}}_j)}{\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_{\mathrm{ik}}^{\left(t\right)}{P}^{\left(t\right)}\left(y_i\right|{\mathrm{\&theta;}}_k)}$                              公式11

参数集计算模块8014c还用于计算参数向量为：

$\frac{\&PartialD;Q}{\&PartialD;\mathrm{\&Theta;}}={\left[\begin{array}{ccccc}\frac{\&PartialD;Q}{\&PartialD;a}& \frac{\&PartialD;Q}{\&PartialD;\mathrm{\&sigma;}}& \frac{\&PartialD;Q}{\&PartialD;c}& \frac{\&PartialD;Q}{\&PartialD;b}& \frac{\&PartialD;Q}{\&PartialD;\mathrm{\&alpha;}}\end{array}\right]}^T$                 公式12

且为

$\frac{\&PartialD;Q}{\&PartialD;a_j}=-\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}P\left({\mathrm{\&theta;}}_j\right|y_i)(\frac{y_i-a_j}{{\mathrm{\&sigma;}}_j^2}-\frac{1}{y_i-a_j})$                      公式13

$\frac{\&PartialD;Q}{\&PartialD;{\mathrm{\&sigma;}}_j}=-\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}P\left({\mathrm{\&theta;}}_j\right|y_i)(-\frac{2}{{\mathrm{\&sigma;}}_j}+\frac{{(y_i-a_j)}^2}{{\mathrm{\&sigma;}}_j^3})$                   公式14

公式15

公式16

公式17

然后参数集计算模块8014c可用于根据公式10更新参数向量。

参数集计算模块8014c用于当公式10中的参数向量不再变换时，EM算法收敛，记此时计算出的参数集为Θ^*；否则，令Θ^(t)＝Θ^(t+1)，然后参数集计算模块8014c用于通过公式1和公式6，由贝叶斯定理可以得到像素i在模型中的后验概率，继续进行计算。

聚类执行模块8016用于使用最大后验概率准则对所述感兴趣区域In1的像素进行聚类。

得到最终的参数集Θ^*后，聚类执行模块8016用于使用最大后验概率准则去对感兴趣区域In1中的每个像素进行聚类；聚类个数等于瑞利混合分布的分量的个数κ。聚类执行模块8016用于把第i个像素归为第j类，如果

            公式18

其中，表示正整数集，K表示聚类个数。

聚类图像获取模块8030用于将聚类结果中属于同一聚类簇的像素的灰度值置为相同的值，且属于不同聚类簇的像素的灰度值不同，获得聚类图像。

具体地，感兴趣区域In1被聚类为K类，每一类称为一个聚类簇。聚类图像获取模块8030将同一聚类簇的像素的灰度值置为相同的值，不同聚类簇的像素的灰度值置为不同的值，且不同聚类簇的像素的灰度值所置的值可与聚类簇中像素的灰度值均值正相关。进一步地，聚类图像获取模块8030可用于若聚类结果中像素属于第k类，可将该像素的灰度值设置为k。由于最后获得的聚类图像中就只有K种灰度值，达到了降维的效果，可以降低计算复杂度。获得的聚类图像如图19所示。

第二图像获取模块8050用于以所述聚类图像的中心像素为直角坐标系的坐标原点，将所述聚类图像变换到极坐标系，获得第二图像I_R。

第二图像获取模块8050用于将把感兴趣区域In1聚类降维后获得的大致为圆形的聚类图像变换到极坐标系，获得的矩形的第二图像I_R，如图21所示，便于后续计算。

本实施例中对感兴趣区域In1进行聚类降维，可降低计算复杂度。

如图9所示，在一个实施例中，该所述系统还包括最大值线获取模块8070。

最大值线获取模块8070用于确定所述第一图像每列的像素中灰度值最大的像素，获得最大值线L^mvl；

具体地，令C＝{1,2,...,W}表示第一图像I的列编号的集合，对于任意i∈C，最大值线获取模块8070用于把I的第i列的具有最大灰度值的像素的行数记作于是对于第一图像I的每一列，可以得到记这里L^mvl即为最大值线。如附图20中线段2201所示。第一图像I中行坐标从左到右递增，列坐标从上到下递增。

如图9所示，在一个实施例中，所述第一区域获取模块8080包括特征提取模块8082、分类模块8084和区域确定模块8086。

特征提取模块8082用于提取所述第一图像I中每列的特征。

第一图像I本身包括大量的像素，为了方便计算，特征提取模块8082用于对第一图像I中每列提取特征，降低计算复杂度。

如图12所示，在一个实施例中，所述特征提取模块8082包括第二区域确定模块8082a、第三图像生成模块8082b、第二区域上边界确定模块8082c特征参数获取模块8082d和特征计算模块8082e。

第二区域确定模块8082a用于计算所述第二图像I_R中每个聚类簇对应的感兴趣区域In1中的像素的灰度值均值，找到灰度值均值最小的预设个数的聚类簇，将所述找到的聚类簇中的像素构成的区域作为第二区域R^D。

由于第二图像I_R是有K类，第二区域确定模块8082a计算出每类的灰度值均值，找到灰度值均值最小的两类，将这两类中所有像素的位置定义为第二区域R^D。当预设个数为2时，能够达到最好的效果，识别正确率显著提高。

第三图像生成模块8082b用于生成与所述第二图像I_R尺寸相同的第三图像I_D，将所述第三图像I_D中属于第二区域R^D的像素的灰度值和不属于第二区域R^D的像素的灰度值分别设置为不同的值。

第三图像生成模块8082b可用于定义第三图像I_D，其大小与第二图像I_R相同，第三图像I_D中属于第二区域R^D的像素的灰度值设为1（如图23中白色区域2301），不属于第二区域R^D的像素的灰度值设为0。

第二区域上边界确定模块8082c用于在所述第三图像I_D中确定第二区域R^D的上边界。

第二区域上边界确定模块8082c可用于从上向下扫描第三图像I_D的第i列(i＝1,2,...,W)，找到属于第二区域R^D的第一个像素（灰度值为1的第一个像素），把该像素的行坐标记作于是，我们可以在第三图像I_D中定义第二区域R^D的上边界，记为 $L^{\mathrm{ubrd}}=\{p_1^{\mathrm{ubrd}},p_2^{\mathrm{ubrd}},...,p_W^{\mathrm{ubrd}}\}.$

特征参数获取模块8082d用于计算所述第三图像I_D每列的底部与所述第二区域R^D的上边界L^ubrd的距离、所述第二区域R^D的上边界与所述最大值线L^mvl的距离以及所述第一图像I中每列对应的所述第二区域的上边界与所述最大值线L^mvl之间的像素的平均灰度值。

特征参数获取模块8082d可用于计算三个自定义的特征参数F^h,F^d,F^v，其中，表示第三图像I_D每列的底部与所述第二区域R^D的上边界L^ubrd的距离。

$F_i^h=H-p_i^{\mathrm{ubrd}},i=\mathrm{1,2},...,W$                           公式20

表示第二区域R^D的上边界L^ubrd与最大值线L^mvl的距离，

$F_i^d=p_i^{\mathrm{ubrd}}-p_i^{\mathrm{mvl}},i=\mathrm{1,2},...,W$                           公式21

表示第一图像I中每列对应的第二区域R^D的上边界L^ubrd与最大值线L^mvl之间的像素的平均灰度值。

$F_i^v=\frac{1}{H-p_i^{\mathrm{ubrd}}+1}\underset{j=p_i^{\mathrm{ubrd}}}{\overset{p_i^{\mathrm{mvl}}}{\mathrm{\&Sigma;}}}I(p_i^{\mathrm{ubrd}},i),i=\mathrm{1,2},...,W$                   公式22

特征计算模块8082e，用于根据所述第三图像I_D每列的底部与所述第二区域R^D的上边界L^ubrd的距离、所述第二区域R^D的上边界与所述最大值线L^mvl的距离以及所述第一图像I中每列对应的所述第二区域的上边界与所述最大值线L^mvl之间的像素的平均灰度值确定所述第一图像每列的特征。

特征计算模块8082e可以用于利用上面三个参数F^h,F^d,F^v，定义F＝{F₁,F₂,...,F_W}为第一图像I每列的特征：

$F=\frac{h_1{F}^h+h_2{F}^d+h_3{F}_v}{255}$                                   公式23

F_i＝{F₁,F₂,...,F_W}，其中h₁,h₂,h₃是权重,优选的,权重h₁,h₂,h₃的值分别为h₁＝5,h₂＝-0.5,h₃＝-1。当h₁＝5,h₂＝-0.5,h₃＝-1时，第一图像I每列的特征能够很好地反映该列是否含有钙化斑块图像和声影的图像，性能最优。

分类模块8084用于根据所述第一图像I中每列的特征将所述第一图像I中的各列分为含有钙化斑块和声影的图像的列以及不含有钙化斑块和声影的图像的列两类。

如图13所示，在一个实施例中，所述分类模块8084包括置信计算模块8084a和分类执行模块8084b。

置信计算模块8084a用于根据所述第一图像I每列的特征，使用置信传播算法在马尔可夫随机场上计算出第一图像每列含有钙化斑块和声影的图像的置信。

如图14所示，所述置信计算模块8084a包括位置集和状态集定义模块8084a1、初始化模块8084a2、迭代模块8084a3和计算模块8084a4。

位置集和状态集定义模块8084a1用于定义马尔可夫随机场的位置集S和状态集L分别为

S＝{1,2,...,W}                          公式19

L＝{-1,+1}

其中，S＝{1,2,...,W}表示该列的位置；L＝{-1,+1}表示该列的状态，若第一图像I的一列的状态为“+1”，表示该列含有钙化斑块及其声影的图像；如果第一图像I的一列的状态为“-1”，则表示该列不含有钙化斑块及其声影的图像。

马尔可夫随机场可由概率图模型表示，如图22所示；图中节点χ₁,...,χ_W是观测变量，表示图像I每列的特征，节点z₁,...,z_W是隐藏变量，表示图像I每列的状态。

初始化模块8084a2对于马尔可夫随机场上的每个隐藏变量z_i，令所有的初始状态的概率满足均匀分布，即隐藏变量z_i的边缘概率为：

P(z_i＝-1)＝P(z_i＝1)＝0.5                          公式24

且隐藏变量z_i的置信初始化为：

$b_i^{\left(0\right)}(z_i=-1)=b_i^{\left(0\right)}(z_i=1)=0.5$                          公式25

局部信息φ(z_i,χ_i)由公式22和公式23决定：

φ_i(z_i,χ_i)＝F_i                          公式26

F_i是第一图像第i列的特征，χ_i表示观测变量χ₁,...,χ_W。

相容函数ψ(z_i,z_j)是一个2×2的矩阵

${\mathrm{\&psi;}}_{i,j}(z_i,z_j)=\left[\begin{array}{cc}0.8& 0.2\\ 0.2& 0.8\end{array}\right]$                       公式27

当节点z_i和z_j是互为邻域，那么从z_i传递到z_j的信息初始化为：

m_i,j(z_j)＝1                                    公式28

迭代模块8084a3用于在第t次迭代中，计算从z_i传递到z_j的信息

$m_{i,j}^{\left(t\right)}\left(z_j\right)=\underset{z_i}{\mathrm{\&Sigma;}}{\mathrm{\&phi;}}_i(z_i,{\mathrm{\&chi;}}_i){\mathrm{\&psi;}}_{i,j}(z_i,z_j){\mathrm{\&Pi;}}_{k\&Element;N\left(i\right)\backslash j}{m}_{k,i}^{\left(t\right)}\left(z_i\right)$                     公式29

并计算节点z_i的置信

$b_i^{\left(t\right)}\left(z_i\right)=k{\mathrm{\&phi;}}_i\left(z_i\right)\underset{j\&Element;N\left(i\right)}{\mathrm{\&Pi;}}{m}_{j,i}^{\left(t\right)}\left(z_i\right)$                          公式30

其中N(i)是节点z_i的邻域。

计算模块8084a4用于当在t+1次迭代中，满足公式31时，迭代算法收敛。

$\frac{1}{W}|\underset{i=1}{\overset{W}{\mathrm{\&Sigma;}}}\left(b^{(i+1)}\left(z_i\right)\right)-\underset{i=1}{\overset{W}{\mathrm{\&Sigma;}}}\left(b^{\left(t\right)}\left(z_i\right)\right)|<\mathrm{\&epsiv;}$                        公式31

其中ε为预设定值。

则每个隐藏变量的置信为，即每列含有钙化斑块和声影的图像的置信为：

b^*(z₁)＝b^(t+1)(z₁),b^*(z₂)＝b^(t+1)(z₂),...,b^*(z_W)＝b^(t+1)(z_W)        公式32

其中，W表示第一图像的列数。

最终第一图像I每列含有钙化斑块和声影的图像的置信如图24所示，图24中的线2401表示每列含有钙化斑块和声影的图像的置信，其中横轴表示第一图像I的列编号，纵轴表示置信，线2402表示置信等于0.5。

分类执行模块8084b用于根据所述第一图像I每列的置信将所述第一图像中的各列分为含有钙化斑块和声影的图像的列和不含有钙化斑块和声影的图像的列两类。

具体地，分类执行模块8084b用于对第一图像I的列进行分类使用的是最大后验准则。具体地，对于第一图像I的第i列的状态z_i，如果满足：

b^*(z_i＝-1)＞b^*(z_i＝+1)                      公式33

那么将第i列分到类“-1”。如果满足

b^*(z_i＝-1)≤b^*(z_i＝+1)                       公式34

那么将第i列分到类“+1”。其中，若该列属于类“+1”表示该列含有钙化斑块和声影的图像，若该列不属于类“-1”则表示该列不含有钙化斑块和声影的图像。

区域确定模块8086用于将所述第一图像I中连续的含有钙化斑块和声影的图像的列构成的区域判断为包括钙化斑块和声影的图像的第一区域R^mrf。

具体地，区域确定模块8086可用于对第一图像I从左向右每列依次扫描，把连续都是类“+1”的列当作同一个区域，记第一区域其中表示第i个第一区域，N_mrf表示这种区域在第一图像I中的个数，N_mrf可以等于1，结果见图25，其中2501表示第一图像I各列的置信，线2502表示置信为0.5的界限。对于任意一个第一区域其最左边列称作左列边界，其最右边的列称作右列边界，其中置信为1的列属于类“+1”，置信为0的列属于类“-1”；类“+1”表示含有钙化斑块和声影的图像，类“-1”表示不含有钙化斑块和声影的图像。

如图15所示，在一个实施例中，所述提取模块8100包括第三区域确定模块8102、过滤模块8104和提取执行模块8106。

第三区域确定模块8102用于根据所述第一区域R^mrf和所述第一图像I的最大值线L^mvl确定第三区域R^cs。

第一区域R^mrf是初步判断包括钙化斑块和声影的图像的区域，但这种判断并不精确，为了精确判断出钙化斑块和声影图像所在位置，需要先根据所述第一区域R^mrf和所述第一图像I的最大值线L^mvl定义第三区域R^cs，再根据第三区域R^cs精确判断钙化斑块及其声影的图像所在的位置。

在一个实施例中该系统还包括上升沿和下降沿处理模块（图中未示出），用于在第一图像I中定义最大值线的上升沿和下降沿。最大值线的上升沿是指最大值线上的一个连续线段；最大值线上的所有上升沿记作其中N_re表示上升沿的个数，表示最大值线上的第i个上升沿，且可以表示成

$L_i^{\mathrm{re}}=\{P_l^{{\mathrm{re}}_i},p_r^{{\mathrm{re}}_i}\}$                                 公式35

其中，分别表示第i个上升沿的最左边点和最右边点在第一图像I中的列坐标；对于i＝1,2,...,N_re，上升沿有以下性质：

$p_{p_l^{{\mathrm{re}}_i}}^{\mathrm{mvl}}\&GreaterEqual;p_{p_l^{{\mathrm{re}}_i}+1}^{\mathrm{mvl}}\&GreaterEqual;...\&GreaterEqual;p_{p_r^{{\mathrm{re}}_i}}^{\mathrm{mvl}}$                                公式36

最大值线的下降沿可表示成 $L^{\mathrm{fe}}=\{L_1^{\mathrm{fe}},L_2^{\mathrm{fe}},...,L_{N_{\mathrm{fe}}}^{\mathrm{fe}}\},L_i^{\mathrm{fe}}=\{p_l^{{\mathrm{fe}}_i},p_r^{{\mathrm{fe}}_i}\},$ 其中N_fe表示下降沿的个数，分别表示第i个下降沿的最左边点和最右边点在第一图像I中的列坐标。对于i＝1,2,...,N_fe，下降沿有以下性质：

$p_{p_l^{{\mathrm{fe}}_i}}^{\mathrm{mvl}}\&le;p_{p_l^{{\mathrm{fe}}_i}+1}^{\mathrm{mvl}}\&le;...\&le;p_{p_r^{{\mathrm{fe}}_i}}^{\mathrm{mvl}}$                                  公式37

然后上升沿和下降沿处理模块还用于筛选上升沿和下降沿。对于某个可能包含钙化斑块图像和钙化斑块的声影图像的区域其左列边界必定和一个上升沿很接近，其右列边界必定和一个下降沿很接近；对于记其左列边界和右列边界为和我们在上升沿集L^re中寻找一个接近上升沿满足条件：

             公式38

我们在下降沿集L^re中寻找一个接近下降沿满足条件

             公式39

对于不满足公式38的上升沿和不满足公式39的下降沿，将其从上升沿集L^re和下降沿集L^fe中去除掉。

第三区域确定模块8102用于通过第一区域定义第三区域为钙化斑块及其声影的图像所在区域，其中表示第一图像I中第i个第三区域，N_cs＝N_mrf表示这种区域的个数。任意一个第三区域的边界是由3条直线和1条曲线组成，如图26所示，线2601包围的区域2602即一个第三区域记最左边的像素和最右边的像素所在的列为和，分别与的左边界、右边界所在的位置是一样的，那么的左边界是一条竖直直线，其列坐标是右边界是一条竖直直线，其列坐标是下边界是一条水平直线，其横坐标是H；上边界是最大值线中的一段，上边界的最左边点是最大值线的第个点，下边界的最右边点是最大值线的第个点。

过滤模块8104用于过滤掉不符合预设约束条件的第三区域。

过滤模块8104用于过滤第三区域需要考虑以下5个约束条件：

约束条件1：所述第三区域的上边界对应的像素中，灰度值大于第一阈值的像素数量与所述第三区域的上边界对应的所有像素的数量的比大于第二阈值。优选地，第一阈值T₁＝200，第二阈值T₂＝0.3。过滤模块8104用于过滤掉不满足该约束条件1的第三区域。

约束条件2：所述第三区域的上边界的最左端和最右端的像素的列坐标的差与行坐标的差的比的绝对值小于第三阈值。

具体地，对于一个第三区域其上边界的最左端和最右端的点的坐标是和给定第三阈值T₃＝1，有

$\left|\frac{d_r^i-d_l^i}{L_{d_r^i}^{\mathrm{mvl}}-L_{d_l^i}^{\mathrm{mvl}}}\right|<T_3$                          公式40

过滤模块8104用于过滤掉不满足该约束条件2的第三区域。

约束条件3：所述第三区域的左边界和右边界的之间的距离小于第四阈值。

具体地，对于一个第三区域其左边界和右边界的之间距离需要小于第四阈值T₄＝300，即过滤模块8104用于过滤掉不满足该约束条件3的第三区域。

约束条件4：所述第三区域中像素的灰度值均值为m_m，所述第三区域左邻的不包含钙化斑块及其声影的图像的区域中像素的灰度均值为m_l，所述第三区域右邻的不包含钙化斑块及其声影的图像的区域中像素的灰度均值为m_r，T₅为第五阈值，则 $\begin{array}{c}{m}_l>m_m,m_r>m_m\\ \frac{1}{2}(m_l+m_r)-m_m>T_5\end{array}.$

具体地，对于一个第三区域考虑与之间的区域，记作与之间的区域，记作的左边界是竖直直线，其列坐标是右边界是一条竖直直线，其列坐标是下边界是一条水平直线，其横坐标是H；上边界是最大值线中的一段，上边界的最左边点是最大值线的第个点，下边界的最右边点是最大值线的第个点；的左边界是竖直直线，其列坐标是右边界是一条竖直直线，其列坐标是下边界是一条水平直线，其横坐标是H；上边界是最大值线中的一段，上边界的最左边点是最大值线的第个点，下边界的最右边点是最大值线的第个点；计算中所有像素的像素均值分别是m_m，m_l，m_r。给定第五阈值T₅＝20，上述3个均值m_m，m_l，m_r必须满足：

m_l＞m_m,m_r＞m_m

                                           公式41

$\frac{1}{2}(m_l+m_r)-m_m>T_5$

过滤模块8104用于过滤掉不满足该约束条件4的第三区域。

约束条件5：过滤模块8104用于对于一个第三区域若以下三个条件B1、B2和B3只要有一个条件满足，该第三区域就保留。

条件B1，所述第三区域对应的上升沿和下降沿的所有像素中最大的行坐标与所述第三区域的上边界中的像素的灰度值均值的差大于第六阈值。

具体地，对于对应的上升沿和下降沿的所有点中，最大的行坐标是令m_ub表示上边界所有的像素的均值，则

$\max (p_{p_l^{{\mathrm{re}}_i}}^{\mathrm{mvl}},p_{p_r^{{\mathrm{fe}}_i}}^{\mathrm{mvl}})-m_{\mathrm{ub}}>T_6$                                     公式42

优选地，第六阈值T₆＝0；

条件B2，用N_ab表示列坐标在和之间的，行坐标在1到之间的像素的个数，表示所述第三区域对应的上升沿最左边点的列坐标，表示所述第三区域对应的下降沿最右边点的列坐标，表示所述第三区域对应的上升沿和下降沿的所有像素中最大的行坐标，T₇表示第七阈值，则

$N_{\mathrm{ab}}/(p_r^{{\mathrm{fe}}_i}-p_l^{{\mathrm{re}}_i})>T_7$                                 公式43

其中优选地，第七阈值T₇＝0.3。

条件B3，用N_ub表示列坐标在和之间的，行坐标在到H之间的像素的个数，表示所述第三区域对应的上升沿最左边点的列坐标，表示所述第三区域对应的下降沿最右边点的列坐标，表示所述第三区域对应的上升沿和下降沿的所有像素中最大的行坐标，T₈表示第八阈值，则

$N_{\mathrm{ub}}/(p_r^{{\mathrm{fe}}_i}-p_l^{{\mathrm{re}}_i})<T_8$                                 公式44

其中优选地，第八阈值T₈＝5。

如果式公式42满足，则该第三区域被保留下来；如果公式42不满足，而公式42满足，则第三区域被保留下来；如果公式42和公式43都不满足，而公式44满足，则该第三区域被保留下来；如果公式42、公式43和公式44都不满足，则该第三区域被丢弃。

提取执行模块8106用于根据过滤后的第三区域提取钙化斑块及其声影的图像。

如图16所示，在一个实施例中，所述提取执行模块8106包括声影图像提取模块8106a、钙化斑块图像提取模块8106b和边界提取模块8106c。

声影图像提取模块8106a用于根据过滤后的第三区域R^cs提取声影图像所在的区域。

由于最大值线是穿过所有的带声影的钙化斑块，而最大值线是任意第三区域的上边界，任意第三区域都含有一部分的钙化斑块图像和全部的声影图像。第三区域的数量等于钙化斑块的数量。钙化斑块和声影的图像之间存在边缘，声影图像提取模块8106a可用于使用图搜索算法（比如深度优先遍历算法、广度优先遍历算法和双向广度优先遍历算法等）来提取钙化斑块和声影之间的边界，该边界也是声影的上边界。声影的下边界是一条水平直线，其行坐标是H；声影的左右边界各是一条竖直直线，其列坐标和的左右边界的列坐标相同。提取声影的结果如图27所示，区域2701就是声影图像所在区域。

钙化斑块图像提取模块8106b用于根据过滤后的第三区域提取钙化斑块所在的区域。

钙化斑块图像提取模块8106b用于对于第i个第三区域其左边界和右边界的列坐标分别为和钙化斑块图像提取模块8106b可用于在和之间的每一列中（假设是第j列），可以找到第k_j行，该行的要求是，从第k_j行到第之间的灰度值全部是1。于是(k_j,j)就是第j列处的钙化斑块的上边界坐标。然后我们可以找到区域钙化斑块的左边界坐标是由组成，右边界坐标是由组成。提取钙化斑块的结果如图27所示，其中区域2702就是钙化斑块图像所在区域。

边界提取模块8106c用于提取钙化斑块及其声影的图像的边界。

如图27，边界提取模块8106c用于在第一图像I中奖钙化斑块图像所在区域和声影图像所在区域分别用第一颜色和第二颜色标识出，并将该第一图像I从极坐标系变换到直角坐标系。然后边界提取模块8106c用于使用Canny（坎尼）边缘算子分别提取第一颜色区域和第二颜色区域的边界。得到结果如图13所示，第一曲线2801包围的是钙化斑块图像，第二曲线2802包围的是声影图像。

上述钙化斑块及其声影的超声图像提取系统，使用图像处理技术自动从血管内超声图像中提取钙化斑块图像及钙化斑块的声影图像，无需人工干预，可自动判断出钙化斑块的位置，自动化程度高，提高了提取钙化斑块图像及钙化斑块的声影图像的效率。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体（Read-Only Memory，ROM）或随机存储记忆体（Random Access Memory，RAM）等。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (17)

1.一种钙化斑块及其声影的超声图像提取系统，其特征在于，所述系统包括： 

血管内超声图像获取模块，用于获取血管内超声图像； 

感兴趣区域确定模块，用于根据所述血管内超声图像确定包含血管组织图像的感兴趣区域； 

第一图像获取模块，用于以所述感兴趣区域的中心像素为直角坐标系的坐标原点，将所述感兴趣区域变换到极坐标系下，获得第一图像； 

第一区域获取模块，用于根据所述第一图像判断包括钙化斑块和声影的图像的区域作为第一区域； 

提取模块，用于根据所述第一区域提取钙化斑块及其声影的图像。 

2.根据权利要求1所述的系统，其特征在于，所述第一区域获取模块包括： 

特征提取模块，用于提取所述第一图像中每列的特征； 

分类模块，用于根据所述第一图像中每列的特征将所述第一图像中的各列分为含有钙化斑块和声影的图像的列以及不含有钙化斑块和声影的图像的列两类； 

区域确定模块，用于将所述第一图像中连续的含有钙化斑块和声影的图像的列构成的区域判断为包括钙化斑块和声影的图像的第一区域。 

3.根据权利要求2所述的系统，其特征在于，所述系统还包括： 

聚类模块，用于对所述感兴趣区域的像素进行聚类； 

聚类图像获取模块，用于将聚类结果中属于同一聚类簇的像素的灰度值置为相同的值，且属于不同聚类簇的像素的灰度值不同，获得聚类图像； 

第二图像获取模块，用于以所述聚类图像的中心像素为直角坐标系的坐标原点，将所述聚类图像变换到极坐标系，获得第二图像。 

4.根据权利要求3所述的系统，其特征在于，所述聚类模块包括： 

概率分布描述模块，用于使用瑞利混合分布模型描述感兴趣区域的每个像素的概率分布； 

参数求解模块，用于使用最大期望算法求解所述瑞利混合模型的参数； 

聚类执行模块，用于使用最大后验概率准则对所述感兴趣区域的像素进行聚类。 

5.根据权利要求4所述的系统，其特征在于，所述感兴趣区域的每个像素的概率分布为： 

其中，y_i表示感兴趣区域的第i个像素的灰度值，θ_j＝{a_j,σ_j}是瑞利分布的参数，具体地，a_j表示第j个瑞利分布在横轴上的平移量，σ_j表示瑞利分布的众数；π_j表示瑞利混合分布中每个分量的权重，K表示瑞利混合分布中混合分量的个数，p(y_i)表示第i个像素的混合概率，p(y_i|θ_j)表示第i个像素属于第j类的概率。 

6.根据权利要求5所述的系统，其特征在于，所述瑞利混合分布的似然函数为： 

其中，N表示感兴趣区域中的像素总数，Θ表示参数集； 

所述参数求解模块包括： 

似然函数最大化模块，用于最大化所述似然函数，即得到目标函数：其中，θ_j＝{a_j,σ_j,c_j,b_j,α}是参数向量，具体地，a_j表示第j个瑞利分布在横轴上的平移量，σ_j表示瑞利分布的众数，c_j和b_j是计算参数，α表示控制变量；π_j是先验概率，P(θ_j|y_i;Θ^(t))是后验概率，P(y_i|θ_j;Θ)是类条件概率密度，Θ^(t)表示第t次迭代中已知的参数集，Θ表示第t次迭代中未知的参数集，Θ^*表示计算获得的最终参数集； 

参数集初始化模块，用于初始化参数集Θ为Θ⁽⁰⁾； 

参数集计算模块，用于根据初始化的参数集Θ⁽⁰⁾计算参数向量，并使用最速下降算法更新所述参数向量，直到最大期望算法收敛，获得最终参数集Θ^*。 

7.根据权利要求3所述的系统，其特征在于，所述系统还包括： 

最大值线获取模块，用于确定所述第一图像每列的像素中灰度值最大的像素，获得最大值线。 

8.根据权利要求7所述的系统，其特征在于，所述特征提取模块包括： 

第二区域确定模块，用于计算所述第二图像中每个聚类簇对应的感兴趣区域中的像素的灰度值均值，找到灰度值均值最小的预设个数的聚类簇，将所述找到的聚类簇中的像素构成的区域作为第二区域； 

第三图像生成模块，用于生成与所述第二图像尺寸相同的第三图像，将所述第三图像中属于第二区域的像素的灰度值和不属于第二区域的像素的灰度值分别设置为不同的值； 

第二区域上边界确定模块，用于在所述第三图像中确定第二区域的上边界； 

特征参数获取模块，用于计算所述第三图像每列的底部与所述第二区域的上边界的距离、所述第三图像每列对应的第二区域的上边界与所述最大值线的距离以及所述第一图像的每列对应的所述第二区域的上边界和所述最大值线之间的像素的平均灰度值； 

特征计算模块，用于根据所述第三图像每列的底部与所述第二区域的上边界的距离、所述第三图像每列对应的第二区域的上边界与所述最大值线的距离以及所述第一图像的每列对应的所述第二区域的上边界和所述最大值线之间的像素的平均灰度值确定所述第一图像每列的特征。 

9.根据权利要求8所述的系统，其特征在于，所述第一图像每列的特征为： 

其中h₁,h₂,h₃是权重，表示所述第三图像每列的底部与所述第二区域的上边界的距离，表示所述第三图像每列对应的第二区域的上边界与所述最大值线的距离，表 示所述第一图像的每列对应的所述第二区域的上边界和所述最大值线之间的像素的平均灰度值。 

10.根据权利要求9所述的系统，其特征在于，h₁＝5,h₂＝-0.5,h₃＝-1。 

11.根据权利要求3所述的系统，其特征在于，所述分类模块包括： 

置信计算模块，用于根据所述第一图像每列的特征，使用置信传播算法在马尔可夫随机场上计算出第一图像每列含有钙化斑块和声影的图像的置信； 

分类执行模块，用于根据所述第一图像每列的置信将所述第一图像中的各列分为含有钙化斑块和声影的图像的列和不含有钙化斑块和声影的图像的列两类。 

12.根据权利要求11所述的系统，其特征在于，所述置信计算模块包括： 

位置集和状态集定义模块，用于定义马尔可夫随机场的位置集S和状态集L分别为其中，S＝{1,2,...,W}表示该列的位置；L＝{-1,+1}表示该列的状态，若第一图像I的一列的状态为“+1”，表示该列含有钙化斑块及其声影的图像；如果第一图像I的一列的状态为“-1”，则表示该列不含有钙化斑块及其声影的图像； 

初始化模块，用于对于马尔可夫随机场上的表示图像I每列的状态的每个隐藏变量z_i，隐藏变量z_i的边缘概率为：P(z_i＝-1)＝P(z_i＝1)＝0.5；且隐藏变量z_i的置信初始化为局部信息φ_i(z_i,χ_i)＝F_i，F_i是第一图像第i列的特征，χ_i是观测变量χ₁,...,χ_W，表示图像I每列的特征；相容函数为 当节点z_i和z_j是互为邻域，那么从z_i传递到z_j的信息初始化为m_i,j(z_j)＝1; 

迭代模块，用于在第t次迭代中，计算从z_i传递到z_j的信息

并计算节点z_i的置信

其中N(i)是节点z_i的邻域； 

计算模块，用于当在t+1次迭代中，满足迭代算法收敛，则每个隐藏变量的置信为：b^*(z₁)＝b^(t+1)(z₁),b^*(z₂)＝b^(t+1)(z₂),...,b^*(z_W)＝b^(t+1)(z_W)；其中ε为预设定值，W表示第一图像的列数。 

13.根据权利要求12所述的系统，其特征在于，所述分类执行模块还用于对于所述第一图像的第i列的状态z_i，b^*(z_i＝-1)＞b^*(z_i＝+1)，那么将第i列分到类“-1”；如果b^*(z_i＝-1)≤b^*(z_i＝+1)那么将第i列分到类“+1”；其中，若该列属于类“+1”表示该列含有钙化斑块及其声影的图像，若该列不属于类“-1”则表示该列不含有钙化斑块和声影图像。 

14.根据权利要求7所述的系统，其特征在于，所述提取模块包括： 

第三区域确定模块，用于根据所述第一区域和所述第一图像的最大值线确定第三区域； 

过滤模块，用于过滤掉不符合预设约束条件的第三区域； 

提取执行模块，用于根据过滤后的第三区域提取钙化斑块及其声影的图像。 

15.根据权利要求14所述的系统，其特征在于，所述预设约束条件包括：满足以下四个条件中的至少一个条件： 

所述第三区域的上边界对应的像素中，灰度值大于第一阈值的像素数量与所述第三区域的上边界对应的所有像素的数量的比大于第二阈值； 

所述第三区域的上边界的最左端和最右端的像素的列坐标的差与行坐标的差的比的绝对值小于第三阈值； 

所述第三区域的左边界和右边界的之间的距离小于第四阈值； 

所述第三区域中像素的灰度值均值为m_m，所述第三区域左邻的不包含钙化斑块及其声影的图像的区域中像素的灰度均值为m_l，所述第三区域右邻的不包含钙化斑块及其声影的图像的区域中像素的灰度均值为m_r，T₅为第五阈值，则m_l＞m_m,m_r＞m_m

16.根据权利要求14或15所述的系统，其特征在于，所述预设约束条件包括满足以下三个条件中的至少一个条件： 

所述第三区域对应的上升沿和下降沿的所有像素中最大的行坐标与所述第三区域的上边界中的像素的灰度值均值的差大于第六阈值； 

其中，N_ab表示列坐标在和之间的，行坐标在1到 之间的像素的个数，表示所述第三区域对应的上升沿最左边点的列坐标，表示所述第三区域对应的下降沿最右边点的列坐标，表示所述第三区域对应的上升沿和下降沿的所有像素中最大的行坐标，T₇表示第七阈值； 

其中，N_ub表示列坐标在和之间的，行坐标在 到H之间的像素的个数，表示所述第三区域对应的上升沿最左边点的列坐标，表示所述第三区域对应的下降沿最右边点的列坐标， 表示所述第三区域对应的上升沿和下降沿的所有像素中最大的行坐标，T₈表示第八阈值。 

17.根据权利要求14所述的系统，其特征在于，所述提取执行模块包括： 

声影图像提取模块，用于根据过滤后的第三区域提取声影图像所在的区域； 

钙化斑块图像提取模块，用于根据过滤后的第三区域提取钙化斑块所在的 区域； 

边界提取模块，用于提取钙化斑块及其声影的图像的边界。