Discrete Wavelets with Applications

122 - JPEG2000 Lossy Transform: 
(9,7) Cohen-Daubechies-Feauveau Filter Pair (CDF97) and Step Quantization

Helmut  Knaust - Department of Mathematical Sciences, University of Texas at El Paso, hknaust@utep.edu
12/1/2010; last edits 6/21/2013

SetDirectory[NotebookDirectory[]];
<<DW.m;

The CDF97 filter pair

P[t_]=1+4t+10t^2+20t^3;
sol=NSolve[P[t]==0];r=sol[[1,1,2]];r1=sol[[2,1,2]];r2=sol[[3,1,2]];a=-1/r;
H[t]=Sqrt[2]Cos[t/2]^4a (Sin[t/2]^2-r);h=Expand[TrigReduce[H[t]]/.{Cos[t]->(Exp[I t]+Exp[-I t])/2}/.{Cos[2t]->(Exp[2I t]+Exp[-2I t])/2}/.{Cos[3t]->(Exp[3I t]+Exp[-3I t])/2}];
H2[t_]=Sqrt[2]Cos[t/2]^4 20/a (Sin[t/2]^2-r1)(Sin[t/2]^2-r2);h2=Chop[Expand[TrigReduce[H2[t]]/.{Cos[t]->(Exp[I t]+Exp[-I t])/2}/.{Cos[2t]->(Exp[2I t]+Exp[-2I t])/2}/.{Cos[3t]->(Exp[3I t]+Exp[-3I t])/2}/.{Cos[4t]->(Exp[4I t]+Exp[-4I t])/2}]];
f={h[[1]],h[[2,1]],h[[4,1]],h[[6,1]]};
Print["The short CDF97 lowpass filter hs= ",hs=Drop[Join[Reverse[f],f],{4}]]
f2={h2[[1]],h2[[2,1]],h2[[4,1]],h2[[6,1]],h2[[8,1]]};
Print["The long CDF97 lowpass filter hl= ",hl=Drop[Join[Reverse[f2],f2],{6}]]
Print["The short CDF97 highpass filter gs= ",gs=Table[(-1)^k,{k,1,7}]*hs]
Print["The long CDF97 highpass filter gl= ",gl=Table[(-1)^(k+1),{k,1,9}]*hl]

GraphicsGrid[{{ListPlot[hl,PlotStyle->AbsolutePointSize[10],ImageSize->400,PlotLabel->"The coefficients of the long CDF97 lowpass filter"],
ListPlot[hs,PlotStyle->AbsolutePointSize[10],ImageSize->400,PlotLabel->"The coefficients of the short CDF97 lowpass filter"]}}]
GraphicsGrid[{{FourierFilterPlot[hl,PlotLabel->"The long CDF97 lowpass filter",ImageSize->400],
FourierFilterPlot[hs,PlotLabel->"The short CDF97 lowpass filter",ImageSize->400]
}}]

The CDF97 wavelet transform matrices

CDF1[12]//Normal
CDF2[12]//Normal
CDF1[12].Transpose[CDF2[12]]//Chop//Normal

Applying the CDF97 wavelet transform

We subtract 128 from all entries of the image.

pic=Import["Images\\Grayscale\\train.png","Data"]-128;
ImagePlot[pic,Scaled->True]
e0=ShannonEntropy[pic]

We define and then apply the CDF97 transfrom n times :

CDF97WLT[pic_,n_:1]:=Module[{pc=pic,r=Dimensions[pic][[1]],c=Dimensions[pic][[2]]},
	If[IntegerQ[n]&&n>=0,
		If[Mod[r,2^n]==0&&Mod[c,2^n]==0,
		Do[pc[[1;;r/2^k,1;;c/2^k]]=CDF2[r/2^k].pc[[1;;r/2^k,1;;c/2^k]].Transpose[CDF1[c/2^k]],{k,0,n-1}],
		Print["To perform ",n," iterations, the number of rows and columns must be divisible by ",2^n,"."]],
	Print["The parameter 'n' must be a non-negative integer."]];Return[pc];]

Intermezzo : Applying the CDF97 transform once reveals the nice edge detection capability of this filter :

ImagePlot[CDF97WLT[pic,1],Scaled->True]

n=3;
p3=CDF97WLT[pic,n];
ImagePlot[p3,Scaled->True]

Step quantization

Now we define the Step Quantization with parameter τ. Larger τ will lead to higher compression.

