{'fraction_non_alphanumeric': 0.036843761472701526, 'fraction_numerical': 0.011512036502858342, 'mean_word_length': 4.30684664625661, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantitative phase and absorption contrast imaging

We present an algorithm for coherent diffractive imaging with phaseless measurements. It treats the forward model as a combination of coherent and incoherent waves. The algorithm reconstructs absorption and phase contrast that quantifies the attenuation and the refraction of the waves propagating through an object. It requires coherent or partially coherent illuminations, and several detectors to record the intensity of the distorted wave that passes through the object under inspection. The diversity of illuminations, obtained by putting masks between the source and the object, provides enough information for imaging. The computational cost of our algorithm is linear in the number of pixels of the image. Therefore, it is efficient for high-resolution imaging. Our algorithm guarantees exact recovery if the image is sparse for a given basis. Numerical experiments in the setting of phaseless diffraction imaging of sparse objects validate the efficiency and the precision of the suggested algorithm.


Introduction
The success of imaging an object with electromagnetic waves is determined as much by the development of new hardware elements as by the progress in designing efficient algorithms. These algorithms solve an inverse problem that, most often, is linear. This is the case, for example, of coherent imaging when the complex amplitudes of the waves are recorded, or the case of incoherent imaging when only the amplitudes squared, or intensities, are recorded. If, however, waves propagate coherently but only the intensities can be measured, then the imaging problem is nonlinear. This is the so called phase retrieval problem where we seek to reconstruct information about an object from intensity measurements of the wave traversing it.
In this context, we consider multiple measurements of a sparse complex-valued object f ∈ C K , which in our case is the complex refractive index. The measurements |b i | 2 ∈ R N of this phase object of size K are of the form where w i is a mask, and • denotes element-wise multiplication. In (1), T is a linear transformation representing the wave propagation operator, which is often approximated by the Fourier transform. Although the available measurements to recover f are intensities, the missing phases are still encoded in the recorded data because the waves propagate coherently and, hence, there is a fixed phase relationship between waves emerging from different points of f . Uniqueness of phase retrieval solution requires more than one diffraction pattern, so N > 1 in (1). The different diffraction patterns can be obtained using different illumination settings or, equivalently, different masks w i in (1), see [2,15,35,6]. In practice, masks are implemented using a spatial light modulator (SLM) or a digital micromirror device (DMD) [43,29,14,55]. We present a two-step algorithm for phase retrieval of sparse objects that mimics the forms in which waves propagate. In the first step, we assume that the intensities add incoherently, and we treat the coherent contribution to the data as a modeling error that is absorbed by a denoising algorithm. We use a Noise Collector [33] as our denoising algorithm.
In the case in which the object is interrogated incoherently, as in for example ghost imaging [21], the modeling error is zero, and this step produces exact support recovery if the data are not too noisy. Moreover, if the noise in the data is very strong, the support is always contained within the true support. For this very reason, only the strong absorbers are recovered in the first step if the object is interrogated coherently. Hence, a second step that takes into account the coherent contribution to the data in (1) is implemented. With this second step, the weak (semi-transparent) absorbers are recovered as well.
Once the support of both the strong and the weakly absorbing structures is known, a straightforward third step that solves the complete problem is implemented for a precise quantitative phase image restoration. This third step uses both the coherent and the incoherent contributions to the data.
This algorithm has applications in, for example, phase-contrast imaging that seeks to visualize semitransparent structures which are otherwise invisible in conventional absorption-contrast images. The idea behind this imaging modality is simple. When a wave is transmitted through an object it is not only absorbed, but also bent inducing phase changes. The phase changes themselves are invisible, so one needs a method to make them visible as brightness variations in the images. Our algorithm works in the far-field regime, where phase-contrast appears as the result of the free propagation of the waves that transforms phase variations due to the presence of an object into detectable intensity variations in the images [9,37,51]. This regime appears to be one of the most relevant and easy to implement for clinical and biomedical research purposes [27,22]. We refer the reader to [3] for an extensive discussion on phase-contrast imaging in clinical and biomedical applications, and to [39] for recent advances in optical phase imaging for investigating cells and tissues in biomedicine. Some methods for phase-contrast imaging that work well in the near-field regime but need specialized optical elements, that are hard to manufacture, are Zernike phase-contrast microscopy [53], analyzer-based methods [11,12], and grating-based methods [40,41].
The proposed algorithm also finds applications in fields such as exploration seismology, where phase errors in the data are often strong, so classical algorithms fail because they try to match the incorrect phase information. Hence, the measured phases are discarded and the inverse problem is formed for intensity-only data [20,1].
The paper is organized as follows. In Section 2, we review some of the most popular algorithms for phase retrieval, and we mention our main contributions. In Section 3, we introduce our model that produces the measurements for absorption and phase contrast imaging. In Section 4, we summarize the construction of the Noise Collector and some of its important properties. In Section 5, we present the algorithm that produces these images. Section 6 shows the numerical experiments. Section 7 summarizes our conclusions.

Related work & main contributions
There exists an important variety of algorithms that serve as approximate inverses of the problem stated in (1). Gerchberg and Saxton [19], and later Fienup [17], introduced a scheme based on iterative projections which is simple to implement and proved to be very flexible in practice. These methods are fast but, due to the absence of convexity of problem (1), they do not always converge to the true solution unless good prior information about the sought object or signal is available. We refer the readers to [30] for more details on these alternating projection techniques.
To guarantee convergence, Chai et al. [7] and Candés et al. [5] proposed a different approach that lifts problem (1) to a higher dimension to make it convex. The corresponding algorithm solves a low-rank matrix linear system using nuclear norm minimization. The algorithm converges to the true solution, even without prior information about the object and independently of its complexity, provided the recorded data is diverse enough.
However, when lifting the problem to a higher dimension, the size increases quadratically in the number of unknowns K, and the solution becomes infeasible for K large.
The authors in [50] follow a similar approach. They reformulate phase recovery as a MaxCut-like semidefinite program, and solve the resulting problem using a block coordinate descent algorithm similar to the Gerchberg and Saxton algorithm [19]. Although with this approach the problem can be solved more efficiently, its size also increases quadratically, so the approach is not suitable for large scale problems. Other interesting algorithms that require matrix lifting are, for example, [46,38].
Inspired by advances in compressed sensing, and in order to improve the efficiency of the existing algorithms, several works explore the idea of sparsity as prior information on the sought objects or signals, with sparsity levels s ≪ K. In [47], for example, the authors propose an efficient local search method that empirically recovers s-sparse objects, with s ∼ O(N 1/3 ). Here N stands for the dimension of the data. However, convergence to the correct solution is not guaranteed. In [24], the authors exploited the fact that the one-dimensional Fourier phase retrieval problem is equivalent to the turnpike reconstruction problem, where the objective is to reconstruct a set of vertices from the set of their pairwise distances. They showed that O(s 4 ) reconstruction times are possible when the set of vertices is s-sparse and s ∼ o(N 1/2 ). In [36], the authors generalize the one-dimensional algorithm [24] to higher dimensions and obtained dramatically faster O(s 2 ) reconstruction times of s-sparse objects, with s ∼ o(N 1/2 ). These algorithms do solve efficiently high-dimensional sparse problems. However, they only work for Fourier phase retrieval settings, i.e. T must be a Fourier transform in (1).
In this paper, we present an efficient and robust phase retrieval algorithm for O(N 1/2 )-sparse objects, whose computational cost grows linearly with the problem size K, so it is suitable for large scale problems. Although for ease of presentation it is introduced for Fourier measurements, it also works for general quadratic measurements without any modification. The algorithm vectorizes the matrix formulation in [7] and [5] to solve an optimization vector problem, and uses a Noise Collector [33] to reduce its dimensionality [34].
The huge reduction of dimensionality is carried out in two steps mimicking the forms in which waves travel. In this way, the algorithm merges coherent and incoherent imaging, as it considers the coherent and incoherent contributions to the data sequentially. With this approach we are able to image transparent or semi-transparent structures that are not visible in the commonly used absorption-based images.
As a by-product, our algorithm can also be used for phase retrieval with partially coherent observations without any modification and without loss of resolution. This might be important for applications where fully coherent sources are hard to produce as in, for example, X-ray phase-contrast imaging [25]. In these cases, only the signal-to-noise ratio (SNR) of the created images is affected. In the extreme case in which the observations are fully incoherent, the transparent structures cannot, of course, be visualized.

Model
When a wave propagates through an object, both its amplitude and phase are altered. Intuitively, the amplitude of the transmitted wave depends on the absorption of the wave, while its phase shift depends on the refraction.
Consider a planar object of finite support and thickness l illuminated perpendicularly by a monochromatic plane wave of wavenumber k = 2π/λ traveling in the z direction, see Figure 1. If the wavelength λ is small compared to its dimension l, the interaction of the wave with the object can be described in terms of integrals of the complex refractive index n(x, y), so its transmissivity is given by Then, the diffracted complex amplitude in the far-field, using the Fraunhofer approximation, is where L ≫ 1 is the distance between the object plane and the detector plane. We use prime symbols for the coordinates in the object plane to avoid confusion with the (x, y) coordinates in the detector plane. We use dimensionless coordinates by scaling the transverse coordinates (x, y) with √ λL. Fraunhofer diffraction occurs when L ≫ a 2 /λ, where a is a characteristic transverse length of the object. Physically, this means that the object is far enough from the sources, so the incident waves are effectively plane waves, i.e., the phase of the waves at each point on the object is the same.
Mathematically, the complex refractive index in (2), is expressed as It is the ratio of the wavenumber in the objectk(x ′ , y ′ ) and the wave number in the vacuum k and, hence, is a measure of how fast the waves travel through the object. The real and imaginary components δ and β determine the refraction and absorption effects of the interaction wave-matter, respectively. For a thin planar object, the phase shift and the absorption coefficient are well approximated by respectively. Here, the integration is taken over the extend of the object l along the direction of wave propagation z. Both (5) and (6) are proportional to the density of electrons at each point (x ′ , y ′ ) of the object and, hence, both allow for reconstructions of electron densities. However, these reconstructions are usually based on wave attenuation only, as phase shifts still remain harder to quantify. However, phase shifts also induce intensity variations on the images, and this effect can reveal important features of the object's structure that are not visible in the attenuation-based images. This is of great importance at high energy ranges where wave attenuation is small, or even negligible. Indeed, the refractive index decrement δ and the extinction coefficient β have a strong dependence upon the energy of the incident wave E, but they behave very differently as the energy increases. In particular, for X-rays, δ decreases approximately as 1/E 2 , while β approximately as 1/E 4 , so at these energies δ is between one and three orders of magnitude larger than β. For example, at E = 30keV , δ ≈ 2.5610 −7 and β ≈ 1.3610 −10 for water.

Weak phase-contrast
If the absorption is negligible so β(x ′ , y ′ ) ≪ δ(x ′ , y ′ ), then the phase shift (5) modulates the detected complex amplitude (3), with transmissivity given by If, in addition, (7) can be approximated as and the phase-contrast is weak. In this case, the diffracted complex amplitude is given by This is called the weak phase object approximation. Integration of the first term gives a delta function (for an infinite aperture) that represents the direct wave that goes through the object without interaction. We do not consider it here, as it is usually blocked by a beam stop [49]. Thus, only the scattered component up to a constant that is set here equal to one, is considered for imaging. Assume that, for imaging purposes, the inspected object is discretized using a grid of K pixels, so its phase φ(x ′ , y ′ ) is well approximated by the size K vector According to (10), if the object is illuminated by a monochromatic coherent plane wave, the complex diffraction pattern b measured at an N -pixel detector is given by the size N discrete Fourier transform of the object's phase given by . . , N , and k = 1, . . . , K. An important special case is when N = K, so the transformation is the classical discrete Fourier transform. If N < K, the case considered here, the Fourier transform is said to be undersampled.
If a known set of N spatially structured patterns or masks illuminate the object, then several complex diffraction patterns are available for imaging. Here, W i is the diagonal matrix W i = diag(w i1 , . . . , w iK ) corresponding to the i-th illumination pattern. With this notation, the s-th component of the data vectors represents the complex field at the s-th detector when i-th illumination impinges on the object. If the detectors can only measure intensities, then the problem is to recover (11) from measurements of the type i = 1, . . . , N, s = 1, . . . , N . Note that, in this case, the unknown vector φ is real. The left side represents the intensity received at the s-th detector when the object is illuminated with the i-th illumination pattern w i = [w i1 , . . . , w iK ] ⊺ . The first term in the right side of this expression is the incoherent contribution to this intensity, and the second term is the coherent contribution. The coherent contribution is characterized by a fixed phase relationship between the waves emerging from different points of the object, while in the incoherent contribution this phase relationship does not exist. The coherent contribution produces the interferences used in coherent imaging to determine the object's structure, and it is the term that makes phase-retrieval non-linear.

Incoherent imaging
According to the previous discussion, if the phases of the waves immediately exiting the object are randomized either by a diffuser or by the medium's inhomogeneities, so the wavefront is totally scrambled, then the second term in (13) is negligible and the inverse problem is linear in the phase-shifts (squared) |φ k | 2 . In this case, the intensities received at all the detectors are equal; the data vector does not depend on s. Thus, to improve the SNR, the collected intensities can be averaged and we can form the linear system Here, each row of the matrix W incoh contains the intensities upon the pixelated object produced by the N masks w i . The linear system W incoh χ = d avg (15) can be easily solved by means of any ℓ 1 -or ℓ 2 -method depending on the number of masks N used for imaging and the sparsity of the vector χ, which is in this case real.

Non-weak phase-contrast
If the absorption is negligible but the phase contrast is not weak but strong, then (7) cannot be approximated by (8) and we should image the transmissivity vector which is now complex. Then, (13) becomes In (17), the first term in the right side is the total intensity c upon the object, which is known, and that does not provide any information for imaging. All the information is contained in the interferences of the scattered wave represented by the second term. At least two approaches are possible to solve the nonlinear problem (17). One is to solve it iteratively for the K complex unknowns t i as in [19,17]. This phase retrieval algorithms work well in practice if good prior information about the object of interest is known; otherwise, convergence to the true solution is not guaranteed. The second option is to reformulate (17) as a linear matrix problem for the K 2 unknowns t i t * j , i, j = 1, . . . , K, as in [7,5], and solve the resulting problem by using nuclear norm minimization. This option guarantees convergence to the unique solution without any prior information about the sought object, but the number of unknowns grows quadratically with the number of unknowns K and, thus, the solution becomes unfeasible for high-resolution images with large K.
Alternatively, one can define the cross-correlation vector of t excluding the real terms |t i | 2 = 1, and solve the huge linear system given T − c represents the N intensities recorded at the detector s minus the total intensity upon the object c. In (19), is a huge matrix of size N N × K(K − 1), where W s,coh is defined in (21) with c si,lm = F sl F * sm w il w * im . This  c sN,12 c sN,13 . . . c sN,1K c sN,21 c sN,23 . . . c sN,2K c sN,31 . . .
matrix models the coherent component of the intensities received at the s-th detector corresponding to the N illumination patterns w i used to filter the object. If the resolution is low, so the number of pixels K is small, one could easily solve (19) using an appropriate solver. However, as in the approach suggested in [7,5], the size of the problem increases quadratically with K, so the search of a solution rapidly becomes prohibitive for large values of K, as well.

Absorption and phase contrast
Absorption can be easily taken into account by using (5) and (6) in (2). If absorption is not negligible, then the transmissivity is given by This is the most general case in which absorption and refraction effects are mixed in the images. In these cases, the transmissivity vector is given by the complex vector with amplitudes different than 1. Then, the problem is to find (23) from measurements of the form i = 1, . . . , N, s = 1, . . . , N . This problem has, of course, the same form as before, but now the incoherent contribution to the intensity is modulated by the absorption. This makes attenuation-based imaging possible with incoherent sources.
The main purpose of this paper is to propose a dimension reduction technique that allows to find the solution of (24) efficiently for large K. This algorithm is presented in Section 5. The set of equations (24) can be written in matrix form as with W incoh and W coh defined in (14) and (20), respectively. Solution of (25) provides complete images, with distributions of both the real and imaginary parts of the complex refractive index (4). However, the bottleneck is still the size of the problem, which is enormous if one wants to form high resolution images. An image with only 1000 × 1000 pixels, amounts to solving a linear system with 10 12 unknowns! In Section 5, we describe the proposed algorithm that makes possible to find the desired solution in polynomial-time if the object to be image is sparse in some appropriate basis, but first we summarize the main properties and the construction of the Noise Collector which is an essential denoising tool in this approach. For simplicity, we will assume that the sought object is sparse in the real space, meaning that only a few pixels in the object plane absorb or bend the waves. We note, though, that many images are naturally compressible by using appropriate sparsifying transforms such as wavelets, or other dictionaries that are directly adapted to the data [42].

The noise collector
The Noise Collector [33] is a denoising algorithm to find the vector χ ∈ C K in from highly incomplete measurement data d = d 0 + e ∈ C N corrupted by additive noise e ∈ C N , where 1 ≪ N < K. Here, A is a general measurement matrix of size N × K, whose columns have unit length. The main result in [33] ensures that we can recover the support of χ by looking at the support of χ τ found as (χ τ , η τ ) = arg min χ,η (τ χ ℓ1 + η ℓ1 ) , subject to Aχ + Cη = d.
Here, τ is an O(1) no-phantom weight, and C ∈ C N ×Σ is a Noise Collector matrix, with Σ = N β , for β > 1 (typically close to one). If the noise e is Gaussian, then the columns of C can be chosen independently and at random on the unit sphere S N −1 . The weight τ > 1 is chosen so it is expensive to approximate e with the columns of T , but it cannot be taken too large because then we lose the signal χ that gets absorbed by the Noise Collector as well. For practical purposes, τ is chosen as the minimal value for which χ = 0 when the data is pure noise, i.e., when d 0 = 0. The key property is that the optimal value of τ does not depend on the level of noise and, therefore, it is chosen in advance, before the Noise Collector is used for a specific task.
It can be shown that if the matrix A is incoherent enough, so its columns are not almost parallel, the minimizer in (27) has no false positives for any level of noise, with probability that tends to one as the dimension of the data increases to infinity. To find the solution of (27), we define the function This strategy finds the minimum in (27) exactly for all values of the regularization parameter λ. Thus, the method is fully automated, meaning that it has no tuning parameters. To determine the exact extremum in (29), we use the iterative soft thresholding algorithm GeLMA [32] that works as follows. Pick a value for the no-phantom weight τ ; for optimal results calibrate τ to be the smallest value for which χ = 0 when the algorithm is fed with pure noise. In our numerical experiments we use τ = 2. Next, pick a value for the regularization parameter, for example λ = 1, and choose step sizes ∆t 1 < 2/ [A | C] 2 and ∆t 2 < λ/ A 1 . Set χ 0 = 0, η 0 = 0, z 0 = 0, and iterate for k 0: where S r (y i ) = sign(y i ) max{0, |y i | − r}. Terminate the iterations when the distance χ k+1 − χ k between two consecutive iterates is below a given tolerance. For more details on the theory and properties of the Noise Collector, we refer the reader to [33,34].

Dimension reduction
Instead of solving the linear system with K 2 variables (25) for non-weak phase objects, with or without absorption, we propose to reduce its dimensionality by constructing a linear problem for only O(K) significant unknowns, and absorb the error corresponding to the contribution of the unmodeled unknowns by using a Noise Collector. Mathematically, we propose to solve the linear system where A is a matrix with O(K) subsampled columns of the matrix [W incoh |W coh ], each column of size N N . In this formulation, χ is a sparse vector that represents the object, η is an unwanted vector with no physical meaning that absorbs the noise, and C is a Noise Collector matrix with N β columns drawn independently and at random on the unit sphere (we use β = 1.5 in our simulations). The algorithm has three steps.
(1) In the first step, we seek the strong absorbing objects. We set A = W incoh so W coh = 0, and solve (31) for χ = [|t 1 | 2 , |t 2 | 2 , . . . , |t k | 2 ] T . The term Cη in (31), where C is a small matrix with N β columns, absorbs the coherent contributions to the intensities that are treated in this step as noise. Since the model we solve is not exact, only the strong absorbing objects are detected.
(2) In the second step, we seek the non-absorbing objects. Since these objects are almost transparent and do not have a significant impact on the recorded intensities, we look for their phases, encoded in the vector χ cross defined in (18). To this end, we first subtract the incoherent contribution to the recorded intensities due to the detected objects in the first step. Then, we set W incoh = 0 and A = (W coh ) sub , where (W coh ) sub is a small subsampled matrix of the huge matrix W coh . It only contains the m (K − 1) columns that correspond to the interactions between the m detected objects in the first step and the other pixels in the image. Since we are not modeling the incoherent contributions of the remaining (K − m) pixels, the system we solve is not exact neither. Hence, we also use a Noise Collector matrix C with N β columns to absorb the noise.
(3) The third step is optional. It is used to obtain more precise quantitative images. Once the strong and weakly absorbing objects are found, we solve the full problem (25) restricted to the recovered support. This is now a small problem that can be solved using an ℓ 2 minimization method that gives very accurate results.
We stress that, for this dimension reduction strategy, it is necessary that the unknown object can be represented as a sparse vector. Otherwise, the (modeling) errors are too big to be absorbed.

Numerical experiments
The simulations shown here illustrate the potential of the proposed algorithm. In this work, a thin object is illuminated with quasi-monochromatic coherent or with partially incoherent sources. The schematic for the imaging setup is shown in Figure 1. The units of the problem are given with respect to the central wavelength λ 0 , as this is the important parameter in the simulations, and all other length scales will be referred to it. The main assumption is the sparsity of the object. We show here reconstructions of small point-like structures. As for other compressed sensing based algorithms the methodology can be used for more complex imaging scenes, as long as there is a sparsifying transform that allows a sparse representation of the scene. In some applications, they use an off-the-shelf transform like the Fourier, Hadamard, wavelet, or curvelet, while in others they find a new transform using dictionary learning. Note that the sparsity is unknown and can scale as √ N N with N N the number of data (intensities) measured. The sources are located on a two dimensional array of size 16000λ 0 × 16000λ 0 , at a distance of 16000λ 0 from the object. There are 21 × 21 point sources, evenly distributed on the source plane. These are used to create N illumination patterns on the imaging plane. The illuminations are created by assigning random amplitudes and phases to each one of the point sources. The resulting illumination patterns w i , for i = 1, . . . , N , are assumed to be known on the imaging plane.
We use 5 × 5 receivers 16000λ 0 downstream to collect the data corresponding to N = 300 different illumination patterns. The receivers are evenly distributed on the measurement plane that is parallel to the object and to the transmitting array (see Figure 1). The receiving array also has an aperture of 16000λ 0 × 16000λ 0 . Wave propagation is modelled using the 3D wave Green's function, with |x − y| the distance between points x and y and k 0 = 2π/λ 0 the wavenumber. In order to form an image, the imaging plane is discretized using 31 × 31 pixels, with pixel size equal to half a wavelength in both directions., i.e. a square of size λ 0 /2 × λ 0 /2. Although we consider a typical transmission setup, same results are obtained for reflection or for more complicated sources and receivers layouts.
We consider absorbing and non-absorbing objects that change the phases of the waves that go through them. We form images of ρ(x ′ , y ′ ) = t(x ′ , y ′ ) − 1. The illumination by the direct wave may be absorbed by adding a column to the Fourier transform matrix F . In our numerical experiments, ρ k = e iπ if there is a strong absorbing object at pixel k, and ρ k = 0.1 e iπ/2 if there is a non-absorbing object. If only intensities are recorded, these non-absorbing objects are very hard to image because the waves that go through them change only slightly.
In the first step of the algorithm, we seek to reconstruct the strong absorbing objects. As explained in section 5, we set A = W incoh , and solve Aχ + Cη = d , The noise collector term Cη absorbs the contributions of ρ * i ρ j for i = j to the data which are neglected in this step and treated as noise. Looking at we observe that the first term is independent of the receiver s. Therefore in this first step we use the total intensity as data (the sum over all the receivers).
To gain some intuition in how the algorithm works, consider that we want to image m strong absorbers |ρ i | = O(1), i = 1, . . . , m and n weak absorbers |ρ j | = O(ε), j = 1, . . . , n. The n weak absorbers have phase contrast that we wish to reconstruct. During the first step of the algorithm we only recover the strong absorbers |ρ i | 2 , i = 1, . . . , m, because the contribution from the weak ones |ρ j | 2 = O(ε 2 ) j = 1, . . . , n, is lost in the noise. true |ρ k | 2  This is illustrated in Figure 2, where we have considered four scatterers, two of them are strong absorbers (m = 2) shown with red squares and the other two are weak (n = 2) shown with white crosses. The ratio between the weak and the strong absorbing objects in |ρ i | 2 is, in this case, of the order of 1/100, making the detection of the weak absorbers very difficult. The left column of Figure 2 are the results for noise-free data, and the right column is the results for data with SNR= 30 dB. In both cases, the locations of the strong absorbing objects are recovered exactly. Moreover, their amplitudes |ρ i | 2 are recovered with a quite good accuracy. In the first and second rows of Figure 2 we display |ρ| 2 as a two-dimensional image, while in the third row we plot |ρ| 2 as a vector. In this third row, we plot the exact |ρ| 2 vector with green circles, and the recovered one with red stars. The black stars are the non-physical unknown η introduced in the algorithm to absorb the contribution to the data due to the cross-terms ρ * i ρ j , with i = j. In both cases, with or without noise in the data, the weakly absorbing objects are not recovered because the neglected contribution to the data of the cross-terms ρ * i ρ j , with i = j , is larger than the contribution of the weakly absorbing objects (these are the second term which is O(1) and third term which is O(ε) in (32)).
To find the weak absorbers we apply the second step of the algorithm. To this end, we first remove from the data the O(1) contributions m k=1 |w ik | 2 |ρ k | 2 from the m strong absorbers already found in step 1. For our model problem with m-strong and n-weak absorbers what remains is Then, for every pixel i = 1, . . . , m detected during the first step we seek for its interactions ρ * i ρ j with all the other K − 1 pixels in the object plane, j = 1, . . . , K, j = i. These are the O(1) and O(ε) contributions to the data (second and third term in (32)). In this case A = (W coh ) sub , where (W coh ) sub contains the m (K − 1) columns that correspond to the interactions between the m detected objects in the first step and all the other pixels in the image. Since we are neglecting the O(ε 2 ) contributions, the system is not exact.
For the example shown in Figure 2, we found m = 2 strong absorbers, so we have 2K − 2 unknowns. The results of this second step are shown in Figure 3. In top row we display the unknown recovered by considering the interactions with the first scatterer, while in the center row we display the unknown recovered by considering the interactions with the second scatterer. In the bottom row we plot the unknown t * i t j ; the green circles represent the true solution, the red stars the unknown recovered by ℓ 1 minimization, and the black stars the non-physical part of the unknown corresponding to the Noise Collector. This second step finds the weak absorbers with amplitudes and phases recovered with good accuracy; see the results in Figures 4 and 5 that summarize these results. In this example we do not need to apply the third step of the algorithm as both amplitudes and phases of the unknown are recovered correctly after the first two steps. Next, we consider in Figure 6 a more challenging example with ten absorbing objects, with only a strong one (m = 1 and n = 9), and SNR = 30dB. As in the previous example, only the strong absorber is found with the first step of the algorithm in the absorption-based image. In the second step, as expected, we also find all the weak ones in the phase contrast image. The last row of Figure 6 shows the reconstructed amplitudes (left) and phases (right) of the weakly absorbing objects plotted as vectors. The exact values are represented with green circles and the reconstructed values with red stars.
We observe from the results in Figure 6 that the phases of the weakly absorbing objects are recovered with much more accuracy than their amplitudes. This is because the error induced by the neglected terms in the second step increases with the number of objects. Indeed, in the second step we only account for the interactions between the strong and the weak absorbers (the O(ε) term in (32) that corresponds to n = 9 contributions here), but we neglect all the interactions between the weak absorbers (this is the last O(ε 2 ) term in (32) which corresponds to n 2 = 81 contributions here). This example is therefore more challenging because the modelling error increases quadratically with the number of weak absorbers. Better results can be obtained by considering, in a third step, the full problem (25) restricted to the recovered support. This third step allows us to recover the unknown X supp = ρ supp ρ * supp accurately using an ℓ 2 minimization method. This is because the locations of all the absorbers both weak and strong are recovered exactly after the first two steps of the algorithm. In the third step the contributions to data from all terms in X are taken into account. Figure 7 shows that this provides a great accuracy in the recovered values of both the amplitudes and the phases. Figure 8, that shows the true and recovered phase distributions, illustrates the potential of the proposed imaging method for imaging both strong and weak absorbers with intensity-only measurements.
true |ρ k | 2 Finally, we present results for the case in which the illuminations are partially coherent. This is an interesting case because in some applications fully coherent illuminations are very hard to obtain. We use the following model to generate the data with 0 ≤ α coh ≤ 1. If α coh = 1, the sources are fully coherent, and if α coh = 0 they are fully incoherent. The parameter α coh that models our uncertainty about the coherence of the illumination is not known when we seek to reconstruct the unknown absorbers.  Figure 9 shows the results when the illumination used for imaging is partially coherent; α coh = 0.5 in this numerical experiment. The left and right columns show the outputs of the first and second steps of the algorithm, respectively. Because the illumination is partially coherent and, thus, the modeling error in the first step is smaller, we observe that the first step recovers the amplitudes of the strong absorbers with great accuracy. The second step is still able to recover the locations of the weak absorbers exactly. However, as expected, we observe that there is a SNR issue, and that if α coh decreases below a certain threshold, we would not be able to image them. This threshold depends on the transparency of these objects, their number, and the noise in the data. For the numerical experiment shown here, the phases of all the absorbers are recovered with the same precision as in the previous experiment (results not shown).

Conclusions
In this paper we have presented a two steps algorithm for phase retrieval for sparse objects. This algorithm is very efficient because its cost is linear in the number of pixels of the image and, thus, it can be employed for high resolution imaging. It guarantees exact recovery if the image is sparse with respect to a given basis, and it can be used, without any modification, when the illumination is partially coherent. In these cases, only the SNR of the created images is affected. Although for ease of presentation this algorithm is introduced for Fourier measurements, it also works for general quadratic measurements without any modification. With this algorithm we are able to image transparent or semi-transparent structures that are not visible in the common used absorption-based images.