in dev matlab guide ~ read.
Image Compression in MATLAB

Image Compression in MATLAB

An implementation of a lossy image compression format (GPJ) in MATLAB that is similar to JPEG. This can be further extended to provide interlaced progressive encoding.


Overview

We will assume all original images are 8 bit grayscale. Conversion from RGB to $\text{Y}'\text{C}_\text{B}\text{C}_\text{R}$ (luma + chroma components) will be skipped.

We will be assuming no options are passed and our quantization matrix is $Q_{i,j} = 1 + (1 + i + j) \cdot \text{quality factor}$.

File Structure

GPJ file structure

File Tag: 7 byte ASCII string ("GPJFile")

Option: 3 byte ASCII string denoting version and options ("000")

Quality: 3 byte ASCII string denoting quality factor used for coding image blocks

Width: 7 byte ASCII string denoting horizontal dimension of image

Height: 7 byte ASCII string denoting vertical dimension of image

M: 1 byte delimiter marker (0x00)

Each block represents an 8px by 8px image block. Blocks are ordered left to right, top to bottom.

Algorithm

  1. Convert each block to [-128, 128].
  2. Perform 2D DCT on each block.
  3. Quantize each DCT coefficient.
  4. Perform entropy coding on each block.

Encode

Generate the compressed GPJ image.

function output = encode(filename, image, factor)  
% Compress an image to GPJ

output = -1;  
fid = fopen(filename, 'w');  
if fid == -1  
  disp('Error: Can not open file.');
  return;
end;

% Initialization
[y_dim, x_dim] = size(image);
Q = QuantMatrix(factor);

% Write file header
fprintf(fid, 'GPJFile');  
fwrite(fid, 0, 'uchar');  
fprintf(fid, '%3d', 0);  
fwrite(fid, 0, 'uchar');  
fprintf(fid, '%3d', factor);  
fwrite(fid, 0, 'uchar');  
fprintf(fid, '%7d', x_dim);  
fwrite(fid, 0, 'uchar');  
fprintf(fid, '%7d', y_dim);  
fwrite(fid, 0, 'uchar');

for y = 1:8:y_dim  
  for x = 1:8:x_dim
    % Convert each block to [-128, 128]
    block = double(image(x:x+7,y:y+7)) - 128;

    % Perform DCT
    dct_block = dct2(block);

    % Quantization
    quantized = round(dct_block./Q);

    % Zigzag coding
    ind = reshape(1:numel(quantized), size(quantized));
    ind = fliplr(spdiags(fliplr(ind)));
    ind(:,1:2:end) = flipud(ind(:,1:2:end));
    ind(ind==0) = [];
    zigzag = quantized(ind);
    clear out;

    % Huffman coding
    counter = 1;
    numzero = 0;
    for(i = 1:length(zigzag))
      if(i >= counter && counter <= length(zigzag))
        if(zigzag(counter) == 0)
          while(zigzag(counter) == 0)
            numzero = numzero + 1;
            if(counter + 1 <= length(zigzag))
              counter = counter + 1;
            else
              break;
            end;
          end;
          fwrite(fid,0,’int8’);
          fwrite(fid,numzero,’int8’);
        end;
        if(abs(zigzag(counter)) <= 127 && zigzag(counter) ~= 0)
          fwrite(fid,zigzag(counter),’int8’);
        end;
        if(abs(zigzag(counter)) > 127)
          fwrite(fid,-128,’int8’);
          fwrite(fid,zigzag(counter),’int16’);
        end;
        counter = counter + 1;
      end;
      numzero = 0;
    end;
  end;
end;

% Close file
fclose(fid);  
output = 0;  
return;  

Generate the quantization matrix.

function Q = QuantMatrix(quant)  
% Generate the quantizer matrix

y = (0:7)';  
Q = [y y y y y y y y];  
Q = Q + Q';  
Q = (Q + 1) * quant + 1;  
return;  

Decode

Undo GPJ image compression.

function img = decode(filename)  
% Decompress GPJ to image

img = -1;  
fid = fopen(filename, 'r');  
if fid == -1  
  disp('Error: Can not open file.');
  return;
end;

% Read the file header
field = fscanf(fid, '%s', 1);  
idTag = sscanf(field, '%s\0', 1);  
field = fscanf(fid, '%s', 1);  
version = sscanf(field, '%d\0', 1);  
if ((idTag ~= 'GPJFile') | (version ~= 0))  
  disp('Error: GPJ file corrupted or wrong version.');
  return;
