Multiwavelength Absolute Phase Retrieval from Noisy Diffractive Patterns: Wavelength Multiplexing Algorithm

We study the problem of multiwavelength absolute phase retrieval from noisy diffraction patterns. The system is lensless with multiwavelength coherent input light beams and random phase masks applied for wavefront modulation. The light beams are formed by light sources radiating all wavelengths simultaneously. A sensor equipped by a Color Filter Array (CFA) is used for spectral measurement registration. The developed algorithm targeted on optimal phase retrieval from noisy observations is based on maximum likelihood technique. The algorithm is specified for Poissonian and Gaussian noise distributions. One of the key elements of the algorithm is an original sparse modeling of the multiwavelength complex-valued wavefronts based on the complex-domain block-matching 3D filtering. Presented numerical experiments are restricted to noisy Poissonian observations. They demonstrate that the developed algorithm leads to effective solutions explicitly using the sparsity for noise suppression and enabling accurate reconstruction of absolute phase of high-dynamic range.


Introduction
Reconstruction of a wavefront phase is on demand in holography [1,2], interferometry [3], holographic tomography [4,5] and is of general interest due to important practical applications in such areas as biomedicine (for instance for studies of living cells [6,7]), 3D-imaging [8], micromachining [9], etc. There are two different approaches for extraction of spatial phase information. The first, which is associated with the holography, implies recording interference patterns of object and reference waves, while the second approach involves only sets of coherent diffraction patterns of object without reference waves. In these intensity measurements the phase information of the registered wavefields is eliminated, and then reconstructed by means of numerical iterative algorithms. The latter problem is well known as the phase retrieval, where "phase" emphasizes that the missing phase defines the principal difficulty of the problem and its difference with respect to holography.
The first phase retrieval algorithm was introduced by R. Gerchberg and W. Saxton [10]. They used squared roots of intensities measured in image plane to replace amplitudes calculated in iterations. Important update was developed by J. Fienup [11], who introduced constrains for the iterative process by the amplitude distribution in Fourier plane, namely, a nonnegativity of the reconstructed image and energy concentration in a certain area. A dataset of recorded intensity patterns which is characterized by diversity of one or several parameters of the optical setup is a crucial moment of the phase retrieval problems. Defocusing and phase modulation are known as effective instruments in order to achieve the measurement diversity sufficient for reliable phase reconstruction (e.g., [12][13][14][15]). Various types of diffraction optical elements can be used for design of phase imaging systems [16]. The review and analysis of this type of the systems as well as further developments of the algorithms can be found in [17].
There are various optimization formulations of conventional Gerchberg-Saxton (GS) algorithm [18][19][20]. Contrary to the intuitively clear heuristic of the GS algorithm, the variational approaches have a stronger mathematical background including image formation modeling, formulation of the objective function (criterion) and finally going to numerical techniques solving corresponding optimization problems. The recent overview [21] is concentrated on the algorithms for the phase retrieval with the Fourier transform measurements and applications in optics. The constrains sufficient for uniqueness of the solution are presented in detail. The fundamental progress for the methods based on the convex matrix optimization is announced in [22], where the novel algorithm "Sketchy Decisions" is presented. This algorithm is effective for high dimensions typical for the phase-lifting methods and is supported by the mathematical analysis and the convergence proof. The Wirtinger flow (WF) and truncated Wirtinger flow (TWF) algorithms presented in [23,24] are iterative complex-domain gradient descent methods applied to Gaussian and Poissonian likelihood criteria, respectively.
The techniques based on the proximity operators developed in [25,26] provide a regularized optimization of the maximum likelihood criteria also for both Gaussian and Poissonian observations. The sparsity based techniques is a hot topic in phase retrieval (e.g., [27]). A sparse dictionary learning for phase retrieval is demonstrated in [28]. The transform domain phase and amplitude sparsity for phase imaging is developed in [29][30][31]. For the phase retrieval these techniques are applied in [32][33][34].
In recent years, there has been an increase in demand for multispectral imaging. Multispectral information by recording object waves with multiple wavelengths that are irradiated from multiwavelength light sources helps to analyze and recognize objects, to clarify color and tissue distributions and dramatically improves the quality of imaging. Multiwavelength digital holography has an enhanced ability for 3D wide-range shape measurements by using multiwavelength phase unwrapping, due to the recording of multiwavelength quantitative phase information. The multiwavelength is an effective instrument enabling good phase and data diversity in phase retrieval.
The multiwavelength absolute phase imaging is much less studied as compared with the standard single-wavelength formulation. These works by the principle of measurements can be separated in two groups. In interferometry and holography scenarios the reference beams are applied to reveal the phase information, e.g., [35][36][37][38][39][40]. After that, the phase unwrapping algorithms for 2D images are applied. Recently, the phase unwrapping with simultaneous processing of multiwavelength complex-exponent observations have been developed based on the maximum likelihood techniques [41,42].
Another group of the multiwavelength absolute phase imaging techniques uses as measurements amplitudes or intensities of coherent diffraction patterns of an object without reference waves. In these intensity measurements the phase information is missing. These formulations are from the class of the multiwavelength phase retrieval problems, e.g., [43][44][45][46][47][48].
The Chinese Remainder Theorems, existing in many modifications, provide methods for absolute phase reconstruction from multiwavelength data for noiseless data [49]. The reformulation of these mathematically strong approaches for more practical scenarios with noisy data and for robust estimation leads to the techniques similar to various forms of maximum likelihood [50,51].
In this paper we study the problem of multiwavelength absolute phase retrieval from noisy object diffraction patterns with no reference beams. The system is lensless with multiwavelength coherent input light beams and random phase masks applied for wavefront modulation. The light beams are formed by light sources radiating all wavelengths simultaneously and a Complementary Metal Oxide Semiconductor (CMOS) sensor equipped by a Color Filter Array (CFA) is used for spectral measurement registration. The developed algorithm is based on maximum likelihood optimization, in this way targeted on optimal phase retrieval from noisy observations. The algorithm valid for both Poissonian and Gaussian noise distributions is presented.
One of the key elements of the algorithm is an original sparse modeling of the multiwavelength complex-valued images based on the complex-domain block-matching 3D filtering. Numerical experiments demonstrate that the developed algorithm lead to the effective solutions explicitly using the sparsity for noise suppression and enabling accurate reconstruction of absolute phase for high-dynamic range of variations.
In this work, we use the approach which can be treated as a development of the Sparse Phase and Amplitude Reconstruction (SPAR) algorithms [32,33] and the maximum likelihood absolute phase reconstruction for multiwavelength observations [42]. This approach has been exploited in our paper [52] for absolute phase retrieval from multifrequency in the different setup with observations obtained for multiple single wavelength experiments. Remind, that in this paper the experiments are multiwavelength and phase information for each wavelength is retrieved by the developed wavelength multiplexing algorithm.

