## 1.

## Introduction

In modern science and technology, phase and wavefield imaging are popular and well-established techniques for high-accuracy measuring, recording, and reconstructing of two- (2-D) and three-dimensional (3-D) objects. The areas of applications are varying from astronomy and engineering to medicine and biology.^{1}^{,}^{2} In engineering, phase and wavefield sensing methods serve for nondestructive testing/control and precise measurements (e.g., see Refs. 3 and 4). In medicine and biology, phase measurements are exploited in microscopy and coherent tomography.

Phase imaging is a unique instrument to study details of internal structure of transparent or semitransparent specimens. While only intensity of light fields can be measured, visualization of phase from intensity observations is an important issue. In phase contrast microscopy, methods of wavefront modulation in the Fourier plane have been developed to resolve the visualization problem (Frits Zernike 1930s, Nobel prize 1953). Despite the revolutionary success of these methods, only qualitative visualization of phase can be achieved in this way, where features of specimens even visible maybe so distorted that accurate measurements and even proper interpretations can be problematic.

Quantitative visualization is targeted on direct phase imaging and precise phase measuring. Roughly speaking, there are two ways to achieve this goal. The first one is holography with measurements given as intensities of the sums of reference and object beams. The second one is phase retrieval, treated as an inverse diffusion imaging, essential alternative to holography, which does not require a reference beam. In modern science and technology, the quantitative phase imaging techniques are fundamentally based on digital data processing.

Let us start from the following general formalization of the phase retrieval problem:

## (1)

$${y}_{s}={|{\mathcal{P}}_{s}\{{u}_{o}\}|}^{2},\phantom{\rule[-0.0ex]{2em}{0.0ex}}s=1,\dots ,L,$$$L$ experiments are assumed in Eq. (1), where $s$ indicates the result for each of them. Equation (1) defines relations between the complex-valued wavefronts at the object plane and the power of the wavefronts at the sensor plane. It is convenient to also introduce a notation for the complex-valued wavefront at the sensor plane

The observations for the noiseless case corresponding to Eq. (1) are of the form

and for noisy observations where $\mathcal{G}$ stands for a generator of random observations.In this paper, we assume that the observations have a Poissonian distribution typical for optics with observations based on photon counting.

Reconstruction of the complex-valued object ${u}_{o}$ (phase and amplitude) from noiseless or noisy observations is phase retrieval problem. Here, phase emphasizes that the object phase is a variable of the first priority while the object amplitude is treated as an auxiliary variable often useful only in order to improve phase imaging.

Note that the term phase retrieval is originated from the following mathematical problem. Let ${u}_{s}$ be the Fourier transform (FT) of ${u}_{o}$, ${u}_{s}=\mathcal{F}({u}_{o})$. If ${u}_{s}$ is given, ${u}_{o}$ can be precisely calculated as the inverse FT, ${u}_{o}={\mathcal{F}}^{-1}({u}_{s})$. Now let the absolute value $|{u}_{s}|$ of ${u}_{s}$ be given and the phase of ${u}_{s}$ is unknown. Is it possible to reconstruct the phase of ${u}_{s}$ and in this way the original ${u}_{o}$ from the amplitude of FT $|{u}_{s}|$? In general, the answer is negative and positive only for the special classes of ${u}_{o}$, in particular, for the so-called minimum phase signals, or provided some restrictions. In this phase retrieval formulation, the term phase points on the phase of the FT, ${u}_{s}$.

In optics, the priority is different. The phase to be retrieved is phase of the object ${u}_{o}$. The image formation operators ${\mathcal{P}}_{s}$ in Eq. (1) depend on optical setups. Various methods are developed in order to make these ${\mathcal{P}}_{s}$ sufficiently different for different $s$ and gain observation diversity, enabling finding ${u}_{o}$ from observations ${\{{y}_{s}\}}_{1}^{L}$. Defocusing of the registered images is one of the popular instruments to get a sufficient phase diversity.^{5}6.7.^{–}^{8} In a recent development, a spatial light modulator (SLM) is exploited for defocusing (e.g., see Refs. 9 and 10). The $4f$-optical configuration with SLM in the Fourier plane for defocus imitation is proposed in Ref. 11 and further studied in Ref. 12.

Random phase modulation of the wavefront is another tool to achieve a desirable phase diversity. It results in observations known as coded diffraction patterns (e.g., see Refs. 1314.–15)

## (5)

$${y}_{s}={|\mathcal{P}\{{\mathcal{M}}_{s}\xb7{u}_{o}\}|}^{2},\phantom{\rule[-0.0ex]{2em}{0.0ex}}s=1,\dots ,L,$$The phases ${\varphi}_{k,l}(s)$ in ${\mathcal{M}}_{s}$ can be generated as deterministic or random. The phase modulation is able to dramatically change the diffraction pattern of $\mathcal{P}\{{u}_{o}\}$ by redistribution of the observed intensities from low to high frequencies.

## 1.1.

### Phase Retrieval Algorithms

There is a growing flow of publications on phase retrieval. Various versions of the Gerchberg–Saxton (GS) techniques are quite universal and applicable for different setups (e.g., see Refs. 7, 8, 1617.–18). These algorithms based on alternating projections between the object and observations planes allow to incorporate any information available for the variables in these planes. Recent developments in this area as well as a review can be seen in Ref. 19.

