CN103745447A - 一种非局部均值滤波的快速并行实现方法 - Google Patents

Abstract

非局部均值滤波算法是一种被广泛应用于抑制图像噪声的算法，该算法根据属于同一图像结构的像素点周围有着相似的邻域结构这一假设，并基于此邻域相似性质构建加权滤波器来抑制图像中的噪声。实验证明该非局部均值滤波算法能在有效地抑制图像中的噪声的同时保持图像中的组织信息，为了有效地抑制图像中的噪声，一般需要较大的搜索窗以引入较多的邻域信息，从而导致大量的计算量和处理时间，影响了其在实际中的应用。为了解决这一问题，本发明提出了一种非局部均值滤波的快速并行实现方法，该方法在原有的以像素为单位的GPU并行的基础上，利用共享存储器特性及非局部均值权重对称性来优化并行操作，显著地提高了非局部均值滤波算法的计算速度。

Description

一种非局部均值滤波的快速并行实现方法

技术领域

本发明涉及一种非局部均值滤波算法在GPU上的快速并行实现方法。

背景技术

图像降噪始终是数字图像处理领域中的一个重要的研究内容，经典的降噪滤波方法有邻域平均值法、中值法以及一些频域滤波方法，这些图像降噪算法一般基于像素的灰度差和梯度等信息，只利用到较小邻域的信息，易导致结构模糊的图像处理结果。而Buades基于从图像中任取一个小窗口，都能够从该图像的一个较大范围内中找到许多与其相似的窗口结构的事实提出了非局部均值滤波算法，这种算法可以充分利用图像中更大范围内的图像信息对噪声抑制，从而能够在不丢失图像细节的前提下有效的抑制图像中的噪声。具体的，非局部均值滤波算法把每个像素替换成其邻近像素乘以权重之后的平均值并利用两个块之间的相似度计算权重，即假设现在处理像素点p(p＝(p_x,p_y))，p搜索窗内的像素点q的权重值等于分别以p，q为中心的比较块进行比较后得到的值，权重值与两个比较块的相似度呈正相关关系。认为X是目标处理图像，Y是待处理的图像，非局部均值滤波算法可用如下公式表示：

$\hat{X}\left(p\right)=\frac{{\mathrm{\Σ}}_{q\∈N_p}\mathrm{\ω}(p,q)Y\left(q\right)}{{\mathrm{\Σ}}_{q\∈N_p}\mathrm{\ω}(p,q)}---\left(1\right)$

$w(p,q)=\exp (-\frac{\underset{(\mathrm{\Δx},\mathrm{\Δy})\∈{[-B,...,B]}^2}{\mathrm{\Σ}}\left|d_{p,q}^{(\mathrm{\Δx},\mathrm{\Δy})}\right|G(\mathrm{\Δx},\mathrm{\Δy})}{h(2B+1)2B+1)})---\left(2\right)$

$d_{p,q}^{(\mathrm{\Δx},\mathrm{\Δy})}=Y(p_x+\mathrm{\Δx},p_y+\mathrm{\Δy})-Y(q_x+\mathrm{\Δx},q_y+\mathrm{\Δy})---\left(3\right)$

