1 Introduction
Compressed sensing (CS) [CandesIntro, B_MathCs] is a technique of recovering a signal from an incomplete measurement based on an assumption that the signal has a sparse representation in a certain domain, for instance in some wavelet basis. In optics, CS has been initially applied for computational ghost imaging [Shapiro2008, Bromberg2009, Sun2012] giving rise to a novel image acquisition technique often referred to as singlepixel camera (SPC) [Baraniuk2008], which allows for capturing images with a sole bucket detector rather than with a highresolution array of detectors. This architecture opens the way for economic electrooptical imaging systems for infrared wavelengths [Gehm2015_EO_IR], as well as for imaging in more exotic ranges of electromagnetic radiation, such as terahertz [Baraniuk_THz2008, Padilla_THz2014] or millimeter waves [Qiao2015_mm_hologr]. Full color imaging [Dinozaury2013], spectral imaging [Gehm2007_spectralimaging], Stokes polarimetric imaging [OL_37_824_Duran], and imaging of threedimensional objects [Busck2004_laser3d, Li2015_laser3d, Javidi2012, Sun2013] have been also demonstrated. The range of research directions related to the use of CS expands rapidly, including lidar imaging [Lidar2012], compressive holography [Lancis2013_hologr, Javidi2013_hologr], sparse subwavelength imaging [Gazit2009_subwave, Eldar2011_subwave]
, compressive pattern recognition
[Nasze_statki_arcxiv], rapid MRI diagnostics [CS_MRI_2008, Weizman2015_MRI] and imaging through scattering media such as a biological tissue [Lancis2015]. Recently, a continuous realtime video recording at Hz with SPC has been also reported [Padget_sci_rep_2015].The common element of CSbased imaging techniques is the use of spatially modulated illumination or aperture. This modulation is frequently achieved with binary spatial light modulators (SLM), for instance with micromirror devices (DMD). The choice of the sampling functions displayed by the SLM depends on the basis, in witch the sampled image has a sparse representation. Indeed, the minimal number of measurements necessary to recover the image is a function of the second power of the coherence between these two bases. As most of the realworld images are compressible in the wavelet domains, the application of sensing matrices incoherent with wavevelets into CSbased imaging systems is essential. Presently, the sensing matrices are usually either based on Hadamard matrices or generated randomly. The Hadamard matrices are both easy to calculate and to display on a binary SLM, however their coherence with commonly used families of wavelets is rather high. The randomly generated patterns are acceptably incoherent with most of the image compression bases, however the necessity of storing the huge sensing matrix in the computer memory during the reconstruction of the image sets considerable limits on the resolution of the imaging system.
These problems are both overcome by sampling the images with noiselet functions [Coifman]. Discrete noislets take values from fourelement complex sets, and importantly for CS with realworld images, are perfectly incoherent with Haar wavelet basis. A unitary noiselet matrix is fast to calculate. A way to calculate the fast noiselet transform (FNT) stems directly from the matrix definition. These properties, which we will overview in more depth in Section 3, decide upon the interest in the use of noiselets for SPC and are our motivation for the present work.
In this paper we focus on the efficient application of noiselet functions for SPC architectures. In a SPC setup with incoherent illumination, light carries a nonnegative real intensity signal. Together with a binary spatial modulation which is respectively represented with real binary functions, the use of more general complex CS measurement matrices is not straightforward. We propose an efficient method which allows to determine complex noiselet coefficients from measurements with nonnegative binary sampling patterns.
2 Compressed imaging
A schematic view of a typical singlepixel sparse imaging setup with incoherent illumination is illustrated in Fig. 1. The system consists of a light source, a SLM, imaging lenses, and a bucket detector. The SLM is used for illuminating the object plane with structured light consisting of a series of binary patterns. Simultaneously, the combined images of the patterns and the object are projected onto a singlepixel detector, which measures their total intensities. Following, the results are digitised, stored and further processed using a PC.
Let
be a vector of length
containing all the measurements captured by the bucket detector. Then, every measurement () within is a dot product of a vector representing the brightness of the consecutive pixels of the th pattern displayed by the SLM and of a vector representing the reflectance (or transmittance) of the corresponding pixels of the image placed in the object plane:(1) 
A matrix , whose rows consist of all the patterns , is called a measurement or a sensing matrix.
In order to reconstruct the original image from the set of measurements , one needs to solve a system of linear equations:
(2) 
In the case of undersampled image sensing, i.e. when , Eq. (2) has an infinite number of solutions. However, the reconstruction of the image is still possible, provided that the image has a sparse representation in a certain compression (sparsity) basis :
(3) 
where is a vector of coefficients of in the basis . The vector can be either literally sparse, containing only a small number of nonzero elements, or all the elements apart from the largest ones can be negligible but still nonzero. In both cases the reconstruction of the image is performed by solving the basis pursuit optimisation problem [CandesIntro]:
(4) 
or, in case of noisy data acquisition, the basis pursuit denoise (BPDN) optimisation problem:
(5) 
where stands for the norm of a vector and represents the level of noise present in the signal. Alternative reconstruction approaches have been also successfully exploited, including minimisation of a certain quasinorm for [Chartrand2007, Chartrand2008] or minimisation of the total variation (TV) [Chan2006_TV, Needell2013_TV, Ma2008_MRI_TV].
The sufficient number of measurements required to collect enough data to reconstruct the original image without distortion satisfies the following inequality [CandesRomberg]:
(6) 
in which is a small constant, is the number of relevant coefficients in the compressed image , and is the mutual coherence of the sensing matrix and the sparsity basis defined as:
(7) 
where (for and ) stand for the row vectors of the matrices and respectively. Therefore, to recover the original image from the least possible number of measurements, it is crucial to choose the sensing matrix in such a manner, that the mutual coherence is kept as small as possible. In other words, the sensing matrix should be almost completely incompressible in the basis .
3 Noiselet matrices
In 2001 a family of functions was introduced, named noiselet functions [Coifman], which is perfectly incoherent with the Haar wavelet basis (mutual coherence between noiselet and Haar orthonormal basis equals ). Since most of the reallife images are well compressible in the Haar wavelet domain, noiselets are then a good candidate for constructing an efficient sensing matrix.
Another advantage of the discrete noiseletbased sensing matrices over e.g. the Gaussian random matrices is that they are defined using a recursive formula based on the Kronecker product (a similar formula was introduced in [Tuma_Hurley_noiselets]):
(8)  
where are unitary matrices whose dimension is a power of two (for ). It is worth mentioning, that Eq. (8) defines the matrices of both onedimensional and twodimensional noiselet transforms. Indeed, the 2D noiselet transform of a matrix takes the form:
(9) 
where denotes vectorisation of matrix obtained by stacking all columns of into a single column vector. The second equality in Eq. (9) is a well known property of the Kronecker product: for matrices , whereas the last equality is a straightforward consequence of the associativity of the Kronecker product, which is recursively used to construct the noiselet matrices (see Eq. (8)).
Therefore, the operation of multiplication by a noiselet matrix is replaceable with a fast transform with a computational complexity of for both onedimensional and twodimensional transforms. Additionally, the huge noiselet matrix is actually never constructed during the evaluation of the matrixvector product and hence does not need to be stored in computer memory.
Ignoring the normalizing factors, the elements of the noiselet matrices take discrete values from one of two 4element sets, depending on the parity of the parameter defining the size of the matrix :
(10)  
where are indices of an element of the matrix .
Owing to these properties, noiselet matrices have an excellent potential of being used as the sensing matrices in CS imaging systems. However, displaying complexvalued patterns with the use of the DMD is impossible, or at least not straightforward. A solution to this problem was proposed recently [zhao2015], suggesting to divide the complexvalued sensing matrix into four separate matrices, namely: 1. the positive real part, 2. the negative real part, 3. the positive imaginary part, and 4. the negative imaginary part, and to perform the measurements with each of these matrices independently. Then, the collected data may be synthesized into an equivalent of a single sequence of measurements obtained with the complexvalued sensing matrix. This method however, suffers from that the number of snapshots taken by the SPC is fourfold larger than the number of samples actually measured, unnecessarily prolonging the data acquisition time. Instead, we want to replace the complex noiselet functions with the same number of real binary functions, but to retain the incoherence of the basis with Haar wavelets.
In conclusion we note that the definition of discrete noiselet transform is very similar to that of WalshHadamard transform, which also shares the same form for onedimensional and twodimensional case, has a similar fast calculation method, and is also commonly used in CS applications especially that it is readily binary. However, the mutual coherence of WalshHadamard matrices with Haar wavelets is much larger.
4 Efficient method of displaying noiseletbased sensing matrices
We propose an efficient method of displaying noiseletbased sensing matrices with the use of a DMD for the purpose of compressive imaging. The method allows for obtaining complexvalued measurement samples by modulating the object with exactly real binary patterns. We note, that in general the DMDs enable displaying also grayscale patterns with multiple intensity levels by flickering of the micromirrors, however the use of binary patterns is the most efficient in terms of the frequency of pattern exposure and the stability of the displayed images.
The proposed procedure takes the following general form:

