Discrete Cosine Transform
(and the JPEG Standard)
Continued

M. E. Mauel
Columbia University

Introduction

The discrete Cosine transfrom (DCT) is commonly used in digital signal processing. It converts a sequence of real numbers to a sequence of real numbers of the same length. (Recall, the DFT converts a sequence of real numbers (or complex) to a sequence of complex numbers.) The use of the DCT for image compression was introduced in 1974, and it plays an important role in JPEG and MPEG compression schemes.

Formally, the DCT of a digital signal, (Subscript[x, i]) == {(Subscript[x, 1]), (Subscript[x, 2]), …, (Subscript[x, N])}, is
	(Subscript[xc, k]) = (∑_(i=1)^N)(Subscript[x, i]) Cos[(π/(2N))(2i-1)(k-1)]
The DCT is exactly equivalent to the discrete Fourier transform (DFT) of 2N - 2 real numbers with even symmetry. (Therefore, we can use tyhe built-in Fourier[…] command to compute the DCT.)

In this notebook, the DCT and the processing of two dimensional graphics is examined. In particular, we explore

I've used two references. First is "Image Compression Using the Discrete Cosine Transform" by Andrew Watson, Mathematica Journal 4 (1994) p. 81. This is a wonderful collection of Mathematica routines from which I have based this lecture. Second is "The JPEG Still Picture Compression Standard" by Gregory Wallace, Communications of the ACM (1991) p. 35. Both articles are available online on our class website; however, you do not need to read these articles in order to understand this notebook.

2D DCT Modes

The JPEG standard is based on the DCT of 8 x 8 blocks of graphic data. The modes of an 8 x 8 block can be calculated easily.

dctMatrix[nSize_:8] := Table[If[m == 0, Sqrt[1/2],1] Sqrt[2/nSize] Cos[(2 n +1) m π/(2 nSize)], {m, 0, nSize-1}, {n, 0, nSize-1}]

