2.1.

Multiwavelength Object and Image Modeling

Let $h_o(x)$, $x\in 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,\lambda }=b_{o,\lambda }\exp (j{\phi }_{o,\lambda })$, $\lambda \in \mathrm{\Lambda }$, where $\mathrm{\Lambda }$ is a set of the wavelengths and ${\phi }_{o,\lambda }=\frac{2\pi }{\lambda }{h}_o(n_{\lambda }-n_{\mathrm{env}})$ are phase delays corresponding to $h_o$, where $n_{\lambda }$ is the refractive index of the object and $n_{\mathrm{env}}$ is the refractive index of the environment, for the air $n_{\mathrm{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 ${\phi }_{o,\lambda }$ take values beyond the basic phase interval $[-\pi ,\pi ]$, then only wrapped phases can be obtained from $u_{o,\lambda }$ for separate $\lambda$.

In what follows, we apply dimensionless relative frequencies ${\mu }_{\lambda }=\frac{{\lambda }^{\prime }(n_{\lambda }-1)}{\lambda (n_{{\lambda }^{\prime }}-1)}$, which replace common notation for object wavefront as

Here ${\lambda }^{\prime }\in \mathrm{\Lambda }$ is a reference wavelength and ${\phi }_o$ is an absolute phase corresponding to this wavelength. The parameter ${\mu }_{\lambda }$ establishes a link between the absolute phase ${\phi }_o$ and wrapped phase ${\psi }_{o,\lambda }$ of $u_{o,\lambda }$, which can be measured at the $\lambda$-channel. The wrapped phase is related with the true absolute phase, ${\phi }_o$, as ${\mu }_{\lambda }{\phi }_o={\psi }_{o,\lambda }+2\pi k_{\lambda }$, where $k_{\lambda }$ is an integer, ${\psi }_{o,\lambda }\in [-\pi ,\pi )$. The link between the absolute and wrapped phase conventionally is installed by the wrapping operator $\mathcal{W}(·)$ as follows:

The image and observation modeling are defined as

The considered multiwavelength phase retrieval problem consists of reconstruction of ${\phi }_o$ and $b_{o,\lambda }$ from the observation $z_{s,\lambda }$ provided that ${\mu }_{\lambda }$, ${\mathcal{P}}_{s,\lambda ,d}\{\}$, and $\mathcal{G}\{\}$ are known. If the phases ${\mu }_{\lambda }{\phi }_o\in [\pi ,-\pi )$, the problem is trivial as estimates of ${\mu }_{\lambda }{\phi }_o$ and $b_{o,\lambda }$ can be found processing data separately for each $\lambda$ by any phase retrieval algorithm applicable for complex-valued object, in particular by the SPAR algorithm.²⁹ However, if the phases ${\mu }_{\lambda }{\phi }_o$ go beyond the range $[\pi ,-\pi )$, 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 ${\phi }_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,²⁶ 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 $\{z_{s,\lambda }\}$ are used jointly for reconstruction of ${\phi }_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.

2.2.

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,\lambda }(x)$ of the mean value $y_{s,\lambda }(x)$ takes a given non-negative integer value $k$, is given as

Another useful noise characteristic is the mean value of photons per pixel. It is defined as $N_{\text{photon},\lambda }=\sum _x{z}_{s,\lambda }(x)/N_{\text{sensor}}$, where $N_{\text{sensor}}$ is the number of sensor pixels. As in case of SNR, smaller values of $\chi$ lead to smaller $N_{\mathrm{photon}}$, i.e., to noisier observations $z_{s,\lambda }(x)$.

3.1.

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_{o,\lambda }^1$ is obtained from the observations $\{z_{s,\lambda }\}$ by the SPAR algorithm²⁸ separately for each wavelength. The main iterations start from the forward propagation (step 1) and follow by the amplitude update for $u_{s,\lambda }^t$ at step 2. The operator ${\mathrm{\varphi }}_1$ derived in Ref. 28 is defined as $u_{s,\lambda }^{t+1/2}={\mathrm{\varphi }}_1(u_{s,\lambda }^t,z_{s,\lambda })$, that means $u_{s,\lambda }^{t+1/2}=w·\frac{u_{s,\lambda }^t}{|u_{s,\lambda }^t|}$. Here the ratio $\frac{u_{s,\lambda }^t}{|u_{s,\lambda }^t|}$ denotes that the variables $u_{s,\lambda }^t$ and $u_{s,\lambda }^{t+1/2}$ have identical phases. The amplitude $w$ is calculated as

The backpropagation is realized in step 3, and the operator ${\mathrm{\varphi }}_2$ is defined by Eq. (9). The absolute phase reconstruction from the wrapped phases of $u_{o,\lambda }^{t+1}$ is produced in step 4 by the APR algorithm.²⁶ The obtained amplitude and phase update $u_{o,\lambda }^{t+1}$ at step 5.

The steps 3 and 4 are completed by the block-matching 3-D (BM3D) filtering.²⁹ In step 3, it is the filtering of complex-valued $u_{o,\lambda }^{t+1/2}$ produced separately for the wrapped phase and amplitude of $u_{o,\lambda }^{t+1/2}$. In step 4, this filtering is applied to the absolute phase ${\phi }_o^{t+1/2}$. 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.³¹ 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²⁹ being modified for various problems demonstrates the state-of-the-art performance.³² The details of BM3D as an advanced image filter can be seen in Ref. 29.

In step 3 of the proposed algorithm, $u_{o,\lambda }^{t+1}=\mathrm{BM}3\mathrm{D}(u_{o,\lambda }^{t+1/2})$, the BM3D is applied to the complex-valued variables $u_{o,\lambda }^{t+1/2}$. It is implemented in this paper as independent filtering of amplitude and wrapped phase

In step 4 of the proposed algorithm, ${\phi }_o^{t+1}=\mathrm{BM}3\mathrm{D}({\phi }_o^{t+1/2})$, the BM3D is applied for filtering of the real-valued variable ${\phi }_o$. In presented experiments, the parameters of the algorithm are fixed for all tests: ${\gamma }_1=1/\chi$, where $\chi$ is the parameter of the Poissonian distribution, ${\gamma }_1/{\gamma }_2=0.2$. The parameters of BM3D filters can be seen in Ref. 28.

The algorithm presented in Table 1 differs from the WM algorithm²⁶ 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 a single wavelength only, while in WM we obtain the observations for all wavelength simultaneously.

4.1.

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 $\mathrm{\Lambda }=[417,\hspace{0.17em}532,\hspace{0.17em}633]\mathrm{nm}$, with the corresponding refractive indexes [1.528, 1.519, 1.515] as taken for BK7 optical glass. The reference wavelength ${\lambda }^{\prime }=417\mathrm{nm}$ then the relative frequencies take values ${\mu }_{\lambda }=[1,\hspace{0.17em}0.7705,\hspace{0.17em}0.6425]$. The pixel sizes of the CMOS camera and SLM are 1.4 and $5.6\mu \mathrm{m}$, respectively. The distance $d$ between the object and CMOS camera is equal to 5 mm.

Fig. 1

Optical setup. R, G, B lasers,– red, green, and blue light sources; $L_1,L_2$, lenses; SLM, spatial light modulator; CMOS, registration camera.

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).³³ 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³⁴

