基于CUDA架构并行设计图像去噪算法

霍迎秋 陈春荣 王雨菲 尹加

摘 要： 针对图像去噪算法中由于数据量大、计算复杂度高导致的实时性低的问题，通过对经典K?SVD图像去噪算法的并行性进行研究分析，设计基于CUDA架构的并行K?SVD图像去噪算法。该算法主要对去噪算法中矩阵拉伸、快速OMP和SVD等部分进行并行设计，采用“共享内存”“归并求和”等策略进行优化。实验结果表明，基于CUDA架构的并行算法比串行算法速度有了显著提高，最高加速比为12倍，极大提高了图像去噪算法的处理速度。

关键词： CUDA； 图像去噪； K?SVD； 图形处理器； 并行优化； 矩阵拉伸

中图分类号： TN911.73?34； TP311 文献标识码： A 文章编号： 1004?373X（2018）22?0053?06

Abstract： In allusion to the low real?time performance caused by mass data quantity and high calculation complexity in the image denoising algorithm， a parallel K?SVD image denoising algorithm based on computer unified device architecture （CUDA） is designed by means of parallelism research analysis of the typical K?SVD image denoising algorithm. In the algorithm， the parallel design of matrix stretching， fast OMP， SVD and other parts in the denoising algorithm is conducted. Strategies like memory sharing and summation by merging are adopted for optimization. The experimental results show that the parallel algorithm based on CUDA has a significant higher speed than the serial algorithm and a speed?up ratio as high as 12， which greatly improves the processing speed of the image denoising algorithm.

Keywords： CUDA； image denoising； K?SVD； graphic processing unit； parallel optimization； matrix stretching

图像处理技术被广泛应用于航空航天[1?3]、CT装置[4?5]、工业工程[6]等多个方面。其中，图像去噪算法能够去除在图像的生成、传输和存储过程中产生的噪声，提高图像的质量。目前主要的图像去噪算法有中值滤波[7]、维纳滤波[8?9]等经典方法。

本文主要研究基于K?SVD[10]（K?Singular Value Decomposition）的图像去噪算法。K?SVD是一种字典学习法[11]，基于K?SVD的图像去噪方法能够有效抑制加性高斯白噪声，在对纹理图像处理中優势明显[12]。但是这种算法计算复杂度高，计算量大，导致算法的实时性较低，难以满足工程中高实时性的实际要求。

图形处理器（Graphic Processing Unit，GPU）的通用计算能力正在高速发展中，极大地提高了计算机图形处理的速度和质量[13]。基于GPU设计并行算法成为近年来研究的热点[14?16]。本文通过利用GPU并行计算能力，提出一种基于CUDA（Computer Unified Device Architecture）架构的并行K?SVD图像去噪算法。使其能够大大降低算法所需时间，提高运行速度，因此在实际工程中具有重要的应用意义。

1 K?SVD原理

Michal Aharon等人在2006年提出了K?SVD算法[11]，它可以对字典进行自适应更新，K指算法进行了K次迭代，每次都使用了SVD算法。此算法可以针对目标图像构建特定的、能够反映目标图像特征的字典，从而在图像去噪、压缩等领域取得较好的效果[17]。

1.1 K?SVD数论基础

当K?SVD算法中每个信号只用到一个原子来近似时，K?SVD算法就变化成K?均值聚类算法。同理，稀疏表示也可以看作是广义的矢量量化（Vector Quantization，VQ），其中的多个原子的线性组合可以表示一个信号，令D表示字典，y表示训练信号，x表示训练信号的稀疏表示向量，Y是N个训练信号的集合，X是Y的解向量集合。

假设字典D是固定的，用MP（Matching Pursuit），OMP（Orthogonal Matching Pursuit）或BP（Basis Pursuit）等算法可以得到字典D，矩阵X是Y的稀疏表示系数，然后根据矩阵X可以求解出更好的训练字典D。

2 串行图像去噪算法

2.1 K?SVD算法技术实现

本文基于C语言来实现K?SVD串行图像去噪算法。在算法实现过程中主要实现图像的读取、图像的拉伸、字典的生成。通过实现快速OMP算法生成稀疏矩阵，实现SVD算法来分解矩阵然后更新字典与稀疏矩阵；通过实现OMP算法生成稀疏矩阵。然后根据生成的稀疏矩阵与字典来更新拉伸矩阵，把拉伸矩阵还原成图片实现图像去噪。

