Discrete Wavelets with Applications - Spring 2009

126 - JPEG

Helmut  Knaust - Department of Mathematical Sciences, University of Texas at El Paso, hknaust@utep.edu
05/06/09

625/8.

78 8

SetDirectory[NotebookDirectory[]];<<DW.m
pic=Import["ING.jpg","Data"];
ImagePlot[pic,Scaled->True]
{r,c}=Dimensions[pic]
pic[[100,100]]
pic=pic[[1;;416,1;;624]]-128;
pic[[100,100]]
ImagePlot[pic,Scaled->True]
ent=ShannonEntropy[pic]
{r,c}=Dimensions[pic]

Padding and splitting the image into 8 x 8 blocks

If necessary, the image is padded with zeros on the right and at the bottom  to ensure that the both dimensions are multiples of eight.
Then the image is split into 8x8 blocks (64 pixels each).

{r,c}=Dimensions[pic];
{rw,cl}={r,c};
pic=Join[pic,Table[0,{i,1,8Ceiling[r/8]-r},{j,1,c}]];
{r,c}=Dimensions[pic];
pic=Transpose[Join[Transpose[pic],Table[0,{j,1,8Ceiling[c/8]-c},{i,1,r}]]];
{r,c}=Dimensions[pic];
P=Partition[pic,{8,8}];
ImagePlot[P[[35,45]],PixelConstrained->40]
Dimensions[P][[1;;2]]

B1=pic[[1;;rw,1;;cl]];

P[[35,45]]

Creating the Discrete Cosine Fourier Transform (DCT) matrix

Below is the 8 x8 orthogonal matrix used by JPEG, both in symbolic and numerial form.

U=1/2Prepend[U1=Table[ Cos[(j-1)(2k-1)π/16],{j,2,8},{k,1,8}],U2=Table[Sqrt[2]/2,{8}]]
U=1/2.Prepend[U1,U2]
TeXForm[%]

(({{(1/(2 (Sqrt[2]))), (1/(2 (Sqrt[2]))), (1/(2 (Sqrt[2]))), (1/(2 (Sqrt[2]))), (1/(2 (Sqrt[2]))), (1/(2 (Sqrt[2]))), (1/(2 (Sqrt[2]))), (1/(2 (Sqrt[2])))}, {(1/2) cos((π/16)), (1/2) cos(((3 π)/16)), (1/2) sin(((3 π)/16)), (1/2) sin((π/16)), -(1/2) sin((π/16)), -(1/2) sin(((3 π)/16)), -(1/2) cos(((3 π)/16)), -(1/2) cos((π/16))}, {(1/2) cos((π/8)), (1/2) sin((π/8)), -(1/2) sin((π/8)), -(1/2) cos((π/8)), -(1/2) cos((π/8)), -(1/2) sin((π/8)), (1/2) sin((π/8)), (1/2) cos((π/8))}, {(1/2) cos(((3 π)/16)), -(1/2) sin((π/16)), -(1/2) cos((π/16)), -(1/2) sin(((3 π)/16)), (1/2) sin(((3 π)/16)), (1/2) cos((π/16)), (1/2) sin((π/16)), -(1/2) cos(((3 π)/16))}, {(1/(2 (Sqrt[2]))), -(1/(2 (Sqrt[2]))), -(1/(2 (Sqrt[2]))), (1/(2 (Sqrt[2]))), (1/(2 (Sqrt[2]))), -(1/(2 (Sqrt[2]))), -(1/(2 (Sqrt[2]))), (1/(2 (Sqrt[2])))}, {(1/2) sin(((3 π)/16)), -(1/2) cos((π/16)), (1/2) sin((π/16)), (1/2) cos(((3 π)/16)), -(1/2) cos(((3 π)/16)), -(1/2) sin((π/16)), (1/2) cos((π/16)), -(1/2) sin(((3 π)/16))}, {(1/2) sin((π/8)), -(1/2) cos((π/8)), (1/2) cos((π/8)), -(1/2) sin((π/8)), -(1/2) sin((π/8)), (1/2) cos((π/8)), -(1/2) cos((π/8)), (1/2) sin((π/8))}, {(1/2) sin((π/16)), -(1/2) sin(((3 π)/16)), (1/2) cos(((3 π)/16)), -(1/2) cos((π/16)), (1/2) cos((π/16)), -(1/2) cos(((3 π)/16)), (1/2) sin(((3 π)/16)), -(1/2) sin((π/16))}}))

The matrix U is indeed orthogonal :

Chop[U.Transpose[U]]

The matrix Z below is a normalization matrix.

Z={{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}}

P[[35,45]]
ShannonEntropy[P[[35,45]]]
U.P[[35,45]].Transpose[U]

ImagePlot[U.P[[35,45]].Transpose[U],PixelConstrained->40]

U.P[[35,45]].Transpose[U]
TeXForm[%]
ShannonEntropy[U.P[[35,45]].Transpose[U]]

Round[(U.P[[35,45]].Transpose[U])/Z]

ShannonEntropy[Round[(U.P[[35,45]].Transpose[U])/Z]]

Applying the Discrete Cosine Fourier Transform (DCT) matrix and Quantization

We apply U and its transform to each block of the image, then divide each entry in the resulting matrix by the corresponding element in Z; finally we round to the nearest integer.
We also list the entropy of the resulting image matrix and the compression factor obtained.

P1=Map[(U.#.Transpose[U]&),P,{2}];
P1[[35,45]]
P2=Round[Map[(#/Z&),P1,{2}]];
P2[[35,45]]
ent2=ShannonEntropy[P2]
Print["Compression Factor: ",Round[ent/ent2,.1]]

Below is a plot of the compressed matrix :

m=Min[P2];M=Max[P2];D1=255/(M-m) (P2-m);D2=Flatten[D1,{{1,3},{2,4}}];
P2a=Flatten[P2,{{1,3},{2,4}}];
ImagePlot[P2a]
ImagePlot[Abs[P2a]]
ShannonEntropy[P2]
Abs[P3[[100;;120,100;;120]]]

D2[[1;;10,1;;10]]

7.67/1.4

ImagePlot[Abs[D2]]

ImagePlot[P2[[35,45]],PixelConstrained->40,Scaled->True]

Decoding

We multiply each block elementwise by the matrix Z, then apply the transpose of U and U itself to the image:

P3=Map[(Z#&),P2,{2}];
P4=Map[(Transpose[U] .#.U&),P3,{2}];

Finally we join all blocks of the partition to obtain an approximation to the original image :

B2=128+Flatten[P4,{{1,3},{2,4}}];
B3=Map[(Min[Max[#,0],255]&),B2,{2}]-128;
ImagePlot[128+B3]
Grid[{{ImagePlot[128+pic],ImagePlot[128+B3]}}]
Print["PSNR: ",Round[PSNR[pic,B3],.1]]

Grid[{{ImagePlot[255-Abs[pic-B3],Scaled->False],ImagePlot[255-Abs[pic-B3],Scaled->True]}}]
Grid[{{ImagePlot[128+pic[[65;;184,161;;280]],PixelConstrained->3],ImagePlot[128+B3[[65;;184,161;;280]],PixelConstrained->3]}}]