%% Haar Edge Detection
% David Ruch and Patrick J. Van Fleet
% Minicourse #4, January 2008 Joint Mathematics Meetings
% San Diego, CA

%% Objective
% In this m-file, we will learn how to perform naive image edge 
% detection using the 2D 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{15});

%% 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 fifteenth image in the list.

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

%% Edge Dectection
% To perform edge detection on the image, we first compute the 2D HWT.

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

%% 
% The next step is to replace the upper left hand corner with a zero 
% matrix of the same dimensions.  The DiscreteWavelets module PutCorner 
% will help here.  PutCorner takes two arguments and replaces the upper 
% left hand corner of the first argument with the second argument.

[r c]=size(A);
Z = zeros(r/2,c/2);
newB = PutCorner(B, Z);
close;
WaveletDensityPlot(newB,1,'DivideLinesThickness',2);

%%
% We are left with the vertical, horizontal, and diagonal differences.  We
% compute the inverse of the modified transform to obtain the edges.

edges = IHWT2D(newB,1);
close;
ImagePlot(edges);

%%
% Sometimes it is easier to look at the image negative to identify the 
% edges.

close;
ImagePlot(255-edges);

%%
% Try the procedure again with different images.

%% Iterating the Process
% You can use the iterated transform as part of the algorithm.  

its = 3; %Do not exceed the maximum number of iterations allowed!! 

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

%%
% The next step is to replace the upper left hand corner with a zero 
% matrix of the same dimensions.  The DiscreteWavelets module PutCorner 
% will help here.  PutCorner takes two arguments and replaces the upper 
% left hand corner of the first argument with the second argument.

Z = zeros(r/2^its,c/2^its);
newB = PutCorner(B, Z);
close;
WaveletDensityPlot(newB,its,'DivideLinesThickness',[2 2 2]);

%%
% Compute the inverse transformation to find the edges.

edges = IHWT2D(newB,its);
close;
ImagePlot(edges);

%% 
% Sometimes it is easier to look at the image negative to identify the
% edges.

close;
ImagePlot(255-edges);

%% To Consider
%
% *What happens to the edge image if you increase the number of iterations?  
%
% *Why do you suppose this happens?
%
% *How could you modify the process to improve the identification of edges?
%

%%
close all;
displayEndOfDemoMessage(mfilename)