2.2 实验结果与分析

本文算法的计算机配制是Intel（R） Core（TM） E7400 2.80 GHz，内存为2 GB。在Windows 7环境下，使用Visual Studio 2010进行研究与测试，其中图片大小为256×256。首先给原始图像加上加性高斯白噪声形成噪声图像，然后使用K?SVD图像去噪算法去噪，最后得到去噪后的图像。

通过对实验结果分析发现，稀疏度对图像去噪效果和算法性能产生很重要的影响。当稀疏度为128时，实验结果如图2所示。加噪后的PSNR值为20.28，去噪之后的PSNR值为28.53。

当稀疏度为512时实验结果如图3所示。加噪后的PSNR值为20.29，去噪之后的PSNR值为28.51。

分析可知，K?SVD图像去噪算法去噪效果非常好，但是在不同稀疏度的情况下算法耗时不同，详细情况如表1所示。

对表1分析发现，串行运行K?SVD算法进行图像去噪非常耗费时间，实时性不高，对算法优化以减少其运行时间就变得非常重要。

3 基于CUDA并行设计图像去噪算法

3.1 去噪算法并行性分析

分析发现，K?SVD算法中用到大量的矩阵运算，如矩阵相乘、矩阵转置、矩阵求逆，矩阵运算可以进行并行优化。K?SVD图像去噪算法主要耗费时间是在快速OMP和SVD部分。下面就快速OMP部分和SVD部分进行并行优化设计可行性分析。

3.1.1 快速OMP并行分析

在快速OMP过程中，取训练信号的一列向量与字典相乘。从得到的向量里面选取最大的[K]个数，再从训练信号中选取对应的列数组成一个矩阵，求该矩阵的广义逆，得到相应列的稀疏矩阵。各列之间没有相互依赖关系，可以进行并行优化。

3.1.2 SVD算法并行分析

在奇异值分解过程中，主要是求矩阵的特征值和特征向量。将矩阵转置与自身相乘得到一个方阵，构造一个值全部为1的向量，方阵与这个向量相乘更新构造的那个向量，将向量的每个元素都除以该向量的最大值，迭代这个过程。研究发现，每个迭代过程之间无关联，能并行优化。

通过以上分析可知，在基于K?SVD的图像去噪算法中用到了大量矩阵运算。在矩阵运算中每个元素都可以独立地进行计算，为后续的并行优化提供了良好的基础和可能性。

3.2 基于CUDA并行设计图像去噪算法

在并行设计过程中，首先对耗时较大的部分进行并行设计，这样才能显著地降低所需要的时间。而且要结合GPU的硬件特点，合理使用硬件资源才能得到好的效果。

优化的重点放在耗时较大的算法上。整体优化过程分矩阵拉伸优化、快速OMP算法优化和[SVD]算法优化。X.m和mat.m用以表示矩阵的行，X.n和mat.n用以表示矩阵的列。

3.2.1 矩阵拉伸优化

1） 读取图像生成一个大小为256×256的图像矩阵，生成字典D。

2） 设计8×8大小的滑动窗口，从原始图像左上角滑动到右下角。将1）中得到的阵拉伸为64×62 001的矩阵。设计一个blocks为（249，249），threads为64的核函数进行矩阵拉伸变换的并行计算。其中，249是对应256-8+1，64对应8×8。从最左边滑到最右边则要滑动249次，同理，从最左边滑到最右边则也要滑动249次。生成的拉伸矩阵共有249×249=62 001列。具体情况如图4所示。

3.2.2 快速OMP算法优化

1） 使用快速OMP算法将第3.2.1中1）得到的字典与第3.2.1中2）拉伸矩阵生成稀疏矩阵。生成稀疏矩阵的函数全部用CUDA核函数实现。构建一个线程数量为D.m的核函数，D.m是字典D的行数，将字典中的一列提取出来。为避免每次循环都需要耗费80 ms时间在数据从CPU到GPU的传输上，直接在显卡端声明一个变量保存每次循环需要提取的第几列。

