一种基于灰阶超声图像的心脏流场速度矢量分析方法

摘要

本发明公开了一种基于灰阶超声图像的心脏流场速度矢量分析方法，其特征在于，包括如下步骤：首先将灰阶超声图像划分为网格单元，在网格单元的每个单元格内选取像素最接近平均灰度值的像素点为控制点，并利用插值算法更新灰阶超声图像灰度值，达到去除斑点噪声的目的，再对去噪后的图像用光流法计算灰阶超声图像的心脏流场运动矢量。该方法所采用的插值算法过程不需要多次递归计算，可以实时进行斑点去除预处理；并且使运算结果更贴近真实情况，更加准确获取涡流等流体细节信息。

说明书

技术领域

本发明涉及医学图像处理及心脏流体力学研究领域，具体涉及一种基于灰 阶超声图像的心脏流场速度矢量分析方法。

背景技术

从心脏超声图像序列中同步提取血液和心肌的动态运动矢量包括两个重要 步骤：去除斑点噪声对图像的影响和运动矢量计算；

(1)超声图像去噪：医学噪声图像去噪的目的是在平滑噪声的同时尽可能 地保持原图像的重要特征。医学超声图像的噪声主要是斑点噪声，对此学术界 提出了一些抑制斑点噪声的方法为后续处理提供高质量图像。斑点噪声的空域 滤波包括均值滤波、中值滤波、Lee滤波、Kuan滤波、Frost滤波、Sigma滤波 及Gamma Map滤波等，这类算法无法克服的矛盾为：一方面增强斑点去噪效果 需要选较大的滤波窗口，另一方面为保持图像的实际分辨率要求所选的窗口较 小，而且窗口选择将大大地影响图像的处理结果。

近年来，出现了许多非线性滤波器，例如各向异性扩散滤波器，SRAD (Speckle reducing anisotropic diffusion)斑点噪声各向异性滤波器，及相干异性 扩散性滤波器等，这些滤波器由于在不同的扩散方向上采用不同的扩散系数， 因此具有增强对比度、保留细节等优点，成为滤波技术的一个亮点，然而需要 多次计算，严重影响速度。

(2)血流运动矢量计算：现有计算心脏超声的方法主要有PIV方法、斑点 跟踪方法、VFM方法和光流法。

PIV方法：PIV根据粒子浓度的高低分成两类，一类是适合于低浓度的 PTV(Particle Tracking Velocimery)方法，另一类是处理高浓度的粒子分布相关法， 这种方法是采用局部区域内粒子分布模型的相似性原理，使该区域内粒子同时 得到确认。其算法由四部分组成：图像粒子的标定、连续两帧图像中粒子对应 关系(采用互相关技术匹配)、粒子速度确定和误对应粒子速度的判定及其消除。

斑点跟踪方法：斑点跟踪一般包括两个方面，匹配和搜索。匹配的准则有 多种，常用的准则有归一化相关函数(NCCF)，平均均方误差函数(MSD)及平均 绝对差函数(MAD)等。搜索方法也有多种，全搜索法、十字搜索法和菱形搜索 法等。

VFM方法：VFM方法是应用到多普勒彩超计算运动矢量的方法，由于算法 先天性的角度依赖问题，使得临床应用很不稳定，而且不能处理灰度超声图像。

光流法：将光流法应用于心脏黑白超声图像研究极为少见，P.Burlina等人 在2011年将现有光流法应用于超声造影图像，但是效果并不是很理想，一方面 没有针对超声图像特点进行高效预处理，另一方面没有根据心脏流体特性(速 度快、流体复杂)对多尺度的权值进行优化。由于光流法与图像梯度关系密切 相关，没有经过平滑的图像几乎不能得到较为准确的运动矢量结果。同时，由 于血流速度较慢，而且血流流场复杂，直接将光流法应用于超声图像很难获得 较为准确的结果。

发明内容

本发明的目的在于：针对于现有技术的不足，提供一种基于灰阶超声图像 的心脏流场速度矢量分析方法，该方法首先对灰阶超声图像进行去除斑点噪声 处理，使心室内血流区域更加平滑，再采用光流法计算血流区域的运动矢量， 从而解决了基于灰阶超声图像的心脏流场速度矢量的分析问题。

为了达到上述发明目的，本发明采用的技术方案是：提供一种基于灰阶超 声图像的心脏流场速度矢量分析方法，其特征在于，包括以下步骤：