$G(\mathrm{\Δx},\mathrm{\Δy})=\left\{\begin{array}{cc}\sqrt{2}& (\mathrm{\Δx},\mathrm{\Δy})=\left(\mathrm{0,0}\right)\\ \frac{1}{\sqrt{{\mathrm{\Δx}}^2+{\mathrm{\Δy}}^2}}& \mathrm{other}\end{array}\right.---\left(4\right)$

这里，Y和

分别表示处理前和处理后的图像，N_p是以p为中心点的搜索窗。w(p,q)表示以p,q为中心点的半径为B两个比较块之间的相似度，G(Δx,Δy)是一个与距离相关的高斯核函数。这样(2B+1)(2B+1)即为比较块B中点的个数，我们可以用公式(2)中的参数h用来控制算法处理的平滑效果。

然而，由于非局部均值滤波算法使用较大的搜索窗以引入较多的邻域图像信息来实现对噪声有效抑制，这同时也带来了对计算量的较大需求，影响了算法在实际中的应用。为了使非局部均值滤波算法更具实用性，需要对这种算法进行加速。

如今，利用并行处理技术来进行算法加速已经成为了一种趋势，而利用GPU进行加速是并行处理技术较为常见的一种。GPU在处理能力和存储器带宽上相较于CPU具有明显的优势，在单精度浮点处理能力上也远远超过CPU，GPU的并行主要是通过粗粒度的block和细粒度的thread联合并行来实现。NVIDIA推出的CUDA是一种将GPU作为数据并行设备的软硬件体系，它是一种使用类C语言进行通用计算的开发环境和软件体系。CUDA为开发人员有效利用GPU的强大性能提供了便利的条件，它被广泛应用于金融、石油、天文学、图像处理等领域。

目前存在着很多基于GPU并行的非局部均值滤波的加速算法，较为经典有效的非局部均值滤波算法的GPU加速算法可阐述如下：

从公式(1)-(4)可以看到，我们可以直接以像素为单位对非局部均值滤波算法进行并行化的GPU加速。出于GPU的每个内核函数所拥有的共享内存、寄存器等GPU硬件限制的考虑，我们可以把这个算法拆分三部分进行循环计算，循环次数为搜索窗的大小|N|(即为搜索窗内像素点个数)，(p_x+i_x,p_y+i_y)是以p为中心点的搜索窗中的某个点的位置，初始化

第一个内核函数计算比较块的像素差异值，可用如下公式表示，它的计算复杂度是O(1)。

$[U_1^{(3i-1)},...,U_4^{(3i-1)}]\left(p\right)=f_{U_1^{(3i-2)}}^{(3i-2)},...,U_4^{(3i-2)}\left(p\right)=\left(\begin{array}{c}|Y(p_x,p_y)-Y(p_x+i_x,p_y+i_y)|\\ U_2^{(3i-2)}\left(p\right)\\ U_3^{(3i-2)}\left(p\right)\\ U_4^{(3i-2)}\left(p\right)\end{array}\right)---\left(5\right)$

第二个内核函数根据第一个内核函数计算出来后的比较块的对应像素差异绝对值根据公式(2)进行比较块相似度的计算，它的计算复杂度是O((2B+1)(2B+1))。

$[U_1^{\left(3i\right)},...,U_4^{\left(3i\right)}]\left(p\right)=f_{U_1^{(3i-1)},...,U_4^{(3i-1)}}^{(3i-1)}\left(p\right)=\left(\begin{array}{c}{U}_1^{(3i-1)}\\ \exp (-\frac{\underset{(\mathrm{\Δx},\mathrm{\Δy})\∈{[-B,...,B]}^2}{\mathrm{\Σ}}{U}_1^{(3i-1)}(p_x+\mathrm{\Δx},p_y+\mathrm{\Δy})G(\mathrm{\Δx},\mathrm{\Δy})}{h(2B+1)(2B+1)})\\ U_3^{(3i-1)}\\ U_4^{(3i-1)}\end{array}\right)---\left(6\right)$

第三个内核函数是用来累加权重和像素和，它的计算复杂度是O(1)。

$[U_1^{(3(i+1)-2)},...,U_4^{(3(i+1)-2)}]\left(p\right)=f_{U_1^{\left(3i\right)},...,U_4^{\left(3i\right)}}^{\left(3i\right)}\left(p\right)=\left(\begin{array}{c}{U}_1^{\left(3i\right)}\left(p\right)\\ U_2^{\left(3i\right)}\left(p\right)\\ U_3^{\left(3i\right)}\left(p\right)+U_2^{\left(3i\right)}\left(p\right)\\ U_4^{\left(3i\right)}\left(p\right)+U_2^{\left(3i\right)}\left(p\right)Y(p_x+i_x,p_y+i_y)\end{array}\right)---\left(7\right)$

此外还有一个内核函数来计算最后的输出图像，表示如下：

$f_{U_1^{\left(I\right)},...,U_4^{\left(I\right)}}^{\left(I\right)}\left(p\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ \frac{U_4^I\left(p\right)}{U_3^I\left(p\right)}\end{array}\right)---\left(8\right)$

此时I＝|N|+1，I可以用来表示数据U₃和U₄针对搜索窗每个元素的内核函数计算的最终循环次数，计算复杂度是O(1)。最终的输出图像为

综上可知，该普通的GPU加速算法算法的计算复杂度是O(|N|((2B+1)(2B+1)+2)+1)。

发明内容

本发明提出了一种非局部均值滤波的快速并行实现方法，该方法在不改变原先算法处理效果的前提下明显提高计算处理速度，具体阐述如下：

在普通的非局部均值滤波的GPU加速算法中，第二个内核函数需要用传统串行计算的模式计算两个比较块的相似度，这样当使用较大的搜索窗处理图像时需要较大的计算量，所以我们的方法的第一个改进点是降低这一部分的计算复杂度。我们首先分析一下两个比较块相似度是如何计算的。我们假设比较块是一个5×5的块（B＝2），其中p点表示中心点，如图1所示。由上面的公式(1)-(4)可知,两个比较块的相似度是通过计算相对应位置的像素的差异值，由于越靠近中心点的像素对相似度的影响越大，所以我们对这些计算出来的差异值乘以与距离相关的系数G(Δx,Δy)（如公式(4)所示）。从图1中我们可以发现，如果是以q点为比较块的中心点，那么图中箭头联系的像素点的像素差异值跟以p点为中心点时是一样的，不同的仅是它们的距离系数。基于这样的思想，我们提出了公式(9)：

$\underset{\mathrm{\Δy}\∈[0,...,B]}{\∪}\underset{\mathrm{\Δx}\∈[-B,...,B]}{\mathrm{\Σ}}|Y(p_x+\mathrm{\Δx},p_y)-Y(q_x+\mathrm{\Δx},q_y)|\×G(\mathrm{\Δx},\mathrm{\Δy})---\left(9\right)$

其中公式中计算两个像素值的差异值为|Y(p_x+Δx,p_y)-Y(q_x+Δx,q_y)|，q＝(q_x,q_y)用来表示搜索窗中的某点。

表示计算垂直方向上不同的中心点，这样如果使用(2B+1)×(2B+1)大小的比较块，基于公式(10)的每一行的计算会有(B+1)个不同结果。公式(11)虽然看起来计算复杂度是O((B+1)(2B+1))，但是由于每行的差异计算结果只要计算一次并放在共享存储器中，其他B次都从共享存储器中读取数据,对于具有较强单精度浮点处理能力的GPU来说，公式(10)的计算开销主要是来自访问显存所花费的时间，考虑到共享存储器的数据存取速度远高于显存数据存取速度，这个公式的计算复杂度可以近似的表示为O(2B+1)。

和普通的GPU加速算法一样，我们同样的把这个算法拆分成三部分进行循环计算。假设输入图像是m×n，设定的大小是m×n×(B+1)，的大小都为m×n，其中i=1,...,|N|，(p_x+i_x,p_y+i_y)为以p为中心点的搜索窗中的某个点的位置，初始化 $U_3^{\left(1\right)}=U_4^{\left(1\right)}=0,$ $U_5^{\left(i\right)}=Y.$

第一个内核函数，用来计算比较块中像素点灰度差异值，同时以比较块中的每一行为单位，根据每一行中的每个像素点与不同中心点的距离乘以不同的距离系数得到每一行所有的可能累加灰度差绝对值，该内核函数的计算复杂度等于O(2B+1)。

$[U_1^{(3i-1)},...,U_5^{(3i-1)}]\left(p\right)=$

$f_{U_1^{(3i-2)},...,U_5^{(3i-2)}}^{(3i-2)}\left(p\right)=\left(\begin{array}{c}\underset{\mathrm{\Δy}\∈[0,...,B]}{\∪}\underset{\mathrm{\Δx}\∈[-B,...,B]}{\mathrm{\Σ}}|Y(p_x+\mathrm{\Δx},p_y)-Y(p_x+i_x+\mathrm{\Δx},p_y+i_y)|G(\mathrm{\Δx},\mathrm{\Δy})\\ U_2^{(3i-2)}\left(p\right)\\ U_3^{(3i-2)}\left(p\right)\\ U_4^{(3i-2)}\left(p\right)\\ U_5^{(3i-2)}\left(p\right)\end{array}\right)---\left(10\right)$

第二个内核函数计算比较块的相似度，即对第一个内核函数的结果，选择相应的累加灰度差绝对值进行累加，并根据公式(2)计算得出两个比较块的相似度，该内核函数的计算复杂度为O(2B+1)。

$[U_1^{\left(3i\right)},...,U_5^{\left(3i\right)}]\left(p\right)=f_{U_1^{(3i-1)},...,U_5^{(3i-1)}}^{(3i-1)}\left(p\right)=\left(\begin{array}{c}{U}_1^{(3i-1)}\left(p\right)\\ \exp (-\frac{\underset{\mathrm{\Δy}\∈[-B,...,B]}{\mathrm{\Σ}}{U}_1^{(3i-1)}(p_x,p_y,\mathrm{\Δy})}{h(2B+1)(2B+1)})\\ U_3^{(3i-1)}\left(p\right)\\ U_4^{(3i-1)}\left(p\right)\\ U_5^{(3i-1)}\left(p\right)\end{array}\right)---\left(11\right)$

第二个改进点是利用权重计算的对称性，很显然w(p,p+Δp)＝w(p+Δq,p)(Δq表示像素p在搜索窗N中的偏移量)。我们利用权重对称性这一事实，当对p位置累加像素和w(p,p+Δq)Y(p+Δq)，我们同时也可以累加w(p-Δq,p)Y(p-Δq)(根据对称性w(p-Δq,p)Y(p-Δq)＝w(p,p-Δq)Y(p-Δq))，这样我们只需要遍历原来搜索窗里的点的一半就可以了，第三个内核函数可用如下公式(12)表示，它的计算复杂度是O(1)。

$[U_1^{(3(i+1)-2)},...,U_5^{(3(i+1)-2)}]\left(p\right)=f_{U_1^{\left(3i\right)},...,U_5^{\left(3i\right)}}^{\left(3i\right)}\left(p\right)=$

$\left(\begin{array}{c}{U}_1^{\left(3i\right)}\left(p\right)\\ U_2^{\left(3i\right)}\left(p\right)\\ U_3^{\left(3i\right)}\left(p\right)+U_2^{\left(3i\right)}(p+\mathrm{\Δq})+U_2^{\left(3i\right)}(p-\mathrm{\Δq})\\ U_4^{\left(3i\right)}\left(p\right)+U_2^{\left(3i\right)}(p+\mathrm{\Δq})U_5^{\left(3i\right)}(p+\mathrm{\Δq})+U_2^{\left(3i\right)}(p-\mathrm{\Δq})U_5^{\left(3i\right)}(p-\mathrm{\Δq})\\ U_5^{\left(3i\right)}\end{array}\right)---\left(12\right)$

像普通的非局部均值滤波算法加速一样，我们同样的通过一个最终的内核函数来计算最后的输出处理图像：

$f_{U_1^{\left(I\right)},...,U_5^{\left(I\right)}}^{\left(I\right)}\left(p\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ \frac{U_4^{\left(I\right)\left(p\right)}}{U_3^{\left(I\right)}\left(p\right)}\end{array}\right)---\left(13\right)$

这里表示数据U₃和U₄针对搜索窗每个元素的内核函数计算的最终循环次数I等于(2T+1)×(T+1)+1(T为搜索窗半径)。最终输出图像为

有益效果：我们的改进算法的累计计算复杂度为O(((2T+1)×(T+1)+1)(2(2B+1)+1))+1)，考虑到((2T+1)×(T+1)+1)近似等于0.5|N|，故此累计计算复杂度近似等于O(|N|(2B+1)+1)，我们可以看到相对于普通非局部均值滤波GPU加速算法，此改进的算法能将计算复杂度大致缩短为原来的