开辟一个X.n个线程格，每个线程格中有64个线程的核函数，每个线程计算对应元素的乘积，其结果保存在共享内存中，最后归并求和得到共享内存数组所有元素的和。

计算上述稀疏矩阵的稀疏度，令为sparse，设计一个线程数为64的核函数。提取出上述过程得到的向量中前sparse个最大的数，并记录下索引值。每个线程找出对应映射的元素里面的最大值，将512个线程中归并找出最大值，获得所有元素中最大值的索引，并将最大值赋值为0。再求次最大值，得到除去最大值之外第二大元素的记录索引。依次类推，得到向量中前sparse个最大值及其索引。

令核函数的线程块大小为X.m，线程数量是sparse，将向量中前sparse个值所对应的列取出，组成一个矩阵。其中有训练信号行数m个线程格，每个线程格中线程都对应取一个元素。该元素的索引为上述选取的最大值对应的索引，得到的元素构成一行，设计核函数的线程块为sparse，线程数为X.m，即要转置矩阵的列，对该行进行矩阵转置。矩阵相乘核函数线程格grid，即第一个矩阵的行与第二个矩阵的列构成的线程格数，线程数为X.m。将转置后的矩阵与原矩阵相乘，每个线程格负责生成相乘之后矩阵中的一个元素。每个线程计算對应元素的乘积，存放到共享内存中，最后归并求和。

对上述相乘得到的矩阵求逆。先构造一个单位矩阵，使用初等变换法求逆。通过一次消元一列的方法让矩阵变成单位阵，每次消元一列用到的核函数线程格和线程数均为sparse。映射元素如图5所示。

消元第一列之后，仅第一行第一列元素非零，第一列其余元素都为零。将求逆得到的矩阵与转置的矩阵使用线程格为grid，线程数为sparse的核函数相乘，相乘得到的矩阵再与字典相应列相乘，令其线程格为sparse，线程数为64。将向量中得到的每个元素放入到对应位置构成稀疏矩阵。将该稀疏矩阵传输到主机端。其中的一行如图6所示。

2） 将稀疏矩阵对应行的对应非零元素所在的列取出，与字典相乘。将训练信号对应列取出，这两个矩阵相减得到待SVD分解的矩阵，将该矩阵传输到显卡端。

3.2.3 SVD算法优化

1） 将第3.2.2节中2）的矩阵转置。

2） 将第3.2.2节中1）的矩阵与原矩阵相乘，核函数开辟的线程格为mat.n×mat.n个，线程数为mat.m个。使用一个线程格为1，线程数为mat.n的核函数构造一个全为1的向量。转置后的矩阵与该向量相乘核函数的线程格和线程数均为mat.n，其结果为一个向量，从中找出最大值，该向量除以得到的最大值。

上述过程循环200次，完成SVD分解。在串行算法中，保持两次最大值不改变从而完成SVD分解。如遵循串行算法的过程，则要不断地将显卡端获得的最大值传输给主机端，比较是否满足退出循环条件。传输的数据不多，但是很耗费时间。统计串行代码中的循环次数，分析得出200次足够得到一个好的特征值和特征向量。当然在稀疏度更大的情况下也要改变循环次数。用分解得到的向量更新字典与稀疏信号。

3） 将训练得到的字典与训练信号稀疏后，用稀疏信号与字典相乘从而更新训练信号，将训练信号还原成一张图片，该图片就是去噪之后的图像。

3.3 实验结果分析

算法运行环境：GPU显卡型号为GTX770，CUDA架构版本是5.5；CPU为Intel Core2 E7400 2.8 GHz，内存2 GB。

在应用OMP算法对原始信号进行稀疏分解的过程中发现得到的稀疏矩阵的稀疏度不可知，可通过经验值估计图像的稀疏度。实验测试发现，稀疏度越大，去噪效果越好，但是算法的运算量也越大。这是因为当稀疏度很大时，相应的图像处理的计算量变大，后续的奇异值分解过程也变得非常复杂。如在稀疏度为4 000时，奇异值分解耗时不小于1 h。在用快速OMP算法对原始信号进行稀疏分解的过程中可以人为地控制稀疏度。在选定稀疏度为128的情况下发现算法整体耗时减少50%，但是图像的去噪效果受到了很大影响。