Multiwavelength Object and Image Modeling
We present the problem of the absolute phase retrieval in the framework of thickness (profile) measurement for a transparent object. A coherent laser beam impinges on surface of the object of interest. A propagation of the coherent beam through the object results in a phase delay between the original input and the output beams. In general, this phase delay is spatially varying with variations proportional to variations of the thickness of the object. If the phase delay and the wavelength of the laser beam are known the thickness of the object can be calculated. However, only wrapped phase delays within the interval of 2π radians can be measured in experiments with a single wavelength. It follows, that only very thin objects can be examined in this way. Use of multiwavelength beams allows to overcome this restriction and to reconstruct absolute phases, and in this way to produce surface measurements for much thicker objects. If the thickness is known the absolute phases for all wavelengths can be defined. Thus, the approach can be focused in two a bit different directions: measurement of the absolute phases for multiple wavelengths and measurement of the object surface. This task is similar to the problem of surface contouring of a reflective object [43]. The model developed below can be applied to both types of objects. Our aim to consider in details a more interesting case of the transparent object with dispersion of the refractive index, because ignoring this dependence can cause significant distortions in the profile of the object being reconstructed [53].
Let h o (x), x ∈ R 2 , be a thickness variation of a transparent object to be reconstructed. A coherent multiwavelength light beam generates corresponding complex-valued wavefronts In the lensless phase retrieval, the problem at hand is a reconstruction of the profile h o (x) from the noisy observations of diffractive patterns, registered at some distance from the object. For the high magnitude variation of h o the corresponding absolute phases ϕ o,λ take values beyond the basic phase interval [−π, π], then only wrapped phases can be obtained from u o,λ .
We demonstrate that the multiwavelength setup allows to reconstruct the profile of the objects which height variations are significantly greater than the used wavelengths as well as to reconstruct the corresponding multiwavelength absolute phases.
The link with the previous notation is obvious: µ λ = λ (n λ − 1) λ(n λ − 1) and Here λ ∈ Λ is a reference wavelength and ϕ o is an absolute phase corresponding to this wavelength.
We are interested in two-dimensional imaging and assume that amplitude, phase and other variables, in what follows, are functions of the argument x given on a regular N × N grid, X ⊂ Z 2 .
The parameter µ λ establishes a link between the absolute phase ϕ o and the wrapped phase ψ o,λ of u o,λ which can be measured at the λ-channel. The wrapped phase is related with the true absolute phase, ϕ o , as µ λ ϕ o = ψ o,λ + 2πk λ , where k λ is an integer, ψ o,λ ∈ [−π, π). The link between the absolute and wrapped phase conventionally is installed by the wrapping operator as follows: (3) W (·) is the wrapping operator, which decomposes the absolute phase µ λ ϕ o into two parts: the fractional part ψ o,λ and the integer part defined as 2πk λ .
The image and observation modeling is defined as where u s,λ is a wavefront propagated to the sensor plane, P s,λ,d {} is an image (diffraction pattern) formation operator, i.e., the propagation operator from the object to the sensor plane, y s,λ is the intensity of the wavefront at the sensor plane and z s,λ are noisy observations of the intensity as defined by the generator G{} of the random variables corresponding to y s,λ , S is a number of experiments, d is a propagation distance. The considered multiwavelength phase retrieval problem consists in reconstruction of ϕ o and b o,λ from the observations z s,λ provided that µ λ , P s,λ,d {} and G{} are known. If the phases µ λ ϕ o ∈ [π, −π), the problem becomes nearly trivial as estimates of µ λ ϕ o and b o,λ can be found processing data separately for each λ by any phase retrieval algorithm applicable for complex-valued object, in particular by the SPAR algorithm [54]. It gives for us the estimates of b o,λ and µ λ ϕ o directly. These estimates of µ λ ϕ o (x) are biased as the intensities y s,λ are insensitive with respect to invariant additive errors in ϕ o . We will represent these estimates in the form The problem becomes nontrivial, much more interesting and challenging when the object phases µ λ ϕ o go beyond the range [π, −π) and the phase unwrapping is embedded in the multiwavelength phase retrieval. This paper is focused on the reconstruction of the absolute phase ϕ o . The light beam is formed by light sources radiating all wavelengths simultaneously and the standard RGB CMOS sensor with CFA is used for measurement registration.
Then, the measurement model Equations (4)-(6) is modified to the following form: where Ω c , c ∈ (R, G, B), are the Bayer sensor areas with the corresponding RGB color filters. The index c of the measured z s,c in Equation (8) is addressed to the sensor pixels with the corresponding color filters.
These Ω c produce a segmentation of the sensor area Ω, such that ∪ c Ω c = Ω, Ω c ∩ Ω c = ∅. In a bit different terms, x ∈ Ω c , c ∈ (R, G, B), defines a subsampling of the sensor output forming images of different colors. The sets Ω c define the structure of CFA (or mosaic filter) used in the sensor, i.e., the coordinates of the filters of the different colors.
The ideal color filters would have equal sensitivity to each of the spectral components and provide a perfect separation of the three RGB sources of the corresponding wavelengths. In reality, the spectral characteristics of the sensor must be measured [55] and taken into account. In particular, the cross-talk between the color channels exist on both sides in sources and sensor color filters [38]. As a model of this cross-talk the following linear convolution equation can be used: where u s,λ are the source wavefronts and u s,c are the sensor wavefronts as they are registered by the sensor pixels.
In general, the number of the spectral components n Λ can be different from the number of the color filters in the sensor.
Here w c,λ are the cross-talk spectral weights formalizing the sensitivity of (R, G, B)-filter outputs to the wavelength of the sources with the index λ. Despite the importance of the cross-talk effects, in this paper for simplicity of presentation we will assume that the color filters of CFA are nearly ideal and the cross-talk problem can be omitted.

