{'fraction_non_alphanumeric': 0.03577003121748179, 'fraction_numerical': 0.009446540062434963, 'mean_word_length': 4.604100227790433, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Image Compression Based on Compressive Sensing: End-to-End Comparison with JPEG

We present an end-to-end image compression system based on compressive sensing. The presented system integrates the conventional scheme of compressive sampling and reconstruction with quantization and entropy coding. The compression performance, in terms of decoded image quality versus data rate, is shown to be comparable with JPEG and significantly better at the low rate range. We study the parameters that influence the system performance, including (i) the choice of sensing matrix, (ii) the trade-off between quantization and compression ratio, and (iii) the reconstruction algorithms. We propose an effective method to jointly control the quantization step and compression ratio in order to achieve near optimal quality at any given bit rate. Furthermore, our proposed image compression system can be directly used in the compressive sensing camera, e.g. the single pixel camera, to construct a hardware compressive sampling system.


I. INTRODUCTION
Multimedia is emerging due to advanced developments of social networks and the mobile Internet. Billions of image and video resources are transmitted through miscellaneous ways on the Internet everyday. Limited by the bandwidth and storage, image and video compression still play a pivotal role in multimedia. Meanwhile, Compressive Sensing (CS) [2]- [4] has been proposed for more than a decade as a method for dimensionality reduction of signals that are known to be sparse or compressible in a specific basis representation. By "sparse" we mean that most coefficients of the representation are zero; "compressible" indicates that the magnitude of the coefficients decays quickly according to a power law, hence the signal can be well approximated by a sparse signal. In the CS paradigm, the signal is projected onto a low-dimension space, resulting in a measurements vector. For a sparse signal, it is possible to exactly reconstruct it from the measurements vector. If the measurements are noisy or the signal is not sparse but compressible, the reconstruction yields an approximation to the original signal. Natural images are inherently compressible in the frequency or wavelet domain and therefore suitable for CS. The past few years saw impressive progress in this field with new reconstruction algorithms [5], [6] achieving better was investigated in [15], [16], and methods of accounting for the quantization effect (including the extreme case of a 1-bit quantizer) in the reconstruction algorithm were studied in [15], [17]- [22]. However this significant body of research was of limited value for our purposes. First, these works assumed a random, or a structurally random sensing matrix, while the sensing matrices suitable for our purposes could be different (see Sec. II-A). Second, most of these works did not assume any lossless coder and therefore did not consider the resulting bit rates. Those that did consider a lossless coder did not study the interaction between the quantizer and the lossless coder and the trade-offs in their design.
In the signal compression field, a "compression system" converts the signal to a bit stream at the encoder and reconstructs the signal from the bit stream at the decoder. The performance of such a system is measured in terms of the reconstruction error vs. bit rate. There are various papers reporting the effect of measurement quantization on the reconstruction error [5], [13], [18]; however, these studies present their results in terms of reconstruction error vs. compression ratio, rather than vs. bit rate. This is because they do not consider the quantization and lossless coding aspects, that is, converting the sequence of quantization codewords into a low rate bit stream. This is the main reason there has been no comparison of performance between CS and classical image compression methods such as JPEG and JPEG2000 -without a bit rate there is no way of meaningful comparison. While we were preparing this work for publication we came across a deep learning inspired paper by Cui et al. [23] (among other deep learning based methods), which also describes an endto-end image compression system based on CS 2 .

B. JPEG, JPEG2000 and Uniqueness of CS-based Image Compression
Since its introduction in 1992, JPEG [24] has been one of the most popular image compression methods in use. JPEG2000 [25], which was released ten years later, has superior compression properties at a wider bit rate/quality range. A 2 The description in [23] was not detailed enough to allow replication and comparison in our lab, but we highlight the main similarities and differences in the system design relative to our work. 1) Cui et al., like us, recognize that a random, or structurally random sensing matrix will not perform well in this application, and they use a neural network to train a sensing matrix. The histogram of measurements shows that most of the measurement energy is concentrated in a just a few measurements, similar to the effect of classical transforms like DCT or WHT, but without the computational advantages of the fast transforms. 2) Cui et al.'s system, like ours, uses uniform quantization and arithmetic coding. Their system does seem to jointly optimize the compression ratio and the quantization step, but it has a neural network sub-module for improving the dequantization. The use of this module results in a modest improvement in PSNR about 0.4 dB.
3) The reconstruction in Cui et al.'s system is done by a neural network, which is trained jointly with the sensing matrix and the dequantization correction sub-network. Consequently, unlike our system, their encoder and decoder seem to be tightly coupled. 4) Cui et al. 's system works on blocks of 32x32 samples, unlike our system which operates on the whole image. A second decoding stage (i.e.a second neural network) is used to eliminate blocking artifacts. This blocking may be because full frame training would require prohibitively large amount of training materials and computational resources. 5) Cui et al. compare their results with other CS-based compression systems to which they added fixed length lossless coding. There is no comparison with standard image compression methods like JPEG or JPEG2000 as performed in this paper. brief overview of these standards is given in Sec. II-F, where their architectures are compared with the proposed framework, dubbed compressive sensing based image compression system (CSbIC).
HEIF [26], [27] (the still-image version of HEVC) and WebP [28] are recently introduced image compression standards based on intra-frame prediction, which claim to have better compression than JPEG. Another recent approach is sparse representation image coding, after JPEG2000, where the wavelet basis is replaced by a learned, over-complete dictionary with constraints on the variance of the representation coefficients [29]- [31]. However, direct comparison of dictionary learning based methods with other methods is difficult, since the dictionaries are learned for a specific class of images.
At the time of its introduction, the higher complexity of JPEG2000 was an impediment for its adoption, but as computing capabilities improved, this was no longer an issue, and this is true for the aforementioned compression enhancements. Nevertheless, JPEG remains the image compression tool of choice in a wide range of applications. It appears that for most conventional applications (web, professional and amateur photography, publishing, etc.), the advantages offered by the more advanced compression schemes, were not significant enough to justify migration to a new standard. This raises the question whether there is need for yet another image compression scheme. In response, a compression scheme based on CS has some unique properties which are radically different from those of any conventional signal compression scheme.
• The encoder is of much lower complexity than the decoder. In fact, a significant part of the encoder processing, namely the measurements generation, can be done in the analog domain via a single-pixel camera or lensless camera [32]- [34]. • The decoder is not fully defined by the encoder. Different reconstruction algorithms can be applied to the same measurements vector, so as more prior information becomes available the reconstruction algorithm can be improved to take it into account [5], [6]. • The loss of a small fraction of the measurements generally results in only a minor degradation in reconstruction quality [16], [18]. This may be used to achieve robustness to channel impairments without the usual overhead of error protection. These unique properties suggest that CSbIC would be very useful in new, non-conventional applications such as media sensor networks, where the sensors need to be inexpensive and have low energy consumption and the network operates at low power and is usually unreliable. In contrast, the control center has ample processing power for decoding, and much prior information may be gathered while the system is running. Therefore, CSbIC is suitable for a class of applications that are different from those that JPEG and JPEG2000 were designed for.