选定稀疏度为128，256，512情况下，测试并行图像去噪算法的性能。分析实验结果可知，在三种稀疏度下并行去噪算法去噪后图片的PSNR值都为29 dB，通过肉眼观察结合PSNR值可知并行算法去噪效果非常好，具体如图7～图9所示。在稀疏度为128的情况下，图像去噪的效果如图7c）所示。其中加噪后PSNR值为20.29，去噪后的PSNR值是28.41。

在稀疏度为256的情况下，图像去噪效果如图8c）所示。其中加噪后的PSNR值为20.29，去噪之后的PSNR值为29.38。

在稀疏度为512时，图像去噪效果见图9c）。其中加噪后的PSNR值为20.29，去噪后的PSNR值为29.13。

随着稀疏度的逐渐增加，图像去噪效果越来越好，但是算法的运算量也越来越大。算法运算量越大越能够发挥GPU并行处理的优势，相应的加速效果越明显。例如：在稀疏度为128时，串行去噪算法整体耗时为827 508 ms，而并行算法耗时76 812 ms。在稀疏度为512时，串行去噪算法整体耗时高达5 952 046 ms，并行去噪算法耗时484 932 ms。详细如表2所示。

在不同稀疏度情况下并行去噪算法的加速比如图10所示。

通过对表2和图10分析可知，随着稀疏度不断变大，并行算法相对于串行算法的加速比不断提升。在稀疏度为128时，加速比为10.8；稀疏度为256时，加速比为11.6；稀疏度为128时，加速比为12.3。分析其主要原因是随着稀疏度不断变大，去噪算法的运算量也逐渐增大。去噪算法的非相关运算量增大，越能够充分发挥GPU大量并行处理单元的优势，因此加速比逐渐增大。综合考量并行去噪算法的去噪效果和运行速度可知，并行去噪算法在保持算法去噪效果的同时，能够极大地加速算法的运行，提高算法的实时性。

4 结 语

本文主要研究了基于K?SVD的图像去噪算法，提出一种基于CUDA架构的并行图像去噪算法。并行算法充分利用了GPU強大的并行计算能力，对基于K?SVD的图像去噪算法进行了并行优化设计。共设计了三组实验对比分析串、并行算法的性能以及去噪效果。分别测试在不同稀疏度下的串行、并行算法的耗时和去噪效果。实验结果表明，在128，256，512三种稀疏度下串、并行算法去噪效果没有区别，去噪后的图片PSNR值都为29 dB。但是并行图像去噪算法相对于串行算法运行速度得到了很大的提高，最高加速比达12，极大地提高了算法的实时处理性能。因此，本文设计的并行图像去噪算法能够替代串行去噪算法，且能够有效提高图像处理的实时性，在航空航天、CT装置、工业、工程等需要进行图像去噪的领域具有广阔的应用前景。

参考文献

[1] 赵育良，许兆林.基于图像处理技术的航空相机镜头焦面自准直定位研究[J].应用光学，2012，33（2）：288?292.

ZHAO Yuliang， XU Zhaolin. Autocollimation location system of aerial camera lens′ focus plane based on image processing technology [J]. Journal of applied optics， 2012， 33（2）： 288?292.

[2] 孙世杰，张凯，孙力，等.基于图像处理的航空仪表自动判读系统设计[J].科学技术与工程，2011，11（6）：1260?1263.

SUN Shijie， ZHANG Kai， SUN Li， et al. Automatic interpretation of aviation instrument based on image processing [J]. Science technology and engineering， 2011， 11（6）： 1260?1263.

[3] 赵志彬，刘晶红.基于图像处理的航空成像设备自动调焦设计[J].液晶与显示，2010，25（6）：863?868.

ZHAO Zhibin， LIU Jinghong. Auto?focusing method for airborne image equipment based on image processing [J]. Chinese journal of liquid crystals and displays， 2010， 25（6）： 863?868.

[4] 郝景宏，姜袁，梅世强，等.基于CT图像处理技术的混凝土损伤演化研究[J].人民长江，2010，41（17）：79?83.