Noisy Observations
The measurement process in optics amounts to count the photons hitting the sensor's elements and is well modeled by independent Poisson random variables (e.g., [26,56,57]). In many applications in biology and medicine the radiation (laser, X-ray, etc.) can be damaging for a specimen. Then, the dose (energy) of radiation can be restricted by a lower exposure time or by use a lower power radiation source, say up to a few or less numbers of photons per pixel of sensor what leads to heavily noisy registered measurements. Imaging from these observations, in particular, phase imaging is called photon-limited.
The probability that a random Poissonian variable z s,λ (x) of the mean value y s,λ (x) takes a given non-negative integer value k, is given by where y s,λ (x) is the intensity of the wavefront at the pixel x.
Recall that the mean and the variance of Poisson random variable z s,λ (x) are equal and are given by y s,λ (x)χ, i.e., E{z s,λ (x)} = var{z s,λ (x)} = y s,λ (x)χ, here E{} is a mathematical expectation, var{} is a variance. Defining the observation Signal-to-Noise Ratio (SNR) as the ratio between the square of the mean and the variance of z s,λ (x), we have SNR = E 2 {z s,λ (x)}/var{z s,λ (x)} = y s,λ (x)χ. Thus, the relative noisiness of observations becomes stronger as χ → 0 (SNR → 0) and approaches zero when χ → ∞ (SNR → ∞). The latter case corresponds to the noiseless scenario: z s,λ (x)/χ → y s,λ (x) with the probability equal to 1.
The parameter χ > 0 in Equation (10) is a scaling factor defining a proportion between the intensity of the observations with respect to the intensity of the input wavefront. This parameter is of importance as it controls a level of the noise in observations. Physically it can be interpreted as an exposure time and as the sensitivity of the sensor with respect to the input radiation.
In order to make the noise more understandable the noise level can be characterized by the estimates of SNR SNR λ = 10 log 10 and by the mean value of photons per pixel: Here N sensor is a number of sensor pixels. Smaller values of χ lead to smaller SNR and N photon , i.e., to noisier observations z s,λ (x).