Contrary to this type of the intuitive heuristic algorithm, the variational formulations have a strong mathematical background and lead to algorithms solving optimization problems. In particular, in Ref. 20 one can find constrains sufficient for uniqueness of the solution and algorithms that are very different from GS, such as the semidefinite programming phase lifting using matrix completion (PhaseLift algorithm)^{21} and the greedy sparse phase retrieval (GESPAR algorithm).^{22}

There are many publications on revisions of the GS techniques using optimization formulations. In particular, the links between the conventional GS and variational formulations are studied in Refs. 23 and 24.

Concerning variational formulations for algorithm design, we wish to note the recent Wirtinger flow and truncated Wirtinger flow (TWF) algorithms.^{14}^{,}^{25} Methodologically, these algorithms are developed for the Poissonian likelihood criterion, i.e., for Poissonian noisy observations. Simulation experiments confirm that these algorithm works precisely provided nearly noiseless observations. However, they are not so efficient for noisy observations.^{26}

Phase retrieval from coded diffraction patterns of the type (5) is of special interest in the recent publications (e.g., see Refs. 14 and 27). The uniqueness of solution for this scenario is proved in the later paper.

A new variational algorithm for phase retrieval from noisy data based on transform domain sparsity for the object phase and amplitude is developed in Ref. 26. Simulation experiments demonstrate that this algorithm enables the accuracy identical to the accuracy of the TWF algorithm for noiseless data and the principal advantage for noisy data. The sparsity concept as a tool for phase retrieval is a topic of the paper,^{28} where an original sparsifying learning transform is developed.

The spatial resolution in phase retrieval is limited by two principal factors: low-pass filtering by the propagation operator $\mathcal{P}$ and by pixel size in the pixelated discrete sensor and modulation phase masks. Due to these factors, high-frequency spatial information is lost in intensity observations, which can be treated as observations of the subsampled true object ${u}_{o}$. Various methods for superresolution imaging allowing compensation of these subsampling effects are of special interest. One of the straightforward approaches to overcome the pixel size limitations is to use a sequence of laterally shifted holograms (e.g., see Refs. 2930.–31). Compressed sensing (CS) or sparse imaging is a computational approach for the restoration of subsampled data based on a special mathematical modeling of ${u}_{o}$. Applications of this sort of the technique in optics can be seen in Refs. 3233.34.35.–36.

Other factors limiting the spatial resolution concern observation errors. First, we need to mention that a Poissonian noise appears due to the measurement process in optics counting the photons hitting the sensor. Second, the use of a digital camera introduces the readout noise usually modeled by a Gaussian distribution and quantization errors. The later ones can be modeled as a uniform distribution random variables. The quantization effects for phase retrieval are studied,^{37} and it is shown that a low-bit quantization may seriously diminish the accuracy of the phase retrieval.

## 1.2.

### Contribution and Structure of This Paper

In this paper, we consider phase retrieval from Poissonian noisy phase coded diffraction patterns with the optical setup shown in Fig. 1. The complex-valued object ${u}_{o}$ to be reconstructed is placed against the SLM applied for phase modulation of the wavefront. A distance between the object and the digital sensor is equal to $2f$, where $f$ is a focal length of this lens, located in the middle between the object and lens. The system is illuminated by a uniform monochromatic coherent laser beam of wavelength $\lambda $ and intensity equal to 1. The forward propagation operator $\mathcal{P}$ in Eq. (5) is calculated as a rescaled FT, provided that the axial Fresnel approximation for the propagation is assumed. In this setup, the object, the SLM, and lens are considered as phase transformers of the wavefronts. The superresolution sparse phase and amplitude reconstruction (SR-SPAR) algorithm proposed in this paper is designed for superresolution phase/amplitude imaging, which is optimal for Poissonian observations. It is shown by computational experiments that high-accuracy superresolution reconstructions can be achieved with spatial resolution going up to a quarter of wavelength. The superresolution GS (SR-GS) algorithm is introduced as a simplified and faster version of SR-SPAR efficient for noiseless data. For the later case, both SR-GS and SR-SPAR demonstrate identical accuracy. SR-SPAR design is based on the methodology developed in Ref. 26 for pixel-resolution phase retrieval (SPAR algorithm).

This paper is organized as follows. In Secs. 1.3 and 1.4, image formation and noisy observation modeling are presented. The sparsity in the complex domain and the SR-SPAR and SR-GS algorithms are given in Sec. 2. Section 3 concerns the experimental study of the proposed algorithms.

## 1.3.

### Image Formation

For wavefront propagation from the object to sensor planes, we use the paraxial Fresnel modeling. It gives the following link between the object wavefront ${u}_{o}$ and the wavefront at the sensor plane ${u}_{s}$ [Ref. 38, Eq. (5.19)]:

## (6)

$${u}_{s}(\xi ,\eta )=\mu {\int}_{-\infty}^{\infty}{\int}_{-\infty}^{\infty}{u}_{o}(x,y){\mathcal{M}}_{s}(x,y)\mathrm{exp}[j\frac{-2\pi}{\lambda f}(x\xi +v\eta )]\mathrm{d}x\text{\hspace{0.17em}}\mathrm{d}y,$$Here ${u}_{o}(x,y)$ and ${u}_{s}(\xi ,\eta )$ are complex-valued distributions of the wavefronts at the object plane [lateral coordinates $(x,y)$] and the sensor plane [lateral coordinates $(\xi ,\eta )$]. $\mathcal{M}(x,y)$ is a complex-valued transmission function of SLM. Using FT, the input–output model (6) can be given in the form

