Imaging With Nature: Compressive Imaging Using a Multiply Scattering Medium

The recent theory of compressive sensing leverages upon the structure of signals to acquire them with much fewer measurements than was previously thought necessary, and certainly well below the traditional Nyquist-Shannon sampling rate. However, most implementations developed to take advantage of this framework revolve around controlling the measurements with carefully engineered material or acquisition sequences. Instead, we use the natural randomness of wave propagation through multiply scattering media as an optimal and instantaneous compressive imaging mechanism. Waves reflected from an object are detected after propagation through a well-characterized complex medium. Each local measurement thus contains global information about the object, yielding a purely analog compressive sensing method. We experimentally demonstrate the effectiveness of the proposed approach for optical imaging by using a 300-micrometer thick layer of white paint as the compressive imaging device. Scattering media are thus promising candidates for designing efficient and compact compressive imagers.

A cquiring digital representations of physical objects -in other words, sampling them -was, for the last half of the 20 th century, mostly governed by the Shannon-Nyquist theorem. In this framework, depicted in Fig. 1(a), a signal is acquired by N regularly-spaced samples whose sampling rate is equal to at least twice its bandwidth. However, this line of thought is thoroughly pessimistic since most signals and objects of interest are not only of limited bandwidth but also generally possess some additional structure 1 . For instance, images of natural scenes are composed of smooth surfaces and/or textures, separated by sharp edges.
Recently, new mathematical results have emerged in the field of Compressive Sensing (or Compressed Sensing, CS in short) that introduce a paradigm shift in signal acquisition. It was indeed demonstrated by Donoho, Candès, Tao and Romberg 2-4 that this additional structure could actually be exploited directly at the acquisition stage so as to provide a drastic reduction in the number of measurements without loss of reconstruction fidelity.
For CS to be efficient, the sampling must fulfill specific technical conditions that are hard to translate into practical design guidelines. In this respect, the most interesting argument featured very early on in [2][3][4] is that a randomized sensing mechanism yields perfect reconstruction with high probability. As a matter of convenience, hardware designers have created physical systems that emulate this property. This way, each measurement gathers information from all parts of the object, in a controlled but pseudo-random fashion. Once this is achieved, CS theory provides good reconstruction guarantees.
In the past few years, several hardware implementations capable of performing such random compressive sampling were introduced [5][6][7][8][9][10][11][12][13] . In optics, these include the single pixel camera 6 , which is depicted in Fig. 1(b), and uses a digital array of micromirrors (abbreviated DMD) to sequentially reflect different random portions of the object onto a single photodetector. Other approaches include phase modulation with a spatial light modulator 10 , or a rotating optical diffuser 13 . The idea of random multiplexing for imaging has also been considered in other domains of wave propagation. CS holds much promise in areas where detectors are rather complicated and expensive such as the THz or far infrared. In this regards, there have been proposals to implement CS imaging procedures in the THz using random pre-fabricated masks 5 , DMD or SLM photo-generated contrast masks on semi-conductors slabs 14 and efforts are also pursued on tunable metamaterial reflectors 15 . Recently, a carefully engineered metamaterial aperture was used to generate complex RF beams at different frequencies 8 .
However, these CS implementations come with some limitations. First, these devices include carefully engineered hardware designed to achieve randomization, via a DMD 6 , a metamaterial 8 or a coded aperture 11 . Second, the acquisition time of most implementations can be large because they require the sequential generation of a large number of random patterns.
In this work, we replace such man-made emulated randomization by a natural multiply scattering material, as depicted in Fig. 1(c). Whereas scattering is usually seen as a time-varying nuisance, for instance when imaging through turbid media 16 , the recent results of wave control in stable complex material have largely demonstrated that it could also be exploited, for example so as to build focusing systems that beat their coherent counterparts in terms of resolution 17,18 . Such complex and stable materials are readily available in several frequency ranges -they were even coined in as one-way physical functions for hardware cryptography 19 . In the context of CS, such materials perform an efficient randomized multiplexing of the object into several sensors and hence appear as analog randomizers. The approach is applicable in a broad wavelength range and in many domains of wave propagation where scattering media are available. As such, this study is close in spirit to earlier approaches such as the random reference structure 20 , the random lens imager 7 , the metamaterial imager 8 , or the CS filters proposed in 21 for microwave imaging. They all abandoned digitally controlled multiplexors as randomizers. Still, we go further in this direction and even drop the need for a designer to craft the randomizer.
Compressive sampling with multiply scattering material has several advantages. First, it has recently been shown that they have an optimal multiplexing power for coherent waves 22 , which consequently makes them optimal sensors within the CS paradigm. Second, these materials are often readily available and require very few engineering. In the domain of optics for example, we demonstrate one successful implementation using a 300 mm layer of Zinc Oxide (ZnO), which is essentially white paint. Third, contrarily to most aforementioned approaches, this sensing method provides the somewhat unique ability to take a scalable number of measurements in parallel, thus with a potential of strongly reducing acquisition time. In practice, if 500 samples are required to reconstruct a given image using CS principles, this imaging framework allows their acquisition at once on 500 independent sensors, whereas state-ofthe-art systems such as the single pixel camera require a sequence of 500 random patterns on the DMD.
On practical grounds, the use of a multiply scattering material in CS raises several ideas that we consider in this study. First, the random multiplexing achieved through multiple scattering must be measured a posteriori, since it is no longer known a priori as in engineered random sensing. This calibration problem has been the topic of recent studies in the context of CS 23 and we propose here a simple least squares calibration procedure that extends our previous work 24,25 . Second, the use of such a measured Transmission Matrix (TM) induces an inherent uncertainty in the sensing mechanism, that can be modeled as noise in the observations. As we show both through extensive simulations and actual experiments, this uncertainty is largely compensated by the use of adequate reconstruction techniques. In effect, the imager we propose almost matches the performance of idealized sub-Nyquist random sensing.