Algorithm Development
We consider the problem of wavefront reconstruction as an estimation of u o ∈ C n from noisy observations {z s,λ }. This problem is rather challenging mainly due the periodic nature of the likelihood function with respect to the phase ϕ o and the non-linearity of the observation model. Provided a stochastic noise model with independent samples, the maximum likelihood leads to the basic criterion function where l(z, |u| 2 ) denotes the minus log-likelihood of a candidate solution for u o given through the observed true intensities |u| 2 and noisy outcome z.
Including the image formation model Equations (4), (7) and (8) and amplitude/phase modeling of the object Equation (1) by the penalty quadratic norms of residual we arrive to the following extended criterion: 1 where γ 1 , γ 2 > 0 are parameters, and || · || 2 2 is the Hadamard norm. Remind, that phase retrieval algorithms may reconstruct the phase µ λ ϕ o only within invariant additive errors. These unknown phase shifts are modelled by δ λ as disturbance parameters.
We say that the object complex exponents of different frequencies are in-phase (or synchronized) if δ λ = 0, λ ∈ Λ and de-phased if some of δ λ = 0. In the synchronized case the complex exponents of the object become equal to 1 provided ϕ o ≡ 0. For proper estimation of ϕ o these phase shifts δ λ should be estimated and compensated.
The first summand in Equation (14) is a version of Equation (13) modified for the multiwavelength CFA. The summation on x ∈ Ω c is produced independently for each c, what separates the observations of different wavelengths. The sets Ω c define the structure of the CFA (or filter mosaic) used in the sensor, i.e., the coordinates of the filters of the different wavelengths.
In the second summand u s,c is approximated by P s,c u o,c . The phase retrieval is a reconstruction of ϕ o and the amplitude of u o,λ from the observations z s,c . The algorithm is designed using minimization of L with respect to u s,c , u o,λ , b o , ϕ and the phase-shifts δ λ .
The following solutions of the optimization problem are used in the developed algorithm.

