CN107730464A - 基于块匹配的图像降噪并行算法 - Google Patents

基于块匹配的图像降噪并行算法 Download PDF

Info

Publication number
CN107730464A
CN107730464A CN201710928330.0A CN201710928330A CN107730464A CN 107730464 A CN107730464 A CN 107730464A CN 201710928330 A CN201710928330 A CN 201710928330A CN 107730464 A CN107730464 A CN 107730464A
Authority
CN
China
Prior art keywords
segment
thread
patch
block
threads
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201710928330.0A
Other languages
English (en)
Other versions
CN107730464B (zh
Inventor
尹皓
刘东权
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Sichuan University
Original Assignee
Sichuan University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Sichuan University filed Critical Sichuan University
Priority to CN201710928330.0A priority Critical patent/CN107730464B/zh
Publication of CN107730464A publication Critical patent/CN107730464A/zh
Application granted granted Critical
Publication of CN107730464B publication Critical patent/CN107730464B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T1/00—General purpose image data processing
    • G06T1/20—Processor architectures; Processor configuration, e.g. pipelining
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20021—Dividing image into blocks, subimages or windows

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Image Analysis (AREA)

Abstract

本发明涉及图像处理领域，特别涉及一种基于块匹配的图像降噪并行算法。本发明提供的算法，根据相邻参考图块寻找相似图块时存在大量重复计算的特性，将同一行的参考图块分配为多个组别，以不同线程块分别处理，在各个线程块中，将各个参考图块的搜索邻域范围划分入多个线程组，并在各线程组中同时进行对应参考图块与相应组邻域范围内相邻图块的相似度计算；从而避免了重复计算，大幅度的加快了算法的实现时间。

Description

基于块匹配的图像降噪并行算法
技术领域
本发明涉及图像处理领域，特别涉及一种基于块匹配的图像降噪并行算法。
背景技术
现有的基于块匹配的图像降噪算法主要包含如下步骤：根据参考图块的大小确定图像中有效的参考图块的范围，根据相似图块搜索窗口的大小n，计算每一个参考图块与其邻域图块的距离，根据阈值T1和相似图块的最大个数N为每一个参考图块构建3-D相似图块栈，默认包含参考图块；对每个参考图块的所有相似图块栈做2-D FFT正变换；再对相似图块中的每一个点做16点1-D复数FFT正变换，进行硬阈值处理，变换域系数小于阈值的置为0，然后做1-D复数FFT逆变换；对参考图块的所有相似图块栈做2-D FFT逆变换；把3-D相似图块栈中的每一个图块还原到图像上对应的位置，因为图块间有重叠，所以需要对每个像素点的值进行加权累加，同时记录下图像中每个像素点的权重；最后一步对所有像素点进行加权平均得到最终的降噪图像。在采用传统方式应用基于块匹配的图像降噪算法时，需要对每个参考块到邻域图块的距离依次进行计算，这种计算方式不仅耗费时间长，而且由于各个参考图块本身是部分重叠图块，因此不同参考图块先后进行计算时，不可避免的出现对同一像素点进行重复计算的问题，造成资源浪费。
发明内容
本发明的目的在于克服现有技术应用基于块匹配的图像降噪算法时，由于依次对每个参考图块进行计算，从而造成计算时间长，且存在同一像素点重复计算的问题，提供一种并行的基于块匹配的图像降噪算法。
为了实现上述发明目的，本发明提供了以下技术方案：
一种基于块匹配的图像降噪并行算法，包含如下步骤：
在为每个参考图块确定相似图块的步骤中，确定待处理图像像素大小width×height，其中，width表示待处理图像x方向上像素，height表示待处理图像y方向上像素；确定待处理图像参考图块图块大小k×k和步长p；确定待处理图像x方向上每行参考图块数量y方向上每列的参考图块数量
以行为基准，将同一行的参考图块分配入多个线程块进行处理，各个线程块负责处理的参考图块不重复；
每个线程块包含Nwarp个线程组；每个线程组包含Nthread个线程，每个线程块负责Npatch个参考图块的处理；Nwarp、Nthread均为1以上自然数，Npatch≤Nthread；即，每个线程块负责Npatch个参考图块包含的k行，(Npatch-1)×p+k列范围的图像数据的相似度计算；
按照线程组数，以待处理的参考图块对应的搜索邻域范围划分为组邻域范围，每个线程组对应一个组邻域范围，各个组邻域范围包含的搜索图块的相对偏移不重复，该相对偏移至基于待对比参考图块的相对偏移，为了方便区分，我们将组邻域范围内的图块又称为搜索图块，每个搜索图块与参考图块大小相同，均是k×k大小、划分方式可以不同，参数p控制参考图块的步长，参数q控制搜索图块的步长；它们在计算中起到的作用不同，搜索相似图块时，参考图块作为比较基准，而此时，参考图块的被比较对象，我们称之为搜索图块，由此，每个图块都可能即当参考图块，又当搜索图块；搜索邻域范围是当前线程块负责计算的邻域范围，其由预设的搜索窗口半径和搜索步长决定，以当前计算的参考图块左上顶点坐标为原点在行，列两个方向上进行搜索；众所周知的，参考图块的邻域范围不能超出图像范围；即，当每一行中的第一个线程块的x方向的邻域范围只有其右侧搜索半径范围内，同样的，本申请中，为了方便并行计算，每个线程块均以其负责计算的第一个参考图块的相对邻域范围为其搜索邻域范围；如，每行中第一个线程块负责该行中第一个参考图块的计算，即，第一个线程块负责的所有参考图块x，y方向上的相对邻域范围均参考第一参考图块，以x方向为例，只包含自身左上顶点右侧搜索半径范围；应注意的是，此处各个组邻域范围包含的搜索图块的相对偏移不重复是指搜索范围不重复，并非搜索图块包含的图像不重复。
各个线程组负责其对应的k行，(Npatch-1)×p+k列图像数据与自身对应的组邻域范围内图像相似度的计算；理所当然的，每个线程块中的各个线程组所负责的参考图块均是相同的(具体的，均负责着Npatch个参考图块包含的k行，(Npatch-1)×p+k行图像数据的相似度计算)，只是其负责的组邻域范围各不相同；而同时，其他线程块则负责不同的参考图块的相应计算；各个线程快的计算过程为：该线程块中各个线程组均负责该线程块对应的(Npatch-1)×p+k列图像数据与本组负责的组邻域范围内各列数据的相似度计算；计算完成后，线程组中每个线程负责整合一个本线程块对应的参考图块与自身对应组邻域范围内图块的相似度计算，(即将一个参考图块包含的k列数据与组邻域范围内任一图块包含的k列数据的距离进行和计算)，将符合要求的计算结果储存在对应的组共享内存中，即各个组共享内存中，存有对应线程块负责的Npatch个参考图块在对应组邻域范围内的符合要求的相似图块。
进一步的，每个线程块中任选一个线程组，由该线程组负责自该线程块对应的所有组共享内存中读取计算结果，以该线程组包含的Nthread个线程，一对一的完成该线程块负责的Npatch个参考图块的相似图块筛选，筛选是指，比较参考图块所有的相似图块，选择相似度最高的Nmax个相似图块并顺序存储，作为该参考图块的3‐D图块栈，该Nmax为4以上的自然数，如可以是8、16、32、64。
进一步的，确定每个线程组负责计算的邻域范围的步骤中，各个线程组顺序循环分配当前搜索邻域范围。
进一步的，确定每个线程组负责计算的邻域范围的步骤中，各个线程组平均分配当前搜索邻域范围。
进一步的，计算时，每一行需要的最小线程块数此时，除最后一个线程块外，每个线程块负责Npatch＝Nthread个参考图块的计算；最后一个线程块负责num_patch_x-(num_block_x-1)*Npatch个参考图块的计算，当其小于Npatch时，需考虑数据不能超过图像的有效坐标范围。
进一步的，对于所有3-D图块栈进行批量FFT正变换，对3-D线性变换得到的变换域系数做硬阈值处理的步骤中，配置一个线程处理参考图块在变换域中一个位置对应的多个个点，判断系数值小于阈值，则置其为零，同时记录变换域系数非零项的个数用于权重计算。
进一步的，对于所有经过阈值处理的3-D图块栈进行批量FFT逆变换，将每一个3‐D图块栈中的图块恢复到其在图像中的原始位置中；配置一个线程块处理一个3-D图块栈，每个线程块中，一个线程处理图块中的一个位置对应的多个点，将相似图块栈中的所有图块还原到图像的原始位置。
进一步的，采用原子操作将每一个3‐D图块栈中的图块恢复到其在图像中的原始位置中。
进一步的，配置一个线程处理图像的一个像素点计算最终降噪图像。
与现有技术相比，本发明的有益效果：本发明提供一种基于块匹配的图像降噪并行算法，根据相邻参考图块寻找相似图块时存在大量重复计算的特性，将同一行的参考图块分配为多个组别，以不同线程块分别处理，在各个线程块中，将各个参考图块的当前搜索邻域范围划分入多个线程组，并在各线程组中同时进行对应参考图块与相应组邻域范围的相似度计算；从而避免了重复计算，大幅度的加快了算法的实现时间，经实践验证，未采用本发明提供的方法进行的基于块匹配的图像降噪算法处理图像的速度通常在数百毫秒，而经本发明方法的基于块匹配的图像降噪算法，处理同样的图像的速度通常可以降低一个数量级，即在数十毫秒，算法优化效果明显。
本发明中还使用常量内存constant保存算法参数，提高了各执行步骤中函数的访存效率；使用纹理内存texture保存含噪图像，进一步的提高了访存效率；引入参考图块patch步长p和搜索步长q的概念，使得用户对相应参数可调，缩短算法执行时间；改变kernel block/thread的配置和shared memory的大小，以适应图像大小，降低kernel调度成本和内存空间开销；同时使用3-D批量FFT变换，z方向固定16点，当相似图块的数量小于16时，z方向的点当作0处理，显著提高3-D图块栈在变换域进行硬阈值滤波的执行效率。
附图说明：
图1为基于块匹配的图像降噪算法流程图。
图2a为算法中构建3-D图块栈的流程图。
图2b为算法中在变换域进行硬阈值滤波处理流程图。
图2c为算法中图块融合流程图。
图3为本发明提供的算法构建3-D图块栈示意图。
图4为本发明提供的算法具体实施例流程图。
具体实施方式
下面结合附图及具体实施例对本发明作进一步的详细描述。但不应将此理解为本发明上述主题的范围仅限于以下的实施例，凡基于本发明内容所实现的技术均属于本发明的范围。
典型的基于块匹配的图像降噪算法包含如图1所示的如下步骤：
步骤1.如图2a所示，根据参考图块的步长p确定图像中有效的参考图块的范围，根据相似图块搜索窗口的大小n及步长q；在其n×n邻域范围内搜索与当前参考图块匹配的相似图块，图块间的相似度可以采用向量间的欧式距离来进行度量。
例如：x水平方向上的图块数量可根据如下公式计算同样的，y水平方向上的图块数量可根据如下公式计算其中，height是图像y垂直方向上的像素数width是图像x水平方向上的像素数，k是预设的参考图像x水平方向上(或y垂直方向)的像素数，p是预设的参考图块的步长，具体的，如果图块大小k＝8，步长p＝3，此时，512×512大小的二维图像在x水平方向可以划分为169个图块；在y垂直方向也可以划分为169个图块，符号表示取最邻近的最大整数以包含图像边缘的最后一个图块。
我们采用每个参考图块的左上角顶点坐标(px,py)来唯一标识该参考图块，例如，本实施例中，图块patch(0,0)，图块patch(3,0)分别表示图像x方向的第1个和第2个图块；图块patch(0,0)包含的像素点为(0:7,0:7)，图块patch(3,0)包含的像素点为(3:10,0:7)，冒号连接的两个数值分别表示图像像素点坐标的起始下标和终止下标。
如果搜索窗口的半径大小n＝19，步长q＝3，对于参考图块图块patch(0,0)，其邻域搜索范围为sx＝(0:3:18),sy＝(0:3:18)，三个数值分别表示邻域图块顶点坐标的(起始下标:步长:终止下标)；由于邻域的搜索范围不能超过图像边缘，从而对于参考图块图块patch(18,0)来说，其x方向上的邻域范围为sx＝(-18:3:18),但是y方向上，由于是在图像边缘处，因此此邻域范围是sy＝(0:3:18)。
参考图块pr与其邻域图块qs间的相似度采用向量间的欧式距离来计算；具体的，本实施例中，一个参考图块是8×8的像素图块，其被看作64×1的向量，向量间的距离可以用2-范数公式来计算。
本发明中，预定义每个参考图块的最大相似图块数量Nmax，通常，Nmax可是是8-128之间的任意数，本实施例中，我们设Nmax＝16，当图块间的距离小于阈值T1时，将当前邻域图块作为相似图块添加到该参考图块的三维相似图块栈(3-D图块栈)中；如果相似图块的数量达已达到16个，则根据比较图块栈中距离值最大的图块判断是否需要替换，同时记录相似图块的个数，具体示意图如图3。
计算每一个参考图块与其邻域图块的距离，根据阈值T1和相似图块的最大个数N为每一个参考图块构建三维相似图块栈(3-D图块栈)，每个参考图块的相似3-D图块栈中包含该参考图块本身。
步骤2.如图2b所示，对所有参考图块的相似图块栈做3-D FFT正变换；对变换域系数进行硬阈值处理，变换域系数小于阈值的置为0，然后做3-D FFT逆变换。具体的，空间域稠密且相关性高的数据在其变换域往往是稀疏的，变换域的很多系数取值为零或接近于零，因此可以在变换域进行阈值处理，对小于阈值的系数置为零以达到降噪的效果。
步骤3.如图2c所示，把3-D图块栈中的每一个图块还原到图像上对应的位置，因为图块间有重叠，所以需要对每个像素点的值进行加权累加，同时记录下图像中每个像素点的权重；最后一步对所有像素点进行加权平均得到最终的降噪图像。
具体的，得到了所有参考图块的3-D图块栈的去噪估计，根据步骤1中得到的相似图块的个数遍历3-D图块栈中的每一个相似图块，将它们融合到其在图像中对应的原始位置。在融合过程中，图块包含的像素点的值需要乘上权重系数，它们分别是
(1)w1：与3-D图块栈在变换域的非零系数的个数相关。非零系数越多，说明3-D图像栈的图块间的相似度越低，所占比重越小；具体的，w1＝1/Nnonzero，Nnonzero为非零系数。
(2)w2：与3-D线性变换的归一化系数相关；具体的，w2＝1/(k×k×N)。
(3)w3：与Kaiser窗相关，位于图块中心的像素点比重大于图块边缘的像素点，每一个像素点在图块中的几何位置是固定的，可以对应到k×k的Kaiser窗指定位置的值；具体的，w3＝wkaiser，其中，wka i ser表示Kaiser窗指定位置的值；
将所有3-D图块栈中的每一个像素点还原到图像中其对应的原始位置。由于图块间存在重叠，所以需要对每一个像素点的取值numerator及其权重值
denominator进行累加，最后再做加权平均得到最终的降噪图像。
而本实施例中，利用GPU对上述算法进行并行处理，具体的，相应的5个kernel函数，它们分别是do_patch_match,create_patch_stack,do_hard_threshold,do_aggregate,do_final；3-D batch(3-D图块栈)FFT的正变换和逆变换使用CUDA平台提供的cuFFT库，为了方便说明，本实施例中均以在x方向(行)的执行步骤为例进行说明，实际执行时，y方向上的执行步骤和x方向完全相同；具体流程如图4所示的以下步骤：
S100：do_patch_match函数运行：对每个参考图块在其邻域内搜索相似图块，并保存该相似图块的标识坐标和个数。
S200：create_patch_stack函数运行：为每个参考图块构建其3-D图块栈。
S300：对每个3-D图块栈进行FFT正变换。
S400：do_hard_rhreshold函数运行：为每个3-D图块栈的变换域系数进行硬阈值滤波。
S500：对每个3-D图块栈进行FFT逆变换。
S600：do_aggregate函数运行：将每一个3-D图块栈中的图像块恢复到其在图像中的原始位置。
S700：do_final函数运行：对图像中的每个像素点做权重归一化处理。
具体的，S100：do_patch_match函数运行步骤中，为每一个参考图块搜索其相似图块，记录搜索到的相似图块顶点的坐标和距离度量值。搜索范围(sx,sy)由搜索窗口半径n和搜索步长q决定。位于同一行的参考图块与其邻域图块计算相似度时，很多计算是重复的，例如，表一给出了参考图块分别为图块(0,0),图块(3,0),图块(6,0)时，其中部分像素需要计算的内容，以参考图块间的步长p＝3，分别计算它们与其在x方向相距3个像素的邻域图块的距离为例，表一中有下划线的1列表示当前参考图块的x坐标，第2列表示邻域图块的x坐标：
表一
从表一中可以看到，对于图块(0，0)和图块(3,0)来说，其中像素列3与像素列6的计算，像素列4与像素列7的计算，像素列5与像素列8的计算，像素列6与像素列9的计算，像素列7与像素列10的计算均为重复，同样的，对图块(3,0)和图块(6,0)来说，像素列6与像素列9的计算，像素列7与像素列10，像素列8与像素列11的计算，像素列9与像素列12的计算，像素列10与像素列13的计算；即，对于每个参考图块来说，由于与其他参考图块有重复的部分，其很多计算是重复的；我们根据图块间距离计算的冗余特性，设计并行算法来提高计算效率。
详细的，首先进行kernel配置：确定待处理图像中参考图块个数，如本实施例中，x水平方向上的图块数量可根据如下公式计算同样的，y水平方向上的图块数量可根据如下公式计算其中，height是图像y垂直方向上的像素数width是图像x水平方向上的像素数，k是预设的参考图像x水平方向上(或y垂直方向)的像素数，p是预设的参考图块的步长，具体的，如果图块大小k＝8，步长p＝3，此时，512×512大小的二维图像在x水平方向可以划分为169个图块；在y垂直方向也可以划分为169个图块，符号表示取最邻近的最大整数以包含图像边缘的最后一个图块。以CUDA平台为例具体说明详细的处理流程，在CUDA平台中，1个线程网格grid分为多个线程块block，1个线程块block包含多个线程thread。当kernel执行时，CUDA把一个个线程块block分配到流多处理器(Streaming Multiprocessor eXtreme,SMX)进行运算，而线程块block中的线程thread又会被分组为线程组warp进行调度执行。线程组warp作为kernel执行时的调度单位包含32个线程，同一个线程组warp中的线程thread执行相同的指令，处理不同的数据。
在本实施例中，1个线程块block在x方向处理Npatch个参考图块patch，本实施例中，Npatch等于线程组warp包含的线程数，实际的kernel线程块block的配置如下：
num_block_y＝num_patch_y＝169
其中，num_block_x表示x方向上需要的线程块数量，num_patch_x表示被处理图像x方向上参考图块数量，本实施例中，num_patch_x＝169，Nthread表示一个线程组中线程数量，本实施例中，Npatch＝Nthread＝32。
考虑到保存输出参数的GPU全局内存(global memory)的并发访问(coalescedaccess)，所以设置num_patch_global_x需保证其x方向保证数据分段对齐。
num_patch_global_x＝num_block_x*Nthread＝192 (5)
num_patch_global_y＝num_block_y；
1个线程块block中的线程thread个数，根据1个线程块block所需的共享内存(sharedmemory)的大小决定；如上所述，因为每一个参考图块patch均会在(sx,sy)二维空间中搜索其相似图块，所以每个参考图块均会被多次访问；为了提高访存效率，把1个线程块block处理的所有32个参考图块patch存放到共享内存(sharedmemory)中，参考图块patch所占用的空间大小为，
sm_patch_size＝k*((Nthread-1)*p+k)*sizeof(uchar) (6)
其中k表示参考图块的行数，也即参考图块patch的高度，本实施例中为8；(Nthread-1)*p+k表示一个线程组warp处理的像素列数(也即32个参考图块patch在x方向的像素个数)，本实施例中为101，sizeof(uchar)表示灰度图像的像素点的数据类型，该区域用于存放参考图块patch的数据。
1个线程块block可以包含若干个线程组warp，1个线程组warp包含32个线程thread；本实施例中，1个线程块block处理位于同一行内的32个参考图块patch，分别搜索每个参考图块邻域范围内的相似图块patch；因此可以把搜索范围按线程组warp的个数进行划分处理，分别在x,y方向搜索参考图块的相似图块；以x方向为例，同一个线程块block中的多个线程组warp分别处理这32个参考图块patch在x方向的1个或多个邻域范围；本实施例中，搜索半径n＝19，搜索步长q＝3，x方向的邻域图块patch相对下标取值为-18:3:18；如果1个线程块block包含7个线程组warp，那么每一个线程组warp处理的x方向的邻域图块patch相对下标(如上文所述，本申请中采用图块patch左上顶角坐标为其相对下标，进而以该相对下标表示该图块patch，表二中为同一行的图块，因此其y坐标均相同，所以表中仅以x坐标唯一表示对应图块)可以如表二所示，y方向的搜索方式类似x方向；
表二
1个线程组warp所需第一共享内存的空间大小sm_distance为，
sm_dis tan ce＝((Nthread-1)*p+k)*sizeof(u int) (7)
其用于保存32个参考图块patch的所有列与其当前搜索邻域列的距离，每个参考图块包含的k列的距离求和就是当前参考图块patch与其搜索邻域范围内内所有图块patch的距离。
第二共享内存的大小为sm_patch_num，该第二共享内存用于保存32个参考图块patch在当前线程组warp对应的邻域内相似图块patch的个数：
sm_patch_num＝Nthread*sizeof(u int) (8)
应注意的是，在第一共享内存，第二共享内存中均为每个线程组划分有相互独立的内存区域。
第三共享内存的大小为sm_patch_stack，该第三共享内存用于保存当前线程组warp对应的邻域内相似图块patch的顶点坐标(x,y)和距离：
sm_patch_stack＝N*Nthread*sizeof(u int 2 float 1 )(9)
这3个共享内存shared memory区域用于保存每个线程组warp计算的中间结果，最后会汇总到全局内存global memory中；根据系统提供的共享内存sharedmemory大小以及线程块block支持的线程thread最大个数计算1个线程块block可以容纳的线程组warp的个数num_warp；1个线程块block中包含的线程thread的个数为，
num_thread_x＝Nwarp*Nthread＝7*32 (10)
num_thread_y＝1；
该do_patch_match函数的输入参数包括d_noisy_tex和params；其中，d_noisy_tex表示含噪图像的纹理(texture)对象，保存在具有片上高速缓存的纹理内存中，寻址延迟小，非线性访问数据时优势非常明显；params表示算法相关的参数，其保存在只读的常量内存(constantmemory)中，这样当所有线程thread访问相同的参数时，存取效率优于global memory。
该do_patch_match函数的输出参数包括：d_patch_stack，其保存每一个参考图块patch的相似图块patch的顶点坐标，一个参考图块patch最多匹配N个相似图块patch，如N可以为8、16、32等；其字节大小为：
N*num_patch_global_x*num_patch_global_y*sizeof(uint2)＝16*192*169*8。
该do_patch_match函数的输出参数还包括d_patch_num，其保存每一个参考图块patch匹配的相似图块patch的实际个数,其占用的字节大小为：
N*num_patch_global_x*num_patch_global_y*sizeof(uint)＝192*169*4。
具体的，该do_patch_match函数的执行步骤如下：
(a)1个线程块block中的所有线程thread将待处理的32个参考图块patch所对应的像素点载入共享内存sharedmemory；
(b)遍历当前搜索邻域范围，1个线程组warp对应1个或多个的组邻域范围，分别计算该32个参考图块patch的所有列与其当前组邻域内列的距离，保存到第一共享内存sm_distance中；
(c)将属于同一参考图块patch中的所有列与属于一个邻域图块中所有列的距离合并为该参考图块patch与该邻域图块的距离；
(d)判断该距离是否小于阈值，如小于，则匹配成功，进而按距离值大小排序后添加到第三共享内存sm_patch_stack和第二共享内存sm_patch_num中；
(e)在同一线程块block中，预留一个线程组warp 0的32个线程thread负责将该线程块block负责处理的参考图块patch的相似图块添加到相应参考图块的3-D图块栈patchstack中，最大相似图块patch的个数为N，对应全局内存global memory的d_patch_stack和d_patch_num；因为1个线程块block处理32个图块patch，所以只需1个线程组warp 0，其包含的32个线程thread分别处理这32个参考图块patch；同时相似图块在存入sm_patch_stack时是按序排列的，所以只需要遍历每个线程组warp得到的中间结果，将最小的N个相似图块patch添加到global memory中。
S200步骤create_patch_stack函数运行时，该create_patch_stack函数将根据步骤100的计算结果进行，将每1个参考图块patch的相似图块patch组装为3-D图块栈patchstack，相似图块patch的最大个数为N，其中包括距离为0的参考图块patch本身。
create_patch_stack函数的kernel配置包括如下内容：CUDA线程块block在x,y方向的个数仍然由公式(5)确定，线程块block中包含的线程thread由下式给出，本步骤中，每个线程块block负责处理一个参考图块patch，本实施例中，一个图块patch的大小为8*8，本步骤中，1个线程thread处理1个像素位置对应的多个点；
num_thread_x＝k＝8 (11)
num_thread_y＝k＝8
其中，num_thread_x表示处理x方向像素点的线程数，num_thread_y表示处理y方向上像素点的线程数。
create_patch_stack函数的输入参数包括d_noisy_tex和params；其中，d_noisy_tex表示含噪图像的纹理(texture)对象，保存在具有片上高速缓存的纹理内存中，寻址延迟小，非线性访问数据时优势非常明显；由于本步骤中每个参考图块patch的相似图块patch不一定是相邻的，因此无法保证全局共享内存global memory的并发访问；params表示算法相关的参数，其保存在只读的常量内存(constantmemory)中，这样当所有线程thread访问相同的参数时，存取效率优于global memory。
相似图块栈d_transformed_stack保存每个参考图块patch的相似图块栈stack，最多N个相似图块patch字节大小为：
N*k*(k/2+1)*num_patch_global_x*num_patch_global_y*sizeof(cucomplex)
其输入数据为实数，字节大小为16*8*10*192*169*4，数据有冗余，主要是为后面的3-D FFT变换做准备，因为实信号的3-D FFT正变换的结果为复数，并且是对称的，x方向的8个变换系数cuFFT函数库只保留了其中的5个，FFT计算结果的字节大小为8*5*16*192*169*8。
S300执行步骤中，使用cuFFT函数库进行批量3-D FFT正变换。cuFFT函数库提供了一系列的函数帮助开发者进行快速FFT变换。需要进行3-D FFT正变换的图块patch stack的个数为num_patch_global_x*num_patch_global_y，3-D数据的维度(z，y，x)＝(N，k，k)，x表示最内层数据连续存放的维度，R2C变换的输入数据为自相似图块栈d_transformed_stack中取得，数据类型是实数，大小为
N*k*k*num_patch_global_x*num_patch_global_y
计算结果也保存在相似图块栈d_transformed_stack中，数据类型是复数，大小为：N*k*(k/2+1)*num_patch_global_x*num_patch_global_y。
S400步骤中，do_hard_threshold函数的执行对3-D线性变换得到的变换域系数做硬阈值处理，如果系数值小于阈值T2则置为零，同时记录变换域系数非零项的个数用于权重计算。
do_hard_threshold函数的kernel配置如下：CUDA线程块block在x,y方向的个数仍然由公式(5)确定，线程块block中包含的线程thread由下式给出，因为每个线程块block负责处理一个参考图块patch，经过3-D FFT变换一个图块patch的大小为8×5，所以1个线程thread处理图块patch中的1个像素位置数据类型为复数cuComplex的多个点。
num_thread_x＝5 (12)
num_thread_y＝k＝8。
do_hard_threshold函数的输入参数为步骤300中3-D FFT正变换的结果，保存在相似图块栈d_transformed_stack中do_hard_threshold函数的输出参数d_transformed_stack保存图块栈的处理结果字节大小为：
N*k*(k/2+1)*num_patch_global_x*num_patch_global_y*sizeof(cucomplex)；
d_stack_nonzero保存每一个参考图块patch对应的图块栈在变换域非零系数的个数的倒数，其字节大小为：
num_patch_global_x*num_patch_global_y*sizeof(float)；
do_hard_threshold函数的执行步骤如下：
(a)根据图块栈中实际相似图块patch的个数对变换域系数做硬阈值处理，将值小于阈值T2的系数赋值为零，并记录非零项的个数nz。
(b)每个线程thread中有一个寄存器register变量nz用以保存当前参考图块的3-D FFT变换系数的非零项的个数，需要将所有线程thread的nz值进行汇总，其倒数作为图块patch聚集的权重保存到global memory d_stack_nonzero中。
在这一步骤中使用CUDA提供的shfl_down函数在1个线程组warp的32个线程thread间交换数据，每个线程组warp中下标为0的线程thread保存求和结果，最后再汇总线程组warp间的求和结果得到图块栈总的系数非零项个数。对于大小为8x5的线程块block，共40个线程thread，需要2个线程组warp。
S500步骤中，执行3-D FFT逆变换，具体的，使用cuFFT函数库进行批量3-D FFT逆变换；需要进行3-D FFT逆变换的图块patch stack的个数为N*num_patch_global_x*num_patch_global_y；
其输入参数为d_transformed_stack，数据类型是复数，大小为N*k*(k/2+1)*num_patch_global_x*num_patch_global_y，计算结果也保存在d_transformed_stack中，数据类型是实数，大小为N*k*k*num_patch_global_x*num_patch_global_y。
S600步骤中，do_aggregate函数运行：将每一个3‐D图块栈中的图像块恢复到其在图像中的原始位置。其中1个线程块block处理一个图块栈，1个线程块block中的1个线程thread处理图块patch中的一个位置对应的多个点，将图块栈中的所有图块patch还原到图像的原始位置。图块patch间会有重叠，因此必须使用原子操作atomicAdd。
do_aggregate函数的kernel配置如下：CUDA线程块block在x,y方向的个数仍然由公式(5)确定，线程块block中包含的线程thread由下式给出，因为每个线程块block负责处理一个参考图块patch，一个参考图块patch的大小为8×8，所以1个线程thread处理图块patch中的1个位置对应的多个点，其处理一个参考图块x方向的线程数num_thread_x和y方向上的线程数num_thread_y分别根据下式得出：
num_thread_x＝k＝8 (13)
num_thread_y＝k＝8
do_aggregate函数的输入参数包括d_transformed_stack中保存图块栈的3-DFFT逆变换的处理结果；其字节大小为:
N*k*(k+2)*num_patch_global_x*num_patch_global_y*sizeof(float)；
其输入参数还包括d_stack_nonzero，其保存每一个参考图块patch对应的图块栈在变换域非零系数的个数的倒数；其字节大小为：
num_patch_global_x*num_patch_global_y*sizeof(float)；
其输入参数还包括S100步骤的输出参数d_patch_stack,保存了所有参考图块的相似图块的顶点坐标，以及d_patch_num保存了相似图块的个数；
do_aggregate函数的输出参数包括：d_numerator和d_denominator，其中d_numerator保存图像中每个像素点的累加和；其字节大小为width*height*sizeof(float)；d_denominator保存图像中每个像素点的权重和；其字节大小为width*height*sizeof(float)。
do_aggregate函数的执行步骤为：遍历相似图块栈中的所有相似图块patch，1个线程thread负责处理一个像素位置对应的多个点，将其还原到图像上对应的原始位置，分别将每个点的累加和以及权重和存放到d_numerator和d_denominator中。
S700步骤中，执行函数do_final，采用1个线程thread处理图像的一个像素点的方式，计算最终降噪图像。
do_final函数的kernel配置为：每个线程块block负责处理32x32个点，所以1个线程thread处理图像中的1个点。线程块block的个数由下面的公式确定：
num_thread_x＝32 (14)
num_thread_y＝32
其中，num_thread_x表示本步骤中每个线程块内负责x方向的线程数；num_thread_y表示本步骤中每个线程块内负责y方向的线程数；width、height分别表示待处理图像x方向、y方向像素数；num_block_x表示x方向上需要的线程块数，num_block_y表示y方向上需要的线程块数。
do_final函数的输入参数即为步骤S600中函数的输出函数；而do_final函数的输出参数为d_denoised，其保存待处理图像处理后的每个像素点信息。do_final函数的执行步骤，1个线程thread处理图像的一个像素点，计算d_numerator/d_denominator得到最终降噪图像。

Claims (9)

1.一种基于块匹配的图像降噪并行算法，其特征在于，包含如下步骤：
在为每个参考图块确定相似图块的步骤中，确定待处理图像像素大小width×height，其中，width表示待处理图像x方向上像素，height表示待处理图像y方向上像素；确定待处理图像参考图块大小k×k和步长p；确定待处理图像x方向上每行参考图块数量y方向上每列的参考图块数量
以行为基准，将同一行的参考图块分配入多个线程块进行处理，各个线程块负责处理的参考图块不重复；
每个线程块包含Nwarp个线程组；每个线程组包含Nthread个线程，每个线程块负责Npatch个参考图块的处理；Nwarp、Nthread均为1以上自然数，Npatch≤Nthread；即，每个线程块负责Npatch个参考图块包含的k行，(Npatch-1)×p+k列图像数据的相似度计算；
按照线程组数，将待处理的参考图块对应的搜索邻域范围划分为组邻域范围，每个线程组对应一个组邻域范围，各个组邻域范围包含的相对偏移不重复；
各个线程块的计算过程为：该线程块中各个线程组均负责该线程块对应的(Npatch-1)×p+k列图像数据与本线程组负责的组邻域范围内各列数据的相似度计算；计算完成后，线程组中每个线程负责整合一个本线程块对应的参考图块与自身对应组邻域范围内图块的相似度计算，将符合要求的图块作为当前参考图块的相似图块，并将结果储存在对应的组共享内存中。
2.如权利要求1所述的算法，其特征在于，每个线程块中任选一个线程组，由该线程组负责自该线程块对应的所有组共享内存中读取计算结果，以该线程组包含的Nthread个线程，一对一的完成该线程块负责的Npatch个参考图块的相似图块筛选，筛选是指，比较当前参考图块所有的相似图块，选择相似度最高的Nmax个相似图块构建3‐D图块栈并顺序存储。
3.如权利要求1所述的算法，其特征在于，确定每个线程组负责计算的邻域范围的步骤中，各个线程组顺序循环分配当前搜索邻域范围。
4.如权利要求1所述的算法，其特征在于，确定每个线程组负责计算的邻域范围的步骤中，各个线程组平均分配当前搜索邻域范围。
5.如权利要求1所述的算法，其特征在于，计算时，每一行需要的最小线程块数此时，除最后一个线程块外，每个线程块负责Npatch＝Nthread个参考图块的计算；最后一个线程块负责num_patch_x-(num_block_x-1)*Npatch个参考图块的计算，当其小于Npatch时，需考虑数据不能超过图像的有效坐标范围。
6.如权利要求1所述的算法，其特征在于，对于所有3-D图块栈进行批量FFT正变换，对3-D线性变换得到的变换域系数做硬阈值处理的步骤中，配置一个线程处理参考图块在变换域中对应的一个位置对应的多个点，判断系数值小于阈值，则置其为零，同时记录变换域系数非零项的个数用于权重计算。
7.如权利要求6所述的算法，其特征在于，对于所有经过阈值处理的3-D图块栈进行批量FFT逆变换，将每一个3‐D图块栈中的图块恢复到其在图像中的原始位置中；配置一个线程块处理一个3‐D图块栈，每个线程块中，一个线程处理图块中的一个位置对应的多个点，将3‐D图块栈中的所有相似图块还原到图像的原始位置。
8.如权利要求7所述的算法，其特征在于，采用原子操作将每一个3‐D图块栈中的相似图块恢复到其在图像中的原始位置中。
9.如权利要求8所述的算法，其特征在于，配置一个线程处理图像的一个像素点计算最终降噪图像。
CN201710928330.0A 2017-10-09 2017-10-09 基于块匹配的图像降噪并行算法 Active CN107730464B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710928330.0A CN107730464B (zh) 2017-10-09 2017-10-09 基于块匹配的图像降噪并行算法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710928330.0A CN107730464B (zh) 2017-10-09 2017-10-09 基于块匹配的图像降噪并行算法

Publications (2)

Publication Number Publication Date
CN107730464A true CN107730464A (zh) 2018-02-23
CN107730464B CN107730464B (zh) 2021-03-30

Family

ID=61209730

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710928330.0A Active CN107730464B (zh) 2017-10-09 2017-10-09 基于块匹配的图像降噪并行算法

Country Status (1)

Country Link
CN (1) CN107730464B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109493295A (zh) * 2018-10-31 2019-03-19 泰山学院 一种非局部哈尔变换图像去噪方法
CN110473146A (zh) * 2019-08-16 2019-11-19 苏州超擎图形软件科技发展有限公司 遥感图像显示方法、装置及存储介质和计算机设备
CN114942942A (zh) * 2022-05-18 2022-08-26 马上消费金融股份有限公司 特征数据的查询及其用户注册查询方法和装置

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101702751A (zh) * 2009-11-03 2010-05-05 上海富瀚微电子有限公司 一种视频降噪处理中的三维块匹配方法
CN101895676A (zh) * 2010-07-07 2010-11-24 上海富瀚微电子有限公司 一种适用于bm3d实时处理的集合方法
CN102572223A (zh) * 2011-12-06 2012-07-11 上海富瀚微电子有限公司 一种用于视频降噪的相似块搜索方法
US20120177119A1 (en) * 2011-01-07 2012-07-12 Sony Corporation Faster motion estimation in an avc software encoder using general purpose graphic process units (gpgpu)
CN102682429A (zh) * 2012-04-13 2012-09-19 泰山学院 一种尺寸自适应块匹配变换域滤波图像去噪方法
CN106445688A (zh) * 2016-09-30 2017-02-22 电子科技大学 一种基于mic计算平台的nlm并行图像增强方法
WO2017100061A1 (en) * 2015-12-07 2017-06-15 Google Inc. Systems and methods for multiscopic noise reduction and high-dynamic range
CN106934775A (zh) * 2017-03-08 2017-07-07 中国海洋大学 一种基于低秩恢复的非局部图像去噪方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101702751A (zh) * 2009-11-03 2010-05-05 上海富瀚微电子有限公司 一种视频降噪处理中的三维块匹配方法
CN101895676A (zh) * 2010-07-07 2010-11-24 上海富瀚微电子有限公司 一种适用于bm3d实时处理的集合方法
US20120177119A1 (en) * 2011-01-07 2012-07-12 Sony Corporation Faster motion estimation in an avc software encoder using general purpose graphic process units (gpgpu)
CN102572223A (zh) * 2011-12-06 2012-07-11 上海富瀚微电子有限公司 一种用于视频降噪的相似块搜索方法
CN102682429A (zh) * 2012-04-13 2012-09-19 泰山学院 一种尺寸自适应块匹配变换域滤波图像去噪方法
WO2017100061A1 (en) * 2015-12-07 2017-06-15 Google Inc. Systems and methods for multiscopic noise reduction and high-dynamic range
CN106445688A (zh) * 2016-09-30 2017-02-22 电子科技大学 一种基于mic计算平台的nlm并行图像增强方法
CN106934775A (zh) * 2017-03-08 2017-07-07 中国海洋大学 一种基于低秩恢复的非局部图像去噪方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
ADITHYA H. K. UPADHYA: "GPU implementation of non-local maximum likelihood estimation method for denoising magnetic resonance images", 《JOURNAL OF REAL-TIME IMAGE PROCESSING VOLUME》 *
EDUARDA MONTEIRO: "REAL-TIME BLOCK MATCHING MOTION ESTIMATION ONTO GPGPU", 《2012 19TH IEEE INTERNATIONAL CONFERENCE ON IMAGE PROCESSING》 *
于沛: "图像处理中块匹配算法的GPU并行化研究", 《中国优秀硕士学位论文全文数据库 信息科技辑》 *
刘思延: "基于自适应阈值的三维块匹配降噪算法研究", 《中国优秀硕士学位论文全文数据库 信息科技辑》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109493295A (zh) * 2018-10-31 2019-03-19 泰山学院 一种非局部哈尔变换图像去噪方法
CN109493295B (zh) * 2018-10-31 2022-02-11 泰山学院 一种非局部哈尔变换图像去噪方法
CN110473146A (zh) * 2019-08-16 2019-11-19 苏州超擎图形软件科技发展有限公司 遥感图像显示方法、装置及存储介质和计算机设备
CN110473146B (zh) * 2019-08-16 2022-12-27 苏州超擎图形软件科技发展有限公司 遥感图像显示方法、装置及存储介质和计算机设备
CN114942942A (zh) * 2022-05-18 2022-08-26 马上消费金融股份有限公司 特征数据的查询及其用户注册查询方法和装置

Also Published As

Publication number Publication date
CN107730464B (zh) 2021-03-30

Similar Documents

Publication Publication Date Title
CN109685152B (zh) 一种基于dc-spp-yolo的图像目标检测方法
KR102243036B1 (ko) 가산기를 이용한 다차원 텐서의 데이터 액세스
Chan et al. GPU-accelerated discontinuous Galerkin methods on hybrid meshes
CN107730464A (zh) 基于块匹配的图像降噪并行算法
CN110516316B (zh) 一种间断伽辽金法求解欧拉方程的gpu加速方法
US8610737B2 (en) Graphic processing unit (GPU) with configurable filtering module and operation method thereof
KR20080042083A (ko) 그래픽 처리 유닛 상에서 콘볼루션 신경망을 트레이닝하는방법
US7675524B1 (en) Image processing using enclosed block convolution
CN104198463B (zh) 拉曼光谱预处理方法及系统
CN104160396B (zh) 在字符串集之中查找最佳匹配字符串的方法和系统
CN109726441B (zh) 体和面混合gpu并行的计算电磁学dgtd方法
CN111539997B (zh) 基于gpu计算平台的图像并行配准方法、系统、装置
EP3093757A2 (en) Multi-dimensional sliding window operation for a vector processor
CN112580662A (zh) 一种基于图像特征识别鱼体方向的方法及系统
CN103390275B (zh) 动态图像拼接的方法
CN113744136A (zh) 基于通道约束多特征融合的图像超分辨率重建方法和系统
EP1760636B1 (en) Ridge direction extraction device, ridge direction extraction method, ridge direction extraction program
GB2583513A (en) Apparatus, system and method for data generation
CN106355545B (zh) 一种数字图像几何变换的处理方法及装置
CN106484532B (zh) 面向sph流体模拟的gpgpu并行计算方法
CN108304845A (zh) 图像处理方法、装置及存储介质
CN107678025A (zh) 海浪波高计算方法和装置、存储介质及处理器
CN108509532B (zh) 一种应用于地图的聚点方法和装置
Younesy et al. Visualization of time-varying volumetric data using differential time-histogram table
CN111027551B (zh) 图像处理方法、设备和介质

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant