一种基于GPU并行加速的三维块匹配降噪方法与流程

本发明涉及图像处理技术领域，具体公开了一种基于gpu并行加速的三维块匹配降噪方法。

背景技术：

数字图像在获取和传输过程中由于受到成像设备和外部环境的影响常常含有大量噪声，影响图像质量。尤其是在医学应用上，由于受低剂量电压电流的影响，ct图像中往往会引入大量的噪声，致使图像质量退化，影响医生的临床病症诊断。

常见的三维块匹配(bm3d)算法凭借其自身结合局部、非局部、多尺度稀疏、自适应的滤波降噪特点，被认为是目前最好的图像去噪算法，但是该算法基于相似图像块协同滤波，算法复杂度高，计算量大，处理大尺寸的ct图像数据时，耗时较长，处理效率较低，无法满足实际工作需求。

因此，需要一种能解决上述问题的方法。

技术实现要素：

为了克服现有技术中存在的缺点和不足，本发明的目的在于提供一种基于gpu并行加速的三维块匹配降噪方法。

为实现上述目的，本发明采用如下方案。

一种基于gpu并行加速的三维块匹配降噪方法，包括：

在cpu端将待处理的图像进行边界对称扩展预处理；

将预处理后的图像数据从cpu主机端发送至gpu的全局存储器中；

创建线程网络grid，采用全局存储器合并访问的模式及共享存储器多次循环利用的加速策略，对相似图像块匹配分组进行并行加速处理；

采用硬阈值协同滤波内核函数并行加速策略获取三维相似矩阵的第一步降噪估计数据；

以第一步降噪估计数据为参考，联合维纳协同滤波内核函数并行加速策略获取第二步降噪估计数据；

将第二步降噪估计数据从gpu发送至cpu主机端，并剔除扩展边界像素以获取降噪后的图像。

进一步地，所述创建线程网络grid，包括：

以图像每个参考块的图像块匹配过程为一个线程块block，搜索窗内每个检测图像块与参考块的相似性匹配过程为每个线程thread进行线程分配；

以一定像素的步长分别从行列方向依次递增选取参考图像块，并根据图像中参考图像块个数确定线程网grid的大小，根据参考图像块搜索窗的图像块个数确定线程block的大小。

进一步地，，所述采用全局存储器合并访问的模式及共享存储器多次循环利用的加速策略，对相似图像块匹配分组进行并行加速处理，包括：

将同一个warp中的所有线程thread执行同一条指令访问全局存储器中连续单元，以获得合并访问模式；

将搜索窗分成大小为32*32的4个分块，并在每个分块中以block(16,16)的线程进行相似性计算获取相似块，其中d为图像块间的距离，定义为两个图像块对应元素数值差的模除以图像块的大小，τthreod为选取的适合的距离阈值；

将搜索窗中像元数据分块循环共享存储器中，并设置threadidx.x<16且threaidx.y<16；

采用最小值归约并行的策略寻找规定数目的最相似图像块。

进一步地，，所述采用最小值归约并行的策略寻找规定数目的最相似图像块，包括：

将参考图像块的相似图像块按照相似距离由小到大排序聚集成三维矩阵；

启用n个线程分别与相似性计算得到的n个距离值d[n]对应；

将第i个线程的值与第(i+n/2)个线程的距离值进行比较并将较小值放在左部分中，较大值放在右部分，则左部分区间为d[0]至d[n/2]，右部分区间为d[n/2]至d[n]；

完成线程并行比较后，将比较线程数减半，对左部分区间距离值进行上述比较，直至多次比较线程数减半比较后左部分区间值缩至d[0]；

以d[0]作为距离值中的最小值，将d[n]起始访问位置后移一位并重复上述步骤寻找最小值，直至找到距离最小的规定数目的相似图像块。

进一步地，，所述采用硬阈值协同滤波内核函数并行加速策略获取三维相似矩阵的第一步降噪估计数据，包括：

指令混合优化加速，将三维相似矩阵的三维正变换、硬阈值滤波、三维逆变换以及加权估计的过程集成于硬阈值协同滤波内核函数中；其中所述三维正变换包括依次进行二维双正交样条小波正变换和一维沃尔什-哈达玛变换；所述三维逆变换包括依次进行一维沃尔什-哈达玛变换和二维双正交样条小波逆变换；

根据参考图像块的大小，选取一定数量的相似图像块，并保持线程网格grid不变和设置线程块block的大小；

采用全局存储器合并访问的模式将选取的相似图像块数据从全局存储器送至共享存储器中以构成三维相似矩阵；

