function task4()

%close all;

% Load image (image size should be a power of 2)

% phantom image
imagesize = 256;
img_nonoise = phantom(imagesize)*255;

% White square on black background
%img_nonoise = zeros(imagesize);
%img_nonoise(2:imagesize/4,4:imagesize/4) = 255;

% Natural images
%img_nonoise = double(imread('cameraman.png'));
%img_nonoise = double(imread('magic.png'));
%img_nonoise = imresize(img_nonoise, [imagesize, imagesize], 'bilinear');
% Or use cameraman.png, lena.png or own images

% add noise
noise = randn(imagesize, imagesize) * 20;
img = img_nonoise + noise;

% compute the mean squared error between original image and noisy image
%mse_noisy = TODO

clim = [0, 255];
figure;
subplot(2,2,1);
imagesc(img,clim);
title('Original');
colormap gray;
axis image;


%-----------------------------------
% DWT with arbitrary wavelet filters
%-----------------------------------

wavelet = 'sym8';

% create high pass filter and low pass filter using wfilters
% lp and hp are for analysis, rlp and rhp for synthesis
%[lp,hp,rlp,rhp] = TODO

% compute 2-D wavelet transformation of the image
%[approx, detailH, detailV, detailD] = mydwt2(img, TODO);

% concatenate the four sets of coefficients to an image:
% |-------------------|
% | approx  | detailH |
% |-------------------|
% | detailV | detailD |
% |-------------------|
waveletcoeffs = [approx detailH; detailV detailD];

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

% compute the inverse 2-D wavelet transformation
%recon = myidwt2(approx, detailH, detailV, detailD, TODO);

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


%--------------
% 2-D denoising
%--------------

figure;
imagesc(img,clim);
title('Original');
colormap gray;
axis image;

% number of levels of decomposition
ndec = 3;

% thresholding weighting factor for denoising - the larger, the stronger
thrweight = 1;
% thresholding mode: 0 - hard thresholding, 1 - soft thresholding
thrmode = 1;

% computation of threshold: globally on all coefficients (la=0) or level-adaptive (la=1)
la = 0;


denoised_UT = denoise_UT(img, wavelet, ndec, thrweight, thrmode, la);

% compute mean squared error between original image and denoised image
%mse_UT = TODO
% compute peak signal to noise ratio
%psnr_UT = TODO

figure;
imagesc(denoised_UT,clim);
colormap gray;
axis image;
title(['Denoised - Global Threshold, PSNR = ' num2str(psnr_UT)]);


% denoising with level-adaptive threshold
la = 1;
denoised_LAT = denoise_UT(img, wavelet, ndec, thrweight, thrmode, la);
% compute the mean squared error between original image and denoised image
%mse_LAT = TODO
% compute peak signal to noise ratio
%psnr_LAT = TODO

figure;
imagesc(denoised_LAT,clim);
colormap gray;
axis image;
title(['Denoised - Level Adaptive Threshold, PSNR = ' num2str(psnr_LAT)]);

% denoising using interscale correlations
denoised_IC = denoise_IC(img, wavelet, ndec, thrweight, thrmode);
% compute the mean squared error between original image and denoised image
%mse_IC = TODO
% compute peak signal to noise ratio
%psnr_IC = TODO

figure;
imagesc(denoised_IC,clim);
colormap gray;
axis image;
title(['Denoised - Interscale Correlations, PSNR = ' num2str(psnr_IC)]);


end


%--------------------------------------------------------------------------
% Estimation of the standard deviation of the noise
%--------------------------------------------------------------------------
function sigma = estimateNoise(waveletcoeffs)

% use the HH detail coefficients of the first level of decomposition
%wHH = TODO

% compute the standard deviation of noise
%sigma = TODO

end


%--------------------------------------------------------------------------
% Denoising using the Universal Threshold
%--------------------------------------------------------------------------
function denoisedimg = denoise_UT(signal, wavelet, maxlevel, thrweight, thrmode, la)

% get the filter coefficients for the requested wavelet
%TODO

% compute multiresolution analysis up until maxlevel
%waveletcoeffs = TODO

% estimate the standard deviation of the noise
%sigma = TODO
coeffsize = size(waveletcoeffs);
denoised = waveletcoeffs; % initialisiere mit ungewichteten koeffizienten
threshold = 0;
if la
    % compute threshold for each level of decomposition
    % initialize all threshold values with zero 
    threshold = zeros(size(waveletcoeffs));
    for i = 1:maxlevel
        % compute the thresholds
        %TODO
    end
else
    % compute (global) universal threshold
    %threshold = TODO
end

% threshold may not be applied to the approximation!
% -> set threshold for area of approximation coefficients to 0
%TODO

% apply threshold to coefficients, using either hard or soft thresholding
switch(thrmode)
    case 0 % hard thresholding
        % all coefficients with modulus lower than threshold are set to 0
        %TODO
    case 1 % soft thresholding
		% all coefficients with modulus lower than threshold are set to 0; all
		% coefficients above threshold are reduced accordingly (in modulus)
        %TODO