dctTensor = Array[(Outer[Times, dctMatrix[8][[#1]], dctMatrix[8][[#2]]])&, {8, 8}];

Dimensions[dctTensor]

dctTensor // Shallow

ArrayPlot[dctTensor[[2,2]],Mesh->False,PlotLabel->"DCT Mode [2,2]", ColorFunction -> GrayLevel]

Show[GraphicsGrid[Map[ArrayPlot[#1,Mesh->False,Frame->False,ColorFunction -> GrayLevel]&,dctTensor,{2}]],PlotLabel->"DCT Modes for an 8x8 Image"]

A DCT 2d Procedure

dctFourier[x_List] := Module[{n,e0, tFactors, xeven, xodd},
n = Length[x];
{xeven, xodd} = Transpose[Partition[x, 2]];
e0 = Exp[-I π/(2 n)];
tFactors = {1} ~ Join~ Table[Sqrt[2] e0^i, {i, 1, n-1}];
Re[tFactors InverseFourier[Join[xeven, Reverse[xodd]]]]]

idctFourier[x_List] :=  Module[{n,e0,tFactors, xt, xfirst, xlast},
n = Length[x];
e0 = Exp[-I π/(2 n)];
tFactors = {1} ~ Join~ Table[Sqrt[2] e0^i, {i, 1, n-1}];
tFactors = Conjugate[N[tFactors]];
xt = Re[Fourier[tFactors x]];
xfirst = Take[xt, n/2];
xlast = Take[xt,-n/2];
Flatten[Transpose[{xfirst, Reverse[xlast]}]]]

dct2D[array_?MatrixQ] := Transpose[dctFourier /@ Transpose[dctFourier /@ array]];

idct2D[array_?MatrixQ] := Transpose[idctFourier /@ Transpose[idctFourier /@ array]];

iSize=(2^6);
image1=Table[N[Sin[((π i 3)/iSize)] Sin[((π j 2)/iSize)]],{i,1,iSize},{j,1,iSize}];
MatrixPlot[image1,Mesh->False,PlotLabel->"Test Image (64x64)", ColorFunction -> GrayLevel]

image1DCT=dct2D[image1];
MatrixPlot[image1DCT,Mesh->False,PlotRange->All,PlotLabel->"DCT of 64x64 Test Sample", ColorFunction -> GrayLevel]

ListPlot[{image1DCT[[All,1]],image1DCT[[All,2]],image1DCT[[All,3]]}, PlotRange -> All, Joined -> True];
ListPlot[{image1DCT[[All,1]],image1DCT[[All,2]],image1DCT[[All,3]]}, PlotRange -> All, PlotStyle -> PointSize[0.02]];
Show[%,%%, PlotLabel -> "DCT Values"]

Blocked DCT 2D

JPEG encoding is performed on blocks of 8 x 8 data. In this section, a general purpose algorithm is defined that partitions a graphic into 8 x 8 blocks and then performs a DCT on the partitions. (This is based on the previously referenced article by Andrew Watson.)

Partition

Clear[a];

alist = Table[a[i], {i, 1, 36}]

Partition[alist, 8]

Blocked DCT Definitions (Function Overloading)

bDCTfactors = N[{1} ~ Join~ Table[Sqrt[2] Exp[-I π i/(2 8)], {i, 1, 8-1}]];
bIDCTfactors = Conjugate[bDCTfactors];
bDCT[list_] := Re[ bDCTfactors InverseFourier[ N[list[[{1,3,5,7,8,6,4,2}]]]]];
bIDCT[list_] := Re[Fourier[bIDCTfactors list]][[{1,8,2,7,3,6,4,5}]];
bDCT[array_?MatrixQ] := Transpose[bDCT /@ Transpose[bDCT /@ array]];
bIDCT[array_?MatrixQ] := Transpose[bIDCT /@ Transpose[bIDCT /@ array]];
bDCT[list_?(Length[#] > 8&)] := Join @@ (bDCT /@ Partition[list,8]);
bIDCT[list_?(Length[#] > 8&)] := Join @@ (bIDCT /@ Partition[list,8]);

?bDCT

?bIDCT

bimage1 = bDCT[image1];

ArrayPlot[bimage1,Mesh->False,PlotRange->All,PlotLabel->"Block DCT of (8x8)x(8x8) Image",ColorFunction -> GrayLevel]

Dylan Block DCT

Get Graphic Data

SetDirectory["~/Desktop/"]

FileNames["*.tif"]

g1 = Import["D_Thomas.tif"];

g1Data = ImageData[g1,"Byte"];

Dimensions[g1Data]

ArrayPlot[g1Data,Mesh->False,AspectRatio->Automatic, ColorFunction -> GrayLevel]

Print["Min/Max = ",Min[g1Data],", ",Max[g1Data]];MatrixPlot[g1Data, ColorFunction -> GrayLevel]

DCT

bg1 = bDCT[g1Data];
Dimensions[bg1]

Print["Min/Max = ",Min[bg1],", ",Max[bg1]];ArrayPlot[bg1,AspectRatio->Automatic,PlotRange->All, ColorFunction -> GrayLevel]

ArrayPlot[bg1,AspectRatio->Automatic, PlotRange -> {-200,200},ColorFunction -> GrayLevel]

big1Data = bIDCT[bg1];

g1Data - big1Data // Chop // Abs // Max

Information Entropy, Compress, and File Size

How much information is stored in our TIFF image? and in our DCT?

TIFF Image

g1Data // Shallow

Print["File size in Finder: ",28 2^10,", Number of Bytes: ",Length[Flatten[g1Data]]]

originalEntropy = Entropy[2,Flatten[g1Data]] // N
originalBytes = Ceiling[%  Length[Flatten[g1Data]]/8 ]
originalCompSize = ByteCount[Compress[g1Data]]

Block DCT of TIFF Image

bg1 // Shallow

Entropy[2,Flatten[bg1]] // N
Ceiling[%  Length[Flatten[bg1]]/8 ]
ByteCount[Compress[bg1]]

The DCT of an image converts 8-bit (grayscale) integers into floating-point representation. As defined, the DCT increases the information content needed by Mathematica to represent the image. 
For JPEG compression, some high-frequency information is discarded and the data is quantized prior to compression.

JPEG Quantization

In the JPEG image compression standard,each DCT coefficient is quantized using a weight that depends on the frequencies for that coefficient.The coefficients in each 8 x 8 block are divided by a corresponding entry of an 8 x 8 quantization matrix,and the result is rounded to the nearest integer.

In general,higher spatial frequencies are less visible to the human eye than low frequencies.Therefore,the quantization factors are usually chosen to be larger for the higher frequencies. (A mode is "quantized" by dividing by the quantization factor. A large number reduces the relative weight of the component.) The following quantization matrix is widely used for monochrome images and for the luminance component of a color image. According to Watson, the "de facto" quantization matrix for JPEG encoding is given by the following.

qLum =
{{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}};

MatrixPlot[qLum,Mesh->False,PlotRange->All,PlotLabel->"JPEG 8x8 Quantization Factors", ColorFunction -> GrayLevel]

Quantized Functions

BlockImage[image_, blocksize_:{8, 8}] :=Partition[image, blocksize] /;And @@ IntegerQ /@ (Dimensions[image]/blocksize)

UnBlockImage[blocks_] := Partition[Flatten[Transpose[blocks, {1, 3, 2}]],{Times @@ Dimensions[blocks][[{2, 4}]]}]

(tImage = Table[i + 8(j-1), {j, 4}, {i, 6}]) // MatrixForm

BlockImage[tImage, {2, 3}] // MatrixForm

UnBlockImage[BlockImage[tImage, {2, 3}]] // MatrixForm

bQuantize[image_, qMatrix_] := Map[(#/qMatrix)&, BlockImage[image, Dimensions[qMatrix]],{2}] //UnBlockImage // Round;

bUnQuantize[image_, qMatrix_] := Map[(# qMatrix)&, BlockImage[image, Dimensions[qMatrix]],{2}] //UnBlockImage

Examples

{Dimensions[g1Data], Dimensions[g1Data]/8}

bg1=bDCT[g1Data];
Print["Min/Max = ",Min[bg1],", ",Max[bg1]];
MatrixPlot[bg1,Mesh->False,AspectRatio->Automatic,Frame->False, ColorFunction -> GrayLevel]

bqg1=bQuantize[bg1,qLum];
Print["Min/Max = ",Min[bqg1],", ",Max[bqg1]];
MatrixPlot[bqg1,Mesh->False,AspectRatio->Automatic,Frame->False, ColorFunction -> GrayLevel]

buqg1=bUnQuantize[bqg1,qLum];
Print["Min/Max = ",Min[buqg1],", ",Max[buqg1]];
MatrixPlot[buqg1,Mesh->False,AspectRatio->Automatic,Frame->False, ColorFunction -> GrayLevel]

lp1=ListPlot[bg1[[Range[1,80],9]],Joined->True,PlotRange->{-10,256},PlotStyle->{Dashing[{0.1,0.1}],Red}];
lp2=ListPlot[buqg1[[Range[1,80],9]],Joined->True,PlotRange->{-10,256},PlotStyle->Blue];
Show[lp1,lp2,PlotLabel->"Comparison: (Un)Quantized Data"]

lp1=ListPlot[bg1[[Range[1,80],9]],Joined->True,PlotRange->{-10,25},PlotStyle->{Dashing[{0.1,0.1}],Red}];
lp2=ListPlot[buqg1[[Range[1,80],9]],Joined->True,PlotRange->{-10,25},PlotStyle->Blue];
Show[lp1,lp2,PlotLabel->"Comparison: (Un)Quantized Data"]

Data Compression: Quantization sets pixels to "zero"

Count[Flatten[bg1], _?(# == 0 &)]

Count[Flatten[buqg1], _?(# == 0 &)]

Histogram[Flatten[bg1],PlotLabel->"DCT Frequencies"]

Histogram[Flatten[buqg1],PlotLabel->"Quantized Frequencies"]

Default Compression

j1Data = g1Data // bDCT // bQuantize[#,qLum]& // bUnQuantize[#, qLum]&;

q1Data = j1Data // bIDCT;

Show[GraphicsRow[(ArrayPlot[#1,PlotRange->{0,255},Mesh->False,AspectRatio->Automatic, ColorFunction -> GrayLevel,Frame->False]&)/@{g1Data,q1Data}]]

Print["Entropy = ",Entropy[2,Flatten[j1Data]] // N,", Orig: ", originalEntropy];
Print["Min bytes = ",Ceiling[N[Entropy[2,Flatten[j1Data]] ]  Length[Flatten[j1Data]]/8 ],", Orig: ", originalBytes];
Print["Compress size = ",ByteCount[Compress[j1Data]],", Orig: ", originalCompSize];

Low Quality Compression

j1Data8 = g1Data // bDCT // bQuantize[#,8 qLum]& // bUnQuantize[#, 8 qLum]&;

q1Data8 = j1Data8 // bIDCT;

Show[GraphicsRow[(ArrayPlot[#1,PlotRange->{0,255},Mesh->False,AspectRatio->Automatic, ColorFunction -> GrayLevel,Frame->False]&)/@{g1Data,q1Data8}]]

Print["Entropy = ",Entropy[2,Flatten[j1Data8]] // N,", Orig: ", originalEntropy];
Print["Min bytes = ",Ceiling[N[Entropy[2,Flatten[j1Data8]] ]  Length[Flatten[j1Data8]]/8 ],", Orig: ", originalBytes];
Print["Compress size = ",ByteCount[Compress[j1Data8]],", Orig: ", originalCompSize];

High Quality Compression

j1Data025 = g1Data // bDCT // bQuantize[#,0.25 qLum]& // bUnQuantize[#, 0.25 qLum]&;

q1Data025 = j1Data025 // bIDCT;

Show[GraphicsRow[(ArrayPlot[#1,PlotRange->{0,255},Mesh->False,AspectRatio->Automatic,ColorFunction -> GrayLevel,Frame->False]&)/@{g1Data,q1Data025}]]

Print["Entropy = ",Entropy[2,Flatten[j1Data025]] // N,", Orig: ", originalEntropy];
Print["Min bytes = ",Ceiling[N[Entropy[2,Flatten[j1Data025]] ]  Length[Flatten[j1Data025]]/8 ],", Orig: ", originalBytes];
Print["Compress size = ",ByteCount[Compress[j1Data025]],", Orig: ", originalCompSize];

Summary

This notebook defines and illustrates the use of the discrete Cosine transform (DCT) for digital signal processing. The DCT transforms a sequence of Real numbers into another sequence of Real numbers. Since it is related to the Fourier transform (DFT), relatively fast algorithms exist that have been implemented both in hardware and software.

When the DCT is applied to two-dimensional images, digital filters and compression technqiues can be applied. The JPEG algorithm is applied to 8 x 8 blocks using a perception-based quantization. The quality of the image and the degree of compression can be adjusted with a single parameter.