CN103544682A - 一种三维超声图像非局部均值滤波方法 - Google Patents

Abstract

本发明公开了一种三维超声图像滤波方法，属于三维图像处理领域。本发明首先载入待滤波的三维超声图像，并将待滤波的目标三维超声图像拆分成基准块，给每个基准块指定搜索区域；然后将搜索区域中的所有与基准块大小相同的图像块作为相似块，并且计算每一相似块与基准块之间相似度；接下来以相似块的加权平均值作为基准块滤波结果，权重为相似度；最后将基准块的滤波结果进行整合得到三维超声图像最终滤波结果。与现有技术相比，本发明采用三维图像直接滤波，不拆分成二维图像帧滤波，利用帧与帧之间的灰度信息，能显著提高滤波效果并能更好的保留图像细节信息,同时本发明提供的滤波方法，适用于GPU加速，从而大大缩短滤波时间。

Description

一种三维超声图像非局部均值滤波方法

技术领域

本发明属于三维图像处理领域，更具体地，涉及一种三维超声图像非局部均值滤波方法。

背景技术

三维超声成像已取得了包括基础理论、成像技术、临床实验等许多方面的研究成果,在临床诊断中越来越受到重视。三维超声成像技术能够提供组织器官三维结构图像，与传统的二维超声相比图像显示直观，便于对医学诊断参数进行精确测量。然而，由于超声成像的相干特性，图像中存在斑点噪声，且在采集、转换和传输三维超声图像过程中，容易受到成像设备与外部环境噪声干扰等影响，导致成像质量下降。质量不佳的三维超声图像给边缘检测、图像分割及图像配准等后续处理工作带来了极大困难。因此，对三维超声图像进行滤波处理，是三维超声图像应用于医学诊断的重要步骤。

目前，三维超声图像滤波方法是将三维超声图像分解成若干帧二维超声图像，利用基于高斯模型的非局部均值滤波方法对其中每一帧进行滤波，再将滤波结果重建成三维图像作为三维超声图像滤波结果。这种方法没有有效利用帧与帧之间的灰度信息，滤波效果存在较为严重的失真。同时，由于三维超声图像的斑点噪声并不符合高斯噪声模型，因此基于高斯模型的非局部均值滤波法方法对于三维超声图像的滤波效果仍然不理想，且存在丢失细节信息较严重等问题，限制了三维超声图像在医学诊断中的应用。目前并没有能对三维超声图像进行直接滤波的非均值滤波方法。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种三维超声图像非局部均值滤波方法，其目的在于通过直接对三维体数据使用非局部均值滤波方法进行滤波，利用帧与帧之间的灰度信息，由此解决目前三维图像滤波方法滤波效果欠佳和丢失细节信息的技术问题。

为实现上述目的，按照本发明，提供了一种三维超声图像非局部均值滤波方法，包括以下步骤：

（1）载入待滤波的目标三维超声图像信息；

（2）将目标三维超声图像拆分为多个基准块；

（3）对于每一个基准块，分别在目标三维超声图像中确定其相应的搜索区域；

（4）对于每一个基准块，逐体素遍历该基准块相应搜索区域中与基准块大小相同的图像块，将这些图像块作为该基准块的相似块，计算基准块之间的相似度；

（5）对于每一个基准块，将其所有相似块的加权平均值作为该基准块的滤波结果，其中相似块的权重为其与基准块的相似度；

（6）将所有基准块滤波结果整合，得到目标三维超声图像最终滤波结果并输出。

优选地，所述三维超声图像非局部均值滤波方法，其相似块与基准块间相似度根据伽马模型模拟贝叶斯分布计算，相似度计算公式为：

$w(B_i,B_j)=\frac{1}{z_i}\exp (-\frac{1}{h^2}{\mathrm{\Σ}}_{p=1}^P\frac{{(u_{\left(B_i\right)}^{\left(p\right)}-u_{\left(B_j\right)}^{\left(p\right)})}^2}{u_{\left(B_j\right)}^{\left(p\right)}}),$

