clc; close all;
X = imread('1.jpg');
X=rgb2gray(X);
Y=X;
Y = imnoise(Y, 'salt & pepper');%新增椒鹽噪聲，也可以改成其他噪聲

A=fspecial('average',3); %生成系統預定義的3X3濾波器  
Z=filter2(A,Y)/255; %均值濾波

%Z=medfilt2(Y,[3,3]);%中值濾波

%A=fspecial('gaussian'); %高斯濾波卷積核
%Z=filter2(A,Y)/255;   %用生成的高斯序列進行濾波



figure; 
subplot(1, 3, 1); imshow(X); title('原影象');
subplot(1, 3, 2); imshow(Y); title('加噪聲影象');
subplot(1, 3, 3); imshow(Z); title('濾波後圖像');

X = double(X);  
Z = double(Z);  

D = Z-X;
MSE = sum(D(:).*D(:))/numel(Z);%均方根誤差MSE
PSNR = 10*log10(255^2/MSE);%峰值信噪比
MAE=mean(mean(abs(D)));%平均絕對誤差

w = fspecial('gaussian', 11, 1.5);  %window 加窗  
K(1) = 0.01;                      
K(2) = 0.03;                      
L = 255;       
Z = double(Z);  
X = double(X);  
C1 = (K(1)*L)^2;  
C2 = (K(2)*L)^2;  
w = w/sum(sum(w));  
ua   = filter2(w, Z, 'valid');%對視窗內並沒有進行平均處理，而是與高斯卷積，  
ub   = filter2(w, X, 'valid'); % 類似加權平均  
ua_sq = ua.*ua;  
ub_sq = ub.*ub;  
ua_ub = ua.*ub;  
siga_sq = filter2(w, Z.*Z, 'valid') - ua_sq;  
sigb_sq = filter2(w, X.*X, 'valid') - ub_sq;  
sigab = filter2(w, Z.*X, 'valid') - ua_ub;  
ssim_map = ((2*ua_ub + C1).*(2*sigab + C2))./((ua_sq + ub_sq + C1).*(siga_sq + sigb_sq + C2));  
MSSIM = mean2(ssim_map);  

display(MSE);%均方根誤差MSE
display(PSNR);%峰值信噪比
display(MAE);%平均絕對誤差
display(MSSIM);%結構相似性SSIM