## (7)

$${u}_{s}(\xi ,\eta )=\mu {\mathcal{F}}_{{u}_{o}\xb7{\mathcal{M}}_{s}}(\xi /\lambda f,\eta /\lambda f),$$Equations (6) to (7) define the forward propagation operator ${\mathcal{P}}_{s}$ in Eq. (2). Note that the model (6) as used in Eq. (5) is discrete-continuous with continuous ${u}_{o}$ and a physical discretization imposed on ${u}_{s}(\xi ,\eta )$ by the pixelated sensor and by the pixelated SLM on ${\mathcal{M}}_{s}$.

It is useful to mention that the forward propagation [Eq. (6)] is valid also for a single lens system, provided that the lens is located in the object/SLM plane and the distance between this plane and the sensor is equal to $f$.

## 1.4.

### Poissonian Observations

The measurement process in optics amounts to count the photons hitting the sensor. This process is well modeled by the independent Poisson random variables in the following form:

where $p({z}_{s}[l]=k)$ is the probability that the random observation ${z}_{s}[l]$ takes integer value $k\ge 0$ and ${y}_{s}[l]$ is the intensity of the wavefront at the pixel $l$ defined by Eq. (1).The parameter $\chi >0$ in Eq. (8) is the scaling factor of the Poisson distribution that can be interpreted as an exposure time. Recall that the mean and the variance of Poisson random variable ${z}_{s}$ are equal and are given by ${y}_{s}\chi $, i.e., $E\{{z}_{s}\}=\mathrm{var}\{{z}_{s}\}={y}_{s}\chi $. Defining the observation signal-to-noise ratio (SNR) as the ratio between the square of the mean and the variance of ${z}_{s}$, we have $\mathrm{SNR}={E}^{2}\{{z}_{s}\}/\mathrm{var}\{{z}_{s}\}={y}_{s}\chi $. It follows that the relative noisiness of observations becomes stronger as $\chi \to 0$ ($\mathrm{SNR}\to 0$) and approaches zero when $\chi \to \infty $ ($\mathrm{SNR}\to \infty $). The later case corresponds to the noiseless scenario, i.e., ${z}_{s}/\chi \to {y}_{s}$ with the probability 1.

The scale parameter $\chi $ is of particular importance for modeling as it allows to control the level of noise in the observations.

## 2.

## Superresolution Sparse Phase Retrieval

## 2.1.

### Sparse Wavefront Representations

The recent (10 to 15 years) sparse approximation and CS techniques state that a signal/image can be sampled at a rate much smaller than what is commonly prescribed by Shannon–Nyquist. The sampling of a signal can indeed be performed as a function of its “intrinsic dimension” rather than according to its cutoff frequency. The modern advanced techniques are able to design parametric regression models for small patches of images and combine them in high-quality reconstructions. What is considered in our paper is a development of this sort of thinking to the complex domain, in particular, for superresolution in phase retrieval.

Image sparsity follows from the commonly observed self-similarity of small fragments (patches) of images, meaning that similar features can be found in patches located in different parts of the image. It follows that an image may admit sparse representations: it can be well approximated by linear combinations of a small number of basic functions. Sparsity has been a hot topic in the last years for various imaging problems (e.g., see Ref. 39).

For the complex domain images, such as the object ${u}_{o}={B}_{o}\text{\hspace{0.17em}}\mathrm{exp}(i{\phi}_{o})$, sparse modeling can be presented in a number of ways essentially different from the methods standard for real-valued variables. This principal difference starts from the fact that ${u}_{o}$ is defined by two variables: phase ${\phi}_{o}$ and amplitude ${B}_{o}$.

The sparse representation can be imposed on complex-valued ${u}_{o}$ directly using complex-valued basic functions or on the following pairs of real-valued variables:

1. Phase $\phi $ (interferometric or absolute) and amplitude ${B}_{o}$;

2. Real and imaginary parts of ${u}_{o}$.

Remember that an interferometric (wrapped) phase is restricted to the interval $[-\pi ,\pi )$, whereas an absolute (unwrapped) phase is different by adding an integer number of $2\pi $ to the interferometric phase. In what follows, we denote the interferometric (wrapped) phase of the object as ${\phi}_{o}$ and the corresponding absolute phase as ${\phi}_{o,\mathrm{abs}}$. We introduce the phase-wrap operator $\mathcal{W}:\mathbb{R}\mapsto [-\pi ,\pi )$, linking the absolute and principal phases as ${\phi}_{o}=\mathcal{W}({\phi}_{o,\mathrm{abs}})$. We also define the unwrapped phase as ${\phi}_{o,\mathrm{abs}}={\mathcal{W}}^{-1}({\phi}_{o})$. Notice that ${\mathcal{W}}^{-1}$ is not an inverse operator for $\mathcal{W}$ because the latter is highly nonlinear, and for signals of dimension two and higher, there is no one-to-one relation between ${\phi}_{o,\mathrm{abs}}$ and ${\phi}_{o}$.

In principle, the absolute phase always can be reconstructed as the interferometric one with the application of unwrapping as postprocessing. However, for objects with phase varying beyond the interval $[-\pi ,\pi )$, the absolute phase sparse modeling brings essential advantage. It is because, wrapped phases are complicated by fringes, making images more difficult for sparse approximation.

The success of the sparsity approach depends on how rich and redundant are the used dictionaries/transforms (sets of basic functions). In this paper, the sparsity analysis and synthesis are based on the recent and proved to be very efficient, block-matching 3D (BM3D) denoising algorithm.^{40}

