% Compute discrete wavelet transform of an image.
%
% Boundaries are taken into account simply by zero extension, leading to
% artefacts near the boundaries.
%
% Arguments:
% h       low-pass wavelet filter, horizontal vector
% g       high-pass wavelet filter, horizontal vector
% im      image, a real-valued MxN square matrix
% D       number of scales used; M and N must both be divisible by 2^D
%
% Returns:
% tr   wavelet coefficients as a transformed image of size NxN
%
% Calls to: wavetrans2Donce.m
%
% Samuli Siltanen June 2016

function tr = wavetrans2D(h,g,im,D)

% Make sure that the convolution kernels h and g are horizontal vectors
h = h(:).';
g = g(:).';

% Record image size and perform checks
tmp = size(im);
row = tmp(1);
col = tmp(2);
if length(tmp)>2
    error('Error: Input image "im" is not monochromatic')
end
if (mod(col,2^D)>0) | (mod(row,2^D)>0)
    error('Error: Image size is not divisible by 2^D')
end

% Initialize transform image
tr = zeros(size(im));

% Loop over transformation levels (scales)
tmpim = im;
for ddd = 1:D
    
    % Form the one-step transform image at the current scale
    [hsmvsm, hdevsm, hsmvde, hdevde] = wavetrans2Donce(h, g, tmpim);
    tmptr = zeros(row/2^(ddd-1),col/2^(ddd-1));
    tmptr(1:end/2,1:end/2) = hsmvsm;
    tmptr((end/2+1):end,1:end/2) = hdevsm;
    tmptr(1:end/2,(end/2+1):end) = hsmvde;
    tmptr((end/2+1):end,(end/2+1):end) = hdevde;
    
    % Update the transform image
    tr(1:row/2^(ddd-1), 1:col/2^(ddd-1)) = tmptr;
    
    % Denote now recursively the low-passed quadrant as tmpim
    tmpim = hsmvsm;
end


