1. Introduction

In lensless imaging, only the magnitudes of Fourier transformation are recorded, and the central computational task, known as “Phase Retrieval”, is to recover underling images from phaseless measurements. It is a challenging inverse quadratic problem. For noiseless case, related multivariable quadratic systems should be solved, the complexity of which is possibly NP-complete [1]. In practice, the measurements are usually contaminated by noise, and one shall seek an approximate solution by solving the reformulated noncovex optimization problem based on the maximum a posteriori (MAP) estimate of noise.

Numerically, researchers have developed numerous iterative algorithms, such as an alternating projection algorithm [2, 3] for the classical phase retrieval problem (Classical phase retrieval refers to retrieval phase from Fourier measurements). Several other heuristic first order algorithms popular in the optics community including [4–7] were proposed to solve the classical phase retrieval problems, and one can refer to [8,9] and the references therein. These methods are based on projections onto a nonconvex set, and it is therefore very difficult to describe and mathematically prove the convergence. Consequently, more recent work has focused on designing fast algorithms with convergence guarantee for such nonconvex minimization problem, e.g. [10–16], with random-masking or oversampling measurements. In the optics community, in order to increase the imaging resolutions and robustness, a popular approach, known as ptychography [17], uses a sequence of phaseless diffraction measurements recorded while translating the specimen with respect to a constant illumination mask. In order to prevent getting trapped into a local minimum for solving conventional nonconvex minimization problem by projection methods, another scheme is to model phase retrieval by convex methods including PhaseLift [18] and PhaseCut [19] based on semi-definite programming(SDP), as well as PhaseMax [20, 21] operating in the original signal space with much less computational cost compared with SDP based convex methods.

An early work for phase retrieval by employing sparse prior of images was the shrink-wrap algorithm [22], which was established by iteratively shrinking the size of support set for unknown objects. It has been successfully applied to several ground breaking X-ray Free Electron Lasers based experiments [23], which also demonstrates that sparse prior should be considered for modeling phase retrieval problem in order to increase the robustness and accuracy. The L¹ norm based phase retrieval method with a compressibility constraint was proposed in [24,25]. SDP-based convex methods combined with L¹ regularization of the lifted matrix were proposed to solve a sparse phase retrieval problem [26, 27]. With the prior knowledge of exact sparsity level, a sparse Fienup method [3] was proposed in [28]. A generic two-step iterative scheme was proposed in [29], which consisted of updating support set, and using damped Gauss-Newton algorithm to solve a nonconvex optimization problem with updated support set. A two-stage solution technique was presented in [30], where first by using arbitrary phase retrieval technique, sampling measurements were recovered and then by compressive sensing the unknown signals were recovered as well. If the transform domain of the unknown objects are corrupted by the additive noise, a probabilistic method based on the generalized approximate message passing algorithm was given in [31]. Hard thresholding based truncated gradient iterations [32] were employed in order to refine a sparse orthogonality-promoting initialization. An L⁰ regularized variational model was established for the sparse phase retrieval problem accompanied by an efficient algorithms based on the alternating direction method of multipliers [33] with adaptive stepsizes. Under the assumption of objects possessing a sparse representation in the transform domain, shearlet and total variation based regularization methods were considered in [34,35]. Dictionary learning methods were proposed to reconstruct the image [36,37] from Gaussian noised measurements with outliers, where Gaussian distributions only provide a limited approximation in most real applications.

Roughly speaking, the intensities are recorded by detectors within a limited time and incident flux, which causes Poisson noise. For linear measurements, in order to remove the Poisson noise, one can adopt the Anscombe root transformation [38,39], where a pointwise nonlinear transformation (variance stabilizing transformation) was introduced such that one just need deal with Gaussian data. The total variation based method was proposed in [40], where Kullback-Leibler (KL) divergence was derived based on the Bayesian framework in the case of Poisson noise. In a similar manner, for phaseless nonlinear measurements of phase retrieval, variational methods based on Tikhonov and total variation regularization were proposed to recover images from Poisson measurements [41, 42]. The recovered images [41] exhibited sharp edges and clean backgrounds reconstructed from very noisy and limited data. However, as such method was built using the local-sparsity regularization term, there exist visible staircase artifacts, and the smaller repetitive features (textures) could not be preserved, particularly for measurements generated with deterministic illuminations in the ptychographic phase retrieval problem. In order to further increase the accuracy and visual quality of recovery results for practical applications, one possible direction is to exploit the redundancy between overlapping image patches to boost the sparsity, such as using dictionary learning based denoising method [43–45], which pursues a sparse representation by a dictionary learned from image patches.

