n = 512;
%f = rescale( load_image('target256', n) );
f = double(imread('target256.pgm'));

Mlist = round( [.01 .05]*n^2 );

% obraz oryginalny
clf;
imshow(f,[]); title('source image');%imageplot(f);
pause;

% aproksymacja nieliniowa w bazie Fouriera
fF = fft2(f)/n;
% T = .3;
% c = fF .* (abs(fF)>T);
% fM = real(ifft2(c)*n);
% imageplot(clamp(fM), ['Non-linear, Fourier, SNR=' num2str(snr(f,fM), 4) 'dB']);

clf;
for i=1:length(Mlist)
    M = Mlist(i);
    fFT = perform_thresholding(fF,M,'largest');
    fM = real( ifft2(fFT)*n );
    %imageplot( clamp(fM), ['M/N=' num2str(M/n^2,2) ', SNR=' num2str(snr(f,fM),3) 'dB'], 1,2,i);
    imshow(clamp(fM),[]);
end
pause;

cR = sort(abs(fF(:)));
if cR(n^2)>cR(1)
    cR = reverse(cR); % be sure it is in reverse order
end

% blad aproksymacji w bazie Fouriera
lw = 2;
err_fft = max( norm(f(:))^2 - cumsum(cR.^2), 1e-10); 
clf;
h = plot(log10(err_fft / norm(f(:))^2)); 
if using_matlab()  
    set(h, 'LineWidth', lw);
end
ax = [1 n^2/50 -2.35 0];
axis(ax);
title('log_{10}( \epsilon^2[M]/||f||^2 )');
pause;

% nieliniowa aproksymacja w bazach falkowych
% Jmin = 0;
% fW = perform_wavelet_transf(f,Jmin,+1);

Jmin = 1;
options.h = compute_wavelet_filter('Daubechies',10);
fW = perform_wavortho_transf(f,Jmin,+1, options);

clf;
for i=1:length(Mlist)
    M = Mlist(i);
    fWT = perform_thresholding(fW,M,'largest');
    fM = perform_wavortho_transf(fWT,Jmin,-1, options);
    imageplot( clamp(fM), ['M/N=' num2str(M/n^2,2) ', SNR=' num2str(snr(f,fM),3) 'dB'], 1,2,i);
end
pause;

cR = sort(abs(fW(:)));
if cR(n^2)>cR(1)
    cR = reverse(cR); % be sure it is in reverse order
end

% blad aproksymacji w bazach falkowych
err_wav = max( norm(f(:))^2 - cumsum(cR.^2), 1e-10);
clf; hold('on');
h = plot(log10(err_fft / norm(f(:))^2), 'r'); 
if using_matlab()  
    set(h, 'LineWidth', lw);
end
h = plot(log10(err_wav / norm(f(:))^2), 'b'); 
if using_matlab()  
    set(h, 'LineWidth', lw);
end
axis(ax);
title('log_{10}( \epsilon^2[M]/||f||^2 )');
legend('Fourier', 'Wavelets');
if using_matlab()  
    box('on');
end
pause; 

% nieliniowa aproksymacja w bazie DCT
fC = dct2(f);

clf;
for i=1:length(Mlist)
    M = Mlist(i);
    fCT = perform_thresholding(fC,M,'largest');
    fM = idct2(fCT);
    imageplot( clamp(fM), ['M/N=' num2str(M/n^2,2) ', SNR=' num2str(snr(f,fM),3) 'dB'], 1,2,i);
end
pause;

cR = sort(abs(fC(:)));
if cR(n^2)>cR(1)
    cR = reverse(cR); % be sure it is in reverse order
end

% blad aproksymacji w bazie DCT
err_dct = max( norm(f(:))^2 - cumsum(cR.^2), 1e-10);
clf; hold('on');
h = plot(log10(err_fft / norm(f(:))^2), 'r'); 
if using_matlab()  
    set(h, 'LineWidth', lw);
end
h = plot(log10(err_dct / norm(f(:))^2), 'b'); 
if using_matlab()  
    set(h, 'LineWidth', lw);
end
axis(ax);
title('log_{10}( \epsilon^2[M]/||f||^2 )');
if using_matlab()  
    box('on');
end
legend('Fourier', 'DCT');
pause;

% nieliniowa aproksymacja w lokalnych bazach DCT
w = 16;

fL = zeros(n,n);

for i=1:n/w
    for j=1:n/w
        seli = (i-1)*w+1:i*w;
        selj = (j-1)*w+1:j*w;
        fL(seli,selj) = dct2( f(seli,selj) );
    end
end

clf;
for u=1:length(Mlist)
    M = Mlist(u);
    fLT = perform_thresholding(fL,M,'largest');
    % inverse
    fM = fLT;
    for i=1:n/w
        for j=1:n/w
            seli = (i-1)*w+1:i*w;
            selj = (j-1)*w+1:j*w;
            fM(seli,selj) = idct2( fM(seli,selj) );
        end
    end
    % display
    imageplot( clamp(fM), ['M/N=' num2str(M/n^2,2) ', SNR=' num2str(snr(f,fM),3) 'dB'], 1,2,u);
end
pause;

cR = sort(abs(fL(:)));
if cR(n^2)>cR(1)
    cR = reverse(cR); % be sure it is in reverse order
end

% blad aproksymacji w lokalnych bazach DCT
err_ldct = max( norm(f(:))^2 - cumsum(cR.^2), 1e-10); 
clf; hold('on');
h = plot(log10(err_fft / norm(f(:))^2), 'r'); 
if using_matlab()  
    set(h, 'LineWidth', lw);
end
h = plot(log10(err_dct / norm(f(:))^2), 'g'); 
if using_matlab()  
    set(h, 'LineWidth', lw);
end
h = plot(log10(err_ldct / norm(f(:))^2), 'b'); 
if using_matlab()  
    set(h, 'LineWidth', lw);
end
h = plot(log10(err_wav / norm(f(:))^2), 'k'); 
if using_matlab()  
    set(h, 'LineWidth', lw);
end
axis(ax);
title('log_{10}( \epsilon^2[M]/||f||^2 )');
legend('Fourier', 'DCT', 'local-DCT', 'Wavelets');
if using_matlab()  
    box('on');
end