jpg.m

jpeg压缩matlab源代码。由均匀量化和huffman编码组成

% This function does a JPEG compression on a grayscale or RGB image
% Usage: [msqe,compr]=jpg(file,scale)
% 'file' must contain an m*n*3 (color) or m*n (grayscale) matrix with
%        RGB values in [0,1]
% 'scale' is the JPG quality factor, higher is better.
% 'msqe' is the mean squared error after compression
% 'compr' is the compression ratio in percentages of the original file
%
% Note that the Huffman compression is not JPEG standard (see PDF)
% Huffman code by Karl Skretting is used.
%
%%%%%%%
% Changes 11-11-2003
% 1. Changed 'true' and '' to 1 and 0. This was undefined in older
% versions of matlab.
% 2. Changed 'if color' to 'if color==1', needed for older matlab versions.
% 3. 21-4-2005 Fixed compression issue. Used 'hufflen' in a wrong way. Now uses a 
% call to 'huff06' to calculate bitlengths. 
%%%%%%
% -Arno Swart, [e-mail]
%%%%%

function [msqe,compr]=jpg(file,scale)

% Add path to Huffman coder.
% Written by Karl Skretting (see .txt)
addpath('./huffman');

% Load file, get statistics
% File contains an m * n * 3 matrix 'im' (color) or m * n matrix (gray)
load(file);
original=im;
width = size(im,2);
height = size(im,1);
xblocks = width/8;
yblocks = height/8;

if length(size(im))==3
    color=1;
    numbits=width*height*8*3;
else
    color=0;
    numbits=width*height*8;
end

% Display original picture.
figure(1);cla;
if color==1
    image(im);
else
    % For grayscale images: scale and give the colormap
    imagesc(im);colormap gray;
end
title(sprintf('Original %s picture.',file));


lumtable  =  lumquant(scale);
if color==1
    chrtable  =  chrquant(scale);
    im = RGB2YUV(im,width,height);   
    
    % Show YUV components
    figure(2);cla
    subplot(2,2,1); imagesc(im(:,:,1));colormap gray;colorbar;
    title('Y component');
    subplot(2,2,2); imagesc(im(:,:,2));colormap gray;colorbar;
    title('U component');
    subplot(2,2,3); imagesc(im(:,:,3));colormap gray;colorbar;
    title('V component');
    
    % Now for the encoding we need values 0-255;
    im = round(im*255);
    % Unfortunately values can get out of range, clip them.
    im(im>255)=255;
    im(im<0)=0;
    
    % Y component uses luminance quantisation
    dctcoefY = jpgtrans(im(:,:,1),width,height,xblocks,yblocks,lumtable);
    % U component uses chrominance quantisation
    dctcoefU = jpgtrans(im(:,:,2),width,height,xblocks,yblocks,chrtable);
    % V component uses chrominance quantisation
    dctcoefV = jpgtrans(im(:,:,3),width,height,xblocks,yblocks,chrtable);
else
    % Now for the encoding we need values 0-255;
    im = round(im*255);
    
    % Y component uses luminance quantisation
    dctcoef = jpgtrans(im,width,height,xblocks,yblocks,lumtable);
end


%
% Now we calculate how many bits a Huffman coding would require
%
if color==1
    bitcount=huffcount(dctcoefY,xblocks,yblocks);
    bitcount=bitcount+huffcount(dctcoefU,xblocks,yblocks);
    bitcount=bitcount+huffcount(dctcoefV,xblocks,yblocks);
else
    bitcount=huffcount(dctcoef,xblocks,yblocks);
end

disp(sprintf('Total bitcount %d',bitcount));
compr = 100-100*bitcount/numbits;
disp(sprintf('Compression: %f',compr));