C. Contributions and Organization of This Work
This paper makes the following contributions:  i) We provide an end-to-end architecture for a CSbIC system, including compressive sensing, quantization, lossless coding and a mechanism to adjust the system parameters to control data rate and reconstruction quality. ii) We address the theoretical issues raised by Goyal et al. [14] in a practical way by using domain-specific knowledge: (i) we employ reconstruction algorithms that do not rely solely on sparsity but also on other properties of natural images, and (ii) we use deterministic sensing matrices that are known to be effective for such signals. iii) Having an end-to-end system enables us to benchmark its performance in terms of quality versus data rate. In order to give the reader a sense of CSbIC compression performance, we compare it with the popular JPEG and JPEG2000 standards. Our goal is not to propose CSbIC as a competitor to these and other standards such as HEIF or WebP. As we pointed out, in conventional applications for which these standards were designed for, the image compression needs are already well met. Rather, the purpose of the comparison to JPEG and JPEG2000 is to give the reader comparative levels of compression that can be expected. We show that our CSbIC system is onpar with JPEG, with a clear advantage in the low data rate range. Please refer to Fig. 1 as an example. iv) We describe various design choices in each of the system components and study the effects that these choices have on the overall system performance.
The system that we describe is far from being fully optimized, and throughout the paper we point out where further improvements may be made. Nevertheless, even at this preliminary stage the performance of our system makes it a viable alternative to the conventional methods (refer to Figs 4-6 for comparison with JPEG). The rest of this paper is organized as follows. Sec. II describes the system architecture and discusses the design choices in each component. Sec. III provides the general framework of reconstruction algorithms. Sec. IV presents the results of our performance testing and Sec. V discusses the implications of our results.

II. SYSTEM ARCHITECTURE
A diagram of the system architecture is given in Fig. 2. Each of the encoding steps is matched by a corresponding decoding step (in reverse order), with the exception of the bit rate/quality control block, which appears only in the encoder. In the following we present a detailed description of each of those processing steps. The upper part of Fig. 2 shows the corresponding block diagram of the JPEG encoder and decoder.
• In JPEG the first stage of the encoder (2D-DCT transform on 8 × 8 blocks) is lossless -the output is of the same dimension, N , as the input. The same is true for JPEG2000, although the transform is different (it is a wavelet transform on the whole image or tile). Therefore the quantizer and lossless coder have to process an Ndimensional signal, and their computational workload is about the same regardless of the compression level. • In CSbIC, the first stage (projection of whole image onto low dimension space) is lossy -the output of the first stage is of dimension M N , and hence the computational workload of the quantizer and lossless coder is much lower, proportional to M . Furthermore, if we do not use a fast transform in the first stage, it becomes a simple multiplication of a vector with dimension N by a matrix of dimension M × N , with O(M N ) complexity. This multiplication can be performed in the analog domain, using a CS camera, or it can be implemented by efficient dedicated hardware. In either case, the complexity of the first stage is also proportional to M . As we will see in Sec. II-C, when the joint bit-rate/quality control block decreases the quality, it does so in part by decreasing M . Therefore, in addition to the usual trade-off between quality and bit rate, CSbIC offers a trade-off between quality and encoder computational complexity.
As we will see later (Sec. IV-C), the choice of reconstruction algorithm offers a similar trade-off in the decoder. Thus, CS-bIC allows dynamic adjustment of quality based on available computational resources.

A. Sensing Matrix and Measurements Generation
We consider monochromatic images of N v ×N h pixels. The pixels of the input image are organized column by column as a pixel vector x ∈ R N (N = N v N h ). The pixel vector x is multiplied by a sensing matrix Φ ∈ R M ×N , M N , yielding the measurements vector Φ is a large matrix, even for small images. Therefore, for practical reasons, Φ is never stored, and the operation (1)   decorrelating an image, and most of the energy of the image is concentrated in the low frequency coefficients. Recently, new sensing matrices were introduced that leverage this property [36], [37]. These matrices are not incoherent with the common sparsity bases of natural images, as classical CS theory would require for guaranteeing robust reconstruction [2], [4], but generally they perform better than the classical sensing matrices in our application of image compression. For example, (1) can be implemented by performing a 2D-DCT on the image pixels and then reordering the resulting coefficients in a "zig-zag" order (similar to the one used in JPEG encoding), and then selecting the first M low-frequency coefficients [36].
In some applications, such as the single-pixel camera and lensless camera [32]- [34], [38], [39], a binary-valued matrix, i.e., a matrix whose entries are only ±1 (or {0, 1}), is more suitable for hardware implementation. In this case we approximate the 2D frequency decomposition by using a 2D Walsh-Hadamard transform (2D-WHT) [40]. Let where ⊗ denotes the Kronecker product [41] and W h , W v are Walsh-Hadamard matrices in sequence order [42]. Similar to the 2D-DCT case, the selected measurements are the first M coefficients of the zig-zag ordered entries of Wx, which can also be computed numerically in an efficient way. We can get the CS theoretical guarantee for successful reconstruction, w.h.p. (with high probability), by replacing the deterministic matrices described above with random matrices whose entries are independent, identically distributed (IID) random variables (RVs) with Gaussian or Rademacher distributions [43], [44]. These matrices are universal, i.e., w.h.p. they are incoherent with any given sparsity basis. Furthermore, the measurements generated by those random matrices are mutually independent and asymptotically normally distributed, which is helpful in the quantization and coding design. Such fully random matrices do not allow fast transform implementation of (1). But similar desired properties and performance guarantees were shown for structurally random matrices (SRM) [45]- [47] where Φx is obtained by applying a randomly selected permutation to x, computing a fast transform, on the permuted vector, and randomly selecting M of the transform coefficients as measurements (the DC coefficient is always selected). We denote these matrices SRM-DCT and SRM-WHT, respectively.

