%% Digital Images and Measures
%
% *Wavelets Workshop*
%
% June 4-7, 2008
%
% University of St. Thomas
%
% *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.  Some basic image manipulation is illustrated as well.  You will
% also learn how to use measures and tools such as cumulative energy, 
% entropy, PSNR, and Huffman coding.

%% 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);

%%
% Here is the first row of A

A(1,:)

%%
%
% There are two additional directives for |ImageRead|.  The first is called
% |PrintInfo| and if set to |True|, will give you the dimensions and type 
% (grayscale, color) of image you have read.  The |PowersOfTwo| option, if 
% set (say to K) will chop off enough rows on the bottom and columns on the
% left to make the dimensions of the image divisible by 2^K.  This is 
% really handy when students pull their own images off the web and we need 
% the dimensions to be divisible by 2^K.

[gray  color] = ImageNames('ImageType','All','ListThumbnails','True');
A = ImageRead(gray{2},'PrintInfo','True');
B = ImageRead(color{3},'PrintInfo','True');

A = ImageRead(gray{3},'PrintInfo','True','PowersOfTwo',6);
[rows cols] = size(A);
str=sprintf('The new dimensions of A are %i x %i.',rows,cols);
disp(str);

%%
% Note that B was processed as a color image.  This means that B returned 
% three matrices - the first matrix is the red portion of the image (with 
% values 0 to 255), the second matrix is the green portion of the image, 
% and the third matrix is the blue portion of the image.
%
% Notation-wise, we have the matrices B(:,:,1) (red), B(:,:,2) (green), 
% and B(:,:,3) (blue), but it might be easier to make the call using 
% |Split3D|:

B=ImageRead(color{3},'PrintInfo','True','PowersOfTwo',6);
[red, green, blue]=Split3D(B);
size(red)

%% Plotting a Digital Image
%
% It is easy to plot digital images read with |ImageRead| in Matlab.
% The command is called |ImagePlot|. 
%
% Once you have an image loaded via ImageRead, here's what you can do:

ImagePlot(A);
B=Make3D(red,green,blue);
ImagePlot(B);

%%
% You can plot individual channels of a color image if you like.  It will
% be grayscale by default, but if you want to see it in its original color,
% simply add the directive |ChannelColor|:

ImagePlot(red);
figure;
ImagePlot(red,'ChannelColor',[1 0 0]); %Plot the red channel
figure;
ImagePlot(green,'ChannelColor',[0 1 0]); %Plot the green channel
figure;
ImagePlot(B);

%% Reading and Plotting Images from the Internet
%
% It is quite easy to read and plot images from the internet.  All you need
% is a connection to the internet and a url for the image.  You can find 
% the url for an image in Internet Explorer by simply right-clicking on it
% and then clicking on Properties:

A=ImageRead('http://math.usf.edu/images/people/faculty/beneteau.jpg','PrintInfo','True','PowersOfTwo',2);
ImagePlot(A); %Plot the image

%%
% *Warning*:  Most of our applications will require the dimensions of the
% image to be divisible by some power of 2.  Use the |PowersOfTwo| 
% directive to appropriately resize internet images.  Images that appear 
% as grayscale are often stored as color images on the internet.  Use the 
% |PrintInfo| directive to check this before using internet images in 
% applications.

%% Manipulating Images
%
% Once you have an image loaded and know how to plot it, you can do all 
% kinds of things with it.  
%
% Convert color to grayscale:

color=ImageNames('ImageType','Color','ListThumbnails','True');
A=ImageRead(color{3});
ImagePlot(A);
figure;
[red, green, blue]=Split3D(A);
gray=(red+green+blue)/3;
ImagePlot(gray);

%%
% Darken the image

ImagePlot(A/2);

%% 
% Image negative

ImagePlot(255-gray);

%%
% Color negative

ImagePlot(255-A);

%%
% Flip it upside down

ImagePlot(flipud(gray));

%%
% What does the transpose do?

