ENGI 398: Engineering Mathematics Project 3: Application of Singular Value Decomposition to Image Compression1 1 Background The singular value decomposition (SVD) is a fundamental concept in linear algebra. The SVD has many applications in science and engineering. Some uses of the SVD can be found in computer graphics, image processing, pattern recognition, noise reduction, information retrieval, robotics, gene expression analysis, computational tomography, and latent semantic indexing (LSI technique incorporated in the Google search engine). The purpose of this project is to illustrate the application of the SVD to the compression of digital images. The idea is that given a digital image, we use the SVD of the image to produce a new image that requires less computer memory for its storage. Of course, we want to achieve this reduction in data storage requirement without incurring appreciable visual degradation of the image. 2 Compression of Digital Images 2.1 The Singular Value Decomposition (SVD) Given an m n matrix A representing an image, the singular value decomposition (SVD) is A = UVT , (1) where U is an mm orthogonal matrix, is an mn matrix, and V is an nn orthogonal matrix. That is, U and V are such that UTU = Imm and VTV = Inn. (2) In Equation (1), the columns of U are the left singular vectors; has the same size as A and contains the singular values in its diagonal; and VT has rows that are the right singular vectors. For an image A, the singular values are given along the diagonal of and are arranged in descending order. The columns of U contain an orthogonal basis for the column space of the image, while the columns of V contain an orthogonal basis for the row space of the image. The singular values are always nonnegative real numbers. If the matrix A is a real matrix, then U and V are also real; this is precisely the case for matrices representing images. Calculating the SVD consists of finding the eigenvalues and eigenvectors of AAT and ATA. The eigenvectors of AAT constitute the columns of U, while the eigenvectors of ATA constitute the columns of 1Based on material developed by Prof. Lige Li and Prof. Diego Maldonado of the Department of Mathematics at Kansas State University, and by Prof. Kristian Sandberg of the Department of Applied Mathematics of the University of Colorado at Boulder. 1 V. Also, the singular values in are the square roots of the eigenvalues from AAT or ATA. The SVD algorithm has been implemented in MATLAB. The syntax is as follows: [U,S,V] = svd(A) For example, take the following 2 3 matrix A = 1 ?2 0 4 ?1 1 . To find the SVD of the above matrix, we type the MATLAB instructions A = [1 -2 0; 4 -1 1]; [U,S,V] = svd(A) and MATLAB will return U = -0.3641 -0.9313 -0.9313 0.3641 S = 4.5106 0 0 0 1.6291 0 V = -0.9066 0.3224 -0.2722 0.3679 0.9198 -0.1361 -0.2065 0.2235 0.9526 which correspond to U = ?0.3641 ?0.9313 ?0.9313 0.3641 , = 4.5106 0 0 0 1.6291 0 , V = 24 ?0.9066 0.3224 ?0.2722 0.3679 0.9198 ?0.1361 ?0.2065 0.2235 0.9526 35 in Equation (1). You can verify that A = UVT . 2.2 The Energy of an Image In many cases it is useful to associate an energy to an image. If we want to compress an image, we usually want to remove as much data as possible without losing too much energy of the image. The difference in energy between the original image and the compressed image is one possible way to define an error measure when compressing images. However, it is important to notice that this measure is not always appropriate as an error measure. In some cases the energy error measure might be small, while the visual error is terrible. 2 One difficult problem in image processing is to find a good way of measuring error so that the error measure agrees with our visual perception. Let a digital image be represented as the m n matrix A and denote its elements as aij , i = 1, 2, . . . ,m and j = 1, 2, . . . , n. We define the Frobenius norm of the matrix A as kAkF =vuut m Xi=1 n Xj=1 |aij |2 . (3) This is not the most common matrix norm, and not always very convenient, but it is related to the so-called L2-norm for functions. The L2-norm (and the l2-norm for sequences) is related to the concept of energy and for this reason the Frobenius matrix norm is sometimes useful in image processing. What is the relation between the energy of an image and the singular values of the matrix representing the image? The answer is given by the following theorem: Theorem. Let an m n matrix A have the singular values 1, 2, . . . , r, where r = rank(A). Then, kAk2 F = r Xi=1 2 i . (4) Proof: Let A have an SVD of the form A = UVT =24 | | | u1 u2 . . . um | | | 35 266666666664 1 0 0 0 0 . . . 0 0 2 0 0 0 . . . 0 . . . ... 0 0 0 r 0 . . . 0 0 0 0 0 0 . . . 0 . . . 0 0 0 0 0 . . . 0 377777777775 266664 ? vT 1 ? ? vT 2 ? ... ? vTn ? 377775 . We will make use of these two facts: First note that for any matrix P expressed as a set of columns, i.e., P =24 | | | p1 p2 . . . pk | | | 35 , we have that kPk2 F = kp1k2 + kp2k2 + . . . + kpkk2. Also, recall that when an orthogonal transformation R is applied to a matrix P, the columns of P will undergo a rotation without changing their lengths. That is, kRPk2 F =

