去噪算法 JAVA 去噪算法属于哪类

转载

mob64ca140f67e3 2024-01-05 23:55:02

文章标签 去噪算法 JAVA denoise NLM non local means matlab 文章分类 Java 后端开发

论文原文：A non-local algorithm for image denoising

该文章2005由Buades等人发表在CVPR上，对于single-image denoise来说，当时基本上是state-of-the-art。

去噪属于图像复原的范畴，通常使用滤波来实现，并且往往是低通（平滑噪声）滤波器。对于单帧图像去噪，使用空间邻域像素来处理，对于多帧图像去噪，则可以考虑时空域相结合的方法，即时间+空间的3DNR方法。

简单的平滑滤波器有均值滤波器、高斯滤波器，算法复杂度低，但会导致图像模糊，双边滤波器是性能较好的非线性滤波器，在去噪的同时，保留了较强的纹理细节，缺点是弱的纹理被滤掉了。非局部均值（Non Local Means）方法，不仅仅考虑了像素的较小邻域，并且邻域点的权重由该点与滤波点相似度计算得到。

1. NLM方法定义

————————————————

NLM去噪后输出图像定义如下：

去噪算法 JAVA 去噪算法属于哪类_NLM

其中I为搜索区域，w(i,j)为权重，由匹配块的相似度决定。

去噪算法 JAVA 去噪算法属于哪类_non local means_02

块的相似度定义如下：

去噪算法 JAVA 去噪算法属于哪类_non local means_03

该值由表示点i和j邻域差值平方卷积高斯核，表征邻域相似度，Z(i)表示权重归一化系数。

2. 算法质量评估

使用该方法时，往往会选择一个搜索窗口和匹配窗口，典型值分别为21x21和7x7。

对于图像去噪、复原类问题，经常使用PSNR（峰值信噪比）来评价质量，单位为dB，峰值信号定义如下，L为图像最高灰度值，对于不同位宽的图像，峰值信号大小不一样。

去噪算法 JAVA 去噪算法属于哪类_去噪算法 JAVA_04

MSE（mean squared error）为均方误差，定义如下：

去噪算法 JAVA 去噪算法属于哪类_denoise_05

3. Matlab实现

Matlab代码实现如下：

(1) 顶层代码

clc;
clear all;
close all;
tic;

imgSrc = imread('E:\00_lemonHe\01_code\matlab\04_digitalImageProcessing\learn\lena.tif');
sigma = 15;
imgRandNoise = imgSrc + uint8(sigma * randn(size(imgSrc)));     %add rand noise
% h = 10 * sqrt(sigma);
h = sigma;
imgDst = NLmeansfilter(double(imgRandNoise), 10, 3, h);
imgDst = uint8(imgDst);
figure;
subplot(2,2,1), imshow(imgSrc), title('original');
subplot(2,2,2), imshow(imgRandNoise), title('noisyImage');
subplot(2,2,3), imshow(imgDst), title('denoising');
subplot(2,2,4), imshow(imgRandNoise - imgDst), title('noisy');

filterGaussian = fspecial('gaussian');      %3x3 gaussian filter
filterAverage = fspecial('average');        %3x3 average filter
imgGaussian = imfilter(imgRandNoise, filterGaussian, 'replicate');
imgAverage = imfilter(imgRandNoise, filterAverage, 'replicate');

figure;
subplot(2,2,1),imshow(imgSrc);
subplot(2,2,2),imshow(imgDst),title(['PSNR = ', num2str(calcPSNR(imgSrc, imgDst, 8))]);
subplot(2,2,3),imshow(imgGaussian),title(['PSNR = ', num2str(calcPSNR(imgSrc, imgGaussian, 8))]);
subplot(2,2,4),imshow(imgAverage),title(['PSNR = ', num2str(calcPSNR(imgSrc, imgAverage, 8))]);

toc;

(2). NLM处理函数

