基于gpu加速的非局部平均滤波实时处理方法
【专利摘要】本发明公开了一种基于GPU加速的非局部平均滤波实时处理方法，主要解决了非局部平均滤波算法难以满足实时性要求这一问题。本方法利用GPU卓越的多线程并行数据处理能力，使用OpenCL平台，在GPU上并行地实现图像的滤波处理。在滤波过程中，对中心图像片和相似图像片进行投影，并且利用拆分核函数、优化访存等方法有效减少了计算量和数据读取时间。与现有的非局部平均滤波方法相比，本发明能够在保证去噪性能的前提下，对VGA(分辨率为640*480)和SVGA(分辨率为600*800)图像进行实时的非局部平均滤波处理，较之传统的串行处理方法，本方法取得了最高达111倍的加速比。
【专利说明】
基于GPU加速的非局部平均滤波实时处理方法
技术领域
[0001] 本发明属于图像处理技术领域，更进一步涉及一种自然图像处理滤波技术领域的 基于图形计算单元GPU(Graphics Processing Unit)加速的非局部平均滤波实时处理方 法。该方法可用于对医学、卫星遥感、工业、军事和视频多媒体等领域的数字图像进行预处 理，缩短预处理时间，达到实时处理的效果。
【背景技术】
[0002] 作为最好的图像去噪方法之一，非局部平均滤波算法能够达到很好的去噪效果。 但是由于算法复杂度高，处理时间长，所以很难在实际应用中被直接使用。自从Buades提出 非局部平均的思想以来，研究者们也提出了各种各样的改进算法。
[0003] 随着大数据时代的到来，人们急切的想寻求一种可以加快数据处理速度的方案， 此时，异构并行计算模式应运而生。GPU专为计算密集型、高度并行化的程序设计，并且在浮 点计算能力方面表现超群，自然而然的，处理机制从多核处理器过渡到GPU，并广泛应用于 各个领域。各行各业的人们都尝试利用这一新兴技术加速数据处理，在图像处理领域就有 很多算法已经在异构平台下成功实现。
[0004] 华中科技大学申请的专利"一种三维超声图像非局部均值滤波方法"（专利申请号 201310428718.6,公布号CN103544682A)中公布了一种三维超声图像非局部平均滤波方法。 该方法将待滤波的三维图像拆分成基准块，给每个基准块指定搜索区域，然后将搜索区域 中的所有与基准块大小相同的图像块作为相似块，并且计算每一相似块与基准块之间相似 度，接下来以相似块的加权平均值作为基准块滤波结果，权重为相似度，最后将基准块的滤 波结果进行整合得到三维超声图像最终滤波结果，该方法应用GPU加速，一定程度上缩短了 滤波时间。但是该方法存在的不足之处是，在利用GPU进行非局部平均滤波的过程中，算法 复杂度高导致运行时间长，GPU数据访存没有优化导致运行时间长，有重复的数据处理过程 导致运行时间长，无法满足对图像进行非局部平均滤波处理的实时性要求。
【发明内容】
[0005] 本发明的目的在于克服上述现有技术中的不足，提出一种基于GPU(GraphicS Processing Unit)加速的非局部平均滤波实时处理方法。通过对图像片进行投影，减小算 法复杂度，优化访存方式，缩短访存时间，并将非局部平均滤波在GPU的实现拆分成两个内 核，减少运算量。
[0006] 为实现上述目的，本发明具体实现步骤包括如下：
[0007] (1)搭建环境；
[0008] (la)调用开放运算语言OpenCL提供的应用程序编程接口API进行平台查询，获取 该平台上可用的设备对象；
[0009] (lb)调用开放运算语言OpenCL提供的应用程序编程接口 API依次创建上下文、创 建命令队列、构建程序对象、创建2个内核对象；
[0010] (2)读取待处理图像；
[0011] 将中央处理器CPU作为主机端读取待处理图像，将读取的待处理图像保存到中央 处理器CPU的内存；
[0012] (3)按照下式，计算半径为Z的一维高斯模板：
[0014]其中，W(r)表示一维高斯模板W上与该模板中心点相距r处的值，W表示一维高斯模 板，r表不一维高斯模板上任意点到中心点的距离，rG (-Z，Z)，G表不属于符号，Z表不一维 高斯模板的半径，e表示自然对数，〇表示高斯模板的标准差；
[0015] (4)拷贝高斯模板和待处理图像；
[0016] (4a)在图形计算单元GPU中申请用于保存高斯模板和待处理图像的显存；
[0017] (4b)将高斯模板拷贝到显存的常量存储器中，将待处理图像拷贝到显存的全局存 储器中；
[0018] (5)执行第一个内核对象；
[0019] (5a)给待处理图像的每一个像素点分配一个工作项；
[0020] (5b)给执行第一个内核对象的每一个工作组分配256个工作项；
[0021] (5c)给第一个内核对象设置4个参数，分别为:待处理图像的长度，待处理图像的 宽度，拷贝进GHJ显存的待处理图像和高斯模板；
[0022] (5d)将图形计算单元GPU全局存储器中的待处理图像缓存到局部存储器中；
[0023] (5e)从缓存到局部存储器中的待处理图像依次读取21X21个像素点，组成搜索 窗；
[0024] (5f)读取缓存到局部存储器中待处理图像的一个小图像片，该小图像片是以搜索 窗中心像素点为中心，大小为dXd的二维搜索窗中心图像片，d = 2 XZ+1，表示二维搜索窗 中心图像片的边长，Z表示高斯模板的半径；
[0025] (5g)读取缓存到局部存储器中待处理图像的一个小图像片，该小图像片是以搜索 窗下半部分每一个像素点为中心，大小为dXd的二维搜索窗相似图像片，d = 2XZ+l，表示 二维搜索窗相似图像片的边长，Z表示高斯模板的半径；
[0026] (5h)按照下式，对二维搜索窗中心图像片和二维搜索窗相似图像片进行投影，得 到一维搜索窗中心图像片和一维搜索窗相似图像片：
[0029]其中，运dxl(；a；>表示以目标像素点(a，b)为中心的二维搜索窗中心图像片IPdxd (a，b) 经过投影后得到的一维搜索窗中心图像片，d = 2XZ+l，Z表示高斯模板的半径，x表示二维 搜索窗中心图像片IPdxd(a，b)的列序号，xG [b_Z，b+Z]，G表不属于符号，E表不求和操作， 表示以像素点（i，j)为中心的二维搜索窗相似图像片ISdxd(i，j)经过投影后得到的 一维搜索窗相似图像片，y表示二维搜索窗相似图像片ISdxd(i，j)的列序号，yG[j-Z，j+Z]; [0030] (5i)按照下式，计算一维搜索窗中心图像片和一维搜索窗相似图像片高斯加权的 欧式距离：
[0031 ] Temp(i, j,a,b) - W ? ISdxi(i)-IPdxi(a) L J 2
[0032]其中，Temp (i , j , a, b)表示一维搜索窗中心图像片iPdrt (_a)和一维搜索窗相似图像 片运dxl(i)高斯加权的欧式距离，W表示高斯模板，?表示卷积操作，|||表示欧氏距离操作；
[0033] (5j)在图形计算单元GPU的全局存储器中创建一个数组，将高斯加权的欧氏距离 保存在该数组里；
[0034] (6)执行第二个内核对象；
[0035] (6a)给待处理图像的每个像素点分配一个工作项；
[0036] (6b)给执行第二个内核对象的每个工作组分配256个工作项；
[0037] (6c)在图形计算单元GPU中申请用于存放滤波后图像的显存；
[0038] (6d)给第二个内核设置4个参数:待处理图像的长度，待处理图像的宽度，拷贝进 GHJ缓存的待处理图像和滤波后图像；
[0039] (6e)将图形计算单元GPU的全局存储器里的待处理图像缓存到局部存储器；
[0040] (6f)从缓存到局部存储器中的待处理图像依次读取21X21个像素点，组成搜索 窗；
[0041] (6g)按照下式，计算搜索窗上像素点的权值：
[0043]其中，《(1，」，&，13)表示搜索窗上像素点（1，」）的权值，6表示自然对数，1 1表示衰减 因子，Temp (i，j，a，b)表不一维搜索窗中心图像片IPdxi (a)和一维搜索窗相似图像片ISdxi (i)_ 高斯加权的欧氏距离；
[0044] (6h)按照下式，计算搜索窗中心像素点经过非局部平均滤波后的像素值：
[0046] 其中，$(a，b):表示搜索窗中心像素点的像素值y(a，b)经过非局部平均滤波后的像 素值，E表示求和操作，（i，j)表示搜索窗上像素点，G表示属于符号，Q(a，b)表示以搜索 窗中心像素点(a，b)为中心的21X21大小的邻域，《(1，」， &，13)表示搜索窗上像素点（1，」） 的权值，y (i，j)表示像素点（i，j)处的像素值；
[0047] (6i)用搜索窗中心像素点滤波后的像素值赋给步骤(6c)申请的滤波后图片的对 应像素点，得到图形计算单元GHJ中的滤波后图像；
[0048] (7)保存结果图像；
[0049]从图形计算单元GPU拷贝滤波后图像，保存于中央处理器CPU内存中。
[0050]本发明与现有的技术相比具有以下优点：
[0051]第一，由于本发明对二维搜索窗中心图像片和二维搜索窗相似图像片进行投影， 得到一维搜索窗中心图像片和一维搜索窗相似图像片，克服了现有技术的算法复杂度高导 致运行时间长，无法满足对图像进行非局部平均滤波处理的实时性要求的问题，使得本发 明算法复杂度明显降低，原本需要对dXd个像素点间进行的运算减少到dXl，运算时间显 著减少，从而满足实时性要求。
[0052]第二，由于本发明将高斯模板拷贝到显存的常量存储器中，将图形计算单元GPU全 局存储器中的待处理图像缓存到局部存储器中，克服了 GPU数据访存没有优化导致运行时 间过长，无法满足对图像进行非局部平均滤波处理的实时性要求的问题，使得本发明访存 时间缩短，工作项不需要直接从全局存储器存取数据，进而缩短了运算时间，满足实时性要 求。
[0053]第三，由于本发明计算一维搜索窗中心图像片和一维搜索窗下半部分相似图像片 高斯加权的欧式距离，克服了现有技术存在的有重复的数据处理过程导致运行时间长，无 法满足对图像进行非局部平均滤波处理的实时性要求的问题，使得本发明避免了大量的重 复数据处理的过程，计算量明显减少，缩短了运算时间，满足实时性要求。
【附图说明】
[0054]图1为本发明的流程图；
[0055]图2为本发明的仿真结果图。
【具体实施方式】
[0056] 下面结合附图对本发明作进一步的详细描述。
[0057] 参照附图1，对本发明的步骤作进一步的详细描述。
[0058] 步骤1，搭建环境。
[0059] 调用开放运算语言OpenCL提供的应用程序编程接口API进行平台查询，获取平台 上可用的设备对象。
[0060] 调用开放运算语言OpenCL提供的应用程序编程接口 API依次创建上下文，创建命 令队列，构建程序对象，创建2个内核对象。
[0061 ] 步骤2,读取待处理图像。
[0062]将中央处理器CPU作为主机端读取待处理图像，将读取的待处理图像保存到中央 处理器CPU的内存。
[0063]步骤3,按照下式，计算半径为Z的一维高斯模板：
[0065]其中，W(r)表示一维高斯模板W上与该模板中心点相距r处的值，W表示一维高斯模 板，r表不一维高斯模板上任意点到中心点的距离，rG (-Z，Z)，G表不属于符号，Z表不一维 高斯模板的半径，e表示自然对数，〇表示高斯模板的标准差。
[0066] -维高斯模板的半径Z的取值范围是1~5个像素。高斯模板的标准差〇的取值范围 是[0，3]。
[0067]步骤4，拷贝高斯模板和待处理图像到显存。
[0068]在图形计算单元GPU中申请显存，用来保存高斯模板和待处理图像。
[0069] 将高斯模板拷贝到显存的常量存储器中，将待处理图像拷贝到显存的全局存储器 中。
[0070] 步骤5,执行第一个内核对象。
[0071] 给待处理图像的每以个像素点分配一个工作项。
[0072]给执行第一个内核对象的每一个工作组分配256个工作项。
[0073]给第一个内核对象设置4个参数，分别为:待处理图像的长度，待处理图像的宽度， 拷贝进GPU显存的待处理图像和高斯模板。
[0074]将图形计算单元GPU全局存储器中的待处理图像缓存到局部存储器中。
[0075]从缓存到局部存储器中的待处理图像依次读取21X21个像素点，组成搜索窗。
[0076] 读取缓存到局部存储器中待处理图像的一个小图像片，该小图像片是以搜索窗中 心像素点为中心，大小为d X d的二维搜索窗中心图像片，d = 2 X Z+1，Z表示高斯模板的半 径。
[0077] 读取缓存到局部存储器中待处理图像的一个小图像片，该小图像片是以搜索窗下 半部分每一个像素点为中心，大小为dXd的二维搜索窗下半部分相似图像片，d = 2XZ+l，Z 表示高斯模板的半径。
[0078] 按照下式，对二维搜索窗中心图像片和二维搜索窗下半部分相似图像片进行投 影，得到一维搜索窗中心图像片和一维搜索窗下半部分相似图像片：
[0081] 其中，I^ixi(a)表示以目标像素点(a，b)为中心的二维搜索窗中心图像片IPdxd(a，b) 经过投影后得到的一维搜索窗中心图像片，d = 2XZ+l，Z表示高斯模板的半径，x表示二维 搜索窗中心图像片IPdxd(a，b)的列序号，xG [b_Z，b+Z]，G表不属于符号，E表不求和操作， 虱xl(i)表示以像素点（i，j)为中心的二维搜索窗下半部分相似图像片ISdxd(i，j)经过投影 后得到的一维搜索窗下半部分相似图像片，y表示二维搜索窗下半部分相似图像片ISdxd(i， j)的列序号，y e [ j_Z, j+Z]。
[0082] 按照下式，计算一维搜索窗中心图像片和一维搜索窗下半部分相似图像片高斯加 权的欧式距离： _ _ 2
[0083] Temp(i9j，a5b) = W? ISdx丨(i)-IPdxi(a) L J 2
[0084] 其中，Temp(i，j，a，b)表示一维搜索窗中心图像片和一维搜索窗下半部分 相似图像片高斯加权的欧式距离，W表示高斯模板，0表示卷积操作，|||表示欧氏距 离操作
[0085]在图形计算单元GPU的全局存储器中创建一个数组，将高斯加权的欧氏距离保存 在该数组里。
[0086]步骤6,执行第二个内核对象。
[0087]给待处理图像的每个像素点分配一个工作项。
[0088]给执行第二个内核对象的每个工作组分配256个工作项。
[0089]在图形计算单元GPU中申请用来存放滤波后图像的显存。
[0090]给第二个内核设置4个参数:待处理图像的长度，待处理图像的宽度，拷贝进GPU缓 存的待处理图像和滤波后图像。
[0091]将图形计算单元GPU的全局存储器里的待处理图像缓存到局部存储器。
[0092]从缓存到局部存储器中的待处理图像依次读取21X21个像素点，组成搜索窗。 [0093]按照下式，计算搜索窗上像素点的权值：
[0095] 其中，《(1，」，8，13)表示搜索窗上像素点（1，」）的权值，6表示自然对数，11表示衰减 因子，1' 61^(1，」，8，13)表示一维搜索窗中心图像片1]^44和一维搜索窗下半部分相似图像 片高斯加权的欧氏距离。
[0096] 衰减因子h的取值范围是[7,20]。
[0097] 按照下式，计算搜索窗中心像素点经过非局部平均滤波后的像素值：
[0099] 其中，Ka，b)表示搜索窗中心像素点的像素值y(a，b)经过非局部平均滤波后的像 素值，E表示求和操作，（i，j)表示搜索窗上像素点，G表示属于符号，Q(a，b)表示以搜索 窗中心像素点(a，b)为中心的21X21大小的邻域，《(1，」， &，13)表示搜索窗上像素点（1，」） 的权值，y (i，j)表示像素点（i，j)处的像素值。
[0100] 用搜索窗中心像素点滤波后的像素值赋给申请的滤波后图片的对应像素点，得到 图形计算单元GHJ中的滤波后图像。
[0101]步骤7,保存结果图像。
[0102]从图形计算单元GPU拷贝滤波后图像，保存于中央处理器CPU内存中。
[0103]下面结合附图2的仿真图对本发明的效果做进一步的描述。
[0104] 1.仿真实验条件：
[0105] 本发明的仿真实验使用的CPU是Intel(R)Xeon(R)E5410,使用的GPU是AMD Radeon R9 270,软件环境是Windows7 32位，Visual Studio 2013，Matlab R2015a。仿真实验算法 去噪效果由PSNR (峰值信噪比）来表征。
[0106] 2 ?仿真效果图：
[0107]图2是使用本发明对不同分辨率图像去噪的效果图。其中，图2(a)是显示绘图阵列 VGA(Video Graphic Array)分辨率(640X480)的灰度图，用Matlab软件加方差为0? 1的高 斯噪声后，PSNR为26.8109的待处理图像，图2(b)是对图2(a)利用本发明方法的去噪后的效 果图，高斯模板半径是3,高斯模板的标准差是1.26,图2(c)是超级显示绘图阵列SVGA (Super Video Graphic Array)分辨率(600X800)的灰度图，用Matlab软件加方差为0.5的 高斯噪声后，PSNR为28.3164的待处理图像，图2(d)是对图2(c)利用本发明方法的去噪效果 图，高斯模板半径是3，高斯模板的标准差是0.86。
[0108] 表1和表2列出了本发明对附图2(a)和附图2(c)分别在5 X 5和7 X 7高斯模板尺寸 下，得到的去噪效果和加速比（时间单位为ms) XPU多线程一项指的是多核CPU采用本发明 进行非局部平均滤波所用时间，GPU多线程一项指的是GPU采用本发明进行非局部平均滤波 所用时间，加速比指的是CPU多线程所用时间与GPU多线程所用时间之比，PSNR指的是滤波 后图像的峰值信噪比，表征滤波效果。
[0109] 表1.对图2(a)处理结果
[0111] 表2.对图2(c)处理结果
[0113] 由表1和表2可见：
[0114] 经本发明提出的方法滤波后图像的峰值信噪比(PSNR)明显提高，且高斯模板大小 为7 X 7时滤波去噪效果更好。经本发明提出的加速方法较传统CPU处理方法运行时间大幅 缩短，取得的最大加速比达111，对VGA和SVGA分辨率的图像能实现实时处理。
[0115] 综上，本发明所提出的基于GPU加速的非局部平均滤波处理方法可有效滤除图像 噪声，且能通过异构并行计算大幅提升算法的处理速度，满足实时应用需求，并且采用的 OpenCL编程环境有很好的可移植性，不仅支持GPU，还支持FPGA、ARM等。
【主权项】
1. 一种基于GPU加速的非局部平均滤波实时处理方法，包括如下步骤： (1) 搭建环境； (la) 调用开放运算语言OpenCL提供的应用程序编程接口API进行平台查询，获取该平 台上可用的设备对象； (lb) 调用开放运算语言OpenCL提供的应用程序编程接口 API依次创建上下文、创建命 令队列、构建程序对象、创建2个内核对象； (2) 读取待处理图像； 将中央处理器CPU作为主机端读取待处理图像，将读取的待处理图像保存到中央处理 器CHJ的内存； (3) 按照下式，计算半径为Z的一维高斯模板：其中，W (r)表不一维高斯模板W上与该模板中心点相距r处的值，W表不一维高斯模板，r 表示一维高斯模板上任意点到中心点的距离，re (-Z，Z)，e表示属于符号，Z表示一维高斯 模板的半径，e表示自然对数，〇表示高斯模板的标准差； (4) 拷贝高斯模板和待处理图像； (4a)在图形计算单元GRJ中申请用于保存高斯模板和待处理图像的显存； (4b)将高斯模板拷贝到显存的常量存储器中，将待处理图像拷贝到显存的全局存储器 中； (5) 执行第一个内核对象； (5a)给待处理图像的每一个像素点分配一个工作项； (5b)给执行第一个内核对象的每一个工作组分配256个工作项； (5c)给第一个内核对象设置4个参数，分别为:待处理图像的长度，待处理图像的宽度， 拷贝进GPU显存的待处理图像和高斯模板； (5d)将图形计算单元GPU全局存储器中的待处理图像缓存到局部存储器中； (5e)从缓存到局部存储器中的待处理图像依次读取21X21个像素点，组成搜索窗； (5f)读取缓存到局部存储器中待处理图像的一个小图像片，该小图像片是以搜索窗中 心像素点为中心，大小为d X d的二维搜索窗中心图像片，d = 2 X Z+1，Z表示高斯模板的半 径； (5g)读取缓存到局部存储器中待处理图像的一个小图像片，该小图像片是以搜索窗下 半部分每一个像素点为中心，大小为dXd的二维搜索窗下半部分相似图像片，d = 2XZ+l，Z 表不尚斯t旲板的半径； (5h)按照下式，对二维搜索窗中心图像片和二维搜索窗下半部分相似图像片进行投 影，得到一维搜索窗中心图像片和一维搜索窗下半部分相似图像片：其中，(a)表示以目标像素点(a，b)为中心的二维搜索窗中心图像片IPdXd (a，b)经过 投影后得到的一维搜索窗中心图像片，d = 2XZ+l，Z表示高斯模板的半径，X表示二维搜索 窗中心图像片IPdxd(a，b)的列序号，xe [b_Z，b+Z]，e表不属于符号，Σ表不求和操作， 瓦ui(i)表示以像素点（i，j)为中心的二维搜索窗下半部分相似图像片ISdxd(i，j)经过投影 后得到的一维搜索窗下半部分相似图像片，y表示二维搜索窗下半部分相似图像片ISdxd(i， 」）的列序号，5^[」__2，」_+2]; (5i)按照下式，计算一维搜索窗中心图像片和一维搜索窗下半部分相似图像片高斯加 权的欧式距离：其中，Temp(i，j，a，b)表示一维搜索窗中心图像片(a)和一维搜索窗下半部分相似 图像片瓦儿)高斯加权的欧式距离，w表示高斯模板，Θ表示卷积操作，HI表示欧氏距离操 作； (5j)在图形计算单元GPU的全局存储器中创建一个数组，将高斯加权的欧氏距离保存 在该数组里； (6)执行第二个内核对象； (6a)给待处理图像的每个像素点分配一个工作项； (6b)给执行第二个内核对象的每个工作组分配256个工作项； (6c)在图形计算单元GRJ中申请用于存放滤波后图像的显存； (6d)给第二个内核设置4个参数:待处理图像的长度，待处理图像的宽度，拷贝进GPU缓 存的待处理图像和滤波后图像； (6e)将图形计算单元GPU的全局存储器里的待处理图像缓存到局部存储器； (6f)从缓存到局部存储器中的待处理图像依次读取21X21个像素点，组成搜索窗； (6g)按照下式，计算搜索窗上像素点的权值：其中，co(i，j，a，b)表示搜索窗上像素点（i，j)的权值，e表示自然对数，h表示衰减因 子，Temp(i，j，a，b)表不一维搜索窗中心图像片IP&^a)和一维搜索窗下半部分相似图像片 高斯加权的欧氏距离； (6h)按照下式，计算搜索窗中心像素点经过非局部平均滤波后的像素值：其中，f(a,b)表示搜索窗中心像素点的像素值y(a，b)经过非局部平均滤波后的像素值， Σ表示求和操作，（i，j)表示搜索窗上像素点，e表示属于符号，Q(a，b)表示以搜索窗中心 像素点(a，b)为中心的21\21大小的邻域，《(1，」，&，13)表示搜索窗上像素点（1，」）的权值， y (i，j)表示像素点（i，j)处的像素值； (6i)用搜索窗中心像素点滤波后的像素值赋给步骤(6c)申请的滤波后图片的对应像 素点，得到图形计算单元GHJ中的滤波后图像； (7)保存结果图像； 从图形计算单元GHJ拷贝滤波后图像，保存于中央处理器CRJ内存中。2. 根据权利要求1所述的基于GPU加速的非局部平均滤波实时处理方法，其特征在于： 步骤(3)中一维高斯模板半径Z的取值范围是1~5个像素。3. 根据权利要求1所述的基于GPU加速的非局部平均滤波实时处理方法，其特征在于： 步骤(3)中所述的高斯模板的标准差σ的取值范围是[〇,3]。4. 根据权利要求1所述的基于GPU加速的非局部平均滤波实时处理方法，其特征在于： 步骤(6g)中所述衰减因子h的取值范围是[7,20]。
【文档编号】G06T5/00GK105931207SQ201610370775
【公开日】2016年9月7日
【申请日】2016年5月30日
【发明人】赖睿, 王维, 孟灵非, 杨银堂, 周慧鑫, 王炳建, 秦翰林
【申请人】西安电子科技大学