%A=imread('luxtest.png');
%[Agrey]=rgb2gray(A);
Agrey = imread('chicken-s.png');
Sorig = size(Agrey);
figure, imshow(Agrey);
P = [8 8]; % side length of DCT cell

% Pad size of image A to multiple of P
s = ceil(size(Agrey)./P).*P;
Agrey(s(1),s(2)) = 0;
Agrey=double(Agrey); % convert to double for DCT later
A2 = zeros(s);
X = zeros(s);
% Quantisation matrix
Q = ones(P); % no quantisation apart from rounding

QJPEG50 = [16 11 10 16 24  40  51  61;
           12 12 14 19 26  58  60  55;
           14 13 16 24 40  57  69  56;
           14 17 22 29 51  87  80  62;
           18 22 37 56 68  109 103 77;
           24 35 55 64 81  104 113 92;
           49 64 78 87 103 121 120 101;
           72 92 95 98 112 100 103 99];
QJPEG90 = clip(round(QJPEG50 * (100-90)/50), 1, 255);
QJPEG10 = clip(round(QJPEG50 * 50 / 10), 1, 255);

Q3 = round((1+log10(1:P(1))' * (1:P(2)))); % stronglyremove high-freq
Q4 = round((1+log10(1+log10(1:P(1))' * (1:P(2))))); % strongly remove high-freq
Q5 = round((linspace(1,P(1)*10, P(1)))' * linspace(1,1,P(2))); % strongly remove one component
Q6 = round((linspace(3,15,P(1))' * linspace(3,15,P(2)))); % remove high-freq, linear.
Q7 = round((linspace(15,3,P(1))' * linspace(15,3,P(2)))); % remove low-freq, linear
Q8 = round((linspace(10,1,P(1))' * linspace(10,1,P(2)))); % remove low-freq, linear
%Q4 = logspace(0,2,P(2))' * logspace(0,2,P(1)); % strongly remove high-freq
%Q5 = logspace(0,1,P)' * logspace(0,1,P); % weakly remove high-freq
%Q6 = linspace(P/4,1,P)' * linspace(P/4,1,P); % weakly remove low-freq
%Q7 = logspace(1,0,P)' * logspace(1,0,P); % weakly remove low-freq
%Q8 = logspace(2,0,P)' * logspace(2,0,P); % strongly remove low-freq
%Q9 = ones(P,P); Q9(1,1)=1000; % strongly remove lowest freq component
%Q10 = ones(P,P); Q10(2:P)=1000; % strongly remove one row of high-freq component
%Q11 = ones(P,P); Q11(1,:)=1000; % strongly remove other row of high-freq component
%Q12 = ones(P,P); Q12(:,P)=1000; % strongly remove one high-freq component
%Q13 = ones(P,P); Q13(1,1)=0;
%Q14 = zeros(P,P); Q14(1,1)=1;
%Q15 = ones(P,P); Q15(1:P)=1000; % strongly remove one row of high-freq component
%Q16 = logspace(0,.5,P)' * logspace(0,.5,P);
%Q = QJPEG50
%Q = Q4
Q=QJPEG50

% iterate through 8x8 blocks
Agrey = Agrey - 128; % Makes the DCT work better
for i = (0:floor(s(1)/P(1) -1))
    for j = (0:floor(s(2)/P(2) -1))
        irange = ( i*P(1):(i+1)*P(1) - 1 ) +1;
        jrange = ( j*P(2):(j+1)*P(2) - 1 ) +1;
        Ap = Agrey(irange,jrange);
        % Transform
        Xp = DCT(Ap);
        % Quantise
        Xp = round(Xp./Q);
        % Save quantised result in X matrix
        X(irange,jrange) = Xp;
        % Now decompress so we can compare....
        % Reverse effect of Q division (mostly)
        Xp = Xp .* Q;
        % Inverse transform into A2 matrix
        A2(irange,jrange) = DCTinv(Xp);
    end
end
Agrey = Agrey + 128; %Undo earlier subtraction
A2 = A2 + 128; %Undo earlier subtraction
Ximg = abs(X);
Ximg = log(Ximg+1-min(min(Ximg)));
Ximg2 = Ximg*255/max(max(Ximg));
figure, imshow(uint8(Ximg2));
Agrey = Agrey(1:Sorig(1),1:Sorig(2));
A2 = A2(1:Sorig(1),1:Sorig(2));
%figure, imshow(uint8(X*255/max(max(X))));
figure, imshow(uint8(A2*max(max(A2))/255));
% calculate fraction of zeros
zerofrac = (s(1)*s(2) - length(find(X))) /  (s(1)*s(2))
entropy = Entropy(X)
figure, imshow(uint8(abs(A2 - Agrey)));



        