function x2 = idwtdaub2(xdaub2, J, Hord)
% Inverse Daubechies wavelet transform in two dimensions
%  x2 = idwtdaub2(xdaub2, J, Hord)
%
% Inputs:
%  xdaub2   array with 2-dimensional Daubechies wavelet transform
%  J        number of stages (resolution levels)
%  Hord     positive integer. Order of Daubechies filter is 2*Hord-1. 
% Output:
%  x2       array with reconstructed 2-dimensional signal
%

sz = size(xdaub2);
nrows = sz(1); % number of rows in xdaub2
ncols = sz(2); % number of columns in xdaub2

% check whether nrows and ncols are divisable with 2^J
if ( rem(nrows, 2^J) ~= 0 ) || ( rem(ncols, 2^J) ~= 0 )
    error('The dimensions of x2 must be multiples of 2^J');
end

% compute Daubechies filter coefficients
H = hdaub(Hord);

x1 = zeros(nrows,ncols);
x2 = zeros(nrows,ncols);

% Inverse transform first column-wise, then row-wise
for j=1:ncols
    x1(:,j) = idwtdaub( xdaub2(:,j), J, H );
end

for i=1:nrows
    x2(i,:) = idwtdaub( x1(i,:), J, H );
end