非局部均值类滤波算法

%ENGLISH
%Load a mat-file containing a photograph of a colorful primate. 
% Display the indexed image using its associated colormap.

%FRANÇAIS
%Charger un fichier mat contenant un photographie d'un primate coloré
%Afficher l'image indexée en utilisant sa colormap associée

%中文
%加载一个.Mat 类型的文件，其中包含一个彩色大猩猩的图片
%通过使用对应的颜色棒来显示索引图像

%% EXEMPLE 1
clc
clear all
close all
load clown  % clown 小丑🤡
Y=X;
map1=map;

figure(1)
image(Y)
%中文
%这里注意 这时显示的图像的颜色并不正确，因为其值是索引值，而并非真正的RGB值，我们后续需要利用索引值来找到对应的RGB值，将RGB值放到图像上才可以得到正确的彩色图像

%ENGLISH
%Attention here that the colors of the image displayed at this time is not correct because the values are index values, not actual RGB
%values. We need to use these index values to find the corresponding RGB values later and place the RGB values on the image in order to
%obtain a correct color image.

%FRANÇAIS
%notez ici que les couleurs affichées à ce stade ne sont pas correct parce que la valeur est la valeur d'index au lieu de la valeur RGB
%réelles. Dans la suite, on devrai utiliser la valeur d'index pour trouver les valeur RGB correspondantes, puis placer ces valeurs RGB sur
%l'image pour obtenir une image coleur correcte

colormap(map1)
colorbar
axis off          % Remove axis ticks and numbers
axis image        % Set aspect ratio to obtain square pixels

%% EXEMPLE 2 -- method 1
clear all;clc;
load mandrill.mat
map2=map;

figure(2)
image(X)
colormap(map2)
colorbar
axis off          % Remove axis ticks and numbers
axis image        % Set aspect ratio to obtain square pixels

%中文
%直接使用gray来变成灰度图

%ENGLISH
%Directly use 'gray' to convert to grayscale image

%FRANÇAIS
%Utilisez directement 'gray' pour convertir en image en niveaux de gris

n=255; %定义灰度图的颜色级别数（这里使用了255级灰度） 也可以64级灰度
map3=gray(n);
figure(3)
image(X);
title('转化为灰度图的mandrill图像')
colormap(map3);
colorbar;
axis off          % Remove axis ticks and numbers
axis image        % Set aspect ratio to obtain square pixels

%% EXEMPLE 2 -- method 2
clear all; clc;
load mandrill.mat
map2 = map;

% Display the original color image
figure(2)
image(X)
colormap(map2)
colorbar
axis image      

%中文
%提取每个RGB值 通过加权得到灰度值  加权公式固定

%ENGLISH
%Extract each RGB value to obtain the grayscale value through weighting. The werghting formula is fixed.
%
%FRANÇAIS
%Extraire la valeur de chaque RGB et obtenir la valeur de gris en pondérant. La formule de pondération est fixée.

R=map2(:,1);  %R=reshape(R,size(X));
G=map2(:,2);  
B=map2(:,3);  

%中文
%把每列单独拿出来，得到RGB对应的值
%ENGLISH
%Take out each column separately to obtain the corresponding values of RGB
%
%FRANÇAIS
%Prenez chaque colonne séparément pour obtenir la valeur correspondantes de RGB.

R = R(X); 
G = G(X); 
B = B(X);  

%中文
%不理解 但是必须这么做，来找到实际颜色值
%ENGLISH
%Dont understand but must do so to find the actual color value
%
%FRANÇAIS
%Je comprend pas mais il faut faire comme ça pour obtenir la valeur réelle

% 使用加权公式计算灰度值
Xrgb = 0.2990 * R + 0.5870 * G + 0.1140 * B;

% 显示生成的灰度图像
figure(3)
imagesc(Xrgb);  % 使用 imagesc 显示灰度图像
colormap(gray); % 使用灰度 colormap 来正确显示灰度值
colorbar;
title('转化为灰度图的mandrill图像')
axis image      % 保持像素为正方形

