{'fraction_non_alphanumeric': 0.031618807194698645, 'fraction_numerical': 0.012811612496055539, 'mean_word_length': 4.793784277879342, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Utilizing the Wavelet Transform's Structure in Compressed Sensing

Compressed sensing has empowered quality image reconstruction with fewer data samples than previously though possible. These techniques rely on a sparsifying linear transformation. The Daubechies wavelet transform is a common sparsifying transformation used for this purpose. In this work, we take advantage of the structure of this wavelet transform and identify an affine transformation that increases the sparsity of the result. After inclusion of this affine transformation, we modify the resulting optimization problem to comply with the form of the Basis Pursuit Denoising problem. Finally, we show theoretically that this yields a lower bound on the error of the reconstruction and present results where solving this modified problem yields images of higher quality for the same sampling patterns.


Introduction
Reducing the number of data samples required to generate a quality image is often beneficial. Compressed sensing [1] has been a remarkable advancement to this end for several imaging systems including Magnetic Resonance Imaging (MRI) [2,3], computed tomography (CT) [4], radio astronomy (RA) [5,6], and optical imaging [7]. Many imaging applications (e.g. MRI, CT, RA) collect samples in the Fourier domain. Conventionally, the high quality of compressed sensing relied on a sparsifying transform such as the wavelet transform [8] or the gradient operator [9]. More recently, learning techniques that use a training set of data to (either explicitly or implicitly) determine a basis for signal reconstruction have been created [10,11,12]. However, learning methods can suffer from over-fitting to the training data, which may take the form of extreme sensitivity to small noise, hallucinations, or an inability to accurately represent the image [13,14]. Wavelets provide an orthogonal computationally efficient transform [15] that, when combined with Fourier sensing, satisfy asymptotic incoherence and multi-level sparsity required for accurate reconstruction [16]. Due to these advantages, the wavelet transform continues to be used in compressed sensing applications [17,18,19,20,21,22,23].
The vast majority of compressed sensing research has considered a general linear transformation as the sparsifying transform. The Discrete Daubechies wavelet transform (DDWT) is an effective choice for compressed sensing [24]. One might hope that when this specific transform is used, the reconstruction algorithm could take advantage of its properties to improve the quality of compressed sensing results. That is the subject of this work. We identify an affine transformation, based on the DDWT, that increases sparsity. And we present a method to adapt this transformation into the existing compressed sensing framework.

Theory
With compressed sensing of noisy imaging systems, one solves the following Basis Pursuit Denoising (BPD) problem: where A is the system matrix, b is the data vector, 2 ≥ 0 is a bound on the noise power [25], and · p is the L p norm. This problem can be solved efficiently with standard algorithms (e.g. the Fast Iterative Shrinkage Threshold Algorithm -FISTA [26,27]).
For many imaging systems, including those discussed in the introduction, A = M F W −1 , where M is a sampling mask, F is the discrete Fourier transform, W is the DDWT, and x is the vector of wavelet coefficients. Once (1) is solved, the image can be recovered with y = W T x. This system matrix A satisfies RIPL [28] which guarantees that x can be recovered robustly [29] and that it solves the corresponding sparse signal recovery problem [30,31].
The DDWT consists of two perfect reconstruction finite impulse response (FIR) filters: a low pass and a high pass filter. There are different types of DDWT, which correspond to different orders of the filters (different numbers of filter coefficients). The coefficients of the DDWT-4 and their spectrums are shown in Fig. 1. The process of applying a wavelet transform is to apply each filter and downsample by 2 (and then concatenate the results). We label the response of each filter/downsample as a bin. This process may be applied recursively to each resulting bin. In the images of Fig.  2, the recursion was applied four times only to the lowest frequency bin. After transforming, the lowest frequency bin is downsampled by 16. Since the composition of FIR low-pass filters is an FIR low-pass filter, the lowest frequency bin in the transform domain is a low-pass filtered (and downsampled) image. The DDWT low-pass filters are not ideal filters (meaning that they have a non-zero transition region and that their support extends beyond the cutoff frequency); therefore, the lowest frequency bin of the transform also includes aliasing. Since the other bins of the Wavelet transform of natural images are sparse, by the sifting property of convolution, the aliasing artifacts present in the lowest frequency bin must also be sparse.
To estimate the lowest frequency bin, since it is not typically sparse for natural images, there is no advantage to utilizing the theory of compressed sensing. Instead, we can rely on the Nyquist-Shannon sampling theorem [32] and collect evenly spaced samples at twice the cutoff frequency. This technique, fully sampling a region centered on the 0 frequency in the Fourier domain, was developed as a heuristic and has been used extensively in compressed sensing [33,34,35,36]. The size and shape of the fully sampled region, though, had not been theoretically justified in these techniques. For example, the SAKE method synthesizes a fully sampled square region of size 80×80 to reconstruct an image of size 200 × 200 [37]. The DISCO method uses a fully sampled spherical region [35,36]. By only collecting the number of samples required to satisfy the Nyquist-Shannon theorem, the total number of samples is reduced over these heuristic techniques.
Let r denote the number of recursion levels of the DDWT. Then the size of the lowest frequency bin, after downsampling, would be (M/2 r ) × (N/2 r ) pixels 2 , where the image is size M × N . Thus, its resolution would be According to the Nyquist-Shannon theorem, this image can be reconstructed accurately if a rectangular region of size (2 r /M ) × (2 r /N ) with evenly spaced samples centered on the 0 frequency is collected; we denote this subset of samples as the Fully Sampled Region (FSR).

Methods
Our approach will be to first reconstruct a low-frequency image and then correct this estimate with high frequency information. One could reconstruct a low frequency image simply by performing an Inverse DFT on the data collected in the FSR. However, doing so leads to ringing (Gibbs phenomenon), which increases the energy in the high frequency bins of the wavelet transform of the image. Instead, we apply a separable Kaiser-Bessel window [38,39] with a parameter of 4.
where F * is the adjoint of the unitary F , and M L and K B are diagonal matrices representing the mask of the fully sampled region and the Kaiser-Bessel kernel, respectively. Then we wish to estimate y by solving the following problem: minimize The key insight is that z = W (y −ŷ L ) will be more sparse (have more values approximately equal to 0) than x for the reasons presented in section 2. Figure 3 shows z for the images of Fig. 2. Indeed, they are now better approximated by a sparse representation.
With this definition of z, (2) is equivalent to minimize where β = b − M Fŷ L . After solving (3) for z , the image can be reconstructed with y = W −1 z +ŷ L . We call this the More Sparse Basis Pursuit Denoising (MSBPD) solution. Problem (3) is also a basis pursuit denoising problem; the difference between it and (1) is the data vector is replaced with β. Since A satisfies RIPL [40], the bound on the error provided in [41] applies. Since the vector z has more elements with value closer to 0, the error that results from solving (3) has a smaller bound than that of solving (1).

Results
All experiments in this study were performed with images of size 512 × 512. All optimization problems were solved using FISTA with line search [27,26] run for 100 iterations. The wavelet transform applied was the DDWT-4, recursively applied to only the lowest-frequency bin r = 4 times. Figure 4 shows the sampling patterns used in the experiments with various sampling percentages with and without the fully sampled center region. The sampling patterns are realizations of a random separable Laplacian distribution [42] with a standard deviation of approximately 20%. This pattern has a higher probability of sampling higher frequencies than a Gaussian distribution with the same standard deviation. This likely better matches the actual power spectral density of natural images and is, therefore, a more optimal sampling pattern for a compressed sensing reconstruction [43]. For the sampling patterns with the fully sampled region, fewer samples from the distribution were included in order to retain approximately the same number of total samples. The value of was chosen independently for each reconstruction by conducting an exhaustive search to find the value that minimized the relative error e = y − y 2 / y 2 . Figure 5 shows reconstructions using BPD and MSBPD. As can be observed in Fig. 5, MSBPD offers the highest quality reconstruction for all sampling percentages. The improvements become more noticeable as the sampling percentage is reduced. Table 1 shows the relative errors for the reconstructions in Fig. 5. In all cases, MSBPD yields the lowest relative error. The most significant improvement comes by including the FSR in the samples; a minor gain is attained by using MSBPD over BPD on this sampling pattern. The improvement is more pronounced as the sampling percentage decreases.  Table 2 shows the Structural Similarity metric (SSIM) comparing the reconstructed images to the originals. This quantification shows the same trends as the relative errors. In Fig. 6, we zoom into the reconstruction of the knee from 8% of data. This figure shows fewer wavelet artifacts in the result of MSBPD than exist in the result of BPD with either sampling pattern. Figure 7 shows the reconstructions with 8% of the data for BPD and MSBPD on all images of Fig. 2. Figure 8 shows the magnitudes of the differences between the reconstructions with 8% of the data shown in Fig. 7 and the original images. The errors are significantly reduced by including the FSR. There is an additional improvement to the reconstruction when using MSBPD. Table 3 shows the relative error for each algorithm on each of the five reconstructions of Fig. 7. As observed previously, the most significant gain is attained by including the FSR. A minor additional gain is attained with MSBPD.  Table 4 shows the Structural Similarity metric (SSIM) comparing the reconstructed images to the originals. Once again, this metric shows the same trends as the relative errors.

Discussion
In this work, we have utilized the structure of the Discrete Daubechies Wavelet Transform to improve image reconstruction results based on compressed sensing. It takes advantage of the prior knowledge that a great deal of images are not sparse in the lowest-frequency bin of the image's wavelet transform.  Note that some images, those that are natively sparse, are indeed sparse in the lowest frequency bin of the DDWT; e.g. images of angiography and some astronomy. There would be no advantage to using MSBPD over BPD for images of that type, and other types of compressed sensing reconstruction algorithms may be more appropriate [44]. When it is the case, however, that the image is not sparse in the lowest frequency bin then MSBPD offers reconstruction with improved quality over BPD even when the sampling pattern includes the FSR.
Due to the increased sparsity, an accurate reconstruction may be possible using the greedy Orthogonal Matching Pursuit [45,46] algorithm to solve the corresponding sparse signal recovery problem: where · 0 is the L 0 penalty and S limits the number of non-zero elements in the vector z. This could make image reconstruction more computationally efficient. The DDWT and the DFT are both radially asymmetric transforms; i.e., a rotation operator and the transform operator do not generally commute. However, in many imaging systems, there is not anything inherently special about the horizontal and vertical directions (e.g. MRI, CT, radio interferometry). Therefore, it may be possible to reconstruct images of higher quality by collecting a circle of data in the low-frequency region (rather than a square) and use a sparsifying transform that is radially symmetric. Symmetric wavelets on the sphere [47] may be such a sparsifying transform. Though this may yield improved quality, the computations required may increase. This side-effect of this possible improvement should be considered for any given application.
With each recursive application of the Wavelet transform, the average sparsity is reduced [16]. Thus, the value of r can be chosen based on the sparsity achievable for the imaging system and subjects of interest. Additionally, the higher the value of r, the more computations are required to implement the wavelet transform. Both aspects should be considered when determining a value of r for a given application.
In this work, we determined the regularization parameter with an exhaustive search. This is only appropriate when ground truth is known and not a practical solution. For some applications, the bound on the noise can be determined with an image of noise; this can be accomplished in MRI, for example, by imaging with an excitation of 0. When a calibration image is unavailable, reconstruction may be accomplished by adapting the MSBPD algorithm to an iterative re-weighting algorithm [48,49,50,51]. Alternatively, an automatic parameter tuning algorithm can be used based on existing training data [52] or with properties of the resulting pareto-optimal set [53]. Some of these methods have the consequence of implicitly altering the objective function [48], which means that the theorems of compressed sensing no longer hold. However, they have been shown heuristically to improve image quality in some cases.
Finally, the method presented in this work can be adapted to situations where information is shared amongst several acquisitions. This can happen with multiple-acquisition MRI [54], parallel MRI [55], or synthetic-aperture radar imaging [56] and may involve joint-sparsity regularization [57]. Moreover, the method can be adapted to deep-learning methods that are regularized with wavelet sparsity or joint sparsity [58].
We leave the investigations of these possible extensions as future work.

Conclusion
In this work, we have presented the MSBPD algorithm. This algorithm utilizes the structure and behavior of the DDWT to justify a sampling pattern and to identify a new sparsifying transform that increases sparsity of most natural images.
Since the system matrix satisfies RIPL, this leads to improved results over solving the standard BPD problem of (1). In experiments, we compared image quality of reconstructions made with BPD, BPD with the FSR, and MSBPD. In There is a significant reduction in the error by including the FSR in the sampling pattern. An additional reduction in error is realized by using MSBPD.
all cases, the most significant gain was attained by including the FSR in the sampling pattern and MSBPD yielded the lowest relative error.

Acknowledgments
ND would like to thank the Quantitative Biosciences Institute at UCSF and the American Heart Association as funding sources for this work. ND is supported by a Postdoctoral Fellowship of the American Heart Association. ND and PL have been supported by the National Institute of Health's Grant No. NIH R01 HL136965.