倍，实现了在不改变原先算法处理效果的前提下明显提高计算处理速度。

附图说明

图1为比较块的相似度计算；

图2为基于非局部均值滤波的CPU串行算法和普通GPU并行算法的计算时间对比，这里处理的低剂量CT图像大小为512×512，比较块大小为9×9；

图3为基于非局部均值滤波的普通的和改进后的GPU加速算法的运算时间对比，这里处理的低剂量CT图像大小为512×512，比较块大小为9×9；

图4为基于非局部均值滤波的的普通的和改进后的GPU加速算法的运算时间对比，这里处理的低剂量CT图像大小为512×512，搜索窗大小为81×81。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

基于改进GPU并行的快速非局部均值滤波算法，包括以下步骤：

步骤1、在GPU中，每一个线程都计算它所对应的像素点与其搜索窗中的某一个位置的像素点的灰度差异值的绝对值。当所有的线程都计算完差异值后，计算以该线程所对应的像素点为中心的比较块（假设半径为B）的中心行的B+1种可能灰度累加值，即根据该中心行在不同的比较块中与该比较块中心点的距离乘以不同的距离系数得到该中心行所有的可能累加灰度差绝对值；

步骤2、在GPU中，每一个线程都计算它所对应的像素点的比较块与其的搜索窗中的某一个位置（这个位置是与步骤1相同的）的像素点的比较块的相似度，即对步骤1的结果，选择相应的累加灰度差绝对值进行累加，根据累加结果计算比较块的相似度，得到一个权重值；