1.
Minimization of L with respect to u s,c leads to the solution Here the ratioũ s,c |ũ s,c | means that the variables u s,c andũ s,c have identical phases.
The amplitudes θ s,c are calculated as for Poissonian noise and as the non-negative solution of the cubic Cardan equation for the Gaussian noise [32].
The input-ouput model of this filtering we denote as u s,c = Φ 1 (ũ s,c , z s,c ), where the operator Φ 1 is defined by the Equations (15)- (18). Let us clarify meaning of the operations in Equation (15). Say c = R, then for all sensor pixels with the red filter, x ∈ Ω R , the amplitudes are updated using the corresponding R measurements, the first line in Equation (15). For all other pixels u s,R =ũ s,R , i.e., the amplitude of theũ s,R is not changed, the second line in Equation (15). In similar way, it works for all colors.
This preserving the signal value if the sensor is not able to provide the relevant observation is proposed and used in [32] for subsampled or undersampled observations. This rule can be interpreted as a complex domain interpolation of the wavefront at the sensor plane as well as a demosaicing algorithm for diffraction pattern observations. It is important to note that this observation processing is derived as an optimal solution for noisy data.
If P s,λ are orthonormal such that ∑ s P * s,λ P s,λ is the identity operator, ∑ s P * s,λ P s,λ = I, the solution takes the form 3. Minimization of L on b o,λ , ϕ o and δ λ (the last summand in the criterion Equation (14)) is the non-linear least square fitting of u o,λ by the parameters b o,λ , ϕ o and δ λ . We simplify this problem assuming that b o.λ |u o,λ |. Then, the criterion for this non-linear least square fitting takes the form: where ψ o,λ = angle(u o,λ ), i.e., the wrapped phase of u o,λ .
In this representation, the phase shifts δ λ are addressed to the wrapped phases ψ o,λ in order to stress that the complex exponent exp(j(ψ o,λ )) can be not in-phase with exp(j(µ λ · ϕ o )) and the variables δ λ serve in order to compensation this phase difference and make the phase modeling of the object by exp(j(µ λ · ϕ o )) corresponding to the complex exponent exp(j(ψ o,λ )).
The assumption b o,λ |u o,λ | is supported in our algorithm implementation by the initialization procedure enabling the high accuracy estimation of the amplitudes |u o,λ | in processing of separate wavelength observations.
The Absolute Phase Reconstruction (APR) algorithm is developed for minimization of L 1 on ϕ o and δ λ . The derivation and details of this algorithm are presented in Appendix.

Algorithm's Implementation
Combining the above solutions, the iterative algorithm is developed of the structure shown in Table 1. We use an abbreviation WM-APR for this algorithm as corresponding to Wavelength Multiplexing Absolute Phase Retrieval (WM-APR).
The initialization by the complex-valued u 1 o,λ is obtained from the observations {z s,c } separately for each wavelength by the SPAR algorithm [32] modified for the data provided by CFA. This initialization procedure is identical to the algorithm in Table 1, where Step 4 is omitted.
The main iterations start from the forward propagation (Step 1) and follows by the amplitude update and the interpolation for the wavefront pixels where there are no observations (Step 2). The operator Φ 1 for this update is defined by Equations (15) The derivation is based on the Nash equilibrium formulation for the phase retrieval instead of the more conventional constrained optimization with a single criterion function as it is in Equation (14). We do not show here this derivation as it is quite similar to developed in [32].

