{'fraction_non_alphanumeric': 0.036909236551851625, 'fraction_numerical': 0.012200382032164643, 'mean_word_length': 4.718816067653277, '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': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Intensity-only optical compressive imaging using a multiply scattering material and a double phase retrieval approach

In this paper, the problem of compressive imaging is addressed using natural randomization by means of a multiply scattering medium. To utilize the medium in this way, its corresponding transmission matrix must be estimated. To calibrate the imager, we use a digital micromirror device (DMD) as a simple, cheap, and high-resolution binary intensity modulator. We propose a phase retrieval algorithm which is well adapted to intensity-only measurements on the camera, and to the input binary intensity patterns, both to estimate the complex transmission matrix as well as image reconstruction. We demonstrate promising experimental results for the proposed algorithm using the MNIST dataset of handwritten digits as example images.


INTRODUCTION
From the perspective of image processing, the goal of compressed sensing (CS) is to reconstruct a highresolution image, which is sparse in either the ambient domain or some transform basis, using few incoherent linear projections [1]. Over the past decade, there has been a tremendous amount of work in the field of CS, including analytical reconstruction guarantees as well as developments of new algorithmic approaches that provide efficient methods of solving the reconstruction task [2,3]. However, to date there have been only a handful of engineering projects where optical imagers based on CS have actually been built. Indeed, performing these incoherent, usually random, projections is a highly non-trivial task, requiring innovative hardware solutions. Amongst such imagers, one can cite, without any claim of completeness, several single-pixel imaging systems [4][5][6], a random lens camera [7], and an imaging setup based on a rotating diffuser [8].
The work presented in this paper is built upon a recently developed optical CS setup [9] that uses a multiply scattering medium to effect the random projection operation. The fundamental difference with this approach and most of the CS systems discussed above is that here the random projections are not designed beforehand and then implemented through sophisticated hardware, as in [10], but are based on the natural randomization properties of coherent light multiply scattering through a layer of opaque material. Here, the word "multiply" refers to the fact that the thickness of the material slab is many times larger than the mean free path, ensuring that the light beam is fully scattered without any remaining ballistic photons at the output. If x is the incoming wavefield (the object to be imaged at the input plane), the scattering operation is well modeled by a simple linear operator H, called the transmission matrix. If y is the output wavefield discretized by receptor pixels, then, in the ideal noiseless case, It has been shown that the transmission matrix of a scattering material is statistically identical to an i.i.d. random matrix with a complex Gaussian distribution [11]. The benefits of using such a system for CS imaging are that one does not have to rely on complex engineering solutions to provide the (pseudo-) randomization, and also that, in theory, only one shot is necessary to obtain any desired number of output features; as opposed to the single-pixel camera which intrinsically requires sequential measurements. There exists, however, an obvious price to pay: the necessity of a precise calibration step. Indeed, to be able to use this system as an imaging device, i.e. to estimate x given measurements y, one must have accurate knowledge of the matrix H. This can be accomplished by sending a series of known images, measuring the corresponding outputs, and performing a least-squares estimate of H. The calibration step is conducted by shaping the input wavefront with a Spatial Light Modulator (SLM), which is only used for calibration and display, and is not part of the direct imaging system.
In this paper, we circumvent one major limitation of the previous proof-of-concept system [12]. Since optical sensors (here, a CCD camera) only measure the field intensity |y| 2 , in [12] the input image is phase-modulated using a phase-only SLM, with relative phases 0, π/4, π/2, and 3π/4. Combining the corresponding four output intensity images, one can easily recover the complex field y using a method known as "phase-stepping holography". Furthermore, such phase-only modulated images have a constant intensity. To obtain an image that is sparse in the spatial domain, one has to make the difference between 2 complex phase-only images which only differs arXiv:1510.01098v2 [cs.CV] 25 Jan 2016 by a sparse number of pixels. Therefore, in order to get the complex measurements corresponding to a single sparse image, 8 intensity measurements are required. This significantly slows down both the calibration and the measurement process. Furthermore, a sufficiently fast continuous-phase SLM is a very expensive device, with limited pixel counts. For example, the SLM used in [12] could only display 32 × 32 images.
Here, we investigate the alternative use of a digital micro-mirror device (DMD) as an SLM, as shown in Fig.  1. This has many advantages: DMDs are cheap, fast, and have high pixel counts. However, the main drawback of these binary intensity modulators is that, without additional hardware, one can no longer use phasestepping to measure the complex output field. Instead of using hardware to measure amplitude and phase, we resort to "phase retrieval" in order to estimate the missing phases from intensity-only measurements |y| 2 . It should be noted that, in this framework, phase retrieval must be applied twice successively; first, for the calibration, and second, for the imaging itself. The success of the second step crucially depends on the first one, as every error in estimating H results in multiplicative noise (also called model error) in the imaging step. It should also be noted that the signal-to-noise ratio is relatively poor, thus we favor Bayesian phase retrieval techniques where noise may be explicitly modeled.
The main contributions of this paper are as follows: • A new Bayesian phase retrieval algorithm known as phase retrieval Swept AMP (prSAMP). prSAMP originates from prGAMP [13] and SwAMP [14] and is designed to work with noisy ill-conditioned transmission matrices.
• The experimental demonstration that prSAMP is efficient both for calibration of the non-sparse measurement matrix H using binary inputs, and for intensity-only CS imaging of sparse inputs.
Although our previous studies [12] demonstrate a proofof concept that CS-based imaging can be made with multiply scattering materials, we believe that this one-shot imager represents a very significant step toward real-life applications of these techniques.

THEORETICAL MODELING
Starting from the idealized model of Eq. (1), we formalize the calibration procedure as in [9]. Given P known binary input images of size N = n 1 × n 2 , X ∈ {0, 1} P ×N , and their corresponding intensity measurements on M output pixels, Y ∈ R M ×P + , M independent phase retrieval problems are solved during the calibration step to estimate the transmission matrix H ∈ C M ×N . Each  [9]. A monochromatic laser at 532 nm is expanded by a telescope and illuminates an SLM, here, a Texas Instruments DLP9500 DMD with 1920×1080 pixels. The light beam carrying the image is then focused on a random medium by means of a microscope lens. Here, the medium is a thick (several tens of microns) opaque layer of Zinc Oxide nanoparticles deposited on a glass slide. The transmitted light is collected on the far side by a second lens, passes through a polarizer, and is detected by an AVT PIKE F-100 monochrome CCD camera. Note that the DMD is only for calibration and display and is not part of the imager itself.
calibration problems is formulated as where (·) m indicates the m-th row of corresponding matrix and (·) T is the transpose operator. The process of recovering a signal from only the magnitude of its projections is the goal of phase retrieval [13,15,16]. Apart from additional noise in the measurements, what makes solving Eq. (2) challenging is using binary input patterns, since most well-known phase retrieval methods work well with complex-valued measurement matrices. We have fixed this issue by mixing the ideas of Swept Approximate Message Passing (SwAMP) [14], which demonstrates good convergence properties over ill-conditioned noisy matrices, with the phase retrieval method prGAMP [13]. The new prSAMP algorithm is explained in the next section.
After calibration, the setup can be used as a generalized CS imager with non-linear (intensity) measurements. In this reconstruction phase, the noiseless model becomes y = |Hx|. We use the same prSAMP method, with different priors, to solve both the calibration and reconstruction tasks.

PRSAMP ALGORITHM
In the context of CS, AMP is an iterative algorithm for the reconstruction of a sparse signal from a set of underdetermined linear noisy measurements y = Hx + w, where w ∼ N (0, σ 2 ) [17]. Although this method originates from loopy belief propagation, it does not suffer from the same computational complexity. AMP has been shown to be effective with a minimal number of measurements while being efficient in terms of computational complexity. Using a Bayesian approach, the main loop of AMP consists of iteratively updating the estimated mean x a and variance x v of the unknown signal until convergence, where • is the element-wise Hadamard product, (·) −1 is understood to be an element-wise reciprocal, (·) t is a time index, (·) * is the conjugate-transpose, and p in is a function based on the desired signal prior which returns both the mean and variance estimate of the unknown signal. We refer the reader to [18] for a detailed description of Bayesian AMP. The calibration and reconstruction phases employ Gaussian and binary priors, respectively [19,20]. From [11], we know that the transmission matrices of scattering mediums appear to be i.i.d. random matrices. Therefore, a Gaussian prior for the calibration phase is a reasonable choice. For the reconstruction phase, two binary priors have been investigated based on global (per-image) and local (per-pixel) sparsity, the details of which are explained in the next section. Generalized AMP (GAMP) [13] is an extension of AMP for arbitrary output channels, i.e. y = q(Hx + w). This adds an output function, p out , which is dependent on the stochastic description of q(·). In Eqs. (4)-(6), the terms (y − ω t ) • (v t−1 + σ 2 ) −1 and −(v t + σ 2 ) −1 indicate p out and p out , respectively, for a Gaussian output channel. One can easily modify these two terms in order to extend the framework to other channels. Following [13], for the phase retrieval problem we have where r 0 = I1(φ) I0(φ) , I 0 and I 1 are 0 th and 1 st order modified Bessel functions of first kind, respectively, and φ = 2y The convergence of both AMP and GAMP has been proved for zero-mean i.i.d. measurement matrices [21], however, they do not necessarily converge for generic matrices [22]. There have been some attempts to prevent divergence of AMP-based methods [14,23,24]. In [14], the authors show that a simple change in the main AMP loop may stabilize AMP significantly. They propose a sequential, or swept, random update of the AMP messages s t , r t , x t a and x t v , instead of their standard parallel calculation. By combining the swept update ordering and the phase retrieval output channel (8)- (9) in the AMP iteration (3)-(7), we create a phase retrieval version of SwAMP, denoted as prSAMP, which we describe in Algorithm 1.

EXPERIMENTAL RESULTS
To investigate the performance of the proposed imaging system, two binary datasets are constructed from the spatially-sparse MNIST handwritten dataset. The first, D 1 , consists of cropped digit images at a resolution of 20 × 20 pixels (N = 400). The second, D 2 , is constructed by rescaling the MNIST dataset to 32 × 32 pixels (N = 1024). Both D 1 and D 2 retain the original MNIST training/testing partition.
For the calibration step, training set D 1 is modified by randomly exchanging 5 × 5 blocks of pixels between digit images. Fig. 2 shows a few samples of these structured patterns. This structured randomization is done to reduce the effect of correlation between the DMD pixels. Additionally, to avoid the possibility of completely zero, or very sparse, lines in X, see Eq. (2), we introduce a fixed number of unstructured i.i.d. Bernoulli random binary patterns to the calibration training set.
The transmission matrix is then estimated from P = αN calibration images, of which the first N are Bernoulli random patterns, for the oversampling ratio α ≥ 1. At the receptor, M samples are randomly selected from a 100 × 100 region of the output image. Fig. 3 (left) shows the performance of the proposed calibration method for varying values of α. In lieu of ground-truth comparisons for transmission matrix estimation, we assess the calibration performance in terms of "dependence," the normalized cross-correlation without mean removal, y y , |Hx| Hx , between observed samples and the predicted output of known input patterns using the estimated transmission matrix. We measure dependence over 400 digits from the testing set of D 1 . We compare the level of achieved dependence between prSAMP and prVBEM [9], a mean-field variational Bayes phase retrieval technique we previously employed for the task of transmission matrix calibration in the context of light focusing.
After calibration, the direct imaging phase can start. As described in Section , the calibration and reconstruction steps are performed using the same prSAMP algorithm with different input priors. During calibration, we assume a complex Gaussian prior since the transmission matrix is modeled as i.i.d. random. However, for reconstruction, a binary prior is required, where , and ρ i indicates the probability of pixel i to be non-zero. We use two strategies to set this parameter. The first is a global approach which sets all ρ i uniformly to the input image sparsity level which we assume is known up to some tolerance. In the second local approach, we empirically calcu-late the per-pixel non-zero probability using the calibration training set, which is a fast off-line process. As the prior calculation must be repeated at each pixel for each sweep of the prSAMP algorithm, we select the simplest possible prior for the sake of computational efficiency. The interested reader may refer to [25] for a more sophisticated method of using learned priors for reconstruction tasks.
We next use the D 2 dataset to study the effectiveness of prSAMP post-calibration reconstruction. We first perform the calibration step to estimate M = N rows of the transmission matrix using α = 5, yielding an average calibration correlation of 97% over 1024 test digits. For reconstruction, we randomly choose 50 images from the test set, with five images for each digit. The correlation of prSAMP reconstructions to the true inputs, using global and local binary priors, are compared in Fig. 3 (right) as a function of the measurement rate M/N . Leveraging the extra information in the local prior provides an average 14.87% increase in reconstruction performance over the global prior. To visually assess the quality of recovered images, Fig. 4 provides one instance from each digit recovered at M = 700, with reconstructions using the local and the global priors. As expected, the local prior provides better subjective quality with fewer spurious isolated pixels.

CONCLUSION
In this study, a phase retrieval compressive imager has been proposed and experimentally evaluated using a simple optical setup. The imager has the potential of providing high resolution images in one shot. We solve the challenging problem of estimating a complex transmission matrix using binary patterns and we solve the phase retrieval problem via swept AMP. Finally, we show that we can estimate the transmission matrix accurately, allowing it to be used for compressive imaging. Further studies are necessary to provide faster calibration methods.