The intensities of the light beams registered on the sensor are calculated as $z_{s,\lambda }=\mathcal{G}\{|\mathcal{A}{\mathcal{S}}_{\lambda ,d}\{M_s\circ u_{o,\lambda }\}{|}^2\}$, $s=1,\dots ,S$, $\lambda \subset \mathrm{\Lambda }$. Here $\mathcal{A}{\mathcal{S}}_{\lambda ,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 ${\mathcal{P}}_{s,\lambda ,d}$ in Eq. (3) is implemented as a combination of the AS propagation $\mathcal{A}{\mathcal{S}}_{\lambda ,d}$ and the modulation phase mask $M_s$.

4.2.

Reconstruction Results

The illustrating reconstructions are presented for two phase objects with the invariant amplitudes equal to 1 and the phases: Gaussian distribution ($100\times 100$) and US Air Force (USAF) resolution test-target ($64\times 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\pi$ rad, what corresponds to about 30 reference wavelengths ${\lambda }^{\prime }$ in variations of the profile $h_o$.

Fig. 2

Wrapped and absolute phases of the investigated objects Gauss and USAF, the reference wavelength ${\lambda }^{\prime }={\lambda }_1$.

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_{\text{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 $\chi$. The RRMSE curves for Gaussian and USAF phases demonstrate a similar behavior and go down for growing $\chi$ numbers, but RRMSE curve for the USAF object goes down more sharply and takes values smaller than 0.1 at $\chi =20$, while RRMSE curve for Gaussian phase takes a close value only at $\chi =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_{\mathrm{photon}}=0.75$ and 1.87, respectively.

Fig. 3

RRMSE and SNR as function of the parameter $\chi$ of the Poissonian distribution: solid triangles and dashed diamonds curves (blue in color images) show RRMSE for the Gaussian and USAF objects, respectively (left $y$-axis); solid circle curve (orange) shows SNR of the observations (right $y$-axis).

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 $\mathrm{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 $\mathrm{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.

Fig. 4

Gaussian object, RRMSE as function of SNR and the experiment number $S$: (a) after 50 and (b) 100 iterations.

In comparison with the WM algorithm²⁶ 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.¹⁸ The conditions of the experiments are: $\mathrm{SNR}=6.5\mathrm{dB}$, $N_{\text{photon}}=1.87$ for the Gaussian object and $\mathrm{SNR}=3.8\mathrm{dB}$, $N_{\text{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.

Fig. 5

Gaussian phase reconstructions $\mathrm{RRMSE}=0.0086$ for $\mathrm{SNR}=6.5$ dB ($N_{\mathrm{photon}}=1.87$), (a) 3-D surfaces and (b) 2-D absolute phases. From left to right: the WD algorithm, and the single wavelength reconstructions for ${\lambda }_1$, ${\lambda }_2$, and ${\lambda }_3$, respectively.

Fig. 6

USAF phase reconstructions $\mathrm{RRMSE}=0.030$ for $\mathrm{SNR}=3.8$ dB ($N_{\mathrm{photon}}=0.75$), (a) 3-D surfaces and (b) 2-D absolute phases. From left to right: the WD algorithm, and the single wavelength reconstructions for ${\lambda }_1$, ${\lambda }_2$, and ${\lambda }_3$, respectively.

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\pi$, 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\times 100\text{pixels}$ images, zero-padded to $1700\times 1700$, this time equals to 12 s.