CN101833757A - 血管图像组织结构的边缘检测及血管内膜检测方法和系统 - Google Patents

Abstract

本发明公开了一种血管图像组织结构的边缘检测及血管内膜检测方法和系统，包括：获取二维超声血管灰度图像；基于血管的结构特点搜索血管的中心位置，设定感兴趣区域；查找感兴趣区域的血管腔与血管内膜，和/或血管中层与血管外膜的初始分割边界，基于所述分割区域的灰度分布函数，建立能量泛函，将所述能量泛函的极值所对应的位置作为相邻区域间的边界位置；计算并显示血管壁参数。本发明将基于灰度分布的能量泛函引入到了组织结构区域边界的检测中，通过分析各个组织结构的灰度分布特点，对组织结构的边缘进行更强的鲁棒性估计，并避免陷入局部极值。

Description

血管图像组织结构的边缘检测及血管内膜检测方法和系统

技术领域

本发明涉及超声图像组织结构的边缘检测方法，具体涉及一种超声血管图像组织结构的边缘检测及其在血管内膜检测中的应用。

背景技术

超声图像中进行组织结构边界分割时通常均是采用灰度极值检测等方法，单一灰度极值检测方法会存在一些的弊端。比如血管内膜中层厚度(Intima-Media Thickness，IMT)检测方法，此IMT检测是基于图像检测血管内膜-中层厚度的一种方法，IMT能够准确反映出动脉管的硬化情况，而颈动脉的IMT检测能够精确的识别、反映出早期的动脉硬化疾病。IMT检测通常包括：感兴趣区域设定、两个分界面的检测和结果计算显示等几个环节构成。如图1所示，给出了一种典型的IMT即内膜中层厚度检测系统的框图，其包括：超声探头101向人体相应检查部位(颈动脉)发射超声波102，成像扫描可以通过ECG控制器110触发控制；接收到的回波信号103经过波束合成器104等信号处理环节，送入图像处理模块105获得人体组织(颈动脉)的解剖结构灰度图像，并保存到电影回放数据存储单元106。图像计算模块108读取图像数据，进行边缘检测和计算，图像计算模块108根据得到的边缘信息，结合已有统计数据模型，进行各项IMT功能参数计算，得到的边缘信息和功能参数送到显示器109进行显示。

针对血管的结构特点，现有技术的IMT检测中利用上述系统在进行边缘检测的时候通常会引入先验的结构知识，如假设血管壁外膜区域为血管壁最亮区域。在边缘的搜索过程用一般都沿垂直于血管壁的方向在每条线上搜索灰度的极大值、灰度梯度极大值以及灰度梯度的二阶导数等局部特性，获得边缘点的位置；然后通过边缘点的拟合、平滑获取边缘的位置。如图2所示，给出了颈动脉血管由内到外分别为：血管腔Lumen 204、血管内膜Intima 203、血管中层Media 202和血管外层Adventitia 201。而我们所要获取的Intima-Media Thickness是指Lumen-Intima界面和Media-Adventitia界面的间的距离；图3为沿图2中205线的灰度变化曲线。上述IMT自动检测的通用方法是基于图像灰度梯度，即从血管中心开始搜索梯度的两个极大值点位置作为两个边界位置(图3所示虚线301，302中间区域)，或者搜索两个灰度极大值的位置来确定分界面。这些方法容易受噪声影响，并且容易陷入局部极值点。也有采用类似区域增长的方法，如基于灰度的相似性获取外膜、内膜等的位置，然后基于血管腔与内膜间边界以及外膜与血管壁中层间边界，通过计算这两个边界间的距离来计算IMT。但是这种简单基于亮度阈值的方法，在边缘亮度不均匀时不容易获得连续边界。

可见，现有技术中进行超声图像组织结构边缘检测，特别是在IMT检测时存在一定的问题，需要进步的改进技术以获得更好的效果。

发明内容

本发明的目的在于提供了一种超声图像组织结构的边缘检测及血管内膜检测方法和系统，其在提高边缘检测的精度的基础上，保证了IMT计算的准确性。

为了解决上述目的，本发明采用如下技术方案：

本发明提供了一种提取血管图像中组织结构的边缘检测方法，所述方法包括：

确定灰度图像中的感兴趣区域；

采用基于灰度分布的能量泛函进行区域分割，将搜索得到的能量泛函极值对应的边缘位置作为相邻分割区域间的边界。

基于上述边缘检测方法，本发明还提供了一种血管内膜自动检测方法，其包括：

获取二维超声血管灰度图像；

基于颈部动脉血管的结构特点搜索血管的中心位置，设定感兴趣区域；

基于灰度分布的能量泛函进行区域分割，将搜索得到的能量泛函极值对应的边缘位置作为相邻分割区域间的边界；

计算并显示血管壁参数。

从上述技术方案可见，本发明将基于灰度分布的能量泛函引入到了组织结构区域边界的检测中，通过分析各个组织结构的灰度分布特点，对组织结构的边缘进行更强的鲁棒性估计，并避免陷入局部极值，减少噪声影响，其特别可以用于到血管内膜自动检测方法，在提高边缘检测的精度的基础上，保证了IMT计算的准确性。

