%% Ayodele Embry %% EE368B
%% PROJECT
%%  Static block size transforms
%%

I = imread('peppers.tif');
img = double(imread('peppers.tif'));
dim = size(img,2);

%figure;
%subplot(2,2,1), imshow(I);
%title(' Original ');

%step =[1, 2, 4, 8, 16, 32];                     %% Step Size
step = 16                     %% Step Size

%   *(1) = 16x16 blocks
%   *(2) = 8x8 blocks
%   *(3) = 4x4 blocks

M = [16, 8, 4]              %% Size of block
B_cnt = dim./M         %% # of blocks per image 
B_opt = size(M,2)      %% # of block sizes available

img_qdct = zeros(dim,dim,B_opt);
img_recon = zeros(dim,dim,B_opt);
img_final = zeros(dim,dim);

lambda = [ 0, 5, 10, 30, 62, 70, 90];
%lambda = 0
lambda_sz = size(lambda, 2);

J = zeros(B_cnt(3),B_cnt(3),B_opt+1);    %% Block Cost Fu nction
D = zeros(B_cnt(3),B_cnt(3),B_opt+1);    %% Block Distortion
R = zeros(B_cnt(3),B_cnt(3),B_opt+1);    %% Block Rate
JT = zeros(lambda_sz,B_opt);

for t=1:B_opt

   B = zeros(M(t),M(t));         % 16x16 block
   B_dct  = zeros(M(t),M(t));   % Block DCT

   for i=1:B_cnt(t)
      for j=1:B_cnt(t)
         row = M(t)*(i-1);
         col = M(t)*(j-1);

         % Create block
         B(:,:) = img(row+1:row+M(t),col+1:col+M(t));

         %% For block, find the DCT, distortion, bitrate, IDCT
         %%   and cost function
         [D(i,j,t), R(i,j,t), B_idct ] = blk_calc(B, step);
         img_recon(row+1:row+M(t),col+1:col+M(t),t) = B_idct;

      end % for j
   end % for i
end % for t


% find the total cost function, rate and distortion for the entire image
%   (4x4, 8x8, 16x16) 
for p=1:3
   DT(p) = sum(sum(D(:,:,p)))/(B_cnt(p)*B_cnt(p))
   RT(p) = sum(sum(R(:,:,p)))/(B_cnt(p)*B_cnt(p))
end

for w=1:lambda_sz
   J = D + lambda(w).*R;
   for p = 1:3
      JT(w,p) = sum(sum(J(:,:,p)))/(B_cnt(p)*B_cnt(p));
   end
end
JT

PSNR = 10*log10(255*255./DT);
%figure;
%plot(RT, PSNR);
%xlabel('Rate[bits/pixel]');
%ylabel('PSNR [dB] - 16x16, 8x8, 4x4, Best');

%figure;
%plot(lambda, JT);
%xlabel('Lambda ');
%ylabel('Lagrangian Cost Function');
%title('Stepsize = 1');

%figure;
%imagesc(img_recon(:,:,1));
%colormap=('gray');
%axis('equal');
%title(' peppers.tif, step = 32, block = 16');

%figure;
%imagesc(img_recon(:,:,2));
%colormap=('gray');
%axis('equal');
%title(' peppers.tif, step = 32, block = 8');

%figure;
%imagesc(img_recon(:,:,3));
%colormap=('gray');
%axis('equal');
%title(' peppers.tif, step = 32, block = 4');


% subplot(2,2,2), imagesc(img_recon(:,:,1));
% colormap=('gray');
% axis('equal');
% title(' step = 2, block = 16');

% subplot(2,2,3), imagesc(img_recon(:,:,2));
% colormap=('gray');
% axis('equal');
%title(' step = 2, block = 8');

%subplot(2,2,4), imagesc(img_recon(:,:,3));
%colormap=('gray');
%axis('equal');
%title(' step = 2, block = 4');