## 3 Algorithm Development  
We consider the problem of wavefront reconstruction as an estimation of \(u_{o}\)\(\in\mathbb{C}^{n}\) from noisy observations \(\{\mathbf{z}_{s,\lambda}\}\). This problem is rather challenging mainly due the periodic nature of the likelihood function with respect to the phase \(\varphi_{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  
\[\mathcal{L}_{0}=\sum_{\lambda\in\Lambda}\sum_{s}\sum_{x\in\Omega}l(z_{s, \lambda}(x),|u_{s,\lambda}(x)|^{2}), \tag{13}\]  
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\). For the Poissonian and Gaussian distributions we have, respectively, \(l(z,|u|^{2})=|u|^{2}\chi-z\log(|u|^{2}\chi)\), and \(l(z,|u|^{2})=\frac{1}{2\sigma^{2}}||u|^{2}-z|^{2}\), \(\sigma^{2}\) stands for the noise variance. In Equation (13) s and \(x\) denote the experiment number and pixels of variables. 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:  
\[\mathcal{L}(u_{s,c},u_{o,\lambda},b_{o},\varphi,\delta_{\lambda})=\sum_{s} \sum_{c\in(R,G,B)}\sum_{x\in\Omega_{c}}l(z_{s,c}(x),|u_{s,c}(x)|^{2})+ \tag{14}\]  
\[\frac{1}{\gamma_{1}}\sum_{s}\sum_{c\in(R,G,B)}||u_{s,c}-\mathcal{P}_{s,c}u_{o, c}||_{2}^{2}+\frac{1}{\gamma_{2}}\sum_{\lambda}||b_{o,\lambda}\exp(j(\mu_{ \lambda}\cdot\varphi+\delta_{\lambda}))-u_{o,\lambda}||_{2}^{2},\]  
where \(\gamma_{1},\gamma_{2}>0\) are parameters, and \(||\cdot||_{2}^{2}\) is the Hadamard norm. Remind, that phase retrieval algorithms may reconstruct the phase \(\mu_{\lambda}\varphi_{o}\) only within invariant additive errors. These unknown phase shifts are modelled by \(\delta_{\lambda}\) as disturbance parameters. We say that the object complex exponents of different frequencies are in-phase (or synchronized) if \(\delta_{\lambda}=0\), \(\lambda\in\Lambda\) and de-phased if some of \(\delta_{\lambda}\
eq 0\). In the synchronized case the complex exponents of the object become equal to 1 provided \(\varphi_{o}\equiv 0\). For proper estimation of \(\varphi_{o}\) these phase shifts \(\delta_{\lambda}\) 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\in\Omega_{c}\) is produced independently for each \(c\), what separates the observations of different wavelengths. The sets \(\Omega_{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 \(\mathcal{P}_{s,c}u_{o,c}\). The phase retrieval is a reconstruction of \(\varphi_{o}\) and the amplitude of \(u_{o,\lambda}\) from the observations \(z_{s,c}\). The algorithm is designed using minimization of \(\mathcal{L}\) with respect to \(u_{s,c}\), \(u_{o,\lambda}\), \(b_{o}\), \(\varphi\) and the phase-shifts \(\delta_{\lambda}\). The following solutions of the optimization problem are used in the developed algorithm. 1. Minimization of \(\mathcal{L}\) with respect to \(u_{s,c}\) leads to the solution \[u_{s,c} = \left\{ \begin{array}{ll}\theta_{s,c}\cdot\frac{\tilde{u}_{s,c} }{|\tilde{u}_{s,c}|},&x\in\Omega_{c},\\ \tilde{u}_{s,c},&x\
otin\Omega_{c},\end{array} \right.\] (15) \[\tilde{u}_{s,c} = \mathcal{P}_{s,c}u_{o,c}.\] (16)Here the ratio \(\dfrac{\tilde{u}_{s,c}}{|\tilde{u}_{s,c}|}\) means that the variables \(u_{s,c}\) and \(\tilde{u}_{s,c}\) have identical phases.