Let us mention the basic steps of this advanced technique known as nonlocal self-similarity sparsity. At the first stage, the image is partitioned into small overlapping square patches. For each patch, a group of similar patches is collected that are stacked together to form a 3-D array (group). This stage is called grouping. The entire 3-D group-array is projected onto a 3-D predefined transform basis. The spectral coefficients obtained as a result of this transform are hard-thresholded (small coefficients are zeroed) and the inverse 3-D 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 patches of the entire image and obtained overlapped filtered patches are aggregated in the final image estimate. This last stage is called aggregation. The details of this algorithm can be seen in Ref. 40.

The links of the BM3D algorithm with the general sparsity concept are revealed in Ref. 41, where it is shown that the grouping operations define the data adaptive analysis and synthesis image transforms (frames) and these transforms combined with the thresholding define the thresholding stage of the BM3D algorithm. It is emphasized that sparsity is achieved mainly due to the grouping, which allows the joint analysis of similar patches and, in this way, to guaranty the sparsity (self-similarity of patches), at least for each of the 3-D groups.

Note that the standard BM3D algorithm, as it is presented in Ref. 40, is composed of two successive stages: thresholding and Wiener filtering. In this paper, we use a simplified version of BM3D, as it is introduced in Ref. 41, including grouping, transforms, and thresholding without Wiener filtering.

In what follows, we exploit, for the complex-valued ${u}_{o}$, the sparsity imposed on phase and amplitude. The variational formulation for reconstruction of a complex-valued ${u}_{o}$, optimal for noisy data, results in the likelihood type criteria and optimization with the constrained imposed on sparsity. It has been shown for various optical problems (e.g., see Refs. 4243.–44) that the algorithms are iterative and the sparsity appeared as BM3D filtering applied separately to estimates of phase and amplitude.

This filtering can be represented in the form

where the filters are applied to the phase and amplitude of ${u}_{o}$.Here ${\widehat{\phi}}_{o}$ and ${\widehat{B}}_{o}$ are sparse approximations of ${\phi}_{o}$ and ${B}_{o}$; phase and amplitude as indices of BM3D are used in order to emphasize that the parameters of BM3D can be different for phase and amplitude; and $t{h}_{\phi}$ and $t{h}_{B}$ are threshold parameters of the algorithms. The phase in Eq. (9) can be interferometric or absolute depending on the sparsity formulation.

The implementation of the sparsity hypothesis in the form of the special filters Eqs. (9)–(10) is in line with the recent concept plug-and-play,^{45}46.^{–}^{47} stating that any efficient filter can serve as a potentially good prior and efficient regularizator in variational design of data processing algorithms.

## 2.2.

### Superresolution SPAR Algorithm

## 2.2.1.

#### Discretization

The computational wavefront restoration is going from the continuous domain wavefront propagation [Eq. (7)] to the corresponding discrete model based on pixelation of the object ${u}_{o}$, thus, we arrive to the discrete modeling of the system linking discrete values of the sensor output (observations) with the discrete values of the object ${u}_{o}$.

Conventionally, the pixels are square of the size ${\mathrm{\Delta}}_{\mathrm{SLM}}\times {\mathrm{\Delta}}_{\mathrm{SLM}}$ and ${\mathrm{\Delta}}_{s}\times {\mathrm{\Delta}}_{s}$ for the SLM and sensor, respectively. In what follows, for simplicity, ${\mathrm{\Delta}}_{\mathrm{SLM}}={\mathrm{\Delta}}_{s}$. A continuous object is discretized by pixels of the size ${\mathrm{\Delta}}_{o}\times {\mathrm{\Delta}}_{o}$. This discretization is necessary both for digital data processing as well as for modeling of wavefront propagation and image formation. Contrary to the pixels of SLM and sensor defined by the corresponding optical-electronic devices, the object pixels are computational, which maybe be taken arbitrary small.

Assuming for a moment that ${\mathrm{\Delta}}_{o}={\mathrm{\Delta}}_{s}={\mathrm{\Delta}}_{\mathrm{SLM}}$, the reconstruction of ${u}_{o}$ from the observations $\{{z}_{s}\}$ is the standard phase retrieval problem with an object resolution dictated by the pixel size of the sensor and the SLM.

Let us term this case pixel-resolution imaging.

If ${\mathrm{\Delta}}_{o}<{\mathrm{\Delta}}_{s}$, we arrive to a much more challenging problem of pixel superresolution or subpixel resolution imaging. Further, if ${\mathrm{\Delta}}_{o}$ is so small that ${\mathrm{\Delta}}_{o}\simeq \lambda $, then it is wavelength resolution. Going further to ${\mathrm{\Delta}}_{o}<\lambda $ we arrive to subwavelength or wavelength superresolution. The superresolution phase retrieval with smaller and very small ${\mathrm{\Delta}}_{o}$ as compared with ${\mathrm{\Delta}}_{s}$ and $\lambda $ is the goal of this paper.

It is convenient to assume that ${\mathrm{\Delta}}_{s}={r}_{s}\xb7{\mathrm{\Delta}}_{o}$, where ${r}_{s}\ge 1$ is an integer pixel superresolution factor. In this case, the SLM pixel ${\mathrm{\Delta}}_{\mathrm{SLM}}\times {\mathrm{\Delta}}_{\mathrm{SLM}}$ covers ${r}_{s}^{2}$ computational object pixels and provides the same modulation phase-shift to all object pixels in this group.