if color==1
    % Y component uses luminance quantisation
    imgY = jpginvtrans(dctcoefY,width,height,xblocks,yblocks,lumtable);
    % U component uses chrominance quantisation
    imgU = jpginvtrans(dctcoefU,width,height,xblocks,yblocks,chrtable);
    % V component uses chrominance quantisation
    imgV = jpginvtrans(dctcoefV,width,height,xblocks,yblocks,chrtable);
    
    % Make a new image
    clear im;
    im(:,:,1)=imgY/255;
    im(:,:,2)=imgU/255;
    im(:,:,3)=imgV/255;
       
    im = YUV2RGB(im,width,height);  
    
    % Unfortunately values can get out of range, clip them.
    im(im>1)=1;
    im(im<0)=0;
    
    figure(3);cla
    image(im);
    title('Reconstructed image');
    msqe = sum(sum(sum((original-im).^2)))/(3*width*height);
else
    % Y component uses luminance quantisation
    im = jpginvtrans(dctcoef,width,height,xblocks,yblocks,lumtable);
    im=im/255;
    figure(3);cla
    imagesc(im);colormap gray
    title('Reconstructed image');
    msqe = sum(sum((original-im).^2))/(width*height);
end


disp(sprintf('Mean Squared Error %f', msqe));

function F = dct(f)

c=ones(1,8);
c(1) = 1./sqrt(2);

