% Lezione 7
% Trasformata di Fourier applicata alla compresisone di immagini
% Immagine Lenna di Playboy
%
vnoise=20;
% Lettura binaria
fid=fopen('lenna.bin','r');
%a=fread(fid,'uint8');
a=fread(fid,'uint8=>single');
fclose(fid);

NN=length(a);
N=sqrt(NN);
a=reshape(a,N,N);
an=a+vnoise*randn(size(a));

figure(1)
set(gcf,'Name','Immagine di Lenna')
subplot(1,2,1)
imagesc(a)
colormap(gray)
axis equal
axis off
subplot(1,2,2)
imagesc(an)
colormap(gray)
axis equal
axis off

fa=fft2(a);
fan=fft2(an);

fa2=fa;
fa2(1,1)=0;
fan2=fan;
fan2(1,1)=0;

figure(3)
clf;
set(gcf,'Name','Lenna - Trasformata di Fourier')
subplot(1,2,1)
imagesc(abs(fa2))
axis equal
axis off
colorbar
subplot(1,2,2)
imagesc(abs(fan2))
axis equal
axis off
colorbar

figure(4)
clf;
set(gcf,'Name','Istogramma Lenna')
subplot(1,2,1)
hist(abs(fa2(:)),500)
axis tight
title('Istogramma FT Lenna')
subplot(1,2,2)
hist(abs(fan2(:)),500)
axis tight
title('Istogramma FT Lenna noise')

[fansort indsort]=sort(abs(fan(:)),1,'descend');
figure(5); clf;
set(gcf,'Name','Lenna denoised')

%cr1=0; cr2=90;
cr1=90; cr2=99;
cr=linspace(cr1,cr2,10);
error=zeros(10,1)';
for i=1:10
    compression=cr(i);
    indmax=round((1-compression/100)*NN);
    far=complex(zeros(NN,1),zeros(NN,1));
    far(indsort(1:indmax))=fan(indsort(1:indmax));
    far=reshape(far,N,N);
    far=ifft2(far);
    far=real(far);
    far=min(far,255);
    far=max(far,0);
    subplot(2,5,i)
    imagesc(far)
    colormap(gray)
    axis equal
    axis off
    title(['Lenna - ' num2str(compression) '%'])
    error(i)=sqrt(sum((a(:)-far(:)).^2)/NN);
    fprintf('Rapp. compr: %f%% - Errore: %f\n',...
        compression,error(i))
end

figure(6)
set(gcf,'Name','Curva di errore')
plot(cr,error)
[errmin indmin]=min(error);
hold on
plot(cr(indmin),errmin,'*r')
xlabel('Rapporto di compressione [%]')
ylabel('RMSE')



% Rimozione rumore da funzione
vnoise=0.3;
NN=2^10;
x=linspace(0,2*pi,NN)';
a=sin(x).*cos(x).^2;
an=a+vnoise*randn(size(a));

figure(11)
set(gcf,'Name','Funzione')
plot(x,a,x,an,'.')

fan=fft(an);

[fansort indsort]=sort(abs(fan(:)),1,'descend');
figure(15); clf;
set(gcf,'Name','Lenna denoised')

%cr1=0; cr2=90;
cr1=99; cr2=99.9;
cr=linspace(cr1,cr2,10);
error=zeros(10,1)';
for i=1:10
    compression=cr(i);
    indmax=round((1-compression/100)*NN);
    far=complex(zeros(NN,1),zeros(NN,1));
    far(indsort(1:indmax))=fan(indsort(1:indmax));
    far=ifft(far);
    far=real(far);
    subplot(2,5,i)
    plot(x,far)
    hold on
    plot(x,a,'r')
    plot(x,an,'.','MarkerSize',1)
    axis tight
    title(['sin(x)cos^2(x) - ' num2str(compression) '%'])
    error(i)=sqrt(sum((a(:)-far(:)).^2)/NN);
    fprintf('Rapp. compr: %f%% - Errore: %f\n',...
        compression,error(i))
end

figure(6); clf;
set(gcf,'Name','Curva di errore')
plot(cr,error)
[errmin indmin]=min(error);
hold on
plot(cr(indmin),errmin,'*r')
xlabel('Rapporto di compressione [%]')
ylabel('RMSE')