JPEG компресия

% JPEG Основен тип компресия

clear;
A = imread('birds.ras'); 
[Height,Width,Depth] = size(A);
N = 8; % Размер на подматрици

% Размерите трябва да са кратни на 8
if mod(Height,N) ~= 0
    Height = floor(Height/N)*N;
end
if mod(Width,N) ~= 0
    Width = floor(Width/N)*N;
end

A1 = A(1:Height,1:Width,:);
clear A
A = A1;

SamplingFormat = '4:2:0';
if Depth == 1 % черно бяло
    y = double(A);
else
    A = double(rgb2ycbcr(A)); % Цветно
    y = A(:,:,1);
    switch SamplingFormat
    case '4:2:0'
        Cb = imresize(A(:,:,2),[Height/2 Width/2],'cubic');
        Cr = imresize(A(:,:,3),[Height/2 Width/2],'cubic');
    case '4:2:2'
        Cb = imresize(A(:,:,2),[Height Width/2],'cubic');
        Cr = imresize(A(:,:,3),[Height Width/2],'cubic');
    end
end

jpgQstepsY = [  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];

QstepsY = jpgQstepsY;
Qscale = 1.5;
Yy = zeros(N,N);
xqY = zeros(Height,Width);
acBitsY = 0;
dcBitsY = 0;

if Depth > 1
    jpgQstepsC = [  17 18 24 47 66 99 99 99;...
                    18 21 26 66 99 99 99 99;...
                    24 26 56 99 99 99 99 99;...
                    47 66 99 99 99 99 99 99;...
                    99 99 99 99 99 99 99 99;...
                    99 99 99 99 99 99 99 99;...
                    99 99 99 99 99 99 99 99;...
                    99 99 99 99 99 99 99 99];
    QstepsC = jpgQstepsC;
    YCb = zeros(N,N);
    YCr = zeros(N,N);
    switch SamplingFormat
    case '4:2:0'
        xqCb = zeros(Height/2,Width/2);
        xqCr = zeros(Height/2,Width/2);
    case '4:2:2'
        xqCb = zeros(Height,Width/2);
        xqCr = zeros(Height,Width/2);
    end
    
    acBitsCb = 0;
    dcBitsCb = 0;
    acBitsCr = 0;
    dcBitsCr = 0;
end

% Изчисляване на y  компонента
for m = 1:N:Height
    for n = 1:N:Width
            
        t = y(m:m+N-1,n:n+N-1) - 128;
        Yy = dct2(t); % N x N 2D DCT of input image
            % quantize the DCT coefficients
        temp = floor(Yy./(Qscale*QstepsY) + 0.5);

        % Calculate bits for the DC difference    
        if n == 1
            DC = temp(1,1);
            dcBitsY = dcBitsY + jpgDCbits(DC,'Y');
        else
                DC = temp(1,1) - DC;
                dcBitsY = dcBitsY + jpgDCbits(DC,’Y’);
                DC = temp(1,1);
        end
            
            % Calculate the bits for the AC coefficients
        ACblkBits = jpgACbits(temp,'Y');
        acBitsY = acBitsY + ACblkBits;
        % dequantize & IDCT the DCT coefficients
        xqY(m:m+N-1,n:n+N-1)= idct2(temp .* (Qscale*QstepsY))+ 128;
    end
end
% If the input image is a color image,
% calculate the bits for the chroma components
if Depth > 1
    if strcmpi(SamplingFormat,'4:2:0')
        EndRow = Height/2;
    else
    EndRow = Height;
    end
    
    for m = 1:N:EndRow
        for n = 1:N:Width/2
            t1 = Cb(m:m+N-1,n:n+N-1) - 128;
            t2 = Cr(m:m+N-1,n:n+N-1) - 128;
            Ycb = dct2(t1); % N x N 2D DCT of Cb image
            Ycr = dct2(t2);
            temp1 = floor(Ycb./(Qscale*QstepsC) + 0.5);
            temp2 = floor(Ycr./(Qscale*QstepsC) + 0.5);
            if n == 1
                DC1 = temp1(1,1);
                DC2 = temp2(1,1);
                dcBitsCb = dcBitsCb + jpgDCbits(DC1,'C');
                dcBitsCr = dcBitsCr + jpgDCbits(DC2,'C');
            else
                DC1 = temp1(1,1) - DC1;
                DC2 = temp2(1,1) - DC2;
                dcBitsCb = dcBitsCb + jpgDCbits(DC1,'C');
                dcBitsCr = dcBitsCr + jpgDCbits(DC2,'C');
                DC1 = temp1(1,1);
                DC2 = temp2(1,1);
            end
            ACblkBits1 = jpgACbits(temp1,'C');
            ACblkBits2 = jpgACbits(temp2,'C');
            acBitsCb = acBitsCb + ACblkBits1;
            acBitsCr = acBitsCr + ACblkBits2;

            % dequantize and IDCT the coefficients
            xqCb(m:m+N-1,n:n+N-1)= idct2(temp1 .* (Qscale*QstepsC))+ 128;
            xqCr(m:m+N-1,n:n+N-1)= idct2(temp2 .* (Qscale*QstepsC))+ 128;
        end
    end
end
%
mse = std2(y-xqY);
snr = 20*log10(std2(y)/mse);
sprintf('SNR = %4.2f\n',snr)

if Depth == 1
    TotalBits = acBitsY + dcBitsY;
    figure,imshow(xqY,[])
    title(['JPG compressed ' '@ ' num2str(TotalBits/(Height*Width)) ' bpp'])
else
    
TotalBits = acBitsY + dcBitsY + dcBitsCb +...
acBitsCb + dcBitsCr + acBitsCr;
c1 = imresize(xqCb,[Height Width],'cubic');
c2 = imresize(xqCr,[Height Width],'cubic');
mseb = std2(A(:,:,2)-c1);
snrb = 20*log10(std2(A(:,:,2))/mseb);
msec = std2(A(:,:,3)-c2);
snrc = 20*log10(std2(A(:,:,3))/msec);
sprintf('SNR(Cb) = %4.2fdB\tSNR(Cr) = %4.2fdB\n',snrb,snrc)
xq(:,:,1) = xqY;
xq(:,:,2) = c1;
xq(:,:,3) = c2;
figure,imshow(ycbcr2rgb(uint8(round(xq))))
title(['JPG compressed ' '@ ' num2str(TotalBits/(Height*Width)) ' bpp'])
end
sprintf('Bit rate = %4.2f bpp\n',TotalBits/(Height*Width))