Using for calculation the fast Fourier transform (FFT), we arrive to the discrete analog of Eq. (7), within an invariant factor $\mu $, in the form

## (11)

$${u}_{s}[k,l]=\mathrm{FFT}\{{u}_{o}[s,t]\xb7{\mathcal{M}}_{s}[s,t]\}\frac{{\mathrm{\Delta}}_{o}^{2}}{\lambda f},$$Here, the variables ${u}_{s}[k,l]$ and ${u}_{o}[s,t]\xb7{\mathcal{M}}_{s}[s,t]$ are sampled with the computational period ${\mathrm{\Delta}}_{o}$. Then, in particular, the modulation function ${\mathcal{M}}_{s}[s,t]$ is a piecewise invariant with ${r}_{s}\times {r}_{s}$ squares of invariant values covering the corresponding pixels of SLM.

The constraint (12) is typical for use of FFT for the calculation of discrete FT. All functions and FFT in Eq. (11) are calculated for the square support ${N}_{D}\times {N}_{D}$, where ${N}_{D}$ is always higher (even much higher) than the pixelated sizes of the object, the SLM, and the sensor.

According to Eq. (5), the discrete diffraction pattern is calculated as ${y}_{s}[k,l]={|{u}_{s}[k,l]|}^{2}$ with the noisy observations obtained according to the Poissonian distribution (8). Note that these computational ${y}_{s}[k,l]$ are given with the computational period ${\mathrm{\Delta}}_{o}$ while the observations are introduced with the sampling period ${\mathrm{\Delta}}_{s}$. Equations (11)–(12) define the discrete forward propagation model of the system shown in Fig. 1. In order to simplify the presentation, we preserve the notation $\mathcal{P}$ for this discrete model initially introduced for the continuous domain variables.

## 2.2.2.

#### SR-SPAR Algorithm

The presented SR-SPAR algorithm is derived from the variational formulation introduced for optimal reconstruction of ${u}_{o}$ from Poissonian observations $\{{z}_{s}[k,l]\}$. The corresponding minus log-likelihood for Poissonian observations according to Eq. (8) is as follows:

## (13)

$$\mathcal{L}=\sum _{s=1}^{L}\sum _{k,l}{[|{u}_{s}[k,l]|}^{2}\chi -{z}_{s}[k,l]\mathrm{log}({|{u}_{s}[k,l]|}^{2}\chi )].$$This criterion should be minimized with respect to ${u}_{o}[k,l]$, provided Eq. (11) linking ${u}_{o}$ and ${u}_{s}$ and restrictions imposed by the sparsity requirements.

The derivation of the algorithm is similar to the technique developed in Ref. 26 for the pixel-resolution phase retrieval. The difference mainly concerns the sampling rates: ${\mathrm{\Delta}}_{o}={\mathrm{\Delta}}_{s}$ in Ref. 26 and in this paper ${\mathrm{\Delta}}_{o}={\mathrm{\Delta}}_{s}/{r}_{s}$ meaning that the observations should be upsampled by a factor ${r}_{s}$.

We present the SR-SPAR algorithm in the form given in Table 1. It is emphasized that SR-SPAR, being based on the minimization of Eq. (13), is optimal, in the statistical sense, for Poissonian observations.

## Table 1

SR-SPAR phase retrieval algorithm.

$\text{Input}:\{{\tilde{z}}_{s}\}$; $s=1,\dots ,L,{x}^{1}$ | |

for $t=1,..,T$; | |

1. | Forward propagation: |

${v}_{s}^{t}=\mathcal{P}\{{\mathcal{M}}_{s}\xb7{x}^{t}\}$, ${v}_{s}^{t}\in {S}_{D}$, $s=1,\dots ,L$; | |

2. | Poissonian noise suppression: |