基于上述方法，本发明还提供了一种超声图像系统，所述图像计算模块包括：区域初始化模块，用于确定灰度图像中的感兴趣区域；组织区域边界自动检测模块，用于基于灰度分布的能量泛函进行区域分割，将搜索得到的能量泛函极值对应的边缘位置作为相邻分割区域间的边界；及组织区域分析模块，用于根据所述边界位置计算组织区域的多项参数。本发明的系统基于上述方法的边缘检测结果，将各个组织结果的参数用于显示，可以通过算法编程或者增加硬件设置的方式将实现上述方法的模块加入到现有的超声系统中。

本发明还提供了一种提取并显示图像中组织结构边缘特征的方法，包括以下步骤：

A、在图像中确定目标组织结构的边界位置；

B、统计所述边界位置中各个边缘点灰度邻域的灰度分布相似性，获取所有边缘点间的相似性的平均值，将所述平均值作为所述边界位置估计的置信度，用以确定所述边界位置的精确性。

上述方法基于目标组织结构不同位置边缘点邻域灰度的相似性定义了一个度量，用于标识所述提取组织结构边缘的一致性，相对于现有技术，更加精确的体现了边缘检测的结果。

附图说明

图1为超声血管内膜自动检测系统流程图；

图2为血管结构示意图；

图3为血管及血管壁灰度分布示意图；

图4为感兴趣区域选取示意图；

图5为血管壁内膜中层自动检测流程图；

图6为内膜中层厚度结果显示示意图。

具体实施方式

本发明提供的提取血管图像中组织结构的边缘检测方法主要是基于超声二维灰度图像，首先查找灰度图像中的感兴趣区域；然后，基于灰度分布的能量泛函进行区域分割，将搜索得到的能量泛函极值对应的边缘位置作为相邻分割区域间的边界。这里所说的灰度图像可以是人体任意组织结构的血管灰度图像，特别是针对颈总动脉、颈内动脉和股动脉等二维超声血管灰度图像，且图像中至少包含血管壁的内膜和外膜。以下将详细说明在具体实现时各个步骤的实现方式。

第一点，能量泛函的建立，此目的就是将区域分布函数映射到一个实数上用于判断边界检测的效果，避免进行边缘检测时陷入局部极小值。

本发明假设图像可以分成前景和背景两部分，各部分灰度分布一致，不同部分灰度分布有较大差异。于是，将待计算的感兴趣区域看作是区域中各个像素点坐标和灰度对应关系的集合，即在感兴趣区域中任意一像素点灰度符合一分布函数的集合，每一集合按照概率统计的思想可以看作是一个样本空间，那么用单个或者多个样本空间灰度的数字特征则可以表示能量泛函，这样就可以考虑各个区域的整个灰度分布情况，避免陷入局部极值。这里的单个或者多个区域灰度的数字特征可以是单个区域的灰度方差、多个区域的灰度方差和、多个区域的灰度平方和、单个区域的灰度分布k阶矩、或相邻区域间灰度分布k阶矩之差的平方和等等，以下以具体公式来体现上述几种能量泛函。

1、多个区域的灰度方差和，如下式：

$F_1\left(c\right)=\underset{i=1,...N}{\mathrm{\Σ}}{D}_{{\mathrm{\Ω}}_i}=\underset{i=1,...,N}{\mathrm{\Σ}}{\∫}_{{\mathrm{\Ω}}_i}{[u\left(l\right)-\stackrel{\‾}{u_{{\mathrm{\Ω}}_i}}]}^2\mathrm{dl}---\left(1\right)$

上述公式中，i表示区域标记数，N表示同时参与能量泛函计算的区域总数，Ω_i表示感兴趣区域中分割的第i个区域，l表示区域Ω_i中每一个像素点的相对坐标，其也可以是绝对坐标，但要相应改变上述公式的表达形式。u(l)表示区域Ω_i中点(x_l，y_l)的灰度，

表示区域Ω_i的灰度均值；

表示方差；F₁(c)表示能量泛函的实数值。

2、针对感兴趣区域中的任意一个区域的灰度方差，如下式所示：

$F_1\left(c\right)={\∫}_{{\mathrm{\Ω}}_i}{[u\left(l\right)-\stackrel{\‾}{u_{{\mathrm{\Ω}}_i}}]}^2\mathrm{dl}---\left(2\right)$

3、针对感兴趣区域中的任意一个区域中各点的灰度均值，如下式所示：

$F_1\left(c\right)\stackrel{\‾}{U_{{\mathrm{\Ω}}_i}}---\left(3\right)$

4、针对感兴趣区域中的任意一个分割区域的灰度分布k阶矩，如下式所示：

$F_1\left(c\right)={\∫}_{{\mathrm{\Ω}}_i}{l}^{k}u\left(l\right)\mathrm{dl}---\left(4\right)$

上式中，k表示阶数。

5、相邻区域间灰度分布k阶矩之差的平方和，如下式所示：

$F_1\left(c\right)=\underset{i=1,...,N}{\mathrm{\Σ}}{({\∫}_{{\mathrm{\Ω}}_i}{l}^{k}u\left(l\right)\mathrm{dl}-{\∫}_{{\mathrm{\Ω}}_{i+1}}{l}^{k}u\left(l\right)\mathrm{dl})}^2---\left(5\right)$

6、相邻区域间灰度均值之差的和，如下式所示，即公式(5)中的k＝0，

$F_1\left(c\right)=\underset{i=1,...,N}{\mathrm{\Σ}}{(\stackrel{\‾}{u_{{\mathrm{\Ω}}_i}}-\stackrel{\‾}{u_{{\mathrm{\Ω}}_{i+1}}})}^2---\left(6\right)$