StepQuantize[pic_,n_,τ_]:=Module[{pc=pic,r=Dimensions[pic][[1]],c=Dimensions[pic][[2]],d=τ},
	Print["τ = ",τ];Print[d];
Do[pc[[1;;r/2^k,c/2^k+1;;c/2^(k-1)]]=Map[(Sign[#]Floor[Abs[#]/d]&),pc[[1;;r/2^k,c/2^k+1;;c/2^(k-1)]],{2}];
pc[[r/2^k+1;;r/2^(k-1),1;;c/2^k]]=
Map[(Sign[#]Floor[Abs[#]/d]&),pc[[r/2^k+1;;r/2^(k-1),1;;c/2^k]],{2}];
pc[[r/2^k+1;;r/2^(k-1),c/2^k+1;;c/2^(k-1)]]=Map[(Sign[#]Floor[Abs[#]/(2d)]&),pc[[r/2^k+1;;r/2^(k-1),c/2^k+1;;c/2^(k-1)]],{2}];d=d/2;,
{k,1,n-1}];
pc[[1;;r/2^n,1;;c/2^n]]=Map[(Sign[#]Floor[2Abs[#]/d]&),pc[[1;;r/2^n,1;;c/2^n]],{2}];
pc[[r/2^n+1;;r/2^(n-1),1;;c/2^n]]=Map[(Sign[#]Floor[Abs[#]/d]&),pc[[r/2^n+1;;r/2^(n-1),1;;c/2^n]],{2}];
pc[[1;;r/2^n,c/2^n+1;;c/2^(n-1)]]=
Map[(Sign[#]Floor[Abs[#]/d]&),pc[[1;;r/2^n,c/2^n+1;;c/2^(n-1)]],{2}];
pc[[r/2^n+1;;r/2^(n-1),c/2^n+1;;c/2^(n-1)]]=Map[(Sign[#]Floor[Abs[#]/(2d)]&),pc[[r/2^n+1;;r/2^(n-1),c/2^n+1;;c/2^(n-1)]],{2}]
;Return[pc];]

τ=20;
pq=StepQuantize[p3,n,τ];
Print["Entropy: ",e1=ShannonEntropy[pq]]
Print["Compression factor: ",cf=Round[e0/e1,0.1]]

This matrix is then encoded using a procedure similar to Huffman encoding. Afterwards the encoded information is sent to the recipient, together with a small overhead (encoding table, value of τ, dimensions of the matrix, etc.).

Decoding by the recipient

The recipient then reverses the quantization (lossy)  and the transform procedure:

InvStepQuantize[pic_,n_,τ_]:=Module[{pc=pic,r=Dimensions[pic][[1]],c=Dimensions[pic][[2]],d=τ},
	Print["τ = ",τ];
Do[pc[[1;;r/2^k,c/2^k+1;;c/2^(k-1)]]=Map[d #&,pc[[1;;r/2^k,c/2^k+1;;c/2^(k-1)]],{2}];
pc[[r/2^k+1;;r/2^(k-1),1;;c/2^k]]=
Map[d#&,pc[[r/2^k+1;;r/2^(k-1),1;;c/2^k]],{2}];
pc[[r/2^k+1;;r/2^(k-1),c/2^k+1;;c/2^(k-1)]]=Map[2d #&,pc[[r/2^k+1;;r/2^(k-1),c/2^k+1;;c/2^(k-1)]],{2}];d=d/2;,
{k,1,n-1}];
pc[[1;;r/2^n,1;;c/2^n]]=Map[d/2 #&,pc[[1;;r/2^n,1;;c/2^n]],{2}];
pc[[r/2^n+1;;r/2^(n-1),1;;c/2^n]]=Map[d # &,pc[[r/2^n+1;;r/2^(n-1),1;;c/2^n]],{2}];
pc[[1;;r/2^n,c/2^n+1;;c/2^(n-1)]]=
Map[d # &,pc[[1;;r/2^n,c/2^n+1;;c/2^(n-1)]],{2}];
pc[[r/2^n+1;;r/2^(n-1),c/2^n+1;;c/2^(n-1)]]=Map[(2d # &),pc[[r/2^n+1;;r/2^(n-1),c/2^n+1;;c/2^(n-1)]],{2}]
;Return[pc];]

p=InvStepQuantize[pq,n,τ];

We now apply the inverse transform three times to see the resulting image. 
	Note: The next command is NOT contained in the package 'DW.m'.

InvCDF97WLT[pic_,n_:1]:=Module[{pc=pic,r=Dimensions[pic][[1]],c=Dimensions[pic][[2]]},
	Do[pc[[1;;r/2^k,1;;c/2^k]]=Transpose[CDF1[r/2^k]].pc[[1;;r/2^k,1;;c/2^k]].CDF2[c/2^k],{k,n-1,0,-1}];
	Return[pc];]

Below is the image the recipient sees :

picq=InvCDF97WLT[p,n];
ImagePlot[128+picq,Scaled->True]

PSNR[pic,picq]

Here is a plot of the difference between the original image and its quantized version (negative) - both unscaled and scaled :

ImagePlot[255-Abs[pic-picq],Scaled->False]
ImagePlot[255-Abs[pic-picq],Scaled->True]

An enlarged detail: the original image is on the left, the recipient' s image is on the right.

Grid[{{ImagePlot[128+pic[[200;;450,130;;330]],Scaled->True,PixelConstrained->2],ImagePlot[128+picq[[200;;450,130;;330]],Scaled->True,PixelConstrained->2]}}]