B. Quantization
The quantizer maps the measurements vector y into a finite sequence of codewords taken from a finite codebook C and the dequantizer maps the codewords sequence into a measurements vector that is an approximation of the original measurements vector. If the sensing matrix is deterministic, the measurements are highly uncorrelated; if it is a SRM the measurements are nearly independent. Hence the advantage of vector quantization [48] over scalar quantization is small and does not justify its added complexity [19]. Therefore, we consider a scalar quantizer that maps each measurement where Q i : R → C is the quantizer of the i th measurement. In this work we use the same quantizer for all measurements, and hence in the following we omit the subscript i from Q i .
The simplest scalar quantizer is the uniform quantizer. We select the "mid-tread" type, defined by where y denotes the largest integer not exceeding y, µ = 1 M M i=1 y i is the mean of the measurements, s is the quantizer's step,Q(y) is the unclipped quantized value, and L is a positive integer that determines the range s(L−0.5) of the actual quantizer Q(y). Consequently there are 2L+1 codewords, C = {−L, · · · , L}. Since the distribution of the measurements is highly non-uniform, the codewords distribution is also not uniform; hence in order to represent codewords effectively by a bit sequence we need to use variable length coding (VLC) in the lossless coder. A uniform quantizer with fixed length coding (FLC) with little data rate penalty is used due to its simplicity. Another reason to use a uniform quantizer is that the reconstruction may be sensitive to the presence of even a few measurements with large errors [16], which is often the case with non-uniform quantizers.
If Q(y) = c and |c| < L, we define the dequantizer by hence the quantization error is bounded by If |c| = L, the quantized measurement is saturated and the quantization error cannot be bounded. Even a small number of saturated measurements can cause severe quality degradation unless they are specially handled in the reconstruction [16].
The simplest way to do it is by not using the saturated measurements at all; attempts to modify the reconstruction algorithm to use those measurements showed little gain over simply discarding them. Another option (not considered in [16]) is to code the value ofQ(y) for each saturated measurement y in some ad hoc manner and transmit it as additional information. In both cases saturated measurements incur a penalty, either in the form of transmitted codewords that are not used, or as ad hoc transmission ofQ(y). Therefore, we select L large enough to make saturation [16] a rare event. In fact, L can be set sufficiently large to eliminate saturation completely, but a very large codebook may have an adverse effect on lossless coding (see Sec. II-D3). We found that a good trade-off is to select L so that quantizer's range s(L−0.5) is about 4 standard deviations of the measurements. However, the system is not very sensitive to this parameterperformance does not change much if the range is 3 or 6 standard deviations. With all the sensing matrices considered, the first measurement y 1 is the DC coefficient, which is the sum of all pixels in the image; y 1 is much larger than the other measurements and is always saturated, therefore it requires special handling: y 1 is excluded when calculating the mean (µ), and the standard deviation of the measurements, and Q(y 1 ) = L is not included in the quantized measurement. Instead,Q(y 1 ) is coded in an ad hoc fashion and transmitted separately. Unless the quantization is very coarse, its effect on the measurements can be modeled as adding white noise, uniformly distributed in [−s/2, s/2], which is uncorrelated with the measurements [48]. Hence the variance of the quantization noise in each measurement is The integral pixel values are generally obtained by sampling an analog signal and rounding the samples values to the nearest integer. Hence the pixel values contain digitization noise with variance of 1/12. This noise appears in the measurement y j with variance where φ j is the j th row of Φ. In the sensing matrices that we consider φ j 2 is constant, φ j 2 = Φ 2 . Clearly, there is no point in spending bits to accurately represent the digitization noise, hence we need to have σ Q ≥ σ D and consequently s ≥ Φ 2 .

C. Image Quality Control
As mentioned earlier, the performance of a compression system is measured in terms of reconstruction quality vs. bit rate. Usually, both of these quantities are determined by a single parameter. By varying this parameter we obtain the receiver operator curve (ROC), which depicts reconstruction quality against bit rate. However, in our case, both the coded image size and the quality of the reconstructed image depend on two parameters, the compression ratio R and the quantizer step size s. One can get the same coded image size with various combinations of these two parameters, but the reconstruction quality varies significantly among those parameter combinations. If the encoded image size is constrained not to exceed b, then the highest possible reconstruction quality Z * (b) is given by where B(R, s) and Z(R, s) are the coded image size and the reconstructed image quality (according to some specified quality criterion) respectively, as functions of the compression ratio and the step size. For each value of b, there is an optimal pair of compression ratio and quantization step, (R * (b), s * (b)), which is the solution of the optimization problem (9). However, getting an analytical solution to this equation, even approximately, turns out to be an extremely difficult problem 3 . It is relatively simple to obtain an approximate expression for B(R, s) 4 , but unfortunately, there is no closed form expression for Z(R, s). The CS literature provides several upper bounds for the mean square reconstruction error [52], but using any of these upper bounds as a substitute to the actual mean square 3 The joint compressive sensing (using sensing matrix to obtain measurements) and data compression (to quantize the measurements) is an unsolved problem in theory to our best knowledge. The most related work can be found in [49]; unfortunately, no closed-form solution has been shown. Most recently, Goldsmith el al. [50], [51] performed some research on the distortion limits of analog-to-digital compression, which can inspire some analysis for this joint analysis of compressive sensing and data compression. However, significant gaps still exist to derive any explicit formula. 4 Since we are using entropy coding, we may approximate the bit rate by the entropy of the measurements, where Hs(y j ) is the entropy of the measurement y j quantized with a step size s, H(y j ) is the differential entropy of y j , and the second inequality is based on the well-known approximation (for small step size), H(y j ) = H S (y j ) − log 2 s. H(y j ) can be computed by estimating the distribution y j over a large number of images, or in the case of stochastic sensing matrices, since the measurements are known to be approximately Gaussian and identically distributed, one could get an analytic expression for H(y j ).
error is a gross oversimplification and the formulas derived for (R * (b), s * (b)) from this "approximation" have no resemblance to the empirical results. Furthermore, the mean square error does not reflect image quality as perceived by the human visual system (this issue is elaborated in Sec. IV-B); hence even if we had an expression for Z(R, s) in the mean square error sense, it would not help us in solving (9) when the objective is to maximize the perceived quality.
Since an analytic solution is not available, we set out to search for an empirical one. For eight widely used images (shown in Fig. 4), we plotted curves of Z(R, s) against B(R, s), where R was kept fixed and s was the independent variable. Plotting these curves in the same graph, for different values of R (see the example in Fig. 11), made it clear that each value of R is optimal only for a particular quantization step, s * (R), that is (R, s * (R)) are the solution to (9) for a particular value of b. For s > s * (R), we could find R < R such that B(R , s) = B(R, s * (R)) but Z(R , s) > Z(R, s * (R)), and similarly, for s < s * (R), we could find R > R such that B(R , s) = B(R, s * (R)) but Z(R , s) > Z(R, s * (R)). Therefore, for each image the function s * (R) was monotonically decreasing, but it was different for different images. A closer examination of the data revealed that in each image R and s * (R) were approximately inversely proportional to each other, with a proportion constant of the same order of magnitude as Φ 2 . This could be expressed by: where C is a constant whose value depends on the type of sensing matrix and varies from image to image, but generally is between 1 to 10. Furthermore, we found that using C = 2.0 is a good choice, albeit suboptimal, for all images. To confirm these findings we ran our tests on 200 images from the Berkeley Segmentation Dataset and Benchmark dataset, i.e., BSDS300 [53] and verified that setting (R, s) using (11) with C = 2.0 is a generally good choice for all natural images. Thus in our tests s * (R) is determined using (11), with C = 2.0. The resulting quantization step is sufficiently fine to enable modeling the quantization noise as uncorrelated, uniformly distributed white noise. It should be emphasized that our choice of C = 2.0 is not optimal for many images, and Eq. (11) itself is an empirical approximation of a very complex relationship between the optimal compression ratio and optimal step size. Consequently, one could often get better results than those presented in this paper by manually adjusting s * (R) for a given image. Therefore, finding a better practical solution for the optimization problem (9) could lead to significant performance improvement in our system.