Sparsity and Block-Matching 3 Dimensions Filtering
The sparsity rationale assumes that there is a transform of image or signal such that it can be represented with a small number of transform coefficients or in a bit different terms with a small number of basis functions [58]. This idea is confirmed and supported by a great success of many sparsity based techniques developed for image or signal processing problems. Overall, the efficiency of the sparsity depends highly on selection of the transforms, i.e., basic functions relevant to the problem at hand. The so-called group-wise sparsity was proposed where similar patches of image are grouped and processed together Equation (22). It is obvious that the similarity inside the groups enhances the sparsity potential.
Within the framework of nonlocal group-wise sparse image modeling, a family of the BM3D algorithms has been developed where the both ideas grouping of similar patches and the transform design are taken into consideration. This type of the algorithms proposed initially for image denoising [54] being modified for various problems demonstrate the state-of-the-art performance [59,60].
Let us recall some basic ideas of this popular BM3D technique. At the first stage the image is partitioned into small overlapping square patches. For each patch a group of similar patches is collected which are stacked together and form a 3D array (group). This stage is called grouping. The entire 3D group-array is projected onto a 3D transform basis. The obtained spectral coefficients are hard-thresholded and the inverse 3D transform gives the filtered patches, which are returned to the original position of these patches in the image. This stage is called collaborative filtering. This process is repeated for all pixels of the entire wavefront and obtained overlapped filtered patches are aggregated in the final image estimate. This last stage is called aggregation. The details of BM3D as an advanced image filter can be seen in [54].
It follows from [61,62], that the above operations including the grouping define the analysis and synthesis transforms which can be combined in a single algorithm. The notation BM3D is used for this filtering algorithm.
In . It is implemented in this paper as independent filtering of amplitude and phase:  )). In real experimental samples, edges and varying features in phase and amplitude often co-exist in the same image pixels. Then, the joint phase/amplitude sparsity is a more efficient instrument for complex-valued object reconstruction. Novel forms of a BM3D based joint phase-amplitude processing are developed in [63], where the high-order singular value decompositions (HOSVD) are exploited for design of sparse basis functions for complex valued variables without separation phase and amplitude.
In Step 4 of the proposed algorithm, ϕ t+1 , the BM3D is applied for filtering of the real-valued ϕ o .
In our experiments the parameters of the algorithm are fixed for all tests. The parameters defining the iterations of the algorithm are as follows: γ 1 = 1/χ, where χ is the parameter of the Poissonian distribution, γ 1 /γ 2 = 0.2. The parameters of BM3D filters can be seen in [32].

Setup of Experiments
In numerical experiments, we model the lensless optical system (Figure 1a), where a thin transparent phase object is illuminated by monochromatic three color (RGB) coherent light beams from lasers. The wavelengths are Λ = [417, 532, 633] nm with the corresponding refractive indexes [1.528, 1.519, 1.515] as taken for BK7 optical glass. The reference wavelength is λ = 633 nm, then the relative wavelengths are µ λ = [1. 48, 1.17, 1]. The pixel sizes of CMOS camera and SLM are 1.4 and 5.6 µm, respectively. The distance d between the object and CMOS camera is equal to 5 mm. Figure 1b illustrates arrangement for each color channel of the CMOS camera [64]. A free propagation of the wavefronts from the object to the sensor is given by the Rayleigh-Sommerfeld model with the transfer function defined through the angular spectrum (AS) [65].
For the proper numerical AS propagation without aliasing effects, we introduce a zero-padding in the object plane of the size obtained from the inequality: which binds the propagation distance d, the number of pixels N in one dimension of the zero-padded object and the pixel size of the sensor ∆x [66]. For the distance d = 5 mm and the object 100 × 100 pixels, the zero-padded object has as a support 1700 × 1700 pixels. The intensities of the light beams registered on the sensor are calculated as Here AS λ,d denotes the AS propagation operator and M s are the modulation phase masks inserted before the object and pixelated as the object, • stands for pixel-wise multiplication of the object and phase masks. These phase masks enable strong diffraction of the wavefield introduced in order to achieve the phase diversity sufficient for reconstruction of the complex-valued object from intensity measurements. As it was described in [34], we use the Gaussian random phase masks. These phase masks are implemented on SLM. It follows that the pixels of these masks and SLM have the same size.
Thus, in our experiments the propagation operator P s,λ,d in Equation (4) is implemented as a combination of the angular spectrum propagation AS λ,d and the modulation phase mask M s .
To produce multiwavelength color images we need the three color components at each pixel location. It can be achieved by three sensors, each of which registers a light of the specific color. This optical setup for absolute phase retrieval from multiwavelength observations has been studied in our paper [52].
To reduce the complexity and the cost, most digital cameras use a single sensor covered by CFAs. In our simulation we assume that the popular RGB Bayer CFA pattern is used. The corresponding sampled color values are termed mosaic images, Figure 1. To render a full-color pattern at the sensor plane from these CFA samples, an image reconstruction process is applied as defined by Equation (15).

