%Chapter 8,  Fig 8.5
%Eric Dubois, updated 2018-12-14
% create two noisy images and compare CPSNR and sCIELAB error
% add two varieties of noise with similar CPSNR but different sCIELAB error
%Data file should be placed in a data directory on the MATLAB path.
%requires function mdsp_scielab_error
clear all, close all;
RGBin = im2double(imread('05.tif'));
figure, imshow(RGBin);
%Assume image is sRGB. Convert to YIQ.
YIQ = rgb2ntsc(RGBin);
%version 1: add lowpass noise to Y, I, and Q
% create a Gaussian filter to lowpass filter the noise in each channel
L = 7;
x = -L:L; y = x;
[X, Y] = meshgrid(x,y);
r0  = 2;
hG = exp(-(X.^2 + Y.^2)/(r0^2)); % unnormalized
hG = hG/(sum(sum(hG))); % normalize for unit DC gain
wt1 = [.2 .05 .05];
for C = 1:3
   YIQ_n1(:,:,C) = YIQ(:,:,C) + wt1(C)*imfilter(randn([size(RGBin,1),size(RGBin,2)]),hG,'same');
end
RGB_n1 = ntsc2rgb(YIQ_n1);
figure, imshow(RGB_n1);

%compute CMSE and CPSNR
err1 = RGBin - RGB_n1;
RMSE1 = sqrt(sum(sum(sum(err1.^2)))/(size(RGBin,1)*size(RGBin,2)));
PSNR1 = 20*log10(1/RMSE1);
%compute sCIELAB error
sCIELABerr1 = mdsp_scielab_error(RGBin,RGB_n1);
sCIELABPSNR1 = 20*log10(100/sCIELABerr1);
%version 2: add highpass Gaussian noise to Y, I and Q
hpG = zeros(size(hG));
hpG(L+1,L+1) = 1;
hpG = hpG - 0.5*hG;
wt2 = [.006 .031 .031];
for C = 1:3
   YIQ_n2(:,:,C) = YIQ(:,:,C) + wt2(C)*imfilter(randn([size(RGBin,1),size(RGBin,2)]),hpG,'same');
end
RGB_n2 = ntsc2rgb(YIQ_n2);

figure, imshow(RGB_n2);
%compute CMSE and CPSNR
err2 = RGBin - RGB_n2;
RMSE2 = sqrt(sum(sum(sum(err2.^2)))/(size(RGBin,1)*size(RGBin,2)));
PSNR2 = 20*log10(1/RMSE2);
sCIELABerr2 = mdsp_scielab_error(RGBin,RGB_n2);
sCIELABPSNR2 = 20*log10(100/sCIELABerr2);