% An implementation of Algorithme Á Trous.
function W = atrous(image, Nscales)
% Returns the wavelet transform of image, divided into Nscales
% layers, plus a final remainder that has all the larger scales
% up to DC. 
%
% The image is treated with periodic boundary conditions, with
% no attempt to ameliorate edge effects. If you are concerned
% about edge effects, try unfolding/mirroring the image or
% padding it in some way.
% 
% The inverse transform is the sum of all layers:
%     image = sum(W, 3)
% CCK 2018-Feb-28


isize = size(image);
Nx = isize(1);
Ny = isize(2);
k1 = zeros(Nx,Ny);
k1(1:3,1:3) = [1,2,1; 2,4,2; 1,2,1]/16;
k1 = circshift(k1, [-1, -1]);

W = zeros(Nx,Ny,Nscales+1);



% Necessary Fourier transforms. Note the F_ prefix notation.
F_image = fft2(double(image));
   % Since image data types are often uint8, double() makes sure
   % we are working with floating point numbers.
F_k1 = fft2(k1);

F_image_blur = F_image;
for j=1:Nscales
   F_image_blur_next = F_image .* (  F_k1.^(2.^(j-1))  );
   W(:,:,j) = real( ifft2( F_image_blur - F_image_blur_next ) );
   F_image_blur = F_image_blur_next;
end
% I could have vectorized the above for loop, but it would save
% practically no time and would greatly increase memory usage.

% The final layer:
W(:,:,Nscales+1) = real( ifft2( F_image_blur ) );