其中，B_i为中心是第i体素的基准块，B_j为中心是第j体素的相似块，w(B_i，B_j)为相似度，Z_i为归一化系数,h为滤波因子,P为图像块B_i或B_j的大小，

为基准块内第p个体素点的图像灰度值，

为相似块内第p个体素点的图像灰度值。

优选地，所述三维超声图像非局部均值滤波方法，其基准块为边长在5到11体素的三维立方体图像。

优选地，所述的三维超声图像非局部均值滤波方法，其相邻基准块中心之间的距离为3到6体素。

优选地，所述的三维超声图像非局部均值滤波方法，其搜索区域为以所述基准块中心为中心的三维立方体图像区域。

优选地，所述的三维超声图像非局部均值滤波方法，其搜索区域边长为11体素。

优选地，所述的三维超声图像非局部均值滤波方法，其相邻基准块其中心之间的距离小于等于基准块边长。

优选地，所述的三维超声图像非局部均值滤波方法，其步骤（6）将基准块滤波结果整合时，当一个体素处于多个基准块内时，将多个基准块中该体素滤波结果的平均值作为该体素最终滤波结果。

优选地，所述的三维超声图像非局部均值滤波方法，其步骤（3）到步骤（5）的并行计算过程使用GPU完成。

优选地，所述的三维超声图像非局部均值滤波方法，其GPU编程环境为CUDA语言。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，采用三维图像直接滤波，而不拆分为二维图像帧滤波，利用了帧与帧之间的灰度信息，避免了现有三维图像滤波方法滤波结果扭曲或者丢失细节严重的问题，能更好的过滤三维超声成像斑点噪声，同时减少了滤波结果失真。

本发明用伽马模型来模拟三维超声图像斑点噪声,根据伽马噪声模型和贝叶斯概率分布计算图像块与基准块之间的相似度，能够取得有效降低三维超声图像斑点噪声有益效果。并使用GPU完成基准块滤波，利用GPU并行计算能力，能大大缩短滤波时间。

对本发明提供的三维超声图像非局部均值滤波方法的各参数进行合理选择，如基准块大小、相邻基准块中心的距离、搜索区域大小等参数，能兼顾滤波效果和滤波计算量，在保证滤波结果质量的前提下缩短滤波时间。相邻基准块中心之间的距离不超过基准块的边长，能保证滤波结果覆盖目标三维超声图像的每一个体素。对于一个体素处在不同基准块中，将所有包括该体素的基准块灰度平均值作为该体素的滤波结果能保持滤波后的图像平滑，灰度过度自然。

附图说明

图1是本发明提供的三维超声图像非局部均值滤波方法流程图；

图2是待滤波的目标三维超声图像：

图3是基于高斯噪声模型的三维超声图像滤波方法和本发明提供的三维超声图像滤波方法搜索区域边长为11体素时滤波结果对比图，其中，图3（a）-1为基准块边长是5体素时基于高斯模型的滤波效果，图3（b）-1为基准块边长是5体素时基于伽马模型的滤波效果，图3（a）-2为基准块边长是7体素时基于高斯模型的滤波效果，图3（b）-2为基准块边长是7体素时基于伽马模型的滤波效果，图3（a）-3为基准块边长是9体素时基于高斯模型的滤波效果，图3（b）-3为基准块边长是9体素时基于伽马模型的滤波效果，图3（a）-4为基准块边长是11体素时基于高斯模型的滤波效果，图3（b）-4为基准块边长是11体素时基于伽马模型的滤波效果。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

对于给定的三维图像，尤其是三维超声图像，如图1所示，本发明提供的非局部均值滤波方法具体步骤如下：

（1）载入大小为W×H×F的三维图像信息。

（2）确定基准块边长和相邻基准块中心之间的距离，将目标三维超声图像拆分为多个基准块。基准块是待滤波的三维图像立方体块，通常是以某体素为中心的一块立方体三维图像区域，该体素为基准块中心，相邻基准块中心之间的距离，就是指一个基准块的中心和在X、Y、Z三个轴向上与其相邻的基准块的中心之间的距离。