R24 | | | p1 p2 . . . pk | | | 35

2 F 3 =

24 | | | Rp1 Rp2 . . . Rpk | | | 35

2 F = kRp1k2 + kRp2k2 + . . . + kRpkk2 = kp1k2 + kp2k2 + . . . + kpkk2 = kPk2 F . Now, let P , VT . Then, kAk2 F = kUVT k2 F = kUPk2 F = kPk2 F (U is orthogonal) = kVT k2 F = k(VT )T k2 F = kVT k2 F = kT k2 F (V is orthogonal) = kk2 F = 2 1 + 2 2 + . . . + 2 r . = r Xi=1 2 i . 2.3 Rank One Decomposition Let A be an m n matrix and assume that it has rank r min(m, n). Without loss of generality, write the SVD of A in the form A = UVT =24 | | | u1 u2 . . . um | | | 35 266666666664 1 0 0 0 0 . . . 0 0 2 0 0 0 . . . 0 . . . ... 0 0 0 r 0 . . . 0 0 0 0 0 0 . . . 0 . . . 0 0 0 0 0 . . . 0 377777777775 266664 ? vT 1 ? ? vT 2 ? ... ? vTn ? 377775 . This in turn can be expressed as A = 1u1vT 1 + 2u2vT 2 + . . . + rurvT r . (5) The preceding expansion is called the rank one decomposition of A. Observe that each term iuivT i in the decomposition is a rank one m n matrix and hence the name. 4 0 2 4 6 8 10 0 100 200 300 400 500 600 Magnitude decay of singular values Singular value Magnitude of singular value Figure 1: Magnitude of singular values of the matrix A given in Equation (6). The rapid decay in the magnitude of the singular values suggests that only a small number of them make a significant contribution to kAkF , the energy of the image. 2.4 Lower Rank Approximation Consider the following 10 10 matrix A = 2666666666666664 87 16 90 63 62 66 59 70 96 75 1 19 1 73 73 39 57 48 6 74 14 42 30 38 19 63 72 11 36 43 82 86 5 1 90 70 51 66 55 63 43 49 69 42 57 40 78 37 26 80 89 82 65 75 63 41 49 14 60 8 73 46 98 79 23 66 19 57 5 95 69 46 55 92 55 84 70 82 57 92 35 45 40 84 93 37 98 67 70 60 17 41 20 37 34 43 81 100 96 25 3777777777777775 . (6) This matrix has singular values {558.63, 136.47, 115.54, 97.31, 86.93, 64.30, 54.28, 41.30, 9.82, 1.56}. The plot shown in Figure 1 reveals a rapid decay of the singular values of A. This suggests the idea of a lower rank approximation to A by keeping only the largest singular values and discarding the smallest ones. In image processing parlance, the contribution to the energy of the image by the largest singular values will 5 be more significant than that of the smallest singular values. Neglecting the smallest singular values, and keeping the largest ones, is the principle upon which image compression is based. A rank k approximation to matrix A, call it Ac, is formed by keeping the first k terms of the rank one decomposition of A, as given in Equation (5). For data compression to be practical, we choose k < r (r = min(m, n) is the rank of the original mn matrix A). Thus, the rank k approximation can be written as Ac = k Xi=1 iuivT i . (7) The rank k approximation to A tells us that the columns of Ac are linear combination of just k column vectors. A similar statement can be made about the rows of Ac. Let Ac represent a compressed version of the image A. We define the relative (Frobenius) error, ", as " = kA ? AckF kAkF . (8) The approximation error, kA?AckF , can be quantified in terms of the singular values. Naturally, this error will be a function of the singular values that were ignored in the approximation. More precisely, kA ? Ack2 F =

r Xi=1 iuivT i ? k Xi=1 iuivT i

2 F =

r X i=k+1 iuivT i

