Homework 2

function O = blockdct(I,N)
% O = BLOCKDCT(I,N) does the 2D DCT on BxB blocks of I.
% With the image processing toolbox, i could use blkproc, as
% in O = blkproc(I,[N N], @dct2), but I don't have the toolbox at home.
% only handles square blocks.

% Adam Berenzweig, February 28, 2001

if ~isa(I,'double')
   I = double(I);
end

[mi,ni] = size(I);
mb = floor(mi/N);	% how many blocks fit down columns
nb = floor(ni/N); % how many blocks fit across rows

O = zeros(size(I));

% define zonal mask
zone = zeros(N,N);
%zone(1:3,1:3) = [1 1 1; 1 1 0; 1 0 0]; % this only keeps 6 coefficients... how do you do 8?
zone(1:4,1:4) = [1 1 1 1; 1 1 0 0; 1 0 0 0; 1 0 0 0];	% i think this is the right 1:8 filter.
%zone = hankel([1 1 1 1 0 0 0 0]);	% keep 10 coefficients
%zone = ones(8,8);	% for testing. should get the original image back.
%zone(1:3,1:3) = [1 1 1; 1 0 0; 0 0 0];	% 1:16 zone filter

for i = 0:(mb-1)
   for j = 0:(nb-1)
      block = I( (N*i+1):(N*(i+1)) , (N*j+1):(N*(j+1)) );
      block = dct(block);	% dct on the columns
      block = dct(block')';	% dct on the rows
      
      % mask it here with zone
      block = block .* zone;
      
      % now do idct to get image back
      block = idct(block);	% down the columns
      block = idct(block')'; % across the rows;
      
      O( (N*i+1):(N*(i+1)) , (N*j+1):(N*(j+1)) ) = block;  % put it back
   end
end
% sheesh.. wasted four hours trying to vectorize this.  who says life is
% too short for writing loops?

% calculate the MSE
% (actually, it's the average least-squares error estimate of MSE
MSE = sum(sum((I-O).^2)) / (mi*ni)

% calculate the SNR
mI = mean(mean(I)); 
vI = sum(sum((I-mI).^2)) / (mi*ni);
SNR = 10*log10(vI/MSE)

% write out to a file
out = uint8(O);
imwrite(out,'output.png','png');