Non Local Means-块匹配MATLAB和GPU实现

matlab 代码

IO接口

[posIdx, weiIdx] = blockMaching(im, par)

输入：

* 图像im

* 匹配参数，window和block尺寸，step等

输出：

* 不同中心块对应的前10个类的索引posIdx

* 前10个类的权重

步骤：

设置块匹配的window大小为5*5，window移动步长为1，block大小为25*25，block的移动步长为2。

提取所有window块

对于3维图像，[d1,d2,d3] = size(im);

方法1.

for d3
for d1
for d2
%抽取5*5的window
end
end
end

从以上代码可知，循环次数为 d1*d2*d3,并且在取window时，存在内存不连续的取数过程。

方法2.

for d3
for window_row
for window_col
%%按列取一维块中不同位置的数据
end
end
end

以上代码可知，循环次数为d3*window_row*window_col,每次取数据都是连续地址。

通过试验对比，在相同的图像尺寸下，方法2比方法1要快至少10倍

//

//

中心块是以步长为2提取的，需要保存中心块索引

两个索引：

当前中心块在所有块内的索引

offInAllBlk=(leftUpCol(j)−1)∗N+leftUpRow(i) offInAllBlk = (leftUpCol(j)-1)*N + leftUpRow(i)

当前中心块在所有中心块中的索引

offInCentralBlk=(j−1)∗N1+i offInCentralBlk = (j-1)*N1 + i

找出当前中心块，在所有块中的搜索范围

计算当前搜素块内所有块和当前中心块的距离，并挑选出前10个距离最大的块、

保存索引块索引

根据距离计算出不同块的权重，利用高斯函数计算并归一化

MATLAB代码

function  [posIdx, weiIdx]   =  blockMaching(im, par)
% 功能 ： 块匹配，计算与中心块相似的块以及权重
% im ： 输入图像
% par ： 参数
% posIdx ： 与中心块相似的块的位置索引
% weiIdx ： 与中心块相似的块的权重索引

%% ------1-------%%
[height, width, ch] = size(im);
searchRadius  =  par.s1; % 搜索窗半径
similarBlkNum =  par.nblk; % 相似块的个数
win           =  par.win; % 块的大小
blockSize     =  win^2; % 每个块的元素个数
step          =  par.step; % 中心窗的步长
hp            =  80/(255^2); % 高斯参数

%% 参数
N          =  height - win + 1; % 块的起始行坐标范围
M          =  width - win + 1; % 块的起始列坐标范围
blockNum   =  N * M; % 块的个数
leftUpRow  =  1:step:N; % 每个块的行起始坐标
leftUpRow  =  [leftUpRow leftUpRow(end)+1:N]; % 添加最后一个块，因为可能最后不一定正好取完
leftUpCol  =  1:step:M; % 每个块的列起始坐标
leftUpCol  =  [leftUpCol leftUpCol(end)+1:M]; % 添加最后一个块，因为可能最后不一定正好取完

%% ------2-------%%
%% 分块
X = divideBlock3D(im, win); % 提取了所有的块，也就是说步长为1个像素
X = X';

%% 查找相似块并计算权重
I      = reshape(1:blockNum, N, M); % 所有块的索引
N1     = length(leftUpRow); % 多少行中心块
M1     = length(leftUpCol); % 多少列中心块中心块
posIdx = zeros(similarBlkNum, N1*M1 ); % 每个中心块的相似块索引
weiIdx = zeros(similarBlkNum, N1*M1 ); % 每个中心块对应的相似块的权重
rmin   = bsxfun(@max, bsxfun(@minus, leftUpRow, searchRadius), 1); % 搜索窗的行最小坐标
rmax   = bsxfun(@min, bsxfun(@plus, leftUpRow, searchRadius), N); % 搜索窗的行最大坐标
cmin   = bsxfun(@max, bsxfun(@minus, leftUpCol, searchRadius), 1); % 搜索窗的列最小坐标
cmax   = bsxfun(@min, bsxfun(@plus, leftUpCol, searchRadius), M); % 搜索窗的列最大坐标

for  i  =  1 : N1
for  j  =  1 : M1

%% ------3-------%%
offInAllBlk   = (leftUpCol(j)-1)*N + leftUpRow(i); % 当前中心块在所有块中的索引
offInCentralBlk  = (j-1)*N1 + i; % 当前中心块在所有中心块中的索引

%% ------4-------%%
idx   = I(rmin(i):rmax(i), cmin(j):cmax(j));  %在由[rmin:rmax, cmin:cmax]确定的搜索窗中搜索相似块
idx   = idx(:);

%% ------5-------%%
dis   = sum(bsxfun(@minus, X(idx, :), X(offInAllBlk, :)).^2, 2); % 搜索窗内的块与中心块的距离
dis   = dis ./ (blockSize*ch); % 归一化

%% ------6-------%%
similarBlkInd = maxN(dis, similarBlkNum); % 找到距离最小的几个块的索引
posIdx(:,offInCentralBlk)  =  idx(similarBlkInd); % 保存相似块索引

%% ------7-------%%
wei   =  exp( -dis(similarBlkInd) ./ hp ); % 高斯
weiIdx(:,offInCentralBlk)  =  wei ./ (sum(wei(:))+eps); % 归一化
end
end

function index = maxN(data, N)
% 功能 ：找到 data 中最大的 N 个数，并返回索引

[~, index] = sort(data);
index = index(1:N);

应用

块匹配在图像处理中是经常用到的，比如非局部均值滤波，BM3D滤波等操作，使用块匹配中的划窗操作可以得到还较好的滤波效果；

通过找出不同搜索块内的主要块，可以对图像进行降维处理，也可以对图像进行分类

GPU并行流程实现

其实现流程如下：

//step1
IndexTableForStep2();//产生所有块索引
IndexTableForStep1();//产生中心块索引
//step2
divideBlock3D();//抽取所有三维块
//step3
calSearchWindowCoord()//得到搜索窗行和列的最小最大坐标
//step4
//计算中心块的索引和权重
for row
{
for col
{
offInAllBlk =();//当前中心块在所有块中的索引
calSimlarBlkIdx();//计算搜索块与中心块的欧式距离
mergeSort();对得到的所有欧式距离进行归并排序
}
}
//step5
normalizationStep1();//高斯形式
normalizationStep2();//归一化

在此输入正文