7、针对感兴趣区域中的单个或者多个区域的灰度分布概率之和，如下式所示：

$F_1\left(c\right)=\underset{i=1,...,N}{\mathrm{\Σ}}\left({\∫}_{{\mathrm{\Ω}}_i}\log P(u\left(l\right);m_i,{\mathrm{\α}}_i)\mathrm{dl}\right)---\left(7\right)$

其中，F₁(c)表示能量泛函的实数值，P(u(l)；m_i，α_i)表示已知的参数为m_i，α_i的灰度分布函数，例如：高斯分布，瑞利(Rayleigh)分布或者偏移瑞利(Rayleigh)分布等；N表示同时计算的区域数量。Ω_i表示感兴趣区域中分割的第i个区域，l表示区域Ω_i中每一个像素点(x_l，y_l)的坐标，u(l)表示区域Ω_i中像素点(x_l，y_l)的灰度，i表示区域标记数，N表示同时参与能量泛函计算的分割区域总数。

针对上述前6个计算公式，可以直接代入图像数据进行计算，然后选择极小值或者极大值作为新的边界位置；而针对上述公式(7)，则区域样本空间为所述感兴趣区域中任一灰度符合一分布函数的像素点的集合，可以符合高斯分布分布，瑞利分布、偏移瑞利分布等概率分布中的一种，则将其相应的分布函数代入上述公式(7)进行计算即可。能量泛函还可以采用如下表现形式：

从上述公式可以看出，能量泛函还可以用一个或多个样本空间的灰度分布概率之和，此样本空间为所述感兴趣区域中任一灰度符合一分布函数的像素点的集合。

此外，在上述能量泛函F₁(c)的基础上，还可以增加利用沿边界的导数进行加权计算所获得的边缘约束条件，如下式：

$F_3\left(c\right)=\mathrm{\β}{\∫}_0^1{\left|c^{\′}\left(q\right)\right|}^2\mathrm{dq}+\mathrm{\γ}{\∫}_0^1\left|c^{\′\′}\left(q\right)\right|\mathrm{dq}---\left(8\right)$

其中，F₃(c)表示加权计算的结果，β，γ表示权重，c(q)：[0，1]→R²表示边界。这样做可以提高边缘的连续性和光滑性。

为了更加准确的确定边界位置，并保证边缘的连续性和图像梯度大小，则对于采用上述F₃(c)+F₁(c)的能量泛函获得的边界位置，还要添加以测地距及灰度梯度的能量泛函取极值做为边界条件，如下式：

$F_2\left(c\right)=\mathrm{\λ}{\∫}_0^{1}g\left(\right|\▿u\left(c\left(q\right)\right)\left|\right)\left|c^{\′}\left(q\right)\right|\mathrm{dq}---\left(9\right)$

其中，F₂(c)表示能量泛函的实数值，λ表示权重，g(x)＝1/(1+x²)，c(q)：[0，1]→R²表示边界。通过公式(9)筛选查找到的相邻区域间的边界位置，保留此能量泛函取极值的像素点，可以提高边缘检测的精度，保证估计结果的鲁棒性。

第二点，确定感兴趣区域及能量泛函的计算起始点。

感兴趣区域可以手动在图像中进行选择，或者根据各个组织结构的特点，利用灰度分布进行自动选择，比如下述针对血管内膜检测所采用的方法。

如图4所示，获得颈动脉血管二维灰度图像，选择颈动脉位置，采用超声高频线阵探头平行于颈动脉走向进行扫描，调整探头角度获得清晰的内膜图像；冻结图像，选取内膜较清楚图像进行计算。感兴趣区域的自动设定过程包括：

第一步，在血管灰度图像中设定一个目标搜索区域，比如手动图4中白色实线框为取定的目标区域，便于快速准确的搜索。

第二步，根据组织结构的灰度特性确定血管内腔区域为起点位置。因为动脉血管腔是一个连续的黑色的区域，其位于灰度极小值的位置，比较容易查找，所以以此作为起点位置。设定血管直径一半大小的窗口，沿扫描线方向进行滤波检测，获取灰度最小点为血管内腔区域，设为血管中心，最小灰度为Min_min，并获取灰度最大值Max_max。

第三步，从起点位置开始沿扫描线分别检测该起点位置两侧血管壁位置，获取目标搜索区域内血管腔区域的长度。在这里，检测方法是基于近似1/5血管直径窗口的滤波和阈值灰度校正获取沿线上扫描区域内血管腔区域长度，近似为血管直径，这里的阈值为

Min_min+(Max_max-Min_min)×α                     (10)

上述α为一常数，根据经验设定。

第四步，以所述血管腔区域的长度(即近似的血管直径)作为参考判断血管的位置走势，若所述长度连续增大，判断为颈动脉球体区域(Carotidbulb)，即颈内动脉(Internal carotid)和颈外动脉(External carotid)分叉区域，将所述血管分叉前一预设长度处设置为感兴趣区域，比如设定感兴趣区域在分叉前1-1.5cm处；否则，将所述目标搜索区域的中间位置设置为感兴趣区域。

上述方法能够初步检测出血管区域大体位置，并设定感兴趣区域的初步位置。该方法把初始感兴趣区域设定在颈动脉血管附近，尽量降低低效的人工操作的同时提高自动检测的重复性，用户还可以对自动搜索的结果进行调整。而感兴趣区域的大小根据实际需要做调整，最优选择是：沿血管壁设定1cm到2cm的长度，垂直血管壁1cm宽度(要求大于1.5毫米)的区域。

