Morlet wavelet based nonergodic random sampling
Incoherent sampling is based on patterns dissimilar to image contents. On the other hand, a sparse wavelet representation of the image could be also found rapidly by probing the image directly with wavelet functions, if the most probable elements of the wavelet representation are known beforehand. What we propose here, is to combine these two contradictory lines of reasoning into a novel sampling scheme which is both random and based on a wavelet representation at the same time. We propose to apply a novel kind of sampling, equivalent to random sampling in the feature space. A feature space is built out of vectors, whose elements correspond to specific features of images. Simple features may be associated with spatial and frequency contents of an image. For instance, a feature space may be constructed using Gabor filters which are defined as Gaussian functions modulated with a linear phase dependence. A two-dimensional Gabor filterhas the following form ,
where is a normalization constant such that . A feature vector is constructed out of a set of dot-products of the image with Gabor filters where the parameters are related to probing a certain location of the image, parameter determines the characteristic scale of the feature, and select the part of the probed spatial spectrum. Gabor functions allow for probing images in the spatial domain and in the frequency domain at the same time with the highest possible resolution 
. In fact, the Fourier representation of the Gabor function is also a Gabor function, and both functions optimize the uncertainty relation for the two-dimensional Fourier transform. In other words, it is not possible to construct narrower probing functions (with smaller variances) in the spatial and frequency domains at the same time. A zero-mean and normalized Gabor filter is known as the Morlet wavelet or Gabor wavelet. In a two-dimensional situation, a Morlet wavelet is equal to
where the constants and assure that the wavelet function is normalized and has zero mean . Parameters are related to the size of the Gaussian envelope, number of periods within the envelope, and the orientation of modulation. A feature vector is obtained by convolving the wavelets with the image .
Taking a measurement with a single-pixel detector consists in probing the measured image with a set of sampling functions . The set of measured dot-products is later used to reconstruct the image . Now, let us probe the feature space with random functions . Since instead of probing the feature space, we propose to probe the image directly with modified sampling functions . These sampling functions are obtained as a convolution of Morlet wavelets with realizations of white Gaussian noise. Some examples of sampling functions and the procedure for their calculation are illustrated in Fig. 1a,b).
Many interesting natural phenomena in physics, biology or artificial intelligence arise on the verge of random and deterministic behavior of a system. In this work, the proposed sampling functionsare calculated by convolving random functions and deterministic Morlet wavelets , and they clearly combine random and deterministic properties. Mathematically, they are zero-mean random matrices with multivariate Gaussian probability density distributions. Each has a distinct power spectrum, dependent on the wavelet used for its generation. As opposed to the wavelets, the sampling functions have a similar random shape at any location, which reflects the property of stationarity. However, their ensemble properties are distinct from the properties of every single realization. This means they do not satisfy the statistical property of ergodicity. This makes them a lot different from the uncorrelated random sampling often used in compressive sensing applications, as well as from the deterministic sampling with noiselet, Walsh-Hadamard or cosine functions as well as from localized functions such as wavelets.
The choice of parameters for the random wavelet-based sampling functions has an obvious influence on the quality of a compressive measurement. By decomposing an image database with images of various content, we found that there exists a common parameter range that may be successfully used to represent most of the images with our sampling functions. This decomposition is not unique and finding an optimal decomposition is challenging from the computational viewpoint. Instead, we have used a simplified approach. We generated a large number of sampling functions with randomly selected parameters , placed them into a rectangular matrix, and decomposed every image into these sampling functions by left-multiplying the image by the pseudoinverse of this matrix. The Moore-Penrose pseudoinverse is a generalization of matrix inverse for rectangular and singular matrices, it finds application in image reconstruction from compressive measurements , and further we also use it as one of the methods for image reconstruction in this paper. In this way, we found a typical distribution of coefficients required to represent a large variety of real-world images. A graphical representation of the decomposition projected onto the parameter space is shown in Fig. 2. As we can see, the interesting part of the feature space spanned with is easily identified from this plot. A method with an adaptive choice of sampling functions could be further developed through a more in-depth analysis of how the decomposition varies for different images. However, in our case the range of obtained parameters was similar for the whole image database and the non-adaptive approach taken in this paper is certainly easier to implement, especially in experimental conditions.
The practical benefit of using the proposed sampling functions becomes clear from a simple comparison with classical sampling methods based on Walsh-Hadamard or noiselet functions. Figure 3 presents results of a simulated measurement from a single-pixel detector conducted at a low compression rate of with the use of three different sampling protocols, including the proposed Morlet wavelet-based random functions. The simulation was performed for two images with different properties, such as spatial frequency spectrum, contrast, or richness of details. The images have been reconstructed by minimizing the total variation norm (TV), which is one of the basic image reconstruction approaches used in CS. Figures 3b,f) show the reconstructions from the measurements with Walsh-Hadamard functions, which do not at all resemble the original images. The average quality of the image recovery, measured by the PSNR criterion (see Methods section), for the Walsh-Hadamard sampling is on the level of 14 dB. This is not surprising at such a low compression rate, although noiselet sampling gives considerably better results (PSNR dB). On the other hand, the Morlet wavelet-based sampling, also measured at the compression rate of only allowed us to reconstruct high quality images with the PSNR of 25 dB on average (see Fig. 3
). We would like to emphasize that these measurements are not adaptive, and the sampling functions are randomly selected from the previously estimated range of values of, as well as that images shown in Fig. 3a,e) were not included in the training database. Moreover, the quality and resolution of these reconstructions are uniform within the entire image areas and it is not possible to notice any characteristic resolution, orientation or region of interest of the images, which is enhanced at the cost of some other property. We think that this impressive result comes from the efficient sampling of the properly selected part of the feature space which maps the images accurately both in the spatial and spatial frequency domains.
Morlet wavelet-based random sampling functions allow to reconstruct images from smaller number of measurements. However, these patterns are neither binary nor orthogonal. We will now discuss the practical consequences of these important limitations and show how to overcome them. In optical single-pixel detectors, light modulation is usually accomplished by using spatial light modulators, such as digital micromirror devices (DMD). The DMDs are capable of displaying binary patterns at a rate of over kHz. Gray-scale modulation could be achieved with time multiplexing or would require using a different spatial modulator, for instance based on liquid crystals. However, both of these approaches offer much lower effective modulation rate. Therefore, we have decided to binarize the real-valued Morlet wavelet-based random sampling functions to retain the high measurement performance with a binary DMD. Some examples of binarized sampling functions are shown in Fig. 1
c). Although binary sampling functions are advantageous for displaying on DMD modulators, they however make the reconstruction of the image more problematic, since the CS algorithms require the basis of sampling functions to be orthogonal. While orthogonalization of a basis of continuous functions is rather straightforward, it cannot be obtained with the binarity constraint placed on the functions. To cope with this problem we use matrices precalculated with the singular value decomposition, which has been explained in the Methods section.
We have analyzed the influence of the binarization procedure on the quality of the reconstruction of compressively measured images, and we have found it negligible in the case of image recovery with the use of CS optimization method. The reconstruction quality, measured with the peak signal-to-noise (PSNR) criterion is shown in Fig. 4 as a function of compression. As an alternative to compressive sensing image reconstruction methods we also make use of a direct recovery from the precalculated pseudoinverse of the measurement matrix. The mathematical details of both kinds of methods are briefly summarized in the Methods section. The CS-based recovery offers better quality of the image reconstruction than the pseudoinverse method, especially when binarized sampling patterns are used for image acquisition. However, it can not be obtained in real time. On the other hand, the pseudoinverse-based recovery with the precalculated pseudoinverse matrix requires only the evaluation of a single matrix-vector multiplication, and therefore is very fast. For images of the size of sampled at the compression rate of a few percent, the reconstruction stage is faster than the measurement with the DMD, and may be also realized on-the-fly in parallel with the measurement.
Our single-pixel camera set-up shown in Fig. 5 includes a state-of-the-art DMD with pixel resolution and maximum sampling rate of 22 kHz. Signal-to-noise ratio of the measurement is enhanced using the technique of differential photodetection [26, 27, 23]. The signals measured by two large-area photodiodes are then collected at the rate of 17 MS/s and digitized with 14-bit resolution using a PC oscilloscope.
A complete measurement at the resolution of takes s and enables us to reconstruct images with a high quality (see Fig. 6a), restricted only by the imperfections of the experimental set-up and the presence of optical and electronic noise. Compressive measurements take proportionally less time, however the choice of the sampling protocol is crucial for the feasibility of reconstructing the images with a reasonable quality. For instance, at the compression rate of %, sampling with a random set of Walsh-Hadamard functions allows to obtain a reconstruction with PSNR on the order of dB on average (see Fig. 6b). At the same time, using the nonergodic Morlet wavelet-based random binary sampling functions leads to the reconstructions with PSNR of over dB at the same compression rate (Fig. 6d), while the PSNR for sampling with discrete noiselet functions is on the order of dB (Fig. 6c). Additionally, an image reconstruction obtained with the proposed sampling using pseudoinverse method yields an improved image quality over Walsh-Hadamard sampling with PSNR of dB. While the noise level is increased in comparison with the iterative CS reconstruction, the pseudoinverse enables rapid image recovery by a single matrix multiplication. The sampling scheme proposed in this paper is clearly superior to both other methods also for different compression rates (see Fig. 7).
In this work, we proposed a novel random sampling method for single-pixel imaging. It utilizes nonergodic and stationary Morlet-wavelet-based random patterns that may also be binarized for use with binary spatial light modulators. These sampling functions are obtained as a convolution of Morlet wavelets with realizations of white Gaussian noise.
The proposed sampling functions have a rich spatial and frequency content. Individually, each is a realization of a multivariate Gaussian noise with a characteristic feature size, orientation and modulation frequency. Combined together, the sampling functions uniformly probe the feature space spanned over these image features. We have selected a subset of the feature space through the analysis of an image database.
We have tested this kind of sampling with a large variety of images, and the proposed method enabled us to reconstruct these images with a good quality at compression rates of just a few percent. Both theoretical and experimental results show that the proposed sampling is a lot better than random Walsh-Hadamard sampling. It is also better than noiselet sampling.
At such low compression rates it is still possible to use the direct and fast reconstruction method based on the pseudoinverse of the measurement matrix. A direct reconstruction based on a precalculated pseudoinverse matrix may be implemented on-the fly in parallel with image acquisition on a multicore processor. CS-based reconstruction with a better quality requires much longer reconstruction times on the order of seconds.
CS-based image recovery: We calculate the singular value decomposition of the measurement matrix and following use the total variation (TV) image recovery method implemented in the NESTA  numerical package. When the measurement matrix (with , and the compression ratio denoted as ) consists of rows with nonorthogonal sampling functions , it is first decomposed with the singular-value decomposition (SVD) into a product of small
square orthogonal matrix, diagonal matrix and rectangular semiorthogonal complex conjugate transposed matrix , i.e. . In effect the TV method operates on orthogonal matrices, as is required to reach convergence. The mathematical model of the measurement (where is the captured image, is the measurement matrix, and is the compressive measurement) is replaced with (where , and ), with a semiorthogonal matrix .
Pseudoinverse-based image recovery: The pseudoinverse of the measurement matrix is calculated through the singular value decomposition . For the Morlet-based random sampling functions, the measurement matrix and its pseudoinverse have been precalculated before the measurement. In effect, image recovery has been based on a simple matrix-vector multiplication , which even for a very large matrix takes a fraction of a second to calculate. Since , it is possible to calculate this expression on-the-fly during the measurement, as subsequent components become available.
Binarization of the Morlet-based random sampling functions: The binarization is based on testing the sign of the real-valued functions with the Heaviside step function , i.e. . Additionally, a constant function has been always included as well to measure the mean value of the image.
Differential photodetection: the two states of the DMD mirrors direct the reflected light at two different angles. Then two photodiodes measure both and at the same time. Their difference is used to eliminate the influence of background light and intensity fluctuations of the source from the measurement .
Peak signal to noise ratio: we use a standard definition of the PSNR for the noisy image and reference image , , where is the mean square error.
Discrete noiselet functions : Noiselet sampling is a lot less popular than Walsh-Hadamard sampling so we include the definition of discrete noiselets. Let denote a Hadamard or noiselet transformation matrix whose rows consist of the basis functions. These matrices may be defined recursively as where denotes the Kronecker product, , and for Walsh-Hadamard matrices, and for noiselet matrices. Apart from the normalization, Hadamard basis consist of binary values , while noiselet basis, depending on , consist of values with when
is an odd power of, and when is an even power of . In the first case the real and imaginary parts of noiselet functions are binary, and in the second the sum and difference of their real and imaginary parts are binary . Two-dimensional transforms are obtained through the Kronecker product of one-dimensional transforms, i.e. .
-  Candes, E. J. & Wakin, M. B. An introduction to compressive sampling. IEEE Sign. Proc. Mag. 25, 21–30 (2008).
-  Duarte, M. et al. Single-pixel imaging via compressive sampling. IEEE Sign. Process. Mag. 25, 83–91 (2008).
-  Bian, L. et al. Multispectral imaging using a single bucket detector. Sci. Rep. 6, 24752 (2016).
-  Li, Z. et al. Efficient single-pixel multispectral imaging via non-mechanical spatio-spectral modulation. Sci. Rep. 7, 41435 (2016).
-  Duran, V., Clemente, P., Fernandez-Alonso, M., Tajahuerce, E. & Lancis, J. Single-pixel polarimetric imaging. Opt. Lett. 37, 824 (2012).
-  Soldevila, F. et al. Single-pixel polarimetric imaging spectrometer by compressive sensing. Applied Physics B 113, 551–558 (2013).
-  Li, J., Li, H., Li, J., Pan, Y. & Li, R. Compressive optical image encryption with two-step-only quadrature phase-shifting digital holography. Opt. Commun. 344, 166–171 (2015).
-  Ramachandran, P., Alex, Z. C. & Nelleri, A. Compressive Fresnel digital holography using Fresnelet based sparse representation. Opt. Commun. 340, 110–115 (2015).
-  Watts, C. M. et al. Terahertz compressive imaging with metamaterial spatial light modulators. Nature Photon. 8, 605–609 (2014).
-  Sun, B. et al. 3D computational imaging with single-pixel detectors. Science 340, 844 (2013).
-  Sun, M. et al. Single-pixel three-dimentional imaging with time-based depth resolution. Nat. Commun. 7, 12010 (2016).
-  Li, L., Xiao, W. & Jian, W. Three-dimensional imaging reconstruction algorithm of gated-viewing laser imaging with compressive sensing. Appl. Opt. 53, 7992 (2014).
-  Durán, V. et al. Compressive imaging in scattering media. Opt. Express 23, 14424–14433 (2015).
-  Baraniuk, R. Compressive sensing. IEEE Sign. Proc. Mag. 24, 118–121 (2007).
-  Romberg, J. Imaging via compressive sampling. IEEE Sign. Proc. Mag. 25, 14–20 (2008).
-  Eldar, Y. C. Sampling theory, Beyond Bandlimted Systems (Cambridge Univ. Press, 2015).
-  Phillips, D. et al. Adaptive foveated single-pixel imaging with dynamic supersampling. Sci. Adv. 3, e1601782 (2017).
-  Sun, M. J., Meng, L. T., Edgar, M. P., Padgett, M. J. & Radwell, N. A Russian dolls ordering of the hadamard basis for compressive single-pixel imaging. Sci. Rep. 7, 3464 (2017).
-  Huo, Y., He, H. & Chen, F. Compressive adaptive ghost imaging via sharing mechanism and fellow relationship. Appl. Opt. 55, 3356 (2016).
-  Aßmann, M. & Bayer, M. Compressive adaptive computational ghost imaging. Sci. Rep. 3, 1545 (2013).
-  Candes, E. & Romberg, J. Sparsity and incoherence in compressive sampling. Inverse Probl. 23, 969 (2007).
-  Coifman, R., Geshwind, F. & Meyer, Y. Noiselets. Appl. Comput. Harm. Anal. 27–44, 10 (2001).
-  Pastuszczak, A., Szczygieł, B., Mikołajczyk, M. & Kotyński, R. Efficient adaptation of complex-valued noiselet sensing matrices for compressed single-pixel imaging. Appl. Opt. 55, 5141–5148 (2016).
-  Daugman, J. G. Uncertainity relation for resolution in space, spatial frequency, and orientation optimized by two-dimensional visual cortical filters. J.Opt.Soc.Am.A 2, 1160 (1985).
-  Zhang, C., Guo, S., Cao, J., Guan, J. & Gao, F. Object reconstitution using pseudo-inverse for ghost imaging. Opt. Express 22, 30063–30073 (2014).
-  Sun, B., Welsh, S. S., Edgar, M. P., Shapiro, J. H. & Padgett, M. J. Normalized ghost imaging. Opt. Express 20, 16892 (2012).
-  Yu, W. et al. Complementary compressive imaging for the telescopic system. Sci. Rep. 4, 5834 (2014).
-  Becker, S., Bobin, J., & Candès, E. J. Nesta: a fast and accurate first-order method for sparse recovery. SIAM J. Imaging Sci. 4, 1–39 (2009).
The authors acknowledge financial support from the Polish National Science Center (grant No. UMO-2014/15/B/ST7/03107).
Author contributions statement
K.C. introduced the Morlet wavelet-based sampling functions and conducted numerical analysis. R.K. proposed the idea of nonergodic random sampling. A.P. constructed the optical set-up and conducted the measurements. All authors reviewed the manuscript.
The files generated in this work are available from the corresponding author on reasonable request.
Competing interests: The authors declare that they have no competing interests.