(* ::Package:: *)

(* DW.m - 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`"]

CEnergy::usage=
	"CEnergy[data] plots the cumulative energy of 'data'. 
	'data' can be a vector, matrix or other array. All entries are rounded to the nearest integer before computation of the cumulative energy. 
	You can specify the same options as in 'ListPlot'."

CFrequency::usage=
	"CFrequency[data] plots the relative cumulative frequency of 'data'. 
	'data' can be a vector, matrix or other array. All entries are rounded to the nearest integer before computation of the cumulative energy. 
	You can specify the same options as in 'ListPlot'."
CDF1::usage=
	"CDF1[n] is one of the two square CDF(9,7) wavelet transform matrices of even dimension 'n' in symbolic sparse array form."
	
CDF2::usage=
	"CDF1[n] is one of the two square CDF(9,7) wavelet transform matrices of even dimension 'n' in symbolic sparse array form."

CDFWLT::usage=
	"CDFWLT[pic,n] applies 'n' iterations of the two-dimensional CDF(9,7) 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."  

ColorChannels::usage=
    "ColorChannels[color_image] returns a vector of the form {rcolor, gcolor,bcolor} where 'rcolor', 'gcolor', and 'bcolor' are the three color channels in color."

D4exact::usage=
	"D4exact[n] is the square Daubechies-4 Wavelet Transform matrix of even dimension 'n' in symbolic sparse array form."

D4::usage=
	"D4[n] is the square Daubechies-4 Wavelet Transform of even dimension 'n' in numerical sparse array form."

D4WLT::usage=
	"D4WLT[pic,n] applies 'n' iterations of the two-dimensional Daubechies 4-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."  

D6exact::usage=
	"D6exact[n] is the square Daubechies-6 Wavelet Transform of even dimension 'n' in symbolic sparse array form."

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

D6Threshold::usage=
	"D6Threshold[pic,\[Alpha]] applies the thresholding algorithm by Gonzalez and Woods to the matrix 'pic' with tolerance '\[Alpha]'. 
	First the Daubechies-6 wavelet transfrom is applied to 'pic' and its blur portion is set to 0. Then the thresholding value is determined. 
	Matrix entries in 'pic' below the threshold value \[Tau] are set to 0. Finally the inverse Daubechies-6 wavelet transform is applied once to visulaize the edges of 'pic'.
	The default value is \[Alpha]=1." 

D6WLT::usage=
	"D6WLT[pic,n] applies 'n' iterations of the 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."  

EnergySpectrum::usage=
	"EnergySpectrum[f(t),t,n] plots the first 'n' elements of the the energy spectrum for the function f(t). The default for 'n' is 6."

FourierFilterPlot::usage=
	"FourierFilterPlot[coeff] plots the absolute value of the complex Fourier expansion of the filter 'coeff'. 'coeff' must be a finite list." 

FourierPlot::usage=
	"FourierPlot[f(t),t,n] plots the function f(t) (dashed) and its 'n'-terms complex Fourier expansion on the interval [-\[Pi],\[Pi]]. The default for 'n' is 6."

FourierSC::usage=
	"FourierSC[f(t),t,terms] lists the first 'terms' of the real Fourier expansion for the function f(t) on the interval [-\[Pi],\[Pi]].
	The solution is given in closed form, so this command may run for a long time or return no results. The default for 'terms' is 6." 

FourierSCPlot::usage=
	"FourierPlot[f(t),t,n] plots the function f(t) (dashed) and its 'n'-terms real Fourier expansion on the interval [-\[Pi],\[Pi]]. The default for 'n' is 6."

FourierS::usage=
	"FourierS[f(t),t,terms] lists the first 'terms' of the complex Fourier expansion for the function f(t) on the interval [-\[Pi],\[Pi]]. 
	The solution is given in closed form, so this command may run for a long time or return no results. The default for 'terms' is 6." 

FrequencyHistogram::usage=
	"FrequencyHistogram[data] plots the absolute frequency distribution of 'data'. 
	'data' can be a vector, matrix, or 3-dimensional array. All entries are rounded to the nearest integer before computation of the histogram. 
	In case of a color image (3-dimensional array), separate histograms are displayed for each color channel."

GrayChannels::usage=
    "GrayChannels[color_image] returns a vector of the form {rgray, ggray, bgray} where  'rgray', ggray', and 'bgray' are the three color channels in graytones."

Haar::usage=
	"Haar[n] is the square Haar Wavelet Transform of even dimension 'n' in numerical sparse array form."

HaarThreshold::usage=
	"HaarThreshold[pic,\[Alpha]] applies the thresholding algorithm by Gonzalez and Woods to the matrix 'pic' with tolerance '\[Alpha]'. 
	First the Haar transfrom is applied to 'pic' and its blur portion is set to 0. Then the thresholding value is determined. 
	Matrix entries in 'pic' below the threshold value \[Tau] are set to 0. Finally the inverse Haar transform is applied once to visualize the edges of 'pic'.
	The default value is \[Alpha]=1." 

HaarWLT::usage=
	"HaarWLT[pic,n] applies 'n' iterations of the two-dimensional Haar 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."  

Huffman::usage=
	"Huffman[list] gives a table with Huffman encoding symbols; the first column of the table lists the original characters;
	 the second column their corresponding Huffman code; the last column the character frequency."

ImagePlot::usage=
	"ImagePlot[data] displays a matrix (array) as a grayscale (or color) image. 
	You can specify the same options as for 'ArrayPlot'. Particularly useful might be 'PixelConstrained\[Rule]x' to magnify the output. 'x' must be a positive integer; its default is 1.
	The option 'Scaled->False'  assumes that the matrix entries lie between 0 and 255. Values below 0 are depicted in 'Black', values above 255 are depicted in 'White'.
	'Scaled->True' (default) scales the display such that the minimum matrix entry corresponds to 0 and the maximum matrix entry corresponds to 255."

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."  

JoinGrayChannels::usage=
	"JoinGrayChannels[rgray,ggray,bgray] joins the three grayscale RGB channels of a color image to recreate the RGB color image. It is essential to input the three matrices in the correct (=RGB) order." 

NTSC::usage=
	"NTSC[image] converts a color image to the NTSC standard grayscale representation."

PicPlot::usage=
	"PicPlot[data,(magnification factor)] is a faster but less sophisticated version of 'ImagePlot', in particular for large color images.
		PicPlot assumes that all data lie between 0 and 255. An integer magnification factor (default is 1) can be given as an optional argument." 

PSNR::usage=
	"PSNR[image1,image2] computes the peak signal to noise ratio of two matrices representing (grayscale or color) images. Both image matrices must have the same dimensions; 
	the matrix elements should lie between 0 and 255."

Quantize::usage=
	"Quantize[pic,percent] quantizes the matrix 'pic'. The matrix entries are rounded to the nearest integer, then all matrix entries with absolute value less than the cutoff value are replaced by zero. The cutoff value is calculated so that the new matrix conserves 'percent' (between 0 and 1) of the energy of the original matrix."

ShannonEntropy::usage=
	"ShannonEntropy[data] computes the entropy of 'data'. 
	'data' can be a vector, matrix or other array. All entries are rounded to the nearest integer before computation of the entropy."

Taylor::usage=
	"Taylor[f (x),x,center,deg] computes the Taylor polynomial of the function f (x) with the specified center and degree." 

TaylorPlot::usage=
	"TaylorPlot[f (x),x,center,rad,degree] plots the Taylor polynomial of the function f (x) with the specified center and 
	degree on the interval of radius 'rad' centered at 'center'."

TransformPlot::usage=
	"TransformPlot[data,nr] displays a transformed matrix (array) as a grayscale image. 'nr' indicates the number of iterations of the transform. 'nr' should be a non-negative integer (default is nr=0).
	You can specify the same options as for 'ArrayPlot'. Particularly useful might be 'PixelConstrained\[Rule]x' to magnify the output. 'x' must be a positive integer; its default is 1.
	The option 'Scaled->False'  (default) assumes that the matrix entries lie between 0 and 255. Values below 0 are depicted in 'Black', values above 255 are depicted in 'White'.
	'Scaled->True' scales the display such that the minimum matrix entry corresponds to 0 and the maximum matrix entry corresponds to 255."



Begin["`Private`"]

Options[CEnergy]=Options[ListPlot];
CEnergy[data_,opts___]:=Module[{dat,ce,max},dat=Round[Flatten[data]];max=Max[dat];
ce=Accumulate[Reverse[Sort[Abs[dat]]]^2/Norm[dat]^2];
ListPlot[ce,opts,Joined->True,PlotStyle->{Red,AbsoluteThickness[2]},PlotRange->{All,{0,1}},
FrameLabel->{None,"Relative Cumulative Energy"},PlotLabel->"CEnergy Plot",Frame->Automatic,Axes->None]]

Options[CFrequency]=Options[ListPlot];
CFrequency[image_,opts___]:=Module[{min,max,flat},flat=Flatten[Round[image]];
min=Min[flat];max=Max[flat];
ListPlot[Transpose[{Table[i,{i,min,max}],Accumulate[Table[Count[flat,i],{i,min,max}]]/Length[flat]}],
opts,PlotRange->{{Min[min,0],Max[max,255]},{0,1}},Joined->True,BaseStyle->{Bold,12},
PlotStyle->{Blue,AbsoluteThickness[2]},FrameLabel->{None,"Relative Cumulative Frequency"},
PlotLabel->"CFrequency Plot",ImageSize->600,Frame->Automatic,Axes->None]]

CDF1[n_]:=Module[{A1,dim=n,h={0.037828455506995484`,-0.02384946501938002`,-0.11062440441842353`,0.37740285561265374`,0.8526986790094037`,0.37740285561265374`,-0.11062440441842353`,-0.02384946501938002`,0.037828455506995484`},
	j={-0.06453888262893842`,-0.04068941760955841`,0.4180922732222122`,0.7884856164056644`,0.4180922732222122`,-0.04068941760955841`,-0.06453888262893842`}},
	A1=SparseArray[Flatten[{Table[Band[{1,k},{dim/2,2dim},{1,2}]->h[[10-k]],{k,1,9}],Table[Band[{dim/2+1,k+2},{dim,2dim},{1,2}]->(-1)^k j[[k]],{k,1,7}]}],{dim,2dim}];
	Transpose[Apply[Plus,Partition[Transpose[Normal[A1]],dim]]]]

CDF2[n_]:=Module[{A2,dim=n,h={0.037828455506995484`,-0.02384946501938002`,-0.11062440441842353`,0.37740285561265374`,0.8526986790094037`,0.37740285561265374`,-0.11062440441842353`,-0.02384946501938002`,0.037828455506995484`},
	j={-0.06453888262893842`,-0.04068941760955841`,0.4180922732222122`,0.7884856164056644`,0.4180922732222122`,-0.04068941760955841`,-0.06453888262893842`}},
	A2=SparseArray[Flatten[{Table[Band[{1,k+1},{dim/2,2dim},{1,2}]->j[[8-k]],{k,1,7}],Table[Band[{dim/2+1,k+1},{dim,2dim},{1,2}]-> (-1)^(k+1) h[[k]],{k,1,9}]}],{dim,2dim}];
	Transpose[Apply[Plus,Partition[Transpose[Normal[A2]],dim]]]]

CDFWLT[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];]

ColorChannels[pic_]:=Module[{},{Map[({#[[1]],0,0}&),pic,{2}],
Map[({0,#[[2]],0}&),pic,{2}],Map[({0,0,#[[3]]}&),pic,{2}]}]

D4exact[n_]:=Module[{h0=(1+Sqrt[3])/(4Sqrt[2]),h1=(3+Sqrt[3])/(4Sqrt[2]),h2=(3-Sqrt[3])/(4Sqrt[2]),h3=(1-Sqrt[3])/(4Sqrt[2])},If[EvenQ[n]&&n>=4,
	Join[SparseArray[{{i_,j_}/;j==2i-1->h3,{i_,j_}/;j==2i->h2,{i_,j_}/;j==2i+1->h1,{i_,j_}/;j==2i+2->h0,{n/2,1}->h1,{n/2,2}->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,{n/2,1}->-h2,{n/2,2}->h3},{n/2,n}]],
	Print["The dimension of the matrix must be a positive even integer \[GreaterEqual] 4."]]]

D4[n_]:=Module[{h0=(1+Sqrt[3])/(4Sqrt[2.]),h1=(3+Sqrt[3])/(4Sqrt[2.]),h2=(3-Sqrt[3])/(4Sqrt[2.]),h3=(1-Sqrt[3])/(4Sqrt[2.])},If[EvenQ[n]&&n>=4,
	Join[SparseArray[{{i_,j_}/;j==2i-1->h3,{i_,j_}/;j==2i->h2,{i_,j_}/;j==2i+1->h1,{i_,j_}/;j==2i+2->h0,{n/2,1}->h1,{n/2,2}->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,{n/2,1}->-h2,{n/2,2}->h3},{n/2,n}]],
	Print["The dimension of the D4 matrix must be a positive even integer \[GreaterEqual] 4."]]]

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

D6exact[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,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}]],
	Print["The dimension of the D6 matrix must be a positive  even integer \[GreaterEqual] 6."]]]

D6[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,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}]],
	Print["The dimension of the D6 matrix must be a positive  even integer \[GreaterEqual] 6."]]]

D6Threshold[pic_,\[Alpha]_:1]:=
Module[{pic1,pic2,\[Tau],\[Tau]0,r,c},{r,c}=Dimensions[pic];pic1=D6WLT[pic,1];pic1[[1;;r/2,1;;c/2]]=Table[0,{i,1,r/2},{j,1,c/2}];
pic2=Drop[Sort[Flatten[Abs[pic1]]],{1,r c/4}];\[Tau]=(Max[pic2]+Min[pic2])/2;\[Tau]0=0;
While[Abs[\[Tau]-\[Tau]0]>\[Alpha],\[Tau]0=\[Tau];s1=Select[pic2,(#<\[Tau]&)];s2=Select[pic2,(#>=\[Tau]&)];
\[Tau]=(Mean[s1]+Mean[s2])/2];Print["Threshold value \[Tau]=",\[Tau]];Transpose[D6[r]].Map[(If[Abs[#]<\[Tau],0,#]&),pic1,{2}].D6[c]]

D6WLT[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]]=D6[r/2^k].pc[[1;;r/2^k,1;;c/2^k]].Transpose[D6[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[EnergySpectrum]=Options[ListPlot];
EnergySpectrum[f_,t_,n_:6,Options___]:=Module[{TE,dat},TE=1/\[Pi] Quiet[NIntegrate[f^2,{t,-\[Pi],\[Pi]}]];ListPlot[Transpose[{Table[i,{i,0,n}],
	dat=100/TE Prepend[Table[1/\[Pi] Quiet[NIntegrate[f Cos[k t],{t,-\[Pi],\[Pi]}]],{k,1,Floor[n]}]^2+Table[1/\[Pi] Quiet[NIntegrate[f Sin[k t],{t,-\[Pi],\[Pi]}]],
	{k,1,Floor[n]}]^2,(1/(Sqrt[2]\[Pi]) Quiet[NIntegrate[f,{t,-\[Pi],\[Pi]}]])^2]}],Filling->Axis,FillingStyle->Orange,PlotRange->All,
	PlotLabel->StringForm["`1` % of the total energy is contained in the first `2` harmonics.",
	NumberForm[Total[dat],{5,2}],n],AxesLabel->{"frequency", "%"},ImageSize->400,Ticks->{Table[i,{i,0,n}],Automatic},PlotStyle->{Red,Thick,PointSize[Large]},Options]]

Options[FourierFilterPlot]=Options[Plot];
FourierFilterPlot[coeff_,options___]:=Module[{f},f=coeff.Table[Exp[ k I t],{k,0,Length[coeff]-1}];
	Plot[Abs[f],{t, 0,\[Pi]},Ticks->{Table[k \[Pi]/4,{k,0,4}],Automatic},PlotStyle->Thick,options]]

FourierS[f_,t_,n_:6]:=Module[{ser},Print[ser=Table[1/(2\[Pi]) Integrate[f Exp[-I k t],{t,-\[Pi],\[Pi]}],{k,-Floor[n],Floor[n]}].Table[Exp[I k t],{k,-Floor[n],Floor[n]}],"\n= ",ComplexExpand[ser]]]

Options[FourierPlot]=Options[Plot];
FourierPlot[f_,t_,n_:6,Options___]:=Module[{g=Table[1/(2\[Pi]) Quiet[N[Integrate[f Exp[-I k t],{t,-\[Pi],\[Pi]}]]],{k,-Floor[n],Floor[n]}].Table[Exp[I k t],{k,-Floor[n],Floor[n]}]},
	Plot[Evaluate[{f,g}],{t,-\[Pi],\[Pi]},PlotStyle->{{Gray,Thick,Dashed},{Thick,Red}},PlotLabel->Floor[n],Ticks->{Table[j \[Pi]/4,{j,-4,4}],Automatic},Options]]

FourierSC[f_,t_,n_:6]:=Module[{},1/(2\[Pi])Integrate[f,{t,-\[Pi],\[Pi]}]+Table[1/\[Pi] Integrate[f Cos[k t],{t,-\[Pi],\[Pi]}],{k,1,Floor[n]}].Table[Cos[k t],{k,1,Floor[n]}]+Table[1/\[Pi] Integrate[f Sin[k t],{t,-\[Pi],\[Pi]}],{k,1,Floor[n]}].Table[Sin[k t],{k,1,Floor[n]}]]

Options[FourierSCPlot]=Options[Plot];
FourierSCPlot[f_,t_,n_:6,Options___]:=Module[{g=1/(2\[Pi])Quiet[N[Integrate[f,{t,-\[Pi],\[Pi]}]]]+Table[1/\[Pi] Quiet[N[Integrate[f Cos[k t],{t,-\[Pi],\[Pi]}]]],{k,1,Floor[n]}].Table[Cos[k t],{k,1,Floor[n]}]+Table[1/\[Pi] Quiet[N[Integrate[f Sin[k t],{t,-\[Pi],\[Pi]}]]],{k,1,Floor[n]}].Table[Sin[k t],{k,1,Floor[n]}]},
	Plot[Evaluate[{f,g}],{t,-\[Pi],\[Pi]},PlotStyle->{{Gray,Thick,Dashed},{Thick,Green}},PlotLabel->Floor[n],Options,Ticks->{Table[j \[Pi]/4,{j,-4,4}],Automatic}]]

Options[FrequencyHistogram]=Options[Histogram];
FrequencyHistogram[image_,opt___]:=Module[{Z},If[Length[Dimensions[image]]==3,
GraphicsArray[{Table[Histogram[Z=Round[Flatten[image[[All,All,i]]]],{1},PlotRange->{{Min[Min[Z],0],Max[Max[Z],255]},All},
ChartStyle->RGBColor[KroneckerDelta[i,1],KroneckerDelta[i,2],KroneckerDelta[i,3]],ChartBaseStyle->EdgeForm[None],
BaseStyle->{Bold,12},AxesOrigin->{0,0},AxesLabel->{"","Frequency"}],{i,1,3}]},ImageSize->900,opt],
If[Length[Dimensions[image]]<=2,Histogram[Z=Round[Flatten[image]],{1},PlotRange->{{Min[Min[Z],0],Max[Max[Z],255]},All},
ChartStyle->Gray,ChartBaseStyle->EdgeForm[None],ImageSize->400,BaseStyle->{Bold,14},AxesOrigin->{0,0},
AxesLabel->{"","Frequency"},opt],Print["The input has to be a vector, a matrix, or a 3-dimensional array."]]]]

GrayChannels[pic_]:=Module[{},Transpose[pic,{2,3,1}]]

Haar[n_]:=Module[{},If[EvenQ[n]&&n>=2,1/Sqrt[2.]Join[SparseArray[{{i_,j_}/;j==2i-1->1,{i_,j_}/;j==2i->1},{n/2,n}],SparseArray[{{i_,j_}/;j==2i-1->1,{i_,j_}/;j==2i->-1},{n/2,n}]],
	Print["The dimension of the Haar matrix must be a positive even integer."]]]

HaarThreshold[pic_,\[Alpha]_:1]:=
Module[{pic1,pic2,\[Tau],\[Tau]0,r,c},{r,c}=Dimensions[pic];pic1=HaarWLT[pic,1];pic1[[1;;r/2,1;;c/2]]=Table[0,{i,1,r/2},{j,1,c/2}];
pic2=Drop[Sort[Flatten[Abs[pic1]]],{1,r c/4}];\[Tau]=(Max[pic2]+Min[pic2])/2;\[Tau]0=0;
While[Abs[\[Tau]-\[Tau]0]>\[Alpha],\[Tau]0=\[Tau];s1=Select[pic2,(#<\[Tau]&)];s2=Select[pic2,(#>=\[Tau]&)];
\[Tau]=(Mean[s1]+Mean[s2])/2];Print["Threshold value \[Tau]=",\[Tau]];Transpose[Haar[r]].Map[(If[Abs[#]<\[Tau],0,#]&),pic1,{2}].Haar[c]]

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

Huffman[lst_]:=Module[{list,tab,tot,new},
	list=Flatten[lst];tab=Sort[Map[({Count[list,#],ToString[#]}&),Union[list]]];
	While[Length[tab]>1,new={tab[[1,1]]+tab[[2,1]],tab[[1]],tab[[2]]};tab=Append[Drop[tab,2],new];
	tab=Sort[tab,#1[[1]]<#2[[1]]&];Length[tab]];
	Print[
	"Bits per pixel with Huffman coding: ",NumberForm[tot=Total[Map[(Length[Drop[Drop[Flatten[Position[tab,ToString[#]]],1],-1]-2]Count[list,#]&),
	Union[list]]]/N[Length[list]]],".\nEntropy: ",ShannonEntropy[lst],"."];
	Print["In the table below, the first column lists the original characters; the second column their corresponding Huffman code; the last column the character frequency"];
	Sort[Map[({ToString[#],StringJoin[Map[ToString[#]&,Drop[Drop[Flatten[Position[tab,ToString[#]]],1],-1]-2]],
		Count[list,#]}&),Union[list]],#1[[3]]>#2[[3]]&]]

Options[ImagePlot] =Union[{Scaled->True},Options[ArrayPlot]];
ImagePlot[A_,opts___]:=Module[{opt,sc},
	opt=Complement[{opts},{Scaled->True},{Scaled->False}];
	sc=Scaled/.{opts}/.Options[ImagePlot];
	If[sc,ArrayPlot[1.-Rescale[A],opt,PixelConstrained->1],
	ArrayPlot[255-Map[(Min[255,Max[0,Round[#]]]&),A,{Length[Dimensions[A]]}],
	opt,PlotRange->{0,255},PixelConstrained->1]]]

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

JoinGrayChannels[rgray_,ggray_,bgray_]:=Module[{},Map[({#,0,0}&),rgray,{2}]+Map[({0,#,0}&),ggray,{2}]+Map[({0,0,#}&),bgray,{2}]]
	
NTSC[pic_]:=Module[{},Floor[0.299 pic[[All,All,1]]+0.587 pic[[All,All,2]]+0.114 pic[[All,All,3]]]]

Options[PicPlot] =Options[ArrayPlot];
PicPlot[A_,pc_:1,opt___]:=Module[{},ArrayPlot[255-A,PlotRange->{0,255},
PixelConstrained->pc,opt]]

PSNR[data1_,data2_]:=Module[{d=Depth[data1]-3,l1=Dimensions[data1][[1]],l2=Dimensions[data1][[2]],err,range},
	PSNR::range="Warning: Maximum of matrix entries exceeds 255, 
	or minimum of matrix entries is below 0.";
	If[Max[{data1,data2}]>255||Min[{data1,data2}]<0,Message[PSNR::range],Null];
	err=1/(l1 l2)/(1+2d) Total[Flatten[(data1-data2)^2]];10. Log[10,255.^2/err]]

Quantize[pic_,pct_]:=
Module[{dat,acc,ce,k=1,cutoff,rd},rd=Round[pic];dat=Reverse[Sort[Round[Flatten[Abs[rd]]]]];
	acc=Accumulate[dat^2/Norm[dat]^2]//N;ce[k_]:=acc[[k]];
	While[ce[k]<pct,k++] ;
	Print["Quantization: At least ",NumberForm[k,DigitBlock->3], " matrix elements with the highest energy (out of ",NumberForm[Length[dat],DigitBlock->3], " elements total) are preserved.
	The rest is set to 0 (all entries with absolute value below ",cutoff=dat[[k]],")."];
	Return[rd/.Table[i->0,{i,-cutoff+1,cutoff-1}]];]

ShannonEntropy[data_]:=N[Entropy[2,Round[Flatten[data]]]]

Taylor[fun_,var_,center_,deg_]:=Module[{f},f=Function[var,fun];
	Table[Derivative[k][f][center]/k!,{k,0,deg}].Table[(var-center)^k,{k,0,deg}]]

Options[TaylorPlot] =Options[Plot];
TaylorPlot[fun_,var_,center_,rad_,deg_,opt___]:=Module[{f,tp,x},f=Function[var,fun];
	tp=N[Table[Derivative[k][f][center]/k!,{k,0,deg}].Table[(x-center)^k,{k,0,deg}]];
	Plot[{f[x],tp},{x,center-rad,center+rad},opt,PlotStyle->{Red,{Blue,Thick}},
	Epilog->{AbsolutePointSize[10],Blue,Point[{center,f[center]}]}]]

Options[TransformPlot] =Union[{Scaled->False},Options[ArrayPlot]];
TransformPlot[A_,n_:0,opts___]:=Module[{opt,k,sc,r=Dimensions[A][[1]],c=Dimensions[A][[2]]},
	k=Ceiling[Abs[n]];
	opt=Complement[{opts},{Scaled->False},{Scaled->True}];
	sc=Scaled/.{opts}/.Options[TransformPlot];
	If[sc,ArrayPlot[1.-Rescale[A],opt,PixelConstrained->1,
		Epilog->Table[{Red,Line[{{c/2^j+1,r},{c/2^j+1,r-r/2^(j-1)}}],Line[{{1,r-r/2^j+1},{c/2^(j-1),r-r/2^j+1}}]},
		{j,1,k}]],
	ArrayPlot[255-Map[(Min[255,Max[0,Round[#]]]&),A,{Length[Dimensions[A]]}],
	opt,PlotRange->{0,255},PixelConstrained->1,
	Epilog->Table[{Yellow,Line[{{c/2^j+1,r},{c/2^j+1,r-r/2^(j-1)}}],Line[{{1,r-r/2^j+1},{c/2^(j-1),r-r/2^j+1}}]},
		{j,1,k}]]]]

End[ ]
EndPackage[ ]