%% ------------------------------------
% @Software: Visual Studio Code
% @version:
% @Author: Shengji He
% @Email: hsjbit@163.com
% @Date: 2020-04-01 16:17:34
% @LastEditors: Shengji He
% @LastEditTime: 2020-04-01 16:44:41
% @Description:
%% ------------------------------------
function metrics = metrics

metrics.calc_snr = @calc_snr;
metrics.calc_psnr = @calc_psnr;
metrics.calc_rmse = @calc_rmse;
metrics.calc_mae = @calc_mae;
metrics.compute_map_ssim = @compute_map_ssim;

end

function gauss_kernel = init_gaussian_kernel(kernel_size, sigma)
gauss_kernel = zeros(kernel_size, kernel_size);
c = 0.5 * (kernel_size - 1);
for i = 1:kernel_size
    for j = 1:kernel_size
        gauss_kernel(i, j) = (1 / (2 * pi * sigma ^ 2)) * exp(-(((i - c) ^ 2) + ((j - c) ^ 2)) / (2 * sigma ^ 2));
    end
end
totalSum = sum(gauss_kernel(:));
gauss_kernel = 1.0 / totalSum * gauss_kernel;
end


%% calc_snr(oracle, image)
function SNR = calc_snr(oracle, image)
num = sum(sum(oracle.*oracle));
den = sum(sum( (oracle - image).^2 ));

SNR = 10 * log10(num / den);

end


% function PSNR = psnr(f1, f2)
function PSNR = calc_psnr(oracle, image)
%计算两幅图像的峰值信噪比
% k = 8;
% %k为图像是表示地个像素点所用的二进制位数，即位深。
% fmax = 2.^k - 1;
% a = fmax.^2;
% e = double(f1) - double(f2);
% [m, n] = size(e);
% e=e(:);
% b = sum(e.^2);
% PSNR = 10*log(m*n*a/b);
[nx, ny] = size(image);
factor = 1/(nx * ny);
num = max(oracle(:)) ^2;

den = factor * sum(sum( (oracle - image).^2 ));

PSNR = 10* log10(num / den);

end


function RMSE = calc_rmse(oracle, image)
[nx, ny] = size(image);
RMSE = sqrt(1.0 / (nx * ny) * sum(sum((oracle - image).^2 )));
end


function MAE = calc_mae(oracle, image)
[nx, ny] = size(image);
MAE = 1.0 / (nx * ny) * sum(sum(abs(oracle - image)));
end


function [ssim_map, mssim] = compute_map_ssim(image1, image2, window_size, sigma)
%Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, "Image
%quality assessment: From error visibility to structural similarity"
%IEEE Transactios on Image Processing, vol. 13, no. 4, pp.600-612,
%Apr. 2004.

% K: constants in the SSIM index formula (see the above reference).
% # Select the parameters C1, C2 to stabilize the calculation
K = [0.01 0.03];
L = 0.5 * (abs(max(image1(:)) - min(image1(:))) + abs(max(image2(:)) - min(image2(:))));

C1 = (K(1)*L)^2;    % 计算C1参数，给亮度L（x，y）用。  
C2 = (K(2)*L)^2;    % 计算C2参数，给对比度C（x，y）用。
fprintf('\nParameters to calculate the SSIM indeces:')
fprintf('K = [%f, %f], L = %f, C1 = %f, C2 = %f\n', K(1), K(2),L, C1, C2);
gauss_kernel = fspecial('gaussian', window_size, sigma);        % 标准偏差1.5，11*11的高斯低通滤波
gauss_kernel = gauss_kernel/sum(sum(gauss_kernel)); %滤波器归一化操作。

% gauss_kernel = init_gaussian_kernel(window_size, sigma);

% image_mean1 = scim.filters.convolve(image1, gauss_kernel);
% image_mean2 = scim.filters.convolve(image2, gauss_kernel);
mu1   = filter2(gauss_kernel, image1, 'valid');
mu2   = filter2(gauss_kernel, image2, 'valid');

mu1_sq = mu1.*mu1;     % 计算出Ux平方值。
mu2_sq = mu2.*mu2;     % 计算出Uy平方值。
mu1_mu2 = mu1.*mu2;    % 计算Ux*Uy值。

sigma1_sq = filter2(gauss_kernel, image1.*image1, 'valid') - mu1_sq;  % 计算sigmax （标准差）
sigma2_sq = filter2(gauss_kernel, image2.*image2, 'valid') - mu2_sq;  % 计算sigmay （标准差）
sigma12 = filter2(gauss_kernel, image1.*image2, 'valid') - mu1_mu2;   % 计算sigmaxy（标准差）
if (C1 > 0 && C2 > 0)
   ssim_map = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))./((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2));
else
   numerator1 = 2*mu1_mu2 + C1;
   numerator2 = 2*sigma12 + C2;
   denominator1 = mu1_sq + mu2_sq + C1;
   denominator2 = sigma1_sq + sigma2_sq + C2;
   ssim_map = ones(size(mu1));
   index = (denominator1.*denominator2 > 0);
   ssim_map(index) = (numerator1(index).*numerator2(index))./(denominator1(index).*denominator2(index));
   index = (denominator1 ~= 0) & (denominator2 == 0);
   ssim_map(index) = numerator1(index)./denominator1(index);
end
mssim = mean2(ssim_map);

end