% Denoising of images by wazlet significance testing.
% CCK 2018-Mar-23, based on denoise.m
function reconstructed = wazdenoise(image, noise, Nscales, Naz, 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 and Naz are passed straight throught to wazlet().
% 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 = wazlet(noise, Nscales, Naz);
W_image = wazlet(image, Nscales, Naz);

% Significance filter layer by layer, and azimuth by azimuth.
for j = 1:Nscales
   for k = 1:Naz
      W_noise_power = W_noise(:,:,j,k) .^ 2; % Layer j of the noise wavelet power spectrum.
      W_image_power = W_image(:,:,j,k) .^ 2; % Layer j of the image wavelet power spectrum.
      % threshold = prctile( reshape(W_noise_power, [1, Nx_noise * Ny_noise]), percent_significance );
         % The above syntax accomplishes the same as the more efficient line below.
         % The W_noise_power(:) notation represents the aray as one-dimensional. Without the (:), 
         % prctile would operate independently on the columns of its argument, returning a vector.
      threshold = prctile( W_noise_power(:), percent_significance );
         % Calculate the threshold power for layer j based on the statistics of the noise image.
      W_image(:,:,j,k) = W_image(:,:,j,k) .* ( W_image_power > threshold );
   end
end

reconstructed = sum( sum( W_image, 3 ), 4 );