${u}_{s}^{t}=\{\begin{array}{c}{b}_{s}^{t}\text{\hspace{0.17em}}\mathrm{exp}[j\xb7\text{angle}({v}_{s}^{t})],{v}_{s}^{t}\in {S}_{S},\\ {v}_{s}^{t},{v}_{s}^{t}\in {S}_{D}\setminus {S}_{S},s=1,\dots ,L;\end{array}$ | |

3. | Backward propagation: |

${x}^{t}=\frac{1}{L}\sum _{s=1}^{L}{\mathcal{M}}_{s}^{*}\xb7{\mathcal{P}}^{-1}\{{u}_{s}^{t}\}$; | |

4. | Phase unwrapping: |

${\phi}_{\mathrm{abs}}^{t}={\mathcal{W}}^{-1}(\text{angle}({x}^{t}))$; | |

5. | Sparse phase and amplitude filtering: |

${\phi}_{\mathrm{abs}}^{t+1}={\mathrm{BM}3\mathrm{D}}_{\text{phase}}({\phi}_{\mathrm{abs}}^{t},t{h}_{\phi})$, | |

${B}^{t+1}={\mathrm{BM}3\mathrm{D}}_{\mathrm{ampl}}(\mathrm{abs}({x}^{t}),t{h}_{a})$; | |

6. | Object wavefront update: |

${x}^{t+1}={B}^{t+1}\text{\hspace{0.17em}}\mathrm{exp}(j{\phi}_{\mathrm{abs}}^{t+1})$; | |

$\text{Output}:{\widehat{\phi}}_{o,\mathrm{abs}}={\phi}_{\mathrm{abs}}^{T+1}$, ${\widehat{B}}_{o}={B}^{T+1}$. |

The inputs ${\tilde{z}}_{s}$ in this algorithm are the observations ${z}_{s}$ upsampled by a factor ${r}_{s}$. We use the zero-order upsampling giving ${\tilde{z}}_{s}$ as piecewise invariant function with the invariant values for computational pixels corresponding to each of the larger size pixels of the sensor. All calculations in the SR-SPAR algorithm are produced for high-resolution variables with the sampling ${\mathrm{\Delta}}_{o}$.

At step 1, the object wavefront estimate ${x}^{t}$ is multiplied by the phase mask ${\mathcal{M}}_{s}$ and propagated by the operator $\mathcal{P}$ to the sensor plane, with the result denoted as ${v}_{s}^{t}$. These wavefronts are calculated for the diffraction area ${N}_{D}\times {N}_{D}$ denoted as ${S}_{D}$.

At step 2, the wavefront is updated to the variable ${u}_{s}^{t}$ by filtering the amplitude of ${v}_{s}^{t}$ according to the given observations ${\tilde{z}}_{s}$. The following formula, as derived in Ref. 26, defines the rule on how the updated amplitude ${b}_{s}$ is calculated:

## (14)

$${b}_{s}=\frac{|{v}_{s}|+\sqrt{{|{v}_{s}|}^{2}+4{\tilde{z}}_{s}{\gamma}_{1}(1+{\gamma}_{1}\chi )}}{2(1+{\gamma}_{1}\chi )}.$$These calculations are pixelwise; ${\gamma}_{1}>0$ is the parameter of the algorithm. This update is produced provided known observation ${\tilde{z}}_{s}$, i.e., for the pixels belonging to the sensor area ${S}_{S}$, ${v}_{s}\in {S}_{S}$. In our modeling, the computational diffraction area ${S}_{D}$ is always equal or larger than the sensor area, ${S}_{S}\subseteq {S}_{D}$. For the area out of the sensor, the wavefront values are preserved, ${u}_{s}={v}_{s}$ for ${v}_{s}\in {S}_{D}\setminus {S}_{S}$.

At step 3, the estimates $\{{u}_{s}^{t}\}$ backpropagate to the object plane and update the object wavefront estimate ${x}^{t+1}$. Here, ${\mathcal{M}}_{s}^{*}$ means a complex conjugate ${\mathcal{M}}_{s}$ and ${\mathcal{P}}^{-1}\{{u}_{s}^{t}\}={\mathrm{FFT}}^{-1}\{{u}_{s}^{t}[k,l]\}\lambda f/{\mathrm{\Delta}}_{o}^{2}$.

The sparsification (filtering on the base of sparse approximations) is produced in step 5. The unwrapping of the phase with reconstruction of the absolute phase in step 4 is necessary only if the range of the object phase goes beyond $2\pi $.

Following to Ref. 26, we introduce also a simplified version of SR-SPAR named the SR-GS. It differs in two points from SR-SPAR in Table 1: the phase unwrapping and BM3D filtering (steps 4 and 5) are omitted and the Poissonian filtering in step 2 is replaced by the rule ${b}_{s}\to \sqrt{{\tilde{z}}_{s}/\chi}$, which corresponds to the amplitude update standard for the GS style algorithms. This later rule follows from the optimal solution (14) provided that the data are noiseless, i.e., $\chi $ is very large.

It was demonstrated in Ref. 42 that using the update shown in step 2 and different for ${v}_{s}^{t}\in {S}_{D}\setminus {S}_{S}$ and ${v}_{s}^{t}\in {S}_{S}$ allows to improve the accuracy of the wavefront reconstruction. In Ref. 48, this effect is interpreted as a self-extrapolation of holograms applied for resolution enhancement.

We make publicly available the MATLAB^{®} democodes^{49} of the developed SR-GS and SR-SPAR algorithms, which can be used to reproduce the experiments presented in this paper as well as for further tests.

## 3.

## Numerical Experiments

Both algorithms, SR-SPAR and SR-GS, were tested for various models of ${u}_{o}$. In what follows, we are restricted mainly to $256\times 256$ phase-objects of invariant amplitude and three types of varying phase: test-images Lena normalized to the interval [0, $\pi /2$]; Gaussian shape absolute phase (phase range 50 radians), and discontinuous shear plane (phase range 65 rad). Respectively, we treat experiments with Lena as interferometric phase imaging, as they do not require phase unwrapping, and experiments with Gaussian and shear plane as absolute phase imaging requiring the unwrapping operation in SR-SPAR. The PUMA algorithm^{50} is used for phase unwrapping in step 4 of SR-SPAR.

The fixed parameters of the experiments are: ${\mathrm{\Delta}}_{s}={\mathrm{\Delta}}_{\mathrm{SLM}}=5.2\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$, $\lambda =632.8\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}$, sensor size, in pixels, $4096\times 4096$, computational diffraction area ${S}_{D}$ of size $5120\times 5120$. The main varying parameters are the computational sampling period ${\mathrm{\Delta}}_{o}$, the Poissonian noise parameter $\chi $, and computational resolution factor calculated with respect to the sensor as ${r}_{s}={\mathrm{\Delta}}_{s}/{\mathrm{\Delta}}_{o}$. It is assumed that ${r}_{s}$ is integer and takes values ${r}_{s}=1$, 8, 16, and 32. The ${r}_{s}=1$ corresponds to the pixel-resolution and larger ${r}_{s}$ mean superresolution of the higher order. The sensor of a much larger size than the object is taken in order to enable a good quality of superresolution with large values of the superresolution factors ${r}_{s}$.

