%% ratedistortion.m
% show image and family of thresholded images
% show the rate distorion curve for the family of images
% use the DFT

% S. Allen Broughton

%% get image
%load clown
load mandrill
%X = zeros(128,128); X(61:80,41:80)=1;

%% find parameters and DFT of X
X= X-mean(X(:));     % eliminate DC coefficient
[m,n] = size(X);
N = m*n;
fX = fft2(X);        % get  2D FFT of image
afX = abs(fX);       % a convenient variable
M = max(abs(fX(:))); % get maximum amplitude;


%% set up "movie"
plevels = [0, (1:9)/1000 (1:9)/100 ,(1:9)/10, 1:9 ];
E = energy(X);
Dist = [ ];
PRC = [ ];
eps = 10^(-5);

%% movie

for pL = plevels 
      PercentThresholdofMaximum = pL;
      cutlevel = pL*M/100;
      J = find(afX <cutlevel);
      fXc = fX; fXc(J) = 0;
      sfXc = fftshift(fXc);
      Xc = real(ifft2(fXc));
      PercentRetainedCoefficients = 100*(N-length(J))/N;
      PercentDistortion  = 100*abs(energy(X-Xc)/E);
      Dist = [PercentDistortion, Dist];
      PRC = [PercentRetainedCoefficients,PRC];
      'PercentThresholdofMaximum, PercentRetainedCoefficients, PercentDistortion ='
      [PercentThresholdofMaximum,PercentRetainedCoefficients,PercentDistortion]
 
      figure(1)
      imagesc(Xc), colormap(gray)
      title(' image at cutlevel ');
      
      figure(2)
      imagesc(log(1+abs(sfXc))), colormap(gray)
      title('shifted log(1+|dft|) at cutlevel')
      
      figure(3)
      subplot(2,1,1)
      plot(log10(Dist+eps),log10(PRC+eps),'o')
      xlabel('log10 of % Distortion')
      ylabel('log10 of % Retained Coeffs')
      subplot(2,1,2)
      plot( Dist, PRC ,'o')
      xlabel('% Distortion')
      ylabel('% Retained Coeffs')
      
      pause(1)
  end;