HAO Jinghong， JIANG Yuan， MEI Shiqiang， et al. Study on damage development of concrete based on CT image processing [J]. Yangtze River， 2010， 41（17）： 79?83.

[5] 田威，党发宁，陈厚群.基于CT图像处理技术的混凝土细观破裂分形分析[J].应用基础与工程科学学报，2012，20（3）：424?431.

TIAN Wei， DANG Faning， CHEN Houqun. Fractal analysis on meso?fracture of concrete based on the technique of CT image processing [J]. Journal of basic science and engineering， 2012， 20（3）： 424?431.

[6] 刘奇，李昌聪，黄韫栀，等.对接管焊缝的双壁双投影工业X射线图像处理研究[J].四川大学学报（工程科学版），2015，47（2）：89?94.

LIU Qi， LI Changcong， HUANG Yunzhi， et al. Image processing research on double?walled double?imaging industrial X?ray radiographs of butt weld [J]. Journal of Sichuan University （Engineering science edition）， 2015， 47（2）： 89?94.

[7] 丁继生，卫伟，杨依忠，等.一种基于FPGA的开关中值滤波算法研究[J].合肥工业大学学报（自然科学版），2016，39（4）：490?493.

DING Jisheng， WEI Wei， YANG Yizhong， et al. Research on switching median filtering algorithm based on FPGA [J]. Journal of Hefei University of Technology， 2016， 39（4）： 490?493.

[8] 汪祖辉，孙刘杰，邵雪，等.一种结合小波变换和维纳滤波的图像去噪算法[J].包装工程，2016，37（13）：173?178.

WANG Zuhui， SUN Liujie， SHAO Xue， et al. An image denoising algorithm combined with wavelet transform and Wiener filtering [J]. Packaging engineering， 2016， 37（13）： 173?178.

[9] 张小波.基于维纳滤波的图像去噪算法研究[D].西安：西安电子科技大学，2014.

ZHANG Xiaobo. Research of image denoising algorithms based on Wiener filter [D]. Xian： Xidian University， 2014.

[10] 焦莉娟，王文剑.一种快速的K?SVD图像去噪方法[J].小型微型計算机系统，2016，37（7）：1608?1612.

JIAO Lijuan， WANG Wenjian. Speeded?up K?SVD image denoising algorithm [J]. Journal of Chinese computer systems， 2016， 37（7）： 1608?1612.

[11] AHARON M， ELAD M， BRUCKSTEIN A. K?SVD： an algorithm for designing overcomplete dictionaries for sparse representation [J]. IEEE transactions on signal processing， 2006， 54（11）： 4311?4322.

[12] ELAD M， AHARON M. Image denoising via sparse and redundant representations over learned dictionaries [J]. IEEE transactions on image processing， 2006， 15（12）： 3736?3745.

[13] 张海军，陈圣波，张旭晴，等.基于GPU的遥感图像快速去噪处理[J].城市勘测，2010（2）：96?98.

ZHANG Haijun， CHEN Shengbo， ZHANG Xuqing， et al. GPU?based denoising to remotely sensing images [J]. Urban geotechnical investigation & surveying， 2010（2）： 96?98.

[14] PARK J Y， CHUNG K S. Parallel LDPC decoding using CUDA and OpenMP [J]. EURASIP journal on wireless communications and networking， 2011（1）： 172?179.

[15] MROZEK D， BROZEK M， MALYSIAK?MROZEK B. Parallel implementation of 3D protein structure similarity searches using a GPU and the CUDA [J]. Journal of molecular modeling， 2014， 20（2）： 2067.

[16] 霍迎秋，秦仁波，邢彩燕，等.基于CUDA的并行K?means聚类图像分割算法优化[J].农业机械学报，2014，45（11）：47?53.

HUO Yingqiu， QIN Renbo， XING Caiyan， et al. CUDA?based parallel K?means clustering algorithm [J]. Transactions of the Chinese Society of Agricultural Machinery， 2014， 45（11）： 47?53.

[17] 张晓阳.基于K?SVD和残差比的稀疏表示图像去噪研究[D].重庆：重庆大学，2012.

ZHANG Xiaoyang. Image denoising via sparse and redundant representations over K?SVD algorithm and residual ratio iteration termination [D]. Chongqing： Chongqing University， 2012.