clear
clc
close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1 Cargamos la imagen y definimos r
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f=imread('cameraman.tif');
f=double(f);
[M N]=size(f);
figure
imagesc(f),colormap(gray)
r=1/200;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2 Calculamos la transformada de Fourier y
%   el maximo y el minimo de su modulo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F=fft2(f)/sqrt(M*N);
aF=abs(F);
minF=min(min(aF));
maxF=max(max(aF));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3 Construimos y dibujamos la funcion h
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h=[];
for T=minF:200
    h=[h nnz(aF>=T)-r*M*N];
end
figure
plot(minF:200, h,[minF,200],[0 0]),xlabel('T')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4 Comprobamos el signo de h desde minF hasta maxF
%   Paramos cuando cambia el signo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h1=1;
T1=floor(minF);
while T1 < maxF && h1 > 0
    h1=nnz(aF>=T1)-r*M*N;
    T1=T1+1;
end
h1
T1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 5 Definimos la matriz de Fourier con el umbral aplicado
%   Dibujamos la imagen comprimida
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
G1=F.*(aF>T1);
g1=ifft2(G1)/sqrt(M*N);
figure
imagesc(g1),colormap(gray)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 6 Comprobamos la razon de compresion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
comp1=nnz(abs(G1))/nnz(aF)
r
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 7 Calculamos el valor umbral de corte 
%   utilizando el metodo de biseccion 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H=@(T) nnz(aF>=T)-r*M*N;
a=minF;b=200;
while b-a > 0.001
    c=(a+b)/2;
    if H(c)*H(b)<0
        a=c;
    elseif H(a)*H(c)<0
        b=c;
    else
        break
    end
end
c
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 7 Dibujamos de nuevo y comprobamos compresion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
G2=F.*(aF>c);
g2=ifft2(G2)/sqrt(M*N);
figure
imagesc(g2),colormap(gray)
comp1=nnz(abs(G2))/nnz(aF)
r