%% EXEMPLE 3
clear all;clc;
load flujet
Z=X;
map4=map;
figure(4)
image(Z)
title('彩色fluid jet 图像')
colormap(map4)
colorbar
axis off         
axis image      

clc
clear
close all;
%% 读取原始图像
Input=imread('lena.bmp'); %read image
Input=double(Input); % convert the image to double type for subsequent processing
[m,n]=size(Input); % get the size of the image
%% 创建高斯噪声
d=30;
Gasnois=d*randn(m,n); %wgn(m,n)  generate gaussian noise with standard 30
STDgn=std(std(Gasnois));% calculate the standard of noise
MEANgn=mean(mean(Gasnois));% calculate the mean of noise

%这里要求高斯噪声的标准差为1.5 均值为0 但是经过上述检验可发现，其标准差和均值并不符合要求，所以在下面添加了几行代码来更加精确的生成高斯噪声
%一维数据进行一次std或mean 二维图像数据进行两次


a=0; 
b=1.5;
IMGgn=a+b.*((Gasnois-MEANgn)/STDgn);%把噪声拿来，减去他的均值，做差除以标准差，因此后面括号里就是均值是0 标准差是1的高斯噪音，本例生成均值=a，标准差=b的高斯噪声

%这个步骤在数学中就是标准化，比如说概率论与数理统计 给一个随机变量 想让它的均值是0 怎么办呢 让它减去原来的均值最后除以原来的标准差


STDYimg=std(std(IMGgn)); % same
MEANYimg=mean(mean(IMGgn));
Input_noise=Input+IMGgn; %add noise to the original image to obtain the final noisy-image 


WSH=2;  %1,2,3,7取均值滤波窗口半径
Input_r=padarray(Input_noise,[WSH,WSH],'replicate');%扩展输入图像矩阵  也就是边界扩展
[m2,n2] = size(Input_r);    %the size of the extended image 

% 显示原始图像,噪声图像和扩边图像
figure(1);
subplot(131); imshow(uint8(Input)); title('原始图像')
subplot(132); imshow(uint8(Input_noise)); title('含噪图像')
subplot(133); imshow(uint8(Input_r)); title('加取均值窗半径长作边界扩展图像');

for i=1:m
     for j=1:n   % 遍历整个图像长宽
         i1 = i+ WSH;  %原图右下角方向设置取均值窗
         j1 = j+ WSH;  %（i1，j1）为取均值窗中心    
         W1= Input_r(i1-WSH:i1+WSH , j1-WSH:j1+WSH); %W1为取均值窗v(Ni)，长宽大小2WSH+1，中心点（i1，j1）
         [mw1,nw1] = size(W1); %中心点（i1，j1）在含噪图像中是点（i,j）。
        d = sum(sum(W1));%中心点（i1,j1）矩形v(Ni)求和
        Output_filt(i,j)=d/(mw1*nw1); %中心点（i1,j1）与点（i,j）的对应 
    end 
end
figure(2);  
subplot(131); imshow(uint8(Input)); title('原始图像')
subplot(132); imshow(uint8(Input_noise)); title('含噪图像')
subplot(133); imshow(uint8(Output_filt));title('降噪后图像');

MSE=sum(sum((Input_noise-Output_filt).^2))/(m*n);% calculate mean square error  between noisy image ande filtered image
%这里是含噪图像减降噪图像然后平方求和
PSNR=10*log10((255*255)/MSE);% calculate peak signal to noise ratio , is used to measure the effect of denoising. The higher the better.
% MSE =  1159   		 PSNR =  17.4899

clc
clear
close all

image = imread('lena.bmp');

load IO.mat;
noisyimage = I0;


[denoisyimage] = non_local_means(image,noisyimage);

function [denoisyimage] = non_local_means(image,noisyimage)

%% step 1: Set the three paramaters of non local mean filtering
searchws  =  5;
similarws =  3;
h = 3;

if ~strcmp(class(image),'double')
    image = double(image)/255;
end

%% step 2: image border expansion (mirror reflection) to avoid edge effects

[H,W] = size(noisyimage);     % obtain the height and width of the image
zeroimage = zeros(H,W);  % use for storing denoised images

