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
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
- Convert each block to [-128, 128].
- Perform 2D DCT on each block.
- Quantize each DCT coefficient.
- 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.
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}$