(* ::Package:: *)

(* DW1.wl - Discrete Wavelets package *) 
(* Author: Helmut Knaust *)
(* Address: Department of Mathematical Sciences, The University of Texas at El Paso, El Paso TX 79968-0514 *)
(* Email: hknaust@utep.edu *)
(* Last edit: 11/12/2015 *)



BeginPackage["DW`"]

D6L::usage=
	"D6[n] is the shifted Daubechies-6 Wavelet Transform of even dimension 'n' in numerical sparse array form."

D6LWLT::usage=
	"D6LWLT[pic,n] applies 'n' iterations of the shifted two-dimensional Daubechies 6-tab Wavelet Transform to the matrix 'pic'. The resulting matrix is not rounded to the closest integer. 
	To perform n iterations, the number of rows and columns of the matrix must be divisible by \!\(\*SuperscriptBox[\"2\", \"n\"]\). The default is n=1."  

HuffmanTree::usage=
	"HuffmanTree[v] displays a Huffman tree of the vector v."

HuffmanTreeStep::usage=
	"HuffmanTreeStep[v,k] displays the kth step in the construction of a Huffman tree of the vector v."

InvCDFWLT::usage=
	"InvD6WLT[pic,n] applies 'n' iterations of the inverse of the two-dimensional (9,7)Cohen-Daubechies-Feaveau (CDF97) Wavelet Transform to the matrix 'pic'. The resulting matrix is not rounded to the closest integer. 
	The default is n=1."  
	
InvD4WLT::usage=
	"InvD4WLT[pic,n] applies 'n' iterations of the inverse of the two-dimensional Daubechies-4 Wavelet Transform to the matrix 'pic'. The resulting matrix is not rounded to the closest integer. 
	The default is n=1."  

InvD6WLT::usage=
	"InvD6WLT[pic,n] applies 'n' iterations of the inverse of the two-dimensional Daubechies-6 Wavelet Transform to the matrix 'pic'. The resulting matrix is not rounded to the closest integer. 
	The default is n=1."  

InvD6LWLT::usage=
	"InvD6LWLT[pic,n] applies 'n' iterations of the inverse of the shifted two-dimensional Daubechies-6 Wavelet Transform to the matrix 'pic'. The resulting matrix is not rounded to the closest integer. 
	The default is n=1."  

InvHaarWLT::usage=
	"InvHaarWLT[pic,n] applies 'n' iterations of the inverse of the two-dimensional Haar Wavelet Transform to the matrix 'pic'. The resulting matrix is not rounded to the closest integer. 
	The default is n=1."  


Begin["`Private`"]

D6L[n_]:=Module[{h0=(Sqrt[2.]/32)(1+Sqrt[10]+Sqrt[5+2Sqrt[10]]),h1=(Sqrt[2.]/32)(5+Sqrt[10]+3Sqrt[5+2Sqrt[10]]),
	h2=(Sqrt[2.]/32)(10-2Sqrt[10]+2Sqrt[5+2Sqrt[10]]),h3=(Sqrt[2.]/32)(10-2Sqrt[10]-2Sqrt[5+2Sqrt[10]]),
	h4=(Sqrt[2.]/32)(5+Sqrt[10]-3Sqrt[5+2Sqrt[10]]),h5=(Sqrt[2.]/32)(1+Sqrt[10]-Sqrt[5+2Sqrt[10]])},
	If[EvenQ[n]&&n>=6,Transpose[RotateLeft[Transpose[Join[SparseArray[{{i_,j_}/;j==2i-1->h5,{i_,j_}/;j==2i->h4,{i_,j_}/;j==2i+1->h3,{i_,j_}/;j==2i+2->h2,{i_,j_}/;j==2i+3->h1,{i_,j_}/;j==2i+4->h0,{n/2-1,1}->h1,{n/2-1,2}->h0,{n/2,1}->h3,{n/2,2}->h2,{n/2,3}->h1,{n/2,4}->h0},{n/2,n}],
	SparseArray[{{i_,j_}/;j==2i-1->-h0,{i_,j_}/;j==2i->h1,{i_,j_}/;j==2i+1->-h2,{i_,j_}/;j==2i+2->h3,{i_,j_}/;j==2i+3->-h4,{i_,j_}/;j==2i+4->h5,{n/2-1,1}->-h4,{n/2-1,2}->h5,{n/2,1}->-h2,{n/2,2}->h3,{n/2,3}->-h4,{n/2,4}->h5},{n/2,n}]]],2]],
	Print["The dimension of the D6 matrix must be a positive  even integer \[GreaterEqual] 6."]]]

D6LWLT[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]]=D6L[r/2^k] . pc[[1;;r/2^k,1;;c/2^k]] . Transpose[D6L[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];]


Options[HuffmanTree]=Options[TreePlot];
HuffmanTree[v_,Options___]:=Module[{g,h,g2,n=0,sym},g=Sort[Union[Table[{sym,Count[v,sym]},{sym,v}]],#1[[2]]<=#2[[2]]&];
	h=Map[Function[t,t->t],g];While[Not[Length[g]==1],n=n+1;g2=Prepend[g,{StringJoin["\[FivePointedStar]",ToString[n]],g[[1,2]]+g[[2,2]]}];
	h=Prepend[h,g2[[1]]->g2[[2]]];h=Prepend[h,g2[[1]]->g2[[3]]];g2=Drop[g2,{2,3}];g=Sort[g2,#1[[2]]<=#2[[2]]&]];
	Print[TreePlot[h,Automatic,h[[1,1]],SelfLoopStyle->None,VertexLabeling->True,BaseStyle->{12,Thick},
	VertexRenderingFunction->({Blue,EdgeForm[Black],Disk[#,0.1],White,
	Text[If[StringQ[#2[[1]]],StringReplace[#2[[1]],"\[FivePointedStar]"~~__->""],#2[[1]]],#1]}&)],Options]]

Options[HuffmanTreeStep]=Options[TreePlot];
HuffmanTreeStep[v_,k_,Options___]:=Module[{g,h,g2,n=0,sym,iteration},g=Sort[Union[Table[{sym,Count[v,sym]},{sym,v}]],#1[[2]]<=#2[[2]]&];
	h=Map[Function[t,t->t],g];Do[n=n+1;g2=Prepend[g,{StringJoin["\[FivePointedStar]",ToString[n]],g[[1,2]]+g[[2,2]]}];
	h=Prepend[h,g2[[1]]->g2[[2]]];h=Prepend[h,g2[[1]]->g2[[3]]];g2=Drop[g2,{2,3}];g=Sort[g2,#1[[2]]<=#2[[2]]&];,{iteration,1,k}];
	TreePlot[h,Automatic,h[[1,1]],SelfLoopStyle->None,VertexLabeling->True,BaseStyle->{14,Thick},Options]]

InvCDFWLT[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[CoDF1[r/2^k]] . pc[[1;;r/2^k,1;;c/2^k]] . CoDF2[c/2^k],{k,n-1,0,-1}];
	Return[pc];]

InvD4WLT[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[D4[r/2^k]] . pc[[1;;r/2^k,1;;c/2^k]] . D4[c/2^k],{k,n-1,0,-1}];
	Return[pc];]

InvD6WLT[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[D6[r/2^k]] . pc[[1;;r/2^k,1;;c/2^k]] . D6[c/2^k],{k,n-1,0,-1}];
	Return[pc];]

InvD6LWLT[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[D6L[r/2^k]] . pc[[1;;r/2^k,1;;c/2^k]] . D6L[c/2^k],{k,n-1,0,-1}];
	Return[pc];]

InvHaarWLT[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[Haar[r/2^k]] . pc[[1;;r/2^k,1;;c/2^k]] . Haar[c/2^k],{k,n-1,0,-1}];
	Return[pc];]

End[ ]
EndPackage[ ]