The sensing matrix is composed of rows chosen randomly from a complexvalued noiselet matrix .

A real binary matrix (further in this paper called a pattern matrix) is defined, whose rows are linear functions of the rows of matrix .

A series of measurements is taken by the SPC using the row vectors of the matrix as the patterns displayed by the DMD:
(11) 
The complexvalued vector is calculated from the measurements and used for reconstructing the image .
Real part 

Imaginary part 

Sum of real and imaginary parts 

Difference of real and imaginary parts 
In order to establish the form of the matrix we exploit several properties of the noiselet matrices. First we note, that for odd values of parameter , real and binary patterns are obtained immediately from the real and imaginary part of the noiselet matrix (see Eq. 10). For even values of the parameter , the real and the imaginary part of a noiselet matrix are triplevalued, however binary patterns are obtained from their sum and difference instead. For better illustration, two examples of a noiselet matrix with either even or odd value of the parameter are presented in Table 1.
However, thus obtained set of patterns suffers from two important shortcomings: 1. the number of realvalued patterns to be displayed by the DMD is twice as large as the number of complexvalued patterns, which they represent and 2. the patterns consist of both positive and negative values. We shall address these issues consecutively.
We begin with the first. This problem may be solved by exploiting the symmetry of the noiselet matrices. Noiselet matrices, similarly to discrete Fourier transform matrices, obey the reflection symmetry of the form:
(12) 
where the symbol denotes the complex conjugate. In other words, each element taken from the upper half of a noiselet matrix (see Fig. 2(a)) is a complex conjugate of the respective mirror element taken from the lower half of the matrix.
Indeed, suppose that Eq. (12) is fulfilled for a noiselet matrix . Then, by induction, for the matrix (see Fig. 2(b)), the elements taken from the left half of the matrix satisfy :
(13) 
where the first and the last equality result from the recursive definition of the matrix (Eq. (8)) and the third equality results from Eq. (12).
Similarly, for the elements from the right half of the matrix :
(14) 
Owing to this property, a pair of patterns consisting of the real and of the imaginary part of a single noiselet (i.e. of a single row vector taken from a noiselet matrix) contains full information about two mirror noiselets taken from the upper and the lower half of the noiselet matrix respectively. Therefore, we propose to choose the sensing matrix in such a manner that it consists of rows randomly picked from the upper half of a noiselet matrix and the mirror rows taken from the lower half of the matrix. Then, the number of unique pattens composed of the real and the imaginary parts (or of their sum and difference) of the rows of matrix equals .
The second problem, concerning patterns consisting of both positive and negative values, is resolved by applying additional rescaling to the patterns. The patterns are already binary, therefore the necessity of displaing negativevalued pixels may be omitted by replacing them with zeros instead. Thus obtained patterns consist of values and only and they are straightforward to be displayed by a DMD. The cost of this rescaling is the necessity of capturing a single additional measurement, independently of the number of patterns or of the size of image . To justify this statement let us first analyse the procedure of restoring the complexvalued measurement vector from the realvalued measurements .
Let us order the rowvectors of the sensing matrix to obey the following formula:
(15) 
The patterns, i.e. the rows of the proposed pattern matrix , take one of two forms, depending on the parity of the parameter defining the size of the sampled image ():
for odd  (16)  
for even  
In both cases, all the patterns stored in the matrix are real and consist of zeros and ones only.
Simultaneously, the row vectors of the sensing matrix are expressed in terms of the patterns by the following equations:
for odd  (17)  
for even  
Finally, the equivalents of the complexvalued measurements are restored from the actual realvalued measurements using the following formulas:
for  
for  
(18) 
and denotes a vector of length with all elements equal to . The term , which occurs in both cases in Eq. (18), is the consequence of the rescaling applied to the patterns in order to obtain only nonnegative values. In physical interpretation, this term represents a measurement of the total intensity of the image without modulating it with any pattern. This single additional measurement is necessary for reverting the scaling during retrieval of . Therefore, exactly measurements with real binary and nonnegative patterns are required in order to restore samples corresponding to measurements taken with a complexvalued and nonbinary noiseletbased sensing matrix.
In practical experimental conditions this number of measurements may be increased. For instance, in order to eliminate the background ‘dark pixel’ signal resulting from the light reflected from the DMD matrix when all the mirrors are in the offstate, it is necessary to measure that dark signal in the calibration stage, and then to subtract it from the measurements. More likely, the intensity of the lightsource may vary with time. An additional detector could be used for the normalization of measurements [Sun2012]. In another approach, a differential measurement with complementary [Yu2015_complementary3D] binary masks increases the signaltonoiseratio, eliminates background dark signal, and accounts for intensity variations but actually doubles the number of measurements.
5 Modification of the fast noiselet transform for realtime generation of the noiseletbased patterns
In the following section we propose a modification of the fast noiselet transform for the purpose of efficient generation of the noiseletbased real binary patterns introduced in Section 4. The onedimensional and twodimensional noiselet transform defined by Eq. (8) operates on complex numbers. We propose a similar procedure, allowing for a direct generation of the real binary patterns using only operations of summation and subtraction on integer variables. Moreover, by utilising the bit representation of a bit integer, a bundle of up to patterns may be generated using the same number of arithmetic operations as in the case of generating only a single pattern (the additional two bits are used in the intermediate calculations). The efficiency of the proposed method allows for generating the patterns during the time of the experiment, without any preparations beforehand.
The modified noiselet transform matrix (where ) relates to the noiselet matrix as follows:
(19) 
and it satisfies a recursive formula similar to that, which defines matrices :
(20)  
We note, that the elements of thus defined matrix belong to a single 4element set, independently on the parameter :
(21) 
To derive the modified noiselet transform , let us rewrite Eq (20) into a more explicit form. We define an auxiliary matrix:
(22) 
Then, the modified noiselet matrix takes a form:
(23)  
where
stands for an identity matrix of size
, and the second equality results from the mixedproduct property of the Kronecker product: . Each factor in the final expression in Eq. (23) of a form for , is a block diagonal matrix with same blocks of the form:(24) 
Therefore, the modified noiselet transform of a vector (i.e. the multiplication of the vector by the matrix ) is equivalent to dividing the vector into blocks, each of length , and multiplying each of them by the matrix (24). The procedure is then repeated times for sizes of the blocks defined by consecutive values of . The product of a single partial vector by the matrix (24) has a simple form:
(25) 
which, by keeping the real and the imaginary parts of the vectors as separate variables, does not require any implementation of complex numbers. Moreover, the only arithmetical operations present in the transform are summations and subtractions applied to the real and the imaginary part of the vector . Therefore, for vector consisting only of integer values, the transform may be implemented purely on integer types, which greatly increases its efficiency.
The real binary patterns introduced in the previous section are obtained directly from the modified noiselet transform. Indeed, if a complex sensing pattern is chosen as the th row of the noiselet matrix , then the pair of its equivalent real binary patterns , defined by Eq. (19) is expressed by the modified noiselets as follows:
(26)  
where
(27)  
indicates the th row of the matrix and stands for the reminder after division of parameter by . The signs in the expressions for , depend on the value of the complex argument in Eq. (19).
Now, let us explain the operation of the packed modified noiselet transform. A single pattern or is efficiently obtained by applying the modified noiselet transform to a unit vector , whose all elements apart from the th one are zeros:
(28)  
The operation required to obtain a bundle of different binary patterns (or with interchanging ) encoded into the consecutive bitplanes of the bit representation of the integer variables is derived from Eq. (28):
(29)  
where is a vector composed of a bundle of unit vectors encoded into the consecutive bitplanes of its integer elements. From Eqs. (28), (29) one may easily determine that the computational complexity of calculating the whole bundle of patterns is identical as in the case of calculating a single pattern. The number of patterns, which may be calculated with a single run of the modified noiselet transform, is limited only by the bitwidth of the integer variables used in the implementation of the transform.
To illustrate the efficiency of the modified noiselet transform, we note that our C++ implementation of the transform generates around bundles of patterns of resolution per second or bundles of patterns of resolution per second on a midrange laptop. For patterns bundled into packages of 24, as is the case in our experiment, this is respectively or patterns per second. This speed already includes all the operations on the computer graphical objects and hardware necessary to display the patterns.
6 Experimental Results
We will now demonstrate the experimental results of reconstructing an image captured by a SPC. We apply binary sampling equivalent to noiselet sampling according to the procedure proposed in Sections 4 and 5.
6.1 Design of the experimental setup
The schematic of our experimental system matches the one shown in Fig. 1. We use a DMD light modulator (TI DLP LightCrafter 4500) integrated with an RGB LED light source and optical lens to project the sampling patterns onto the object plane. The DMD consists of square micromirrors organised into a diamond grid, i.e. a square grid rotated by with respect to the boundaries of the DMD. Out of this array, we use a square subarea of pixels oriented along the edges of the pixels, so that the boundaries of the area form straight lines. All the patterns are transformed in order to be displayed in this area, which ensures accurate projection of the patterns into a square grid, without distortions resulting from the diamond pixel layout of the device. A similar approach was previously reported in [Lancis_dual_mode_2016]. The singlepixel detector consists of a photodiode integrated with an onchip transimpedance amplifier (TI OPT101P) with peak sensitivity wavelength of nm. The analog signal is digitized with a bit A/D converter (NI USB6003 100kS/s multifunction DAQ) and streamed to a PC via USB port. Data acquisition is controlled with a LabView routine.
The DLP displays binary patterns consisting of arbitrarily selected bitplanes of 24bit RGB images. These images may be either first stored in the internal flash memory or streamed via HDMI video port. The latter solution is more convenient, since it does not impose memory restrictions and allows for a more flexible choice of patterns, by generating them in real time during the measurements instead of preparing and uploading them beforehand. We have developed a dedicated fast routine in C++ based on modified noiselet transform (see Section 5) to generate noiseletbased patterns (Eq. (26)) bundled together in packages of (with one bit left for synchronization) into the bitrepresentations of RGB images. Such images are transferred through HDMI port to the DLP and then displayed in sequence, bit by bit. Our system is capable of displaying binary patterns at the maximum rate of Hz ( bits Hz of the video rate). However, the actual speed of displaying the patterns in the experiment is set to Hz in order to preserve high signal to noise ratio and to ensure stability of the data acquisition.
The accuracy of the measurements is an important issue for the SPC imaging. Since all the noiseletbased binary patterns of the same resolution have the same total brightness, the standard deviation of the measurements taken with different patterns is usually at least two orders of magnitude smaller than their average value. Therefore, in order to increase the signal to noise ratio and to reduce the influence of both the background signal from the dark pixels and the light intensity variations over time, we use the technique of differential measurements, involving displaying both the patterns and their binary negations. In a different SPC design, similar effect could be obtained without doubling the number of measurements, by introducing a second detector to measure synchronously the total intensity of light reflected from the DMD.
Finally, we recover the image from the measured data by solving the BPDN problem (see Eq. (5)) with the use of the SPGL1 package [SPGL1art2011]. In the worst case scenario of recovering an image from measurements with 50% elements of the noiselet basis, the optimisation takes approximately s for an image with resolution (in Matlab, using a PC with a single eightcore processor). In the case of reconstructing the image from the entire noiselet basis, a straightforward approach of calculating the inverse noiselet transform may be applied, reducing the time of recovery to s.
6.2 Results
The original image used for the experiment is presented in the Fig. 3(a). The image is well compressible in the Haar wavelets, as shown in Fig. 3(b), with only 24% of nonzero Haar coefficients. Further lossy compression is also possible, with the peak signal to noise ratio (PSNR, see Eq. (30)) of the compressed image on the order of dB when 8% of the largest Haar coefficients are preserved.
The image is sampled with the noiseletbased real binary patterns of the resolution . The recovery of the image is repeated using several different values of the sampling ratio from the range between 5% and 100%, where 100% corresponds to the sensing matrix consisting of all possible noiselets from a single noiselet matrix. Several examples of reconstructed images obtained with different values of the sampling ratio are presented in Figs. 3(c)(f). Additionally, in Fig. 3(g) we present how the sampling ratio influences the accuracy of reconstruction of the image. To measure the quality of the reconstruction, we use the PSNR defined as:
(30) 
where represents distribution of brightness in the reconstructed image, stands for the peak brightness of the original image, and is the mean squared error of the reconstructed image as compared to the original one. Before calculating the PSNR, both images have been registered and normalised in order to match their intensity levels. We note that PSNR is not a deterministic function of the sampling ratio  it depends on the specific patterns used in the measurement or simulation. Therefore, in Fig. 3(g) we show also the standard deviation of PSNR calculated over a number of randomly chosen sets of patterns with the same value of the sampling ratio.
We define the experimental noise as the difference between the measured and calculated values of . The standard deviation of the experimental noise is on the order of
of the mean value of the measured signal and it is by approximately two orders of magnitude lower than the peaktopeak amplitude of the measurements taken with different patterns. By introducing additive Gaussian white noise to the theoretical model of the measurement, we obtain a good agreement between the simulated compressive measurement and the actual experimental results (see Fig.
3(g)).7 Conclusions
We have proposed theoretically and validated experimentally an efficient method of using complexvalued and nonbinary noiselet functions for object sampling in singlepixel cameras with binary spatial light modulators and incoherent illumination. Minimal mutual coherence of discrete noiselets and Haar wavelets makes this pair of bases an essential choice for the sensing and compression matrices in compressed sensing with singlepixel detectors. Indeed, most realworld images are compressible in the Haar basis. The proposed method allows to determine noiselet coefficients from binary sampling measurements. Moreover, we have proposed a modification to the complex fast noiselet transform, which enables computationallyefficient generation of the binary noiseletbased patterns using only operations of summation and subtraction on integer variables. Further acceleration is obtained by utilising the bit representation of a bit integer to calculate a bundle of up to patterns without any additional computational cost as compared with generating only a single pattern. The efficiency of the proposed method allows for generating patterns in real time on a PC or even on a singleboard computer.
Funding Information
Funding. We acknowledge financial support from the National Science Centre, Poland grant UMO2014/15/B/ST7/03107 and European Union Seventh Framework Programme (FP7/20072013, grant agreement no 316244).
Comments
There are no comments yet.