%% Haar Transform 2D
%
% *Wavelets Workshop*
%
% June 4-7, 2008
% 
% University of St. Thomas
%
% *Catherine Beneteau*, *Caroline Haddad*, *David Ruch*, *Patrick Van Fleet*

%% Objective
% In this m-file, we will explore the two-dimensional discrete Haar
% wavelet transform.

%% Conventions
% This m-file uses the DiscreteWavelets Toolbox.  Help is available on all
% functions in the toolbox by accessing the Help menu and clicking on the
% DiscreteWavelets Toolbox.  
%
% Demos are also available.  Click the Demos tab under Help to access all
% demos.

%% Available Images
% The DiscreteWavelets packages comes with 18 grayscale images.  You can
% see information about these images (name, size, etc.) by issuing the 
% command ImageList.

ImageList('ImageType','GrayScale');

%% 
% You can also get a look at these images by using the command
% ShowThumbnails.  

ShowThumbnails('ImageType','GrayScale');

%% 
% No matter where you installed the DiscreteWavelets package on your
% computer, you can retrieve the absolute path and file name for each 
% included image.  The command ImageNames produces a list of file names.

gray = ImageNames('ImageType','GrayScale');
disp(gray{1});

%% Loading and Plotting Images
% Once you have the list of file names, it is very easy to load and plot
% images.  Let's load the first image in the list.

A = ImageRead(gray{1});
close;
ImagePlot(A);

%% Using HWT2D
% In the cell below, we compute and plot the two-dimensional HWT of image
% A.  We compute the transform of N[A] to speed up the computation.

B = HWT2D(A,1);
close;
WaveletDensityPlot(B,1,'DivideLinesThickness',2);

%% Using IHWT2D
% In the cell below, we compute the inverse transform of the image.

B = HWT2D(A,1);
WaveletDensityPlot(B,1,'DivideLinesThickness',2)

origA = IHWT2D(B,1);
figure;
ImagePlot(origA);

%% Iterating the Process
% Just as with the one-dimensional HWT, we can iterate the two-dimensional
% HWT.  In this case, the HWT2D is applied the approximation (blur) 
% portion of the transform.
% 
% Warning:  Make sure to know the maximum amount of iterations you can 
% compute before proceeding!!  You can use ImageList to display the 
% maximum number of iterations that can be performed on each of the 
% included images.

ImageList('ImageType','GrayScale');

%%
% Here we compute four iterations of the Haar wavelet transformation.

its = 4;
B = HWT2D(A, its);
close;
WaveletDensityPlot(B,its,'DivideLinesThickness',[2 2 2 2]);

%%
% You can plot various portions of the transform.  Here is the vertical
% portion of the third iteration.

WaveletDensityPlot(B,its,'Iteration',3,'Region','Vertical');

%%
% Here is the horizontal portion of the first iteration.

close;
WaveletDensityPlot(B,its,'Iteration',1,'Region','Horizontal');

%%
% Here is the blur.

close;
WaveletDensityPlot(B,its,'Region','Blur');

%%
% We can "blow up" the blur.

close;
WaveletDensityPlot(B,its,'Region','Blur','Magnification',8)

%% Writing Your Own Modules - MyHWT2D1
%
% Let's write a function to compute one iteration of the two - dimensional 
% Haar wavelet transformation.  We will call the routine |MyHWT2D1|.  This 
% function will use the module |MyHWT1D1| or |HWT1D1|.  In the case of the 
% former, make sure you have it loaded before you proceed.
%
% Our goal is to compute B = W A W^T.  The routine |MyHWT1D1| computes W v 
% where v is a vector.  So to compute WA, we apply |MyHWT1D1| to each column 
% of A.  
%
% You can write a loop to apply MyHWT1D1 to each column of A.  (Hint: The 
% jth column of A is A(:,j).)
%
% Once you have C = WA, all you need is B = C W^T.  But B^T = W C^T.  
% Can you modify the code to compute WA  in conjunction with MyHWT1D1 to 
% complete the task?  (Hint: B^T = W C^T.)
%
% To create your function, find the DiscreteWavelets folder and inside of
% it open the MyFiles folder.  In this folder, look for MyHWT2D1.m  Open
% this file and put your code there.  Save and exit.
%

%% Testing the Function
%
% Here is some code to test your function:

gray = ImageNames('ImageType','GrayScale','ListThumbnails','True');
A = ImageRead(gray{1});
B1 = MyHWT2D1(A);
B2 = MyHWT2D1(A);

x = max(max(abs(B2 - B1)));
if (x==0)
   disp('The function works.');
else
   disp('The function does not work.');
end;
%% Writing Your Own Modules - MyIHWT2D1
%
% Writing a function to compute for |MyIHWT2D1| uses the same ideas as 
% |MyHWT2D1|.  The only thing that is different is that we use |MyIHWT1D1|
% instead of |MyHWT1D1|.  
%
% To create your function, find the DiscreteWavelets folder and inside of
% it open the MyFiles folder.  In this folder, look for MyIHWT2D1.m  Open
% this file and put your code there.  Save and exit.
%

%% Testing the Function
%
% Here is some code to test your function:

gray = ImageNames('ImageType','GrayScale','ListThumbnails','True');
A = ImageRead(gray{1});
B = MyHWT2D1(A);
newA = MyIHWT2D1(B);

x = max(max(abs(round(A-newA))));
if (x==0)
   disp('The function works.');
else
   disp('The function does not work.');
end;
%%
close all;
displayEndOfDemoMessage(mfilename)