D. Lossless Coding
The lossless encoder encodes the codeword sequence generated by the quantizer as well as some miscellaneous information (e.g., µ, s, and the ad hoc representation of saturated measurements) as a bit sequence. The lossless decoder exactly decodes the codeword sequence and the miscellaneous information from the bit sequence.
1) Coded Numbers Format : Various types of numbers are coded by the lossless encoder (and decoded by the lossless decoder). Each type is encoded in a different way: Unbounded signed or unsigned integers are integers whose maximal magnitude is not known in advance. They are represented by byte sequences, where the most significant bit (MSB) in each byte is a continuation bit. It is clear in the last byte of the sequence and set in all other bytes. The rest of the bits are payload bits which represent the integer. The number of bytes is the minimal number that has enough payload bits to fully represent the unsigned or unsigned integer.
Real numbers, which are natively stored in single or double precision floating point format [54] are coded as pairs of unbounded signed integers representing the mantissas and the exponents in the floating point format.
Bit arrays are zero padded to a length which is a multiple of 8 and coded as a sequence of bytes, 8 bits per bytes.
Bounded unsigned integer arrays are arrays of unsigned integers, each of which may be represented by a fixed number of bits. An array of n integers, each of which can be represented by b bits, is encoded as a bit array of bn bits.
Each of these number formats can be easily decoded. Note that these formats are byte aligned for simpler implementation.
2) Entropy Coding : The codewords ±L represent saturated measurements. We merge these two labels into a single label L, thus we have 2L codewords: −L + 1, . . . , L.
Let p c , c ∈ C be the probability of a measurement to be quantized to c. If the codewords {Q(y i )} M i=1 are IID random variables, then a tight lower bound on the expected number of bits required to transmit these codewords is the entropy rate: Arithmetic coding (AC) [54] represents the codeword sequence by a bit sequence, the length of which can get arbitrarily close to the entropy rate for a large M . We use AC to encode the codewords sequence.
Since the probabilities p c , c ∈ C are not known a priori, they need to be estimated and sent to the receiver, in addition to the AC bit sequence. These probability estimates can be obtained in two ways: For SRMs, the measurements are approximately normally distributed [45], hence for |c| < L, p c is the normal probability of the quantization intervals [cs + µ − 0.5s, cs + µ + 0.5s), and p L = 1 − |c|<L p c . Thus, all that needs to be sent to the receiver is the estimated standard deviation of the measurements, which is coded as a real number. For deterministic sensing matrices, it is necessary to compute a histogram of the quantized measurements sequence, use it to determine the probabilities, and then code the histogram and send it to the receiver along with the AC bit sequence. Sending the histogram is an overhead, but it is small in comparison to the gain achieved by arithmetic coding.
In natural images the magnitudes of the coefficients of 2D-DCT or 2D-WHT decay quickly, hence measurements generated using a deterministic sensing matrix are not identically distributed, which violates the assumptions under which AC is asymptotically optimal. In order to handle this problem we partition the codeword sequence into sections, and for each section we compute a histogram and an AC sequence separately. This, of course, makes the overhead of coding the histograms significant. In the following we describe how the histograms are coded efficiently and how to select a locally optimal partition of the codewords sequence.
3) Coding of Histograms : In order to be efficient, the code of histograms of short codeword sequences should be short as well. Fortunately, such histograms often have many zero counts, which can be used for efficient coding. A histogram is coded in one of three ways: Full histogram: A sequence of 2L unbounded unsigned integers, containing the counts for each codeword. This method is effective when most counts are non-zero.
Flagged histogram: A bit array of 2L bits indicates for each codeword whether the corresponding count is non-zero, and a sequence of unbounded unsigned integers contains the non-zero counts. This method is effective when a significant share of the counts is zero.
Indexed histogram: A bounded integer indicates the number of non-zero counts, an array of bounded integers contains the indices of the non-zero counts, and a sequence of unbounded unsigned integers contains the non-zero counts. This method is effective when most of the counts are zero. In the extreme case of a single non-zero count the AC bit sequence is of zero length, hence this histogram coding is effectively a run length encoding (RLE).
The histogram is coded in these three ways and the shortest code is transmitted. A 2-bit histogram format selector (HFS) indicates which representation was chosen. The HFSs of all sections of the codeword sequence are coded as a bit array. Thus, each section of the codeword sequence is represented by the HFS, the selected histogram representation and the AC bit sequence.