在感兴趣区域确定的基础上，需要选择一个能量泛函的计算起点。根据上述能量泛函建立的原理，在感兴趣区域中涉及两个以上区域计算时，需要有一定的边界将感兴趣区域预先分割成相应的区域，然后再套用上述公式进行计算。而这样的边界可以是感兴趣区域中任意一条连续边界，也可以是通过检测获得的初始分割边界，比如，可以利用基于区域灰度分布一致性获取边缘的方法，也可以利用基于灰度梯度和亮度值等的边缘检测方法来进行初始分割边界的检测。

假设，采用感兴趣区域中的任意一条连续的边界进行初始分割来套用能量泛函公式的话，可以采用以下过程进行：

以感兴趣区域中的任意一个像素点所在的任意连续边界作为初始分割边界，建立分割区域间能量泛函，并计算当前像素点对应的能量泛函；

搜索能量泛函的极值，并将该位置作为相邻分割区域的边界位置；

在所述感兴趣区域内，连接所有扫描线上的所述边界位置，获得所述相邻分割区域间的边界。

采用上述方法则需要将每一扫描线上的像素点都进行计算，因此，为了减少计算量，可以限定一个能量泛函的计算范围，比如采用初始分割边界进行感兴区域的初始划分，然后利用能量泛函纠正初始分割边界的精确度。此方法可以采用如下过程：

在感兴趣区域中搜索初始边缘，获得初始分割边界；

由初始分割边界将所述感兴趣区域分割成至少两个分割区域，计算获得分割区域间的能量泛函，取能量泛函极值点作为所述相邻分割区域间的边界。

上述初始分割边界的查找过程可以采用如下方法：

首先，从感兴趣区域中的一个预定位置开始进行边缘搜索，计算整个区域上的第一灰度均值，比如在血管各腔边界检测时，此预定位置选择血管腔的一侧。

其次，计算灰度小于所述第一灰度均值的像素点的第二灰度均值Gm_ean；

然后，从所述预定位置开始，计算能量

的最小值位置n，其中g_i表示第i个像素点的灰度；

最后，将整个感兴趣区域中所有位置n以内判断为所述预定位置所在的分割区域，则获得一条由多个位置n构成的初始分割边界。这一方法可以应用到各种组织结构的边界检测中，特别适合应用于血管腔的边缘检测。在血管腔的边缘检测时感兴趣区域内主要包含血液区域和血管壁组织区域，这两个区域具有不同的灰度分布，因此从血管腔一侧开始边缘检测，计算整个区域上灰度均值，然后计算灰度小于该均值的点的灰度均值，从血管腔区域开始，计算能量的最小值位置n，其中g_i表示第i个点的灰度，把位置n以内判断为血管腔区域，这样获得一个初始边缘。在初始边缘的基础上进行能量泛函计算可以减小计算量，更快接近极值点。

第三点，为了提高检测精度，避免不相关数据对分割结果的影响，在进行边缘检测之前，需要对感兴趣区域内每条扫描线上的点灰度进行校正处理。感兴趣区域的灰度校正环节用于规范化感兴趣区域内的灰度分布剔除无效数据的影响，其可以采用更简单的方式进行灰度校正环节，比如沿一定方向搜索灰度极大值点为校正点，把校正点以外的点的灰度赋为校正点的灰度值乘以一个系数以进行边界检测；或者只对所述校正点以内的区域进行边界检测，比如，沿扫描线方向搜索灰度极大值或极小值，并将所述极大值点或极小值作为校正点；计算区域的灰度平均值；利用所述灰度极大值或极小值和灰度平均值的加权计算结果，对扫描线所述校正点以内的区域进行灰度校正。下面以血管腔边缘检测为例详细说明。

通常颈动脉外膜为整个感兴趣区域内最亮的位置，如图3。沿每条扫描线搜索灰度的极大值Grey_Max，即外膜位置；计算血管腔初始区域内的平均灰度Grey_LumenMean。为了避免极大值点后面数据对分割结果的影响，采用(1-α)*Grey_LumenMean+α*Grey_Max对扫描线上灰度极大值点后面的区域进行灰度校正。参数α在检测不同边界时选择不同的数值：在进行血管腔与内膜边缘检测时采用较小α，减少外膜区域内灰度与内膜区域灰度间的差异；在进行中层与外膜边界检测时采用较大α，增大极大值点后区域与内膜间差异，保证灰度极大值点后区域灰度与外膜的一致性。为了增加血管腔与内膜间界面分割的稳定性，本发明在搜索到Grey_Max后，沿扫描线向血管腔位置搜索一个局部极大值Grey_IntimaMax(内膜位置)来代替当前校正点Grey_Max。在进行中层与外膜边界检测时仍然使用Grey_Max作为校正位置。通过灰度校正处理，保证当前搜索区域内只包含血管腔内血液区域和血管壁区域两个连续区域，而且血管壁区域只包含内膜、中层和外膜三个部分。