In this paper, we propose a novel dictionary learning model to reconstruct high-quality images from Poisson phaseless measurements, which can deal with both real and complex-valued images. It consists of three terms, i.e. a representation term by an orthogonal dictionary, an L⁰ pseudo norm of the coefficient matrix, and a Kullback-Leibler divergence term to fit phaseless Poisson data. The pseudo L⁰ norm is used as a sparse regularization of the coefficient matrix with both isotropic and anisotropic forms. Meanwhile, in order to solve the proposed nonconvex nondifferential optimization model, fast iterative algorithms are designed with the convergence guarantee. Since an orthogonal dictionary is employed, subproblems w.r.t. the dictionary and coefficients both have closed-form solutions. Numerous experiments for coded diffraction and ptychographic patterns are performed, which demonstrate that our proposed methods can recover images with sharp edges, clean background compared with the phase retrieval methods without any regularization. Moreover, these methods are capable of producing higher quality results by preserving the texture features compared with the total variation based method [41]. Additionally, by unitizing the introduced anisotropic L⁰ pseudo norm, the proposed methods further increase the quality for reconstructing complex-valued images.

The rest of this paper is organized in the following way. The proposed dictionary learning model is given in section 2. An alternative minimization method is given in section 3. Then a proximal alternating linearized minimization method is further designed in section 4 with the theoretical convergence guarantee. Numerous experiments are performed in section 5 to demonstrate the efficiency of proposed methods. We conclude the paper in section 6.

2. Proposed model

We consider phase retrieval for a 2-dimensional (2D) image in a discrete setting, i.e., an underlying object u : Ω = {0, 1, · · · , n − 1} → ℂ is of size n with n = n₁ × n₂, in which we represent a 2D object in terms of a vector of size n by a lexicographical order. For the classical phase retrieval problem, the measured data are only the magnitudes of the discrete Fourier transform of u, i.e., |ℱu|², where | · |² denotes the pointwise square of the absolute value of a vector, ℱ : ℂⁿ → ℂⁿ denotes the discrete 2D Fourier transform (DFT). In this paper, we consider a general phase retrieval problem [18], i.e.

We will first give a brief review of dictionary learning based denoising for the linear inverse problem, and regularized methods for denosing noisy phaseless measurements contaminated by Poisson noise in the following subsection.

2.1. Review of dictionary learning and Poisson denoising

2.1.1. Dictionary learning and denoising

Given a matrix Y = [Y(0), Y(1), · · · ,Y(d − 1)] ∈ ℂ^l,d, whose columns represent 1-D signals, the target of dictionary learning (or sparse coding) [46] is to find a “sparse” linear representation for the matrix Y by a learned dictionary D = [D(0), D(1), · · · , D(c − 1)] ∈ ℂ^l,c, i.e. Y ≈ Dα, such that the coefficient matrix α has more zeros (“sparse”) entries. Obviously, the dictionary D and coefficient matrix α ∈ ℂ^c,d are coupled and both needed in order to derive an optimal sparse representation. Roughly speaking, it can be reformulated as below:

The well-known work named K-SVD denoising [43] is to first produce redundant signals by sliding a small window over an entire image, and then to solve a sparse coding problem [47]. Given the noisy image f, the K-SVD denoising model can be given as follows:

2.1.2. Regularized Poisson denoising for phase retrieval

Thibault and Guizar-Sicairos [42] proposed the following model in coherent diffractive imaging with the Tikhonov regularization min_u∈ℂⁿ λ‖∇u‖² + ℬ(|𝒜u|², f), with the Kullback-Leibler divergence denoted as follows,

2.2. Proposed dictionary learning phase retrieval (DicPR)

Assuming that the measurements f = Poisson(|𝒜u|²) ∈ ℝ^m in (1) are corrupted by Poisson noise, we establish a minimization problem driven by the patch based sparsity [43–45] of the underlying images, referred to as “DicPR”, which reads,

Instead of directly using K-SVD [43], we have proposed an orthogonal dictionary learning model such that subproblems involving D and α both have closed forms. Moreover, the denoising results with the orthogonal dictionary seem comparable with the nonorthogonal one [48]. Due the the existence of the data fitting term for phaseless measurements and L₀ regularization of the coefficient matrix, the proposed model is a nonconvex and nonsmooth optimization problem, which is difficult to design the corresponding algorithm with the convergence guarantee. In the following sections, we first give a simple Alternating Minimization Method (AMM). Furthermore, a convergent algorithm called Proximal Alternating Linearized Minimization method (PALM) will be provided.