步骤3、在GPU中，每一个线程都累加它所对应的像素点在步骤2计算出来的权重值，同时也累加权重值乘以像素的值；

循环遍历搜索窗里的像素点，每一个像素点都执行上述三个步骤；

遍历搜索窗里的所有位置，每一次都执行上述三个步骤；

步骤4、根据最后一次步骤3得出的累加权重值和累加像素和计算出处理后的像素值。

我们的GPU的并行是每一个线程对应于处理图像中的一个像素点，上述的四个步骤对应四个不同的内核函数，我们现结合具体的例子阐述这四个内核函数所做的工作。现以线程对应的像素点为p，计算p的搜索窗里的像素点q(q＝(p_x+i_x,p_y+i_y))与p的相似度，比较块的半径为B为例。

在步骤1中，我们首先计算图像中的像素点p与q的像素点差异值，取绝对值，放在共享存储器里，等待所有的线程都处理完毕。然后计算以p为中心点，半径为B的行的差异值累加值，根据p与比较块的中心点的垂直距离(共有B+1种不同)乘以不同的距离系数得到每一行所有的可能累加灰度差绝对值。

在步骤2中，我们计算像素点p与像素点q的相似度，首先根据步骤1的的结果，累加p所在的比较块与q的比较块的像素差异值，然后根据公式(14)计算q对p的权重wieght，其中sum表示p与q的比较块的像素差异值的累加值。

