(Much debt is owed to David Martin for parts of this assignment)

For those wanting to do some assignments in C/C++ see the last question. Even if you prefer to work in C/C++, it is highly recomended that you explore Matlab if you are not already familiar with it. If you expect that you might join/collaborate with the vision group, doing some of the assignments in C/C++ and using the software mentioned below is highly recomended.

Matlab is installed on the UA CS linux machines (the path is /usr/local/bin). Normally, to begin, you just type "Matlab" at the shell prompt. For example, you may want to use the machines gr01 through gr08 in GS 920, either on-site, or remotely. If your DISPLAY environment is properly set, the default behavior is the GUI interface. Otherwise, you will get the the command line interface. This GUI is required for much of this assignment. You may find that the GUI is slow to start up, and painfully slow to use over a slow connection. For a faster startup (but without some of the features needed for some parts of this assignment), you can try

        matlab -nodisplay -nojvm

Matlab is available for personal use to UA faculty, staff and students. It can be downloaded from http://sitelicense.arizona.edu/matlab/ by any University of Arizona faculty, staff or student with a netid. The site-license includes MATLAB, Simulink, plus 48 toolboxes. This license is available for one year. At the end of the year term, partners will assess continuance of the project.

IMPORTANT NOTE. Even if you develop programs on some other machine, they must work on our refernece machines, gr01 through gr08.

1. Become familiar with Matlab.

   Read this short Matlab tutorial and be aware of this longer Matlab primer . The complete documentation for Matlab is also available on the web. All Matlab commands are well documented with Matlab's help system. For example, if you want to know how the svd function works, type help svd. You can also get help on all the built in operators and even the language itself. Typing help gives a list of help topics.

2. Reading and Displaying Images

   Place a color JPEG (or TIFF) image in your working directory and load it into a variable using imread. Type help imread to find out how to do this. Create a new figure using figure(1) and display the image in it using imshow.

   Tip: You may want to turn on paging (more on) for reading help pages.

   Tip: Make sure you put a semicolon at the end of your imread command. The semicolon at the end of a line prevents Matlab from displaying the result of an expression. Most of the time, you want the semicolon.

   What data structure is used to represent the image? Type whos to get a list of active variables along with their types and sizes. You can see that your image (if it is an RGB image) is a 3D array of bytes.

   What is the range of values contained in the image array? Use the min and max functions to do this.

   Hint: These commands by default operate on only one dimension of their argument. Convert your image into a 1D array (a vector) using the (:) notation when you pass it to min and max.

   Convert the image to grayscale using rgb2gray. What now is the representation of the image? Use size or whos to find this out. You will see that we now have a 2D array, or a matrix.

3. Visualizing Matrices

   Convert the grayscale image into a double-precision type, and scale the values so that they lie in the range [0,1].

   Hint: Use the double function to do the type conversion, and then divide by the maximum allowable value for a byte, which is 255.

   Create a new figure using figure(2), and use imagesc to display the grayscale image in it. Why does it look so strange?! Type colorbar to find out. You can see that 0 maps to blue, 0.5 to green, and 1 to red. This is the default jet colormap that is useful for data visualization. For this example, we might prefer the grayscale colormap. Type colormap(gray) to do this.

   Tip: Type help gray to see what default colormaps are available. Note that you can also make your own colormaps.

   Tip: The imagesc command takes a second argument that lets you specify the range of values. By default, imagesc scales the data to use the full colormap, but this is not always what you want. In this example, we could add the argument [0 1] to the imagesc command so that the gray values are displayed faithfully.

   The image is probably distorted, i.e. the pixels aren't square. Use the axis command to fix this.

   Tip: When you display matrices as images, you usually want the axes to be scaled so that the pixels are square and the image not distorted. See help axis to determine how to do this. The axis image command is a useful one. You can also turn the display of the axes on and off with the axis command.

4. Files and Paths

   The interactivity of Matlab is great for debugging and experimenting, but usually one wants to type code into a file. Create a file hw1.m and put the commands in it to read and display your image in figure(1), convert the image to grayscale, convert it to double-precision scaled to [0,1], and then display the grayscale image in figure(2) with the colorbar, grayscale colormap, and square pixels. Now type hw1 at the Matlab prompt to execute the commands in the file. A file of this sort is known as a script.

   Tip: Matlab has pwd, ls, and cd commands like the shell.

   Tip: If your script is in some other directory than pwd, then you can add that other directory to Matlab's search path with the addpath command. The current directory is in Matlab's search path by default.