3. Algorithms 1: Alternating minimization method (AMM)

Since the variables u, D, α in (5) are coupled together, the overall iterative algorithm of AMM to solve the problem consists of three steps:

3.1. Subproblem w.r.t. the variable u

In order to avoid directly calculating the gradient of the KL divergence (which is possibly unbounded), the alternating direction of multiplier method (ADMM) is adopted to solve the above subproblem. By introducing an auxiliary variable z, we have:

For the second subproblem of (12), we have

We remark that alternatively one can consider a gradient descent type algorithm with adaptive stepsizes such as Wirtinger flow [10], or the nonlinear conjugate gradient method or Newton type algorithm [49] if the number of measurements are sufficient. The ADMM can work pretty well with very few measurements [41], and has low computation cost compared to Newton type algorithm. In the experimental tests, very few iterations are adopted.

3.2. Subproblems w.r.t. variables D and α

For the second subproblem of (8) w.r.t. D, we have

For the third subproblem of (8) w.r.t α, since D^*D = I, we have

Algorithm 1:. AMM for “DicPR” in (5)

View Table | View all tables in this article

In order to accelerate the computation, in the first iterations (three iterations in the numerical experiments), we use a truncated version for (20) where only parts of all columns of U and V are selected corresponding to the few largest singular values.

Assumption 1 The sequence ${\left\{u^k\right\}}_{k=0}^{\infty }$ is bounded if and only if ${\left\{\left|𝒜u\right|\right\}}_{k=0}^{\infty }$ is bounded.

In order to guarantee the uniqueness of the solution for phase retrieval, oversampling or multiple measurements are required. For example, for the coded diffraction patter and ptychographic pattern, 𝒜^*𝒜 is a diagonal matrix as (15). One readily sees that Assumption 1 holds, if we assume that the diagonal elements are all non-zeros. 𝒜^*𝒜 is usually invertible.

Theorem 1 Under Assumption 1, there exists an accumulation point Z̃ of {Z^k} such that lim_k→∞ ϒ(Z^k) = ϒ(Z̃).

Proof: Readily one knows that the objective functional values of the iterative sequences generated by Algorithm 1 are non-increasing, i.e. ϒ(u^k+1, D^k+1, α^k+1) ≤ ϒ(u^k, D^k, α^k). Hence, there exists a positive constant C independent with k, such that ℋ(u^k, D^k, α^k)+ ℱ(u^k)+ℐ_K (D^k)+ τ‖α^k ‖₀ ≤ C. First we have ℱ(u^k) ≤ C and readily one has {|𝒜u^k |} is bounded. By Assumption 1, the boundedness of {u^k} is proved. Since (D^k)^*D^k = I, it implies the boundedness of {D^k}. Since ℋ(u^k, D^k, α^k) ≤ C, we have ‖α^k −(D^k)^*R(u^k)‖ ≤ C by the orthogonality of the dictionary, which leads to the boundedness of {α^k}. There must exist a subsequence ${\left\{Z^{k_n}\right\}}_{n=0}^{\infty }\subset \left\{Z^k\right\}$ with the limit point Z̃ = (ũ, D̃, α̃), s.t. lim_n→∞ Z^k_n = Z̃. Therefore we conclude the theorem.

For AMM algorithm, we only prove that it converges to the accumulation point of the iterative sequences. Generally speaking, for a nonconvex problem, the AMM algorithm may cycle indefinitely without converging [52]. Although numerically AMM for the proposed model works well, it is still challenging to study the convergence mathematically. Inspired by the proximal linearized technique, a proximal type splitting algorithm will be further given, where the sufficient decrease of the objective function is guaranteed due to the introduction of the additional proximal terms such that the convergence of such algorithm can be derived following [50].

4. Algorithm 2: Proximal alternating linearized minimization (PALM)

We adopt the proximal alternating linearized minimization method (PALM) [50] to solve (5). First we give the partial derivatives of

For u–subproblem in (26), one can readily derive the ADMM following Algorithm 1-1, where only the solver for u in Step 2 is slightly different, and we just list Algorithm 2-1 directly.

Algorithm 1-1:. ADMM for u–subproblem in (8)