$\mathrm{weight}=\exp (-\frac{\mathrm{sum}}{h(2B+1)(2B+1)})---\left(14\right)$

在步骤3中，我们累加步骤2得到的q的权重值到p所对应的U₃(U₃存放p的搜索窗中的所有点的权重值)，同时累加q的权重值乘以q的像素到p所对应的U₄(U₄存放p的搜索窗中的所有点的累加像素和)，同时我们利用权重计算的对称性，累加p搜索窗中的像素点s(s＝(p_x-i_x,p_y-i_y))。

遍历p的搜索窗里的所有点，每一次都执行上述三个步骤；

在步骤4中，我们利用根据公式计算出处理后的像素点p的值，

5．效果评估准则

在同一台机器上比较非局部均值滤波算法在CPU上的串行时间以及使用两种GPU并行加速后的运算时间（包括普通的和改进的），所使用的实验计算机环境的配置参数如下：

1）硬件：

CPU：Inter(R)Core(TM)i7-3770CPU3.40GHz

内存：8GB

显卡：NVIDIAGeForceGTX680，其中流处理器：1536个，显存频率：6008MHz，显存带宽：192GB/S，显存容量：2GB，显存位宽：256bit

2）软件

操作系统：Win764位

Matlab：R2011a

CUDA：4.0

