% Andrew Puryear
% 362 Project - compress_jpega
%   Matlab file takes images and returns COLA-JPEG files. Allows for variable compression.

function [ imw, imsize] = compress_jpega( X , halfMax)

% discreatize to 256 values and normalize
im = round(X/max(X(:))*256);
hopsize=8;
xpixel=size(im,2);
ypixel=size(im,1);

% create COLA window
wind=[ones(1,10)*.25;
.25, .25, ones(1,6)*3/4, .25, .25;
.25, .75, ones(1,6), .75, .25;
.25, .75, ones(1,6), .75, .25;
.25, .75, ones(1,6), .75, .25;
.25, .75, ones(1,6), .75, .25;
.25, .75, ones(1,6), .75, .25;
.25, .75, ones(1,6), .75, .25;
.25, .25, ones(1,6)*3/4, .25, .25;
ones(1,10)*.25];

% decompose image.

% top left corner
imw = [zeros(10,2),[zeros(2,8);im(1:8,1:8).*wind(3:10,3:10)]];

% top edge
m=1;
for k=7:hopsize:xpixel
    if k+9>xpixel
        imw(1:10,m*10+1:m*10+10)=[zeros(2,10);[im(1:8,k:xpixel).*wind(3:10,1:xpixel-k+1),zeros(8,9-xpixel+k)]];
    else
        imw(1:10,m*10+1:m*10+10)=[zeros(2,10);im(1:8,k:k+9).*wind(3:10,1:10)];
    end
    m=m+1;
end

% left edge
m=1;
for k=7:hopsize:ypixel
    if k+9>ypixel
        imw(m*10+1:m*10+10,1:10)=[zeros(10,2),[im(k:ypixel,1:8).*wind(1:ypixel-k+1,3:10);zeros(9-xpixel+k,8)]];
    else
        imw(m*10+1:m*10+10,1:10)=[zeros(10,2),im(k:k+9,1:8).*wind(1:10,3:10)];
    end
    m=m+1;
end

% middle
n=1;
for k=7:hopsize:xpixel
    m=1;
    for p=7:hopsize:ypixel
        if k+9>xpixel & p+9>ypixel           
            imw(m*10+1:m*10+10,n*10+1:n*10+10)=[[im(k:xpixel,p:ypixel).*wind(1:xpixel-k+1,1:ypixel-p+1),zeros(xpixel-k+1,9-ypixel+p)];zeros(9-xpixel+k,10)];
        else
            if k+9>xpixel
                imw(m*10+1:m*10+10,n*10+1:n*10+10)=[im(p:p+9,k:xpixel).*wind(1:10,1:xpixel-k+1),zeros(10,9-xpixel+k)];
            elseif p+9>ypixel
                imw(m*10+1:m*10+10,n*10+1:n*10+10)=[im(p:ypixel,k:k+9).*wind(1:ypixel-p+1,1:10);zeros(9-ypixel+p,10)];
            else    
                imw(m*10+1:m*10+10,n*10+1:n*10+10)=im(p:p+9,k:k+9).*wind;
            end
        end
        m=m+1;
    end
    n=n+1;
end
imsize=[size(im,1),size(im,2)];
imw=round(imw);

% Make Q table
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nBits = zeros(8,8);
for i=1:8
 for j=1:8
%  nBits(i,j) = 8 - (i^2 + j^2)^0.5;
   d = i^2 + j^2;
  nBits(i,j) = 8 -  ( (8*d) / (d + halfMax));
 end
end
ma = max(nBits(:));
mn = min(nBits(:));
nBits = round( 8* (nBits - mn) ./ (ma - mn));
compressionFactor = 512/sum(sum(nBits))
coefScaleFactor = 2 .^ (nBits - 8);
q = [[coefScaleFactor;coefScaleFactor(size(coefScaleFactor,1),:);coefScaleFactor(size(coefScaleFactor,1),:)],[coefScaleFactor(:,size(coefScaleFactor,2));0.0039;0.0039],[coefScaleFactor(:,size(coefScaleFactor,2));0.0039;0.0039]];



for i=1:10:size(imw,1)
    for j=1:10:size(imw,2)
        dctCoef(i:i+9,j:j+9) = round(dct2(imw(i:i+9,j:j+9)) .* q) ./ q;
    end
end
imw=dctCoef;