View Table | View all tables in this article

Algorithm 2-1:. ADMM for u–subproblem in (26)

View Table | View all tables in this article

For the D–subproblem in (27), we need to solve the following problem:

For the α–subproblem in (28), one can directly get the closed form represented by the hard thresholding as

Therefore we can give an overall PALM for (5) as below.

In the following parts, we will provide convergence analysis for Algorithm 2. For simplicity, the analysis is conducted for real-valued images, and it can be easily generalized to the complex-valued images, where one only needs to reformulate the related norms and operators of variables w.r.t. their corresponding real and imaginary parts respectively [41].

Algorithm 2:. PALM for “DicPR” in (5)

View Table | View all tables in this article

The existence of the critical point of proposed model is given below.

Lemma 4.1 The critical point of (10) exists with the finite functional value.

Following Lemma 5 in [50] and the lower semi-continuity of ℱ, ℐ_D and ‖ · ‖₀, one can readily prove it, and details are omitted here. We need the assumption of the boundedness of u^k.

Assumption 2 (Boundedness) The sequences {u^k} generated by (26) is bounded.

Based on above assumption, the boundedness of the iterative sequences (u^k, D^k, α^k) can be proved, and see details in Lemma .1 in the appendix.

Lemma 4.2 (Lipschitz-Gradient) Under Assumption 2, for the functional ℋ(u, D, α), we have

The proof for the first item is trivial. By (25) and Assumption 2, one can readily finish the proof for the left parts, based on the boundedness of iterative sequences proved by Lemma .1 in the appendix. Therefore, we omit the details. Furthermore, We assume that stepsizes are big enough.

Assumption 3 (Step Sizes) ${\lambda }^{+}:=\text{min}\left\{c^k-{\lambda }_u^{+},d^k-{\lambda }_D^{+},e^k-{\lambda }_{\alpha }^{+}\right\}>0$.

Then we can consider the convergence of Algorithm 2. Readily one knows that ϒ(u, D, α) is semi-algebra [50,51] with the data term ℬ(·, ·) in (4), and finally we can get the final convergence theorem as follows.

Theorem 2 Let Assumptions 2, 3 hold and min_k e^k > 1/2. The sequence {u^k, D^k, α^k} generated by Algorithm 2 converges to a critical point of (6).

Based on Lemma .2 and Lemma .3 in the appendix, one can readily finish the proof following Theorem 1 in [50], and therefore, we omit the details here.

Remark 4.1 Although in this paper, we focus on the orthogonal dictionary, the proposed PALM can also be applied to the case with non-orthogonal dictionary, and further speedup of PALM should be investigated as in [51]. We leave them as future work.

At the end of this section, we analyze the computational complexity of Algorithm 1 and Algorithm 2 which have similar cost. Assuming that the operator 𝒜 generates Fourier masked measurements which means fast Fourier transformation can be adopted, the computation complexity for Algorithm 1-1 and Algorithm 2-1 is of O(n log(n)T_in + ln) after T_in inner iterations. The complexity for Step 3 and 4 is O(l³ + l²d). Hence the complexity of Algorithm 1 and Algorithm 2 is O(n log(n)T_in + ln + l³ + l²d)T_out) after T_out outer iterations. In our experiments, we set d ≈ n, l ≪ n and as a result the complexity is about O(n(log(n)T_in + l²)T_out).

5. Numerical experiments

All the tests are performed on a laptop with Intel I7-5600U2.6GHZ, and 16GB RAM, and the codes are implemented in MATLAB. The dictionary is initialized by discrete cosine transform. The model parameters η and τ for Algorithm 1 and Algorithm 2 are selected by hand. In the inner iterative algorithm e.g. Algorithm 1-1 and Algorithm 2-1, we set default parameter r = 1 × 10⁻³, and the default inner iteration number to five. For the left parameters for Algorithm 2, one can use dynamic schemes to update the parameters as [50,51]. In our experiments, for simplicity, these parameters are fixed. We stop Algorithm 1 and Algorithm 2 after a given maximum outer iteration number T to guarantee the convergence, which will be specified in the following subsections. All the other needed parameters will also be addressed in the following subsections if we do not give or use the default values.