对三维矩阵进行二维双正交样条小波正变换和一维沃尔什-哈达玛变换，并在变换域中进行硬阈值滤波，硬阈值滤波后再通过一维沃尔什-哈达玛变换和二维双正交样条小波逆变换得到图像块的第一步降噪估计数据，

对硬阈值滤波后图像块中灰度值进行加权平均，将加权平均值赋给图像块的单个像素，并在加权平均时引入凯瑟窗系数进行加权优化以获取第一步图像降噪值；

同时，将三维相似矩阵进行二维双正交样条小波变换的四个滤波器系数lpd,hpd,lpr,hpr存储在常数存储器中，并利用寄存器在每个线程中定义私有变量存储中间结果。

进一步地，，所述以第一步降噪估计数据为参考，联合维纳协同滤波内核函数并行加速策略获取第二步降噪估计数据，包括：

指令混合优化加速，将三维相似矩阵的三维正变换、维纳滤波、三维逆变换以及加权估计的过程集成于维纳协同滤波内核函数中；其中所述三维正变换包括依次进行二维离散余弦正变换和一维沃尔什-哈达玛变换；所述三维逆变换包括依次进行一维沃尔什-哈达玛变换和二维离散余弦逆变换；

根据参考图像块的大小，选取一定数量的相似图像块，并保持线程网格grid不变和设置线程块block的大小；

采用全局存储器合并访问的模式将选取的相似图像块数据从全局存储器送至共享存储器中以构成三维相似矩阵；

以第一步降噪估计数据为参考，将原始三维相似矩阵进行二维离散余弦正变换和一维沃尔什-哈达玛变换；进行维纳滤波，并对维纳滤波后的三维相似矩阵进行一维沃尔什-哈达玛变换和二维离散余弦逆变换获取第二步降噪估计数据；

对维纳滤波后图像块中灰度值进行加权平均，将加权平均值赋给图像块的单个像素，并在加权平均时引入凯瑟窗系数进行加权优化以获取第二步图像降噪值；

同时，利用寄存器将二维离散余弦变换系数存储在每个线程中定义的私有变量中。

进一步地，，边界对称扩展预处理时，首先分别对左右两侧的边界列像元进行对称扩展，然后再分别对上下两侧的边界行像元进行对称扩展，并且边界扩展的像元宽度由搜索窗的半径大小确定。

本发明的有益效果：提供一种基于gpu并行加速的三维块匹配降噪方法，通过合并访问的模式将每个gpu线程块所需要的数据一次从全局存储器读取到共享存储器中，并采用循环利用共享存储器的策略，从而提高数据访问速度和减少重复访问延时，极大提升了算法整体性能，提高了计算效率；同时配合硬阈值协同滤波内核函数和维纳协同滤波内核函数能有效去除图像中的噪声，有利于对大尺寸图像实时降噪处理。

附图说明

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

图2为本发明实施例图像进行边界对称扩展预处理的示意图。

图3为本发明实施例线程网格分配的示意图。

图4为本发明实施例共享存储器使用机制的示意图。

图5为本发明实施例最小值归约并行排序的示意图。

图6为本发明使用头模的原始ct图像的示意图。

图7为本发明去噪后头模的ct图像的示意图。

图8为本发明使用体模的原始ct图像的示意图。

图9为本发明去噪后体模的ct图像的示意图。

具体实施方式

为了便于本领域技术人员的理解，下面结合实施例及附图对本发明作进一步的说明，实施方式提及的内容并非对本发明的限定。

一种基于gpu并行加速的三维块匹配降噪方法,如图1所示，包括：

为了防止ct图像边界像素越界不能进行降噪处理，将待处理的图像进行边界对称扩展预处理。如图2所示，进行对称扩展预处理时需先对噪声方差、相似块个数、搜索窗半径等参数进行设置，而边界扩展的像素宽度由搜索窗半径大小确定，一般取值为16个像素。

完成图像边界对称扩展预处理后，将扩展预处理后的图像数据从主机端发送至gpu的全局存储器中；

创建线程网络grid时需进行线程分配，如图3所示，考虑不同参考块的相似块匹配过程的互不相关性，以图像每个参考块的图像块匹配过程为一个线程块block，搜索窗内每个检测图像块与参考块的相似性匹配过程为每个线程thread。然后在图像中沿行列方向以每3个像元的步长依次递增选取参考图像块，当图像中参考图像块个数为(m*n)时，则线程网格grid的大小为(m，n)，检测参考块搜索窗内的32*32个图像块进行相似性匹配，线程块block的大小恒定设为(32,32)。

