Multi-frequency phase retrieval from noisy data

Multi-frequenc y phase retrie v al from noisy data Vladimir Katko vnik and Karen Egiazarian Laboratory of Signal Processing, T ampere Uni versity of T echnology , P .O. Box 553 FI-33101 T ampere, Finland e-mails: vladimir .katkovnik@tut.fi, karen.egiazarian@tut.fi Abstract —The phase r etrieval from multi-frequency intensity (power) observations is considered. The object to be recon- structed is complex-valued. A nov el algorithm is presented that accomplishes both the object phase (absolute phase) retriev al and denoising for Poissonian and Gaussian measurements. The algorithm is derived from the maximum likelihood formulation with Block Matching 3D (BM3D) sparsity priors. These priors result in two filtering: one is in the complex domain for complex- valued multi-fr equency object images and another one in the real domain f or the object phase. The algorithm is iterative with alternating projections between the object and measurement variables. The simulation experiments are produced for Fourier transform image formation and random phase modulations of the object, then the observations are random object diffraction patterns. The results demonstrate the success of the algorithm for reconstruction of the complex phase objects with the high- accuracy performance even for very noisy data. I. I NTRODUCTION Phase retriev al is a problem concerning reconstruction of the phase from intensity (power) measurements of complex- valued variables. It is inherently ill-posed due to the lack of the phase in observations. The fields of application are varying from physics and engineering to medicine and biology , e.g. [1]-[3]. W e consider the problem of 2 D imaging from the multi- frequency observations. The multi-frequency complex-valued object model is of the form u o, λ = b o, λ exp( j μ λ ϕ ) , λ ∈ Λ , (1) where u o, λ ∈ C N × N , ϕ is the object absolute phase in radians, μ λ > 0 are dimensionless relative frequencies and Λ =[ λ o , λ 1 , ..., λ n λ − 1 ] is a set of periods (wav elengths) of the multi-frequency observations. In what follows, the amplitude, phase and other v ariables are functions of the argument x given on a regular 2 D grid, X ⊂ Z 2 . In the model (1), μ λ is a frequency depending scale para- meter of the object phase which establishes a link between the absolute phase ϕ and the wrapped phase ψ λ of u o, λ which can be measured at the λ -channel. The wrapped phase is related with the true absolute phase ϕ as μ λ ϕ = ψ λ +2 π k λ , where k λ is an integer , ψ λ ∈ [ − π , π ) . The link between the absolute and wrapped phase is installed as: ψ λ = W ( μ λ ϕ ) ≡ mo d( μ λ ϕ + π , 2 π ) − π . where W ( · ) is the wrapping operator, which decomposes the absolute phase μ λ ϕ into two parts: the fractional part ψ λ ∈ [ π , − π ) and the integer part defined as 2 π k λ . The frequencies, i.e. the phase factors μ λ a r eg i v e nb u tt h e phase ϕ and the amplitude b o, λ should be retriev ed from the intensity observations, which are for the noiseless y s, λ and noisy z s, λ observations hav e the form: u s, λ = P s, λ u o , y s, λ = | u s, λ | 2 , (2) z s, λ = G{ | u s , λ | 2 } , λ ∈ Λ , s =1 , ..., S , (3) where S is a number of experiments. The image formation operator P s, λ , usually depending on λ , is also known. The observations are giv en for each λ ∈ Λ . If the phases μ λ ϕ are interferometric, i.e. μ λ ϕ ∈ [ π , − π ) , then the problem becomes trivial as the estimates of μ λ ϕ and b o, λ can be found separately for each λ by one of the phase retriev al algorithms applicable for complex domain imaging. The problem becomes much more challenging, when the phases μ λ ϕ are non-interferometric, μ λ ϕ / ∈ [ π , − π ) because the data processing produced independently for each λ giv es estimates only for the wrapped phases ψ λ . Aggregation of these estimates in order to reconstruct the absolute phase ϕ is a crucial issue of the multi-frequency phase retriev al problem. A. Related works For a single frequency , μ λ =1 and b o, λ = b o in (1), the multi-frequency problem becomes the con ventional phase retriev al: to reconstruct the object u o ∈ C N × N from noiseless y s or noisy z s intensity observations. A flow of publications old and recent concern various formulations and algorithms for solution of this con ventional setup starting from the classic Gerchberg-Saxton (GS) style algorithms (e.g. [4], [5]) based on alternating projections between the complex-v alued u o and u s . In this iterations the amplitudes in u s are tuned to the given observations. The back projections are modified accordingly to the prior information on the object such as a support size and shape, amplitude of u o , etc. The revie w and analysis of the GS algorithms as well as further developments can be seen in [6]. The variational formulations of the phase retrieval have a stronger mathematical background and go to solving optimiza- tion tasks [7], [8]. The maximum likelihood approaches using the transform domain phase/amplitude sparsity are proved to be effecti ve for noisy data. It is a base of the Sparse Phase Amplitude Retrie val (SP AR) algorithms [9]-[12]. The multi-frequency phase retrieval is a much less studied problem. The lower synthetic frequencies (beat-frequencies) appear in the differences of the wrapped phases Δ λ . λ  = W ( μ λ ϕ ) − W ( μ λ  ϕ )= W (( μ λ − μ λ  ) ϕ ) . The problem can be solved in the straightforward manner provided that μ λ − μ λ  is small and ( μ λ − μ λ  ) ϕ ∈ [ π , − π ) , then ϕ can be easily calculated from Δ λ . λ  . Howe ver , the noise lev el in these synthetic frequency estimates gro ws proportionally to 1 / ( μ λ − μ λ  ) [13]. It is a drawback of this approach. The various forms of t he synthetic frequency techniques are pro- posed for multi-frequency scenarios (e.g. [14], [15]). The 2 D phase unwrapping algorithms, i.e. absolute phase retrie val with simultaneous processing of multi-frequency noisy complex- exponents, hav e been de veloped in terms of the maximum likelihood formulation [16], [17]. The Chinese Remainder Theorem provides a wide class of the methods with a good the- oretical background [18]. Reformulation of these approaches to noisy data and for robust estimation leads to the techniques similar to the maximum likelihood methods [19], [20]. In this paper we propose the approach and deriv e the algorithm which can be treated as development of the SP AR algorithm [11]-[12] and the maximum likelihood absolute phase reconstruction from multi-frequency observations [17]. B. Contribution of this paper W e present a novel multi-frequency phase retrie val algo- rithm that accomplishes both phase unwrapping and denoising from intensity observations. The noise models account for Poissonian and Gaussian distributions. The proposed algorithm is iterative with alternating projections between object and measurement variables. The algorithm is deriv ed from the maximum likelihood formulation with the Block Matching 3D (BM3D) sparsity priors. These priors result in the two types of filtering: one is in the complex domain for multi-frequency object images and another one is in the real domain for the object absolute phase (object shape). One of the paramount issues of the algorithm is the absolute phase reconstruction by aggregation of the object multi-frequency estimates. This aggregation is achiev ed by the least square solution including compensation of the in variant phase-shifts in the λ -channel object estimates. The simulation experiments are produced for Fourier trans- form image formation with the random phase modulations of the object. The observ ations are noisy random object diffraction patterns. The results demonstrate the success of the algorithm for reconstruction of the high complexity phase objects with the precise performance even for very noisy observations. II. A LGORITHM DEVELOPMENT A. Pr oblem formulation Reconstruction of u o ∈ C N × N from noisy { z s, λ } is rather challenging mainly due the periodic nature of the likelihood function with respect to the phase ϕ and the non-linearity of the observ ation model. Provided a stochastic noise model with independent samples, the maximum likelihood leads to the basic criterion L 0 =  λ ∈ Λ  S s =1  x l ( z s, λ ( x ) , | u s, λ ( x ) | 2 ) , where l ( z, | u | 2 ) denotes the minus log-likelihood of a candi- date solution for u o giv en through the observed true intensities | u | 2 and noisy outcome z . For the Poissonian and Gaussian distributions we hav e, respectiv ely , l ( z, | u | 2 )= | u | 2 χ − z log( | u | 2 χ ) , χ > 0 is the noise scaling parameter [11] and l ( z, | u | 2 )= 1 2 σ 2 || u | 2 − z | 2 , σ 2 stands for the noise variance. W e introduce the following criterion including the image formation model P s, λ for u s, λ and the complex-v alued expo- nent modeling of the multi-frequency object u o, λ : L ( u s, λ , u o, λ , b o , ϕ , δ λ )=  λ ,s,x l ( z s, λ ( x ) , | u s, λ ( x ) | 2 )+ 1 γ 1  λ ,s || u s, λ − P s, λ u o, λ || 2 2 + (4) 1 γ 2  λ || b o. λ exp( j ( μ λ · ϕ + δ λ )) − u o, λ || 2 2 , where γ 1 , γ 2 > 0 are parameters, and | |·| | 2 2 is the Hadamard norm. The maximum likelihood in (4) is penalized by the quadratic residual function for correspondence of u s, λ to u o, λ . The estimates of μ λ ϕ ( x ) in (4) can be biased due to the fact that in the phase retrieval the solution is defined only within an inv ariant additiv e phase shift. W e model this biasedness by the parameters δ λ in v ariant with respect to x and additi ve to the phases μ λ ˆ ϕ . W e say that the object complex exponents u o, λ are in-phase (synchronized) if δ λ =0 , for all λ ∈ Λ , and out-of-phase otherwise. The phase shifts δ λ should be compensated for a proper estimation of the absolute phase ϕ . The criterion L is minimized with respect to u s, λ , u o, λ , b o , ϕ as well as with respect to the phase-shifts δ λ . 1) Minimization with respect to u s, λ concerns the first two summands in L . The problem is additiv e and can be obtained separately for each x , s and λ .T h e corresponding analytical solutions derived for both the Poissonian and Gaussian distributions can be seen in [11]. These solutions optimal in the maximum likelihood sense define u s, λ as functions of noisy observation z s, λ and projections P s u o, λ of u o, λ on the sensor . 2) Minimization with respect to u o, λ goes to the last two summands of the criterion. It is the quadratic problem. If P s, λ are orthonormal such that  s P ∗ s, λ P s, λ = I , where I is the identity operator and P ∗ s, λ is Hermitian adjoint for P s, λ , the solution is of the form ˆ u o, λ =  s P ∗ s, λ u s, λ + γ 1 / γ 2 b o, λ exp( j μ λ ϕ ) 1+ γ 1 / γ 2 . . (5) 3) Minimization on b o, λ , ϕ and δ λ concerning the last summand in L is a non-linear least square fitting of the frequency dependent u o, λ . Minimization on ϕ is an absolute phase estimation, i.e. the phase unwrapping problem. The solution is composed from the following two successive stages A and B . (A) Phase synchr onization . Let λ  ∈ Λ be a reference channel. Define the estimates for the phase shift δ λ for (4) in the following way ˆ δ λ = W ( ˆ δ λ  · μ λ / μ λ  + δ λ , λ  μ λ ) , (6) where ˆ δ λ  is a hypothetical value of unknown δ λ  and δ λ , λ  is a shift between the weighted wrapped phases of the reference and λ -channels δ λ , λ  = median x ( W ( ψ o, λ ( x ) / μ λ − ψ o, λ  ( x ) / μ λ  )) . (7) One of the goals is to reduce the number of unknown phase shifts δ λ , λ ∈ Λ ,t oas i n g l e δ λ  , i.e. to the phase- shift in the reference channel. It was tested by Monte-Carlo modeling that δ λ , λ  is a good estimate for δ λ / μ λ − δ λ  / μ λ  , then ˆ δ λ = δ λ + ( ˆ δ λ  − δ λ  ) μ λ / μ λ  . Inserting this ˆ δ λ in (4) instead of δ λ we arri ve to L 1 ( b o, λ , ϕ , δ λ )= (8)  λ || b o. λ exp( j ( μ λ ( ϕ + Δ ϕ )) − | u o, λ | exp( j ( ψ o, λ − δ λ )) || 2 2 , Δ ϕ =( ˆ δ λ  − δ λ  ) / μ λ  . The accurate synchronization of the exponents exp( j ( μ λ ( ϕ + Δ ϕ )) is achie ved in (8) as a result of compensation by ˆ δ λ the phase shifts existing in the wrapped phases ψ o, λ .I f Δ ϕ  =0 the phase ϕ is estimated within an inv ariant Δ ϕ . It is not impose any restrictions as by default the phase ϕ in the phase-less measurements of the phase retriev al can be estimated within an arbitrary inv ariant shift only . (B) Absolute phase r etrieval (APR). In order to simplify the presentation assume that the complex exponents of u o, λ are perfectly in phase, i.e. Δ ϕ =0 ,a n d b o. λ = | u o, λ | , then the absolute phase is reconstructed as ˆ ϕ = (9) arg min ϕ  λ | u o, λ | 2 [1 − cos( μ λ · ϕ − ˆ ψ o, λ )] , where ˆ ψ o, λ − ˆ δ λ is denoted as ˆ ψ o, λ , i.e. the δ λ unknown in (8) are replaced by their estimates ˆ δ λ . For calculation of ˆ ϕ we use the approach proposed in [17]. Due to the periodicity of the cosines in (9) μ λ · ϕ = ˆ ψ o, λ +2 π k λ , where k λ are integers. By summation ov er λ we obtain ϕ =  λ ( ˆ ψ o, λ +2 π k λ )  λ μ λ and substituting this ϕ into (9) we arrive to ˆ ϕ =a r g m i n k ∈ [0 ,Q )  λ | u o, λ | 2 [1 − (10) cos( μ λ (  λ ˆ ψ o, λ +2 π k )  λ μ λ − ˆ ψ o, λ )] ,k =  λ k λ , ˆ b o, λ = | u o, λ | . (11) It shows that the multiv ariable optimization on k λ is reduced to the scalar optimization on the integer k .I f μ λ = p λ /q λ , where ( p λ ,q λ ) are coprime integer , then the criterion in (10) is a periodic function of k with the synthetic period Q equal to the nominator of  λ μ λ .I t makes clear why we restrict the optimization in (10) to the interval [0 ,Q ) . W e obtain from the above formula for ϕ that the maximal range of the absolute ϕ which can be reconstructed in this approach depends on the spectrum of the multiple frequencies and restricted by 2 π ( n Λ + Q ) /  λ μ λ , where n Λ is a number of the frequencies in Λ . Rational approximations are used for μ λ in these calculations if μ λ are not rational. B. Algorithm’s implementation Using the above solutions the iterative algorithm is dev el- oped of the structure shown in T able I. The initialization by the complex-valued u 1 o, λ is obtained from the observations { z s, λ } by the SP AR algorithm [11] separately for each λ .T h e main iterations start from the forward propagation (Step 1) and follows by the amplitude update for u t s, λ at Step 2. The operator Φ 1 for this update is obtained by minimization of L on u s, λ . It can be seen in [11] for Poissonian and Gaussian noise models. The back propagation in Step 3, the operator Φ 2 , is defined by Eq.(5). The absolute phase reconstruction from the wrapped phases of u t +1 o, λ is produced in Step 4 by the ARP algorithm as defined by the solutions (10)-(11). The obtained amplitude and phase update u t +1 o, λ at Step 5. The number of iteration is fixed in this implementation of the algorithm. The steps 3 and 4 are completed by the BM3D filtering. 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 . These BM3D filters are derived for the considered phase retriev al problem from the group-wise sparsity priors for the filtered variables. This technique is based on the Nash equilibrium formulation for the phase retriev al instead of the more conv entional constrained optimization with a single criterion function as it is, for instance, in (4). W e do not show this deriv ation as it is quite similar to the deriv ation presented in [11]. TA B L E I M UL TI -F REQUENCY A BSOLUTE P HASE R ETRIEV AL (MF-APR) A LGORITHM Input : { z s, λ } , s =1 , ..., S , λ ∈ Λ , Initialization : u 1 o, λ , λ ∈ Λ Main iterations: t =1 , 2 , ..., T 1. Forward propagation : u t +1 / 2 s, λ = P s, λ u t o, λ , s =1 , ..., S , λ ∈ Λ ; 2. Noise suppression and update of u t s, λ : u t s, λ = Φ 1 ( u t +1 / 2 s, λ , z s, λ ) ; 3. Backward propagation and filtering: u t +1 / 2 o, λ = Φ 2 ( u t s ) , u t o, λ = BM 3 D ( u t +1 / 2 o, λ ); 4. Absolute phase retrieval and filtering: { ϕ t +1 / 2 , b t +1 o. λ } = AP R ( u t o, λ ) , ϕ t +1 = BM 3 D ( ϕ t +1 / 2 ); 5. Object update: u t +1 o, λ = b t +1 o, λ exp( j ϕ t +1 μ λ ) , λ ∈ Λ ; Output : ϕ T +1 ,b T +1 o, λ , u T +1 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 [11]. The MA TLAB demo-codes of the algorithm are publicly av ailable 1 . III. N UMERICAL EXPERIMENTS In numerical experiments we model a lensless optical sys- tem where a thin flat transparent phase object is illuminated by monochromatic three color RGB (red-green-blue) coherent light beams from lasers or LEDs. The wa vefronts just behind the object are u o, λ = b o, λ exp( j 2 π · h · ( n λ − 1) / λ ) , where h ≥ 0 is an object thickness (height, shape, profile), λ is the wav elength, and n λ is a refracti ve index of the object material. W e recalculate the wav elengths including the effect of the refractive index as λ → λ / ( n λ − 1) ∈ Λ , Λ = [500 , 600 , 700] × 10 − 9 m. Then, the formula 2 π · h/ λ defines the phase delay in propagation of the coherent light-beam through the object. In notation used in (1), the relativ e frequencies are μ λ = λ ref / λ , where λ ref is a reference wav elength, and ϕ =2 π h/ λ ref is the reference absolute phase in radians. The propagation of the wav efronts to the sensor is gi ven by the Fourier transform [21]. The intensities of the light beam impinging on the sensor parallel to the object plane are calculated as z o, λ = G {|F { M s ◦ u o, λ }| 2 } , s =1 , ..., S , λ ∈ Λ . Here M s are the modulation phase-masks inserted next to the object,  ◦  stands for the pixel-wise multiplication of the object and phase-mask transition functions. The phase-masks M s enable strong diffraction of the wavefield propagation and are introduced in order to achiev e the phase di versity sufficient for improved reconstruction of the complex-valued object. Similar to [22], we use the four phase values [0 , π / 2 , π , 3 / 2 π ] prescribed to the pixels of the modulation phase- masks randomly . The accuracy of the object reconstruction is characterized by the Relative Root-Mean-Square Error ( RRM S E ) criteria calculated as RM S E divided by the root of the mean square power of the signal. In these criteria the biasedness of the phase estimate is corrected by the mean value of the error between the estimate and the true value. W e show the results obtained for Poissonian observa- tions and very noisy data. The noisiness of observ ations is characterized by Signal-to-Noise Ratio (SNR) SN R = 10 log 10 ( χ 2  s, λ || u s, λ || 2 2 /  s, λ || χ ·| u s, λ | 2 − z s, λ || 2 2 ) dB and by the mean number of photons per sensor pixel, N photon . The illustrating results are presented for two phase objects with the in v ariant amplitude equal to 1 : Gaussian ( 100 × 100 ) and U.S. Air Force (USAF) resolution target ( 612 × 612 ). For the giv en Λ and the reference wav elength λ ref = λ 1 we have μ 1 =1 , μ 2 =0 . 0833 , μ 3 =0 . 7143 . The rational approximations for these μ λ are μ 1 =1 , μ 2 =5 / 6 , μ 3 = 5 / 7 and the sum  λ μ λ = 107 / 42 . Thus, Q = 107 and the upper bound for the range of ϕ is calculated as 2 π ( n Λ + 1 http://www .cs.tut.fi/sgn/imaging/sparse/ Q ) /  λ μ λ  272 rad s . The number of experiments for each wa velength is S =6 . The peak-value of the Gaussian phase for the reference wa velength is ϕ = 242 radians, what is close to the upper limit 272 ra ds . T o make the problem more difficult we consider a very noisy scenario SN R =1 2 . 5 dB and N photon =7 . The achieved phase reconstructions are shown in Fig.1 as 3 D/ 2 D images. The proposed algorithm produced a nearly perfect reconstruction with RRM S E =0 . 00529 while coun- terpart reconstructions are completely failed for all wav e- lengths. These latter absolute phase reconstructions are ob- tained by the SP AR algorithm as the wrapped phase estimates unwrapped further by the PUMA algorithm [23] for each wa velength separately . The true wrapped phase pattern in this test is so complex and so irregular (see. Fig.2) that it is just impossible to unwrap it by any modern 2D algorithms. It is not surprising that PUMA, which is one of the best algorithms in the field, failed. The proposed MF-APR algorithm is able to resolve the problem ev en for the noisy data. Fig. 1. Gaussian absolute phase reconstructions, 3D/2D images: left-top, by the proposed MF-APR algorithms with the best RRMSE value, others are absolute phases obtained from the wrapped phases giv en by the SP AR algorithm and unwrapped by the PUMA algorithm. The phase reconstructions for the USAF test-object are s h o w ni nF i g . 3a s 3 D / 2 D phase images for noisy observations with SN R =8 . 74 dB and N photon =3 . 74 . The proposed algorithm demonstrates a quite accurate reconstruction with RRM S E =0 . 063 . The counterparts produced as above for each wa velength separately are failed. Fig. 2. The 3 D true absolute phase of the Gaussian phase object, and the corresponding 2 D wrapped phase, the wa velength λ ref = λ 1 . The peak-value of the phase for this test with the reference wa velength λ  = λ 1 is ϕ =3 0 rads. It is not so high as is for the Gaussian phase image but this test images is discontinuous binary with noisy observations what defines the complexity of the problem at hand. The synthetic wa velength approach discussed in Section I-A is not effecti ve for the considered tests as the wrapped phase differences are not small enough in order to ( μ λ − μ λ  ) ϕ ∈ [ π , − π ) . Thus, it cannot be a v aluable alternati ve to the proposed algorithm. For our experiments we use MA T - LAB R2015b and the computer with the processor Intel(R) Core(TM) i7-4800MQ@ 2.7 GHz. The complexity of the algorithm is characterized by the computational time. For 100 × 100 and 612 × 612 images, the computation time is 3 . 3 and 14 . 6 sec. per iteration. The APR algorithm takes, respectiv ely , 0 . 2 sec. and 8 . 5 sec. of this time. IV . C ONCLUSION The multi-frequency absolute phase retrieval from intensity observations is considered. This paper introduces a variational approach to object phase and amplitude reconstruction from noisy intensity observations. The maximum likelihood crite- rion used in the dev eloped optimization approach defines the intention to reach statistically optimal solutions. The phase re- triev al is an ill-possed in verse problem, where the observation noise is amplified and transferred to estimates of phase and amplitude. The sparsity dev eloped for modeling of varying variables is one of the key instruments for regularization of this in verse problem. The block-wise sparsity implemented by BM3D filters applied to amplitude/phase of complex-v alued variables and the real-valued object absolute phase. One of the key original elements of the algorithm is the weighted least square solution for aggregation of the multi-frequency estimates developed for both the absolute phase reconstruction from the observed multi-frequency wrapped phases and the phase-shifts compensation. The ef fectiv eness of the dev eloped algorithm is demonstrated by simulation experiments for the coded diffraction pattern scenario and very noisy Poissonian observations. The MA TLAB demo-codes of the algorithm are publicly av ailable. V. A CKNOWLEDGMENT This work is supported by Academy of Finland, project no. 287150, 2015-2019. R EFERENCES [1] A. W alther , “The question of phase retrieval in optics, ” Journal of Modern Optics, 10, no. 1, pp. 41-49, (1963). [2] R. P . Millane, “Phase retrieval in crystallography and optics, ” JOSA A 7 , no. 3, pp. 394-411 (1990). [3] J. C. Dainty and J. R. Fienup, “Phase retrie val and image reconstructio n for astronomy , ” Image Recovery: Theory and Application, pp. 231-275 (1987). [4] R. W . Gerchberg and W . O. Saxton, A practical algorithm for the determination of phase from image and diffraction plane pictures, Optik , vol. 35, pp. 237–246 (1972). [5] J. R. Fienup, Phase retrieval algorithms: a comparison, Appl. Opt ., vol. 21, pp. 2758–2769 (1982). [6] C. Guo, S. Liu, Shi, J. T . Sheridan, ”Iterativ e phase retrieval algorit hms. I: Optimization, ” Appl. Opt, vol. 54, no.15, pp. 4698-4708 (2015). Fig. 3. USAF phase reconstructions, 3D/2D images: left-top, by the proposed MF-APR algorithms with the best RRMSE value, others are absolute phases obtained from the wrapped phases given by the SP AR algorithm and unwrapped by the PUMA algorithm. [7] Y . Shechtman, Y . C. Eldar, O. Cohen, H. N. Chapman, J. Miao, and M. Segev , ”Phase retriev al with application to optical imaging: a contemporary overview , ” IEEE Signal Processing Magazine , pp. 87-109 (2015). [8] E. Pauwels, A. Beck, Y . C. Eldar , S. Sabach, ”On Fienup methods for regularized phase retrieval, ” arXiv:1702.08339v1, 2017. [9] V . Katkovnik and J. Astola, “ High-accuracy wav efield reconstruction: decoupled inverse imaging with sparse modeling of phase and amplitude, ” J. Opt. Soc. Am. A, vol. 29, pp. 44 – 54, 2012. [10] V . Katkovnik and J. Astola, ”Sparse ptychographical coherent diffra cti ve imaging from noisy measurements, ” J. Opt. Soc. Am. A , vol. 30, pp. 367-379, 2013. [11] V . Katkovnik, ”Phase retriev al from noisy data based on sparse approx- imation of object phase and amplitude, ” 2017, [12] V . Katkovnik and K. Egiazarian, ”Sparse super-resolution phase retr ie val from phase-coded noisy intensity patterns, ” Optical Engineering 56(9), 094103 (2017). [13] C.E. T owers , D.P . T owers, J.D.C. Jones, ”Absolute fringe order calcu la- tion using optimized multi-frequency selection in full-field profilometr y,” Optics and Lasers in Engineering 43, pp. 788–800 (2005). [14] P . Bao, G. Situ, G. Pedrini, and W . Osten, “Lensless phase microscopy using phase retrieval with multiple illumination wa velengths, ” Appl. Opt . 51(22), pp. 5486–5494 (2012). [15] C. Falldorf, S. Huferath-von Luepke, C. von Kopylo w , and R. B. Bergmann , ”Reduction of speckle noise in multiwav elength contouring, ” Appl. Opt. , vol. 51, Issue 34, pp. 8211-8215 (2012). [16] V . Katkovnik and J. Bioucas-Dias, "Multi-frequency Phase Unwrapping from Noisy Data: Adapti ve Least Squares Approach”, in Proceedings of the AIP Conference, vol. 1236, pp. 471-, Locarno, Switzerland (2010). [17] J. Bioucas-Dias and G. V aladão, "Multifrequency absolute phase esti- mation via graph cuts”, 17th European Signal Processing Conference, EUSIPCO’2009, Glasgow (2009). [18] K. Falaggis, D. P . T owers, and C. E. T o wers, Algebraic solution for phase unwrapping problems in multiwavelength interferometry , Applied Optics , vol. 53, No. 17, pp. 3737-3747 (2014). [19] W . W ang, X. Li, W . W ang, and Xiang-Gen Xia, ”Maximum likelihood estimation based robust Chinese Remainder Theorem for real numbers and its fast algorithm, ” IEEE T ransactions on Signal Processing , vol. 63, no. 13 (2015). [20] X. Li, Xiang-Gen Xia, and H. Huo, "T o wards robustness in residue number systems." IEEE Tr ansactions on Signal Processing 65.6, pp. 1497-1510 (2017). [21] J. W . Goodman, Introduction to Fourier optics, Roberts and Company Publishers (2005). [22] Y . Chen and E. J. Candès, ”Solving random quadratic systems of equations is nearly as easy as solving linear systems, ” Communications on Pur e and Applied Mathematics 70(5), pp. 822–883 (2017). [23] J. M. Bioucas-Dias and G. V aladão, " Phase unwrapping via graph cuts," IEEE T rans. Image Process. , vol. 16, no. 3, pp. 698–709, (2007).