Set image patch size to 8 × 8 empirically such that l = c = 64, and the number of patches $d={\left(\sqrt{n}-7\right)}^2$ for square images with n pixels. In this paper we have given a framework for phase retrieval with arbitrary linear operator 𝒜. However, it it more practical to consider the Fourier type transform, and we will show the performance on Fourier masked measurements involved with two types of patterns: Coded diffraction pattern (CDP) with random masks and ptychographic phase retrieval with deterministic masks generated by zone plate lens. Especially, ptychographic phase retrieval is a very promising technique to generate high resolution images with large field of view compared with the traditional diffraction imaging and meanwhile it requires less temporal and spatial coherence. However, the related phase retrieval problem is more challengeable.

The ground truth images are provided in Fig. 1, where four real-valued images with resolution 512 × 512 are put in Figs. 1(a)–1(d), and a complex-valued image with resolution 256 × 256 is put in Figs. 1(e)–1(g). Given a ground truth u, the noisy measurement is generated as f(j) = Poisson(|(𝒜u_δ)(j)|²) ∀j ∈ Ω̃, with u δ = δu at peak level δ (Noise level increases as peak level δ decreases). We measure the quality of the reconstructed image ũ by signal-to-noise ratio (SNR) (The SNR of noise free image is + ∞. The reconstructed image with larger SNR usually means higher visual quality) SNR (ũ, u) = −20 min_{ς∈{ς ∈ℂ: |ς|=1}} log(‖ςũ − u ‖/‖ũ‖), with the ground truth image u, where ς denotes the phase factor for the trivial ambiguities. In order to measure the sparsity of coefficient matrix α ∈ ℂ^c,d or ℝ^c,d, we introduce the sparsity level (Smaller sparsity values means sparser of the data with a higher level) $S(\mathit{\alpha })=\frac{{\Vert \mathit{\alpha }\Vert }_0^{\mathit{iso}}}{c\times d}\times 100\%$.

 

Fig. 1 Ground truth images. Real-valued images of resolution 512 × 512 in (a): Peppers; (b): Fingerprint; (c): Barbara; (d): House; A complex-valued image “Goldballs” of resolution 256 × 256 with magnitude, real and imaginary parts in (e), (f) and (g), respectively.

Download Full Size | PDF

We compare our proposed methods with ADMM for phase retrieval algorithm without any regularization [9] and the total variation (TV) based Poisson noise removal method [41]. Denote ADMM for phase retrieval [9], TV based Poisson denoising method [41], our proposed Algorithm 1 and Algorithm 2 with isotropic L ⁰ norm by “PR”, “TVPR”, “ALGI” and “ALGII”, respectively. We denote Algorithm 1 with anisotropic L ⁰ term for complex-valued images by “ALGI_aniso”. Since our proposed algorithms solve the nonconvex optimization problem, they are initialized by the results of “PR” in order to increase the robustness and convergence speed.

5.1. CDP for real-valued and complex-valued images

For the first pattern, the octanary CDP is explored, and specifically each element of I_j in (14) takes a value randomly among the eight candidates, i.e., {$\pm \sqrt{2}/2$, $\pm \sqrt{2}\mathbf{i}/2$, $\pm \sqrt{3}$, $\pm \sqrt{3}\mathbf{i}$}. Set the number of masks K = 2, 4 as [41] for real-valued and complex-valued images, respectively.

5.1.1. Real-valued Images

We will show the performance of our proposed algorithms from noisy measurement on four real-valued images shown in Figs. 1(a)–1(d), with peak level δ ∈ {5.0 × 10⁻³, 1.0 × 10⁻²}. Set η = 8.0 × 10⁻⁶, τ = 4.5 × 10⁻⁴ in Algorithm 1 and Algorithm 2. Set stepsizes c^k = 10, d^k = 50, and e^k = 1.5, 1.2 for δ = 5.0 × 10⁻³, 1.0 × 10⁻² respectively in Algorithm 2. Both algorithms stop after T = 100 iterations. All four images share the same parameters with the others.