优选方案，基准块边长在5到11体素之间，相邻基准块中心之间的距离一般不超过基准块边长。在目标三维图超声图像中选定一个体素，以该体素为中心的基准块，作为起始基准块，根据相邻基准块中心之间的距离确定目标三维超声图像中所有基准块，从而将目标三维超声图像拆为大小确定的基准块。相邻基准块中心之间的距离不超过基准块边长，使得所有基准块能全面覆盖目标三维超声图像，基准块之间允许有重叠。

（3）对于每一个基准块，分别确定相应的搜索区域。基准块与搜索区域的相对位置关系是固定的，一旦确定，在滤波过程中不能变更。舍弃其搜索区域溢出图像范围的基准块。

优选方案，搜索区域是以基准块中心为中心的三维立方体图像区域，边长为11体素。

（4）对于每一个基准块，逐体素遍历该基准块相应搜索区域中与基准块大小相同的图像块，将这些图像块作为该基准块的相似块，计算基准块之间的相似度。

超声图像的斑点噪声可采用高斯模型、伽马模型或其他数学模型来模拟，由于三维超声图像斑点噪声更符合伽马分布的噪声模型，本发明优选伽马模型模拟三维超声图像斑点噪声，模型数学表达式为：u(x)＝ν(x)+ν^γ(x)η(x)，

其中,ν(x)是理论上无噪声的三维图像,ν(x)是实际得到的三维超声图像,η(x)符合均值为0的高斯分布。实验表明，伽马模型较高斯模型能更好的模拟三维超声图像，且当γ=0.5的伽马分布噪声模型更适合于对数压缩的超声图像，此时伽马模型的数学表达式为：

$u\left(x\right)=\mathrm{\ν}\left(x\right)+\sqrt{\mathrm{\ν}\left(x\right)}\mathrm{\η}\left(x\right).$

对于基准块B_i，根据贝叶斯概率分布公式可定义如下：

$\hat{v}\left(B_i\right)=\frac{{\frac{1}{{\mathrm{\Δ}}_i}\mathrm{\Σ}}_{j=1}^{\left|{\mathrm{\Δ}}_i\right|}v\left(B_j\right)p\left(u\left(B_i\right)\right|v\left(B_j\right))}{\frac{1}{{\mathrm{\Δ}}_i}{\mathrm{\Σ}}_{j=1}^{\left|{\mathrm{\Δ}}_i\right|}p\left(u\left(B_i\right)\right|v\left(B_j\right))},$

其中，

是理论上图像块的B_i无噪声灰度值的估计，v(B_i)和v(B_j)分别表示无噪声图像块B_i和B_j,u(B_i)和u(B_j)分别表示以体素i和j为中心的实际得到的三维超声图像块，p(u(B_i)|v(B_j))是u(B_i)对无噪声的未知块v(B_j)的概率密度函数。因为v(B_j)是未知的，采用u(B_j)来代替v(B_j)。即得到：

$\hat{v}\left(B_i\right)=\frac{{\mathrm{\Σ}}_{j=1}^{\left|{\mathrm{\Δ}}_i\right|}u\left(B_j\right)p\left(u\left(B_i\right)\right|u\left(B_j\right))}{{\mathrm{\Σ}}_{j=1}^{\left|{\mathrm{\Δ}}_i\right|}p\left(u\left(B_i\right)\right|u\left(B_j\right))},$

假设:u(x)|ν(x)～N(ν(x)，ν(x)σ²),则:

$p\left(u\left(x\right)\right|\mathrm{\ν}\left(x\right))\∝\exp (-\frac{{(u\left(x\right)-\mathrm{\ν}\left(x\right))}^2}{2\mathrm{\ν}\left(x\right){\mathrm{\σ}}^2}).$

对于每个基准块B_i，其与相似块B_i的相似度的距离公式为：

$\mathrm{dp}(u\left(B_i\right),u\left(B_j\right))={\mathrm{\Σ}}_{p=1}^P\frac{{(u_{\left(B_i\right)}^{\left(p\right)}-u_{\left(B_j\right)}^{\left(p\right)})}^2}{u_{\left(B_j\right)}^{\left(p\right)}},$