It is natural to also measure the superresolution with respect to the wavelength $\lambda $ as the ratio ${r}_{\lambda}={\mathrm{\Delta}}_{o}/\lambda ={\mathrm{\Delta}}_{s}/({r}_{s}\lambda )$. Then, ${r}_{s}=1$ gives the pixel-resolution with ${r}_{\lambda}=8.21$, i.e., the sensor and SLM pixels are about eight times larger than the wavelength $\lambda $. For ${r}_{s}=8$, we obtain the wavelength resolution with ${r}_{\lambda}=1.03$, the higher values ${r}_{s}=16$, 32 correspond to the subwavelength resolution with computational pixels ${\mathrm{\Delta}}_{o}$ smaller than the wavelength with the wavelength resolution factors ${r}_{\lambda}=0.515$ and 0.257.

Note that according to the restriction (11), smaller ${\mathrm{\Delta}}_{o}$ (larger ${r}_{s}$) assumes that the lens with a smaller focal distance should be used in the considered optical setup. For the introduced set ${r}_{s}$, we obtain the following focal distances $f=[54.7,6.8,3.4,1.7]\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$, respectively.

In our experiments, the phase modulation masks ${\mathcal{M}}_{s}(k,l)=\mathrm{exp}[j{\varphi}_{k,l}(s)]$ are random with the Gaussian independent zero-mean phase values, ${\varphi}_{k,l}(s)\sim \mathcal{N}(0,\pi /4)$.

The accuracy of the wavefront reconstruction is characterized by RMSE criteria calculated independently for amplitude and phase. The object phase image can be estimated at least within an invariant global phase-shift ${\phi}_{\text{shift}}$. It is estimated using as reference the phase of the true object. This correction of the phase is done only for the calculation of the criteria and for result imaging and is not used in the algorithm iterations.

In what follows, we produce calculations for noisy and nearly noiseless data with the Poissonian scale parameter $\chi $ taking values in the internal [1, 1000]. The smallest $\chi $ results in the noisiest observations. The corresponding SNR is calculated in dB as

## (15)

$$\mathrm{SNR}=10\text{\hspace{0.17em}}{\mathrm{log}}_{10}({\chi}^{2}\sum _{s=1}^{L}{\Vert {y}_{s}\Vert}_{F}^{2}/\sum _{s=1}^{L}{\Vert {y}_{s}\chi -{z}_{s}\Vert}_{F}^{2})\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{dB}.$$For the superresolution experiments, we use the objects with a fixed number of computational pixels of size ${\mathrm{\Delta}}_{o}={\mathrm{\Delta}}_{s}/{r}_{s}$, thus larger ${r}_{s}$ means a smaller physical size of the object. The successful superresolution imaging, in particular the wavelength resolution, requires a sensor size being much larger than the object size.

## 3.1.

### Modulation Phase Mask and Sparsity

Let us start from qualitative observations concerning the effects of the basic ingredients of the considered setup and the developed algorithm. Figure 2 shows the reconstruction of the object phase for noiseless data ($\chi =1000$), provided that only a single experiment is produced, $L=1$. The left image shows the true phase. The next completely destroyed image in Fig. 2 is obtained from the experiment with no phase modulation and without the sparse modeling for phase and amplitude. Thus, the diffraction pattern is a squared amplitude of FT of the object complex exponent. The third image is obtained from a single experiment, where the phase modulation is employed and no BM3D filters are used (SR-GS algorithm). This modulation makes the main features of the phase distribution at least visible but quite noisy. This noise is a result of the used phase modulation. The fourth image is obtained by the SR-SPAR algorithm, i.e., with phase modulation and BM3D filtering for the phase and amplitude. It shows nearly perfect reconstruction of the object phase.

For a larger number of experiments ($L>1$), the accuracy of the phase reconstruction with the phase modulation improves quickly both for processing with SR-SPAR and with SR-GS. Figure 2 and above comments are given for the pixel-resolution imaging, ${r}_{s}=1$, and for nearly noiseless observations.

## 3.2.

### Superresolution for Interferometric Phase

The reconstruction results for the Lena phase test-image with the superresolution factor ${r}_{s}=8$ are shown in Figs. 3 and 4. The cross sections for phase and amplitude are shown for middle horizontal lines, where the red (solid) and blue (dotted) curves correspond to the reconstructions and true images, respectively. The reconstructions in Fig. 3 are of the high accuracy. They are obtained for the low level noise ($\chi =1000$, $\mathrm{SNR}=60\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{dB}$). The results in Fig. 4 are much worse but we need to take into account that the observations are very noisy ($\chi =1$, $\mathrm{SNR}=30.5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{dB}$) for superresolution with the factor ${r}_{s}=8$. Thus, we can treat these results as of acceptable quality. It can be noted that the low level noise reconstruction is, visually, nearly perfect. The accuracy of reconstruction is good for both phase and amplitude. The experiments produced for higher order resolution (not shown) demonstrate a noticeable degradation of results for ${r}_{s}=16$ and fail completely for ${r}_{s}=32$.

## 3.3.

### Absolute Phase Imaging

