JPEG компресия

Компресия със загуби:

Теория

Пространствена област

Цветното изображението най-често се представя като съвкупност от три двумерни матрици. Всеки елемент от матрица е стойността за съответния пиксел на един от условно наречените основни цветове Червено, Зелено, Синьо. Тези цветове отговарят на спектралната чувствителност на окото.

Цветно (RGB) изображени представено като три измерна матрица

Визуализиране на отделните матрици (цветни канали)

Понякога при обработка на избраженията е по-удачно използването на други цветови схеми (YCrCb, HSV). В случая при компресиране се използва най-често YCrCb:

• Y - яркостен сигнал
• Cr, Cb - хроматични съставящи

${\begin{aligned}Y'&=&16&+&(65.481\cdot R'&+&128.553\cdot G'&+&24.966\cdot B')\\C_{B}&=&128&+&(-37.797\cdot R'&-&74.203\cdot G'&+&112.0\cdot B')\\C_{R}&=&128&+&(112.0\cdot R'&-&93.786\cdot G'&-&18.214\cdot B')\end{aligned}}$

A = double(rgb2ycbcr(A)); % 

Поради това, че окото е по чувствително към измененията на яркоста, хроматичните съставящи могат да се предават с по-малка резолюция и да се компресират с по-ниско качество без това да окаже видимо влияние върху качеството на избражението.(4:2:2,4:2:0...)

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');

Преминаване в честотна област

Изображението представено в пространствената област е със сравнително висока ентропия. (няма връзка между стойностите на съседни пиксели). С помощта на двумерна DCT трансформация изображението може да се изследва в честотната област.

При преминаване в честотната област изображението се разделя на блокове от 8х8 пиксела. Всеки блок се обработва поотделно.

N=8; 
for m = 1:N:Height
  for n = 1:N:Width
           
    t = y(m:m+N-1,n:n+N-1) - 128;
    Yy = dct2(t); 
    ...

изображнение в пространствена област (ляво) и в честотна област (дясно)

Коефициентите на DCT трансформацията са с висока сходимост, а повечето от тях са пренебрежимо малки. В замисимост от желаното качество се избират тегловни (квантуващи) коефициенти, на които се разделят DCT коефициентите.

След преминаване в честотната област лесно може да се премахне информацията за честотите, към които човешкото възприятие е по-малко чувствително.

• Високите честоти отговарят на резки преходи, контрастни линии, дребни обекти.
• Честоти с малки коефициенти (малка аплитуда) остават незабелижими в цялостното изображение

Квантуващата матрица е генирирана в съотвествие качеството, което желаем. При квантуването коефициентите с малка стойност се нулират нула. Квантуваща матрица за яркостния сигнал при качество 50%

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];

temp = floor(Yy./(Qscale*QstepsY) + 0.5);

Зиг-заг пренареждане на AC коефициентите

• DC - компонентата (постоянната съставна, първия елемент от трансформираната матрица) се кодира диференциално (за всички блокове от по 8х8), като разликите са Hoffman кодиране
• AC - коефициенти се пренареждат зигзаг образно, след това се прилага отстраняване на нулите, а накрая кодиране на хофман (или аритметична компресия), на ненулевите коефициенти, и статестическа компресия на нулевите.

референтно изображение

Лена

Matlab реализация

% 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 = 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(uint8(xqCb),[Height Width],'cubic');
c2 = imresize(uint8(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))

A

function y = ZigZag(x)
% y = ZigZag(x)
% returns the zigzag scan addresses of an 8 x 8 array
y(1) = x(1,1);y(2) = x(1,2);y(3) = x(2,1);y(4) = x(3,1);
y(5) = x(2,2);y(6) = x(1,3);y(7) = x(1,4);y(8) = x(2,3);
y(9) = x(3,2);y(10) = x(4,1);y(11) = x(5,1);y(12) = x(4,2);
y(13) = x(3,3);y(14) = x(2,4);y(15) = x(1,5);y(16) = x(1,6);
y(17) = x(2,5);y(18) = x(3,4);y(19) = x(4,3);y(20) = x(5,2);
y(21) = x(6,1);y(22) = x(7,1);y(23) = x(6,2);y(24) = x(5,3);
y(25) = x(4,4);y(26) = x(3,5);y(27) = x(2,6);y(28) = x(1,7);
y(29) = x(1,8);y(30) = x(2,7);y(31) = x(3,6);y(32) = x(4,5);
y(33) = x(5,4);y(34) = x(6,3);y(35) = x(7,2);y(36) = x(8,1);
y(37) = x(8,2);y(38) = x(7,3);y(39) = x(6,4);y(40) = x(5,5);
y(41) = x(4,6);y(42) = x(3,7);y(43) = x(2,8);y(44) = x(3,8);
y(45) = x(4,7);y(46) = x(5,6);y(47) = x(6,5);y(48) = x(7,4);
y(49) = x(8,3);y(50) = x(8,4);y(51) = x(7,5);y(52) = x(6,6);
y(53) = x(5,7);y(54) = x(4,8);y(55) = x(5,8);y(56) = x(6,7);
y(57) = x(7,6);y(58) = x(8,5);y(59) = x(8,6);y(60) = x(7,7);
y(61) = x(6,8);y(62) = x(7,8);y(63) = x(8,7);y(64) = x(8,8);

B