第四点，为提高边缘精度，本发明还可以包含亚象素的搜索方式，即在获得边界点位置后，提取边界点附近，沿扫描线前、后共5个点的灰度。对这五个点的灰度进行三次多项式拟合或者样条拟合，获得拟合曲线。在上述边界点前、后各两个像素范围内，搜索梯度最大值位置作为新的边界点的位置。在上述边缘检测结果由于图像的不均匀，会产生边缘的毛刺，需要采用平滑方法去掉毛刺。通常采用的边界平滑方法为均值平滑或者曲线拟合(分别拟合两个界面以修正局部不准确位置)。常用的拟合方法为线性拟合(因为两个界面近似线性)、二次多项式拟合，样条拟合等等。先采用分段多项式拟合方法分别拟合两个界面，剔除离拟合出来曲线较远的点，然后采用均值滤波的形式进一步平滑边缘轮廓。在检测出边缘后允许用户手动对边界位置进行调整。

在上述边缘检测方法的基础上，本发明还增加了对采用上述方法确定的边界给予置信度的测量。如下所示：

在图像中确定目标组织结构的边界位置后，统计所述边界位置中各个边缘点灰度邻域的灰度分布相似性，获取所有线数据间的相似性的平均值；将所述平均值作为所述边界位置估计的置信度，用以确定所述边界位置的精确性。该参数能够充分反应检测结果的准确性和可靠性，而且能够同时反映图像质量的好坏。图像噪声大的时候该参数较低，图像噪声小而且边缘上灰度分布一致时，边缘检测置信度较高。若所述置信度小于一阈值，则手动进行边缘拟合，或者还可以将相应的置信度计算结果直接输出以作参考。此阈值可以预先根据实际情况设置，或者直接输出估计出的置信度。

上述所有线数据间的相似性的平均值的计算采用如下公式：

${\mathrm{\κ}}_2=\frac{\underset{i=1,...N}{\mathrm{\Σ}}(x_i-\stackrel{\‾}{x})(y_i-\stackrel{\‾}{y})}{\sqrt{\left(\underset{i=1...N}{\mathrm{\Σ}}\left({x_i-\stackrel{\‾}{x})}^2\right)(\underset{i=1...N}{\mathrm{\Σ}}{(y_i-\stackrel{\‾}{y})}^2\right)}}---\left(11\right)$

其中，κ₂表示所有线数据间的相似性的平均值；i＝1，...，N表示灰度邻域R_neib内的所有像素点；x_i，y_i分别表示不同边缘点邻域

内的像素点的灰度；

表示

的平均灰度。这里的灰度邻域R_neib，可以定义为每条边界上的边缘点获取每个边缘点前后各3个点的灰度做为该点的灰度邻域R_neib。

还可以采用另一公式表示：

${\mathrm{\κ}}_3=\frac{\underset{i=1...N}{\mathrm{\Σ}}\left(x_i{y}_i\right)}{\sqrt{\left(\underset{i=1...N}{\mathrm{\Σ}}{x}_i^2\right)\left(\underset{i=1...N}{\mathrm{\Σ}}{y}_i^2\right)}}---\left(12\right)$

其中，κ₂表示所有线数据间的相似性的平均值；i＝1，...，N表示灰度邻域R_neib内的所有像素点；x_i，y_i分别表示不同边缘点邻域

内的像素点的灰度。这里边缘点邻域灰度的相似性可以是灰度的相关性度量等各种计算两个区域内灰度分布相似性的度量，如公式(11)和(12)；本发明通过定义一个度量对目标结构边缘不同部分区域分别计算个区域内边缘点的灰度分布相似性，然后把这几个区域产生的相似性系数加权生成整个边缘的一致性度量。

由于，在边缘拟合的计算过程中，采用分段多项式对边界点进行拟合。则还可以基于给定的距离阈值ω，对距离拟合出来曲线距离大于阈值ω的点，采用拟合出来的界代替原先边界位置，则统计在拟合过程中位置未改变的点的数量，除以边缘总点数即为获得边缘轮廓的置信度κ₁，也就是位置未发生改变的像素点占整个所述边界位置总像素点的比例值，该度量用于显示准确检测到边界的点的扫描线的数量。然后，将此比例值和相似性的平均值取平均作为所述边界位置估计的置信度。

本发明在置信度κ₁的基础上，还针对边界上的边缘点计算不同边缘点的灰度邻域R_neib的灰度分布的相关性，如κ₂和κ₃，然后，结合上述两个参数κ＝(κ₁+κ₂)/2或者κ＝(κ₁+κ₃)/2作为总的置信度。通常现有技术中只采用置信度κ₁的方式，这样只能显示检测边缘的光滑性，不能从图像的角度深入反映检测结果的正确性，而本发明的方法从边缘的连续性和边缘点的自相似性给出了置信度的综合度量，能较好的反应边缘检测的结果。上述置信度的综合度量方法不仅可以度量上述边缘检测的结果，而且还可以用于度量其他目标组织结构的检测结果。

基于上述边缘检测方法，本发明将其应用于血管内膜的自动检测上，给出了一种血管内膜自动检测方法，其针对血管壁中各部分组织灰度分布的差异，基于各区域灰度分布进行各区域边缘检测，能有效克服噪声的影响。其具体过程如图5所示：

400，输入颈动脉超声图像，其可以是颈总动脉、颈内动脉和股动脉等二维超声血管灰度图像；

402，基于颈部动脉血管的结构特点搜索血管的中心位置，设定感兴趣区域；

404，在感兴趣区域内，基于灰度分布的能量泛函进行区域分割，将所述能量泛函的极值所对应的位置作为相邻分割区域间的边界位置，比如基于区域灰度分布搜索血管腔和血管内膜的边界；

406，利用上述步骤404的原理，在感兴趣区域内搜索血管中层和血管外膜的界面；