Reconstruction Results
The illustrating reconstructions are presented for two phase test-objects with the invariant amplitudes equal to 1 and the phases: "truncated Gaussian" and "logo" of Tampere University of Technology (TUT), 100 × 100 pixels both. The wrapped and absolute phases of these test-objects are shown in Figure 2. The TUT logo and truncated Gaussian phases are very different, the first one is discontinuous binary and the second one is piecewise-smooth continuous. The both test-objects are taken with a high peak-value of the absolute phase equal to 35π rad, that corresponds to about 35 reference wavelength λ in variations of the profile h o . We demonstrate the performance of the WM-APR algorithm for this quite difficult test-objects provided very noisy Poissonian observations. The noisiness of observations is illustrated by values of SNR (Equation (11)) and N photon (Equation (12)).
The accuracy of the object reconstruction is characterized by Relative Root-Mean-Square Error (RRMSE) calculated as RMSE divided by the root of the mean square power of the true signal: In this criterion the phase estimate is corrected by the mean value of the error between the estimate and true absolute phase. It is the standard approach to evaluation of the accuracy for the phase retrieval, where the phase is estimated within an invariant error.
Note, that values of RRMSE are identical for accuracy of the phase and the profile h o . It happens because these variables are different only by the invariant factors, which are the same in nominator and denominator of RRMSE. Figure 3 shows the performance of the developed WM-APR algorithm for different noise levels by a dependence of RRMSE from the noisiness parameter χ. The RRMSE curves for truncated Gaussian and logo TUT are given for different numbers of experiments S = 7 and S = 4, respectively. They demonstrate a similar behavior: go down for growing χ and take values about 0.1 at χ = 40, what corresponds to the very noisy observed data with a low value of SNR = 5.6 dB and small mean photon number N photon = 1.42.
RRMSEs for the logo TUT and Gaussian phases are shown in Figure 4 as functions of the experiment number S and SNR after 300 iterations. Nearly horizontal areas (dark blue in color images) correspond to high-accuracy reconstructions with small values of RRMSE, for other areas RRMSE values are much higher and the accuracy is not so good.
As it can be seen from these maps, reconstruction results for logo TUT demonstrate good small RRMSE values in wide ranges of SNR and S, while truncated Gaussian has such small RRMSE values only for S ≥ 7 and SNR > 5. It can be explained by a complex structure of the wrapped phase of the truncated Gaussian phase. It follows from Figure 4, that some improvement in the accuracy can be achieved at the price of the larger number of experiments S.
In what follows, we provide the images of the reconstructed absolute phases obtained for the iteration number T = 300, S = 4 and S = 7 for logo TUT and truncated Gaussian objects, respectively. The high iteration number T = 300 is governed by complexity of reconstruction from the CFA mosaic observations, where each pixel of the sensor captures only one red, green or blue color-channel information.
It means that the color information is subsampled with the subsampling ratio 50% for green and 25% for red and blue channels. The missed color-information is restored for each pixel of the sensor (demosaiced) in the iterations of the algorithm as defined by Equations (15)-(17). Figures 5 and 6 illustrate the absolute phase reconstruction (3D/2D images) obtained in two different ways: by the proposed WM-APR algorithm and by the SPAR algorithm reconstructing the wrapped phases for the separate wavelength channels following by the phase unwrapping using the PUMA algorithm [42]. The conditions of these experiments are: SNR = 6.4 dB, N photon = 1.78 for both test-objects.
The multiwavelength WM-APR algorithm demonstrates a clear strong advantage in reconstruction of the absolute phases while the processing based on separate wavelength channels completely failed. The accuracy of the WM-APR reconstructions is very high despite of the very noisy observations. The wrapped phase for the true truncated Gaussian phase should have a form of circular fringes for the parts of the image corresponding to the Gaussian peak. Instead of it in Figure 2, we can see quite irregular structures which look like destroyed by random errors or moiré patterns. This disturbance of the ideal circular fringes is an aliasing effect due to discretization (sampling) of the continuous Gaussian phase with a low sampling rate, which does not satisfy the Nyquist-Shannon conditions. It happens because the peak value of the truncated Gaussian phase is very high and the number of samples (100 × 100) is comparatively small. This aliasing disturbance of the Gaussian phase is so strong that even advanced 2D unwrapping algorithms, in particular, the PUMA algorithm, failed to reconstruct the absolute truncated Gaussian phase from the precisely noiseless sampled data shown in Figure 2. The impressive high-accuracy performance of the proposed WM-APR algorithm in this extra difficult situation is amazing as demonstrating a strong robustness of this algorithm with respect to errors caused by aliasing effects.
The reconstruction of the absolute phase for the binary logo TUT is also very difficult, in particular, due to huge differences of the phase values for adjacent pixels in the areas of the discontinuity of the phase. The sampling (pixelation) of this absolute phase also leads to the strong aliasing which can be noticed in 2D image in Figure 2. Again the unwrapping of the true noiseless wrapped phase by the 2D unwrapping algorithms is not possible, and then, the fail of the single wavelength based algorithms in Figure 6 is not surprising. Again, the amazing success of the WM-APR algorithm demonstrates that it is robust with respect to the sampling effects and is able to achieve the high-accuracy reconstruction in very difficult noisy scenarios.
In conclusion, we wish to note that the pair-wise beat-frequency synthetic wavelengths, one of the conventional methods for phase retrieval from multiwavelength observations (e.g., [43,46]), do not work in the considered scenario as the differences between the RGB wavelengths are not small.   From left to right: WM-APR algorithm, and phase reconstructions obtained separately for λ 1 , λ 2 , λ 3 wavelengths, respectively, followed by the 2D phase unwrapping.
For simulations we use MATLAB R2016b on a computer with 32 Gb of RAM and CPU with 3.40 GHz Intel R Core TM i7-3770 processor, 108.8 GFLOPS. The computation complexity of the algorithm is characterized by the time required for processing. For 1 iteration, S = 4, and 100 × 100 images, zero-padded to 1700 × 1700, this time equals to 12 s.