其中，P代表块B_i的大小，即体素的个数。三维图像块相似块与基准块的相似度为：

$w(B_i,B_j)=\frac{1}{z_i}\exp (-\frac{1}{h^2}{\mathrm{\Σ}}_{p=1}^P\frac{{(u_{\left(B_i\right)}^{\left(p\right)}-u_{\left(B_j\right)}^{\left(p\right)})}^2}{u_{\left(B_j\right)}^{\left(p\right)}}),$

其中，B_i为中心是第i体素的基准块，B_j为中心是第j体素的相似块，w(B_i,B_j)为相似度，Z_i为归一化系数,h为滤波因子,P为图像块B_i或B_j的大小，

为基准块内第p个体素点的图像灰度值，为相似块内第p个体素点的图像灰度值。归一化系数的计算公式为：

$Z_i={\mathrm{\Σ}}_{j\∈S}\exp (-\frac{1}{h^2}{\mathrm{\Σ}}_{p=1}^P\frac{{(u_{\left(B_i\right)}^{\left(p\right)}-u_{\left(B_j\right)}^{\left(p\right)})}^2}{u_{\left(B_j\right)}^{\left(p\right)}}),$

其中，S为搜索区域，j为搜索区域内的体素。

（5）将基准块与相似块灰度的加权平均值作为基准块的滤波结果，其中基准块的权重为权重最大值，相似块的权重为与基准块的相似度。具体计算公式为：

NL(v)(B_i)＝Σ_j∈Iw(B_i，B_j)u(B_j)

其中，NL(v)(B_i)为滤波后的基准块，B_i为中心是第i体素的基准块，B_j为中心是第j体素的相似块，u(B_j)为以体素j为中心的相似块实际三维超声图像。

（5）将所有基准块的滤波结果整合，得到目标三维图像滤波结果，并输出。当基准块中心之间距离小于基准块边长时，基准块之间出现重叠，同一体素位于多个基准块中，将多个基准块中该体素滤波结果的平均值作为该体素最终滤波结果。

总体而言，本发明由于采用了伽马噪声模型，较之现有的三维图像的滤波方法，滤波效果显著提高，但由于所处理的对象是三维体数据，计算量非常庞大。本发明提供的三维图像滤波方法，其滤波效果和计算量受到以下因素影响：基准块大小、搜索区域大小、相邻基准块中心之间的距离。基准块大小越大，得到的滤波结果越光滑，丢失的细节信息越多，计算量越大，反之，基准块大小越小，得到的滤波结果越粗糙，保留的细节信息越多，计算量越小；搜索区域大小越大，滤波效果越好，计算量越大，反之，搜索区域大小越小，滤波效果越差，计算量越小；相邻基准块中心之间距离不超过基准块边长，能保证基准块全面覆盖目标三维图像，保证滤波效果，相邻基准块中心之间距离越小，滤波结果越光滑，计算量越大，反之，相邻基准块中心之间距离大，滤波结果越粗糙，计算量越小。

由于处理对象是三维体数据，计算量庞大，现有基于CPU的滤波方法处理能力相对较弱，本发明提供的三维图像非局部均值滤波方法通过选择上述参数，兼顾滤波效果和滤波时间。另外，优选方案，采用GPU并行处理，加速滤波过程，缩短滤波时间，具体步骤如下：

（1）为目标三维超声图像数据分配内存，并将数据读入内存；

（2）初始化GPU计算环境，为GPU分配显存空间用于内核函数输入和输出；

（3）将三维超声图像数据由内存传送到步骤（2）中分配的显存空间；

（4）根据基准块大小、搜索区域大小、相邻基准块中心之间的距离分配线程结构；

（5）调用GPU进行并行计算各基准块的滤波结果，并将滤波结果存储在显存中；

（6）将步骤（5）中存储在显存中的基准块滤波结果传送到内存空间；

（7）调用CPU，将内存中基准块滤波结果进行整合，得到最终滤波结果，将最终滤波结果输出到文件，退出GPU计算环境，并将计算过程中分配的内存和显存释放。

以下为实施例：