5.1视觉评估

通过对一幅图像利用三种不同的非局部均值滤波算法现实方式（CPU串行、普通GPU并行、改进GPU并行）得到的处理结果一致，利用CPU串行实现的非局部均值滤波算法运行时间非常长，普通GPU并行大幅度地减低非局部均值滤波算法的运行时间，而改进的GPU并行实现方式相比于普通GPU并行方式又再一次地减低了算法的运行时间。

5.2量化评估

为了量化的验证本发明方法的有效性，我们分别采用CPU串行方法，普通GPU加速算法，改进的GPU加速算法对一副512×512的低剂量CT图像进行处理，计算的并行化没有改变算法的处理原理，不同的方法能够得到同样的处理结果。

首先我们在不同搜索窗尺寸下，固定比较块大小为9×9，对比CPU串行算法和普通GPU算法的计算时间。实验结果如图2所示。从图中我们可以看到基于GPU的并行计算能够大幅度降低运算时间，相对于原来的基于CPU的串行算法能够获得超过一百倍的加速。

接下来，我们继续比较在不同尺寸的搜索窗设置下，普通GPU并行算法和改进后的GPU加速算法的运行时间，比较块大小为9×9，实验结果如图3所示。从图3我们可以发现当搜索窗尺寸比较大时，加速倍数比较接近等于2B+1=2×4+1=9。

最后，我们再比较一下当图像尺寸为512×512，搜索窗大小为81×81，比较块大小改变时，普通GPU加速算法和改进后的GPU加速算法的运行时间，实验结果如图4所示。由于当比较块较大时，会导致图像模糊，体现不出非局部均值滤波算法的去噪效果，所以在这里我们比较块最大就取到15×15。从图4中我们很清楚可以看到加速比对着比较快的增大而增加，满足前面分析得出的2B+1数值。

Claims (3)

1.一种非局部均值滤波的快速并行实现方法，其特征在于，以像素为单位，每一个图形处理单元GPU中的线程计算图像中的一个像素点，对非局部均值滤波算法进行并行化的GPU加速，包括以下步骤：

步骤1、在GPU中，每一个线程都计算它所对应的像素点与其搜索窗中的某一个位置的像素点的灰度差异值的绝对值；当所有的线程都计算完差异值后，假设比较块变径为B，计算以该线程所对应的像素点为中心的比较块的中心行的B+1种可能灰度累加值，即根据该中心行在不同的比较块中与该比较块中心点的距离乘以不同的距离系数得到该中心行所有的可能累加灰度差绝对值；

步骤2、在GPU中，每一个线程都计算它所对应的像素点的比较块与其的搜索窗中的步骤1中的某一个位置的像素点的比较块的相似度，即对步骤1的结果，选择相应的累加灰度差绝对值进行累加，根据累加结果计算比较块的相似度，得到一个权重值；

步骤3、在GPU中，每一个线程都累加它所对应的像素点在步骤2计算的权重值，同时也累加权重值乘以像素的值；

遍历搜索窗里的所有位置，每一次都执行上述三个步骤；

步骤4、根据最后一次步骤3得出的累加权重值和累加像素和计算出处理后的像素值。

2.如权利要求1所述的方法，其特征在于，在步骤1中，把计算出来的像素点灰度差异值保存在共享存储器中，这样可以减少访问显存的次数，由于访问共享存储器的速度远高于访问显存的速度，因此这样可以大幅降低数据传输时间；然后根据比较块的中心行在不同的比较块中与该比较块中心点的距离乘以不同的距离系数得到该中心行所有的可能累加灰度差绝对值。

3.如权利要求1所述的方法，其特征在于，在步骤3中，利用权重计算的对称性，当对搜索窗中的某一位置累加时，同时也累加该位置的对称位置。

Images