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:
experiments are assumed in Eq. (1), where 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
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 (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 be the Fourier transform (FT) of , . If is given, can be precisely calculated as the inverse FT, . Now let the absolute value of be given and the phase of is unknown. Is it possible to reconstruct the phase of and in this way the original from the amplitude of FT ? In general, the answer is negative and positive only for the special classes of , 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, .
In optics, the priority is different. The phase to be retrieved is phase of the object . The image formation operators in Eq. (1) depend on optical setups. Various methods are developed in order to make these sufficiently different for different and gain observation diversity, enabling finding from observations . Defocusing of the registered images is one of the popular instruments to get a sufficient phase diversity.56.7.–8 In a recent development, a spatial light modulator (SLM) is exploited for defocusing (e.g., see Refs. 9 and 10). The -optical configuration with SLM in the Fourier plane for defocus imitation is proposed in Ref. 11 and further studied in Ref. 12.
The phases in can be generated as deterministic or random. The phase modulation is able to dramatically change the diffraction pattern of by redistribution of the observed intensities from low to high frequencies.
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 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 . 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 . 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.
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 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 , where 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 and intensity equal to 1. The forward propagation operator 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.
For wavefront propagation from the object to sensor planes, we use the paraxial Fresnel modeling. It gives the following link between the object wavefront and the wavefront at the sensor plane [Ref. 38, Eq. (5.19)]:
Here and are complex-valued distributions of the wavefronts at the object plane [lateral coordinates ] and the sensor plane [lateral coordinates ]. is a complex-valued transmission function of SLM. Using FT, the input–output model (6) can be given in the form
Equations (6) to (7) define the forward propagation operator in Eq. (2). Note that the model (6) as used in Eq. (5) is discrete-continuous with continuous and a physical discretization imposed on by the pixelated sensor and by the pixelated SLM on .
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 .
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:
The parameter 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 are equal and are given by , i.e., . Defining the observation signal-to-noise ratio (SNR) as the ratio between the square of the mean and the variance of , we have . It follows that the relative noisiness of observations becomes stronger as () and approaches zero when (). The later case corresponds to the noiseless scenario, i.e., with the probability 1.
The scale parameter is of particular importance for modeling as it allows to control the level of noise in the observations.
Superresolution Sparse Phase Retrieval
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 , 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 is defined by two variables: phase and amplitude .
The sparse representation can be imposed on complex-valued directly using complex-valued basic functions or on the following pairs of real-valued variables:
Remember that an interferometric (wrapped) phase is restricted to the interval , whereas an absolute (unwrapped) phase is different by adding an integer number of to the interferometric phase. In what follows, we denote the interferometric (wrapped) phase of the object as and the corresponding absolute phase as . We introduce the phase-wrap operator , linking the absolute and principal phases as . We also define the unwrapped phase as . Notice that is not an inverse operator for because the latter is highly nonlinear, and for signals of dimension two and higher, there is no one-to-one relation between and .
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 , 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 , the sparsity imposed on phase and amplitude. The variational formulation for reconstruction of a complex-valued , 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
Here and are sparse approximations of and ; 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 and 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,4546.–47 stating that any efficient filter can serve as a potentially good prior and efficient regularizator in variational design of data processing algorithms.
Superresolution SPAR Algorithm
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 , 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 .
Conventionally, the pixels are square of the size and for the SLM and sensor, respectively. In what follows, for simplicity, . A continuous object is discretized by pixels of the size . 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 , the reconstruction of from the observations 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 , we arrive to a much more challenging problem of pixel superresolution or subpixel resolution imaging. Further, if is so small that , then it is wavelength resolution. Going further to we arrive to subwavelength or wavelength superresolution. The superresolution phase retrieval with smaller and very small as compared with and is the goal of this paper.
It is convenient to assume that , where is an integer pixel superresolution factor. In this case, the SLM pixel covers 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 , in the form
Here, the variables and are sampled with the computational period . Then, in particular, the modulation function is a piecewise invariant with 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 , where 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 with the noisy observations obtained according to the Poissonian distribution (8). Note that these computational are given with the computational period while the observations are introduced with the sampling period . 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 for this discrete model initially introduced for the continuous domain variables.
The presented SR-SPAR algorithm is derived from the variational formulation introduced for optimal reconstruction of from Poissonian observations . The corresponding minus log-likelihood for Poissonian observations according to Eq. (8) is as follows:
This criterion should be minimized with respect to , provided Eq. (11) linking and 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: in Ref. 26 and in this paper meaning that the observations should be upsampled by a factor .
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.
SR-SPAR phase retrieval algorithm.
|, , ;|
|2.||Poissonian noise suppression:|
|5.||Sparse phase and amplitude filtering:|
|6.||Object wavefront update:|
The inputs in this algorithm are the observations upsampled by a factor . We use the zero-order upsampling giving 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 .
At step 1, the object wavefront estimate is multiplied by the phase mask and propagated by the operator to the sensor plane, with the result denoted as . These wavefronts are calculated for the diffraction area denoted as .
At step 2, the wavefront is updated to the variable by filtering the amplitude of according to the given observations . The following formula, as derived in Ref. 26, defines the rule on how the updated amplitude is calculated:
These calculations are pixelwise; is the parameter of the algorithm. This update is produced provided known observation , i.e., for the pixels belonging to the sensor area , . In our modeling, the computational diffraction area is always equal or larger than the sensor area, . For the area out of the sensor, the wavefront values are preserved, for .
At step 3, the estimates backpropagate to the object plane and update the object wavefront estimate . Here, means a complex conjugate and .
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 .
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 , 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., is very large.
It was demonstrated in Ref. 42 that using the update shown in step 2 and different for and 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® democodes49 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.
Both algorithms, SR-SPAR and SR-GS, were tested for various models of . In what follows, we are restricted mainly to phase-objects of invariant amplitude and three types of varying phase: test-images Lena normalized to the interval [0, ]; 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 algorithm50 is used for phase unwrapping in step 4 of SR-SPAR.
The fixed parameters of the experiments are: , , sensor size, in pixels, , computational diffraction area of size . The main varying parameters are the computational sampling period , the Poissonian noise parameter , and computational resolution factor calculated with respect to the sensor as . It is assumed that is integer and takes values , 8, 16, and 32. The corresponds to the pixel-resolution and larger 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 .
It is natural to also measure the superresolution with respect to the wavelength as the ratio . Then, gives the pixel-resolution with , i.e., the sensor and SLM pixels are about eight times larger than the wavelength . For , we obtain the wavelength resolution with , the higher values , 32 correspond to the subwavelength resolution with computational pixels smaller than the wavelength with the wavelength resolution factors and 0.257.
Note that according to the restriction (11), smaller (larger ) assumes that the lens with a smaller focal distance should be used in the considered optical setup. For the introduced set , we obtain the following focal distances , respectively.
In our experiments, the phase modulation masks are random with the Gaussian independent zero-mean phase values, .
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 . 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 taking values in the internal [1, 1000]. The smallest results in the noisiest observations. The corresponding SNR is calculated in dB as
For the superresolution experiments, we use the objects with a fixed number of computational pixels of size , thus larger 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.
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 (), provided that only a single experiment is produced, . 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 (), 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, , and for nearly noiseless observations.
Superresolution for Interferometric Phase
The reconstruction results for the Lena phase test-image with the superresolution factor 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 (, ). The results in Fig. 4 are much worse but we need to take into account that the observations are very noisy (, ) for superresolution with the factor . 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 and fail completely for .
Absolute Phase Imaging
SR-SPAR phase imaging for shear plane phase distribution 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 . 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 . 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 . Finally in Fig. 7, we show the attempt to get reconstruction for the superresolution factor . These results are definitely negative, and the phase reconstruction failed.
Phase imaging for Gaussian phase distribution with maximum value of about 50 rad is shown in Figs. 8 and 9. The reconstructions are obtained for quite noisy data and the superresolution factors and . For (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 (Fig. 9). The errors in the phase reconstruction are obvious and quite large. The attempt to get reconstruction for the superresolution factor 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 and corresponding to the wavelength and half wavelength superresolution.
More on Subwavelength Resolution
Let us demonstrate a few interesting tests on subwavelength imaging with corresponding to with the object size .
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 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 . It is a demonstration of the subwavelength resolution.
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 . 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: , , and . 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, and images this time is as follows: s; SR-SPAR without phase unwrapping ; SR-SPAR with phase unwrapping .
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.
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.