随后，采用全局存储器合并访问的模式及共享存储器多次循环利用的加速策略，对相似图像块匹配分组进行并行加速处理。如图4所示，使同一个warp中的所有线程thread都执行同一条指令访问全局存储器中连续的单元时，获得最有利的访问模式。其中一个warp束有32个线程，因此将搜索窗分成大小为32*32的4个分块，使其满足全局存储器合并访问优化的要求，尽可能的达到全局存储器带宽通信峰值。并且在每个分块中以block(16,16)的线程分配方式进行并行相似性计算获取相似块。

进行相似性计算时，以参考图像块pk为中心，半径为r的搜索窗内逐像素滑动检测窗(大小为n*n)，计算判断检测窗选中的图像块与参考图像块pk的相似性，并将相似的图像块聚集得到参考图像块pk的三维相似图像块组。相似性判断公式为：

式中，d为图像块间的距离，定义为两个图像块对应元素数值差的模除以图像块的大小，τthreod为选取的适合的距离阈值，满足d＜τthreod时，即认为两个图像块相似，反之，则不然。

完成相似性并行匹配运算后，将分块数据拷贝至共享存储器。由于全局存储器每次访问都有高达400～600个时钟周期的延迟，且图像块匹配过程中像元数据又存在着大量的重复访问，因此将搜索窗内的像元数据拷贝到高速缓存的共享存储器过程中尽可能减少对全局存储器的访问次数。而受gpu硬件条件的限制，线程块block的最大线程数为1024，将搜索窗中32*32个检测块的像元数据一次性从全局存储器中拷贝到共享存储器时，像元总个数超出了线程最大数，为了有效解决这个问题，采用将搜索窗中像元数据分块循环拷贝到共享存储器中的策略进行并行加速。同时为了避免由分块重叠造成的图像块重复匹配的情况，在分块数据拷贝至共享存储器后，设置threadidx.x<16且threaidx.y<16。

采用最小值归约并行的策略寻找规定数目的最相似图像块，如图5所示，根据实际情况选取一定数目的最相似的图像块参与后续的处理，即需在相似的图像块中选出距离最小的相似图像块。首先参考图像块的相似图像块按照相似距离由小至大排序聚集成三维矩阵，启用n个线程分别与相似性计算得到的n个距离值d[n]对应；将第i个线程的值与第(i+n/2)个线程的距离值进行比较并将较小值放在左部分中，较大值放在右部分，则左部分区间为d[0]至d[n/2]，右部分区间为d[n/2]至d[n]；此时n/2个线程并行比较后，最小值一定在左部分区间为d[0]至d[n/2]。完成n/2线程并行比较后，比较线程数减半成n/4，对左部分区间同样进行上述比较，则最小值所在的左部分区间也缩小为d[0]至d[n/4]；重复上述步骤，直至左部分区间缩小为d[0]，则此时d[0]为所有距离值中的最小值。每次归约并行排序找到最小值后，将距离数组d的起始访问位置后移一位，数组长度变为n-1，然后循环重复上述归约运算过程最终即可找到规定数目的最相似图像块，直至找到距离最小的规定数目的相似图像块。以此每次归约比较线程数都缩小为原来的一半，极大减少了排序的时间。

不管是浮点计算指令、加载指令还是分支指令均占用指令处理带宽，每个sm的指令处理带宽都是有限的。因此采用指令混合优化加速，将三维相似矩阵的三维正变换(二维双正交样条小波正变换+一维沃尔什-哈达玛变换)、硬阈值滤波、三维逆变换(一维沃尔什-哈达玛变换+二维双正交样条小波逆变换)以及加权估计的过程集成统一到一硬阈值协同滤波内核函数中，以减少不必要中间变量的反复加载、拷贝指令，节约不必要的耗时。

进行线程分配，参考图像块的大小为8*8，硬阈值协同滤波选择16个相似图像块参与处理，故将线程块block的大小具体设置为(64,16)，线程网格grid保持不变，并同样采用合并访问的方式将选取的最相似图像块数据从全局存储器拷贝到共享存储器中以构成三维相似矩阵

三维矩阵进行二维双正交样条小波(bior)正变换和块间一维沃尔什-哈达玛变换，并在变换域中进行硬阈值滤波，滤波后再通过一维沃尔什-哈达玛变换和二维双正交样条小波(bior)逆变换得到组中所有图像块的第一步降噪估计值过程为其中为组中每个图像块的二维bior变换和块间的径向一维沃尔什-哈达玛变换；为硬阈值滤波后三维矩阵的块间径向一维沃尔什-哈达玛变换和每个图像块的二维bior逆变换；γ为硬阈值滤波处理；