图2是一个胎儿的三维超声图像，大小为428×354×209（体素），该图像较为模糊，噪点多，难以辨认，对该图像使用本发明的方法进行滤波处理。所使用的计算系统为：CPU，Intel(R)Core(TM)i3-21203.3GHz；GPU型号是NVIDIA GeForce GTX660Ti，操作系统是32位Windows7，编程环境CUDA语言。

当基准块大小为5×5×5体素到11×11×11体素之间时，基于高斯模型的滤波效果为如图3(a)-1至3(a)-4,本发明提供的三维超声滤波方法其效果为图3(b)-1至3(b)-4。从图3中可看出，在条件相同时，基于伽马模型的非均值滤波效果明显优于基于高斯模型的非均值滤波效果。

利用GPU的并行处理能力进行加速的效果见表1。表1的结果显示，利用GPU能显著缩短算法运行时间，在计算密度越高的条件下效果越发显著。由于采用GPU加速，使本发明提供的滤波方法在实际应用于三维超声体数据去噪时，其运行时间可以被接受。

表1 GPU并行处理加速结果

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

Claims (10)

1.一种三维超声图像非局部均值滤波方法，其特征在于，包括以下步骤：

（1）载入待滤波的目标三维超声图像信息；

（2）将目标三维超声图像拆分为多个基准块；

（3）对于每一个基准块，分别在目标三维超声图像中确定其相应的搜索区域；

（4）对于每一个基准块，逐体素遍历该基准块相应搜索区域中与基准块大小相同的图像块，将这些图像块作为该基准块的相似块，计算这些图像块与所述基准块之间的相似度；

（5）对于每一个基准块，将其所有相似块的加权平均值作为该基准块的滤波结果，其中相似块的权重为其与基准块的相似度；

（6）将所有基准块滤波结果整合，得到目标三维超声图像最终滤波结果并输出。

2.如权利要求1所述的三维超声图像非局部均值滤波方法，其特征在于，所述相似块与基准块间相似度根据伽马模型模拟贝叶斯分布计算，相似度计算公式为：

$w(B_i,B_j)=\frac{1}{z_i}\exp (-\frac{1}{h^2}{\mathrm{\Σ}}_{p=1}^P\frac{{(u_{\left(B_i\right)}^{\left(p\right)}-u_{\left(B_j\right)}^{\left(p\right)})}^2}{u_{\left(B_j\right)}^{\left(p\right)}}),$

其中，B_i为中心是第i体素的基准块，B_j为中心是第j体素的相似块，w(B_i,B_j)为相似度，Z_i为归一化系数,h为滤波因子,P为图像块B_i或B_j的大小，

为基准块内第p个体素点的图像灰度值，

为相似块内第p个体素点的图像灰度值。

3.如权利要求1或2所述的三维超声图像非局部均值滤波方法，其特征在于，所述基准块为边长在5到11体素的三维立方体图像。

4.如权利要求1至3中任意一项所述的三维超声图像非局部均值滤波方法，其特征在于，所述相邻基准块中心之间的距离为3到6体素。

5.如权利要求1至4中任意一项所述的三维超声图像非局部均值滤波方法，其特征在于，所述搜索区域为以所述基准块中心为中心的三维立方体图像区域。

6.如权利要求1至5中任意一项所述的三维超声图像非局部均值滤波方法，其特征在于，所述搜索区域边长为11体素。

7.如权利要求1到6中任意一项所述的三维超声图像非局部均值滤波方法，其特征在于，所述相邻基准块其中心之间的距离小于等于基准块边长。

8.如权利要求7所述的三维超声图像非局部均值滤波方法，其特征在于，步骤（6）将基准块滤波结果整合时，当一个体素处于多个基准块内时，将多个基准块中该体素滤波结果的平均值作为该体素最终滤波结果。

9.如权利要求1到8中任意一项所述的三维超声图像非局部均值滤波方法，其特征在于，所述步骤（3）到步骤（5）的并行计算过程使用GPU完成。

10.如权利要求9所述的三维超声图像非局部均值滤波方法，其特征在于，所述GPU编程环境为CUDA语言。

Images