[SOLVED] Math551 - lab11

Description

Goal: In this lab, you will learn to use the singular value decomposition for image compression.

Getting started: You will need to download the file TarantulaNebula.png and put it in your working directory.

 

Tasks

1. Open an M-file called m and add the following commands

t=linspace(0,2*pi,100); X=[cos(t);sin(t)]; figure(1);

plot(X(1,:),X(2,:),’b’); axis equal

this will create a 2 × 100 matrix X = [x₁ x₂ ··· x₁₀₀] whose columns are unit vectors pointing in various directions. The plot will show a blue circle corresponding to the endpoints of these vectors. Variables: X

2. Now in the M-file, define a variable A holding the matrix

A

and a variable AX holding the product A*X. Note that this product equals

AX = [Ax₁ Ax₂ ··· Ax₁₀₀],

so the columns of AX show where the matrix A maps each of the columns of X. Add the following lines to your M-file after the definition of AX.

figure(2); plot(AX(1,:),AX(2,:),’b’); axis equal

Now when you run the M-file you should see a blue circle in Figure 1 and a blue ellipse in Figure 2. This shows that A maps vectors that end on the unit circle to vectors ending on an ellipse. Variables: A, AX

3. Add the following line to your M-file and re-run the file.

[U,S,V] = svd(A)

Your output should look like this

U =

-0.9571	0.2898
0.2898 S =	0.9571
2.3028	0
0 V =	1.3028
-0.9571	-0.2898
-0.2898	0.9571

The function svd computes the singular value decomposition of A:

A = USV^T

where U = [u₁ u₂] and V = [v₁ v₂] are orthogonal matrices and

S

is a non-negative diagonal matrix. Moreover

[Av₁ Av₂] = AV = USV^TV = US  ,

showing that

Av₁ = σ₁u₁                  and           Av₂ = σ₂v₂.

You can check this fact in the command window as follows.

>> A – U*S*V’ ans =

1.0e-15 *

0.8882          0.4441

-0.4441          0.2220

>> U’*U ans =

1.0000                    0

• 0000

>> V’*V ans =

• 0

0           1

Variables: U, S, V

4. Plot the columns of V on Figure 1 by adding the following to your M-file.

figure(1); hold on; quiver(0,0,V(1,1),V(2,1),0,’r’); quiver(0,0,V(1,2),V(2,2),0,’g’); hold off;

The vector corresponding to the first column will be plotted in red and the second column in green.

5. Plot the columns of US on Figure 2 by adding the following to your M-file.

figure(2); hold on; quiver(0,0,S(1,1)*U(1,1),S(1,1)*U(2,1),0,’r’); quiver(0,0,S(2,2)*U(1,2),S(2,2)*U(2,2),0,’g’); hold off;

The vector corresponding to the first column will be plotted in red and the second column in green.

6. Now we’ll work on a bigger matrix. Add the following commands to your M-file and run it. (If you want, you can create a new cell by entering two percent signs (%%) at the beginning of a line and just run the cell from now on.)

%% (start a new cell)

Img = double(imread(’TarantulaNebula.png’)); figure(3); image(Img); colormap(gray(256));

This code loads an image as a real matrix and displays it on the screen. Each entry in the matrix corresponds to a pixel on the screen and takes a value somewhere between 0 (black) and 255 (white).

Variables: Img

7. Perform the singular value decomposition of Img, save the output in matrices UImg, SImg, and VImg.

Variables: UImg, SImg, VImg

8. Add the following code and run the M-file.

figure(4); semilogy(diag(SImg));

This shows the singular values (the diagonal entries of the S matrix) for our image matrix (the scale of the y axis is logarithmic). Notice that the diagonal entries of S are ordered so that σ₁ ≥ σ₂ ≥ σ₃ ≥ ···. One way of writing the SVD is

A

where r is the rank of A and u_i and v_i are the ith columns of U and V respectively.

If for some k < r, σ_k+1 is very small compared to σ₁, we should expect

A

to be a good approximation of A. (This is called a truncated SVD.) This idea can be used for image compression as follows. Instead of storing the whole m × n matrix A, we instead store the m × k and n × k matrices

C = [u₁ u₂ ··· u_k]             and              D = [σ₁v₁ σ₂v₂ ··· σ_kv_k].

If k(m + n) is much smaller than mn, then storing C and D will take much less space than storing A. Moreover, if we wish to display the image, we can reconstruct A (approximately) as A ≈ CD^T.

9. Add the following code to your M-file. %% (start a new cell) k = 10;

% compressed image

• = UImg(:,1:k);
• = VImg(:,1:k)*SImg(1:k,1:k);

% compression percentage pct = 1 – (numel(C)+numel(D))/numel(Img);

figure(5); image(C*D’); colormap(gray(256)); title( sprintf( ’%.2f%% compression’, 100*pct ) );

figure(6); image(Img); colormap(gray(256)); title( ’Original’ );

This code compresses the image as described above using k = 10 singular values, then displays the reconstructed image along with the original. It also shows the percent decrease in storage size required for the compressed image. Try several values of k to see how the quality and compression percentage vary with k. Set k so that the compression percentage is 80.23% before turning in your M-file.

Variables: k, C, D