ImagePlot(gray');

%%
% The previous cell is basically a rotation of the image counterclockwise
% 90 degrees.   In the cell below, can you write code to rotate the image 
% clockwise 90 degrees?

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

%% Converting Color Images to YCbCr Space
%
% DiscreteWavelets contains a command for converting from RGB color space
% to YCbCr space.  The command is called |RGBToYCbCr|.  Here is an example
% call:

color=ImageNames('ImageType','Color'); %Get names of all color images included with the toolbox
A=ImageRead(color{5}); %Read a color image
ImagePlot(A); %Plot the image
[R,G,B]=Split3D(A); %Split A into the three color channels
figure
ImagePlot(R,'ChannelColor',[1 0 0]); %Plot the red channel
figure
ImagePlot(G,'ChannelColor',[0 1 0]); %Plot the green channel
figure
ImagePlot(B,'ChannelColor',[0 0 1]); %Plot the blue channel

%%
% Now convert to YCbCr space and plot the resulting images

B=RGBToYCbCr(A,'DisplayMode','True'); %Do the color space conversion
[Y,Cb,Cr]=Split3D(B); %Split B into the three channels
figure
ImagePlot(Y,'Title','Y Channel'); %Plot the Y channel
figure
ImagePlot(Cb,'Title','Cr Channel'); %Plot the Cr channel
figure
ImagePlot(Cr,'Title','Cb Channel'); %Plot the Cb channel

%%
% Note that we have used the directive |DisplayMode| set to |True|.  This
% maps the values of the Y, Cb, Cr channels from [0,1], [-1/2,1/2],
% [1/2,1/2] to integer intervals in the range [0,255].

%%
% You can convert back to RGB space using the command |YCbCrToRGB|:

A=YCbCrToRGB(B,'DisplayMode','True'); %Do the color space conversion
[R,G,B]=Split3D(A); %Split A into the three channels
figure
ImagePlot(R,'ChannelColor',[1 0 0]); %Plot the red channel
figure
ImagePlot(G,'ChannelColor',[0 0 1]); %Plot the green channel
figure
ImagePlot(R,'ChannelColor',[0 0 1]); %Plot the blue channel
figure
ImagePlot(A);

%% Exercises
%
% 
%% 
% *Exercise 1* 
%
% In the cell below, find and plot the grayscale image of the clown that
% comes with the DiscreteWavelets Toolbox.

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

%%
% *Exercise 2*
% In the cell below, find and plot the red, green, and blue color channels
%(when reading the image, name the three matrices red, green, and blue) of
% the image of the goats that comes with the DiscreteWavelets Toolbox.

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

%%
% *Exercise 3*
%
% In an example above, we converted a color image to grayscale by simply
% averaging the channels.  The National Television System Committee 
% suggests using a weight average of these three channels for this 
% conversion.  The NTSC map is
%
% gray = .299*red + .587*green +  .114*blue
%
% In the cell below, load and plot the legos image that comes with the
% package.   Next create two grayscale images from this color image. Build
% the first by averaging the three channels and construct the second using
% the NTSC map.  Call the first image gray1 and the second gray2.  Plot 
% both images - do you notice a difference? 
%
% Repeat the exercise with some other color images included with the
% package or one from the internet.

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

%%
% *Exercise 4*
%
% Gamma correction is a simple image processing tool that can be used to
% lighten or darken an image.  The process is quite simple.  Suppose A is a
% grayscale image matrix.  We first divide all elements of A by 255 to 
% obtain values in the interval [0, 1] and then raise each value to some 
% positive exponent |gamma|.  We then multiply each element by 255 and 
% round the result to the nearest integer.
%
% In Matlab we can do "incorrect" mathematical operations.  If A is
% our image matrix, then A/255 makes some sense, but (A/255).^gamma 
% actually raises each element of A to the gamma power.  In the cell 
% below, using the |round| function, write code that will perform gamma 
% correction on the image from Exercise 1.  Call the output B.  Try several
% values of |gamma| (I would suggest naming a variable gamma and then 
% changing its values) and plotting the result.  
%
% Can you characterize what different values of gamma do to the image?
%
% Feel free to load other grayscale images in the package and perform gamma
% correction on them.

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

%%
% *Exercise 5*
%
% In an example above we computed the negative of a color image.  In the
% original image, each pixel was represented by a linear combination of 
% vectors representing red, green, and blue.  For the image negative, what
% vectors are used to form the linear combination for each pixel?  What are
% the associated colors?

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


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

%% Entropy
%
% The DiscreteWavelets command for entropy is |Entropy|.  The function 
% takes either a vector or matrix as input.  Here are some example calls:

v = 1:16
str=sprintf('The entropy of vector v is %f.', Entropy(v));
disp(str);
A = eye(8);
str=sprintf('The entropy of the 8x8 identity matrix is %f.', Entropy(A));
disp(str);

%% Exercises
%

%% 
% *Exercise 1 (Challenge)*
%
% Under the MyFiles folder, open the m-file MyEntropy.m.  In this file,
% create a function that computes the entropy of either a vector or a
% matrix.  Useful Matlab commands are |histc|, |uniq|, |length|, and
% |log2|.
%
% Use the cell below to test your code.

v = round(rand(1,10)*10);
Entropy(v)
MyEntropy(v)

A = round(rand(8,8)*10);
Entropy(A)
MyEntropy(A)

%% 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)