function Bits = jpgDCbits(dc,C)
% Bits = jpgDCbits(dc,C)
% Computes the exact number of bits to entropy code
% the DC differential coefficient using the JPEG default
% Huffman tables.
% Input:
% dc = predicted DC differential value
% C = ’Y’ for luma or ’C’ for chroma
% Output:
% # of bits for the DC differential
%
% CodeLengthY contains the lengths of codes for the
% luma DC differential for size categories fro 0 to 11 inclusive.
% Code length equals the length of the Huffman code for
% the size category plus the binary code for the value.
% CodeLengthC is the corresponding table for the chroma.

CodeLengthY = int16([3 4 5 5 7 8 10 12 14 16 18 20]);
CodeLengthC = int16([2 3 5 6 7 9 11 13 15 18 19 20]);
switch C
case 'Y'
    CodeLength = CodeLengthY;
case 'C'
    CodeLength = CodeLengthC;
end
if dc == 0
    Bits = double(CodeLength(1));
else
    Bits = double(CodeLength(round(log2(abs(dc))+0.5)+1));
end

C

function Bits = jpgACbits(x,C)
% y = jpgACbits(x,C)
% Computes the exact number of bits to Huffman code
% an N x N quantized DCT block
% using JPEG Huffman table for AC coefficients.
% DCT coefficients are zigzag scanned and
% the RL/Amplitude pairs are found. Then,
% using the JPEG RL/Amp table, number of bits to
% compress the block is determined.
% Input:
% x = N x N quantized DCT block
% C = ’Y’ for luma or ’C’ for chroma
% Outputs:
% Bits = total bits to encode the block
%
% RLCy is the JPEG run-length/size-category Huffman table
% for luma AC coefficients &
% RLCc is the corresponding Huffman table for the chroma component
% Each entry in the table is the code length for the RL/size
% category which equals the RL/size Huffman code plus the
% binary code of the actual non-zero coefficient amplitude.
% RLC{1}(1) is the EOB code corresponding to 0/0
% RLC{16}(1) is the Huffman code for RL/size = 16/0
%

RLCy = cell(1,15);
RLCy{1}=int16([4 3 4 6 8 10 12 14 18 25 26]);
RLCy{2} = int16([5 8 10 13 16 22 23 24 25 26]);
RLCy{3} = int16([6 10 13 20 21 22 23 24 25 26]);
RLCy{4} = int16([7 11 14 20 21 22 23 24 25 26]);
RLCy{5} = int16([7 12 19 20 21 22 23 24 25 26]);
RLCy{6} = int16([8 12 19 20 21 22 23 24 25 26]);
RLCy{7} = int16([8 13 19 20 21 22 23 24 25 26]);
RLCy{8} = int16([9 13 19 20 21 22 23 24 25 26]);
RLCy{9} = int16([9 17 19 20 21 22 23 24 25 26]);
RLCy{10} = int16([10 18 19 20 21 22 23 24 25 26]);
RLCy{11} = int16([10 18 19 20 21 22 23 24 25 26]);
RLCy{12} = int16([10 18 19 20 21 22 23 24 25 26]);
RLCy{13} = int16([11 18 19 20 21 22 23 24 25 26]);
RLCy{14} = int16([12 18 19 20 21 22 23 24 25 26]);
RLCy{15} = int16([13 18 19 20 21 22 23 24 25 26]);
RLCy{16} = int16([12 17 18 19 20 21 22 23 24 25 26]);
%
RLCc = cell(1,15);
RLCc{1}=int16([2 3 5 7 9 12 15 23 24 25 26]);
RLCc{2} = int16([5 8 11 14 20 22 23 24 25 26]);
RLCc{3} = int16([6 9 13 15 21 22 23 24 25 26]);
RLCc{4} = int16([6 9 12 15 21 22 23 24 25 26]);
RLCc{5} = int16([6 10 14 20 21 22 23 24 25 26]);
RLCc{6} = int16([7 11 17 20 21 22 23 24 25 26]);
RLCc{7} = int16([8 12 19 20 21 22 23 24 25 26]);
RLCc{8} = int16([7 12 19 20 21 22 23 24 25 26]);
RLCc{9} = int16([8 14 19 20 21 22 23 24 25 26]);
RLCc{10} = int16([9 14 19 20 21 22 23 24 25 26]);
RLCc{11} = int16([9 14 19 20 21 22 23 24 25 26]);
RLCc{12} = int16([9 14 19 20 21 22 23 24 25 26]);
RLCc{13} = int16([9 18 19 20 21 22 23 24 25 26]);
RLCc{14} = int16([11 18 19 20 21 22 23 24 25 26]);
RLCc{15} = int16([13 18 19 20 21 22 23 24 25 26]);
RLCc{16} = int16([11 17 18 19 20 21 22 23 24 25 26]);
%
switch C
case 'Y'
    RLC = RLCy;
case 'C'
    RLC = RLCc;
end
x1 = ZigZag(x); % zigzag scan the 8 x 8 DCT blocks
%
k = 2; Count = 0; Bits = double(0);
while k <= 64
    if x1(k) == 0
        Count = Count + 1;
        if k == 64
            Bits = Bits + double(RLC{1}(1));
            break;
        end
    else
        if Count == 0
            RL = Count;
            Level = round(log2(abs(x1(k)))+0.5);
            Bits = Bits + double(RLC{RL+1}(Level+1));
        elseif Count >= 1 && Count <=15
            RL = Count;
            Level = round(log2(abs(x1(k)))+0.5);
            Bits = Bits + double(RLC{RL+1}(Level));
            Count = 0;
        else
            Bits = Bits + double(RLC{16}(1));
            Count = Count - 16;
        end
    end
    k = k + 1;
end

Литература

K. S. Thyagarajan, Still Image and Video Compression with MATLAB

DCT - discrete cosine transformation - дискретна косиносова трансформация