408，计算并显示IMT。这里IMT计算时，边界面间厚度的定义方式可以为：计算某个层面中一个点到对面界面点的最短距离，作为当前点位置的厚度；或者计算某个层面中一个点位置的法线与对面界面的交点，以这两个点作为两个对应点，两点间距离作为该处厚度；或者计算扫描线与拟合出来的边界所成角度，以扫描线上内膜-中层厚度乘以所得角的余弦的绝对值作为IMT。上述步骤的具体操作请见上述相关说明。上述血管内膜的自动检测方法中，步骤402后还可以通过查找感兴趣区域的血管腔与血管内膜，和/或血管中层与血管外膜的初始分割边界来确定参与能量泛函的计算范围，用以减少计算量。此初始分割边界将感兴趣区域分割成血管腔、血管内膜分割区域和/或血管中层、血管外膜分割区域，然后再进行能量泛函更快接近极值点。

通常采用上述方法确定的血管感兴趣区域主要包含血液区域和血管壁组织区域，这两个区域具有不同的灰度分布，另外，血管壁区域只包含内膜、中层和外膜三个部分，即三个分割区域。采用上述公式(1)进行血管腔和内膜间的边界检测时，进行双区域能量泛函，搜索得到能量泛函最小的边缘位置为血管腔和内膜间的边界。假设初始边界c把感兴趣区域L分成inside(L)和outside(L)，设c₁为内部区域inside(L)的灰度均值，设c₂为外部区域outside(L)的灰度均值，u(l)为区域中点(线上)(x_l，y_l)的灰度。则确定边界c后，能量泛函定义为：

F₁(c)＝∫_nside(L)u(l)-c₁)²dl+∫_outside(L)u(l)-c₂)²dl                  (13)

上述公式(13)是基于两个区域，首先计算获得一个边界，然后基于这个边界在剩下区域内基于灰度分布确定另一侧边界。同理，还可以基于多个区域灰度分布能量组合的范函同时检测感兴趣区域内的多个边界，如下式所示：

F₁(c)＝∫_Region(1)(u(l)-c₁)²dl+∫_Region(2)(u(l)-c₂)²dl+∫_Region(3)(u(l)-c₃)²dl        (14)

其中：Region(i)，i＝1，2，3表示感兴趣区域内的三个区域，c_i，i＝1，2，3分别表示Region(i)，i＝1，2，3内的灰度均值。从上述公式(13)和(14)可见边界在两区域边缘位置时，上述能量函数F₁(c)最小，所以采用这种方法对初始边界依赖小，受噪声影响小。

在进行中层与外膜间边界检测时，在已知血管腔与内膜的边界情况下，以去掉血管腔后剩下区域作为感兴趣区域。由于在剩下的感兴趣区域内，经过灰度校正后只包含内膜、中层和外膜三个部分，而且血管外膜的亮度(灰度分布)显著高于其它两层组织，因此可继续采用基于灰度分布的方法，进行中层和外膜间的边缘检测。即采用上述搜索方法，获取血管壁中层与外膜间边界如图3中302位置。当获取的颈动脉图像较深，或者图像分辨较低时，直接基于图像象素的IMT边缘检测结果往往达不到精度要求。此时，在进行边缘检测之前，对获取的感兴趣区域内每条扫描线进行加密采样，比如采用分段多项式拟合方法对每条扫描线上点的灰度分布进行拟合，并插值出更小象素点的灰度。基于插值出来的扫描线上点的灰度分布，进行边缘检测。

本发明提出的血管内膜自动检测方法，与现有技术不同，其基于血管腔、血管内膜与血管外膜灰度分布的不同，采用计算基于灰度分布的能量泛函来搜索区域边界，具有更强的鲁棒性。在边缘检测过程中，首先基于血管位置特征，设计感兴趣区域的自动检测方案获取感兴趣区域；在感兴趣区域内，利用血管腔灰度分布与整个血管壁区域灰度分布的差异，基于灰度分布的能量泛函获取血管腔与血管壁内膜间的界面；然后在剩下区域内，利用血管内膜和中层与血管壁外膜的灰度分布差异，获取血管壁中层与血管壁外膜间的分界面。上述边缘检测方法在采用基于全局的搜索，因此有较强的稳定性，避免陷入局部极小值。基于分割出来的边界计算IMT(内膜-中层)厚度，并对计算结果做统计分析，可以给出给出了新的有效显示IMT的方式，提高精确度。

基于上述方法，本发明还提供了一种超声图像中组织结构特征的提取及显示系统。如图1所示，本发明可以通过软件编程或者添加相关硬件结果的方式在现有的系统中增加实现上述方法的各个功能模块。本发明的图像计算模块108包括：区域初始化模块，用于确定灰度图像中的感兴趣区域；组织区域边界自动检测模块，用于基于灰度分布的能量泛函进行区域分割，将搜索得到的能量泛函极值对应的边缘位置作为相邻分割区域间的边界；及组织区域分析模块，用于根据所述边界位置计算组织区域的多项参数。为了存储和读取超声图像，在上述系统中还设置有图像数据存储和读取模块。上述各个功能模块的原理可以参见上述方法的相关说明。

上述组织区域分析模块包括：内膜中层厚度计算单元，用于根据不同位置的所述边界位置计算各个位置的内膜中层厚度，并建立内膜中层厚度随位置变化的曲线予以输出显示；置信度计算单元，用于利用边缘点邻域灰度的相似性计算边缘检测结果的置信度，予以输出显示。根据置信度的大小可以随意通过操作控制模块107调整相应的边界位置进行重新计算。