Conclusions
The multiwavelength absolute phase retrieval from noisy intensity observations is considered. The system is lensless with multiwavelength coherent input light beams and random phase masks applied for wavefront modulation. The light beams are formed by light sources radiating all wavelengths simultaneously and CMOS sensor equipped with a color filter array (CFA) is used for spectral measurements. The developed algorithm is based on maximum likelihood optimization in this way targeted on optimal phase retrieval from noisy observations. The algorithm for the Poissonian and Gaussian noise distribution is presented. The sparse modeling enables a regularization of the inverse imaging problem and efficient suppression of random observation errors by BM3D filtering of phase and amplitude. The effectiveness of the developed WM-APR algorithm is verified by simulation experiments demonstrated in this work only for the Poissonian noise. The high-accuracy performance of the algorithm is achieved for a high-dynamic range of absolute phase variations while the algorithms operating with the single wavelength observations separately failed. Physical experiments are our priorities for a further work. We consider nonidealities of SLM as a phase modulator and the AS propagation operator as the main challenges for the accuracy of profile and absolute phase reconstruction.

Conflicts of Interest:
The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.
In general case, when δ λ =δ λ , the in-phase situation is not achieved.
Going back to the criterion (Equation (A3)) we note that with this compensation it takes the form L 1 (ϕ o , δ λ ) = ∑ λ |||u o,λ |[exp(j(µ λ · ϕ o )) − exp(j(ψ o,λ −δ λ ))]|| 2 2 (A5) and using (Equation (A4)) it can be rewritten as Thus, the accurate compensation of the phase shifts δ λ in the wrapped phases ψ o,λ is achieved while the absolute phase ϕ o can be estimated within an unknown but invariant phase shift ∆ϕ o . This is not essential error as in the phase retrieval setup, ϕ o can be estimated only within an invariant shift.
The equivalence of the calculations (Equations (A1) and (A2)) to the phase manipulations shown in Equation (A4) cannot be proved analytically but careful numerical experiments confirm that they are quite precise even for noisy data.

2.
Minimization of Equation (A5) on ϕ o is a solution of the problem The calculations in Equation (A8) are produced pixel-wise for each pixel x independently. Note thatδ λ as well as the solutionφ o both depend on the unknown invariantδ λ .

3.
Minimization of Equation (A8) on δ λ leads to the following scalar numerical optimization: Finalization of estimation. Whenδ λ is found the optimal values forφ o are calculated aŝ In our algorithm implementation, the solutions for Equations (A6) and (A9) are obtained by the grid optimization.
It has been noted in our experiments that an essential improvement in the accuracy of the absolute phase reconstruction is achieved due to optimization on δ λ .