ABSTRACT Title of dissertation: IMAGE REPRESENTATION AND COMPRESSION VIA SPARSE SOLUTIONS OF SYSTEMS OF LINEAR EQUATIONS Alfredo Nava Tudela, Doctor of Philosophy, 2012 Dissertation directed by: Professor John J. Benedetto Department of Mathematics We are interested in nding sparse solutions to systems of linear equations Ax = b, where A is underdetermined and fully-ranked. In this thesis we examine an implementation of the orthogonal matching pursuit (OMP) algorithm [11], an algorithm to nd sparse solutions to equations like the one described above, and present a logic for its validation and corresponding validation protocol results. The implementation presented in this work improves on the performance reported in previously published work [2] that used software from SparseLab [19]. We also use and test OMP in the study of the compression properties of A in the context of image processing. We follow the common technique of image blocking used in the JPEG and JPEG 2000 standards [3, 22, 25]. We make a small modi cation in the stopping criteria of OMP that results in better compression ratio vs image quality as measured by the structural similarity (SSIM) and mean structural similarity (MSSIM) indices which capture perceptual image quality [26]. This results in slightly better compression than when using the more common peak signal to noise ratio (PSNR) [30]. We study various matrices whose column vectors come from the concatenation of waveforms based on the discrete cosine transform (DCT), and the Haar wavelet. We try multiple linearization algorithms and characterize their performance with respect to compression. An introduction and brief historical review on the topics of information theory, quantization and coding, and the theory of rate-distortion leads us to compute the distortion D properties of the image compression and representation approach presented in this work. A choice for a lossless encoder is left open for future work in order to obtain the complete characterization of the rate-distortion properties of the quantization/coding scheme proposed here. However, the analysis of natural image statistics is identi ed as a good design guideline for the eventual choice of . The lossless encoder is to be understood under the terms of a quantizer ( ; ; ) as introduced in [6]. IMAGE REPRESENTATION AND COMPRESSION VIA SPARSE SOLUTIONS OF SYSTEMS OF LINEAR EQUATIONS by Alfredo Nava Tudela Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park in partial ful llment of the requirements for the degree of Doctor of Philosophy 2012 Advisory Committee: Professor John J. Benedetto, Chair/Advisor Professor Radu V. Balan Professor Wojciech Czaja Professor Ramani Duraiswami Professor Stuart Vogel c Copyright by Alfredo Nava Tudela 2012 Preface During her excellent Scienti c Computing course in Spring 2010, Professor Dianne P. O?Leary o ered a variety of articles on which to base our nal course projects. Among them was an article from SIAM Review by Alfred M. Bruckstein, David L. Donoho, and Michael Elad [2] that caught my attention and on which I based my work. However, it was not until the Advanced Scienti c Computing year- long, two-course sequence that followed that I fully developed what would become the core of this dissertation, further inspired by this article. During this period, Professors Radu Balan and Manuel Tiglio o ered valuable input that shaped some of the results presented in this work. Enthused by my discovery of the eld of Sparseland, I talked to Professor John J. Benedetto, my long-time friend and advisor, about switching gears on what had been until then my dissertation work|trying to improve Magnetic Resonance Imaging|and instead exploring further image representation and compression by means of sparse solutions of systems of linear equations. Dr. Benedetto agreed to let me explore my newfound interest. As they commonly say, the rest is history. I am deeply indebted to all these wonderful professors. Especially to you John, who patiently over the years believed in me and, through our regular weekly meetings, saw me come to the end of this long and winding road. I would also like to thank Professor Konstantina Trivisa, director of the Ap- plied Mathematics & Statistics, and Scienti c Computation Program at the Univer- sity of Maryland, College Park, for persuading me to graduate while she was still director of the program. Her energy was a huge motivation to nish. ii Dedication To Kelly and Luc a. Thank you for your love and support. iii Table of Contents List of Tables vi List of Figures vii 1 Introduction 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Finding sparse solutions . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Image representation and compression . . . . . . . . . . . . . . . . . 3 2 Finding sparse solutions 6 2.1 Orthogonal Matching Pursuit . . . . . . . . . . . . . . . . . . . . . . 6 2.2 An OMP implementation . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2.1 Least-squares approximation by QR decomposition . . . . . . 8 2.2.2 Implementation ne tuning and speedup . . . . . . . . . . . . 10 2.3 OMP validation protocol and validation results . . . . . . . . . . . . 12 2.3.1 Theoretical results that motivate and justify the protocol . . . 12 2.3.2 Validation protocol . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3.3 Validation results . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.4 OMP testing protocol and results . . . . . . . . . . . . . . . . . . . . 18 2.4.1 Reproducing previous work . . . . . . . . . . . . . . . . . . . 21 3 Image representation and compression 23 3.1 Elementary image representation concepts . . . . . . . . . . . . . . . 23 3.2 Image compression via sparsity . . . . . . . . . . . . . . . . . . . . . 25 3.2.1 Choosing A: \One-dimensional" basis elements . . . . . . . . 26 3.2.2 Choosing A: \Two-dimensional" basis elements . . . . . . . . 29 3.2.3 Some properties of [DCT1 Haar1] and [DCT2;j Haar2;j] . . . . 31 3.3 Image database . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.4 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.5 E ects on image reconstruction from choice of map ci . . . . . . . . . 41 3.5.1 Results for A = [DCT1 Haar1] . . . . . . . . . . . . . . . . . . 41 3.5.2 Results for A = [DCT2;j Haar2;j] . . . . . . . . . . . . . . . . 47 3.6 Normalized bit-rate vs tolerance . . . . . . . . . . . . . . . . . . . . . 49 3.7 Image reconstruction and error estimation . . . . . . . . . . . . . . . 54 3.8 PSNR and MSSIM comparison . . . . . . . . . . . . . . . . . . . . . 59 3.9 Other matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.9.1 \Rotations" of the basis elements of DCT2;3 . . . . . . . . . . 69 4 Quantization and coding 74 4.1 Image quantization and coding . . . . . . . . . . . . . . . . . . . . . 81 5 New directions 97 iv A Basic frame de nitions 99 Bibliography 100 v List of Tables 2.1 Algorithm performance. The speedup column refers to the speedup from the immediately previous implementation. To compute the total speedup from rst to last implementations multiply all speedup values together. Total speedup from ompQR to ompQRf3 is 2.068, which means we doubled the speed of our implementation for the generic matrix version of our code. We used Matlab version R2010b Service Pack 1 to run \experiment.m" which performs many OMP calls on randomized input. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.1 List of images used. Image Peppers was converted to a grayscale image from the 24-bit color original using Photoshop CS3. Image Barbara is not in the USC-SIPI database but we included since it is widely used by the image processing community. . . . . . . . . . . . . 37 3.2 Performance results for A = [DCT1 Haar1] for each c1, c2, and c3. For each image in our test database, we linearized each 8 by 8 sub-image using either c1, c2, or c3. In all cases, the PSNR value was larger using c2; and in all cases, except for image Barbara|although minimally|, the normalized bit-rate was smaller. Both of these measures make c2 a better choice over c1 or c3. The values correspond to runs of OMP with an ?2 termination criteria, and a tolerance = 32. . . . . . . . . 43 3.3 Total variation averages Vavg(I; ci) and maximums Vmax(I; ci) . . . . 43 3.4 Performance results for A = [DCT2;1 Haar2;1], A = [DCT2;2 Haar2;2], and A = [DCT2;3 Haar2;3] with corresponding vectorization functions c1, c2, and c3, respectively. In all cases, the PSNR and normalized bit- rate values were almost identical, with minor di erences. Matrix A = [DCT2;3 Haar2;3] performs slightly better on average. Mismatching function ci with matrix A = [DCT2;j Haar2;j] when i 6= j, results in degraded performance. The values correspond to runs of OMP with an ?2 termination criteria, and a tolerance = 32. . . . . . . . . . . . 48 3.5 Performance results for A = [DCT2;3(0; 6) : : :DCT2;3(5; 6)]. In all cases both the PSNR and normalized bit-rate values were better than the values obtained for A = [DCT2;3 Haar2;3], i.e., larger and smaller, respectively. See Table 3.4 for comparison. The values correspond to runs of OMP with an ?2 termination criteria, and a tolerance = 32. 73 vi List of Figures 2.1 Suppose that after k = 3 iterations of the main loop, OMP has chosen, in the following order, columns a5, a2, and a7 from matrix A. We form sub-matrix A(3) = (a5 a2 a7), and its QR decomposition A(3) = Q(3)R(3), which we use to solve the least-squares problem kA(3)~x bk22 = 0, with ~x 2 R 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2 OMP behavior for a matrix A with (A) = 0:3713, which corresponds to k0 = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3 OMP behavior for a matrix A with (A) = 0:3064, which corresponds to k0 = 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.4 The three modal behaviors, dependent on 0, observed for the matrix A used in Figure 2.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.5 The three modal behaviors, dependent on 0, observed for the matrix A used in Figure 2.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.6 Reproduction of results in Section 3.3.1 of [2] for OMP. Our imple- mentation ompQRf3 is slightly better at recovering with higher prob- ability the sparsest solution to Ax = b when compared to SolveOMP, an implementation publicly available at [19]. . . . . . . . . . . . . . . 22 3.1 Top row: The rst six waveforms of the DCT1 (a), and the Haar1 (b) normalized bases for R64. Bottom row: Full 2D representation using the inverse map of c2 de ned in Section 3.4 for the DCT1 (c), and Haar1 (d) bases for R64. White corresponds to the maximum value achieved by the basis element, black to the minimum. The intermediate shade of gray corresponds to 0. . . . . . . . . . . . . . . 28 3.2 Full natural 2D representation for the DCT2 (a), and Haar2 (b) bases for R64. White corresponds to the maximum value achieved by the basis element, black to the minimum. The intermediate shade of gray corresponds to 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.3 Images used for our compression algorithms based on sparse image representation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.4 Two possible ways to vectorize a matrix. (a) Concatenate from top to bottom one after the other, from left to right, the columns of the matrix; or (b) rst ip every other column, and then concatenate the columns as before. This is what functions c1 and c2 do, respectively. . 39 3.5 How to vectorize a matrix following a zigzag pattern. A 4 by 4 ma- trix is shown with its 16 elements enumerated from left to right, top to bottom. The path of arrows shows the new induced order after traversing the matrix in a zigzag pattern to obtain a vector. . . . . . 40 3.6 Total variation for image Barbara for the vectorization function (a) c1, (b) c2, and (c) c3. Blue represents low values, red high values. . . 44 3.7 Total variation for image Boat for the vectorization function (a) c1, (b) c2, and (c) c3. Blue represents low values, red high values. . . . . 44 vii 3.8 Total variation for image Elaine for the vectorization function (a) c1, (b) c2, and (c) c3. Blue represents low values, red high values. . . . . 45 3.9 Total variation for image Peppers for the vectorization function (a) c1, (b) c2, and (c) c3. Blue represents low values, red high values. . . 45 3.10 Total variation for image Stream for the vectorization function (a) c1, (b) c2, and (c) c3. Blue represents low values, red high values. . . . . 45 3.11 Normalized bit-rate vs tolerance: \One-dimensional" basis elements. We can observe that for all images the best normalized bit-rate for a given tolerance is obtained for matrix A = [DCT1 Haar1] which combines both the DCT1 and Haar1 bases for R64. . . . . . . . . . . . 52 3.12 Normalized bit-rate vs tolerance. We show results to compare the performance of A = [DCT1 Haar1] and ~A = [DCT2;3 Haar2;3]. . . . . 53 3.13 Peak Signal-to-Noise Ratio vs tolerance. We observe three typical behaviors for all images. For large values of the tolerance, about > 40, the PSNR of all images is above the PSNR value for the idealized error distribution marked by the black dashed line. This behavior is also observed for very small values of the tolerance, about < 3. Then for values between these two extreme behaviors, all images conform very closely to the idealized error distribution, a fact that is expected from the least-squares approximation at the core of the OMP algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.14 Normalized bit-rate vs MSSIM, PSNR . . . . . . . . . . . . . . . . . 59 3.15 Peak Signal-to-Noise Ratio vs Mean Structural Similarity . . . . . . . 60 3.16 Normalized bit-rate and corresponding MSSIM vs tolerance. In this graph we have plotted together the best normalized bit-rate obtained by combining the DCT and Haar bases, and the corresponding value of the MSSIM index for a given tolerance. The normalized bit-rate graphs are on the bottom left, and the MSSIM index values are above these. This gure combines results for all images. . . . . . . . . . . . 61 3.17 Normalized bit-rate vs MSSIM. Comparison of di erent termination criteria for OMP. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.18 Barbara: 0 = 32, PSNR = 36.9952 dB, MSSIM = 0.9444, normalized bit-rate = 0.1863 bpp, termination criteria: k k2. . . . . . . . . . . . 64 3.19 Barbara: 0 = 0:94, PSNR = 32.1482 dB, MSSIM = 0.9462, normal- ized bit-rate = 0.1539 bpp, termination criteria: k kMSSIM . . . . . . 64 3.20 Boat: 0 = 32, PSNR = 36.6020 dB, MSSIM = 0.9210, normalized bit-rate = 0.1608 bpp, termination criteria: k k2. . . . . . . . . . . . 65 3.21 Boat: 0 = 0:92, PSNR = 34.1405 dB, MSSIM = 0.9351, normalized bit-rate = 0.1595 bpp, termination criteria: k kMSSIM . . . . . . . . . 65 3.22 Elaine: 0 = 32, PSNR = 36.5155 dB, MSSIM = 0.9096, normalized bit-rate = 0.1682 bpp, termination criteria: k k2. . . . . . . . . . . . 66 3.23 Elaine: 0 = 0:90, PSNR = 35.6288 dB, MSSIM = 0.9168, normalized bit-rate = 0.1686 bpp, termination criteria: k kMSSIM . . . . . . . . . 66 3.24 Peppers: 0 = 32, PSNR = 36.8674 dB, MSSIM = 0.8983, normalized bit-rate = 0.1024 bpp, termination criteria: k k2. . . . . . . . . . . . 67 viii 3.25 Peppers: 0 = 0:89, PSNR = 35.4309 dB, MSSIM = 0.9080, normal- ized bit-rate = 0.1011 bpp, termination criteria: k kMSSIM . . . . . . 67 3.26 Stream: 0 = 32, PSNR = 36.4686 dB, MSSIM = 0.9622, normalized bit-rate = 0.3050 bpp, termination criteria: k k2. . . . . . . . . . . . 68 3.27 Stream: 0 = 0:95, PSNR = 34.1398 dB, MSSIM = 0.9634, normal- ized bit-rate = 0.2853 bpp, termination criteria: k kMSSIM . . . . . . 68 3.28 Full natural 2D representation for (a) DCT2;3(0; 6), (b) DCT2;3(1; 6), (c) DCT2;3(2; 6), (d) DCT2;3(3; 6), (e) DCT2;3(4; 6), and (f) DCT2;3(5; 6), bases for R64. White corresponds to the maximum value achieved by the basis element, black to the minimum. The intermediate shade of gray corresponds to 0. . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.1 Schematic diagram of a communication system, from [16]. . . . . . . 78 4.2 Schematic diagram of a general transform coding system. . . . . . . . 82 4.3 (a) MSSIM vs Normalized bit-rate, (b) PSNR vs Normalized bit-rate. 89 4.4 PSNR vs bit-rate: (a) Normalized bit-rate results for A = [DCT1 Haar1] prior to any coding, and (b) bit-rate coding performances published in [14] for image Lena: Said and Pearlman?s SPIHT algorithm [15], Embeded Coding and the Wavelet-Di erence-Reduction compression algorithm (\new algorithm") [23], and Shapiro?s EZW algorithm [18]. 90 4.5 Histograms for image Barbara. A = [DCT2;3 Haar2;3], tolerance = 32. 91 4.6 Histograms for image Boat. A = [DCT2;3 Haar2;3], tolerance = 32 . 92 4.7 Histograms for image Elaine. A = [DCT2;3 Haar2;3], tolerance = 32 93 4.8 Histograms for image Peppers. A = [DCT2;3 Haar2;3], tolerance = 32 94 4.9 Histograms for image Stream. A = [DCT2;3 Haar2;3], tolerance = 32 95 4.10 Histograms for uniform random input. A = [DCT2;3 Haar2;3], toler- ance = 32 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 ix Chapter 1 Introduction 1.1 Motivation In recent years, interest has grown in the study of sparse solutions of under- determined systems of linear equations because of their many and potential appli- cations [2]. In particular, these types of solutions can be used to describe images in a compact form, provided one is willing to accept an imperfect representation. We are interested in studying this approach to image compression because of the ever-increasing volume of images in use by many multimedia channels. The basic idea is the following. Suppose that we have a full-rank matrix A 2 Rn m, where n < m, and that we want to nd solutions to the equation, Ax = b; (1.1) where b is a given \signal." Since the matrix A is full-rank, and we have more unknowns than equations, we have an in nite number of solutions to Equation 1.1. What if from all possible solutions we could nd x0, the \sparsest" one, in the sense of having the least amount of nonzero entries? Then, if the number of nonzero entries in x0 happens to be less than the number of nonzero entries in b, we could store x0 instead of b, achieving a representation of the original signal b in a compressed 1 way. There are many questions that arise to this approach for image representation and compression. For example, is there a unique \sparsest" solution of Equation 1.1? How do you nd such solutions? What are the practical implications for this approach to image compression? We note that this idea can be framed in the context of signal transform com- pression techniques. For example, the JPEG and JPEG 2000 standards have at their core transformations that result in di erent representations of the original im- age which can be truncated to achieve compression at the expense of some acceptable error [3, 22, 25]. 1.2 Finding sparse solutions The orthogonal matching pursuit (OMP) algorithm [2, 11] is one of the existing techniques to nd sparse solutions of systems of linear equations, as in Equation 1.1, where A 2 Rn m is a full-rank matrix, and n < m. This is one of many greedy algorithms that attempt to solve the general problem, (P 0) : minx kxk0 subject to kAx bk2 < : (1.2) Here, kxk0 = #fj : jxjj > 0g is the \zero-norm" of vector x, which counts the number of nonzero entries in x. A greedy algorithm approach is necessary in the solution of the optimization problem de ned in Equation 1.2 because this is an NP- complete problem [12]. Moreover, it can be proven that under certain circumstances 2 there is a unique sparsest solution to (P 0); and, under those same circumstances, OMP is then guaranteed to nd it [2]. We have implemented OMP using a QR matrix decomposition [20] in one of the steps of the algorithm in a way that takes advantage of calculations performed in previous steps, resulting in a signi cant speedup. For details, see Section 2.2. Our implementation improves on the convergence to the sparsest solution of Equation 1.2 as compared to results previously published and obtained using software from [19], see Section 2.4. Moreover, compared to our initial naive implementation of OMP as de ned in [2], we modi ed certain aspects of the algorithm that resulted in further speedup. Armed with a tested and validated implementation of OMP (Section 2.3), we proceeded to study image representation and compression as explained before. 1.3 Image representation and compression To make a practical implementation of the compression technique described in the introduction, we looked for inspiration in the JPEG and JPEG 2000 standards [25, 22]. Also, in order to test our image representation and compression approach, we selected a set of ve commonly used test images, four of them from [24]. All images in our image database are 512 by 512, 8-bit depth, grayscale im- ages. We proceeded to partition each image in subsets of 8 by 8 non-overlapping sub-images to process them individually. Partitioning a signal to work with more manageable pieces is a common technique [29]. 3 Subsequently, we vectorize each 8 by 8 sub-image into a vector b 2 R64 to be used as a right hand side in Equation 1.1. There are many ways to do this vectorization and we explored three of them in detail, see Sections 3.4 and 3.5. To complete the setup, we needed a matrix A. We decided initially on A = [DCT1 Haar1] 2 R64 128, where DCT1 represents a basis of one-dimensional discrete cosine transform waveforms, and Haar1 a basis of Haar wavelets, respectively. See Section 3.2.1 for details. That is, we concatenated two bases of R64, where b lives. Our initial thought for choice of basis elements drew from one-dimensional waveforms given that we are trying to approximate a vector|a one-dimensional object|, however we also considered bases for R64 built from tensor products of the one-dimensional waveforms mentioned above to capture the two-dimensional nature of the underlying problem, the representation of an image, an inherently two-dimensional object. We constructed in speci c matrices A = [DCT2;j Haar2;j] de ned in Section 3.2.2. We have compared the compression properties of A = [DCT1 Haar1] to those of B = [DCT1] and C = [Haar1] that only use the one-dimensional discrete co- sine transform waveforms, or the one-dimensional Haar wavelets, respectively, and found that combining bases results in a representation x0 for each sub-image that requires fewer nonzero entries. We found similar compression results for the two- dimensional bases. A comparison between the one-dimensional and two-dimensional concatenation of bases showed that up to a range of tolerance values between ap- proximately 3 and higher the two-dimensional basis elements perform better than the one-dimensional ones. See Section 3.6. 4 Also, as part of our research we have measured the quality of the reconstruc- tion from these compressed representations by way of the peak signal-to-noise ratio (PSNR) [22, 30], the structural similarity index (SSIM), and the mean structural sim- ilarity index (MSSIM) [26], all as functions of the tolerance chosen when solving (P 0), see Sections 3.7 and 3.8. This led us to a novel modi cation of the termination criteria in the OMP algorithm in order to achieve a desired PSNR or MSSIM value for the resulting decompressed image. There is still some work left to do on how to modify OMP to optimize for a desired MSSIM value. Even though this is a promising compression technique, the full potential of this approach to image representation and compression cannot be assessed until a bit-stream c is produced to be able to quantify the net compression ratio, among other criteria. That is, we still need to tackle the elements found in complete image compression standards that draw on the information-theoretical paradigms raised by the work of Shannon [16], that are incorporated into working JPEG and JPEG 2000 implementations, and in general information transmission/storage systems. Along these lines, we gave a brief review and historical introduction to the top- ics of information theory, quantization and coding, and the theory of rate-distortion, see Chapter 4. The structure of the image compression approach that we have cho- sen is amenable to a ne computation of its distortion D properties, for details see Section 4.1 Finally, we o ered a few new directions where future research can lead into in Chapter 5, and in that spirit we brie y explored variations on the matrices we used in the bulk of our work. For details see Section 3.9. 5 Chapter 2 Finding sparse solutions 2.1 Orthogonal Matching Pursuit In the rst half of this project, we are interested in implementing and validating one of the many greedy algorithms (GAs) that attempt to solve (P 0). The general idea is as follows. Starting from x0 = 0, a greedy strategy iteratively constructs a k-term approximation xk by maintaining a set of active columns|initially empty| and, at each stage, expanding that set by one additional column. The column chosen at each stage maximally reduces the residual ?2 error in approximating b from the current set of active columns. After constructing an approximation including the new column, the residual error ?2 is evaluated; if it now falls below a speci ed threshold, the algorithm terminates. Orthogonal Matching Pursuit (OMP)|a GA for calculating the solution to (P 0): Task: Calculate the solution to (P 0) : minx kxk0 subject to kAx bk2 < . Parameters: We are given the matrix A, the vector b, and the threshold 0. Initialization: Initialize k = 0, and set The initial solution x0 = 0. 6 The initial residual r0 = b Ax0 = b. The initial solution support S0 = Supportfx0g = ;. Main Iteration: Increment k by 1 and perform the following steps: Sweep: Compute the errors (j) = minzj kzjaj r k 1k22 for all j using the optimal choice z j = a T j r k 1=kajk22. Update Support: Find a minimizer j0 of (j): 8j =2 Sk 1; (j0) (j), and update Sk = Sk 1 [ fj0g. Update Provisional Solution: Compute xk, the minimizer of kAx bk22 subject to Supportfxg = S k. Update Residual: Compute rk = b Axk. Stopping Rule: If krkk2 < 0, stop. Otherwise, apply another iteration. Output: The proposed solution is xk obtained after k iterations. This algorithm is known in the literature of signal processing by the name orthogonal matching pursuit (OMP), and this is the algorithm we have implemented and validated. OMP solves, in essence, (P 0) for = 0, a given positive threshold. See Equation 1.2 in Section 1.2 for details. 2.2 An OMP implementation For a given matrix A 2 Rn m, if the approximation delivered by OMP has k0 zeros, the method requires O(k0mn) ops in general; this can be dramatically 7 better than the exhaustive search, which requires O(nmk0k20) ops. 2.2.1 Least-squares approximation by QR decomposition We would like to make the following observations about the OMP algorithm described in Section 2.1. The step that updates the provisional solution seeks to minimize kAx bk22, subject to Supportfxg = S k. This is equivalent to solving the least-squares approximation problem min~x kA(k)~x bk22 for the matrix A (k) that results from using only the k active columns of A de ned by Sk, and ~x is the vector in Rk whose i-th entry corresponds to the column of A that was chosen during the i-th iteration of the main loop. See Figure 2.1. n 1 1 1 a5 a2 a7 A = Q n 3 n-3 Q2 = x n-3 3 R 0 0 3 R2 R1 (3) (3) (3) Q1(3) (3) (3) (3) Figure 2.1: Suppose that after k = 3 iterations of the main loop, OMP has chosen, in the following order, columns a5, a2, and a7 from matrix A. We form sub-matrix A(3) = (a5 a2 a7), and its QR decomposition A(3) = Q(3)R(3), which we use to solve the least-squares problem kA(3)~x bk22 = 0, with ~x 2 R 3. For the case when A is a relatively small matrix, we can solve this problem, for example, by factorizing A(k) = Q(k)R(k) with the QR-algorithm, and then observing 8 that A(k) = Q(k)R(k) = Q(k)1 R (k) 1 + Q (k) 2 R (k) 2 = Q (k) 1 R (k) 1 + 0 = Q (k) 1 R (k) 1 ; (2.1) where|using Matlab notation| Q(k)1 = Q (k)(:; 1:k); Q(k)2 = Q (k)(:; k+1:n); R(k)1 = R (k)(1:k; :); and R(k)2 = R (k)(k+1:n; :): Then, from Equation 2.1, we have A(k)~x0 = b, Q (k) 1 R (k) 1 ~x0 = b ) Q(k)T1 Q (k) 1 R (k) 1 ~x0 = Q (k)T 1 b , R(k)1 ~x0 = Q (k)T 1 b , ~x0 = (R (k) 1 ) 1Q(k)T1 b; where ~x0 2 Rk is the solution to the equivalent minimization problem described above, and the inverse of R(k)1 exists because A (k) is full-rank. Finally, when OMP returns successfully after k0 iterations, we embed ~x0 2 Rk0 in 0 2 Rm \naturally" to obtain the solution x0 2 Rm to the initial least-squares approximation problem kAx bk22 subject to the nal active column set S k0 . The natural embedding refers 9 to setting the j-th entry of 0 2 Rm equal to the i-th entry in ~x0 2 Rk0 if during the i-th loop of the main algorithm, OMP chose the j-th column of A. 2.2.2 Implementation ne tuning and speedup We went through a series of code iterations to speedup our original imple- mentation ompQR, initially done from a simplistic reading of the OMP algorithm described in Section 2.1. We also had a generic implementation and a couple of dedicated implementations. In Table 2.1 we show the speedup results for the generic version, which can take any full-rank matrix A as input. The dedicated implementa- tions exploited the structure of known input matrices used during OMP testing with further speedup gains, as in resorting to the FFT as part of the internal calculations, for example. Algorithm Seconds Speedup ompQR 617.802467 | ompQRf 360.192118 1.715 ompQRf2 308.379138 1.168 ompQRf3 298.622174 1.032 Table 2.1: Algorithm performance. The speedup column refers to the speedup from the immediately previous implementation. To compute the total speedup from rst to last implementations multiply all speedup values together. Total speedup from ompQR to ompQRf3 is 2.068, which means we doubled the speed of our implementation for the generic matrix version of our code. We used Matlab version R2010b Service Pack 1 to run \experiment.m" which performs many OMP calls on randomized input. The rst improvement came from computing krkkj cos( j)j during the Sweep portion of the algorithm. In this case j is the angle between aj and the residue 10 rk 1. This number re ects how good an approximation zjaj to the residue is, and it is faster to compute than (j). During this step we also kept track of the best approximation to the residue so that during the Update Support stage we could update Sk more e ciently as compared to what was done in ompQR. Finally, we do the sweep only on the set of columns that have not been added to the support set, resulting in further time gains on the Sweep stage whenever k > 1. All these changes where incorporated into ompQRf. For the next round of improvements, we stop building Ak at each iteration as explained in Section 2.2.1. Rather, we initialize Q = In and R = ;, where In 2 Rn n is the identity, and for consistency we de ne A0 = In ; = ;. Subsequently, we update Q and R each time we add a column vector aj of A in the following way. Suppose that at step k > 0 we have a QR decomposition of Ak 1 = QR, and that column ajk is chosen from the Update Support step. Set w = (a T jkQ) T and let H be a Householder re exion such that Hw = v, where v = (#; : : : ;#; 0; : : : ; 0)T has n k zeros after the rst k entries. Then, since HT = H, H2 = In, and HR = R, it is easy to see that QHT RjHTw = Q HTRjH2w = QRjQw = Ak 1jQQTajk = Ak 1jajk = Ak: Therefore, if we set Q0 = QHT, and R0 = RjHTw , we would have found a QR decomposition of Ak = Q0R0 as a function of Q, R, and ajk . The implementation of this update results in faster code compared to the implementation that com- 11 putes a QR decomposition of Ak from scratch for each k. This new approach was implemented in ompQRf2. A nal time improvement came simply from allocating all required variables as opposed to have them grow dynamically as needed. This was implemented in the nal version ompQRf3. 2.3 OMP validation protocol and validation results In this section we present the validation protocol that we followed to verify the correctness of our OMP implementation. 2.3.1 Theoretical results that motivate and justify the protocol The following results provide the foundation for the validation protocol that we adopted. This protocol can be used to validate any OMP implementation. Given a matrix A 2 Rn m with n < m, we can compute its mutual coherence de ned as follows [2]. De nition 1. The mutual coherence of a given matrix A is the largest absolute nor- malized inner product between di erent columns from A. Denoting the k-th column in A by ak, the mutual coherence is given by (A) = max 1 k;j m; k 6=j jaTk ajj kakk2 kajk2 : (2.2) The mutual coherence gives us a simple criterion by which we can test when a 12 solution to Equation 1.1 is the unique sparsest solution available. In what follows, we assume that A 2 Rn m, n < m, and rank(A) = n. Lemma 1. If x solves Ax = b, and kxk0 < 12 1 + 1= (A) , then x is the sparsest solution. That is, if y 6= x also solves the equation, then kxk0 < kyk0. This same criterion can be used to test when OMP will nd the sparsest solution. Lemma 2. For a system of linear equations Ax = b, if a solution x exists obeying kxk0 < 12 1 + 1= (A) , then an OMP run with threshold parameter 0 = 0 is guaranteed to nd x exactly. The proofs of these lemmas can be found or are inspired by results in [2]. In light of these lemmas, we can envision the following roadmap to validate an implementation of OMP. We have a simple uni ed theoretical criterion to guarantee both solution uniqueness and OMP convergence. The following theorem simply uni es the previous lemmas into one statement. Theorem 3. If x is a solution to Ax = b, and kxk0 < 12 1 + 1= (A) , then x is the unique sparsest solution to Ax = b, and OMP will nd it. In light of this result, we can establish the following protocol to validate any implementation of OMP. 2.3.2 Validation protocol Given a full-rank matrix A 2 Rn m, with n < m, compute (A), and nd the largest integer k smaller than or equal to 12 1 + 1= (A) . That is, nd the integer 13 k = 1 2 1 + 1= (A) . Then, build a vector x with exactly k nonzero entries and produce a right hand side vector b = Ax. This way, you have a known sparsest solution x to which to compare the output of any OMP implementation. Pass A, b, and 0 to OMP to produce a solution vector xOMP = OMP(A;b; 0). If OMP terminates after k iterations (or less), and kAxOMP bk < 0, for all possible x and 0 > 0, then the OMP implementation would have been validated. 2.3.3 Validation results Call A = 12 1 + 1= (A) , the constant dependent on A that guarantees the results of Theorem 3 for matrix A. To test our implementation, we ran two experiments involving two random matrices. 1. A1 2 R100 200, with entries in the Gaussian distribution N(0; 1), i.i.d., for which its mutual coherence turned out to be (A1) = 0:3713, corresponding to k = 1 = b A1c. 2. A2 2 R200 400, with entries in the Gaussian distribution N(0; 1), i.i.d., for which its mutual coherence turned out to be (A2) = 0:3064, corresponding to k = 2 = b A2c. We rst note that, with probability 1, Ai, (i = 1; 2), is a full-rank matrix [2]. Second, we would like to mention that for full-rank matrices A of size n m, the mutual coherence satis es (A) p (m n)=(n (m 1)), with the equality being sharp [21]. We used these results to guide us into obtaining matrix A2 for which 14 k = 2 = b A2c > 1. For each matrix Ai, (i = 1; 2), we chose 100 compatible vectors with k nonzero entries whose positions were chosen at random, and whose entries were in the Gaus- sian distribution N(0; 1), i.i.d.. Then, for each such vector x, we built a corresponding right hand side vector b = Aix. Each of these vectors would then be the unique sparsest solution to Aix = b, and OMP should be able to nd them. Finally, given 0 > 0, if our implementation of OMP were correct, it should stop after k steps (or less), and if xOMP = OMP(Ai;b; 0), then kb AixOMPk2 < 0. We ran these experiments for twelve values of 0 equal to 10, 1, 10 1, 10 2, 10 4, 10 6, 10 8, 10 10, 10 12, 10 14, 10 15, and 10 16. For each of these values of 0 we built 100 vectors as described above, with their respective right hand side vectors, both of which were fed to OMP together with the tolerance 0 being tested. We kept track of how many iterations it took OMP to stop, and the value of the norm of the residue kb AixOMPk2 at the end of each run. We mention that our implementation of OMP had as stopping condition that either the residue would be less than the tolerance 0 given, or that n iterations of the main loop would have been executed. Figure 2.2 shows the summary of the results for matrix A1. It contains two graphs, the top graph represents the average of the norm of the residue over the 100 experiments executed for a given tolerance, versus the 12 tolerances chosen. The red line represents the identity in this case. The second graph is the same but for the average number of iterations it took OMP to stop vs the tolerances chosen. The 15 10!20 10!15 10!10 10!5 100 105 10!16 10!14 10!12 10!10 10!8 10!6 10!4 10!2 100 102 !0 ||b !A x om p|| 2 Average residue norm vs tolerance 10?20 10?15 10?10 10?5 100 105 10?1 100 101 102 ?0 # of iteration s Average # of iterations vs tolerance Figure 2.2: OMP behavior for a matrix A with (A) = 0:3713, which corresponds to k0 = 1. 16 red line in this case is the expected number of iterations k = b A1c at stop time. Figure 2.3 is the same as Figure 2.2, but for matrix A2. 10?20 10?15 10?10 10?5 100 105 10?16 10?14 10?12 10?10 10?8 10?6 10?4 10?2 100 102 ?0 ||b ?A x om p|| 2 Average residue norm vs tolerance 10?20 10?15 10?10 10?5 100 105 10?1 100 101 102 103 ?0 # of iteration s Average # of iterations vs tolerance Figure 2.3: OMP behavior for a matrix A with (A) = 0:3064, which corresponds to k0 = 2. One can observe in both cases that there are three modal behaviors of OMP. The rightmost points in each graph correspond to tolerances 0 that are \too large". For them, OMP converges, but it does not have to do much work necessarily, since the default initial solution x = 0 is already close to the right hand side b. The typical behavior corresponds to points in the middle of the graph, they represent the cases when OMP converges in exactly k iterations to the sparsest solution within machine 17 precision. And, nally, the leftmost points represent when OMP fails to converge because the tolerances 0 are too close to machine precision, basically trampling OMP e orts to converge due to roundo and truncation errors. In Figure 2.4 we exempli ed each of the three modal behaviors with three values of 0 typical of each mode. The gure contains three graphs, the top graph is for 0 = 10, the middle graph is for 0 = 10 6, and the bottom graph is for 0 = 10 16. Each of the graphs shows the individual results for each of the 100 experiments run for each tolerance 0. This gure corresponds to matrix A1. In Figure 2.5 we have the same graphs but for matrix A2. We can conclude then that the validation protocol and results con rm that our implementation of OMP is correct. This implementation will return a solution x = OMP(A;b; 0) to Ax = b, within machine precision, whenever the tolerance 0 10 14, and provided kxk0 A. 2.4 OMP testing protocol and results For the rst part of our testing protocol, we set out to reproduce a portion of an experiment described in [2]. The second part deals with studying the image compression properties of multiple matrices A which will be described in more detail in Chapter 3. 18 0 20 40 60 80 100 10?16 10?14 10?12 10?10 10?8 10?6 10?4 10?2 100 102 Experiment ID ||b ?A x om p|| 2 ?0 = 10 0 20 40 60 80 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Experiment ID # of iteration s ?0 = 10 0 20 40 60 80 100 10?18 10?16 10?14 10?12 10?10 10?8 10?6 Experiment ID ||b ?A x om p|| 2 ?0 = 1e?06 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Experiment ID # of teration s ?0 = 1e?06 0 20 40 60 80 100 10?18 10?17 10?16 10?15 10?14 10?13 10?12 10?11 Experiment ID ||b ?A x om p|| 2 ?0 = 1e?16 0 20 40 60 80 100 0 10 20 30 40 50 60 70 80 90 100 Experiment ID # of iteration s ?0 = 1e?16 Figure 2.4: The three modal behaviors, dependent on 0, observed for the matrix A used in Figure 2.2. 19 0 20 40 60 80 100 10?16 10?14 10?12 10?10 10?8 10?6 10?4 10?2 100 102 Experiment ID ||b ?A x om p|| 2 ?0 = 10 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Experiment ID # of iteration s ?0 = 10 0 20 40 60 80 100 10?16 10?15 10?14 10?13 10?12 10?11 10?10 10?9 10?8 10?7 10?6 Experiment ID ||b ?A x om p|| 2 ?0 = 1e?06 0 20 40 60 80 100 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 Experiment ID # of teration s ?0 = 1e?06 0 20 40 60 80 100 10?16 10?15 10?14 10?13 10?12 10?11 Experiment ID ||b ?A x om p|| 2 ?0 = 1e?16 0 20 40 60 80 100 0 20 40 60 80 100 120 140 160 180 200 Experiment ID # of iteration s ?0 = 1e?16 Figure 2.5: The three modal behaviors, dependent on 0, observed for the matrix A used in Figure 2.3. 20 2.4.1 Reproducing previous work In a SIAM Review article by A. M. Bruckstein, D. L. Donoho, and M. Elad [2], the following experiment is presented. We set out to reproduce the portion corresponding to OMP. Consider a random matrix A of size 100 200, with entries independently drawn at random from a Gaussian distribution of zero mean and unit variance, N (0; 1). It can be proven that, with probability 1, every solution for the system Ax = b with less than 51 entries is necessarily the sparsest one possible, and, as such, it is the solution of (P 0), for any > 0. By randomly generating such su - ciently sparse vectors x (choosing the nonzero locations uniformly over the support in a random way, and their values from N (0; 1)), we generate vectors b. This way, we know the sparsest solution to Ax = b, and we shall be able to compare this to the results given by OMP. Since we set to reproduce the results that pertain to OMP in Figure 2 of page 56 of [2], we considered cardinalities in the range of 1 to 70|even though we knew that, with probability 1, only those solutions with cardinality equal to or less than 51 were uniquely the sparsest ones possible|, and we conducted 100 repetitions and averaged the results to obtain the probability of the algorithm nding the solution with which we had generated the right hand side b. Comparing our results with results obtained by published OMP implementations, e.g., like the ones available at SparseLab [19], Figure 2.6 shows that our implementation of OMP reproduces the published experiment, and it performs slightly better than the software found 21 0 10 20 30 40 50 60 70 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Cardinality of the solution Probability of succes s A is 100x200 ompQRf3 SolveOMP Figure 2.6: Reproduction of results in Section 3.3.1 of [2] for OMP. Our implemen- tation ompQRf3 is slightly better at recovering with higher probability the sparsest solution to Ax = b when compared to SolveOMP, an implementation publicly avail- able at [19]. in [19]. 22 Chapter 3 Image representation and compression Image compression plays a central role in modern multimedia communications. Compressed images arguably represent the dominant source of the Internet tra c today, and multiple applications ranging from medical records, to publishing|both in print and online media|, to military imagery, use them. 3.1 Elementary image representation concepts For our purposes an image is a two dimensional sequence of sample values, I[n1; n2]; 0 n1 < N1; 0 n2 < N2; having nite extents, N1 and N2, in the vertical and horizontal directions, respec- tively. The term \pixel" is synonymous with an image sample. The rst coordinate, n1 is understood as the row index, while the second coordinate, n2, is identi ed as the column index of the sample or pixel. The ordering of the pixels follows the canonical ordering of a matrix?s rows and columns [22]. The sample value, I[n1; n2], represents the intensity (brightness) of the image at location [n1; n2]. The sample values will usually be B-bit signed or unsigned 23 integers. Thus, I[n1; n2] 2 f0; 1; : : : ; 2 B 1g for unsigned imagery, I[n1; n2] 2 f 2 B 1; 2B 1 + 1; : : : ; 2B 1 1g for signed imagery. In many cases, the B-bit sample values are best interpreted as uniformly quantized representations of real-valued quantities, I 0[n1; n2], in the range 0 to 1 (unsigned) or 12 to 1 2 (signed). Letting round( ) denote rounding to the nearest integer, the relationship between the real-valued and integer sample values may be written as I[n1; n2] = round(2 BI 0[n1; n2]): (3.1) This accounts for the sampling quantization error which is introduced by rounding the physically measured brightness at location [n1; n2] on a light sensor to one of the allowed pixel values [22, 29]. We will use this framework to represent grayscale images, where a pixel value of 0 will represent \black", and a value of 2B 1 \white". The value of B is called the depth of the image, and typical values for B are 8, 10, 12 and 16. Color images are represented by either three values per sample, IR[n1; n2], IG[n1; n2], and IB[n1; n2] each for the red, green, and blue channels, respectively; or by four values per pixel IC [n1; n2], IM [n1; n2], IY [n1; n2], and IK [n1; n2], each for the cyan, magenta, yellow, and black channels commonly used in applications for color printing, respectively. We will restrict ourselves to grayscale images given that it is always possible to 24 apply a compression system separately to each component in turn [22]. 3.2 Image compression via sparsity In this section, we give an overview of image compression via sparsity. The basic idea is that if Ax = b, b is dense|that is, it has mostly nonzero entries|, and x is sparse, we can achieve compression by storing wisely x instead of b. In speci c, suppose we have a signal b 2 Rn that usually requires a description by n numbers. However, suppose that we can solve problem (P 0), which we now recall from Section 1.2, viz., (P 0) : minx kxk0 subject to kAx bk2 < ; and whose solution x 0 has k nonzeros, with k n, then we would have obtained an approximation b^ = Ax 0 to b using k scalars, with an approximation error of at most . Thus, by increasing we can obtain better compression at the expense of a larger approximation error. We will characterize this relationship between error and compression, or equivalently, error and bit-rate per sample, later on. The choice of matrix A is clearly central to our approach to compression. Inspired by the JPEG and JPEG 2000 standards that use at their core the discrete cosine transform (DCT) and a wavelet transform [25, 22], respectively, we decided to incorporate both transforms in some capacity in our choices of matrix A. 25 3.2.1 Choosing A: \One-dimensional" basis elements Given that the signal that we are going to process comes in the form of a vector b, an inherently one-dimensional object, a rst approach is to consider the one-dimensional DCT waveforms, and any one-dimensional wavelet basis for L2[0; 1]. For the choice of the wavelet basis we opt for the Haar wavelet plus its scaling function, the identity on [0 1]. More speci cally, we know that the one-dimensional DCT-II transform [1, 13], Xk = N 1X n=0 xn cos N n+ 1 2 k ; k = 0; : : : ; N 1; (3.2) has at its core a sampling of the function fk;N(x) = cos x+ 12N k on the reg- ularly spaced set of points S(N) = si 2 [0 1) : si = iN ; i = 0; : : : ; N 1 . With this, we can de ne the vector wk;N = p 2 sgn(k) (fk;N(s0); : : : ; fk;N(sN 1))T 2 RN , which we call generically a \DCT waveform (of wave number k, and length N )." In speci c, we will use DCT waveforms with N = 64 for the one-dimensional compres- sion approach. This is because we subdivide each image in our study database into collections of 8 by 8 non-overlapping sub-images, which are then transformed into vectors b 2 R64, and subsequently compressed, as described above. See Section 3.4. This collection of DCT waveforms is a basis for R64, and we arrange its elements column-wise in matrix form as DCT1 = (w0;64 : : :w63;64) 2 R64 64. Note that all column vectors of DCT1 have the same ?2 norm. The corresponding basis of R64 based on the Haar wavelet is built in the 26 following way. Consider the Haar wavelet?s mother function [10, 14], (x) = 8 >>>>>>>>< >>>>>>>>: 1 if 0 x < 1=2, 1 if 1=2 x < 1, 0 otherwise, (3.3) and its scaling function, (x) = 8 >>>< >>>: 1 if 0 x < 1, 0 otherwise. (3.4) For a natural number n 2 N, build the set of functions Hn = fx 7! n;k(x) = 2 n=2 (2nx k) : 0 k < 2ng; and de ne H 1 = fx 7! (x)g. Note that #Hn = 2n for n 0, #H 1 = 1, and therefore # Sn j= 1Hj = 1 + Pn j=0 2 j = 2n+1. Since 64 = 26, a value of n = 5 will produce 64 functions to choose from in H(n) = Sn j= 1Hj. For each function h 2 H(n = 5), create vector vh;64 2 R64 by sampling h, as before, on the set of points S(N = 64). That is, vh;64 = (h(s0); : : : ; h(s63))T. Observe that kvh;64k2 = kw0;64k2 for all h 2 H(5). We drop the N in v ;N or w ;N when clear from the context. Note that we can order the elements of H(5) in a natural way, namely, h0 = , h1 = 0;0, h2 = 1;0, h3 = 1;1, etc. It is easy to see that the set of vectors fvhjg 63 j=0 is a basis of R64. As with the DCT waveforms, we arrange column-wise these vectors in matrix form as Haar1 = (vh0 : : :vh63) 2 R 64 64. Then, for the one-dimensional 27 approach, we de ne A = [DCT1 Haar1] 2 R64 128, the concatenation of both bases. In Figure 3.1 we can visualize the columns of DCT1 and Haar1 in two di erent ways. A rst method, given a column vector u = (u1; : : : ; u64)T of either basis, is to plot the map j 7! uj. A second method, slightly more elaborate, is to show the two-dimensional mapping of u to c 12 (u), where the invertible map c2 is de ned in Section 3.4. We choose c2 because the ordering of the entries of u induced by this invertible map results in optimal compression performance. See Table 3.2. (a) (b) (c) (d) Figure 3.1: Top row: The rst six waveforms of the DCT1 (a), and the Haar1 (b) normalized bases for R64. Bottom row: Full 2D representation using the inverse map of c2 de ned in Section 3.4 for the DCT1 (c), and Haar1 (d) bases for R64. White corresponds to the maximum value achieved by the basis element, black to the minimum. The intermediate shade of gray corresponds to 0. 28 3.2.2 Choosing A: \Two-dimensional" basis elements Another way to choose a basis for R64 results from taking into account the intrinsic two-dimensional nature of an image and create basis elements that re ect this fact. In speci c, consider the two-dimesional DCT-II used in the JPEG standard [25], Xk1;k2 = N1 1X n1=0 N2 1X n2=0 xn1;n2 cos N1 n1 + 1 2 k1 cos N2 n1 + 1 2 k2 ; (3.5) and, as before, consider the family of functions indexed by k1 and k2, gk1;k2;N1;N2(x; y) = cos x+ 1 2N1 k1 cos y + 1 2N2 k2 ; sampled on all points (x; y) 2 S(N1) S(N2), where in our case, N1 = N2 = 8, and k1; k2 2 f0; : : : ; 7g. (a) (b) Figure 3.2: Full natural 2D representation for the DCT2 (a), and Haar2 (b) bases for R64. White corresponds to the maximum value achieved by the basis element, black to the minimum. The intermediate shade of gray corresponds to 0. 29 Note that gk1;k2;N1;N2(x; y) = fk1;N1(x)fk2;N2(y), and therefore the image of S(8) S(8) under gk1;k2;8;8 can be naturally identi ed with the outer product wk1;8 wk2;8; k1; k2 = 0; : : : ; 7; (3.6) modulo the constants p 2 sgn(k1) and p 2 sgn(k2), that make kwk1k2 = kwk2k2, respec- tively. See Figure 3.2(a) for a graphical representation of each and all of these outer products. There is a total of 64 such outer products, and for each of them we can obtain a vector ewk1;k2;j = c 1 j (wk1;8 wk2;8), where c 1 j is the inverse map of any of the bijections c1, c2, or c3 that take a vector into a matrix, de ned in Section 3.4. It is easy to see that the vector columns of the matrix DCT2;j = (ewk1;k2;j) 2 R 64 64; k1; k2 = 0; : : : ; 7; (3.7) form a basis for R64. The ordering of the column vectors of DCT2;j follow the lexicographical order of the sequence of ordered pairs (k1; k2) for k1; k2 = 0; : : : ; 7 when k2 moves faster than k1, i.e., (0; 0); (0; 1); : : : ; (0; 7); (1; 0); : : : ; (7; 7). In a similar fashion, we can build Haar2;j. In speci c, consider rst the set of functions H(2), which contains 8 functions. Sampling hk 2 H(2) on S(8), we obtain vector vhk;8 2 R 8. Given k1; k2 2 f0; : : : ; 7g we then compute vector evk1;k2;j = c 1 j vhk1 ;8 vhk2 ;8 2 R64; 30 with which, for all k1; k2 2 f0; : : : ; 7g, and following the same order for the ordered pairs (k1; k2) mentioned above, we can create Haar2;j = (evk1;k2;j) 2 R 64 64; k1; k2 = 0; : : : ; 7: (3.8) For a visual representation of the outer products vhk1 ;8 vhk2 ;8, see Figure 3.2(b). Finally, we can de ne A = A(j) = [DCT2;j Haar2;j] 2 R64 128, the concate- nation of both bases. 3.2.3 Some properties of [DCT1 Haar1] and [DCT2;j Haar2;j] From the de nitions of DCT1, Haar1, DCT2;j, and Haar2;j, with j = 1; 2; or 3 (from now on, we won?t make explicitly clear that the index j runs through 1, 2, and 3 as we will assume that this is always the case when it appears,) it is easy to see that they are matrices whose column vectors are pairwise orthogonal and with a common norm, i.e., kak2 = c for any column vector a of any of these matrices. Also, it is easy to see in this case that c = 8. This means that if we were to consider 1 8DCT1, 1 8Haar1, 1 8DCT2;j, and 1 8Haar2;j, separately, they would all be orthogonal matrices, i.e., matrices whose columns are pairwise orthonormal. We present now some properties that are derived from these facts. Lemma 4 (Parseval?s identity). Let U be an orthogonal basis for Rn such that 8u 2 U kuk2 = c. Then 8w 2 Rn X u2U (wTu)2 = c2kwk22: (3.9) 31 Proof. Let w 2 Rn. Since U is a basis for Rn, there exists a linear combination of elements of U = fuk 2 Rn : k = 1; : : : ; ng such that X k xkuk = w; where xk 2 R: Let j 2 f1; : : : ; ng, then we have that X k xkuk = w =) X k xkuTj uk = u T j w; , xjkujk22 = u T j w; since U is an orthogonal basis, , xj = uTj w kujk22 : Hence P k uTkw kujk22 uk = w, but 8j kujk2 = c by assumption. Hence P k u T kw uk = c 2w. From this equation, premultiplying by wT, we obtain that X k uTkw uk = c 2w =) X k uTkw w Tuk = c2wTw; , X k (wTuk)2 = c2kwk22: Lemma 5 (Union of two disjoint orthogonal bases). Let U and V be two orthogonal bases for Rn such that 8u 2 U;v 2 V kuk2 = kvk2 = c, and U \ V = ;. Let W = U [ V , then W is a tight uniform frame with frame bounds A = B = 2c2. Proof. Since 8w 2 W either w 2 U or w 2 V , we must have kwk2 = c. This takes care of uniformity, by De nition 7. See Appendix A for more basic frame de nitions. 32 Let b 2 Rn, then by Parseval?s identity (Lemma 4, Equation 3.9) we have X k (bTuk)2 = c2kbk22; if U = fuk 2 R n : k = 1; : : : ; ng, and X k (bTvk)2 = c2kbk22; if V = fvk 2 R n : k = 1; : : : ; ng: Hence X w2W (bTw)2 = X u2U (bTu)2 + X v2V (bTv)2 X w2U\V (bTw)2 (3.10) = c2kbk22 + c 2kbk22 0; since U \ V = ;; = 2c2kbk22: Let A = B = 2c2. Then from Equation 3.10, Akbk22 = P k(b Tvk)2 = Bkbk22. This establishes the required frame condition for a tight frame [4, 7]. From the proof of Lemma 5, we obtain the following result. Lemma 6 (Upper bound). Let U and V be as in Lemma 5, except that U \ V 6= ;. Then the following inequality holds, 8b 2 Rn X w2U[V (bTw)2 2c2kbk22; and the inequality is tight whenever U \ V 6= U; V . 33 Proof. From Equation 3.10 we have that, for 8b 2 Rn, X w2U[V (bTw)2 = X u2U (bTu)2 + X v2V (bTv)2 X w2U\V (bTw)2; X u2U (bTu)2 + X v2V (bTv)2; since X w2U\V 6=; (bTw)2 0, = c2kbk22 + c 2kbk22; by Lemma 4 = 2c2kbk22: Finally, assume that U \ V 6= U; V , or equivalently U \ V ( U; V . This implies that spanfU \ V g ( Rn, therefore there is a nonzero vector ~b 2 Rn such that ~b ? spanfU \ V g. This implies that X w2U\V 6=; (~bTw)2 = 0; and therefore P w2U[V (~b Tw)2 = 2c2k~bk. We have a similar result for the lower bound in the general case. Lemma 7 (Lower bound). Let U and V be as in Lemma 5, except that U \ V 6= ;. Then the following inequality holds, 8b 2 Rn c2kbk22 kbk 2 2 0 B @2c2 max k~bk2=1 ~b2Rn X w2U\V (~bTw)2 1 C A X w2U[V (bTw)2: Moreover, the second inequality is tight. 34 Proof. Let 0 6= b 2 Rn, then X w2U\V (bTw)2 = kbk22 X w2U\V b kbk2 T w 2 ; kbk22 max k~bk2=1 ~b2Rn X w2U\V ~b k~bk2 T w !2 ; (3.11) kbk22 max k~bk2=1 ~b2Rn X w2U ~b k~bk2 T w !2 ; = kbk22c 2; by Lemma 4: If b = 0 then P w2U\V (b Tw)2 kbk22c 2 trivially since both sides of the inequality are 0. Therefore Inequality 3.11 holds for all b 2 Rn, and we have c2kbk22 kbk 2 2 max k~bk2=1 ~b2Rn X w2U\V ~b k~bk2 T w !2 X w2U\V (bTw)2 (3.12) Combining Equation 3.10 and the Inequalities 3.12 we have that, for 8b 2 Rn, X u2U (bTu)2 + X v2V (bTv)2 X w2U\V (bTw)2 = X w2U[V (bTw)2; ) 2c2kbk22 X w2U\V (bTw)2 = X w2U[V (bTw)2; by Lemma 4; ) c2kbk22 kbk 2 2 0 B @2c2 max k~bk2=1 ~b2Rn X w2U\V (~bTw)2 1 C A X w2U[V (bTw)2: Let b = arg maxb P w2U\V b kbk T w 2 , then Inequality 3.11 becomes an equality and by Equation 3.10 the claim of tightness follows. 35 Combining Lemmas 5, 6, and 7 we obtain the following result, Theorem 8. Let U and V be orthogonal bases for Rn such that 8u 2 U kuk2 = c, and 8v 2 V kvk2 = c. Then W = U [ V is a uniform frame with frame bounds A = 2c2 max kbk2=1 b2Rn X w2U\V (bTw)2; and B = 8 >>< >>: 2c2 if U \ V 6= U; V c2 otherwise : Moreover, if U \ V = ; then W is a tight frame with frame bounds A = B = 2c2. Proof. The only thing left to prove is that B = c2 whenever U \ V = U or U \ V = V . In either case we have that U = V , and the claim follows from Parseval?s identity (Lemma 4, Equation 3.9.) Note that in this case, A = c2 since maxkbk2=1 P w2U\V (b Tw)2 = c2, for the same reason. In subsequent sections, we will conduct a detailed study of the compression properties of both [DCT1 Haar1] and [DCT2;j Haar2;j], for j = 1; 2; 3, as well as matrices derived from them. 3.3 Image database To carry out our experiments and test image compression via sparsity, as well as the properties of the matrices described in Section 3.2 for that purpose, we selected 5 natural images, 4 of them from the University of Southern California?s Signal & Image Processing Institute (USC-SIPI) image database [24]. This database has been widely used for image processing benchmarking. The images are described 36 in Table 3.1 and shown in Figure 3.3. All images are 512 by 512, 8-bit grayscale images, which means they are composed of 5122 = 262,144 pixels that can take integer values from 0 (black) to 255 (white). USC-SIPI le Name Size Notes n/a Barbara 512x512 pixels, 8-bit boat.512 Boat 512x512 pixels, 8-bit elaine.512 Elaine 512x512 pixels, 8-bit 4.2.07 Peppers 512x512 pixels, 24-bit color image Converted to 8-bit grayscale image. 5.2.10 Stream 512x512 pixels, 8-bit Table 3.1: List of images used. Image Peppers was converted to a grayscale image from the 24-bit color original using Photoshop CS3. Image Barbara is not in the USC-SIPI database but we included since it is widely used by the image processing community. 3.4 Methodology Following the approach to image processing at the core of the JPEG image compression standard [25], we subdivide each image in our database in 8 by 8 non- overlapping squares that will be treated individually. Since we need to generate a right hand side vector b to implement our compression scheme via sparsity, cf. Section 3.2, a sub-image Y 2 R8 8 of size 8 by 8 pixels needs be linearized into a vector y 2 R64 to play the role of b. There are many ways to do this, but we tested only three possible approaches. The rst one consisted of concatenating one after the other the columns of Y to form y, we shall call this method c1, which can be thought of as a bijection c1 : R8 8 ! R64 that maps Y 7! y the way just described above. See Figure 3.4(a). 37 (a) Barbara (b) Boat (c) Elaine (d) Peppers (e) Stream Figure 3.3: Images used for our compression algorithms based on sparse image representation. 38 C1 1 2 3 4 1 2 3 4 (a) C2 1 2 3 4 1 2 3 4 (b) Figure 3.4: Two possible ways to vectorize a matrix. (a) Concatenate from top to bottom one after the other, from left to right, the columns of the matrix; or (b) rst ip every other column, and then concatenate the columns as before. This is what functions c1 and c2 do, respectively. Yet another method would be to reverse the ordering of the entries of every even column of Y and then concatenate the columns of the resulting matrix. That is, from Y rst obtain the matrix Y0, where Y 0i;2j = Y8+1 i;2j, and Y 0 ; and Y ; are the entries of Y0 and Y, respectively. Then concatenate one after the other the columns of Y0 to form y0. Call this method c2, also a bijection c2 : R8 8 ! R64 that maps Y 7! y0 the way just described. See Figure 3.4(b). Finally, a third method, c3 : R8 8 ! R64, would traverse an 8 by 8 sub-image Y in a zigzag pattern from left to right, top to bottom, along the diagonal that goes from its top left to its bottom right corners and map it to a vector y00 2 R64, in a similar way as the zigzag mapping shown in Figure 3.5, which depicts this mapping but for a 4 by 4 image that is mapped into a vector in R16. These are the three methods we consider to vectorize a matrix Y representing an 8 by 8 sub-image from any of the images in our test database. Whether we generate a signal vector b 2 R64 by either setting b = c1(Y), b = c2(Y), or 39 C3 1 2 5 9 . . . 12 15 16 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Figure 3.5: How to vectorize a matrix following a zigzag pattern. A 4 by 4 matrix is shown with its 16 elements enumerated from left to right, top to bottom. The path of arrows shows the new induced order after traversing the matrix in a zigzag pattern to obtain a vector. b = c3(Y), we still need to designate a matrix A 2 R64 128 as well to complete the setup. This is where we use [DCT1 Haar1] and [DCT2;j Haar2;j], with j = 1; 2; 3. That is, we will consider A equal to any of these matrices, de ned in Sections 3.2.1 and 3.2.2, respectively. Once we have chosen a way to vectorize a sub-image, say ci|with i either 1, 2, or 3|and a matrix A, we proceed in the following way. Given a tolerance 0 > 0, and an image I that has been partitioned in 8 by 8 non-overlapping sub-images, say fYlg|where l = 1; : : : ; (512=8)2 for the images in our database|we obtain the approximation to yl = ci(Yl) derived from the OMP algorithm, i.e., from the sparse xl = OMP(A;yl; 0), we compute ~yl = Axl, where A = [DCT1 Haar1], or A = [DCT2;j Haar2;j], for j equal to 1, 2, or 3. Using ~yl we can then reconstruct a sub-image by setting eYl = c 1i (~yl). Finally, we can rebuild and approximate the original image I by pasting to- gether, in the right order and position, the set of sub-images feYlg, and form that way eI, the approximate image reconstruction of I. This new image eI is necessarily 40 an approximation, and not the original image I, because in the process we have introduced an error by setting the tolerance 0 > 0, and not 0 = 0. Recall that k~y yk2 < 0, and therefore we will likely have eI 6= I. Since kxlk0 kylk0, and more likely kxlk0 kylk0, storing wisely the set fxlgMl=1, where M is the number of sub-images partitioning I, will provide a compressed representation of image I. How to e ciently create a compressed representation of image I using fxlgMl=1, and the e ects of the choice of tolerance 0, map ci : R8 8 ! R64, and matrix A on such representation will all be addressed in subsequent sections. 3.5 E ects on image reconstruction from choice of map ci Given an image I in our database, we can follow and apply to it the method- ology described in Section 3.4, and obtain at the end of this process a recon- structed image eI from it. In this section we explore the e ects of the choice of map ci : R8 8 ! R64 on the characteristics of image eI for the di erent choices of matrix A that we have selected to study. We summarize the empirical observations obtained from the experiments de- scribed below, and suggest a possible explanation for them without further proof. 3.5.1 Results for A = [DCT1 Haar1] We conducted the following experiments. We set A = [DCT1 Haar1], and chose a tolerance of = 32. Then, for each image I in our database|and each index i = 1; 2;3|we chose map ci : R8 8 ! R64 for the step in the methodology 41 mentioned before that converts a sub-image Yl, from a partition in 8 by 8 sub- images of I, into a vector yl = ci(Yl). The collection of vectors fylgl2L is used to eventually build image eIi, which depends on the choice of map ci this way. See Section 3.4. The characteristics that we used to measure the impact of the choice of map ci are the normalized bit-rate, and the peak signal-to-noise ratio (PSNR). See De ni- tions 3 and 4, respectively. A smaller value for the normalized bit-rate is better than a bigger one given that this implies less bits are necessary to represent the image. A larger value for the PSNR is better than a smaller one as this means the delity of the representation is higher. Table 3.2 summarizes the results of the experiments. From these, we can conclude that for matrix A = [DCT1 Haar1], the choice of c2 over c1 or c3 produces better results. This could be because, if Y = (Yi;j)i;j=1;:::;8 is a natural image, on average jY8;j Y8;j+1j < jY8;j Y1;j+1j for j = 1; : : : ; 7, which makes yc2 = c2(Y) change more slowly than yc1 = c1(Y). By analogy to the be- havior of the DFT, this must translate into needing less column vectors from A to describe, within a certain error , the signal yc2 compared to the number of columns needed to approximate signal yc1 to the same error tolerance. What could explain the superiority of the ordering induced by c1 over c3? For this, we introduce the concept of total variation for a vector v. De nition 2 (Total Variation). Let v 2 Rn. The total variation of v = (v1; : : : ; vn)T is the quantity, V (v) = n 1X i=1 jvi+1 vij: (3.13) 42 Image/Function PSNR (dB) Normalized bit-rate (bpp) Barbara c1: 36.8996 0.1833 c2: 36.9952 0.1863 c3: 36.8470 0.2338 Boat c1: 36.5791 0.1812 c2: 36.6020 0.1608 c3: 36.5615 0.2205 Elaine c1: 36.5003 0.1763 c2: 36.5155 0.1682 c3: 36.4877 0.1885 Peppers c1: 36.7936 0.1193 c2: 36.8674 0.1024 c3: 36.7648 0.1372 Stream c1: 36.4423 0.3161 c2: 36.4686 0.3050 c3: 36.4400 0.3504 Table 3.2: Performance results for A = [DCT1 Haar1] for each c1, c2, and c3. For each image in our test database, we linearized each 8 by 8 sub-image using either c1, c2, or c3. In all cases, the PSNR value was larger using c2; and in all cases, except for image Barbara|although minimally|, the normalized bit-rate was smaller. Both of these measures make c2 a better choice over c1 or c3. The values correspond to runs of OMP with an ?2 termination criteria, and a tolerance = 32. Image Averages Maximums c1 c2 c3 c1 c2 c3 Barbara 684.9177 624.9761 850.5647 3793 3894 4681 Boat 559.6167 471.2769 625.9233 2991 2102 4005 Elaine 552.3401 487.6824 479.7908 2564 1570 2082 Peppers 423.7781 356.2764 405.8694 2707 1940 2430 Stream 978.5122 844.6191 971.4321 3563 2826 2869 Table 3.3: Total variation averages Vavg(I; ci) and maximums Vmax(I; ci) 43 For an image I in our database, create from it the partition of non-overlapping 8 by 8 sub-images fYlgl2L. Since our images are 512 by 512 pixels, we will have 64 64 = 4096 sub-images Yl. That is, #L = 4096. Set i = 1; 2; or 3, and for each and everyone of those sub-images Yl compute the total variation V (ci(Yl)). Then plot in a 64 by 64 image the results to visualize the aggregate of all these values. (a) (b) (c) Figure 3.6: Total variation for image Barbara for the vectorization function (a) c1, (b) c2, and (c) c3. Blue represents low values, red high values. (a) (b) (c) Figure 3.7: Total variation for image Boat for the vectorization function (a) c1, (b) c2, and (c) c3. Blue represents low values, red high values. The images obtained this way are shown in Figures 3.6 through 3.10. Also, 44 (a) (b) (c) Figure 3.8: Total variation for image Elaine for the vectorization function (a) c1, (b) c2, and (c) c3. Blue represents low values, red high values. (a) (b) (c) Figure 3.9: Total variation for image Peppers for the vectorization function (a) c1, (b) c2, and (c) c3. Blue represents low values, red high values. (a) (b) (c) Figure 3.10: Total variation for image Stream for the vectorization function (a) c1, (b) c2, and (c) c3. Blue represents low values, red high values. 45 Table 3.3 contains the maximum value, Vmax(I; ci) = max l2L V (ci(Yl)); (3.14) and the average value, Vavg(I; ci) = 1 #L #LX l=1 V (ci(Yl)); (3.15) for each image I in our database, and each and all of the vectorization functions ci : R8 8 ! R64 that we considered. We observe that usually, but not always, that the smallest value of Vavg predicts the largest PSNR, and the smallest value of Vmax predicts the smallest normalized bit-rate. However, as observed, this is not always the case. If we had only considered c1 and c3 for linearizing functions, the prediction would have failed both for the PSNR and normalized bit-rate measures for image Elaine, for example. This could be because the inherent properties of the matrix chosen to perform the compression, A = [DCT1 Haar1] in this case, may have an impact as well. This is showcased in the following section, where we will see that the total variation of the di erent vectors resulting from di erent choices of linearizing functions ci seems to have less of an impact in the PSNR and normalized bit-ratio characteristics of the reconstructed image when A = [DCT2;3 Haar2;3]. 46 3.5.2 Results for A = [DCT2;j Haar2;j] Proceeding in a similar fashion as in Section 3.5.1, we set A = [DCT2;j Haar2;j] for each j = 1; 2; and 3, and for each such instance of A, we perform the following experiment. Assume that A = [DCT2;j Haar2;j], and set the vectorization function to match the same ordering that was used to create A. This means we pick the vectorization function to be cj. We compute yl = cj(Yl), where|as before|sub- image Yl comes from the partition fYlgl2L of an image I in our image database. Then, continuing with the compression methodology described in Section 3.4, we obtain xl = OMP(A;yl; 0), setting 0 = 32 for this experiment. Finally, from the set of vectors fxlgl2L we obtain ~yl = Axl, and use f~ylgl2L to eventually obtain the reconstructed image eI of our original image I. Again, as in Section 3.5.1, we asses the e ects of the choice of linearizing function ci by the values of PSNR and normalized bit-rate resulting from this representation of I by eI. We show the summary of the results of this experiment for all matrices A considered above, and all images in our database, in Table 3.4. We point out that choosing ci, with i 6= j, when A = [DCT2;j Haar2;j] results in worse values of both PSNR and normalized bit-rate than when i = j, as described above. We report only the results where i = j. Remarkably, any choice of A = [DCT2;j Haar2;j] with a matching vectorization function cj performs better than when A = [DCT1 Haar1] for the normalized bit-rate metric, and almost better for the PSNR metric except for image Stream by 0.0008 47 of a dB. Also, on average, the vectorization order imposed by c3 is slightly better than those by either c1 or c2, although the di erence is practically imperceptible to the human eye. The normalized bit-rate gures all coincide. Image/Function PSNR (dB) Normalized bit-rate (bpp) Matrix Barbara c1 : 37.0442 0.1634 [DCT2;1 Haar2;1] c2 : 37.0443 0.1634 [DCT2;2 Haar2;2] c3 : 37.0443 0.1634 [DCT2;3 Haar2;3] Boat c1 : 36.6122 0.1541 [DCT2;1 Haar2;1] c2 : 36.6120 0.1541 [DCT2;2 Haar2;2] c3 : 36.6120 0.1541 [DCT2;3 Haar2;3] Elaine c1 : 36.5219 0.1609 [DCT2;1 Haar2;1] c2 : 36.5219 0.1609 [DCT2;2 Haar2;2] c3 : 36.5220 0.1609 [DCT2;3 Haar2;3] Peppers c1 : 36.8780 0.0955 [DCT2;1 Haar2;1] c2 : 36.8780 0.0955 [DCT2;2 Haar2;2] c3 : 36.8780 0.0955 [DCT2;3 Haar2;3] Stream c1 : 36.4678 0.2957 [DCT2;1 Haar2;1] c2 : 36.4676 0.2957 [DCT2;2 Haar2;2] c3 : 36.4677 0.2957 [DCT2;3 Haar2;3] Table 3.4: Performance results for A = [DCT2;1 Haar2;1], A = [DCT2;2 Haar2;2], and A = [DCT2;3 Haar2;3] with corresponding vectorization functions c1, c2, and c3, respectively. In all cases, the PSNR and normalized bit-rate values were almost identical, with minor di erences. Matrix A = [DCT2;3 Haar2;3] performs slightly better on average. Mismatching function ci with matrix A = [DCT2;j Haar2;j] when i 6= j, results in degraded performance. The values correspond to runs of OMP with an ?2 termination criteria, and a tolerance = 32. 48 3.6 Normalized bit-rate vs tolerance Following the methodology described in Section 3.4, suppose that we have an image I from our database and let fYlgl2L be a partition of I in #L sub-images of size 8 by 8. Let yl = ci(Yl) for some i = 1; 2; or 3, where ci : R8 8 ! R64 is one of the three maps de ned in Section 3.4. See Figures 3.4 and 3.5. Using OMP, with a full-rank matrix A 2 R64 128, and a tolerance 0 > 0, obtain xl = OMP(A;yl; 0), and compute the approximation to yl given by ~yl = Axl. We know that k~yl ylk2 < 0, and we can count how many nonzero entries there are in xl, namely kxlk0. With this, we can de ne the normalized bit-rate. De nition 3 (Normalized Bit-Rate). Given the context above, the normalized bit- rate measured in bits per pixel (bpp) for image I|given compression matrix A, and tolerance 0|is the number nbr(I;A; 0) = P l kxlk0 N1N2 ; (3.16) where image I is of size N1 by N2 pixels (N1 = N2 = 512 in our case). This is how to interpret this de nition. Suppose that it takes a binary digit or \bit" (0 or 1, for example) to represent a nonzero coordinate in vector xl, and that we|for the moment|ignore how to keep track of the index l, then we will need kxlk0 bits to store or transmit xl. Then, the total number of bits to represent image I is P l kxlk0. The average number of bits per pixel (bpp) is then obtained by dividing this quantity by N1N2. This is the normalized bit-rate. 49 Compare De nition 3 with more traditional forms of compression measures. From [22], the purpose of image compression is to represent the image I of size N1 by N2 and depth B with a string of bits, called the compressed \bit-stream", denoted c. The objective is to keep the length of c, say length(c), as small as possible. In the absence of any compression, we require N1N2B bits to represent the image sample values. The compression ratio is then de ned as cr(I; c) = N1N2B length(c) ; (3.17) and the compressed bit-rate, expressed in bpp, as br(I; c) = length(c) N1N2 : (3.18) The compression ratio is a dimensionless quantity that tells how many times we have managed to reduce in size the original representation of the image, while the compressed bit-rate has bpp units and tells how many bits are used on average per sample by the compressed bit-stream c to represent the original image I. We note that from Equations 3.16 and 3.18 we must have that nbr(I;A; 0) br(I; c) if the bit-stream c is derived from the sparse representation induced by A and 0. The reason for this is two-fold. Firstly, it is unlikely that the coordinates in each of the resulting xl vectors will be realistically represented by only one bit, and secondly, because c would have to somehow include a coding for the indices l for each xl, necessarily increasing the bit count some more. These particular issues, 50 i.e., the number of bits to represent the entries of vectors xl, and the coding of the indices l for each of those vectors into a nal compressed bit-stream c, will be addressed in Section 4. Notwithstanding the above remarks, there is still value in using the normalized bit-stream measure to quantify and plot normalized bit-rate vs tolerance graphs to gauge the compression properties of various compression matrices. Given the results in Table 3.4, we look into the compression properties of matri- ces A = [DCT1 Haar1], and ~A = [DCT2;3 Haar2;3]. We compare these properties for both matrices relative to each other, and to the compression properties of B and C, which are formed from the DCT1 or the Haar1 submatrices of matrix A, respectively. We plot for all images in our database their respective normalized bit-rate vs tol- erance graphs. We took 0 2 T = f2kg11k=0 [ f3; 5; 6; 7; 24; 40; 48; 56; 80; 96; 112g and for each image I in our image database we obtained the corresponding normalized bit-rates nbr(I;A; 0), nbr(I;B; 0), and nbr(I;C; 0) to obtain the plots in Figure 3.11. In Figure 3.12 we compare A = [DCT1 Haar1] with ~A = [DCT2;3 Haar2;3]. We observe that up to a tolerance I , dependent on image I, matrix ~A performs better for tolerance values I . That is, the value of the normalized bit-rate is smaller when performing compression utilizing matrix ~A. For values of I compression with matrix A results in better normalized bit-rate values. We shall see in Section 3.7 that for values of = 32, and smaller, the quality of the image reconstruction is very good. We note from Figure 3.12 that, for all images in our database, I < 32. This means that, for most practical cases, the use of [DCT2;3 Haar2;3] results in slightly smaller normalized bit-rate values than when using [DCT1 Haar1]. 51 (a) Barbara (b) Boat (c) Elaine (d) Peppers (e) Stream Figure 3.11: Normalized bit-rate vs tolerance: \One-dimensional" basis elements. We can observe that for all images the best normalized bit-rate for a given tolerance is obtained for matrix A = [DCT1 Haar1] which combines both the DCT1 and Haar1 bases for R64. 52 (a) Barbara (b) Boat (c) Elaine (d) Peppers (e) Stream Figure 3.12: Normalized bit-rate vs tolerance. We show results to compare the performance of A = [DCT1 Haar1] and ~A = [DCT2;3 Haar2;3]. 53 From the results shown in Figure 3.11, we can see that the DCT1 basis elements perform better compression for any given tolerance than when using the Haar1 basis elements, except for image Stream: when the tolerance is approximately less than 3, the Haar1 basis elements actually result in a smaller normalized bit-rate value. Moreover, and more importantly, combining both the DCT1 and Haar1 bases results in better compression than if either basis is used alone. The same is true for DCT2;3 and Haar2;3. This is very encouraging. However, what can be said of the quality of the reconstructed images? For what range of the tolerances tested is the image quality acceptable? What does an \acceptable" image quality mean? To answer these and related questions, we need to address the issues of image reconstruction and error estimation. 3.7 Image reconstruction and error estimation For the work in this and further sections, unless otherwise noted, we will work with A = [DCT1 Haar1] and the vectorization function c2. Given a 512 512 image I in our database, we can proceed to compress it using the methodology described in Section 3.4. If I is broken down in 8 8 non-overlapping sub-images fYlgl=1;:::;4096, we can obtain for each of them a cor- responding reconstructed sub image eYl = c 12 (~yl), and from those reconstruct an approximation eI to I. Here, ~yl = Axl, xl = OMP(A;yl; 0), and yl = c2(Yl), as before. The compression would come from storing wisely fxlg. We can summarize this procedure with the following notation: eI = rec(I;A; 0). 54 How do we assess the quality of eI when compared to the original image I? We introduce two error estimators|three in actuality, two of them being related. Traditionally, the signal processing community has relied on the peak signal- to-noise ratio, or PSNR [22, 30]. De nition 4 (PSNR). The Peak Signal-to-Noise Ratio between two images eI and I is the quantity, measured in dB, PSNR(eI; I) = 20 log10 0 @ maxI q mse(eI; I) 1 A ; where maxI is the maximum possible value for any given pixel in I, maxI = 2B 1 typically; and mse(eI; I) = 1N1N2 P i;j eI[i; j] I[i; j] 2 is the mean square error between both images. Here N1 and N2 represent the dimensions of I, and I[i; j] represents the value of the pixel at coordinates [i; j] in image I|similarly for eI[i; j]. In our case N1 = N2 = 512, and maxI = 255. PSNR has the advantage that it is easy to compute and has widespread use, but it has been criticized for poorly correlating with perceived image quality [26, 27]. However, in recent years extensive work on other error estimators that take into account the human visual system have arisen. In particular, we present and de ne the structural similarity and mean structural similarity indices [26]. De nition 5 (SSIM). Let eI and I be two images that have been decomposed in L L non-overlapping sub images feYlg and fYlg, respectively. Then the Structural Similarity index for two corresponding sub-image vectorizations, say ~yl = c2(eYl) 55 and yl = c2(Yl), is de ned as follows SSIM(~yl;yl) = (2 ~yl yl + C1)(2 ~ylyl + C2) ( 2~yl + 2 yl + C1)( 2 ~yl + 2yl + C2) ; where yl and yl represent the mean and standard deviation of yl, respectively; and similarly for eyl. The term ~ylyl is the correlation between ~yl and yl. The values C1 and C2 are two small constants. For our purposes, we used the default values of L = 11, C1 = 0:01, and C2 = 0:03 used in [26] when assessing the SSIM of an image in our database and its reconstruction. We used a value of L = 4 when we modi ed OMP to use internally the SSIM as a stopping criteria. More on this later. From the above de nition, we can see that the SSIM index is a localized quality measure that can be represented on a plane that maps its values. It can take values from 0 to 1 and when it takes the value of 1 the two images are identical. In practice, we usually require a single overall quality of measure for the entire image. In that case we use the mean SSIM index to evaluate the overall image quality. De nition 6 (MSSIM). Let eI and I be two images, where the former is the ap- proximation and the later is the original. Then the Mean Structural Similarity index is MSSIM(eI; I) = 1 M MX l=1 SSIM(~yl;yl); where ~yl and yl are the image contents at the l-th local sub-image, and M is the number of local sub-images in the image. 56 Finally, we take a look at the relationship between the size of the sub-image and the tolerance, and how this a ects the quality of the approximation. We analyze the idealized error distribution in which all pixels of the approximation are c units apart from the original. Consider an L L sub image that has been linearized to a vector y of length L2. Assume that the OMP approximation within has distributed the error evenly, that is, if x = OMP(A;y; ) and ~y = Ax, then kAx yk2 < , k~y yk22 < 2; , L2X j=1 ~y(j) y(j) 2 < 2; , L2c2 < 2; , c < L : (3.19) That is, if we want to be within c units from each pixel, we have to choose a tolerance such that c = =L. We note that the least-squares approximation at the core of OMP approx- imates the idealized error distribution. This can be seen in Figure 3.13 where the black dashed line represents this idealized error approximation. For tolerances > 40, we can see that the PSNR for all images considered is above this idealized error distribution. This can be explained by noting that, for example, for = 2048, we would have from Equation 3.19 that c = 2048=8 = 256, but the maximum pixel value is only 255. Therefore, unless the original image I is just a white patch, the initial value of the OMP approximation being an all black image, there are 57 matching pixels in the original and the approximation image eI = rec(I;A; 2048) that are less than 256 units apart. This would necessarily, by De nition 4, im- ply PSNR(eI; I) > 0, a value above the value of the PSNR for the idealized error distribution when = 2048, which is a small negative value. Figure 3.13: Peak Signal-to-Noise Ratio vs tolerance. We observe three typical behaviors for all images. For large values of the tolerance, about > 40, the PSNR of all images is above the PSNR value for the idealized error distribution marked by the black dashed line. This behavior is also observed for very small values of the tolerance, about < 3. Then for values between these two extreme behaviors, all images conform very closely to the idealized error distribution, a fact that is expected from the least-squares approximation at the core of the OMP algorithm. On the other hand, for really small tolerances, about < 3, we observe that the PSNR value for all images jumps again above the PSNR for the idealized error model. This is a happy case when roundo error actually helps. What happens is that for such small tolerances, the roundo to the closest integer for all entries in ~yl = Axl when we form the sub image approximation eYl = c 12 (~yl), coincides with 58 the true value of the pixels in the original sub image Yl. Again, by De nition 4, this increases the value of PSNR(eI; I) compared to the case where roundo would not have taken place. 3.8 PSNR and MSSIM comparison Now that we have some tools to asses the quality of a reconstruction, how do they compare? Figure 3.14: Normalized bit-rate vs MSSIM, PSNR In Figure 3.14 we have plotted the normalized bit-rate versus both error indices MSSIM and PSNR. The rst thing that we observe is that the sensitivity for PSNR varies more dramatically than the sensitivity for MSSIM over the range of tolerances chosen. From Figure 3.15 we can observe that for the range of 20 to 40 dB in PSNR, the MSSIM index ranges from about 0.33 to 0.98. Since a value of 1 in MSSIM corresponds to two identical images, we can focus on values of PSNR no greater than 40 dB in our analysis. Also in Figure 3.15 we corroborate the criticism that 59 Figure 3.15: Peak Signal-to-Noise Ratio vs Mean Structural Similarity has been addressed to PSNR as a measure of image quality. For example, for image Stream at 20 dB we have an MSSIM value of 0.33, whereas for image Peppers we have an MSSIM value of 0.51. A similar wide range between 0.69 (Elaine) and 0.86 (Stream) for MSSIM is observed for 30 dB in PSNR. It is not until 40 dB that we have a much smaller range of MSSIM values, 0.95 (Peppers) to 0.98 (Stream). Therefore, if SSIM and MSSIM capture more accurately the human visual system?s perception of image quality, then the PSNR index is shown to be not so good at it until after values larger than or equal to 35 dB. We therefore drop from the rest of our analysis the PSNR index, other than for an occasional reference point or comparison, and focus on the SSIM and MSSIM indices. We retake the questions at the end of Section 3.6 and answer them with Figure 3.16. From it, if we were to consider desirable values of MSSIM to be above 60 or equal to 0.9, we would see that this would correspond to a tolerance of less than 32 for the Peppers and of 48 for Stream, all other tolerances for the rest of the images falling in between these two values. Figure 3.16: Normalized bit-rate and corresponding MSSIM vs tolerance. In this graph we have plotted together the best normalized bit-rate obtained by combining the DCT and Haar bases, and the corresponding value of the MSSIM index for a given tolerance. The normalized bit-rate graphs are on the bottom left, and the MSSIM index values are above these. This gure combines results for all images. This means that if we wanted all images to have an MSSIM index of 0.9 or better, we would have to pick a tolerance no larger than = 32. This tolerance corresponds, according to Equation 3.19, to a distance of no more than 32/8 = 4 units per pixel between the reconstructed image and the original, on average. Under these circumstances we would achieve a normalized bit-rate of 0.102 to 0.305 bits per pixel. But how good would images with MSSIM greater than or equal to 0.9 actually look? Moreover, what if we could modify OMP as to guarantee a certain 61 minimum MSSIM quality level? It turns out that this modi cation is possible. Consider the following change in the termination condition for the OMP algo- rithm from kAx bk2 < 0 to kAx bkMSSIM MSSIM(c 12 (Ax); c 1 2 (b)) > 0, where 0 is a desired minimum MSSIM index value to achieve in each individual sub image of the reconstruction of I. When we make this change, and recompute the normalized bit-rate vs MSSIM graphs, we obtain the plots shown in Figure 3.17. In this gure, we observe that changing the termination condition for OMP leads to an improvement in the normalized bit-rate without sacri cing image quality. Or, to see this from the opposite perspective, given a normalized bit-rate, we can achieve a better image quality index MSSIM when we use the new stopping criteria. As we shall see in the pictures below, this change redistributes the work that OMP performs more evenly across the image. Finally, to address the question of how good the images actually look, we let the images speak for themselves, and let you|the reader|be the judge. See Figures 3.18 through 3.25. All gures consist of two images, the reconstruction from the original, to the left, and the SSIM index map to the right. The SSIM map represents the localized quality of the image reconstruction. Lighter values are values closer to 1 (\white" = 1), whereas darker values are values closer to 0 (\black" = 0). For each image we obtained a reconstruction for 0 = 32 for the ?2 stopping criteria, and a reconstruction with value of 0 close to the MSSIM from the former for the MSSIM stopping criteria, for comparison purposes. 62 (a) Barbara (b) Boat (c) Elaine (d) Peppers (e) Stream Figure 3.17: Normalized bit-rate vs MSSIM. Comparison of di erent termination criteria for OMP. 63 (a) Barbara (b) SSIM Figure 3.18: Barbara: 0 = 32, PSNR = 36.9952 dB, MSSIM = 0.9444, normalized bit-rate = 0.1863 bpp, termination criteria: k k2. (a) Barbara (b) SSIM Figure 3.19: Barbara: 0 = 0:94, PSNR = 32.1482 dB, MSSIM = 0.9462, normalized bit-rate = 0.1539 bpp, termination criteria: k kMSSIM . 64 (a) Boat (b) SSIM Figure 3.20: Boat: 0 = 32, PSNR = 36.6020 dB, MSSIM = 0.9210, normalized bit-rate = 0.1608 bpp, termination criteria: k k2. (a) Boat (b) SSIM Figure 3.21: Boat: 0 = 0:92, PSNR = 34.1405 dB, MSSIM = 0.9351, normalized bit-rate = 0.1595 bpp, termination criteria: k kMSSIM . 65 (a) Elaine (b) SSIM Figure 3.22: Elaine: 0 = 32, PSNR = 36.5155 dB, MSSIM = 0.9096, normalized bit-rate = 0.1682 bpp, termination criteria: k k2. (a) Elaine (b) SSIM Figure 3.23: Elaine: 0 = 0:90, PSNR = 35.6288 dB, MSSIM = 0.9168, normalized bit-rate = 0.1686 bpp, termination criteria: k kMSSIM . 66 (a) Peppers (b) SSIM Figure 3.24: Peppers: 0 = 32, PSNR = 36.8674 dB, MSSIM = 0.8983, normalized bit-rate = 0.1024 bpp, termination criteria: k k2. (a) Peppers (b) SSIM Figure 3.25: Peppers: 0 = 0:89, PSNR = 35.4309 dB, MSSIM = 0.9080, normalized bit-rate = 0.1011 bpp, termination criteria: k kMSSIM . 67 (a) Stream (b) SSIM Figure 3.26: Stream: 0 = 32, PSNR = 36.4686 dB, MSSIM = 0.9622, normalized bit-rate = 0.3050 bpp, termination criteria: k k2. (a) Stream (b) SSIM Figure 3.27: Stream: 0 = 0:95, PSNR = 34.1398 dB, MSSIM = 0.9634, normalized bit-rate = 0.2853 bpp, termination criteria: k kMSSIM . 68 3.9 Other matrices We have devoted a lot of attention to speci c matrices A to solve problem (P 0), see Equation 1.2, which we reproduce below for convenience, (P 0) : minx kxk0 subject to kAx bk2 < : Namely, we have studied [DCT1 Haar1], and [DCT2;j Haar2;j], for j = 1; 2; 3, de ned in Sections 3.2.1 and 3.2.2, respectively. But there are many other matrices we could potentially use, an in nite number in fact. We brie y present here some results for other matrices that could be of interest but for which further study would be desirable. 3.9.1 \Rotations" of the basis elements of DCT2;3 Consider the matrix DCT2;3 de ned in Section 3.2.2, and recall that its columns are obtained from the outer products of DCT waveforms, see Section 3.2.1, and Equations 3.6 and 3.7. We recall that the two-dimensional DCT-II transform plays a fundamental role in the de nition of the column vectors of matrix DCT2;3, in speci c see Equation 3.5, which we reproduce below for convenience, Xk1;k2 = N1 1X n1=0 N2 1X n2=0 xn1;n2 cos N1 n1 + 1 2 k1 cos N2 n1 + 1 2 k2 : As pointed out in Section 3.2.2, at the core of this transform we nd the family 69 of functions indexed by k1 and k2, gk1;k2;N1;N2(x; y) = cos x+ 1 2N1 k1 cos y + 1 2N2 k2 ; sampled on all points (x; y) 2 S(N1) S(N2), where in our case, N1 = N2 = 8, and k1; k2 2 f0; : : : ; 7g. Recall that S(N) = si 2 [0 1) : si = iN ; i = 0; : : : ; N 1 . Now observe that given that we are fundamentally sampling the function gk1;k2;N1;N2 : R 2 ! R that takes a vector v = (x; y)T 2 R2 and maps it to gk1;k2;N1;N2(x; y) 2 R, we could transform v rst by applying to it a rotation in the R2 plane, and then apply to the transformed vector the function gk1;k2;N1;N2 . Let ( ) be the counter-clockwise rotation by angle in the R2 plane, and identify it with its matrix representation, i.e., ( ) = 0 B B @ cos( ) sin( ) sin( ) cos( ) 1 C C A : Then, we have the following construction for a new matrix A. Let n;N = 2 n N , with n = 0; : : : ; N 1, be N equidistant points in 0; 2 . For each ordered pair (k1; k2) 2 f0; : : : ; 7g f0; : : : ; 7g, and setting vector vi;j = i 8 ; j 8 T 2 S(8) S(8), we can form the matrix, Wk1;k2( n;N) = 0 B B B B B B @ gk1;k2;8;8 ( n;N)v0;0 : : : gk1;k2;8;8 ( n;N)v0;7 ... ... gk1;k2;8;8 ( n;N)v7;0 : : : gk1;k2;8;8 ( n;N)v7;7 1 C C C C C C A 2 R8 8: 70 Using this matrix, form the column vector, wk1;k2;j( n;N) = cj Wk1;k2( n;N) ; where the vectorization functions cj were de ned in Section 3.4, see Figures 3.4 and 3.5. Choosing j = 3, we can compare wk1;k2;3( n;N) with ewk1;k2;3 in Equation 3.7. Basically, when n = 0, we have wk1;k2;3( 0;N) = wk1;k2;3(0) = ewk1;k2;3. We recall that fewk1;k2;3 : k1; k2 = 0; : : : ; 7g form the columns of the DCT2;3 submatrix in A = [DCT2;3 Haar2;3], hence the name for this section. Now, nally, traversing k1 and k2 in lexicographical order for the sequence of ordered pairs (k1; k2) with k1; k2 = 0; : : : ; 7 and k2 moving faster than k1, i.e., (0; 0); (0; 1); : : : ; (0; 7); (1; 0); : : : ; (7; 7), we can form a new matrix, DCT2;3(n;N) = w0;0;3( n;N) : : :w7;7;3( n;N) : Suppose that we choose N = 6, and that we concatenate the six matrices fDCT2;3(n; 6)g5n=0, to form the matrix below, A = [DCT2;3(0; 6) : : :DCT2;3(5; 6)]: (3.20) This is a matrix that uses \rotations"|as introduced above|of some of the column vectors of [DCT2;3 Haar2;3] that could be used in problem (P 0), for example. We have run experiments for comparison purposes for a tolerance = 32 using 71 (a) (b) (c) (d) (e) (f) Figure 3.28: Full natural 2D representation for (a) DCT2;3(0; 6), (b) DCT2;3(1; 6), (c) DCT2;3(2; 6), (d) DCT2;3(3; 6), (e) DCT2;3(4; 6), and (f) DCT2;3(5; 6), bases for R64. White corresponds to the maximum value achieved by the basis element, black to the minimum. The intermediate shade of gray corresponds to 0. 72 A as de ned in Equation 3.20. The results are shown in Table 3.5. For a visual representation of the basis elements of A, see Figure 3.28. Image PSNR (dB) Normalized bit-rate (bpp) MSSIM Barbara 37.1180 0.1325 0.9464 Boat 36.6579 0.1321 0.9234 Elaine 36.5671 0.1359 0.9118 Peppers 36.9145 0.0859 0.9001 Stream 36.5089 0.2482 0.9636 Table 3.5: Performance results for A = [DCT2;3(0; 6) : : :DCT2;3(5; 6)]. In all cases both the PSNR and normalized bit-rate values were better than the values obtained for A = [DCT2;3 Haar2;3], i.e., larger and smaller, respectively. See Table 3.4 for comparison. The values correspond to runs of OMP with an ?2 termination criteria, and a tolerance = 32. 73 Chapter 4 Quantization and coding In 1948 Claude E. Shannon published a seminal paper in two parts in The Bell Systems Technical Journal named \A mathematical theory of communication" [16]. This work marked the dawn of both information theory and coding theory [8, 9]. In it, Shannon de ned precisely what \information" meant mathematically. He drew from work done by his colleagues at Bell Labs, Harry Nyquist and Ralph Hartley, and he certainly was aware of the work on the topic by Norbert Wiener, with whom he had taken a class at MIT. Hartley wrote the following equation to quantify \information" in a discrete setting, H = n log s; where H is the amount of information, n is the number of symbols transmitted, and s is the size of a given alphabet from which the symbols are drawn [5]. Shannon pushed this notion by identifying the amount of information with en- tropy, as Norbert Wiener commented on a ve paragraph review in Physics Today, September 1950, of the book that Shannon and Weaver published together in 1949. The book, \The mathematical theory of communication" [17], contained in a single tome the two-part original paper by Shannon and an expanded and slightly more technical essay that Warren Weaver, the director of natural sciences for the Rock- 74 efeller Foundation, had written earlier for Scienti c American in 1949 [28] about Shannon?s work. More speci cally, in the case of a discrete information source, Shannon repre- sented this as a Markov process, and asked if we could \de ne a quantity which will measure, in some sense, how much information is ?produced? by such a process, or better, at what rate information is produced?" He equated this quantity to entropy, H = nX i=1 pi log2 pi; (4.1) where in Equation 4.1 we have supposed that we have a set of n possible events whose probabilities of occurrence are p1; p2; : : : ; pn. To develop some intuition about this de nition, we review and compare the material in [8] with Shannon?s work [16]. Assume that we are given a random variable X on the nite set f1; 2; : : : ; ng with probability distribution p. The elements X(1) = x1; X(2) = x2; : : : ; X(n) = xn are distinct and p(x1); p(x2); : : : ; p(xn) are nonnegative real numbers with p(x1) + p(x2) + : : :+ p(xn) = 1: We are using pi = p(xi) as a shorthand for prob(X = xi). Now think of the following. The smaller the probability p(xi), the more uncertain we are that an observation of X will result in xi. Thus we can regard 1=p(xi) as a measure of the uncertainty of xi. The smaller the probability, the larger the uncertainty. 75 Shannon thought of uncertainty as information (contrary to common intu- ition.) This is because if an event has probability 1, there is no information gained in asking the outcome of such an event given that the answer will always be the same. Moreover, he also thought the following three properties as desirable for a function that could quantify information: 1. H should be continuous in the pi. 2. If all the pi are equal, pi = 1n , then H should be a monotonic increasing function of n. With equally likely events there is more choice, or uncertainty, when there are more possible events. 3. If a choice can be broken down into two successive choices, the original H should be the weighted sum of the individual values of H. Shannon proved that the only H satisfying these three assumptions is of the form H = K nX i=1 pi log pi: Hence, if we de ne the uncertainty of xi to be log2 1 p(xi) = log2 p(xi); measured in bits, the entropy of the random variable X is de ned to be the expected value H(X) = nX i=1 p(xi) log2 1 p(xi) = nX i=1 p(xi) log2 p(xi); 76 of the uncertainty for the random variable X, i.e., the entropy of X will measure the information gained from observing X. This reasoning establishes Equation 4.1. Moreover, given a communication channel|like a pair of twisted copper for telephony|Shannon identi ed a number called the capacity of the channel and proved, in a nonconstructive way, that arbitrarily reliable communication is possible at any rate below the channel capacity. For example, he de ned the capacity C of the discrete noiseless channel as C = lim T!1 logN(T ) T ; where N(T ) is the number of allowed signals of duration T . The link between entropy and channel capacity was stated by Shannon in the following theorem. Theorem 9 (Shannon, Fundamental Theorem for a Noiseless Channel [16]). Let a source have entropy H (bits per symbol) and a channel have a capacity C (bits per second). Then it is possible to encode the output of the source in such a way as to transmit at the average rate CH symbols per second over the channel where is arbitrarily small. It is not possible to transmit at an average rate greater than CH . However, in reality, all communication channels exhibit some noise. That is, one cannot expect a communication channel to transmit with 100% delity the source message. Shannon included this fact in his schematic diagram for a general communication system, which we reproduce below. Notwithstanding the limitation imposed by noise in the channel, Shannon 77 Figure 4.1: Schematic diagram of a communication system, from [16]. also proved that we can be as sure as we want (but not absolutely sure!) of the information in a received message while transmitting at any rate below channel capacity. The conclusions of Shannon?s theorem are not without some trade-o . In order to communicate reliably at a rate close to channel capacity, we must increase the complexity of our communication scheme [8, 10]. To a large extent, coding theory has been the quest for these good codes. \Quantization", the extraordinary invited paper on the occasion of the 50th anniversary of \A mathematical theory of communication" surveying this quest from an engineering perspective by Robert Gray and David Neuho [6], describes the history of the theory and practice of quantization up until 1998. In this paper the term quantization encompasses transform coding as one of several approaches to exploit redundancy in the source signal. The topic of redundancy in messages was of great interest to Shannon and he studied it in [16] as well. In transform encoding, the source samples are collected into a vector of, say, dimension n that is multiplied by an orthogonal matrix (an orthogonal transform) and the resulting transform coe cients are scalar quantized. A scalar quantizer can be de ned as consisting of a set of intervals or cells 78 S = fSi R : i 2 Ig that form a partition of the real line, where the index set I is ordinarily a collection of consecutive integers beginning with 0 or 1, together with a set of reproduction values or points or levels C = fyi 2 R : i 2 Ig so that the overall quantizer q is de ned by q(x) = yi for x 2 Si, which can be expressed concisely as q(x) = X i2I yi1Si(x); (4.2) where the indicator function 1S(x) is 1 if x 2 S and 0 otherwise [6]. More generally, a vast class of memoryless quantizers can be described as follows. A quantizer of dimension k, a positive integer, takes as input a vector x = (x1; : : : ; xk)T 2 A Rk. Memoryless refers to a quantizer which operates independently on successive vectors. The set A is called the alphabet and it is often called the support of the source distribution. If k = 1 the quantizer is scalar, otherwise it is vector. The quantizer consists then of three components|a lossy encoder : A ! I, where the index set I is an arbitrary countable set, usually taken as a collection of consecutive integers, a reproduction decoder : I ! A^, where A^ Rk is the reproduction alphabet, and a lossless encoder : I ! J , an invertible mapping (with probability 1) into a collection J of variable-lenght binary vectors that satis es the pre x condition, that is, no vector in J can be the pre x of any other vector in the collection [6]. Alternatively, a lossy encoder is speci ed by a partition S = fSi R : i 2 Ig of A; a reproduction decoder is speci ed by a codebook C = f (i) 2 A^ : i 2 Ig of points, codevectors, or reproduction codewords, also known as the reproduction 79 codebook; and the lossless encoder can be described by its binary codebook J = f (i) : i 2 Ig containing binary or channel codewords. The quantizer rule is the function q(x) = ( (x)) or, equivalently, q(x) = (i) whenever x 2 Si [6]. This type of quantizers are known as \vanilla" vector quantizers or block- source codes. In the literature the term code is used as a generic substitute for quantizer. The instantaneous rate of a quantizer applied to a particular input is the normalized length r(x) = 1k l( ( (x))) of the channel codeword, the number of bits per source symbol that must be sent to describe the reproduction. If all binary codewords have the same length, we talk of a xed-length or xed-rate quantizer. To measure the quality of the reproduction, we assume the existence of a nonnegative distortion measure d(x; x^) which assigns a distortion or cost to the reproduction of input x by x^. Ideally, one would like a distortion measure that is easy to compute, useful in analysis, and perceptually meaningful in the sense that small (large) distortion means good (poor) perceived quality. No single distortion measure accomplishes all three goals [6]. However, d(x; x^) = kx x^k22 satis es the rst two. We also assume that d(x; x^) = 0 if and only if x = x^. In this light we say that a code is lossless if d(x; ( (x))) = 0 for all inputs x, lossy otherwise. Finally, the overall performance of a quantizer applied to a source is charac- 80 terized by the normalized rate R( ; ) = E[r(X)] = 1 k E[l( ( (X)))] = 1 k X i l( (i)) Z Si f(x) dx; (4.3) and the normalized average distortion D( ; ) = 1 k E[d(X; ( (X)))] = 1 k X i Z Si d(x;yi)f(x) dx: (4.4) Every quantizer ( ; ; ) is thus described by a rate-distortion pair (R( ; ); D( ; )). The goal of compression system design is to optimize the rate-distortion trade-o [6]. In light of the results by Shannon, compression system design will also have to take into account the characteristics of the communication channel in managing the rate-distortion trade-o . Also, from the de nitions above, it is clear that knowledge of the probability density function of the source messages is very important. 4.1 Image quantization and coding Under the perspective from the previous introduction, we return to the topic of image quantization and coding. The image standards JPEG and JPEG 2000 are framed in the transform coding paradigm, and contain two more steps besides their respective discrete cosine transform (DCT) and the Cohen-Daubechies-Feauveau 81 5/3 (lossless) and 9/7 (lossy) biorthogonal wavelet transforms at the core of their compression engines. Both have a quantization and an coding step, as described in the literature, with their respective \inverses", to complete their de nitions [3, 22, 25]. Note how the use of the words \quantization" and \coding" refer to and from the previous introduction, respectively. The general schematic for transform coding, under which both JPEG and JPEG 2000 techniques|as well as our approach|fall under, is described in the schematic drawing shown in Figure 4.2. T Q E T' Q' E' b b' Storage/ Transmission Figure 4.2: Schematic diagram of a general transform coding system. The transform T is meant to exploit the redundancy in the source signal b and decorrelate it. It has an inverse T 1, or at the very least a left inverse T 0 such that T 0Tb = b. In our approach we have A = T 0, and T is de ned via OMP by Tb = OMP(A;b; 0). In this case, we have kT 0Tb bk2 < . Therefore, our compression scheme is lossy. Q is a non-invertible scalar quantizer that will be applied to the coe cients of vector Tb, and Q0 is its reproduction function. Finally, we have an invertible lossless encoder E, in the introduction, with E 0 = E 1. The function composition QT is equivalent to the lossy encoder , and the function composition T 0Q0 corresponds to the reproduction decoder from the introduction 82 above. In order to describe the overall performance of our quantizer (QT;E; T 0Q0), we would need to characterize the rate-distortion pair (R(QT;E); D(QT; T 0Q0)). For the scope of the present work, we have studied the error sensitivity for Ax = y, in terms of A 2 Rn m and x 2 Rm. Assume that the columns aj of A = (aj) all have the same norm kajk2 = c. Let ~x = x + , where = ( 1; : : : ; m)T. What can we say of the distance between ~y = A~x and y = Ax? k~y yk2 = kA~x Axk2 = kA(x + ) Axk2 = kA k2 = mX j=1 aj j 2 = X j 6=0 aj j 2 X j 6=0 kaj jk2 = X j 6=0 kajk2j jj = X j 6=0 cj jj ; kajk2 = c 8j X j 6=0 ck k1 ; k k1 = max j j jj = ck k1k k0 ; k k0 = #fj : j jj > 0g: (4.5) Hence, the error k~y yk2 is bound by ck k1k k0. Moreover, the value of k k0 is linked to the sparsity of x, because in our case the error comes from scalar quantizing the entries of x. That is, if|for example|~x = round(x), where ~x is the vector whose entries are exactly those of x but rounded to the closest integer, then 83 necessarily k k0 = k~x xk0 kxk0 (4.6) Hence, in the case where a scalar quantization scheme satis es (4.6), Equation 4.5 gives k~y yk2 ck k1kxk0; (4.7) with A = (aj), ~y = A~x, y = Ax, and c = kajk2 for all j. Observe that, in particular, when ~x = round(x), we have k k1 = 1=2. From Equation 4.7 we see that the error in the reconstruction due to scalar quantization is linked to the size of c, k k1, and the sparsity of x. We are tempted to make c as small as possible, or modify the quantization scheme to make k k1 smaller to reduce the reconstruction error due to scalar quantization. The problem is that the magnitude of c is linked to the norm of x, which a ects the number of bits we would need to represent x. To see the link between c and the norm of x, consider the following example. We will come back to the issue of the representation of x shortly. Assume that kajk2 = c for all j, where A = (aj), and aj are the columns of A. Suppose that Ax = y and x = (x1; 0; : : : ; 0)T, i.e. y = x1a1. Keeping y xed, as a given right hand side, what happens to kxk2 = jx1j as we modify c? Let a > 0 and consider Aa = (aaj). Then Aax = x1(aa1) = a(x1a1) = ay, hence Aa(a 1x) = y. It is easy to see by linearity from this example that the general case, where Ax = y and x = (x1; : : : ; xm)T, also implies Aa(a 1x) = y. The norm of the columns of Aa is kaajk2 = akajk2 = ac, a > 0. Making ac small implies making 84 the parameter a small, since c is xed by the choice of A, and therefore making ka 1xk2 = a 1kxk2 big. Hence, if a vector x0 solves Ax = y, then the norm of a 1x0, which solves Aax = y, will increase as a ! 0. It is easy to see that if x0 is the sparsest vector that solves Ax = y, then a 1x0 is the sparsest vector that solves Aax = y, and kx0k0 = kax0k0. From the above discussion we conclude that the norm of a 1x0 has an inversely proportional relationship with the size of ac. Therefore, if we are to use a nite predetermined number of bits to represent a 1x0, the solution of Aax = y, we need to choose a carefully. Generalizing Equation 4.7 to the class of matrices Aa, a > 0, the error bound due to the scalar quantization of the entries of x becomes k~y yk2 ack k1kxk0; (4.8) where A = A1 = (aj), kajk2 = c for all j, ~y = Aa~x, y = Aax, and ~x = x + with k k0 kxk0. From the results of Section 3.2.3 we can see that the magnitude of the co- ordinates of x are bound by a multiple of kyk2. This also has an impact on how many bits would be needed to represent x. Therefore the choice of a has to take into consideration the maximum value that the value of kyk2 will impose on the magnitude of the coordinates of x. Finally, let us recall that our image compression approach comes from the 85 OMP solution x0 to problem (P 0 0 ) for a given matrix A = (aj) whose column vectors satisfy kajk2 = c, a vector b, and a tolerance 0 > 0 (such that kAx0 bk2 < 0.) Then, choosing a > 0, and following the description of T at the beginning of this section, if we set x0 = Tb = OMP(aA;b; 0) for a given signal b and a tolerance 0 > 0, aA = T 0, Q a scalar quantizer that satis es Inequality 4.6, and Q0 its corresponding reproduction function, the triangle inequality and Inequality 4.8 give d( ( (b));b) = kT 0Q0QTb bk2 = kT 0Q0QTb T 0Tb + T 0Tb bk2 = kaA~x0 aAx0 + aAx0 bk2 kaA~x0 aAx0k2 + kaAx0 bk2 = ack k1kx0k0 + 0; where = ~x0 x0. This inequality would give us a footing in the computation of the normalized average distortion D( ; ), see Equation 4.4. From the de nition of D( ; ), it is clear that we need to know something about the probability density function of the input sources, i.e., the statistics of the 8 8 linearized sub images into which each image is partitioned, if we are to compute D. As a surrogate of such knowledge, we can observe the distribution of the coe cients for each of the vectors resulting from the analysis of the images in our database, and their corresponding histograms. This is what Figures 4.5 through 4.9 show. For each image I in our database, we used matrix A = [DCT2;3 Haar2;3] with 86 a tolerance of = 32 to compute its statistics. On the x-axis of each sub gure we have matched column ai at position i to integer i for each of the columns of A = (ai). Hence, positions 1 to 64 correspond to the DCT waveforms, and positions 65 to 128 to the Haar waveforms. A sub gure labeled (a) corresponds to the histogram for the frequency each column vector of A is chosen. For example, since there are 4096 sub images of size 8 8 in a 512 512 image, column vector a1 will be chosen 4096 times since it corresponds to the constant vector, which computes the mean or DC component for each sub image. Sub gure (b) corresponds to a partial representation of the distribution of the coe cients that multiply each and every column vector whenever such vector is chosen in the representation/approximation of an input b. For example, suppose that column a74 was multiplied by a coe cient a74 = 3:2310 to obtain the representation of some input b = a74a74 + r within a tolerance of = 32; then we would have plotted point (74; 3:2310) in sub gure (b). We have not plotted the coe cients for a1 since they correspond to the DC components of the sub images of I, which vary between 0 and 255. We note that all images in our database have a similar structure. For comparison purposes, we obtained the histogram and the distribution of coe cients for a randomly generated image with a uniform distribution on [0 255], see Figure 4.10. The rst thing we note is that unlike for the natural images in our database, all column vectors|except a1 and a65, which correspond to the constant vectors (one for the DCT and one for the Haar waveforms)|are chosen about the same number of times regardless of their position; and the distribution of the values of the coe cients is uniform as well. It is also clear that this is the image that 87 requires the most nonzero coe cients to be reconstructed within the tolerance = 32 chosen. This is consistent with the de nition of information by Shannon: the more the uncertainty in a source, the more there is information in it. Regarding the computation of the rate R( ; ) for our image quantizer, we would have to choose a lossless encoder , which we have not done here. On the other hand, we have plotted the distortion as measured by the MSSIM index and the PSNR versus the idealized normalized bit-rate, de ned in Section 3.6. This bit-rate is unattainable in practice, as observed in that section, but this gives us nonetheless an idea of an upper bound in the rate-distortion trade-o , and lets us compare how much room we have to select a lossless encoder to complete the implementation of a quantizer using our sparse image representation approach. Figure 4.3 shows the MSSIM and PSNR versus normalized bit-rate trade-o . Figure 4.4(a) shows the PSNR versus normalized bit-rate trade-o for image Lena, which we computed to compare with Figure 4.4(b) which shows results for that image for three published fully implemented quantizers. We observe that there is enough room to pick an encoder that could compete with these implementations. 88 (a) (b) Figure 4.3: (a) MSSIM vs Normalized bit-rate, (b) PSNR vs Normalized bit-rate. 89 (a) (b) Figure 4.4: PSNR vs bit-rate: (a) Normalized bit-rate results for A = [DCT1 Haar1] prior to any coding, and (b) bit-rate coding performances published in [14] for image Lena: Said and Pearlman?s SPIHT algorithm [15], Embeded Coding and the Wavelet-Di erence-Reduction compression algorithm (\new algorithm") [23], and Shapiro?s EZW algorithm [18]. 90 (a) (b) Figure 4.5: Histograms for image Barbara. A = [DCT2;3 Haar2;3], tolerance = 32. 91 (a) (b) Figure 4.6: Histograms for image Boat. A = [DCT2;3 Haar2;3], tolerance = 32 92 (a) (b) Figure 4.7: Histograms for image Elaine. A = [DCT2;3 Haar2;3], tolerance = 32 93 (a) (b) Figure 4.8: Histograms for image Peppers. A = [DCT2;3 Haar2;3], tolerance = 32 94 (a) (b) Figure 4.9: Histograms for image Stream. A = [DCT2;3 Haar2;3], tolerance = 32 95 (a) (b) Figure 4.10: Histograms for uniform random input. A = [DCT2;3 Haar2;3], tolerance = 32 96 Chapter 5 New directions From the discussion in Chapter 4, it is clear that to complete an implemen- tation of a quantizer based in the sparse image representation we would have to choose a lossless encoder that would allow us to match or beat the results shown in Figure 4.4(b). Such an encoder would have to be able to take advantage of the sparsity of the solution x0 to problem (P 0) : minx kxk0 subject to kAx bk2 < ; for a given right hand side b and tolerance > 0. A combination of techniques such as Hu man encoding and run-legth encoding come to mind, see for example [10] for details. Said and Pearlman?s SPIHT algorithm [15], and Shapiro?s EZW algorithm [18] would be sources of inspiration as well. As noted in Chapter 4, knowledge of natural image statistics is essential to the computation of rate R and distortion D quantities. This knowledge can also be used to optimize the design of the scalar quantizer in the setting of transform encoding by allocating more cells where most of the coordinates xi of solution vector x0 = (x1; : : : ; xm)T for (P 0) would fall, as illustrated by Figures 4.5 through 4.9 Moreover, as pointed out in Section 3.9 at the end of Chapter 3, a study of other matrices A other than the matrices studied in this work would be desirable. A 97 promising approach would be to study the properties of Gabor systems, for example. Finally, given that the NP-completeness of the solution to problem (P 0) in the general case [12] does not preclude a priori the existence of fast and e cient algorithms for the computation of sparse solutions to a restricted class of matrices A, it would be most desirable to discover such algorithms and identify the class of matrices for which they would exist, if at all. 98 Appendix A Basic frame de nitions De nition 7 (Frame). A frame for a Hilbert space H is a sequence of vectors fxig H for which there exists constants 0 < A B < 1 such that, for every x 2 H, Akxk2 X i jhx; xiij 2 Bkxk2: (A.1) The constants A and B are known respectively as lower and upper frame bounds. We will often presume without stating it that A is the greatest lower frame bound and B is the least upper frame bound. A frame is called a tight frame if the optimal upper and lower frame bounds are equal; A = B. A frame is a Parseval frame if A = B = 1. A uniform frame is a frame in which all the vectors have equal norm [7]. Equation A.1 is usually called the frame condition. 99 Bibliography [1] William L. Briggs and Van Emden Henson. The DFT, an owner?s manual for the discrete Fourier transform. SIAM, 1995. [2] Alfred M. Bruckstein, David L. Donoho, and Michael Elad. From sparse solu- tions of systems of equations to sparse modeling of signals and images. SIAM Review, 51(1):34{81, 2009. [3] Charilaos Christopoulos, Athanassios Skodras, and Touradj Ebrahimi. The JPEG 2000 still image coding system: an overview. IEEE Transactions on Consumer Electronics, 46(4):1103{1127, November 2000. [4] R. J. Du n and A. C. Schae er. A class of nonharmonic Fourier series. Trans- actions of the American Mathematical Society, 72(2):341{366, 1952. [5] James Gleick. The information: a history, a theory, a ood. Pantheon Books, 2011. [6] Robert M. Gray and David L. Neuho . Quantization. IEEE Transactions on Information Theory, 44(6):2325{2383, October 1998. [7] Deguang Han, Keri Kornelson, David Larson, and Eric Weber. Frames for undergraduates, volume 40 of Student Mathematical Library. American Math- ematical Society, 2007. [8] W. Cary Hu man and Vera Pless, editors. Handbook of coding theory, volume 1, chapter 1. Elsevier Science B. V., 1998. [9] W. Cary Hu man and Vera Pless. Fundamentals of error-correcting codes. Cambridge University Press, 2010. [10] St ephane G. Mallat. A wavelet tour of signal processing. Academic Press, 1998. [11] St ephane G. Mallat and Zhifeng Zhang. Matching pursuits with time-frequency dictionaries. IEEE Transactions on Signal Processing, 41(12):3397{3415, De- cember 1993. [12] Balas Kausik Natarajan. Sparse approximate solutions to linear systems. SIAM Journal on Computing, 24(2):227{234, 1995. [13] K. R. Rao and P. Yip. Discrete cosine transform: algorithms, advantages, applications. Academic Press Professional, Inc., San Diego, CA, USA, 1990. [14] Howard L. Resniko and Raymond O. Wells, Jr. Wavelet analysis. The scalable structure of information. Springer-Verlag, 1998. Corrected 2nd printing. 100 [15] A. Said and W. A. Pearlman. A new, fast, and e cient image codec based on set partitioning in hierarchical trees. IEEE Transactions on Circuits and Systems for Video Technology, 6(3):243{250, June 1996. [16] Claude E. Shannon. A mathematical theory of communication. The Bell System Technical Journal, 27:379{423, 623{656, July, October 1948. [17] Claude E. Shannon and Warren Weaver. The mathematical theory of commu- nication. The University of Illinois Press: Urbana, 1949. [18] Jerome M. Shapiro. Embedded image coding using zerotrees of wavelet coe - cients. IEEE Transactions on Signal Processing, 41(12):3445{3462, December 1993. [19] SparseLab. Seeking sparse solutions to linear systems of equations: Software tools. http://sparselab.stanford.edu/. [20] G. W. Stewart. Introduction to matrix computations. Academic Press, 1973. [21] T. Strohmer and R. W. Heath. Grassmanian frames with applications to coding and communication. Appl. Comput. Harmon. Anal., 14(3):257{275, 2003. [22] D. S. Taubman and Michael W. Marcellin. JPEG 2000: image compression fundamentals, standards and practice. Kluwer Academic Publishers, 2001. [23] J. Tian and R. O. Wells, Jr. A lossy image codec based on index coding. In IEEE Data Compression Conference, DCC ?96, page 456, 1996. [24] Viterbi School of Engineering, University of Southern California. The USC-SIPI image database. http://sipi.usc.edu/database/. [25] Gregory K. Wallace. The JPEG still picture compression standard. Communi- cations of the ACM, 34(4):30{44, 1991. [26] Zhou Wang, Alan C. Bovik, Hamid R. Sheikh, and Eero P. Simoncelli. Image quality assessment: From error measurement to structural similarity. IEEE Transactions on Image Processing, 13(1):1{14, 2004. [27] Andrew B. Watson, editor. Digital images and human vision. The MIT Press, 1993. [28] Warren Weaver. The mathematics of communication. Scienti c American, 181(1):11{15, July 1949. [29] Mladen V. Wickerhauser. Adapted wavelet analysis from theory to software. A K Peters, Ltd., 1996. [30] Wikipedia. Peak signal-to-noise ratio. http://en.wikipedia.org/wiki/Peak signal- to-noise ratio. 101