5. Manipulating Matrices

   In this section, we will work with the converted image that you should have displayed in figure(2).

   Tip: Array indices in Matlab start at 1, not at 0.

   Tip: As in standard mathematical notation, the first index of a matrix is the row, and the second index the column. When viewing a matrix as an image, this means that the first index is the y direction, and the second the x direction. As is usual with image, the origin is at the top left corner, the positive x axis points right, and the positive y axis points down.

   Get the width and height of the image using size. The size function, as is common in Matlab, can return multiple arguments.

   Hint: To store both the width and height into variables in one go, try [h,w]=size(im). You could alternately do h=size(im,1) and w=size(im,2).

   Write a couple of nested for loops to set every 5th pixel, horizontally and vertically, to 1. This should set 1/25th of the pixels to white in a square lattice pattern.

   Hint: Use the colon operator to define the limits of the for loops. See help colon and help for to see how to do this. Specifically, you want the minval:interval:maxval form.

   Set all the same pixels to 0, making the ones that you just set to white become black. Do it this time without using any for loops.

   Hint: You can index arrays in Matlab with vectors as well as with scalars, so im(Y,X)=0 will set multiple entries of im to zero when either X or Y are vectors, such as the vectors returned by the colon operator.

   Set all the pixels whose values are greater than 0.9 to zero. The command to do this is im(find(im>0.9))=0. It would be good for you to understand why this works, but if this is just too confusing, move on for now.

   Hint: The (im>0.9) expression evaluates to a boolean matrix, which is true where the condition holds. The find function returns a vector containing the indices of the true values. This vector is then used to index the image and set all the values to zero. Why does this last step work when the matrix is 2D and the indexing is 1D? Matlab lets you treat matrices as 1D vectors too, linearizing the matrix in column-major order.

   Tip: Matlab makes it easy to write "vectorized" expressions without having to write for loops or if statements. For example, this will add all the values of the image:

   sum(im(:))

   The following will count the number of values greater than 0.9:

   numel(find(im>0.9))

   So will this:

   sum(sum(im>0.9))

   This will halve only those values greater than 0.9 (note the use of the .* operator to do element-size multiplication of matrices):

   im = im - 0.5*im.*(im>0.9);

   And so will this:

   mask = (im>0.9);
im = im.*~mask + im.*mask*0.5;

   This will set 100 unique random pixels to zero:

   p = randperm(numel(im));
im(p(1:100)) = 0;

   See help elmat for a list of interesting matrix manipulation and creation routines.

   Extract a rectangular region of the image into another variable. Use notation like im2 = im(colmin:colmax,rowmin:rowmax) to do this. Use the getrect function to visually select the rectangle you want from figure(2). Create a new figure with figure(3), and then display the extracted region in that figure using imagesc. Fix the axes so that the pixels are square.

   Tip: getrect returns floating-point values that you must convert to integers before using them as indices. Use the round function for this.

   Tip: In addition to getrect there are also getline and getpts functions. Also, the ginput function is an extremely useful function to know. It allows you to retrieve mouse clicks and key presses from a figure window.

6. Writing Images

   Write the image to a file called out.jpg using the imwrite function. Use some independent image viewer like display, xv or a web browser to verify that this worked. I recomend learning about the ImageMagick suite of tools for converting, and displaying images (do "man convert", and a "man display" to find out more). Also the program import can be used to get a screen-shot in linux.

7. Functions

   Go back to your hw1.m file and add the code to set every 5th horizontal/vertical pixel to 0, extract a user-selected image region using getrect from figure(2), display the new image in figure(3), and write the image to an image file.

   You can define a new Matlab function by simply putting code in a file (as we have just done) and placing a function declaration at the top. Write the following as the first line of your hw1.m file:

   function [imregion] = hw1(infile,outfile)

   Tip: The name of the file and the name of the function should always match. A file can export only one function, but the file may contain internal functions.

   Modify your code so that the imread and imwrite commands use the infile and outfile variables instead of hard-coded filenames. Also make sure that the extracted image region is stored in a variable called imregion so that it gets returned by the function. If the hw1.m file is in your pwd or in your Matlab path, then you can execute this function by typing hw1('mug.jpg','mugcrop.jpg'), or cropped=hw1('mug.jpg','mugcrop.jpg') if you want the cropped image returned in a variable.

   Tip: If you want to return more than one variable from a function, add it to the list in square brackets in the function's declaration. Simply setting these variables inside the function will cause them to be returned.

   Tip: If you want a Matlab function to modify a variable, then you have to duplicate it in the input and output lists. Matlab behaves like a functional language in this regard.

   Tip: Functions can also have optional input and output arguments. To determine the actual number of input and output arguments a function was called with, look at the special variables nargin and nargout.