4) Partitioning into AC Sections:
The partition of the codeword sequence into sections is preformed using a greedy algorithm (Algorithm 1). We begin by partitioning the sequence into RLE sections, and in each iteration we merge up to m consecutive sections, so that the merging yields the greatest possible reduction in total coded sections length; m is a constant which was set to 4 in our experiments. When the algorithm starts, all the bits are spent on HFSs and histogram representation, and none on AC bit sequences. As the algorithm progresses and sections are merged, more bits are spent on AC, and the histograms become fewer in number, but having more non-zero counts.
Determining the effect of merging several consecutive AC sections requires repeated calculation of the the coded length of AC sections. In principle this could be done by performing arithmetic coding and checking the resulting length. However, doing that in each iteration for each potential sections merge would be prohibitively computationally intensive. Instead we use the tight lower bound of (12) as an approximation to the section length. Let {h s (c), c ∈ C} be the histogram of section s, where h s (c) is the count for codeword c, and let be the number of codewords in the merged section. The codeword probability is approximated by p c ≈ h s (c)/M s , c ∈ C. Hence, using (12) and rounding up the result to the nearest multiple of 8 (for byte alignment), the estimated AC coded length of section s iŝ where a is the least integer not smaller than a. The total section code length, including the histogram is estimated bŷ where L H (s) is the number of bits used for the histogram coding and the last term on the right hand side is the bits used for the HFS. 5) Complexity of Algorithm 1: To analyze the amount of computations, we note that J r − m ≤ J r+1 < J r for 0 ≤ r < r * and therefore J r ≤ M − r and r * ≤ M − J r * ≤ M − 1. We also note that Consider the steps of the algorithm: Assuming that the arrays are organized in a tree structure, which allows insertion and deletion of elements at a fixed complexity, each iteration of step 2b and 2c is of complexity O(1) and the complexity of step 2d is O(m(m + 1)/2). Since the iterations are repeated at most (M − 1) times, the expression gives the complexity of Algorithm 1 as Since M N and the complexity of computing the measurements is O (N log N ), assuming a sensing matrix based on a fast transform, the complexity of performing the entropy coding is significantly lower than that of computing the measurements.
A note about the significance of choice of m: This is a greedy algorithm which stops when it cannot make an improvement by merging up to m adjacent sections. Increasing m will make the algorithm consider more options in each iteration, and therefore potentially find a better solution. Thus it appears that choosing m is a trade-off between complexity and better compression. However, if the algorithm chooses to merge more than two sections at a time it skips an iteration, therefore choosing a larger m may sometimes improve performance. In our simulations we used m = 4. i. Let S 0 j,k be the section obtained by merging sections S 0 j , . . . , S 0 k and compute the histogram resulting from this merge: , c ∈ C and compute the values ofL(S 0 j,k ) using the histogram and Eq. (16). ii. Using the above histogram, compute the gain in bit rate obtained by merging sections S 0 j , . . . , S 0 k : Let m be a small positive integer. for r = 0, 1, . . . do a. Let (j * , k * ) = arg max (j,k) g r j,k . b. If g r j * ,k * ≤ 0 break out of the loop. c. Merge sections S r j * , . . . , S r k * to get partition (r + 1), in which Let P r+1 = {(j, k)|1 ≤ j < k ≤ min(J r+1 , j + m − 1). For each pair (j, k) ∈ P r+1 : let S r+1 j,k be the section obtained by merging sections S r+1 j , . . . , S r+1 k and compute the histogram resulting from this merge: h r+1 j,k = k i=j h r i (c), c ∈ C, and compute the values ofL(S r+1 j,k ) and the the merging gain g r+1 j,k = k i=j l r+1 j −L(S r+1 j,k ) using the histogram and Eq. (16). In most cases the values will be copied directly from corresponding values of h r j,k (c), c ∈ C,L(S r j,k ) and g r j,k . Actual calculation is needed only if j < j * ≤ k < k * , that is for m(m + 1)/2 entries. end for 3. Let r * be the final value of r, J * = J r * be the final number of sections and S * j = S r * j , h * j (c) = h r * j (c), j = 1, . . . , r * , c ∈ C be the final sections and corresponding histograms. 6) Coded Image Structure: The data structure of the coded image is shown in Fig. 3. The coded image begins with a header specifying the image parameters (size, bits per pixel, color scheme, etc.) and the encoder's algorithmic choices (e.g., sensing matrix type, compression ratio, quantizer parameters). The body of the code consists of some global image parameters (µ, s, L,Q(y 0 ), etc.) followed by the entropy coded quantized measurements, which are the bulk of the transmitted data: the HFSs, the histograms and the AC sequences for each section. If the values ofQ(y) are transmitted for saturated measurements (corresponding to terms equaling L in the codeword sequence), these are coded as an array of unbounded signed integers.

E. Decoding
Decoding is done in reverse order of encoding, as follows: • The numerical parameters and bit arrays in the coded image are parsed. • The quantization codewords are recovered by arithmetic decoding. • The unsaturated quantized measurements are computed using (5). For the saturated measurements (those having a codeword of L), if values ofQ(y) are transmitted, they are used. Otherwise, the values of the quantized saturated measurements are set to zero. • The sensing matrix is determined according to the algorithmic choices in the coded image. If there was no ad hoc transmission ofQ(y) for saturated measurements, the rows corresponding to these measurements are set to zero. In practice this is done by replacing the original sensing matrix Φ by DΦ, where D is a M ×M diagonal matrix whose diagonal elements are zero for saturated measurements and one for unsaturated measurements. • The image is reconstructed using the sensing matrix and the quantized measurements, as described in detail in Sec. III.

