Difference between revisions of "JPEG компресия"

From Ilianko
Line 8: Line 8:
 
*Cr, Cb - хроматични съставящи
 
*Cr, Cb - хроматични съставящи
  
−
<math>
+
<math>
 
\begin{align}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{align}
 
\begin{align}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{align}
 
</math>
 
</math>
 +
 
Поради това, че окото е по чувствително към измененията на яркоста, хроматичните съставящи могат да се предават с по-малка резолюция и да се компресират с по-ниско качество без това да окаже видимо влияние върху качеството на избражението.
 
Поради това, че окото е по чувствително към измененията на яркоста, хроматичните съставящи могат да се предават с по-малка резолюция и да се компресират с по-ниско качество без това да окаже видимо влияние върху качеството на избражението.
  

Revision as of 11:16, 30 June 2014

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

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

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

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

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

Матрица пространствена област - space domain 

|10|20|30|255...
|20|10|35..
|255...
|...

малка връзка между стойността на всеки пиксел

Всеки цвят се обработва по отделно. Възможно е и пъвро да се преобразува в друга друга светова област. (HSV). Зеления/Яркостния сигнал могат да се обработват в нормалната си резолюция, а останалите компонетни се децемират два пъти (4:2:2,4:2:1...)

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

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

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

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

3. Зиг-заг пренареждане

  • 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 = 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))


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

Retrieved from "https://ilianko.com/index.php?title=JPEG_компресия&oldid=5176"