Theoretical background
In its simplest form, CS may be understood as a way to solve an underdetermined linear inverse problem. Let x be the object to image, understood as a N 3 1 vector, and let us suppose that x is only observed through its multiplication y by a known measurement matrix H, of dimension M 3 N, we have y 5 Hx. Each one of the M entries of y is thus the scalar product of the object with the corresponding row of H. When there are fewer measurements than the size of the object, i.e. M , N, it is impossible to recover x perfectly without further assumptions, since the problem has infinitely many solutions. However, if x is known to be sparse, meaning that only a few of its coefficients are nonzero (such as stars in astronomical images), and provided H is sufficiently random, x can still be recovered uniquely through sparse optimization techniques 1 .
In a signal processing framework, the notion of structure may also be embodied as sparsity in a known representation 1 . For example, most natural images are not sparse, yet often yield near-sparse representations in the wavelet domain. If the object x is known to have some sparse or near-sparse representation s in a known basis B (x 5 Bs), then it may again be possible to recover it from a few samples, by solving y 5 HBs, provided H and B obey some technical conditions such as incoherence [1][2][3][4]26 .
In practice, when trying to implement Compressive Sensing in a hardware device, fulfilling this incoherence requirement is nontrivial. It requires a way to deterministically scramble the information somewhere between the object and the sensors. Theory shows that an efficient way to do this is by using random measurement matrices H or HB 2-4 . Using such matrices, it can indeed be shown 26 that the number of samples required to recover the object is mostly governed by its sparsity k, i.e. the number of its nonzero coefficients in the given basis. If the coefficients of the M 3 N measurement matrix are  independent and identically distributed (i.i.d.) with respect to a Gaussian distribution, perfect reconstruction can be achieved with only O(k log(N/k)) measurements 27 . Furthermore, many algorithms are available, for instance Orthogonal Matching Pursuit (OMP) or Lasso 26,28 , which can efficiently perform such reconstruction under sparsity constraints.
Using natural complex media as random sensing devices Our approach is summarized in Fig. 1(c) and its implementation in an optical experiment is depicted in Fig. 2. The coherent waves originating from the object and entering the imaging system propagate through a multiply scattering medium. Within the imager, propagation produces a seemingly random and wavelength-dependent interference pattern called speckle on the sensors plane. The speckle figure is the result of the random phase variations imposed on the waves by the propagation within the multiply scattering sample 29 . Scattering, although the realization of a random process, is deterministic: for a given input, and as long as the medium is stable, the interference speckle figure is fully determined and remains constant. In essence, the complex medium acts as a highly efficient analog multiplexer for light, with an input-output response characterized by its transmission-matrix 24,25 . We highlight the fact that the multiple scattering material is not understood here as a nuisance occurring between the object and the sensors, but rather as a desirable component of the imaging system itself. After propagation, sensing takes place using a limited number M , N of sensors.
Let x and y denote the N 3 1 and M 3 1 vectors gathering the value of the complex optical field at discrete positions before and after, respectively, the scattering material. It was confirmed experimentally 24,25 that any particular output y m can be efficiently modeled as a linear function of the N complex values x n of the input optical field: where the mixing factor h mn [C corresponds to the overall contribution of the input field x n into the output field y m . All these factors can be gathered into a complex matrix [H] mn 5 h mn called the Transmission Matrix (TM), which characterizes the action of the scattering material on the propagating waves between input and output. The medium hence produces a very complex but deterministic mixing of the input to the output, that can be understood as spatial multiplexing. This linear model, in the ideal noiseless case, can be written more concisely as:

y~Hx:
As can be seen, each of the M measurements of the output complex field may hence be considered as a scalar product between the input and the corresponding row of the TM. If multiply scattering materials have already been considered for the purpose of focusing, thus serving as perfect ''opaque lenses'' 17,18 , the main idea of the present study is to exploit them for compressive imaging. In other wavelength domains than optics, analogous configurations may be designed to achieve CS through multiple scattering. For instance, a collection of randomly packed metallic scatterers could be used as a multiply scattering media from the microwave domain up to the far infrared, and the method proposed here could allow imaging at these frequencies with only a few sensors. A similar approach could be used to lower the number of sensors in 3D ultrasound imaging using CS through multiple scattering media.
In our optical experimental setup, we used a Spatial Light Modulator (an array of N 5 1024 micromirrors, abbreviated as SLM) to calibrate the system and also to display various objects, using a monochromatic continuous wave laser as light source.
During a first calibration phase, which lasts a few minutes and needs to be performed only once, a series of controlled inputs x are emitted and the corresponding outputs y are measured. The TM can be estimated through a simple least-squares error procedure, which generalizes the method proposed in 24,25 , as detailed in the supplementary material below. In short, this calibration procedure benefits from an arbitrarily high number of measurements for calibration, which permits to better estimate the TM. It is important to note here that the need for calibration is the main disadvantage of this go through a scattering material (ii) that efficiently multiplexes the information to all M sensors (iii). Provided the transmission matrix of the material has been estimated beforehand, reconstruction can be performed using only a limited number of sensors, potentially much lower than without the scattering material. In our optical scenario, the light coming from the object is displayed using a spatial light modulator. technique, compared to the more classical CS imagers based on pseudo-random projections, which have direct control on the TM. However, this calibration step here involves only standard leastsquares estimation of the linear mapping between input and output of the scattering material 24,25 . In our experimental setup, the whole calibration is performed in less than 1 minute. While we here rely on optical holography to extract complex amplitude from intensity measurements, the TM measurement can be implemented in a simplified way for other types of waves (RF, acoustics, Terahertz), where direct access to the field amplitude is possible. It may not be so straightforward in practical situations when only the intensity of the output is available, and where more sophisticated methods 30 would be required.
After calibration, the scattering medium can be used to perform CS, using this estimated TM as a measurement matrix. Note that, in our experiment, the same SLM used for calibration is then used as a display to generate the sparse objects. This approach is not restrictive as any sparse optical field or other device capable of modulating light could equivalently be used at this stage. As demonstrated in our results section, using such an estimated TM instead of a perfectly controlled one does yield very good results all the same, while bringing important advantages such as ease of implementation and acquisition speed. Hence, even if the proposed methodology does require the introduction of a supplementary calibration step, this step comes at the cost of a few mandatory supplementary calibration measurements rather than at the cost of performance. This claim is further developed in our results and methods sections.
For a TM to be efficient in a CS setup, it has to correctly scramble the information from all of its inputs to each of its outputs. It is known that a matrix with i.i.d Gaussian entries is an excellent candidate for CS 31 and the TM of optical multiple scattering materials were recently shown to be well approximated by such matrices 22 . The rationale for this fact is that the transmission of light through an opaque lens leads to a very large number of independent scattering events. Even if the total transmission matrix that links the whole input field to the transmitted field shows some non-trivial meso-scopic correlations 32 , recent studies proved that these correlations vanish when controlling/measuring only a random partition of input/output channels 22 . In our experimental setup, the number of sensors is very small compared to the total number of output speckle grains and we can hence safely disregard any mesoscopic correlation.
Several previous studies 24, 25 have shown on experimental grounds that TMs were close to i.i.d. Gaussians by considering their spectral behavior, i.e. the distribution of their eigenvalues. As a consistency check, we also verified that our experimentally-obtained TMs are close to Gaussian i.i.d., through a complementary study of their coherence, which is the maximal correlation between their columns with values between 0 and 1. Among all the features that were proposed to characterize a matrix as a good candidate for CS [33][34][35] , coherence plays a special role because it is easily computed and because a low coherence is sufficient for good recovery performance in CS applications [36][37][38][39] , even if it is not necessary 40 . In Fig. 3(a), we display one actual TM obtained in our experiments. In Fig. 3(b), we compare its coherence with the one of randomly generated i.i.d. Gaussian matrices. The similar behavior confirms the results and discussions given in 22,24 , but also suggests that TMs are good candidates in a CS setup, as will be demonstrated below.

Results and discussion
During our experiments, we measured the reconstruction performance of the imaging system, when the image to reconstruct is composed of N 5 32 3 32 5 1024 pixels, using a varying number M of measurements. In practice, we use a CCD array, out of which we select M pixels. These are chosen at random in the array, with an exclusion distance equal to the coherence length of the speckle, in order to ensure uncorrelated measurements. Details of the experiments can be found in the methods section below. For each sparsity level k between 1 and N, a sparse object with only k nonzero coefficients was displayed under P 5 3 different random phase illuminations [Since our SLM can only do phase modulation, we used a simple trick as in 41 to simulate actual amplitude objects, based on two phase- Reconstruction of the sparse objects was then achieved numerically using the M 3 P measurements only. The TM used for reconstruction is the one estimated in the calibration phase. In order to demonstrate the efficiency and the simplicity of the proposed system, we used the simple Multichannel Orthogonal Matching Pursuit algorithm 42 for MMV reconstruction. It should be noted that more specialized algorithms may lead to better performance and should be considered in the future.
Examples of actual reconstructions performed by our analog compressive sampler are shown on Fig. 4. As can be seen, near-perfect reconstruction of complex sparse patterns occur with sensor density ratios M/N that are much smaller than in classical Shannon-Nyquist sampling (M 5 N). An important feature of the approach is its universality: reconstruction is also efficient for objects that are sparse in the Fourier domain.
The performance of the proposed compressive sampler for all sampling and sparsity rates of interest is summarized on Fig. 5, which is the main result of this paper. It gives the probability of successful reconstruction displayed as a function of the sensor density M/N and relative sparsity k/M. Each point of this surface is the average reconstruction performance for real measurements over approximately 50 independent trials. As can be seen, this experimental diagram exhibits a clear ''phase transition'' from complete failure to systematic success. This thorough experimental study largely confirms that the proposed methodology for sampling using scattering media indeed reaches very competitive sampling rates that are far below the Shannon-Nyquist traditional scheme.
The phase transition observed on Fig. 5 appears to be slightly different from the ones described in the literature 31,43 . The main reason for this fact is that this diagram concerns reconstruction under P 5 3 Multiple Measurement Vectors (MMV) instead of the classical Single Measurement Vector (SMV) case. This choice, which proves important in practice, is motivated by the fact that MMV is much more robust to noise than SMV 44 . In order to compare our experimental performance to its numerical counterpart, we performed a numerical experiment whose 50% success-rate transition curve is represented by the dashed green line. The transmission matrix is taken as i.i.d Gaussian. The measurement matrix is estimated with the same calibration procedure as in the physical experiment. Each measurement, during calibration and imaging, is contaminated by additive Gaussian noise of variance 3%. Performance obtained in this idealized situation is close to that obtained in our practical setup, for this level of additive noise.

Conclusion
In this study, we have demonstrated that a simple natural layer of multiply scattering material can be used to successfully perform compressive sensing. The compressive imager relies on scattering theory to optimally dispatch information from the object to all measurement sensors, shifting the complexity of devising CS hardware from the design, fabrication and electronic control to a simple calibration procedure.
As in any hardware implementation of CS, experimental noise is an important issue limiting the performance, especially since it impacts the measurement matrix. Using baseline sparse reconstruction algorithms along with standard least-squares calibration techniques, we demonstrated that successful reconstruction exhibits a clear phase transition between failure and success even at very competitive sampling rates. The proposed methodology can be considered to be a truly analog compressive sampler and as such, benefits from both theoretical elegance and ease of implementation.
The imaging system we introduced has many advantageous features. First, it enables the implementation of an extremely flat imaging device with few detectors. Second, this imaging methodology can be implemented in practice with very few conventional lenses, as in 45 for instance. This is a strong point for implementation in domains outside optics where it is hard to fabricate lenses. Indeed, the concept presented here can directly be used in other domains of optics such as holography, but also in other disciplines such as THz, RF or ultrasound imaging. Third, similarly to recent work on metamaterials apertures, non-resonant scattering materials work over a wide frequency range and have a strongly frequency-dependent response. Fourth, unlike most current compressive sensing hardware, this system gives access to many compressive measurements in a parallel fashion, potentially speeding up acquisition. These advantages come at the simple cost of a calibration step, which amounts to estimate the Transmission Matrix of the scattering material considered. As we demonstrated, this can be achieved by simple input/output mapping techniques such as linear least-squares and needs to be done only once.
While conventional direct imaging can be thought as an embarrassingly parallel process that does not exploit the structure of the scene, in contrast most current CS hardware (such as the single pixel camera) require a heavily sequential process that does take into account the structure of the scene. Our approach borrows from the best of both acquisition processes, in that it is both embarrassingly parallel and takes into account the structure of the scene.