function task3()

% create phantom image.
imagesize = 256;
img = phantom(imagesize);
% Later: load the 'magic.png' image, and resize it to a power of 2, for simplicity
% img = double(imread('magic.png'));
% img = imresize(img, [imagesize, imagesize], 'bilinear');

% Later: add noise to the image
% noise = randn(imagesize, imagesize) * 0.05;
% TODO

figure;
subplot(2, 2, 1);
imagesc(img);
title('Original');
colormap gray;
axis image;

% Compute a single step of the 2-D Haar transform of the image
% [approx, detailH, detailV, detailD] = TODO

% compose an image from the four sets of coefficients, aligned as follows:
% |-------------------|
% | approx  | detailH |
% |-------------------|
% | detailV | detailD |
% |-------------------|
% waveletcoeffs = [TODO];

subplot(2, 2, 2);
imagesc(waveletcoeffs);
title('Wavelet Coefficients')
colormap gray;
axis image;

% compute the inverse 2-D Haar transform
% recon = TODO;

subplot(2, 2, 3)
imagesc(recon);
title('Reconstruction')
colormap gray;
axis image;
subplot(2, 2, 4)
imagesc(recon-img);
title('Difference to original image')
colormap gray;
axis image;

% Compute ndec-many levels of decomposition using the 2-D Haar wavelet transform
ndec = 3;
% waveletcoeffs = TODO;

% ----------- processing ---------------

% compute the number of pixels in the original image that are not zero
% TODO
% compute the number of wavelet coefficients that are not zero (later: higher
% than the chosen threshold, respectively)
% TODO
% set all coefficients that are lower than the chosen threshold to zero
% TODO
% compute the compression rate if only values that are not zero (later: higher
% than the chosen threshold, respectively) have to be stored/transmitted?
% compression = TODO

% -------- ende processing ------------

figure;
imagesc(waveletcoeffs);
colormap gray;
axis image;

% compute the reconstruction from the (modified) wavelet coefficients
% recon = TODO;
figure;
imagesc(recon);
colormap gray;
axis image;
title(['compression = ' num2str(compression)]);

% compute the mean quadratic error between original and reconstructed image
% mse = TODO;

figure;
imagesc(recon-img);
colormap gray;
axis image;
title(['difference, mse = ' num2str(mse)]);

end

%--------------------------------------------------------------------------
% 2-D DHT
%--------------------------------------------------------------------------
function [approx, detailH, detailV, detailD] = dht2(signal)

% definition of high pass and low pass filter for Haar analysis
lp = 1/sqrt(2) * [ 1 1];
hp = 1/sqrt(2) * [-1 1];

% ------------rows-------------
% 1-D DHT along the rows
% Periodically supplement the signal for one element
% sig = TODO;

% Compute convolutions with the high pass and low pass filter
% approximation = TODO;
% detail = TODO;

% downsample the coefficients by a factor of 2
% approximation = approximation(:, TODO);
% detail = detail(:, TODO);

% ------------columns-------------
% Operate first on approximation coefficients:
  % Periodically supplement the signal for one element
  % sig = [TODO];
  
  % Compute convolutions with the high pass and low pass filter
  % (hint: lp' is the transposed vector of lp)
  % approx = TODO;
  % detailH = TODO;
  
  % downsample the coefficients by a factor of 2
  % approx = TODO;
  % detailH = TODO;

% Operate now on the detail coefficients:
  % Periodically supplement the signal for one element
  % sig = TODO;

  % Compute convolutions with the high pass and low pass filter
  % (hint: lp' is the transposed vector of lp)
  % detailV = TODO;
  % detailD = TODO;

  % downsample coefficients by a factor of 2
  % detailV = TODO;
  % detailD = TODO;

end

%--------------------------------------------------------------------------
% 2-D IDHT
%--------------------------------------------------------------------------
function [recon] = idht2(approx, detailH, detailV, detailD)

% definition of high pass and low pass filter for Haar synthesis
rlp = 1/sqrt(2) * [1  1];
rhp = 1/sqrt(2) * [1 -1];


% ------------columns-------------
% compute the inverse DHT along the columns by combining approx and detailH
% into approx, and detailV and detailD into detail
%
% upsample the approximation and detail coefficients by a factor of 2 along the
% columns
% usapprox = zeros(TODO);
% usapprox(TODO) = TODO;
% usdetail = zeros(TODO);
% usdetail(TODO) = TODO

% periodically supplement both sets of coefficients by one element
% usapprox = TODO;
% usdetail = TODO;

% convolve the columns with high pass and low pass filters. Add the results.
% approx = TODO;

% ------------rows-------------
% upsample the coefficients by a factor of 2 along the rows
% usapprox = TODO;
% usdetail = TODO;

% periodically supplement both sets of coefficients by one element
% usapprox = TODO;
% usdetail = TODO;

% convolve the columns with high pass and low pass filters. Add the results.
% recon = TODO;

end

%--------------------------------------------------------------------------
% 2-D filter bank for multiresolution analysis
%--------------------------------------------------------------------------
function waveletcoeffs = mra2D(signal, ndec)

% initialize the wavelet coefficients with the signal
% if ndec = 0, then the wavelet decomposition is the signal itself
% waveletcoeffs = TODO;
sizex = size(signal,1);
sizey = size(signal,2);
for i=1:ndec

    % compute dht2 for the current approximation coefficients
%    [TODO] = TODO;

	% replace the approximation coefficients of the previous decomposition
	% level by the coefficients obtained at the current level
%     waveletcoeffs(TODO) = [TODO];
    
    % reduce the index limits by a factor of 2
%     sizex = TODO;
%     sizey = TODO;

end


end

%--------------------------------------------------------------------------
% 2-D filter bank for multiresolution synthesis
%--------------------------------------------------------------------------
function recon = mrs2D(waveletcoeffs, ndec)

sizex = size(waveletcoeffs,1)/(2^ndec);
sizey = size(waveletcoeffs,2)/(2^ndec);
for i=ndec:-1:1
    
    % compute inverse Haar transform
%     waveletcoeffs(TODO) = TODO
    
    % the number of approximation coefficients doubles after each iteration
%     sizex = TODO;
%     sizey = TODO;
    
end

% return result
%recon = TODO;

end
