% Project 1
% Andy Doran and Brent Hendricks
%
% function pict = dcomp2(cname,name,xdim)
%
% cname = 'compressed image file' (as in 'bitsLib')
% name = 'input image' (as in 'lenna.256')
% xdim = x dimension of input image (usually 256 or 512)
%
% pict is the decompressed image.  We load in the input image (name)
% in order to compare against decompressed image.
%
% This program loads in the compressed bitstream, parses of the header,
% and decompresses the image based on a known bitmask.
%

function pict = dcomp2(cname,name,xdim)

temp = 16;
qpict = zeros(xdim,xdim);
unquant = zeros(xdim,xdim);
Fmax = zeros(xdim/temp,xdim/temp);
Fmin = zeros(xdim/temp,xdim/temp);
xmin = 32;
dcmp = [1 .75 .5 .375 .25 .125 .0625];

%
% Read in from file
%
fid2 = fopen(cname,'rb');
scale=fread(fid2,1,'double');
maskNum=fread(fid2,1,'uchar');
[msk,dcomp] = mask(dcmp(maskNum));
for i = 1:(xdim/temp)    % 1:16
  for j = 1:(xdim/temp)  % 1:16
    if (msk(i,j) ~= 0)
      i1 = (i-1)*temp+1;
      i2 = i*temp;
      j1 = (j-1)*temp+1;
      j2 = j*temp;
      Fmin(i,j) = bread(fid2, -8, 1 );
      Fmax(i,j) = bread(fid2, -8, 1 );
      qpict(i1:i2,j1:j2)=reshape(bread(fid2, msk(i,j),(xdim/temp)*(xdim/temp) ), xdim/temp,xdim/temp);
   end
  end
end
fclose(fid2);

%
% Reverse scale
%
for i = 1:(xdim/temp)    % 1:16
  for j = 1:(xdim/temp)  % 1:16
    if (msk(i,j) ~= 0)
      i1 = (i-1)*temp+1;
      i2 = i*temp;
      j1 = (j-1)*temp+1;
      j2 = j*temp;
      unquant(i1:i2,j1:j2)=dquant(qpict(i1:i2,j1:j2),Fmax(i,j),Fmin(i,j),msk(i,j));
    end
  end
end

pict = unquant*scale/127;
size(pict);

% Plot quantized transformed image
figure(3);
image(pict);
colormap(gray(xdim));
txt = [num2str(name) ' after Haar wavelet transform and Quantization'];
title(txt)

% Get image range
imax = round(log2(xdim));
imin = round(log2(xmin));

% Call inverse transform application
for i = imin:imax
  b = 2^i;                % Get number of pixels to work on
  pict(1:b,1:b) = itrans(pict(1:b,1:b)')';  % Work on cols of image
  pict(1:b,1:b) = itrans(pict(1:b,1:b));    % Work on rows of image
end

% Plot decoded image
figure(4)
image(pict);
colormap(gray(xdim));
txt = [num2str(name) ' after Inverse Haar wavelet transform'];
title(txt)

% Get input comparison image
in = freadbin(name,xdim,xdim);
a = in(:);  % Get column vector of initial image
b = pict(:);  % Get column vector of final image

% Compute Mean Squared Error
MSE = sum([a-b].^2)/xdim^2