for u=0:7
    for v=0:7
        F(u+1,v+1) = c(u+1)*c(v+1)* sum(sum(f.*(cos((2*[0:7]'+1)*u*pi/16)*cos((2*[0:7]+1)*v*pi/16))))/4;
    end
end


function f = invdct(F)

c=ones(1,8);
c(1) = 1./sqrt(2);

for x=0:7
    for y=0:7
        f(x+1,y+1) = sum(sum(F.*(c'*c).*(cos((2*x+1)*[0:7]'*pi/16)*cos((2*y+1)*[0:7]*pi/16))))/4;
    end
end

%
% Standard JPG luminance quantisation table
%
function qtable  =  lumquant(scale)

    qtable=[    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];
    qtable = round(qtable./scale);
    qtable(qtable<1) = 1;
    qtable(qtable>255) = 255;
    
%
% Standard JPG chrominance quantisation table
%
function qtable  =  chrquant(scale)
    qtable = [  17, 18, 24, 47, 99, 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];
    
    qtable = round(qtable./scale);
    qtable(qtable<1) = 1;
    qtable(qtable>255) = 255;    

%
% Do zigzag ordering on matrix x, return vector x
%
function x = zigzag(x)

    zigzag= [   0,  1,  5,  6, 14, 15, 27, 28,
                2,  4,  7, 13, 16, 26, 29, 42,
                3,  8, 12, 17, 25, 30, 41, 43,
                9, 11, 18, 24, 31, 40, 44, 53,
                10, 19, 23, 32, 39, 45, 52, 54,
                20, 22, 33, 38, 46, 51, 55, 60,
                21, 34, 37, 47, 50, 56, 59, 61,
                35, 36, 48, 49, 57, 58, 62, 63]+1;

    % convert matrices to column
    zigzag = zigzag(:);
    x = x(:);

    % do zigzag ordering
    x(zigzag) = x;

%
% Do 'inverse' zigzag ordering on vector x, return matrix x
%
function x = dezigzag(x)

    zigzag= [   0,  1,  5,  6, 14, 15, 27, 28,
                2,  4,  7, 13, 16, 26, 29, 42,
                3,  8, 12, 17, 25, 30, 41, 43,
                9, 11, 18, 24, 31, 40, 44, 53,
                10, 19, 23, 32, 39, 45, 52, 54,
                20, 22, 33, 38, 46, 51, 55, 60,
                21, 34, 37, 47, 50, 56, 59, 61,
                35, 36, 48, 49, 57, 58, 62, 63]+1;
    % work on columns
    zigzag = zigzag(:);
    x = reshape(x(zigzag),8,8);

%    
% Do a zero RLE on acoef, return arrays zerocounts and nonzeros.
%
function [zerocounts,nonzeros]=zerorle(acoef);   

	% Easiest is to first determine the nr of zeros at the end.
	k=31;
	while acoef(k)==0
        k=k-1;
        if k==0;
            zerocounts(1)=0;
            nonzeros(1)=0;
            break;
        end
	end
	
	curzerocount=0;
	l=1;
	for i=1:k
        if(acoef(i)==0)
            curzerocount=curzerocount+1;
            % Exception: 16 zeros
            if curzerocount==16;
                zerocounts(l)=15;
                nonzeros(l)=0;
                curzerocount=0;
                l=l+1;
            end
        else
            zerocounts(l)=curzerocount;
            nonzeros(l)=acoef(i);
            l=l+1;
            curzerocount=0;
        end
	end
	zerocounts(l)=0;
	nonzeros(l)=0;

%
% Convert from YUV to RGB color space.
%
function x=YUV2RGB(x,width,height)
    rgbmat=[1,0,1.402; 1,-0.34414,-0.71414;1,1.772,0];
       
    % apply transformation to all RGB tripples
    for i=1:width
        for j=1:height
            x(j,i,:) = rgbmat*(squeeze(x(j,i,:)) -[0;0.5;0.5]);
        end
    end
    
     
%
% Convert from RGB to YUV color space.
%
function x=RGB2YUV(x,width,height)
    yuvmat = [  0.299   ,0.587      ,0.114; 
                -0.1687 ,-0.3313    ,0.5; 
                0.5     ,-0.4187    ,-0.0813];

    % apply transformation to all RGB tripples
    for i=1:width
        for j=1:height
            x(j,i,:)=yuvmat*(squeeze(x(j,i,:)))+[0;0.5;0.5];
        end
    end

%    
% JPEG procedure for a color component
%
function dctcoef = jpgtrans(im,xsize,ysize,xblocks,yblocks,qtable)
	prevdc=0;
	%JPG scan order is columns first
	for j=1:xblocks
        for i=1:yblocks
            % Extract 8x8 block
            imblock{i,j} = im((i-1)*8+1:(i*8),(j-1)*8+1:(j*8));
            % Shift 128 down
            imblock{i,j} = imblock{i,j} - 128;
            % Do a DCT
            dctcoef{i,j} = dct(imblock{i,j});
            % Quantize
            dctcoef{i,j} = round(dctcoef{i,j}./qtable);
            
            % Differential code the DC
            temp = dctcoef{i,j}(1,1);
            dctcoef{i,j}(1,1) = dctcoef{i,j}(1,1) - prevdc;      
            prevdc=temp;  
            
            % Zigzag
            dctcoef{i,j} = zigzag(dctcoef{i,j});
        end
	end
    
%
% Inverse JPEG procedure
%
function img = jpginvtrans(dctcoef,xsize,ysize,xblocks,yblocks,qtable)   
	prevdc=0;
	for j=1:yblocks
        for i=1:xblocks
            % De-zigzag
            dctcoef{i,j} = dezigzag(dctcoef{i,j});
            
            %Un-differential code the DC
            dctcoef{i,j}(1,1) = dctcoef{i,j}(1,1) + prevdc;   
            prevdc=dctcoef{i,j}(1,1);
            
            % De-Quantize
            dctcoef{i,j} = dctcoef{i,j}.*qtable;
            
            % Do an inverse DCT
            imblock{i,j} = invdct(dctcoef{i,j});
            
            % Shift 128 up
            imblock{i,j} = imblock{i,j} + 128;   
            
            %Rebuild the image
            img((i-1)*8+1:i*8,(j-1)*8+1:j*8) = imblock{i,j};
        end
	end
    
%
% Count the nr. bits that a Huffman coding would take.
%
function bitcount=huffcount(dctcoef,xblocks,yblocks)
	l=1;
	bitcount=0;
	for j=1:yblocks
        for i=1:xblocks
            % The dccoef's are treated separately.
            dccoef(l)=dctcoef{i,j}(1);
            l=l+1;
            % Get the ac coef's
            accoef = dctcoef{i,j}(2:32);
	
            [zerocounts,nonzeros]=zerorle(accoef);  
            [dummy,res]=huff06({zerocounts,nonzeros},8,1);
            bitcount = bitcount + res(3,3);
        end
    end
    [dummy,res]=huff06({dccoef},8,1);
    bitcount = bitcount+ res(2,3);