% Test of wazdenoise.m
% CCK 2018-Mar-23, based on test2.m

% Set parameters for applying noise and for denoising
sigma_noise = 20; % standard deviation of the noise
Nscales = 8;      % number of wavelet scales
Naz = 8;          % number of azimuths
SNR = 8;          % signal-to-noise ratio
percentile = 99.9;  % 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
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, sqrt(variance/SNR), 2048, 2048);

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

figure(3,'position',[1,1,2*768,768]) % How does it look??
   % For some reason, octave's position keyword has the vertical window size and then the horizontal size.
%figure();
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
axis tight
title('Residuals')