F. Module Comparison of CSbIC with JPEG and JPEG2000
We now compare the architectures of CSbIC, JPEG and JPEG2000 and consider the aspects that may lead to performance differences. Fig. 2 compares the architecture of CSbIC with that of JPEG side by side. The main common points and differences are: • JPEG, JPEG2000 and CSbIC begin with a linear projection of the image onto a different space. However: -JPEG2000 may partition the image into tiles of varying sizes, which are processed separately. This is equivalent to using a block-diagonal projection matrix, where each block corresponds to a tile. In JPEG the tiles (referred to as blocks) are of fixed size of 8 × 8 pixels. In CSbIC the projection is done on the whole image, which was adequate for the image sizes that we experimented with. -In both JPEG and JPEG2000, each tile/block is projected on a space of the same dimension, N , hence there is no data loss in this operation, whereas in CSbIC the projection is lossy since it is on a Mdimensional space, M N . -In JPEG the projection is a 2D-DCT with the output organized in zig-zag order, while in JPEG2000 it is a 2D-wavelet transform. In CSbIC, a 2D-DCT based projection is one of several options. -JPEG uses block to block prediction of the DC coefficient in order to deal with the issue of the DC coefficient being much larger than the other coefficients. In JPEG2000 this is done by subtracting a fixed value from all pixels before the projection. In contrast, CSbIC takes care of this issue in the quantization stage. The effect of these different methods is similar and has little impact on performance.
• CSbIC uses a simple uniform scalar quantizer with an identical step size for all measurements. JPEG uses the same type of quantizer, but the step size is selected from a quantization table and is different for each coefficient, resulting in quantization noise shaping, which may give JPEG an advantage at higher quality/data rate (Fig. 4). JPEG2000 also performs quantization noise shaping through varying step size, and in addition, its quantizer is not exactly uniform-the quantization interval around zero (the "dead-zone") is larger than the other quantization intervals, effectively forcing small wavelet coefficients to zero and reducing the amount of bits spent on coding them. • The bit rate and quality trade-off in JPEG and JPEG2000 is controlled by tuning the operation of a single module -the quantizer. In contrast, in CSbIC this trade-off is controlled by jointly tuning two different modules: The projection, or measurement capturing module is tuned by changing the compression ratio, and the quantizer is tuned by changing the quantization step. • In JPEG, entropy coding is based on Huffman coding and RLE. JPEG2000 uses arithmetic coding, with a sophisticated adaptive algorithm to determine the associated probabilities. CSbIC uses arithmetic coding, which is known to be better than Huffman coding, but instead of using adaptive estimation of the probabilities, the codewords are partitioned into sections and for each section a histograms of codewords is computed and sent as side information. The overhead of the transmitted histograms may be a disadvantage of CSbIC relative to JPEG2000. • In JPEG and JPEG2000, the decoder generates the image from the dequantized transform coefficients by an inverse 2D-DCT or wavelet transform, respectivelya simple linear operation that does not rely on any prior information not included in the received data. In contrast, the CS reconstruction in CSbIC is an iterative, non-linear optimization, which relies on prior assumptions about the structure of the image (e.g., sparsity) or some recently developed deep learning based algorithms.

III. IMAGE RECONSTRUCTION
Since the seminal works of Candès et al. [2] and Donoho [3], various reconstruction algorithms have been developed. The early reconstruction algorithms leveraged the property of natural images of being compressible when projected by a suitable sparsity operator D: where f denotes the projected vector and it is usually forced to be sparse. D can be a pre-defined basis (DCT or wavelet), or learned on the fly [55]. Another popular sparsity operator is the Total Variation (TV). Recently, better results were obtained by algorithms such as D-AMP (denoising based approximate message passing) [6] and NLR-CS (nonlocal low-rank compressive sensing) [5], which exploit established image denoising methods or the natural images property of having a low rank on small patches. While using different projection operators D, most of these algorithms computex, the estimated signal, by solving the minimization problem where p can be 0 ( f 0 denotes the number of non-zero components in f ) or 1, using f 1 as a computationallytractable approximation to f 0 . Alternatively, p can stand for the nuclear norm * to impose a low rank assumption [5], [56]- [58].
Problem (20) is usually solved iteratively, where each iteration consists of two steps. Beginning with an initial guess, the first step in each iteration projects the current guess onto the subspace of images, which are consistent with the measurements, and the second step denoises the results obtained in the first step [59]. The denoising based image compressive sensing reconstruction algorithm (D-AMP), proposed in [6], integrates  various denoising algorithms [60] into the approximate message passing (AMP) [61] framework and leads to the state-ofthe-art reconstruction results using Gaussian sensing matrices. AMP extends iterative soft-thresholding by adding an extra term to the residual known as the Onsager correction term. By exploiting the nonlocal self similarity of image patches and low-rank property, Dong et al. [5] proposed the nonlocal lowrank regularization based compressive sensing reconstruction algorithm, i.e., NLR-CS. NLR-CS is an iterative optimization algorithm based on the alternating direction method of multipliers (ADMM) [62] framework. In each iteration, after the updating of x, instead of using various denoising algorithms as in D-AMP, NLR-CS imposes low rank property of similar patches to achieve a better image. 5 In our experiments we have used an improved variant of TV known as Generalized Alternating Projection Total Variation (GAP-TV) [65], as well as D-AMP and NLR-CS.

IV. PERFORMANCE TESTING
Our initial test material consisted of 8 widely used monochrome images of 256x256 8-bit pixels (Fig. 4). These images were processed in a variety of test conditions, specified by encoder and reconstruction parameters. The outcome of processing an image in a particular test condition is the data rate, as expressed by the size of the coded image file (in bytes), and the reconstructed image quality, measured by structural similarity (SSIM) [35]. We preferred SSIM over Peak-Signalto-Noise-Ratio (PSNR) as a measure of quality because SSIM better matches the quality perception of human visual system (refer to Fig. 5 and [1] for PSNR results, which provide similar observations to SSIM). Fig. 6 shows a few examples where higher SSIM is clearly consistent with higher image quality. Results are presented as points in SSIM vs. coded image size diagram. By connecting the points corresponding to test cases with different rate control parameters (but identical in all other parameters) we get performance curves which can be compared with each other to determine the better operating parameters. Performance curves for JPEG or JPEG2000 were obtained using the MATLAB implementation of those standards (the "imwrite" function). We confirmed our results by repeating our tests on 200 images from the BSDS 300 dataset [53]. Fig. 7 shows SSIM average and standard deviation, over the 200 images, vs. compressed file size. The 2D-DCT sensing matrix with NLR-CS reconstruction performs better than JPEG for compressed file size less than 4000 bytes. It is also more robust to variations among images, as shown by its significantly lower standard deviation. All the results presented in this section are supported by performance diagrams of additional images in [1], along with the source code.

