%% Measures: Cumulative Energy (CE) and Peak Signal to Noise Ratio (PSNR)
%
% *Wavelets Workshop*
%
% June 2 - 5, 2009
%
% University of South Florida
%
% *Catherine Beneteau*, *Caroline Haddad*, *David Ruch*, *Patrick Van Fleet*

%% Objectives
% The purpose of this notebook is to give you a brief introduction to the
% DiscreteWavelets Toolbox and show you how to use it to load 
% images.  You will learn how to use measures and tools such as cumulative
% energy and PSNR.

%% Help on the DiscreteWavelets Toolbox
%
% Help for the toolbox is available by clicking on Help and then Product
% Help (or press F1) and then clicking on the DiscreteWavelets Toolbox.
% Several demos and examples are available as well by clicking on the Demos
% tab on the Help menu.

%% Image Basics
%
% The DiscreteWavelets Toolbox comes with 18 grayscale images and 9 color
% images for you to use.  There are three functions available to tell you more about these images. 
%
% The first function is called |ImageList|.  This function can tell you the
% names and sizes of the digital images in the Toolbox.

ImageList('ImageType','GrayScale')
ImageList('ImageType','Color')
ImageList()

%%
%
% The second function is called |ImageNames|.  This function returns a list
% of absolute path names to each image.  Since the Toolbox can be installed
% in multiple locations, this functions takes the work out of locating the 
% images.

color=ImageNames('ImageType','Color')

%%
%
% You'll notice a list of 9 path/file names.  If you want to use 
% building.png in an application (the second name in the list), you simply 
% use |color{2}|.

color{2}

%%
% From |ImageList|, you probably noticed that the images were quite large.
% Sometimes it's handy to work with smaller images.  If you set the 
% directive |ListThumbnails| to |True|, then the file names returned are 
% those of the thumbnails.

color=ImageNames('ImageType','Color','ListThumbnails','True');
color{2}

%%
%
% If you wish to view thumbnails of all images included in the Toolbox, use
% the |ShowThumbnails| command.  You can specify the |ImageType| as either
% |GrayScale| or |Color|.  If no image type is given the function shows 
% thumbnails of all included images.

ShowThumbnails('ImageType','GrayScale')
ShowThumbnails('ImageType','Color')

%% Loading a Digital Image
%
% DiscreteWavelets comes with a command called |ImageRead|.  This command
% essentially needs one argument - the location, either on a computer or 
% on the internet, of a digital image.  In the case of a grayscale image, 
% |ImageRead| returns a matrix (named A in the commands below) containing 
% the gray scale intensity values for each pixel.  
%
% Here are some examples.

gray = ImageNames('ImageType','GrayScale');
A = ImageRead(gray{1});
[rows cols]=size(A);
str=sprintf('The dimensions of A are %i x %i.',rows,cols);
disp(str);

%%
%
% You can read in the thumbnail images too.
%

gray = ImageNames('ImageType','GrayScale','ListThumbnails','True');
A = ImageRead(gray{1});
[rows cols]=size(A);
str=sprintf('The dimensions of A are %i x %i.',rows,cols);
disp(str);
%% Cumulative Energy
%
% One of the important skills students develop in this course is the
% ability to write modules or functions.  The cumulative energy function is 
% an excellent starting point for this development. 
%
% To compute the cumulative energy of vector v, we sort the absolute values
% of the elements largest to smallest and then square the components of the
% resulting vector.  Call this vector y.  We find the kth element of the 
% cumulative energy vector by summing the first k elements of y and 
% dividing the result by the norm of v squared.
%
% The first two steps of this process are easy.  Let's use a small vector
% as an example.

v = [-3 -2 0 0 5 -8 3 1 1 2]
y = fliplr(sort(abs(v)))
y = y.^2

%%
% Next comes the cumulative sums.  We could write a loop for k = 1 through
% 10 and an inner loop j = 1 to k to perform the task, but we instead 
% encourage our students to take advantage of built-in functions that 
% accomplish tasks.  In most cases, these built-in functions are much 
% faster than attempting to extract individual elements from a vector or 
% list.  For cumulative sums, the built-in command |cumsum| is well-suited
% for our needs. 

x = cumsum(y)
ce = x/(v*v')

%%
%
% We can plot the cumulative energy using the Matlab command |plot|.
% Note the x-axis is the element number of ce.

plot(ce);

%% The Cumulative Energy Function
%
% In Matlab, it is very easy to write modules or functions.  In the My
% Documents folder, you will find a copy of the DiscreteWavelets Toolbox -
% the folder is called DiscreteWavelets.  In this folder, you should find a
% folder called MyFiles.  In this folder, open the m-file called MyCE.m
%
% Using this file, we will write a function for computing the cumulative
% energy of a vector (or matrix!) together.

%%
% 
% The cell below creates a vector of 20 random integers in the range 
% 0,...,10, and then calls the |MyCE| function and the DiscreteWavelets 
% command |CE| to compute the cumulative energy.  Don't execute the cell 
% until you have completed and executed the cell above.

v=round(rand(1,20)*10)
MyCE(v)
CE(v)

%% Exercises
%
%%
% *Exercise 1*
%
% The |MyCE| function above will not work on matrices.  The function is 
% expecting a vector.  How can you modify the |MyCE| function so that it will
% work on matrices?  (Hint: Check the command |reshape| under Help.)

%%
% Place your answer here.  To type an answer, start a line with %.  Matlab
% commands are entered as usual.

%%
% *Exercise 2*
%
% Load the grayscale image of the chess pieces from the DiscreteWavelets
% package, find its cumulative energy using the |MyCE| function, and then 
% plot the cumulative energy of the image.  How many elements of the image 
% constitute 90% of the energy?  How many zeros are in the image?

%%
% Place your answer here.  To type an answer, start a line with %.  Matlab
% commands are entered as usual.

%%
% *Exercise 3*
%
% Describe the plot of the cumulative energy of a (nonzero) constant
% vector.  (Hint: If you wish, you can generate a constant vector using the
% |ones| command - see Help for more information.)

%%
% Place your answer here.  To type an answer, start a line with %.  Matlab
% commands are entered as usual.
%% Peak Signal-To-Noise Ratio
%
% The DiscreteWavelets command for peak signal-to-noise ratio is |PSNR|.  
% The input values are two matrices of equal size.

img = ImageNames('ImageType','GrayScale','ListThumbnails','True');
A = ImageRead(img{4});
[rows cols] = size(A);

A1 = A + round(rand(rows,cols)*20 - 10);
A2 = A + round(rand(rows,cols)*80 - 40);

ImagePlot(A);
figure;
ImagePlot(A1);
figure;
ImagePlot(A2);

psnr1=PSNR(A,A1);
psnr2=PSNR(A,A2);

str=sprintf('The PSNR of A and A1 is %f and the PSNR of A and A2 is %f.',psnr1,psnr2);
disp(str);

%% Exercises
%

%%
% *Exercise 1*
%
% Write a function to compute the PSNR of two matrices.  Use the m-file
% MyPSNR that can be found in the MyFiles folder.  Useful Matlab commands
% are |sum| and |log10|.
%
% Here is some code to test your module :

mypsnr1 = MyPSNR(A, A1)
mypsnr2 = MyPSNR(A, A2)