function [output] = NLmeansfilter(input,t,f,h)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %  https://cn.mathworks.com/matlabcentral/fileexchange/13176-non-local-means-filter
 %  input: image to be filtered
 %  t: radius of search window
 %  f: radius of similarity window
 %  h: degree of filtering
 %
 %  Author: Jose Vicente Manjon Herrera & Antoni Buades
 %  Date: 09-03-2006
 %  Modify: lemonHe
 %  Modify Date: 09-15-2017
 %
 %  Implementation of the Non local filter proposed for A. Buades, B. Coll and J.M. Morel in
 %  "A non-local algorithm for image denoising"
 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Size of the image
[m n]=size(input);
% Memory for the output
output=zeros(m,n);

% Replicate the boundaries of the input image
input2 = padarray(input,[f f],'replicate');

% Used kernel
% kernel = make_kernel(f);
% kernel = kernel / sum(sum(kernel));
kernel = fspecial('gaussian', [2*f+1 2*f+1], 1);    %gaussian filter kernel
hwait=waitbar(0,'计算中...');                       %由于计算太慢，此处加入进度bar
for i=1:m
    for j=1:n
        value = 100 * i / m;
        waitbar(i/m, hwait, sprintf('计算中:%3.2f%%',value));
        i1 = i + f;
        j1 = j + f;
        W1= input2(i1-f:i1+f , j1-f:j1+f);      %输入图像的2f+1领域
        wmax = 0;
        average = 0;
        sweight = 0;
        rmin = max(i1-t,f+1);
        rmax = min(i1+t,m+f);
        smin = max(j1-t,f+1);
        smax = min(j1+t,n+f);
        %遍历21x21的search领域
        for r = rmin:1:rmax
            for s = smin:1:smax
                if(r==i1 && s==j1)
                    continue;
                end
                W2= input2(r-f:r+f , s-f:s+f);          %search window中的2f+1领域
                d = sum(sum(kernel.*((W1-W2).^2)));     %两个2f+1邻域的高斯加权欧式距离
                w = exp(-d/(h^2));                      %像素点(r-f,s-f)的权重系数Z(i)
                if w>wmax
                    wmax = w;
                end
                sweight = sweight + w;                  %除了像素点(i,j)外的所有点的权值和
                average = average + w * input2(r,s);    %除了像素点(i,j)外的所有点的加权值
            end
        end
        average = average + wmax*input2(i1,j1);         %search window的加权和
        sweight = sweight + wmax;                       %search window的权值和   
        if sweight > 0
            output(i,j) = uint8(average / sweight);
        else
            output(i,j) = input(i,j);
        end
    end
end
close(hwait);

(3). 高斯滤波窗口函数

function [kernel] = make_kernel(f)

kernel=zeros(2*f+1,2*f+1);
for d=1:f
    value= 1 / (2*d+1)^2;
    for i=-d:d
        for j=-d:d
            kernel(f+1-i,f+1-j)= kernel(f+1-i,f+1-j) + value;
        end
    end
end
kernel = kernel ./ f;

(4) PSNR计算函数

function [imgPSNR] = calcPSNR(imgSrc, imgDst, bitWidth)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% imgSrc: original image
% imgDst: dst image
% 
% Author: lemonHe
% Date: 20170920
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(size(imgSrc,3) == 1)
    imgMSE = sum(sum((double(imgDst) - double(imgSrc)).^2)) / (size(imgSrc,1) * size(imgSrc,2));
    imgPSNR = 10 * log10((2^bitWidth-1)^2 / imgMSE);
else
    imgSrcGray = rgb2ycbcr(imgSrc);
    imgDstGray = rgb2ycbcr(imgDst);
    imgMSE = sum(sum((double(imgDstGray(:,:,1)) - double(imgSrcGray(:,:,1))).^2)) / (size(imgSrc,1) * size(imgSrc,2));
    imgPSNR = 10 * log10((2^bitWidth-1)^2 / imgMSE);
end