name = 'boat';
n0 = 512;
f = rescale( load_image(name,n0) );

imageplot( f, 'Image f');
pause;

% transformacja Fouriera z normalizacją
F = fft2(f) / n0;

disp(strcat(['Energy of Image:   ' num2str(norm(f(:)))]));
disp(strcat(['Energy of Fourier: ' num2str(norm(F(:)))]));
pause;

L = fftshift(log( abs(F)+1e-1 ));

% liniowa aproksymacja w bazie Fouriera
M = n0^2/64;

m = sqrt(M);
F = fftshift(fft2(f));
F1 = zeros(n0,n0);
sel = (n0/2-m/2:n0/2+m/2)+1;
F1(sel,sel) = F(sel,sel); 
fM = real( ifft2(fftshift(F1)) );
% display
clf;
imageplot(clamp(fM), ['Linear, Fourier, SNR=' num2str(snr(f,fM), 4) 'dB']);
pause;

% porownanie profili
clf;
subplot(2,1,1);
plot(f(:,n0/2));
axis('tight'); title('f');
subplot(2,1,2);
plot(fM(:,n0/2));
axis('tight'); title('f_M');
pause;

% nieliniowa aproksymacja w bazie Fouriera
T = .2;
F = fft2(f) / n0;
FT = F .* (abs(F)>T);
fM = real( ifft2(FT)*n0 );
clf;
imageplot(clamp(fM), ['Non-linear, Fourier, SNR=' num2str(snr(f,fM), 4) 'dB']);
pause;

m = sum(FT(:)~=0);
disp(['M/N = 1/'  num2str(round(n0^2/m)) '.']);
pause;

% nieliniowa aproksymacja w bazie Fouriera
F = fft2(f);
a = sort(abs(F(:))); 
if a(1)<a(n0*n0)
	a = a(n0^2:-1:1);
end
T = a(M+1);
F = F .* (abs(F)>T);
fM = real( ifft2(F) );
% display
clf;
imageplot(clamp(fM), ['Non-linear, Fourier, SNR=' num2str(snr(f,fM), 4) 'dB']);

% porownanie profili
clf;
subplot(2,1,1);
plot(f(:,n0/2));
axis('tight'); title('f');
subplot(2,1,2);
plot(fM(:,n0/2));
axis('tight'); title('f_M');
pause;

% aproksymacje w bazach falkowych
Jmin = 0;

fw = perform_wavelet_transf(f,Jmin,+1);

% aproksymacja liniowa w bazach falkowych
q = sqrt(M);
fw = perform_wavelet_transf(f,Jmin,+1);
fw1 = zeros(n0,n0);
fw1(1:q,1:q) = fw(1:q,1:q);
fM = perform_wavelet_transf(fw1,Jmin,-1);
% display
clf;
imageplot(clamp(fM), ['Linear, Wavelets, SNR=' num2str(snr(f,fM), 4) 'dB']);
pause;

% porownanie profili
clf;
subplot(2,1,1);
plot(f(:,n0/2));
axis('tight'); title('f');
subplot(2,1,2);
plot(fM(:,n0/2));
axis('tight'); title('f_M');
pause;

% nieliniowa aproksymacja w bazach falkowych
T = .2;
fwT = fw .* (abs(fw)>T);
fM = perform_wavelet_transf(fwT,Jmin,-1);
clf;
imageplot(clamp(fM), strcat(['Approximation, SNR=' num2str(snr(f,fM),3) 'dB']));
pause;

% nieliniowa aproksymacja w bazach falkowych
%fw = perform_wavelet_transf(f,Jmin,+1);
a = sort(abs(fw(:))); 
if a(1)<a(n0*n0)
	a = a(n0^2:-1:1);
end
T = a(M+1);
fw1 = fw .* (abs(fw)>T);
fM = perform_wavelet_transf(fw1,Jmin,-1);
% display
clf;
imageplot(clamp(fM), ['Non-linear, wavelets, SNR=' num2str(snr(f,fM), 4) 'dB']);
pause;

% porownanie profili
clf;
subplot(2,1,1);
plot(f(:,n0/2));
axis('tight'); title('f');
subplot(2,1,2);
plot(fM(:,n0/2));
axis('tight'); title('f_M');