块匹配算法GPU并行化

来源：互联网发布：数据结构与算法的关系编辑：程序博客网时间：2024/06/05 18:39

1. 介绍

在《块匹配算法》中介绍了块匹配算法及优化策略，但这些方法都是在 CPU 端执行，由于块匹配算法的计算量很大，因此会耗费很多时间。本文侧重于块匹配算法的并行实现。

本文实验的平台为：

CPU：i7 4790k
GPU：Nvidia GTX980
Matlab 2014a
CUDA 6.5

2. 提取图像块

本文的实验以 512x512x3 的彩色 lena 图为例，每个图像块的大小为 5x5x3，块与块之间的间隔为 1 个像素。为了便于计算，将一个图像块拉成一个向量，所有的图像块就变成一个矩阵。

[height, width, ch] = size(im); % 图像的行、列和页数blockSize = win * win * ch; % 块的元素个数N         =  height - win + 1; % 块的起始行坐标范围M         =  width - win + 1; % 块的起始列坐标范围blockNum  =  N * M; % 块的个数block     =  zeros(blockSize, blockNum, 'single');% 依次提取块k = 0;for m = 1:ch    for i  = 1:win        for j  = 1:win            k    =  k+1;            blk  =  im(i : height-win+i, j : width-win+j, m);            blk = blk'; % 更改取块的方向，变成行优先取            blk  =  blk(:);            block(k,:) =  blk';        end    endend

上述提取块的方式用了优化技巧：即每次所有块中的一个元素，这样循环的次数只有 5x5x3，而如果每次提取一个块，如下所示，则循环次数为 508x508，因此效率有明显差异。

[height, width, ch] = size(im); % 图像的行、列和页数blockSize = win * win * ch; % 块的元素个数N         =  height - win + 1; % 块的起始行坐标范围M         =  width - win + 1; % 块的起始列坐标范围blockNum  =  N * M; % 块的个数block     =  zeros(blockSize, blockNum, 'single');k = 0;for i = 1:N    for j = 1:M     k = k+1;     blk = im(i : i+win-1, j : j+win-1, :);     blk = permute(blk, [2,1,3]);     blk = blk(:);     block(:, k) = blk';    endend

优化前后的时间对比为：

优化前：1.2682s优化后：0.1680s

加速比为：7.5 倍

虽然上述CPU端的程序使用优化的方式能获得更好的效果，但是在GPU段并不一定如此，为了获得更好的并行性，在本文中使用优化前的方案，即每个线程提前一个图像块，因此可以开辟 508 * 508 个线程，另外，涉及到数据的反复读取，可以使用共享内存。

并行的运行时间为：

GPU并行运行时间为：4.6756ms

相较于未优化的本加速比为 270，对于优化后的加速比为 36。

3. 方案1-local

如下图所示，每个查询块搜索局部的窗，块之间有重叠，其中查询块的间隔为两个像素。需要说明的是：搜索窗中的块之间是有重叠的，块与块之间的间隔为1个像素。

这里写图片描述

为了便于计算，可以提前计算出每个查询块的搜索范围，由四个变量决定：
+ 搜索窗的行最小坐标
+ 搜索窗的行最大坐标
+ 搜索窗的列最小坐标
+ 搜索窗的列最大坐标

matlab 上的核心代码如下所示：

for  i  =  1 : N1    for  j  =  1 : M1                   offInAllBlk   = (leftUpRow(i) - 1) * M + leftUpCol(j); % 当前中心块在所有块中的索引(行优先)        offInCentralBlk  = i * M1 + j; % 当前中心块在所有中心块中的索引(行优先)        idx   = I(rmin(i):rmax(i), cmin(j):cmax(j));  %在由[rmin:rmax, cmin:cmax]确定的搜索窗中搜索相似块        idx = idx';        idx   = idx(:);        dis   = sum(bsxfun(@minus, X(idx, :), X(offInAllBlk, :)).^2, 2); % 搜索窗内的块与中心块的距离        dis   = dis ./ (blockSize*ch); % 归一化        similarBlkInd = maxN(dis, similarBlkNum); % 找到距离最小的几个块的索引        posIdx(:,offInCentralBlk)  =  idx(similarBlkInd); % 保存相似块索引        wei   =  exp( -dis(similarBlkInd) ./ hp ); % 高斯        weiIdx(:,offInCentralBlk)  =  wei ./ (sum(wei(:))+eps); % 归一化    endend

在GPU端开辟的线程数为508*508，也就是查询块的个数。因此每个线程计算量会很大，在此对不同的窗进行测试对比：

窗的半径 4 8 10 16 20 25 匹配块数 81 289 441 1225 1681 2601 cpu时间 3.6635s 6.1905s 7.5098s 13.7732s 19.5151s 27.8922s gpu时间 222ms 714ms 1085ms 2569ms 3869ms 5838ms 加速比 16.5 8.7 6.9 5.4 5 4.8

块匹配的这个过程与 Kmeans 中计算样本与聚类中心的欧式距离很相似，计算量也很相近：

Kmeans：508x508个样本，80个聚类中心，样本长度 75，迭代 14 次，计算时间为：183.562ms（优化前），25.213ms（使用共享内存优化后）
本文：508x508个查询块，81个近似块，每个块的大小75，执行1次，计算时间为：222ms

但为何在《Kmeans 的 CUDA 并行实现》中可以实现 345 倍的加速比而在此处只有几倍。两者相比，相差几十倍。最主要的原因是：Kmeans 中计算的是 508x508个样本与固定的80个聚类中心的欧式距离；而块匹配中每个查询块匹配的近邻块不完全相同，共享内存的使用不方便，内存的读取也会不连续。下一步会对其进行进一步的优化。

考虑到上述加速效果不明显，为此对其进行优化，按照《块匹配算法》中的优化方式，首先提取计算每个块的均值和方差，然后按照下式的方式进行计算权重，以此可以避免大量的重复计算。

ω (B i k, B j) = ⎧ ⎩ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 1 Z i k e - ∥ ∥ ∥ u ( B ) - u ( B j ) ∥ ∥ ∥ 2 2 2 β δ 2 | N i |; i f μ 1 < u ( B i k ) ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ u ( B j ) ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ < 1 μ 1 a n d σ 21 < var ( u ( B i k ) ) var ( u ( B j ) ) 0 o t h e r w i s e < 1 σ 2 1

上述方案的意思是说：如果两个块的均值和方差的比值相差较大，那么就肯定不是相似块。

如果设置参数为：

μ 1 = 0.95, μ 2 = 1.05, σ 1 = 0.5, σ 2 = 1.5

优化后
的重构时间为：

窗的半径 4 8 10 16 20 25 匹配块数 81 289 441 1225 1681 2601 cpu时间 3.6635s 6.1905s 7.5098s 13.7732s 19.5151s 27.8922s gpu时间 169ms 481ms 689ms 1408ms 1978ms 2766ms 加速比 21.6 12.8 10.9 9.8 9.9 10.1

相较于优化前加速2倍左右，由于窗较小时较少的计算量不多，所以优化的加速效果小；而窗较大时，较少的计算量多，优化的加速效果明显。

如果设置参数为：

μ 1 = 0.96, μ 2 = 1.04, σ 1 = 0.95, σ 2 = 1.05

优化后的重构时间为：

窗的半径 4 8 10 16 20 25 匹配块数 81 289 441 1225 1681 2601 cpu时间 3.6635s 6.1905s 7.5098s 13.7732s 19.5151s 27.8922s gpu时间 125ms 298ms 386ms 730ms 996ms 1358ms 加速比 29.4 20.7 19.5 18.9 19.5 20.5