Reconstructed results with different noise levels are put in Fig. 2–Fig. 3, and zoom-in parts of corresponding results are put in Fig. 4–Fig. 5. Readily one can see that “PR” generates very noisy results by observing the first rows of Fig. 2–Fig. 3 and the second rows in the zoom-in results of Fig. 4–Fig. 5, where the edges, background and the repetitive structures are contaminated severely. By TV regularization method and our proposed dictionary learning methods, all the recovery images have sharper edges and cleaner background. However, by observing the zoom-in results in the third rows of Fig. 4–Fig. 5, “TVPR” generates images with visible stair case artifacts, and some important texture information cannot be kept, while “ALGI” and “ALGII” can both produce high quality images and one does not find any visible staircase artifacts in the zoom-in parts. By observing Figs. 4(m), 4(q) or Figs. 5(m), 5(q), “ALGI” and “ALGII” generate cleaner background than “TVPR”. Moreover, by Figs. 4(n)–4(p), 4(r)–4(t) or Figs. 5(n)–5(p), 5(r)–5(t), the textures are preserved pretty well compared with the results by “TVPR”. In order to qualify the improvement of our proposed methods, the corresponding SNRs are put in Table 1 and Table 2, where the maximum in the same row is marked in bold font. Inferred from these two tables, SNRs by “TVPR”, “ALGI” and “ALGII” are about double of those by “PR”, and SNRs by “ALGI” and “ALGII” increase averagely about 1.5dB and 2dB, respectively. The increase or improvement with the proposed methods compared with “PR” and “TVPR” is more obvious on the image with more texture information e.g. “Fingerprint” or “Barbara” than that e.g. “Peppers” with piecewise smooth features, which demonstrates that its advantage is to handle images with sophisticated texture features.

 

Fig. 2 CDP with δ = 5.0 × 10⁻³. First row: PR; Second row: TVPR; Third row: ALGI; Fourth row: ALGII.

Download Full Size | PDF

 

Fig. 3 CDP with δ = 1.0 × 10⁻². First row: PR; Second row: TVPR; Third row: ALGI; Fourth row: ALGII

Download Full Size | PDF

 

Fig. 4 Zoom-in images of Fig. 2 for CDP with δ = 5.0 × 10⁻³. First row: Ground truth; Second row: PR; Third row: TVPR; Fourth row: ALGI; Fifth row: ALGII

Download Full Size | PDF

 

Fig. 5 Zoom-in images of Fig. 3 for CDP with δ = 1.0 × 10⁻². First row: Ground truth; Second row: PR; Third row: TVPR; Fourth row: ALGI; Fifth row: ALGII.

Download Full Size | PDF

Table 1. SNRs for CDP with δ = 5.0 × 10⁻³.

View Table | View all tables in this article

Table 2. SNRs for CDP with δ = 1.0

View Table | View all tables in this article

By observing the zoom-in parts, the results in the fourth rows of Fig. 4 and Fig. 5 seem more smoothing by “ALGI” than those in the fifth rows by “ALGII”. On the other hand, the results by “ALGII” contain a bit more features than by those by “ALGI”. In order to investigate the performance differences of “ALGI” and ”ALGII”, the sparsity of the coefficient matrix α is put in Table 3, and one can readily see that the coefficient matrix is sparser by “ALGI” than by “ALGII”, which possibly leads to those above behaviors. Anyway, both two algorithms produce comparable results with cleaner background and well-preserved textures. Hereafter, we only provide the results by “ALGI” since it has fewer parameters and seems more suitable for practical use.

Table 3. Sparsity S(α) for CDP on real-valued images.

View Table | View all tables in this article

We also report the runtimes of all compared algorithms in Table 4, where one can readily see that the proposed dictionary learning algorithms are much slower than the other two algorithms, since the inner loops for the subproblem of u are needed in Algorithm 1-1 and Algorithm 2-1. In the future, we will explore much faster algorithm without inner loops.

Table 4. Runtimes (in seconds) for CDP with δ = 5.0 × 10⁻³.

View Table | View all tables in this article

5.1.2. Complex-valued Image

More experiments on the complex-valued image “Goldballs” in Fig. 1 for CDP with peak level δ ∈ {0.08, 0.1} are performed, and please see the reconstructed results in Fig. 6. Notice that for complex valued images, we introduce two kinds of definitions for L⁰ pseudo norm in (7), and therefore we have two versions of PALM, i.e. “ALGI” and “ALGI_aniso”. Set the parameters η = 3.0 × 10⁻⁵, τ = 1.0 × 10⁻³ for “ALGI” and η = 2.5 × 10⁻⁵, τ = 8.0 × 10⁻⁴ for “ALGI_aniso”. Set iteration number in the outer loop as T = 50. By observing Fig. 6, from noisy measurements one can only derive noisy results by “PR”. By “TVPR”, “ALGI” and “ALGI_aniso”, the recovery results are almost noise free and also have very clean background. “TVPR” can only recover the large scale structures at the left top corner, but cannot keep the smaller repetitive structures. Obviously our proposed methods can produce better results, where both the large and small structures can be preserved well. Related SNRs are put in Table 5, about 1dB, 2dB increases are gained by “ALGI”, and “ALGI_aniso” compared with “TVPR”.

 

