% Test of denoise.m
% CCK 2018-Feb-28

% Set parameters for applying noise and for denoising
sigma_noise = 20; % standard deviation of the noise
Nscales = 6;      % number of wavelet scales
SNR = 8;          % signal-to-noise ratio
percentile = 99;  % significance threshold for keeping wavelet coefficients (%)

% Load up an interesting image
image_big = double( imread('Ansel_Adams_and_camera.jpg') );

% Downsample and crop
step = 4; % downsample factor
image_filt = conv2(image_big, ones(step)/step^2, 'valid');
image = image_filt(1:step:3072, 1:step:2048); 

isize = size(image);
Nx = isize(1);
Ny = isize(2);


%{
figure(1) % The true image
hold off
imagesc(image)
colormap(gray)
axis equal
title('True Image')
%}

% Add noise
%{
noise1 = normrnd(0, 30, Nx, Ny);
image_noisy = image + noise1;
%}

% Add noise
variance = var(image(:));
   % image(:) looks like a vector (1D array) to var.
noise1 = normrnd(0, sqrt(variance/SNR), Nx, Ny);
image_noisy = image + noise1;

%{
figure(2) % The noisy image
hold off
imagesc(image_noisy)
colormap(gray)
axis equal
title('Noisy Image')
%}

% Create separately a sample of noise from the same distribution
noise2 = normrnd(0, 30, 2048, 2048);

% Wavelet denoising
foo = denoise(image_noisy, noise2, Nscales, percentile);

figure(3) % How does it look??
hold off
imagesc([image, image_noisy, foo])
colormap(gray)
axis equal
axis tight
title('True (left); Noisy (middle); Reconstructed (right)')

figure(4) % The residuals
hold off
imagesc(image_noisy - foo)
colormap(gray)
axis equal
title('Residuals')