①将灰阶超声图像在横向和纵向方向等距离划分为网格单元；

②计算网格单元中每个单元格的平均灰度值P_mn，并选取单元格内最接近 平均灰度值的像素点作为单元格的控制点；其中，平均灰度值的的计算公式为：

$P_{\mathrm{mn}}=\frac{\underset{i,j}{\Sigma }{p}_{\mathrm{ij}}}{\mathrm{CellW}\times \mathrm{CellH}}$

其中，P_mn表示第m行n列单元格的平均灰度值；CellW和CellH分别为单 元格的宽度和高度；(m-1)·CellW≤i≤m·CellW，(n-1)·CellH≤j≤n·CellH；

③将步骤②中获取的控制点集进行网格剖分，在每个剖分的单元格内插值 计算所有像素的灰度值，并更新灰阶超声图像的灰度值，得到去除斑点噪声的 灰阶超声图像；

④将步骤③中得到的去除斑点噪声的灰阶超声图像通过高斯图像金字塔 分解方法分解成高度为n的高斯金字塔多分辨率图像；其中，n的计算公式为：

width和height分别是灰阶超声图像的宽度和高度；pyramid_space是下采样 因子，等于2；

⑤计算高斯金字塔多分辨率图像各层图像权重η_x，权重的计算公式为：

${\eta }_x=1-e^{-k{(x+t)}^2}$

其中，x＝0，1，2，…n-1；k＝0.05；t＝0.6；

⑥计算灰阶超声图像的心脏流场运动矢量V₀；

首先计算高斯金字塔多分辨率图像中第i层图像的运动矢量V_i：如果i＝n-1， 利用光流法计算第i层图像的初始运动矢量V′_i，且V_i＝V′_i；如果i＜n-1，则利 用光流法计算第i层图像的初始运动矢量V′_i，再根据V′_i计算第i层图像的运动 矢量V_i，计算公式如下：

V_i＝η_i+1·V_i+1+V′_i

其中，0≤i≤n-2；

重复上述步骤，直到i等于0，计算结果为V₀，即为灰阶超声图像的心脏流 场运动矢量。

优选地，步骤①中，在灰阶超声图像中横向和纵向方向上的划分间距均为7 个像素；所述网格单元为7×7的网格。

优选地，步骤②中，所述CellW为7，CellH为7。

优选地，步骤③中，采用Delaunay方法对控制点集进行三角形网格剖分， 在每个剖分的三角形单元格内采用分段三次Hermite方法插值计算所有像素的 灰度值。

优选地，步骤④中，将去除斑点噪声的灰阶超声图像进行高斯图像金字塔 分解的具体步骤为：

将去除斑点噪声的灰阶超声图像与高斯算子进行卷积，下采样，得到一个 分辨率为去除斑点噪声灰阶超声图像分辨率1/2的第一层图像；再将第一层图像 与高斯算子卷积，下采样，得到一个分辨率为去除斑点噪声灰阶超声图像分辨 率1/4的第二层图像，重复进行n次上述操作，直至得到高度为n的高斯金字 塔多分辨率图像。

优选地，所述高斯算子为5×5的模板，如下所示：

  1   1   2   1   1   1   2   4   2   1   2   4   8   4   2   1   2   4   2   1   1   1   2   1   1 。

综上所述，本发明提供的基于灰阶超声图像的心脏流场速度矢量分析方法 具有以下有益效果：

①本发明通过选取控制点并利用插值算法更新图像灰度值以达到快速去除 斑点噪声的目的，使心室内血流区域更加平滑；在去除斑点噪声过程中，所采 用的插值算法过程不需要多次递归计算，可以实时进行斑点去除预处理。

②本发明通过对去除斑点噪声的灰阶超声图像进行高斯金字塔分解，得到 高斯金字塔多分辨率图像，更便于准确计算该图像的运动矢量；由于低分辨率 的灰阶超声图像能够很好地反应总体血流流向，光流算法将低分辨率的运动矢 量赋予较大权重再进行计算，使运算结果更贴近真实情况，并且在表现总体运 动规律的同时，还获得了局部涡流信息。

附图说明

图1为本发明-实施例的流程图。

图2为本发明-实施例图像采用传统的光流法处理后的结果图。

图3为本发明-实施例图像采用本发明的方法处理后的结果图。

具体实施方式

下面结合具体实施例对本发明的具体实施方式做详细地描述：