%% Huffman Codes
%
% The DiscreteWavelets Toolbox includes two functions that are useful for
% creating and analyzing Huffman codes.
%

%% MakeHuffmanCodes
%
% The DiscreteWavelets function |MakeHuffmanCodes| can be used to generate 
% Huffman codes for a list of integers, a string, or a matrix of integers. 
% Note that integer input must be nonnegative.
%
% The routine returns five pieces of information.  The first is a list of 
% unique integers or characters that appeared in the input, the second is
% list of the relative frequencies of each unique value, and the third is 
% a list of Huffman codes for each unique value.  The fourth output is the
% original bitstream length of the input and the last return value is the
% new bitstream length.
%

[uniq, freq, codes, origlen, newlen] = MakeHuffmanCodes('pitterpatter')

%%
% We can also find the Huffman codes of an image:

img = ImageNames('ImageType', 'GrayScale', 'ListThumbnails', 'True');
A = ImageRead(img{1});
ImagePlot(A);
[uniq, freq, codes, origlen, newlen] = MakeHuffmanCodes(A);
str=sprintf('The original bitstream length is %i and the new bistream length is %i.',origlen,newlen);
disp(str);

%% HuffmanTree
%
% For small input, we can use |HuffmanTree| to visualize the output of
% |MakeHuffmanCodes|.  The input of |HuffmanTree| is the first three 
% outputs of |MakeHuffmanCodes|.  There are several graphics options for 
% |HuffmanTree|.  See Help for more information.

[uniq, freq, codes, origlen, newlen] = MakeHuffmanCodes('tennessee');
figure;
HuffmanTree(uniq,freq,codes);

%% Exercises
%
% *Exercise 1*
%
% Load the surfer image (small version) from the DiscreteWavelets package
% and find the bits per pixel that results from the Huffman coded version
% of the image.

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

%%
%
% *Exercise 2 (Challenge)
%
% Suppose you pass a string |s| to |MakeHuffmanCodes|.  Write a function 
% that will take the codes returned by |MakeHuffmanCodes| and |s| and 
% return the new bitstream.  Call the function |MakeBitstream|.  You will 
% have to add this m-file to the MyFiles folder.
%
% For example:

s = 'boohoo';
figure;
[uniq, freq, codes, origlen, newlen] = MakeHuffmanCodes(s);
HuffmanTree(uniq, freq, codes);

%%
%
% Using the tree above and |s|, we see that the new bitstream for 
% "boohoo" is 01110011.

%%
close all;
displayEndOfDemoMessage(mfilename)