Fig. 6 CDP for complex-valued image with δ ∈ {0.08, 0.1}. First row: δ = 0.08; Second row: δ = 0.1. From left to right: reconstructed images by PR in (a) and (e), TVPR in (b) and (f), ALGI in (c) and (g) and ALGI_aniso in (d) and (h).

Download Full Size | PDF

Table 5. SNRs of CDP for complex-valued image

View Table | View all tables in this article

One also notices that “ALG_aniso” produces the recovery results with higher accuracy than “ALGI”. In order to investigate the improvements by the anisotropic norm, we show the learned dictionaries in Fig. 7, and the sparsity of coefficient matrix α in Table 6. Especially in Figs. 7(d) and 7(h) the imaginary parts of learned dictionary by “ALGI_aniso” seem to have more features than those in Figs. 7(b) and 7(f) by “AlGI”. Meanwhile the real parts of learned dictionaries by the anisotropic version algorithm are quite close to the imaginary parts, which is consistent with the similarity of structure of real and complex parts of ground truth in Figs. 1(f) and 1(g). Therefore, it produces better dictionaries by the anisotropic style in (7). Moreover, due to such better dictionaries learned by “ALGI_aniso”, the representation is much sparser, i.e. the overall sparsity level S(ℜ(α)) + S (ℑ(α)) by “ALGI_aniso” is much higher than by “ALGI”, since S (ℑ(α)) is much smaller by “ALGI_aniso” than by “ALGI”, and corresponding sparsity levels of the real part are almost same in Table 6. It demonstrates that the anisotropic version algorithm can help to improve the image qualities for complex-valued images by learning the better dictionary and sparser coefficient matrix.

 

Fig. 7 Learned dictionaries for CDP corresponding to the results in Fig. 6. First row: δ = 0.08; Second row: δ = 0.1. From left to right: Real and complex parts of learned dictionaries by ALGI in the first two columns, and by ALGI_aniso in the third and fourth columns.

Download Full Size | PDF

Table 6. Sparsity for CDP on complex-valued images with different L⁰ pseudo norms.

View Table | View all tables in this article

5.2. Ptychographic phase retrieval (PtychoPR) for complex-valued image

5.2.1. Fixed sliding distances

For PtychoPR, a zone plate lens and responding illumination mask are employed as in [13]. Set sliding distance SlidDist = 16. The number of frames is 16 × 16 and the frame size is of 64 × 64, and in such setting total number of measurements m = 16n. Set peak level δ ∈ {0.2, 0.5}. Here only two inner iterations for PALM are sufficient. Set outer iteration number T = 50, τ = 5.0 × 10⁻³, η = 3.0 × 10⁻² for δ = 0.2 and τ = 8.0 × 10⁻³, η = 3.0 × 10⁻² for δ = 0.5. The reconstructed results are shown in Fig. 8, and the corresponding SNRs are shown in Table 7. Inferred from Figs. 8(a) and 8(e), it seems very blurry in the results of “PR” from noisy measurements, since the phaseless data are generated by structured deterministic illumination, and corrupted low frequency parts are worse than high frequency parts. By “TVPR”, at higher noise level δ = 0.2 in the first row of Fig. 8, it can produce results with sharp edges for large scale features and clean background, but cannot preserve the smaller features at all. For the case with peak level δ = 0.5, “TVPR” can produce pretty good recovery results for both smaller and larger scale features. The proposed “ALGI” and “ALGI_aniso” can recover the smaller features very well especially at peak level δ = 0.2. Inferred from Table 7, SNRs are increased about 3.5dB, 4dB for “ALGI” and “ALGI_aniso”, respectively compared with “TVPR” at peak level δ = 0.2; SNRs are increased about 1.5dB, 1.7dB for “ALGI” and “ALGI_aniso” respectively compared with “TVPR” at peak level δ = 0.5. Such gains in SNRs show that our proposed algorithms can produce more accurate results. Similarly to the previous subsection for complex-valued image, the anisotropic version algorithm “ALGI_aniso” has better performances than isotropic version algorithm.

 

Fig. 8 PtychoPR for complex-valued image with δ ∈ {0.2, 0.5}. First row: δ = 0.2; Second row: δ = 0.5. From left to right: reconstructed images by PR in (a) and (e), TVPR in (b) and (f), ALGI in (c) and (g) and ALGI_aniso in (d) and (h).