上述组织区域分析模块还包括：方差统计单元，用于根据内膜中层厚度计算单元的输出结果，统计内膜中层厚度的分布方差，并叠加在所述内膜中层厚度的变化曲线中予以输出显示；平均值统计单元，用于根据内膜中层厚度计算单元的输出结果，计算内膜中层厚度的平均值，并叠加在所述内膜中层厚度的变化曲线中予以输出显示。如图6所示，显示厚度的变化曲线的测量结果502，正常人厚度统计平均厚度结果504以及统计方差范围500等信息，也可以在该坐标系内同时显示过去诊断结果。这样能直观给出定量的结果显示。同时在图像中用彩色对厚度进行编码，不同颜色在图像中显示不同的IMT厚度。在已知正常人的平均厚度和分布的方差时，可基于测量结果采用统计分析(T检验)给出患者的IMT是否正常的概率。

本发明根据颈动脉超声图像的特点，提出了一种新的IMT自动检测技术，首先自动计算选取感兴趣区域，然后采用基于灰度分布方法检测获取两个界面位置。本发明采用的边缘检测泛函能够有效的避免超声图像噪声的影响，并且有效检测到弱边缘。把基于灰度分布的能量泛函和其它基于边缘的能量泛函结合，能够获取精确的检测结果，同时本发明系统给出了结果的直观显示方法。

上述各具体步骤的举例说明较为具体，并不能因此而认为是对本发明的专利保护范围的限制，本发明的专利保护范围应以所附权利要求为准。

Claims (23)

1.一种提取血管图像中组织结构的边缘检测方法，其特征在于，所述方法包括：

确定灰度图像中的感兴趣区域；

基于灰度分布的能量泛函进行区域分割，将搜索得到的能量泛函极值对应的边缘位置作为相邻分割区域间的边界。

2.根据权利要求1所述的方法，其特征在于，所述区域分割过程包括：

以一像素点所在的任意连续边界作为初始分割边界，建立分割区域间能量泛函，并计算当前像素点对应的能量泛函；

在每一扫描线上，搜索能量泛函的极值，并将该位置作为相邻分割区域的边界位置；

在所述感兴趣区域内，连接所有扫描线上的所述边界位置，获得所述相邻分割区域间的边界。

3.根据权利要求1所述的方法，其特征在于，所述区域分割过程包括：

搜索初始分割边界；

由初始分割边界将所述感兴趣区域分割成至少两个分割区域，计算获得分割区域间的能量泛函，取能量泛函极值点作为所述相邻分割区域间的边界。

4.根据权利要求3所述的方法，其特征在于，所述初始分割边界的边缘检测过程包括：

从感兴趣区域中的一个预定位置开始进行边缘搜索，计算整个区域上的第一灰度均值；

计算灰度小于所述第一灰度均值的像素点的第二灰度均值G_mean；

从所述预定位置开始，计算能量

的最小值位置n，其中g_i表示第i个像素点的灰度；

将整个感兴趣区域中所有位置n以内判断为所述预定位置所在的分割区域，则获得一条由多个位置n构成的初始分割边界。

5.根据权利要求1所述的方法，其特征在于，所述能量泛函的实数值为：单个或多个样本空间的数字特征，此样本空间为所述感兴趣区域中灰度符合一分布函数的像素点的任一集合。

6.根据权利要求5所述的方法，其特征在于，所述单个或多个样本空间灰度的数字特征包括：单个区域的灰度方差、单个区域的灰度均值，多个区域的灰度方差和、相邻区域间灰度均值之差的和、单个区域的灰度分布k阶矩、或相邻区域间灰度分布k阶矩之差的平方和。

7.根据权利要求1所述的方法，其特征在于，所述能量泛函的实数值为：单个或多个样本空间的点的灰度分布概率之和，此样本空间为所述感兴趣区域中灰度符合一分布函数的像素点的任一集合。

8.根据权利要求7所述的方法，其特征在于，所述能量泛函的实数值表示为：

$F_1\left(c\right)=\underset{i=1,...,N}{\mathrm{\Σ}}\left({\∫}_{{\mathrm{\Ω}}_i}\log P(u\left(l\right);m_i,{\mathrm{\α}}_i)\mathrm{dl}\right)$

其中，F₁(c)表示能量泛函的实数值，P(u(l)；m_i，α_i)表示已知的参数为m_i，α_i的灰度分布函数，Ω_i表示所述感兴趣区域中分割的第i个分割区域，l表示区域Ω_i中每一个像素点(x_l，y_l)的坐标，u(l)表示区域Ω_i中像素点(x_l，y_l)的灰度，i表示区域标记数，N表示同时参与能量泛函计算的区域总数。

9.根据权利要求8所述的方法，其特征在于，所述灰度分布函数为：高斯分布分布，瑞利分布、偏移瑞利分布概率分布中的一种。

10.根据权利要求6或7所述的方法，其特征在于，所述能量泛函中还包括：利用沿边界的导数进行加权计算所获得的边缘约束条件。

11.根据权利要求10所述的方法，其特征在于，所述方法还包括：以测地距及灰度梯度的能量泛函取极值为边界条件，筛选查找到的相邻分割区域间的边界位置，保留此能量泛函取极值的像素点。

12.根据权利要求1所述的方法，其特征在于，所述方法还包括：对所述感兴趣区域进行灰度校正的步骤，其包括：