A. Effect of The Choice of Sensing Matrix
The performance with the deterministic sensing matrices, 2D-DCT and 2D-WHT, was always significantly better than with the SRMs with the same fast transforms, SRM-DCT and SRM-WHT, respectively (Fig. 8), regardless of the reconstruction algorithm that was used. Within each of those groups (deterministic matrices and SRMs), the sensing matrices based on DCT generally yielded better performance than the those based on WHT (Fig. 8 right). However, (i) the difference, in terms of SSIM for the same file size, is generally smaller than the performance difference between deterministic and structurally random matrices, (ii) the difference between SRM-DCT and SRM-WHT is generally much smaller than the difference between 2D-DCT and 2D-WHT, and (iii) the magnitude of the difference varies with the image and the reconstruction algorithm that is used. In particular, with GAP-TV the performance difference between DCT and WHT based matrices was most pronounced, whereas with D-AMP WHT performance was very close to that of DCT, and in a few cases the WHT based matrices slightly outperformed the DCT based ones.
An obvious reason for these differences is the fraction of the signal energy that is captured by the measurements with each sensing matrix. The DC measurement is the same, up to a scaling factor, in all four matrices. Therefore, when speaking about signal energy and measurements we refer only to the AC component and exclude the DC component. Since SRMs randomize and whiten the signal prior to applying the transform, their measurements capture about M/N of the signal energy. The 2D-DCT measurements capture a much larger fraction of the signal energy because most of the signal energy is concentrated in the low frequencies. Since 2D-WHT is a crude approximation of spectral decomposition, its low "frequency" components capture less energy than those of 2D-DCT, but much more than the M/N of the SRM matrices. This argument is certainly true, but it is only a part of the explanation. Consider the sensing matrix given bỹ where Φ is a M × N 2D-DCT or 2D-WHT deterministic sensing matrix and Θ is given by whereΘ is a random orthonormal where R is a diagonal matrix whose diagonal entries are IID random variables that have the values of ±1 with equal probability; H is the DCT transform matrix of order M − 1, with rows scaled to have unit norm; and S is a random selection matrix, i.e., a matrix in which each row and each column contains one entry of 1 and the rest are zero. The measurements computed byΦ are obtained by a random orthonormal transformation (ROT) of the measurements computed by Φ, hence we denote the sensing matrices Φ ROT-2D-DCT and ROT-2D-WHT, respectively. The ROT matrices capture the same components of the signal as the corresponding deterministic matrices and the addition of the orthonormal transformation Θ should not matter to any of the reconstruction algorithms. Furthermore, the quantization step, given by Eq. (11), is the same and so are the quantization noise and digitization noise, given by Eq. (8). Nevertheless, the performance of the ROT matrices, while significantly better than the SRM matrices, was not as good as the deterministic ones. The difference in SSIM between the performance curves of a deterministic matrix and the corresponding ROT matrix was about a third or a quarter of the difference between the performance curves of the same deterministic matrix and the corresponding SRM matrix (see example in Fig. 9). Thus, multiplying the measurements vector Φx by Θ, causes a significant performance degradation. By its definition,Θ is a DCT-based SRM with a compression ratio of 1, which uses local randomization, that is, the input signal is randomized by randomly toggling the signs of its entries. This randomization method is different from the global randomization used in the SRM-DCT sensing matrix, where the input signal is randomized by randomly permuting its entries, but both variants have a similar effect on the distribution of the measurements: when a vector is multiplied by either one of these matrices, the distribution of the entries of the result is approximately zero mean Gaussian or a mixture of zero mean Gaussians with similar covariances [45]. Figure 10 shows the effect of multiplication byΘ on the AC measurements. The measurements generated by the deterministic matrices have a narrowly peaked distribution with a long tail, while the measurements generated by the ROT matrices have a much wider, Gaussian like distribution, with a shorter tail. As a result, the entropy coding needs significantly more bits per measurement for coding the ROT measurements. The distribution of the measurements generated by SRM-DCT and SRM-WHT has also a Gaussian-like shape, similar to that of the ROT generated measurements. Therefore, we may conclude that the difference between the SRM and the corresponding deterministic matrices has two components: one component is the different fraction of signal energy captured by each method, and is represented by the improvement from the SRM matrices performance to that of the ROT matrices. A second component relates to the different distributions of the measurements in the deterministic case and it is represented by the improvement from the ROT matrices performance to that of the corresponding deterministic matrices.
One can notice in Figure 10 that the histogram of 2D-WHT is a little wider than that of 2D-DCT. Thus, the advantage of 2D-DCT over 2D-WHT seems also to be not only because of the amount of signal energy captured in each case but also because of differences in the measurements distribution. For all the three algorithms, SRM-DCT performs similar to SRM-WHT, since the captured fraction of signal energy, and the measurements distribution are similar in both cases.

B. Analysis of Joint Quality Control
The effect of quantization on performance can be seen in Fig. 11. In general, as the quantization step is decreased, the compressed file size and the reconstructed image quality are increased. However, the rate of increase is not constant. For each compression ratio there seems to be three distinct regions for the quantization step. In the lower quality/data rate region, decreasing the quantization step results in a significant gain in quality with little increase in the compressed file size; in the higher quality/data rate region the opposite is true: decreasing the quantization step increases the file size significantly with hardly any quality gain; and in the narrow intermediate region between those two we get moderate increase in quality and file size as the quantization step is decreased. The intermediate region occurs at higher quantization step for lower compression ratios. The highest possible reconstruction quality of Eq. (11) is the maximum of all the curves depicted in Fig. 11 (as well as any potential curve of intermediate compression ratio). It is clear from the figure that at each bit rate b, the maximum is achieved in a different curve and b is at the center of the intermediate region of that curve.
As mentioned earlier in Section II-C, after extensive analysis of the eight images shown in Fig. 4, we determined empirically that the maximum is achieved when the compression ratio and quantization step satisfy (11). While the optimal value for the constant C in the right hand side varies from picture to picture, we found that C = 2.0 gave good results for all images and we used it in our experiments on 200 images from the BSDS dataset [53].