end

% wavelet synthesis
%denoisedimg = TODO


end

%----------------------------------------
% Denoising using interscale correlations
%----------------------------------------
function denoisedimg = denoise_IC(signal, wavelet, maxlevel, thrweight, thrmode)

% get the filter coefficients for the requested wavelet
%TODO

% compute multiresolution analysis up until maxlevel
%waveletcoeffs = TODO

sigma = estimateNoise(waveletcoeffs)
coeffsize = size(waveletcoeffs);
denoised = waveletcoeffs; % initialisiere denoised coefficients with original coefficients
threshold = zeros(size(waveletcoeffs));
mask = zeros(size(waveletcoeffs));
% iterate over all decomposition levels and compute thresholds
for i = 1:maxlevel

    % horizontal details
    % compute x- and y- array ranges of horizontal detail coefficients
    % prev: coefficients of level l+1; cur: coefficients of level l
%    x_prev = TODO;
%    y_prev = TODO;
%    x_cur = TODO;
%    y_cur = TODO;
    % get coefficients of level l+1
%    hprev = TODO;
    hpreve = [hprev hprev(:,1); hprev(1,:) hprev(1,1)]; % periodic supplement
    % interpolate to larger resolution of current level l
%    hprevi = TODO;
    % convolve hprevi with box filter
    hprevi = conv2(TODO, [0.25 0.25; 0.25 0.25], 'valid');

	% get coefficients of level l
    hcur = denoised(x_cur, y_cur);

	% compute m_l according to Eqn. 7 of the worksheet
%    mask(x_cur,y_cur) = TODO

    % vertical details: perform steps analogously to horizontal detail coefficients
    % compute x- and y- array ranges of horizontal detail coefficients
    % prev: coefficients of level l+1; cur: coefficients of level l
%    x_prev = TODO;
%    y_prev = TODO;
%    x_cur = TODO;
%    y_cur = TODO;
%    TODO

	% compute m_l according to Eqn. 7 of the worksheet
%    mask(x_cur, y_cur) = TODO

    % diagonal details: perform steps analogously to horizontal detail coefficients
    % compute x- and y- array ranges of horizontal detail coefficients
    % prev: coefficients of level l+1; cur: coefficients of level l
%    x_prev = TODO;
%    y_prev = TODO;
%    x_cur = TODO;
%    y_cur = TODO;
%    TODO
    
    % bilde produkt
%    mask(x_cur,y_cur) = TODO

    % computation of the thresholds - use now sigma.^2 instead of sigma
%    threshold(1:coeffsize(1)/(2^(i-1)),1:coeffsize(2)/(2^(i-1))) = TODO

end
% threshold may not be applied to the approximation!
% -> set threshold for area of approximation coefficients to 0
%TODO

% figure; imagesc(threshold); axis image;
figure; imagesc(mask); colormap gray; axis image;

% apply threshold to coefficients, using either hard or soft thresholding
switch(thrmode)
    case 0 % hard thresholding
        % all coefficients with mask modulus lower than threshold are set to 0
%        TODO
    case 1 % soft thresholding
		% all coefficients with mask modulus lower than threshold are set to 0;
		% all coefficients above threshold are reduced accordingly relative to
		% threshold (in modulus)
%        TODO
end

% wavelet synthesis
%denoisedimg = TODO


end

%--------------------------------------------------------------------------
% 2-D wavelet transform
% use your dwt method here from previous exercises, but adapt it such that
% - low pass and high pass filters are now passed as arguments hp and lp
% - filter length may be longer than 2, i.e. periodically supplement by more
%   elements
%--------------------------------------------------------------------------
function [approx, detailH, detailV, detailD] = mydwt2(signal, lp, hp)
	% TODO
	approx = 0;
	detailH = 0;
	detailV = 0;
	detailD = 0;
end


%--------------------------------------------------------------------------
% 2-D inverse wavelet transform
% use your idwt method here from previous exercises, but adapt it such that
% - low pass and high pass filters are now passed as arguments hp and lp
% - filter length may be longer than 2, i.e. periodically supplement by more
%   elements
%--------------------------------------------------------------------------
function [recon] = myidwt2(approx, detailH, detailV, detailD, rlp, rhp)
	% TODO
	recon = 0;
end

%--------------------------------------------------------------------------
% 2-D Multiresolution Analysis
% use your mra method here from previous exercises, but adapt it such that low
% pass and high pass filters are now passed as arguments hp and lp
%--------------------------------------------------------------------------
function waveletcoeffs = mra2(signal, ndec, lp, hp)

% initialize wavelet coefficients with signal
	waveletcoeffs = 0;
	% TODO
end

%--------------------------------------------------------------------------
% 2-D Multiresolution Synthesis
% use your mrs method here from previous exercises, but adapt it such that low
% pass and high pass filters are now passed as arguments hp and lp
%--------------------------------------------------------------------------
function recon = mrs2(waveletcoeffs, ndec, rlp, rhp)

	% TODO
	recon = 0;
end