Download Full Size | PDF

Table 7. SNRs of PtychoPR for complex-valued image corresponding to Fig. 8

View Table | View all tables in this article

5.2.2. Variable sliding distance

In order to further study the robustness of the proposed algorithms, more experiments with different sliding distances (collecting less data by increasing the sliding distances) are done. Set peak level δ = 0.2. Same parameters as the previous tests are used. Performances for PtychoPR with different number of measurements by increasing SlidDist are shown in Fig. 9, where m/n = 12.25, 9, and 7.56 when the sliding distances SlidDist = 18, 20, 22 respectively. On can see that in the first column of Fig. 9, the results of “PR” are not only blurry, but also contain some visible structured artifacts. “TVPR” can remove such artifacts, and recover some edges of large scale feature. Our proposed algorithms can further recover most of smaller features in Figs. 9(c) and 9(d). In an extreme case with SlidDist = 22 shown in the third row of Fig. 9, the proposed algorithms can also recover some parts of smaller features. We also put the corresponding SNRs in Table 8, where one can see the obvious gains at about 1.3dB, 1.6dB for “ALGI” and “ALGI_aniso” respectively compared with “TVPR”.

 

Fig. 9 PtychoPR for complex-valued image with different sliding distances (peak level δ = 0.2). From first row to the last row: SlidDist= 18, 20, 22 respectively. From left to right: reconstructed images by PR in (a), (e) and (i), TVPR in (b),(f) and (j), ALGI in (c), (g) and (k), and ALGI_aniso in (d), (h) and (l).

Download Full Size | PDF

Table 8. SNRs of PtychoPR with different sliding distances for complex-valued image corresponding to Fig. 9

View Table | View all tables in this article

5.3. Convergence study

To check the convergence of proposed algorithms, we monitor the histories of SNRs and the successive errors of iterative solution u^k and D^k w.r.t. the iteration number k, which are defined as $\frac{\Vert u^k-u^{k-1}\Vert }{\Vert u^k\Vert }$ and $\frac{\Vert D^k-D^{k-1}\Vert }{\Vert D^k\Vert }$. We show the histories of errors and SNRs in Fig. 10 for CDP on “Peppers”, which implies the proposed algorithms are stable and convergent, which is consistent with the provided theories. It seems that the dictionary converges fast by “ALGII” than by “ALGI”. By inferred from Figs. 10 (c) and 10(f), SNR first increase, and then get stable, which demonstrates that the proposed algorithms are quite robust. Moreover, when the noise level increases, more iterations are needed. It is very interesting and important to give an optimal iteration condition, and we leave it as a future study.

 

Fig. 10 Convergence histories for successive errors and SNRs on “Peppers”. First row: δ = 5.0 × 10⁻³; Second row: δ = 1.0 × 10⁻². Histories of successive error of dictionaries, the iterative solution u and SNRs from the left to right rows respectively.

Download Full Size | PDF

5.4. Performances w.r.t parameters

Finally we provide some experiments to test the parameter impact, especially on the parameters τ and η which balance sparsity term and data fitting term. We show the impact of “ALGI” on “Peppers” with δ = 5.0 × 10⁻³. Particularly, we select (η, τ) ∈ {η⁰ × 2^−x, η⁰ × 2^−x+1, · · · , η⁰ × 2^x−1, η⁰×2^x}×{τ⁰×2^−y, τ⁰×2^−y+1, · · · , τ⁰×2^y−1, τ⁰×2^y}, with x = y = 4, η₀ = 8.0×10⁻⁶, τ₀ = 4.5 × 10⁻⁴, and plot the corresponding SNRs of the reconstructed images in Fig. 11. One can see that too large or smaller parameters decrease the quality of reconstructed results. In the future, an automatic schemes for optimal parameter selection should be investigated for the best results.

 

Fig. 11 SNRs for ALGI w.r.t different parameters τ and η.

Download Full Size | PDF

6. Conclusion

In this paper, we have proposed a novel orthogonal dictionary learning based phase retrieval model to denoise the phaseless measurements contaminated by Possion noise, and further designed two efficient algorithms. The conducted experiments demonstrate that both algorithms are capable of producing higher quality results, especially preserving the texture features compared with the TV denoising method. An incoherent dictionary based method and further speedup of the proposed algorithms should be investigated in the future.