% Denoising of images by wavelet significance testing.
% CCK 2018-Feb-28
function reconstructed = denoise(image, noise, Nscales, percent_significance)
% image and noise are 2d arrays, not necessarily of the same size.
% It is assumed that the input image contains additive noise, 
% of which the noise array is a representative sample. 
% It is particularly helpful to have a large noise array, to get good statistics.
% Nscales is passed straight throught to atrous().
% percent_significance is the percentile threshold for keeping wavelet coefficients.
% The result is a reconstructed image, hopefully with less noise.
% No effort is made to prevent artifacts at the edges of the image.


nsize = size(noise);
Nx_noise = nsize(1);
Ny_noise = nsize(2);

% Calculate wavelet transforms and wavelet power spectra.
W_noise = atrous(noise, Nscales);
W_image = atrous(image, Nscales);

% Significance filter layer by layer.
for j = 1:Nscales
   W_noise_power = W_noise(:,:,j) .^ 2; % Layer j of the noise wavelet power spectrum.
   W_image_power = W_image(:,:,j) .^ 2; % Layer j of the image wavelet power spectrum.
   threshold = prctile( reshape(W_noise_power, [1, Nx_noise * Ny_noise]), percent_significance );
      % Calculate the threshold power for layer j based on the statistics of the noise image.
   W_image(:,:,j) = W_image(:,:,j) .* ( W_image_power > threshold );
end

reconstructed = sum( W_image, 3 );