8. Plotting

   Explore the plot command. Plot the sin function over the domain -pi:pi. Use the linspace command to define the domain x and then do plot(x,sin(x)). Use the hold on command to plot another function on the same graph, such as cos. Use a different color, e.g. plot(x,cos(x),'r'). The runing of hw1.m should produce a plot along these lines.

   Use the hist command to plot a histogram of the pixel values in a grayscale image. Use the h=hist(y,x) form so that you can control the bin centers and so that you can plot the result using plot(x,h). You may want to plot the log of the bin counts.

   Your program should now (three times consecutively) prompt the user to select rectangular regions from the gray scale version of the first image (figure 2) for which you will plot a histogram. It will be instructive to compare the histograms of different parts of the image, and thus previous plots should stay visible while additional ones are added.

   Your program should now prompt the user to select a point from the first image (figure 1) for which you will plot a horizontal trace of the image through that point. Compare with the "apricot" image in the lecture notes (the "apricot" image is available here ).

   It is instructive to think about histograms and traces of different parts of the image and different images. I recommend running your program on several images. However, for grading purposes, only the functionality described need to be implemented. In particular, there is no need to formally report the results of your experimentation.

9. Playing with Linear Algebra

   Matlab is a great tool to for experimenting with linear algebra.

   Use the fact that inv() inverts a matrix to solve for X=(x,y,z):

      3*x + 4*y +   z = 9
      2*x -   y + 2*z = 8
        x +   y -   z = 0 
   

   Verify that your "solution" works. Make sure that you give the answer in what you hand in.

   If there are more equations than unkowns, then, in the general case, "classically" the problem is over constrained and there is no solution. However, in this course, we will often be assuming that sych equations are approximations and have errors due to noise or other reasons, and that an exact solution cannot be found regardless. Thus we will want to find the "best" solution. This is konwn as solving the equations in the least squares sense. The solution for AX=b, where A has more rows than columns, is given by X=inv(A'*A)A'b, where inv(A'A)A' is known as the Moore-Penrose inverse of A. Use this to solve for X=(x,y,z) in:

      3.0*x + 4.0*y +  1.0*z = 9
      3.1*x + 2.9*y +  0.9*z = 9
      2.0*x - 1.0*y +  2.0*z = 8
      2.1*x - 1.1*y +  2.0*z = 8
      1.0*x + 1.0*y -  1.0*z = 0 
      1.1*x + 1.0*y -  0.9*z = 0 
   

   What is the magnitude of the error vector, AX - b?

   Recall that an eigenvector of a matrix A is a vector v, so that Av=kv, for some scalar constant k. If A is real and symmetric, then A has real eigenvalues and eigenvectors. Note that for a random matrix, R, R*R' is symmetric. (Try it!). The the matlab function eig() gives you eigenvectors and eigenvalues. Use these hints to give the TA a 3 by 3 matrix A, and a vector v, so that Av=kv, where k is a constant. In particular, the TA needs to see the output of Av./v which should be a vector of three elements that are all the same. Note that you will have to use the form of eig() that gives you two outputs (the default is just to give you one---you need to looks at the documentation).

10. Documenting Functions

    When you create a new function, it should always be documented so that help returns something informative. The convention in Matlab is to place the help message in comments after the function declaration (the comment character is "%"). For an example, you can look at some of the code in the Matlab library. The jet command that we used above to set the colormap exists as some .m file on the system. Type which jet to locate it, and then use the type command to view the source code. You see that the first line of the comment contains a one-line description of the command. All subsequent contiguous comment lines are included in the help message. Document your hw1.m function in this style.

11. Using C/C++ instead. Doing the assigments in C/C++ is relatively straightforward with appropriate library support. There are many options for this that you can use, including the UA CS vision group library. (If you expect that you might join/collaborate with the vision group, you will want to get comfortable with this software.) If you want to experiment with this software, then have a look at this example program (information about the library is in comments). To run the program you will also need a matrix input file called matrix.txt and a file named image.tiff. .

    A good warm-up project is to implement the least squares computation above. In increasingly simplicity, but decreasing enlightenment, you could: 1) Use basic matrix operations, including get_matrix_inverse(); 2) Shorten things using get_MP_inverse(); 3) Shorten things even further with least_squares(). Doing all three might be a good way to play with the library.

What to Hand In