JPEG компресия
From Ilianko
% 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))
C
function 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