2 F = r X i=k+1 2 i (by Equation (4)). (9) 2.5 Data Compression In Section 2.4 we introduced a measure of the error incurred by compression (cf. Equation (8)). We also want a way to measure how efficient our compression is in terms of information we need to store. Let nA denote the number of memory units (e.g., bytes) we need to represent an image A, and let nAc denote the number of memory units we need to represent a compressed version of the image. Here, Ac denotes the compressed version of the image A. For an image array A of size m n there are nA = mn (10) numbers that need to be stored. However, using a rank k approximation as given by Equation (7), we need only store the first k singular values (the k largest entries in ) the first k columns of U, each of which has m entries the first k rows of VT , each of which has n entries 6 for a total number nAc = k(m + n + 1). (11) We can now measure the compression efficiency as C = nA nAc . (12) Hence, C = 1 means no compression at all and C = 10 means that we only need 10% as much memory to store Ac compared to A. If C = 10 we say that we have a compression ratio of 10:1. 3 Handling Images with MATLAB As mentioned in Section 2.2, a digital image is an array of picture elements or pixels. In MATLAB an image is represented by an array of nonnegative real numbers. These numbers represent pixel intensities. A monochrome (grayscale) image is represented as a pixel array of dimension mn, where m is the number of rows and n is the number of columns. For an RGB (color) image, the array consists of 3 grayscale mn arrays (red, green, and blue) stacked into an array of dimension m n 3. MATLAB can read and write image files in a variety of formats, including JPEG, TIFF, GIF, BMP, and PNG. In this project you will retrieve images from, and store them as, JPEG files. In the JPEG file format, pixel values are integers ranging from 0 to 255, i.e., 8-bit unsigned integers. The appropriate data type inMATLAB is uint8. In a grayscale image, 0 is black and 255 is white, thus rendering 256 shades of gray. To display a grayscale image in MATLAB, type A = imread(filename.jpg); % reads image in filename.jpg % and stores it in array A image(A); % displays array A as image colormap(gray(256)); % uses gray scale axis(image); % fits data to window frame axis off; % removes tick marks For a color image you can simply try B = imread(filename.jpg); image(B); axis(image); axis off; To extract the red, green, and blue color layers of image B, type Bred = B( :, :, 1); Bgreen = B( :, :, 2); Bblue = B( :, :, 3); You can create a color image by stacking monochrome images as shown below: C = zeros(m,n,3); C( :, :, 1) = Cred; 7 Original image Rank 8 approximation Rank 32 approximation Rank 128 approximation Figure 2: An image and different rank approximations. C( :, :, 2) = Cgreen; C( :, :, 3) = Cblue; The command A = imread(filename.jpg) reads a JPEG file and converts the result to an image array of MATLAB data type uint8. If you need to do arithmetic on this array you must first convert it using the function double. MATLAB provides a function with the same name (uint8) to convert a real (double) array back to data type uint8. Suppose you write a = uint8(x). If x is a real number greater than or equal to 255, then a will be 255. If x is a real number less than or equal to zero, then a will be 0. If 0 < x < 255 then a is the uint8 representation of x rounded down or floored to the nearest integer. Figure 2 shows what you can do with the SVD. The original grayscale image was a 512512 array, thereby requiring 262 144 bytes to store the image array. With the original image, rank 8, 32, and 128 approximations are shown. The approximations require 8 200, 32 800, and 131 200 bytes, respectively, with corresponding 32:1, 8:1, and 2:1 compression ratios. Notice that a rank approximation as low as 8 reveals some structure of the image, which has rank 512! It is amazing how very cool stuff can be done with nice mathematics. Enjoy your project! 4 MATLAB Project In this project you will be asked to analyze digital images using the SVD. Before you start the exercises, familiarize yourself with theMATLAB functions: imread, image, imwrite, svds, uint8, and double. For example, type help imread at the MATLAB prompt to learn more about the imread command. 8 4.1 Simple Image Manipulations 1. Download the files goldhill.jpg and krusty.jpg from the course web page. Use the command imread to read them into MATLAB and store them in matrices. For example, A = imread(goldhill.jpg); B = imread(krusty.jpg); Use the command size to find the dimensions of A and B. Without displaying the images, how can you tell whether the images are monochrome (black and white) or RGB (color)? Take a look at some of the entries for A and B, say A(10,288) or B(150,22,1). What do these numbers mean? Now use image to display A and B. 2. Create a subimage of krusty.jpg. You will do this by creating a submatrix of B, the matrix in which krusty.jpg is stored. Define the matrix subB as the submatrix of B given by subB = B(1:200,100:300,:); This will store in subB the submatrix of B consisting of the rows 1 through 200, columns 100 through 300, in each layer of B. Save this matrix subB as a JPEG image called kface.jpg as indicated below: imwrite(subB,kface.jpg); This will save the image stored in subB as kface.jpg in the current directory. Display the image kface.jpg by typing image(subB); to show the image stored in subB. Then, try C = imread(kface.jpg); image(C); 3. Download the file goldengate.jpg from the course web page. Read the file into MATLAB and display it using image: D = imread(goldengate.jpg); image(D); You should see an RGB image of the Golden Gate Bridge of San Francisco. Find the dimension of matrix D by typing size(D) and check that MATLAB returns ans = 108 250 3 The last 3 indicates that this is an RGB or color image, and each color plane (red, green, and blue) is a 108 by 250 array. Extract the red, green, and blue color layers by doing the following: 9 Dred = D(:,:,1); Dgreen = D(:,:,2); Dblue = D(:,:,3); Use the command subplot to show the original color image and its three color planes in the same figure. Remember that each color plane is a grayscale image. Now do the following: DLayer = D; % define DLayer as the whole image DLayer(:,:,2) = zeros(108,250); % zero out green plane DLayer(:,:,3) = zeros(108,250); % zero out blue plane image(DLayer) What is the visual result? Modify the preceding instructions to get the green and blue layers of the image. 4. Download the file toucan.jpg from the course web page. Store the image in matrix E and display it using image. Perform simple arithmetic operations on the image. If necessary, convert to double all uint8 data before you do arithmetic operations. You can then convert back to 8-bit unsigned integers using the uint8 function. Just for fun, display the result of the following operations: image(E); image(2*E); image(E.2); Comment on the visual effects that each of the above operations had on the original image. Generate a random matrix R of the same size as E, with values uniformly distributed between 0 and 255: R = uint8(255*rand(size(E))); Display the image of the random matrix R: image(R); How about the matrix R + E? Type image(R + E); Comment on your observations of the images of R and R + E. 4.2 Compression of Grayscale Image Using the SVD 1. Download the file barbara.jpg from the course web page. Store the image in matrix A and then display it using image. 2. Compute the 100 largest singular values of A. Type Ad = double(A); [U,S,V] = svds(Ad,100); This will give an approximation of the matrix A by keeping its 100 largest singular values. 10 3. Show the image of the above approximation by doing the following: A100 = U*S*V; A100 = uint8(A100); figure; image(A100); 4. Find the size of A and use Equation (10) to calculate the amount of memory needed to store A. Keep in mind that each entry of A occupies 8 bits or 1 byte. Then use Equation (11) to calculate the amount of memory required to store the rank 100 approximation. Determine the compression ratio for the rank 100 approximation. According to your visual perception, is A100 an acceptable approximation of A?What happens if you keep the 200 largest singular values? How much memory does A200 take up? Is A200 an acceptable approximation? What happens if you keep the 70 largest singular values? How much memory does A70 take up? Is A70 an acceptable approximation? 4.3 Compression of Color Image Using the SVD 1. In Section 4.2 you compressed a black and white image using different rank approximations. In this exercise you will compress an RGB or color image. Download the file track.jpg from the course web page. Store the image in matrix B and then display it using image. 2. Approximate B with rank 8, 32, and 128 matrices. For RGB images do it layer by layer; simply extend the ideas presented in Section 4.2 for monochrome images. Display the resulting approximations and comment on the quality of the images. Use subplot to show the original and approximate images in the same figure. 4.4 Analysis of Compressed Images 1. The SVD of an Image Download the file bridge.jpg from the course web page. Read it into MATLAB and display it. From the SVD of the image, display a rank 1 approximation of the image. What can you say about the columns in the rank 1 approximation of the image in terms of linear dependence/independence? Does your image appear to agree that you indeed produced a rank 1 approximation? Next, reconstruct your image using the SVD of rank 2, 4, 8, 16, 32, 64, and 128. Using the command subplot, plot all these approximations along with the original image in the same window for an easy comparison. 2. The Decay of the Singular Values Make a plot that shows how the singular values decay along the diagonal of the matrix in the singular value decomposition A = UVT . How do the singular values decay for the image bridge.jpg? 11 3. Compression Error Using the SVD, approximate the image bridge.jpg with a rank 64 matrix. Compute the relative error using the Frobenius norm as defined in Equation (8). To find the Frobenius norm of an array M, type the MATLAB instruction norm(M,fro). Compare the relative errors between the images and discuss the results. 4. Compression Based on Visual Perception How much can you compress bridge.jpg using the SVD such that your compressed image is almost perfect based on your visual perception? How large did the rank need to be before it looked natural to your eyes? 5. Illustration of the Frobenius Norm Reconstruct bridge.jpg using two different methods: A rank 100 approximation. Using all but the first two singular values of the image. Compute the relative error in the Frobenius norm for both of the approximation. Does the result of the norm roughly agree with the error based on visual perception? 6. Compression Based on the Frobenius Norm For the image bridge.jpg, how large does the rank need to be so that the relative error (in the Frobenius norm) is less than 10%? 5 Project Report Your project report should include: Brief explanation of the tasks to be performed in each of the proposed exercises. Listings of MATLAB programs and figures. Answers to the questions in the MATLAB Project section and conclusions. Remember the following: Prepare your report using a word processor (MicrosoftWord or similar) and turn your report in on the due date. 12