C. Effect of the Reconstruction Algorithm
We found that the more modern reconstruction algorithms, NLR-CS and D-AMP, give a better performance than GAP-TV (Fig. 4). In most cases NLR-CS is a little better than D-AMP, especially at low-quality/low bit rate settings. We note however, that both NLR-CS and D-AMP require careful manual parameter tuning and data normalization in order to perform well, so it is possible that the results would be somewhat different with better tuning.
The performance gains of D-AMP, and NLR-CS come at the cost of a much higher complexity. The reconstruction times of GAP-TV, D-AMP and NLR-CS were 0.6, 35 and 147 seconds receptively, with the same compression rate on the same Intel i7 CPU with 24G memory. The algorithms were implemented in Matlab. We expect that an optimized implementation in a compiled programming language would reduce the computation time by one or two orders of magnitude, and a significant further improvement can be achieved by using a GPU. However, it appears that even with optimized implementation the last two algorithms would run considerably slower than GAP-TV, because they are inherently more complex. Therefore, the choice of reconstruction algorithm can present a trade-off between quality and execution speed. Using different reconstruction algorithms, an implementation can offer a "fast mode" for quick browsing and an "accurate mode" for scrutinizing blown up images.

D. Effect of Quantization
In order to get a better insight into the performance of our system, we isolate the component of the error which is due to quantization by running tests where the measurements vectors are reconstructed without quantization. Of course, in this case we cannot draw the usual performance curves of SSIM vs. compressed image file size, because without quantization (and entropy coding), there is no meaningful "compressed image file". Instead, we draw curves of quality, with and without quantization, vs. the compression ratio (CSr). In Fig. 12 we plot the SSIM of reconstructed images processed by CSbIC (including quantization), with that of images reconstructed directly from the unquantized measurements vectors. As expected, the results without quantization are consistently better than the results with quantization. In SSIM terms, the total decoding error is the difference between 1, which represents perfect reconstruction, and the SSIM value of the image decoded by the CSbIC. The error component due to quantization is the gap between curves with and without quantization. This gap varies with the reconstruction algorithms and compression ratio, both in absolute terms and relative to the total decoding error. GAP-TV is not that sensitive to quantization error, with a SSIM decrease of 0.02-0.07 due to quantization. NLR-CS is the most sensitive, with SSIM decrease of 0.04-0.10. The behavior of D-AMP is not uniform: at compression ratios below 0.10 the SSIM decrease is similar to NLR-CS, but at higher compression ratios the SSIM difference decreases, down to 0.01-0.04 at a compression ratio of 0.25.

E. Comparison with JPEG and JPEG 2000
CSbIC was tested with various sensing matrices and various reconstruction algorithms. As might be expected, the best combination depends on the particular image as well as on the quality or data rate operating point. However, in general the best performance was achieved with 2D-DCT+NLR-CS or 2D-WHT+D-AMP. In almost all test images CSbIC is decidedly better than JPEG at low bit rates (Fig. 4). A visual inspection shows (Fig. 6) that at low bit rates JPEG images have blocking artifacts, while in CSbIC the artifacts are more subtle. In half of the images CS compression is better than JPEG for all quality levels. In the rest of the images JPEG is better for the higher quality cases, with the crossover happening at SSIM of 0.7-0.9 (0.8-0.9 for two out of the four images).
At this stage, the performance CSbIC is inferior to that of JPEG2000, but at the low quality regions, below about 0.7 SSIM, their performance is quite similar (Fig. 4).

V. DISCUSSION
Both JPEG and CS leverage the inherent compressibility of natural images, but in different ways: In JPEG this is done at the encoder -the 2D DCT decomposition on 8 × 8 pixel blocks essentially represents the signal according to a sparsity basis; in classical CS the compressibility assumption comes into play only at the reconstruction. JPEG also exploits another property of natural images, namely that the high magnitude coefficients are the low frequency ones, therefore JPEG allocates more bits to their accurate representation. This domain-specific insight was built into CSbIC by using the deterministic sensing matrices described in Sec. II-A, which improved performance significantly in comparison to the SRMs of classical CS.
JPEG achieves lower bit rates solely through quantization, while CSbIC does it by quantization, and in addition, by discarding most of the transform coefficients during the measurements selection. This is possible because CS reconstruction recovers the missing coefficients by invoking various structural assumptions about the original image. This gives CSbIC a significant advantage at low bit rates, where excessive quantization results in poor performance of JPEG. Another possible explanation for that is that JPEG processes 8 × 8 pixel blocks nearly independently. At low bit rates, only a small number of DCT coefficients are adequately represented in each JPEG block, which may be too little for blocks with much detail and causing the blocking artifacts which were observed. In contrast, all CS measurements contribute equally to the description of the whole image.
CSbIC yields better performance with certain types of sensing matrices and reconstruction algorithms. However, the less successful candidates also performed decently, and they may be the preferred choice in some situations. Binary-valued sensing matrices may be preferred for capturing measurements in the analog domain, and faster reconstruction may be useful when speed is important, for example, for quick image browsing.
JPEG is the culmination of many years of research into the properties of natural images. Wherever we could we leveraged this body of knowledge in CSbIC. However, much improvement can still be made, both in the quantization and entropy coding aspects of the system and in its CS aspects. The latter includes the sensing matrix design, the reconstruction algorithm, and the derivation of theoretical guarantees for successful reconstruction. The following are a few examples of such potential improvement in the area of reconstruction.
Recently, deep convolutional neural networks have been used, with good results, to solve the inverse problem (i.e., the reconstruction) in CS [34], [66]- [69]. This approach could dramatically reduce the reconstruction computation time, and po-tentially lead to improved performance, especially in domainspecific applications where the specifics of the domain may be difficult to express analytically. Such networks need to be trained on a very large image database. The quality of the reconstruction depends on the size and variability of the content of the database, which makes it difficult to make a fair comparison between the neural nets methods and the analytical reconstruction methods discussed in this paper. However, this is one of the most promising directions for improving the system. As an additional point, the snapshot compressive video imaging systems [70], [71] along with the recently developed algorithms [67], [72]- [74] and the theoretical guarantees derived in [75] demonstrate a promising way for video coding based on CS. This will be one of our future directions and can also be used for hyperspectral image [76], [77] compression.