SR-SPAR phase imaging for shear plane phase distribution $(256\times 256)$ with the maximum value of about 65 rad is shown in Figs. 5Fig. 6–7. In Fig. 5, the results are shown for the resolution factor ${r}_{s}=8$. Visually, the obtained 3-D surface is very close to the true one, thus, it is not necessary to show it. The 2-D images in these figures show the wrapped phase and amplitude reconstructions. As it is seen from the wrapped phase, the errors in the reconstruction can be noticed. Nevertheless, the quality of this superresolution reconstruction is very good. In Fig. 6, the similar results are shown for the much higher resolution factor ${r}_{s}=16$. The shown 3-D surface is of a quite acceptable quality while the wrapped phase and amplitude reconstruction, shown as 2-D images, definitely demonstrate that these results are of a much lower quality than those achieved for ${r}_{s}=8$. Finally in Fig. 7, we show the attempt to get reconstruction for the superresolution factor ${r}_{s}=32$. These results are definitely negative, and the phase reconstruction failed.

Phase imaging for Gaussian phase distribution $(256\times 256)$ with maximum value of about 50 rad is shown in Figs. 8 and 9. The reconstructions are obtained for quite noisy data $\chi =100$ and the superresolution factors ${r}_{s}=8$ and ${r}_{s}=16$. For ${r}_{s}=8$ (Fig. 8), the quality of the reconstruction is very good and the 3-D phase reconstruction is very close to the true Gaussian phase object. The situation becomes much worse for the superresolution factor ${r}_{s}=16$ (Fig. 9). The errors in the phase reconstruction are obvious and quite large. The attempt to get reconstruction for the superresolution factor ${r}_{s}=32$ failed and we do not show these images.

In conclusion of this section, we wish to note that we are talking about very high levels of the pixel superresolution ${r}_{s}=8$ and ${r}_{s}=16$ corresponding to the wavelength and half wavelength superresolution.

## 3.4.

### More on Subwavelength Resolution

Let us demonstrate a few interesting tests on subwavelength imaging with ${r}_{s}=32$ corresponding to ${r}_{\lambda}=0.257$ with the object size $128\times 128$.

In Figs. 10 and 11, the reconstructions for the two-peak phase object are shown. The four squares clearly seen in the amplitude reconstructions are patterns of four pixels of SLM each covering $32\times 32$ computational pixels of the object. These images confirm that both developed algorithms SR-GS and SR-SPAR are able to reconstruct two pointwise phase peaks separated by the distance equal to $0.257\lambda $. It is a demonstration of the subwavelength resolution.

## 3.5.

### Parameters of the SR-SPAR Algorithm

The performance of the SR-SPAR algorithm essentially depends on its parameters. Optimization can be produced for each magnitude/phase distribution and noise level. However, in our experiments, the parameters are fixed. The image patches in BM3D are of size $8\times 8$. The group size is limited to 39 patches. The step size among the neighboring patches is equal to 3. The transforms DCT (for patches) and Haar (for the group length) are used for 3-D group data processing in BM3D. In the shown results as an initial guess for the iterative SR-GS and SR-SPAR algorithm, we use an image with the invariant amplitude equal to 1.3 and zero phase.

The parameters defining the iterations of the algorithm are as follows: ${\gamma}_{1}=1/\chi $, $t{h}_{a}=4.0$, and $t{h}_{\phi}=4.0$. The number of the iterations is fixed to 50.

For our experiments, we use MATLAB R2015a and a computer with the processor Intel(R) Core(TM) i7-4800MQ@ 2.7 GHz.

The complexity of the algorithm is characterized by the time required for processing. For 50 iterations, $L=12$ and $256\times 256$ images this time is as follows: $\mathrm{SR}\text{-}\mathrm{GS}\simeq 2500$ s; SR-SPAR without phase unwrapping $\simeq 3300\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$; SR-SPAR with phase unwrapping $\simeq 1500\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$.

## 4.

## Conclusion

Computational superresolution phase retrieval is considered for phase-coded intensity observations. The proposed algorithm is derived as an optimal solution for Poissonian noisy observations. One of the essential instruments of the algorithm is a sparsity hypothesis applied to both phase and amplitude. The efficiency of the algorithm is confirmed by simulation experiments. It is shown that high level superresolution can be achieved with the pixel superresolution factor up to 32, i.e., the pixel size of the reconstructed object in 32 times smaller than the pixel size of the sensor and the SLM. In comparison with the wavelength, the superresolution up to one-quarter of the wavelength is demonstrated.

## References

## Biography

**Vladimir Katkovnik** received his PhD and DSc degrees in technical cybernetics from Leningrad Polytechnic Institute (LPI) in 1964 and 1974, respectively. From 1964 to 1991, he was an associate professor and then a professor in the Department of Mechanics and Control Processes, LPI. Since 2003, he has been with the Department of Signal Processing, Tampere University of Technology, Finland. He published six books and over 350 refereed journal and conference papers. His research interests include stochastic image/signal processing, nonparametric estimation, computational imaging, and computational phase imaging.

**Karen Egiazarian** received his MSc degrees from Yerevan State University, in 1981, his PhD from Moscow State University, Russia, in 1986, and the DTech degree from Tampere University of Technology (TUT), Finland, in 1994. He is a professor leading the Computational Imaging Group, Signal Processing Laboratory, TUT. He has authored about 650 refereed journal and conference papers. His research interests include computational imaging, sparse coding, and image and video restoration. He serves as an associate editor for the *IEEE Transactions of Image Processing* and is the editor-in-chief of the *Journal of Electronic Imaging*.