在三维正变换后，噪声往往集中在变换域系数较小值处，真实细节信息集中在变换域系数较大值处，因此硬阈值滤波将小于阈值参数的变换域系数置0，其他系数保留不变，过程如下式所示：

式中，x为三维矩阵三维正变换后的变换系数，为硬阈值滤波的收缩参数，σ为估计的噪声偏差。

经过硬阈值协同滤波后，每个图像块的每个像元都得到一个估计值，但对于某一个像元i，可能会出现在多个图像块内，从而拥有多个估计值，需要对这些有重叠的块的像元估计值进行加权平均从而得到像元i的基础降噪估计值，其基础权值由硬阈值滤波后三维矩阵中的非零变换域系数的个数决定，

在加权平均过程中，为了进一步减少边界影响，加入凯瑟窗系数wkaiser进行加权聚集，则基础降噪估计图像中任一像元i的值为：

式中，

同时，由于在对三维相似矩阵进行图像块的二维双正交样条小波变换时，小波变换的四个滤波器系数lpd,hpd,lpr,hpr是常量，且需要被每个线程频繁反复访问，因此将其存储在拥有缓存加速的常数存储器中，充分利用寄存器加快访问速度，节约运算耗时。寄存器和共享存储器位于gpu芯片内部，分别是访问速度最快的两个存储器。二维双正交样条小波变换过程涉及迭代运算，需要定义变量存储中间变换结果，同时充分利用寄存器，在每个线程中定义私有变量存储中间结果，以提高迭代运算效率。

以第一步降噪估计图像为参考，对原始含噪图像进行维纳滤波降噪估计。首先同样采用指令混合优化加速，将三维相似矩阵的三维正变换(二维离散余弦正变换+一维沃尔什-哈达玛变换)、维纳滤波、三维逆变换(一维沃尔什-哈达玛变换+二维离散余弦逆变换)以及加权估计的过程集成统一到一维纳协同滤波内核函数中。

进行线程分配优化时，维纳滤波选择32个相似图像块参与处理，故将线程块block的大小具体设置为(64,32)，线程网格grid保持不变，并同样合并访问拷贝到共享存储器，构成新的三维相似矩阵同时从原始含噪图像中获取原始图像的三维相似块矩阵对两个三维矩阵和均分别进行二维离散余弦正变换(dct)和一维沃尔什-哈达玛变换，然后以基础估计图像的三维正变换矩阵计算经验维纳系数

接着利用经验维纳系数对原始含噪图像的三维矩阵进行维纳滤波处理，处理完成后再通过逆变换即可得到组中所有图像块的降噪估计

利用类似第一步的原则对重叠像元估计加权平均，不同的是其中基础权值由经验维纳系数决定，则第二步降噪估计图像中任一像元i的值为：

同时，由于图像块的二维离散余弦变换中，变换系数涉及开销极大的三角函数的计算，因此将其预先计算出来，并充分利用寄存器将其存储在了每个线程定义的私有变量中。随后将第二步降噪估计数据从gpu发送至cpu主机端，并剔除扩展边界像素以获取降噪后的图像。

本发明提供的一种基于gpu并行加速的三维块匹配降噪方法，考虑到线程间的图像块检测中存在着大量的数据访问重复，且全局存储器的每次访问都有高达400～600个时钟周期的延迟，采用将每个gpu线程块所需要的数据一次从全局存储器读取到共享存储器，并重复利用共享存储器的策略，以节约重复访问延时，提升计算效率。而且数据从全局存储器读取到共享存储器时均采用合并访问的模式，能提高数据的访问速度，尽可能的达到全局存储器带宽的峰值，能将基于三维块匹配算法进行gpu加速，极大提升了算法整体性能，节省了处理时间，计算速度比传统串行算法提高了近80倍；同时配合硬阈值协同滤波内核函数和维纳协同滤波内核函数，不仅有效去除了图像中的噪声，还有利于对于大尺寸ct图像的实时降噪处理。

更具体地，图6是本发明使用的头模包含噪声原始ct图像，扫描部位是人体的口腔牙齿，图7是采用本发明方法去噪后的头模ct图像；图8是本发明是使用体模包含大量噪声的原始ct图像，图9是采用本发明方法去噪后体模的ct图像。通过对比图6、7和图8、9的实验结果，从目视效果上可以看出本发明方法有效去除了原始图像中的噪声，图像整体清晰度较好，边缘细节信息也保持良好。从处理速度上，原始串行算法处理幅面为2048*2048大小的ct图像需要1.5个小时，而本发明方法处理相同图像的耗时则仅需要69秒，速度提升了将近80倍。

以上内容仅为本发明的较佳实施例，对于本领域的普通技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，本说明书内容不应理解为对本发明的限制。