% Expand the boundaries of the image
expandimage = padarray(noisyimage,[searchws + similarws,searchws + similarws],'symmetric');  

figure(1)
subplot(131);  imshow(image); title('original image');
subplot(132);  imshow(uint8(noisyimage)); title('Noisy image');
subplot(133);   imshow(uint8(expandimage));  title('Expand image ')

%% step 3: Main loop
kernel = make_kernel(similarws);
kernel = kernel/sum(sum(kernel));

noDenoiseCount = 0;

for rowIdx=1:H
   for colIdx=1:W
       centerRow = rowIdx + searchws + similarws;
       centerCol = colIdx + searchws + similarws;
       w1 = expandimage(centerRow-similarws:centerRow+similarws,centerCol-similarws:centerCol+similarws); %相似窗口，长宽大小2WSH+1，中心点（i1，j1）
       [hw1,ww1] = size(w1);
       maxWeight = 0;
       weightedSum = 0; % weighted values of  various similar windows
       totalWeight = 0; % sum of weights
        
       searchRowMin = centerRow - searchws;
       searchRowMax = centerRow + searchws;
       searchColMin = centerCol - similarws;
       searchColMax = centerCol + similarws;

       L = 1;
       for r = searchRowMin : L :searchRowMax
           for c = searchColMin : L : searchColMax
               if(r == centerRow && c==centerCol) 
                   continue;
               end;

               w2 = expandimage(r-similarws:r+similarws , c-similarws:c+similarws);
               distance = sum(sum(kernel .* (w1-w2).*(w1-w2)));
               similarityWeight = exp(-distance/h.^2);
               weightedSum = weightedSum + similarityWeight * expandimage(r,c);
               totalWeight = totalWeight + similarityWeight;

               maxWeight = max(maxWeight,similarityWeight);
           end
       end

       weightedSum = weightedSum + maxWeight * expandimage(centerRow,centerCol);
       totalWeight = totalWeight + maxWeight;

       if totalWeight > 0
           denoisyimage(rowIdx,colIdx) = weightedSum / totalWeight;
       else
           denoisyimage(rowIdx,colIdx) = noisyimage(rowIdx,colIdx);

           noDenoiseCount = noDenoiseCount +1;
       end
   end
end

figure(4);
subplot(131);imshow(image);title('original image')
subplot(132); imshow(uint8(noisyimage)); title('noisy image')
subplot(133); imshow(uint8(denoisyimage));title('denoisy image') 

%% Calculate PSNR 

% can be encapsumated as a function

cumulative = sum(sum((double(noisyimage) - (denoisyimage)).^2)); % 2 sum ---> row and column
MSE = cumulative / prod(size(noisyimage));                   % divide by the total number of pixels  --->   obtain the mean square error of each pixel (MSE)
PSNR  = 10 * log10(255*255/MSE)


%% Calculate SNR
% can be encapsumated as a function

[row,col,nchannel] = size(noisyimage);

SNR = 0;
if nchannel == 1 % gray image
    signalenergy = sum( sum(( noisyimage - mean( mean( noisyimage ))) .^2));
    noiseenergy  = sum( sum(( noisyimage - denoisyimage ).^2));
    SNR = 10 * log10(signalenergy/noiseenergy)

elseif nchannel == 3 % coloful image
        for i = 1:3
            signalenergy = sum(sum((noisyimage(:,:,i)-mean(mean(noisyimage(:,:,i)))).^2));
            noiseenergy  = sum(sum((noisyimage(:,:,i)-denoisyimage(:,:,i)).^2));
            SNR = SNR + 10 * log10(signalenergy/noiseenergy);
        end
        SNR = SNR / 3
end

end

%% subfunction

function [kernel] = make_kernel(input)
kernel = zeros(2*input+1, 2*input+1);

for d = 1 : input
    value = 1 / (2*d+1)^2 ;
    for i = -d : d
        for j = -d : d
            kernel(input+1-i,input+1-j) = kernel(input+1-i,input+1-j) + value;

        end
    end
end

kernel = kernel./input;

end