end  
field = fscanf(fid, '%s', 1);  
quant = sscanf(field, '%d\0');  
field = fscanf(fid, '%s', 1);  
x_dim = sscanf(field, '%d\0');  
field = fscanf(fid, '%s', 1);  
y_dim = sscanf(field, '%d\0');  
fseek(fid, 32, 'bof');

% Initialization
img = zeros(y_dim, x_dim);  
Q = GenerateQ(quant);  
for y = 1:8:y_dim  
  for x = 1:8:x_dim
    % Huffman decoding
    outval = 1;
    num = 1;
    clear val1 val2 zigzag;
    while(num <= 64);
      val1 = fread(fid,1,’int8’);
      num = num + 1;
      if(val1 == 0)
        val2 = fread(fid,1,’int8’);
        num = num + val2;
          for(i = outval:outval+val2-1)
            zigzag(i) = 0;
          end;
          outval = outval + val2;
        end;
        if(val1 == -128)
          val2 = fread(fid,1,’int16’);
          zigzag(outval) = val2;
          outval = outval + 1;
        end;
        if(abs(val1) <= 127 && val1 ~= 0)
          zigzag(outval) = val1;
          outval = outval + 1;
        end;
      end;

      % Zigzag decoding
      rezig = izigzag(zigzag);

      % Dequantization 
      rezig = rezig.*Q;

      % Perform inverse DCT
      rezig = idct2(rezig);

      % Convert each block to [0, 255]
      rezig = rezig + 128;
      image(x:x+7,y:y+7) = rezig;
    end;
  end;
end;

% Close file and display image
fclose(fid);  
img = uint8(imcrop(img, [1, 1, x_dim - 1, y_dim - 1]));  
return;  

Generate the inverse zigzag (credit Alexey Sokolov).

function output = izigzag(input)  
% Generate matrix from zigzag

% Initialization
h = 1;  
v = 1;  
vmin = 1;  
hmin = 1;  
output = zeros(8, 8);  
i = 1;

while ((v <= 8) & (h <= 8))  
  if(mod(h + v, 2) == 0)
    if(v == vmin)
      output(v, h) = input(i);
      if(h == 8)
        v = v + 1;
      else
        h = h + 1;
      end;
      i = i + 1;
    elseif((h == 8) & (v < 8))
      output(v, h) = input(i);
      v = v + 1;
      i = i + 1;
    elseif((v > vmin) & (h < 8))
      output(v, h) = input(i);
      v = v - 1;
      h = h + 1;
      i = i + 1;
    end;  
  else
    if((v == 8) & (h <= 8))
      output(v, h) = input(i);
      h = h + 1;
      i = i + 1;
    elseif(h == hmin)
      output(v, h) = input(i);
      if(v == 8)
        h = h + 1;
      else
        v = v + 1;
      end;
      i = i + 1;
    elseif ((v < 8) & (h > hmin))
      output(v, h) = input(i);
      v = v + 1;
      h = h - 1;
      i = i + 1;
    end;
  end;
  if((v == 8) & (h == 8))
    output(v, h) = input(i);
    break;
  end;
end;  

Analysis

For a sample image, we can vary the quality factor and trade-off file size for PSNR.

Comparison of compressed images with different quality factors

Quality factor File size PSNR
1 257,730 bytes 39.6460
2 210,908 bytes 34.4921
4 141,556 bytes 30.2566
8 86,552 bytes 27.0617
16 51,214 bytes 24.5101
32 29,725 bytes 22.4995

A better quantization table can be constructed by having finer steps for low frequency components and coarser steps for high frequency components.

$\begin{bmatrix}13 & 25 & 33 & 35 & 55 & 60 & 69 & 70 \\\\25 & 37 & 35 & 37 & 60 & 65 & 70 & 71 \\\\33 & 35 & 37 & 39 & 65 & 70 & 71 & 72 \\\\35 & 37 & 39 & 41 & 70 & 75 & 72 & 73 \\\\55 & 60 & 65 & 70 & 75 & 80 & 73 & 74 \\\\60 & 65 & 70 & 75 & 80 & 85 & 74 & 75 \\\\69 & 70 & 71 & 72 & 73 & 74 & 75 & 76 \\\\70 & 71 & 72 & 73 & 74 & 75 & 76 & 77 \end{bmatrix}$

comments powered by Disqus