实施例

参见附图1，本发明提供的基于灰阶超声图像的心脏流场速度矢量分析方法 包括以下步骤：

①将大小为520×512的灰阶超声图像(如图2和图3所示)在横向和纵向 方向上等间距划分为7×7的网格单元；其中，横向和纵向间距可以相同也可以 不同，在此，横向和纵向划分间距都为7个像素；

②计算网格单元中每个单元格的平均灰度值P_mn，并选取单元格内最接近平 均灰度值的像素点作为单元格的控制点；其中，平均灰度值的的计算公式为：

$P_{\mathrm{mn}}=\frac{\underset{i,j}{\Sigma }{p}_{\mathrm{ij}}}{\mathrm{CellW}\times \mathrm{CellH}}$

其中，P_mn表示第m行n列单元格的平均灰度值；CellW和CellH分别为单 元格的宽度和高度；(m-1)·CellW≤i≤m·CellW，(n-1)·CellH≤j≤n·CellH；

③采用Delaunay方法对步骤②中获取的控制点集进行三角形网格剖分，在 每个剖分的三角形单元格内采用分段三次Hermite方法插值计算所有像素的灰 度值，并更新灰阶超声图像的灰度值，得到去除斑点噪声的灰阶超声图像；

④将步骤③中得到的去除斑点噪声的灰阶超声图像进行高斯图像金字塔分 解，具体步骤为：

该灰阶超声图像的宽度为512，高度为520，根据公式

$n=\frac{\log \frac{\min (\mathrm{width},\mathrm{height})}{15}}{\log \mathrm{pyramid}\_\mathrm{space}}$

计算得到n＝5；再将去除斑点噪声的灰阶超声图像与高斯算子进行卷积，下 采样，得到一个分辨率为去除斑点噪声灰阶超声图像分辨率1/2的第一层图像； 再将第一层图像与高斯算子卷积，下采样，得到一个分辨率为去除斑点噪声灰 阶超声图像分辨率1/4的第二层图像，重复进行5次上述操作，直至得到高度 为5的高斯金字塔多分辨率图像；其中，所述高斯算子为5×5的模板，即

  1   1   2   1   1   1   2   4   2   1   2   4   8   4   2   1   2   4   2   1   1   1   2   1   1 ；

⑤计算由步骤④得到的高斯金字塔多分辨率图像各层图像权重η_x，权重 的计算公式为：

${\eta }_x=1-e^{-k{(x+t)}^2}$

其中，x＝0，1，2，3，4；k＝0.05；t＝0.6；由此可计算出，第0、1、2、3、 4层的运动矢量权重{η₀，η₁，η₂，η₃，η₄}分别为{0.9512，0.6376，0.2865，0.0863， 0.0174}；

⑥计算灰阶超声图像的心脏流场运动矢量V₀；

首先计算高斯金字塔多分辨率图像中第i层图像的运动矢量V_i：对于第i层 图像，将用光流法(Optical Flow)计算的运动矢量称为初始运动矢量，用V′_i表示；初始运动矢量V′_i与上一层运动矢量加权求和以后的矢量称为“运动矢量”， 用V_i表示；计算过程如下：

如果i＝n-1，利用光流法计算第i层图像的初始运动矢量V′_i，同时也是第i 层图像的最终运动矢量V_i；如果i＜n-1，则利用光流法计算第i层图像的初始运 动矢量V′_i，再根据V′_i计算第i层图像的运动矢量V_i，计算公式如下：

V_i＝η_i+1·V_i+1+V′_i

其中，0≤i≤n-2；

重复上述步骤，直到i等于0，计算结果为V₀，即为灰阶超声图像的心脏流 场运动矢量。

实验结果对照图

图2是采用传统的光流法处理后的心脏灰阶超声图像的区域截图；图3为 采用本发明的方法对心脏灰阶超声图像进行处理后的区域截图(即本发明实施 例所述)；从两个图像的对比可知，采用本发明的算法更能观察到流体细节涡流 信息，如图3黑色圆圈标注区域所示，它更能准确地表现血流流体湍流状态， 而这些流体特征区域也正是心脏功能研究的重点。

虽然结合具体实施例对本发明的具体实施方式进行了详细地描述，但并非 是对本专利保护范围的限定。在权利要求书所限定的范围内，本领域的技术人 员不经创造性劳动即可做出的各种修改或调整仍受本专利的保护。