The mathematics behind JPEG

©2015 by Gerald Teschl  <Gerald.Teschl@univie.ac.at>, <http://www.mat.univie.ac.at/~gerald/>

The purpose of this notebook is to visualize some of the main mathematical ideas behind the JPEG compression algorithm.

The idea behind the image compression in JPEG is the fact that storing the individual pixels is not adapted to the perception of the human eye. In fact, the human eye is more sensitive to graduate changes in intensity and color than to sharp transitions. Hence one tries to switch from the standard basis to a new basis where the graduate changes correspond to the lower coefficients and the faster changes to the higher coefficients. Dropping some of the higher coefficients results in a data reduction which is hardly visible. Nevertheless it is important to emphasize that the dropped coefficients are lost and cannot be restored, that is, JPEG is not a lossless compression.

The discrete cosine transform

The basis of choice is the discrete cosine transform with basis vectors given by:

u[n_,k_]:=Table[(Sqrt[(2/n)])Cos[(((2j-1)(k-1)π)/(2n))],{j,n}]

Note that the first vector is constant while the consecutive ones oscillate faster and faster:

Manipulate[ListPlot[Table[(Sqrt[(2/n)])Cos[(((2j-1)(k-1)π)/(2n))],{j,n}],PlotRange->{-(Sqrt[(2.1/n)]),(Sqrt[(2.1/n)])}],{k,1,n,1},Initialization:>{n=12;}]

Since our image has two spatial direction we will keep this structure and work with matrices rather than vectors. To get a basis for matrices we take the crosswise product of the above basis vectors (i.e. their tensor product):

u[{n1_,n2_},{k1_,k2_}]:=Table[u[n1,k1][[j1]]u[n2,k2][[j2]],{j1,n1},{j2,n2}]

Here is a graphical vizualization of these basis vectors

Manipulate[ArrayPlot[Table[(2/n)Cos[(((2j1-1)(k1-1)π)/(2n))]Cos[(((2j2-1)(k2-1)π)/(2n))],{j1,n},{j2,n}],ImageSize->250,PlotRange->{-(2/n),(2/n)},Frame->False],{k1,1,n,1},{k2,1,n,1},Initialization:>{n=12;}]

Since these basis vectors are orthogonal the expansion coefficients are just the projections and can be computed via the scalar product. The mapping from the original coordinates to the new coordinates is known as the discrete cosine transformation (DCT) and is built into Mathematica as FourierDCT[ ] and the inversion can be done with FourierDCT[ ,3].

Applications to image compression

Our data compression will be done by dropping a given percentage of the highest coefficients of the DCT. To this end the following function will truncate a given matrix setting p percent of its coefficients to zero (it will never drop the last one to avoid errors):

Manipulate[MatrixForm[truncate[Table[1,{10},{10}],p]],{p,0,99},Initialization:>{truncate[data_,p_]:=
Module[{i,j},
{i,j}=Map[Max[#,1]&,Floor[Dimensions[data](Sqrt[((100-p)/100)])]];
PadRight[Take[data,i,j],Dimensions[data],0]
];}]

We are going to use the following grayscale image to play with:

img=255-ExampleData[{"TestImage","Boat"},"Data"];
ArrayPlot[img,ImageSize->200,Frame->False]

Now we take our image, compute the expansion coefficients with respect to the DCT and throw away a given percentage of the coefficients as above. This results in a compression of the original image data. The resulting image together with the compressed one is visualized.

Manipulate[Grid[{Map[Style[#,"Label"]&,{"original",StringJoin[ToString[p],"% compressed"]}],{oplot,ArrayPlot[Map[Max[0,Min[#,255]]&,FourierDCT[truncate[dct,p],3],{2}],ImageSize->400,PlotRange->{0,255},Frame->False]}}],{{p,0,"compression percentage"},0,99},Initialization:>{
img=255-ExampleData[{"TestImage","Boat"},"Data"];
oplot=ArrayPlot[img,ImageSize->400,PlotRange->{0,255},Frame->False];
dct=FourierDCT[img];
truncate[data_,p_]:=
Module[{i,j},
{i,j}=Map[Max[#,1]&,Floor[Dimensions[data](Sqrt[((100-p)/100)])]];
PadRight[Take[data,i,j],Dimensions[data],0]
];}]

This is the main idea of the JPEG algorithm. However, in the real algorithm the color information is represented as one luma (the intensity) and two chroma numbers. Each of them is treated separately but the color information is compressed more as the eye is more sensitive to intensities. Moreover, the image is split into 8x8 blocks which are treated individually. Also the truncation is done in a more sophisticated way: the coefficients are stored with decreasing accuracy. Finally the coefficients are additionally compressed using a lossless Huffman encoding.

You could also take a portion of the above image:

img = 255 - ExampleData[{"TestImage", "Boat"}, "Data"][[120 ;; 220, 120 ;; 220]];
ArrayPlot[img, ImageSize -> 200, PlotRangePadding -> None, Frame -> False]

Our you could take a color image and convert it to grayscales:

img = Map[255 - {.3,.59,.11}.# &, ExampleData[{"TestImage", "Lena"}, "Data"], {2}][[200 ;; 299, 250 ;; 349]];
ArrayPlot[img, ImageSize -> 200, PlotRangePadding -> None, Frame -> False]