Multiwavelength surface contouring from phase-coded noisy diffraction patterns: wavelength-division optical setup

Abstract. We propose an algorithm for absolute phase retrieval from multiwavelength noisy phase coded diffraction patterns. A lensless optical system is considered with a set of successive single wavelength experiments (wavelength-division setup). The phase masks are applied for modulation of the multiwavelength object wavefronts. The algorithm uses the forward/backward propagation for coherent light beams and sparsely encoding wavefronts, which leads to the complex-domain block-matching three-dimensional filtering. The key-element of the algorithm is an original aggregation of the multiwavelength object wavefronts for high-dynamic-range absolute phase reconstruction. Simulation tests demonstrate that the developed approach leads to the effective solutions explicitly using the sparsity for noise suppression and high-accuracy object absolute phase reconstruction from noisy data.


Introduction
We consider lensless diffraction imaging where information on the object phase and amplitude in question is given by intensities of phase-coded diffraction patterns. The problem of the complex-valued object reconstruction from intensity measurements is well known as the phase retrieval, where "phase" emphasizes that this missing variable missing in observations defines the principle difficulties of the problem in particular, with respect to holography and interferometry, where reference beams allow to preserve phase information. The phase retrieval is from the class of the challenging illposed problems especially difficult for noisy observations. The phase/observation diversity sufficient for reliable reconstruction of the object phase from the intensity patterns is a crucial moment of the problem. Defocusing and phase modulation are known as effective instruments for solving the problem (e.g., Refs. [1][2][3][4]. Diffraction optical elements are used to design the diffraction phase imaging systems. 5 The Gerchberg-Saxton (GS) algorithms 6,7 are from the most popular in the field of phase retrieval. These iterative algorithms originally proposed the noiseless data use, alternating projection between the complex-valued object u o and complex-valued wavefronts u s at the sensor plane. In these iterations, the amplitudes in u s are replaced by square roots of the corresponding items of the measurements. The backprojection of the u s to the object u o is modified according to the prior information on the object, e.g., support size and shape, amplitude value, etc. The GS algorithms exist in various modifications. 8,9 The review and analysis of this type of the algorithms as well as further developments can be found in Ref. 10.
Many publications concern variational formulations for the GS-type algorithms. In particular, the links between the conventional GS and variational techniques are studied in Refs. 11-13. Contrary to the intuitively clear heuristic of the GS algorithms, the variational approaches usually have a strong mathematical background including image formation modeling, formulation of the objective function (criterion), and finally going to numerical techniques solving corresponding optimization tasks.
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/color light sources helps to analyze and recognize objects, to clarify color and tissue distributions of an object and dramatically improve the quality of imaging. Multiwavelength digital holography has an enhanced ability for three-dimensional (3-D) wide-range shape measurements using multiwavelength phase unwrapping, due to the recording of quantitative phase information with multiple wavelengths. In phase retrieval, the multiwavelength is an effective instrument enabling good phase/data diversity.
The multiwavelength phase retrieval is much less studied as compared with the standard single-wavelength formulation. These works by the principle of measurements can be separated into two groups. In the first one, the absolute phase is estimated from phase measurements obtained in some or another way. This scenario is typical for interferometry/holography where reference beams are applied to reveal the phase information, e.g., Refs. 14-17. The phase unwrapping algorithms for two-dimensional (2-D) images with simultaneous processing of multiple noisy complexexponent observations have been developed based on the maximum-likelihood techniques. 18 Another group of the techniques uses amplitudes or intensities (powers) as measurements. These formulations are from the class of the multiwavelength phase retrieval problems, e.g., Refs. [19][20][21][22]. The Chinese Remainder Theorems provide a class of the methods with a good theoretical background. 23 The reformulation of these approaches for more practical scenarios with noisy data and for robust estimation leads to the techniques similar to various forms of the maximum likelihood. 24,25 An approach to multiwavelength phase retrieval has been proposed in our recent papers. 26,27 The phase retrieval in the wavelength multiplexing (WM) setup is studied in the first of these two papers. A broadband light source radiates all wavelengths (RGB in our tests) simultaneously and a CMOS sensor equipped with a color filter array (CFA) registers spectral measurements. In the second paper, 27 we consider a wavelength division (WD) setup with successive and separate three RGB wavelength exposures registered by a broadband CCD sensor without CFA. The algorithms for these two setups are quite different, in particular because the WM algorithm requires interpolation of the RGB data subsampled by CFA for pixels where there are no observations of some of the wavebands. The WD algorithm deals with the observations obtained for all pixels for each of the RGB bands separately.
This paper is a further development of the WD setup 27 with more details concerning the approach and the algorithm as well as the simulation study performed for the complexvalued objects and parameters of the experiments different from those in Ref. 26. The main results of this paper can be summarized as an extension of the conference paper 27 with demonstration of the efficiency of the proposed algorithm, which is computationally simpler, much more accurate and faster than the algorithm for WM in Ref. 26. The experimental study of the algorithm is restricted to simulation tests.
The algorithms presented in Refs. 26 and 27 as well as in this paper can be treated as a development of the sparse phase and amplitude reconstruction (SPAR) algorithm 28 and the maximum likelihood absolute phase reconstruction for multiwavelength observations. 18 The paper is organized as follows. The multiwavelength object and image formation modeling as well as the Poissonian noise observations are discussed in Sec. 2. The development of the algorithm for the WD optical setup is a topic of Sec. 3. Simulation tests and results are discussed in Sec. 4.

Image Formation Model
The notation and the image formation model presented in this section follow those in our paper. 26 We present these results here for completeness of the presentation and in order to make clearer relations of the problems considered in Ref. 26 and in this paper.

Multiwavelength Object and Image Modeling
Let h o ðxÞ, x ∈ R 2 , be a profile of the transparent 2-D object to be reconstructed. A coherent multiwavelength light beam generates corresponding complex-valued wavefronts u o;λ ¼ b o;λ expðjφ o;λ Þ, λ ∈ Λ, where Λ is a set of the wavelengths and φ o;λ ¼ 2π λ h o ðn λ − n env Þ are phase delays corresponding to h o , where n λ is the refractive index of the object and n env is the refractive index of the environment, for the air n env ¼ 1, what is assumed in our tests.
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;λ for separate λ.
Here λ 0 ∈ Λ is a reference wavelength and φ o is an absolute phase corresponding to this wavelength. The parameter μ λ establishes a link between the absolute phase φ o and 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 Wð·Þ as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 3 2 6 ; 4 4 2 (2) Wð·Þ 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 are defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 3 2 6 ; 3 1 7 z s;λ ¼ Gfy s;λ g; s ¼ 1; : : : ; S; λ ∈ Λ; where u s;λ is a wavefront propagated to the sensor plane, P s;λ;d fg is an image (diffraction pattern) formation operator, i.e., the propagation operator from the object to the sensor plane, including in particular random phase masks used for wavefront modulation, d is a propagation distance, y s;λ is the intensity of the wavefront at the sensor plane, and z s;λ are noisy observations as defined by the generator Gfg of the random variables corresponding to y s;λ , and S is a number of experiments. The considered multiwavelength phase retrieval problem consists of reconstruction of φ o and b o;λ from the observation z s;λ provided that μ λ , P s;λ;d fg, and Gfg are known. If the phases μ λ φ o ∈ ½π; −πÞ, the problem is 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. 29 However, if the phases μ λ φ o go beyond the range ½π; −πÞ, the problem becomes nontrivial and challenging since a phase unwrapping must be embedded in the multiwavelength phase retrieval.
We are focused on the reconstruction of the absolute phase φ o from a set of single wavelength experiments produced in a successive manner, which are registered by a wide-band gray-scale CCD sensor, instead of using all wavelengths simultaneously with CFA and demosaicing processing as it is done in the WM algorithm derived in Ref. 26.
The proposed WD algorithm demonstrates a faster convergence than the WM algorithm, 26 but makes the experimental realization more complex as it requires separate experiments and measurements for each wavelength. A single wavelength is used for each experiment and the obtained observations fz s;λ g are used jointly for reconstruction of φ o ðxÞ and h o ðxÞ.
Note also that the proposed approach assumes that the refractive indexes of the object are known for the considered set of the wavelengths. Only in this case, the problem of the absolute phase retrieval and contouring can be resolved.

Noisy Observation
For noise simulation, we use independent Poisson random variables since a measurement process in optics amounts to count the photons hitting the sensor's elements (e.g., Ref. 30). The presented form of the noisy observations corresponds to Refs. 26 and 27.
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 as where y s;λ ðxÞ is the intensity of the wavefront at the pixel x. The parameter χ > 0 in Eq. (6) can be interpreted as an exposure time and as sensitivity of the sensor with respect to the input radiation. Defining the observation signal-to-noise ratio (SNR) as the ratio between the square of the mean and variance of z s;λ ðxÞ, we have SNR ¼ E 2 fz s;λ ðxÞg∕varfz s;λ ðxÞg ¼ 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. Another useful noise characteristic is the mean value of photons per pixel. It is defined as N photon;λ ¼ P x z s;λ ðxÞ∕N sensor , where N sensor is the number of sensor pixels. As in case of SNR, smaller values of χ lead to smaller N photon , i.e., to noisier observations z s;λ ðxÞ.

Development of Algorithm
We consider the problem of the absolute phase retrieval as an estimation of the object wavefront u o ∈ C n from noisy observations fz s;λ g. The problem is challenging mainly due to the periodic nature of the likelihood function with respect to the phase φ o and the nonlinearity of the observation model.
Following Ref. 26 and using the same notation, the maximum likelihood function leads to the criterion function E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 3 2 6 ; 7 5 2 where γ 1 ; γ 2 > 0 are regularization parameters, and k · k 2 2 stands for the Hadamard norm.
Here we introduce invariant additive errors δ λ , which model phase errors for different wavelengths that occur due to uncertainty in the phase retrieval problem with respect to invariant errors in the phase appearing for each wavelength. The criterion Eq. (7) is minimized with respect to u s;λ , u o;λ , b o , and φ o as well as with respect to the phaseshifts δ λ .
The algorithm is composed from the following steps. where I stands for the identity operator.
Here the operator P Ã s;λ;d is Hermitian adjoint for P s;λ;d . If P s;λ;d are orthonormal such that P s P Ã s;λ P s;λ;d is the identity operator, P s P Ã s;λ;d P s;λ;d ¼ I, then the solution is simplified to the form E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 3 2 6 ; 2 3 7û 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 exponents exp½jðψ o;λ Þ can be out-ofphase 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;λ Þ.
Minimization of L 1 , with respect to φ o and the phaseshifts δ λ , is implemented as minimization on the grids of these variables. The derivation of this algorithm and its details can be seen in Ref. 26, where it was done for the WM optical setup. This algorithm was named absolute phase reconstruction, and it is shown in step 4 of Table 1 with the abbreviation APR.

Algorithm's Implementation
Using the above solutions, the iterative algorithm is developed of the structure shown in Table 1. The initialization by the complex-valued u 1 o;λ is obtained from the observations fz s;λ g by the SPAR algorithm 28 separately for each wavelength. The main iterations start from the forward propagation (step 1) and follow by the amplitude update for u t s;λ at step 2. The operator Φ 1 derived in Ref. 28 The backpropagation is realized in step 3, and the operator Φ 2 is defined by Eq. (9). The absolute phase reconstruction from the wrapped phases of u tþ1 o;λ is produced in step 4 by the APR algorithm. 26 The obtained amplitude and phase update u tþ1 o;λ at step 5.
The steps 3 and 4 are completed by the block-matching 3-D (BM3D) filtering. 29 In step 3, it is the filtering of complex-valued u tþ1∕2 o;λ produced separately for the wrapped phase and amplitude of u tþ1∕2 o;λ . In step 4, this filtering is applied to the absolute phase φ tþ1∕2 o . These BM3D filters are derived from the groupwise sparsity priors for the filtered variables. This technique 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 Eq. (7). We do not show here this derivation as it is quite similar to that developed in Ref. 28. The sparsity rationale assumes that there is a transform of image/ 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 basic functions. 31 This idea is confirmed and supported by the great success of many sparsity-based techniques developed for image/signal processing problems. Within the sparse theory, a family of the BM3D algorithms has been developed where both ideas of grouping similar patches and the transform design are taken into consideration. This type of the algorithms proposed initially for image denoising 29 being modified for various problems demonstrates the state-of-the-art performance. 32 The details of BM3D as an advanced image filter can be seen in Ref. 29. In step 3 of the proposed algorithm, u tþ1 o;λ ¼ BM3Dðu tþ1∕2 o;λ Þ, the BM3D is applied to the complex-valued variables u tþ1∕2 o;λ . It is implemented in this paper as independent filtering of amplitude and wrapped phase The algorithm presented in Table 1 differs from the WM algorithm 26 by steps 2 and 3. The main iterations in Ref. 26 also start from the forward propagation and follow by the amplitude update and the interpolation for the CFA sensor pixels where there are no observations (step 2). In the WD algorithm, there are no updates and interpolation. In step 3 in Ref. 26, the wavefront pixels not given in observations preserve values obtained in step 2 and used for backpropagation. This manipulation with wavefront reconstruction for RGB pixel is not given in observations from the principal difference with the WD algorithm, presented here in Table 1. The other steps of the algorithms for WM and WD, as it is in this paper, concerning filtering and forward-and backward propagation are identical.
Comparing the presented WD algorithm developed for the WD setup with that developed in Ref. 26 for the WM setup, we may note that the faster convergence and better accuracy are obtained for WD, but it requires a larger number of observations, because each experiment gives the data of

5.
Object wavefront update: o;λ a single wavelength only, while in WM we obtain the observations for all wavelength simultaneously.

Proposed Setup
In simulation, we model a lensless optical system (Fig. 1), where a thin transparent phase object is illuminated by monochromatic three color (RGB) coherent light beams from lasers or LEDs. 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 λ 0 ¼ 417 nm then the relative frequencies take values μ λ ¼ ½1; 0.7705; 0.6425. The pixel sizes of the 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. The free propagation of the wavefronts to the sensor is given by the Rayleigh-Sommerfeld model with the transfer function defined through the angular spectrum (AS). 33 For the proper numerical AS propagation without aliasing effects, the zero-padding of the object must be applied. Its size is obtained from the inequality 34 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 6 3 ; 3 9 7 N ≥ λ · d Δx 2 ; which binds the propagation distance d and number of pixels N in one-dimension of the zero-padded object for a given values of pixel size of the sensor Δx and λ. For the distance d ¼ 5 mm, Δx ¼ 1.4 μm, λ ¼ 633 nm, and object 100 × 100 pixels the zero-padded object has a support of 1700 × 1700 pixels. The intensities of the light beams registered on the sensor are calculated as z s;λ ¼ GfjAS λ;d fM s ∘ u o;λ gj 2 g, s ¼ 1; : : : ; S, λ ⊂ Λ. 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 the pixelwise multiplication of the object and phase masks. These phase masks enable strong diffraction of the wave-field and are introduced in order to achieve the phase diversity sufficient for reconstruction of the complex-valued object from intensity measurements. As it was described in Ref. 33, we use the Gaussian random phase masks. Thus, in our simulations the propagation operator P s;λ;d in Eq. (3) is implemented as a combination of the AS propagation AS λ;d and the modulation phase mask M s .

Reconstruction Results
The illustrating reconstructions are presented for two phase objects with the invariant amplitudes equal to 1 and the phases: Gaussian distribution (100 × 100) and US Air Force (USAF) resolution test-target (64 × 64). The absolute and wrapped phases of these test-objects are shown in Fig. 2. The Gaussian and USAF phases are very different, the first one has a smooth continuous shape while the second one is discontinuous binary. Both phases are taken with the high peak-value equal to 30π rad, what corresponds to about 30 reference wavelengths λ 0 in variations of the profile h o .
As a result of this high peak-value of the absolute phases, the corresponding wrapped phases are very complex with fringes overlapping for the Gaussian phase and lack of height information in the USAF test-phase. Due to this complexity, the unwrapping is not possible using single frequency wrapped phases only. We show that the proposed multiwavelength WD algorithm is quite successful and is able to reconstruct the absolute phase even from very noisy data.
We demonstrate the performance of the WD algorithm for the very difficult scenarios with very noisy Poissonian observations. The noisiness of observations is characterized by SNR and by the mean number of photons per sensor pixel, N photon .
The accuracy of the object reconstruction is characterized by the relative root-mean-square error (RRMSE) criteria calculated as RMSE divided by the root of the mean square power of the signal   Figure 3 shows the performance of the WD algorithm with respect to different noise levels characterized by the parameter χ. The RRMSE curves for Gaussian and USAF phases demonstrate a similar behavior and go down for growing χ numbers, but RRMSE curve for the USAF object goes down more sharply and takes values smaller than 0.1 at χ ¼ 20, while RRMSE curve for Gaussian phase takes a close value only at χ ¼ 50. Nevertheless, in both cases, the reconstructions are nearly perfect even for very noisy observed data with SNR values as low as 3.8 and 6.5 dB and very small photon numbers N photon ¼ 0.75 and 1.87, respectively.
RRMSEs for the Gaussian phase are shown in Fig. 4 as functions of the experiments number S and SNR. Number of experiments S is in range of 1 to 10, and SNR from 1 (noisiest case) to 25 dB (noiseless case). Nearly horizontal (dark blue) areas 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. For 50 iterations (left image), the high accuracy can be achieved starting from S ¼ 8 even for very noisy data with SNR ¼ 4. For the number of experiments smaller than S ¼ 4, the algorithm fails to converge even for low noise levels due to a not sufficient diversity of the observed diffraction patters. The smallest number of experiments that can be used for reconstruction is S ¼ 4. However, even for S ¼ 4 good results can be obtained only after 100 iterations (right image). It follows from Fig. 4 that some improvement in the accuracy can be achieved at the price of the larger number of experiments S and the larger number of iterations.
In comparison with the WM algorithm 26 such small values for RRMSE can be obtained only after a larger number of iterations, about T ¼ 300, mainly because of 75%, 50%, and 75% of CFA pixels lacking intensity measurements for red, green, and blue wavelengths, respectively. This information is recovered by the WM algorithm but at the price of the larger number of iterations.  In what follows, we provide the images of the reconstructed absolute phase obtained for S ¼ 4 and the iteration number T ¼ 100. Figures 5 and 6 show absolute phase reconstructions (3-D/2-D images) obtained by the WD algorithm and by the SPAR algorithm reconstructing the wrapped phases for the separate wavelengths following by the phase unwrapping by the PUMA phase unwrapping algorithm. 18 The conditions of the experiments are: SNR ¼ 6.5 dB, N photon ¼ 1.87 for the Gaussian object and SNR ¼ 3.8 dB, N photon ¼ 0.75 for the USAF object. The WD algorithm demonstrates a strong ability to reconstruct the absolute phases while the single wavelength-based approach completely failed. The accuracy of the WD reconstruction is very high despite a high level of the noise.
The wrapped phase pattern for such high peak-value of the object phases is very complex and irregular (see Fig. 2) so that it is not possible to unwrap it by modern 2-D unwrapping algorithms, but the proposed algorithm is able to resolve the problem even for such complex case. Especially it is challenging task to reconstruct objects such as USAF with big differences between adjacent pixels exceeding 2π, but the WD algorithm successfully reconstructs both objects with high reconstruction quality.
For simulations, we use MATLAB R2016b on a computer with 32 GB of RAM and CPU with a 3.40 GHz Intel ® Core TM i7-3770 processor. The computation complexity of the algorithm is characterized by the time required for processing. For 1 iteration, S ¼ 4, and 100 × 100 pixels images, zero-padded to 1700 × 1700, this time equals to 12 s.

Conclusion
The multiwavelength absolute phase retrieval and surface contouring from noisy intensity observations is considered. The observations are recorded by a broadband monochromatic sensor in the successive manner for each wavelength in the WD optical setup. The maximum likelihood criterion used in the developed variational derivation of the phase retrieval algorithm defines the general intention to reach statistically optimal estimates. The phase retrieval is an illposed inverse problem where the observation noise is amplified and transferred to phase and amplitude as variables of optimization. The sparse modeling enables a regularization of this inverse problem and efficient suppression of these random errors by BM3D filtering of phase and amplitude. The efficiency of the developed algorithm is demonstrated by simulation tests for the coded diffraction pattern WD scenario.

Disclosures
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.