x = imread('DHSquare.jpg'); 
% Original image is a multiple of 8 
x=double(x);
r=x(:,:,1); g=x(:,:,2); b=x(:,:,3);
xgray=0.2126*r+0.7152*g+0.0722*b;
n=1; 
figure(3);imagesc(uint8(xgray));colormap(gray)
title( 'Original');
sizeImage=size(xgray);
Xf=zeros(sizeImage(1));
for p = [1 2 4 10 16 32];
    for k=0:(sizeImage(1)/8)-1
        for l=0:(sizeImage(2)/8)-1
        c=xgray(8*k+1:8*k+8,8*l+1:8*l+8);
        c=dctV2(p(n),c);
        Xf(8*k+1:8*k+8,8*l+1:8*l+8)=c;    
        end
    end
figure(p(n));imagesc(uint8(Xf));colormap(gray)
title('Compressed p = 1' );
end


function Z=dctV2(p,X) 
[m,n]=size(X);
for i=2:n
    for j=1:n
        C(i,j)=sqrt(2/n)*cos((i-1)*(j-1/2)*pi/n);
    end
end
C(1,1:n)=ones(1,n)/sqrt(n);
Q = p*8./hilb(8);
Xc=X-128;
Y=C*Xc*C';
Yq =round(Y./Q);
Ydq = Yq.*Q;
Xdq = C'*Ydq*C;
Z = Xdq+128;
end