沿扫描线方向搜索灰度极大值或极小值，并将所述极大值点或极小值作为校正点；

计算区域的灰度平均值；

利用所述灰度极大值或极小值和灰度平均值的加权计算结果，对扫描线所述校正点以内的区域进行灰度校正。

13.利用权利要求1所述边缘检测的一种血管内膜自动检测方法，其特征在于，所述方法包括：

获取二维超声血管灰度图像；

基于颈部动脉血管的结构特点搜索血管的中心位置，设定感兴趣区域；

基于灰度分布的能量泛函进行区域分割，将搜索得到的能量泛函极值对应的边缘位置作为相邻分割区域间的边界；

计算并显示血管壁参数。

14.根据权利要求13所述的方法，其特征在于，所述感兴趣区域的自动设定过程包括：

在所述血管灰度图像中设定一个目标搜索区域；

根据组织结构的灰度特性确定血管内腔区域为起点位置；

从所述起点位置开始沿扫描线分别检测该起点位置两侧血管壁位置，获取所述目标搜索区域内血管腔区域的长度；

以所述血管腔区域的长度作为参考判断血管的位置走势，若所述长度连续增大，则判定为血管分叉，将所述血管分叉前一预设长度处设置为感兴趣区域；否则，将所述目标搜索区域的中间位置设置为感兴趣区域。

15.根据权利要求14所述的方法，其特征在于，所述起点位置两侧血管壁的位置的检测过程基于阈值灰度校正，此阈值为下式所示：

Min_min+(Max_max-Min_min)×α

其中，Min_min表示所述血管腔区域的灰度最小值，Max_max表示所述目标搜索区域的灰度最大值，α为一常数。

16.一种提取并显示图像中组织结构边缘特征的方法，其特征包括：

A、在图像中确定目标组织结构的边界；

B、统计所述边界位置中各个边缘点灰度邻域的灰度分布相似性，获取所有这些边缘点间相似性的平均值，将所述平均值作为所述边界位置估计的置信度，用以确定所述边界位置的精确性。

17.根据权利要求16所述的方法，其特征在于，所述步骤B中还包括：

对所述边界位置进行边缘拟合，统计所述拟合过程中位置发生改变的像素点占整个所述边界位置总像素点的比例值；

所述步骤B中将所述比例值和相似性的平均值之和作为所述边界位置估计的置信度。

18.根据权利要求16所述的方法，其特征在于，所述所有线数据间的相似性的平均值的计算采用如下公式：

${\mathrm{\κ}}_2=\frac{\underset{i=1...N}{\mathrm{\Σ}}(x_i-\stackrel{\‾}{x})(y_i-\stackrel{\‾}{y})}{\sqrt{\left(\underset{i=1...N}{\mathrm{\Σ}}\left({x_i-\stackrel{\‾}{x})}^2\right)(\underset{i=1...N}{\mathrm{\Σ}}{(y_i-\stackrel{\‾}{y})}^2\right)}}$

其中，κ₂表示所有线数据间的相似性的平均值；i＝1，...，N表示灰度邻域R_neib内的所有像素点；x_i，y_i分别表示不同边缘点邻域内的像素点的灰度；

表示的平均灰度。

19.根据权利要求16所述的方法，其特征在于，所述所有线数据间的相似性的平均值的计算采用如下公式：

${\mathrm{\κ}}_3=\frac{\underset{i=1...N}{\mathrm{\Σ}}\left(x_i{y}_i\right)}{\sqrt{\left(\underset{i=1...N}{\mathrm{\Σ}}{x}_i^2\right)\left(\underset{i=1...N}{\mathrm{\Σ}}{y}_i^2\right)}}$

其中，κ₃表示所有线数据间的相似性的平均值；i＝1，...，N表示灰度邻域R_neib内的所有像素点；x_i，y_i分别表示不同边缘点邻域内的像素点的灰度。

20.根据权利要求16所述的方法，其特征在于，所述步骤A包括：

确定灰度图像中的感兴趣区域；

采用基于灰度分布的能量泛函进行区域分割，将搜索得到的能量泛函极值对应的边缘位置作为相邻分割区域间的边界。

21.一种超声图像系统，其特征在于，所述图像计算模块包括：

区域初始化模块，用于确定灰度图像中的感兴趣区域；

组织区域边界自动检测模块，用于基于灰度分布的能量泛函进行区域分割，将搜索得到的能量泛函极值对应的边缘位置作为相邻分割区域间的边界；及

组织区域分析模块，用于根据所述边界位置计算组织区域的多项参数。

22.根据权利要求21所述的系统，其特征在于，所述组织区域分析模块包括：

内膜中层厚度计算单元，用于根据不同位置的所述边界位置计算各个位置的内膜中层厚度，并建立内膜中层厚度随位置变化的曲线予以输出显示；

置信度计算单元，用于利用边缘点邻域灰度的相似性计算边缘检测结果的置信度，予以输出显示。

23.根据权利要求22所述的系统，其特征在于，所述组织区域分析模块还包括：

方差统计单元，用于根据内膜中层厚度计算单元的输出结果，统计内膜中层厚度的分布方差，并叠加在所述内膜中层厚度的变化曲线中予以输出显示；

平均值统计单元，用于根据内膜中层厚度计算单元的输出结果，计算内膜中层厚度的平均值，并叠加在所述内膜中层厚度的变化曲线中予以输出显示。

Images