PPR: Plug-and-play regularization model for solving nonlinear imaging inverse problems

Signal Processing 162 (2019) 83–96

Contents lists available at ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

PPR: Plug-and-play regularization model for solving nonlinear imaging inverse problems Baoshun Shi a,b, Qiusheng Lian a,b,∗, Xiaoyu Fan a,b a b

School of Information Science and Engineering, Yanshan University, Qinhuang Dao, 066004, Hebei province, China Hebei Key Laboratory of Information Transmission and Signal Processing, Qinhuang Dao, 066004, Hebei province, China

a r t i c l e

i n f o

Article history: Received 9 October 2018 Revised 6 March 2019 Accepted 7 April 2019 Available online 8 April 2019 Keywords: Nonlinear imaging inverse problem Computational imaging Plug-and-play priors Deep denoiser Sparse coding

a b s t r a c t The problem of recovering an image of interest from nonlinear measured data is challenging. To address this nonlinear imaging inverse problem, we propose a novel Plug-and-Play Regularization (PPR) approach that can exploit multiple priors. The underlying image and its filtered image by a denoiser should have similar structures. To exploit this similarity, we enforce the similarity of their sparse coefficients with respect to a tight frame. We formulate a PPR-based nonlinear imaging optimization problem and solve it by using the alternating optimization strategy that consists of filtering step, sparse coding step and image updating step. To avoid the finely-tuned regularization parameter, the epigraph concept is employed in the image updating step. Multiple priors, including the priors employed by the denoiser and the sparsity in the sparsifying transform domain, can be utilized in the proposed PPR model. Under the coded diffraction imaging scenario, we show that the proposed algorithm of exploiting a deep denoiser can achieve higher quality images, compared to the previous imaging algorithms. © 2019 Published by Elsevier B.V.

1. Introduction The imaging technique in a nonlinear sampling system aims to recover the image or object of interest from nonlinear measured data. It arises in various science and engineering fields, including optical imaging [1], holograph [2], X-ray imaging [3], scatter imaging [4], crystallography [5], to name just a few. A typical imaging technique that recovers the image from nonlinear measured data is phase retrieval (PR). In the PR regime, the phase information of the diffraction pattern is lost in the measured data; therefore, recovering the image from the recorded intensity or magnitude diffraction pattern is a fundamental nonlinear imaging inverse problem. In the past several decades, various efforts to address this problem were developed. The PR algorithms can be classified into the following types: alternating projection (AP) methods [6–9], methods based on semidefinite program (SDP) [10,11], nonlinear compressed sensing (NLCS) methods [12,13], Wirtinger gradient descent methods [14,15] and regularized PR methods [16–19]. The AP methods start from an initial guess, and project the estimated image onto the measurement domain constraint and object domain constraint sets alternatively. The measurement domain constraint set guarantees ∗

Corresponding author. E-mail addresses: [email protected] (B. Shi), [email protected] (Q. Lian), [email protected] (X. Fan). https://doi.org/10.1016/j.sigpro.2019.04.013 0165-1684/© 2019 Published by Elsevier B.V.

that the reconstructed image matches with the measured data. The object domain constraint set is flexible, and one can formulate this constraint set by using the priors such as sparsity, non-negative support, and object intensity ranging. The AP methods are simple. However, since the formulated PR problems are non-convex, the AP methods often suffer from the stagnation. Translating the nonconvex problem into a convex one and solving the corresponding convex problem by using the convex relaxation technique are the core ideas of the SDP-based methods [10,11]. Unfortunately, SDPbased methods are not suitable for large scale imaging problems. The NLCS method that addresses the nonlinear compressed sensing problem is the extension of the linear compressed sensing method [12,13]. Indeed, the PR problem belongs to a NLCS problem. Therefore, the NLCS methods can be utilized for solving PR problems. Recently, the Wirtinger gradient descent methods are popular in the PR regime. An elaborated initial guess is designed by this class of methods, and then the image is estimated by using a gradient descent method whose gradient direction is the Wirtinger gradient direction. The provable convergence of these algorithms leads them to become popular ones [14,15]. Due to the lack of priors, the Wirtinger gradient descent methods suffer from the low reconstruction quality when the measured data is highly incomplete. The regularized PR methods can overcome this limitation by using some reasonable regularization terms. The cost functions of these PR methods consist of two terms: a data fidelity term that ensures the consistency of the image with the measured data and

84

B. Shi, Q. Lian and X. Fan / Signal Processing 162 (2019) 83–96

a regularization term that imposes desirable properties onto the unknown image. The prior information, such as sparsity, smoothness and non-local similarity, can be utilized for the image recovery via the regularization term [16–19]. Moreover, the data fidelity term can be adjusted based on the measurement noise type. The aforementioned imaging methods of exploiting priors have beaten the ones without utilizing priors. Prior information plays an important role in image reconstruction methods. How to explore prior information to improve the imaging quality is still an important issue. Recently, the Plug-and-Play Prior (PPP) framework that can exploit image priors implicitly by using a denoiser was suggested to address this issue. The PPP framework can leverage the advanced denoiser to improve the quality of reconstructed images. A class of PPP frameworks is designed based on the structures of the solvers. Venkatakrishnan et al. [20] proposed the PPP framework firstly, and they showed that the subproblem with a regularization term in the alternating direction method of multipliers (ADMM) [21] solver could be solved by a suitable denoiser. The effectiveness of this method was validated for solving linear imaging problems, but how to demonstrate its convergence is a problem. To cope with this problem, Chan et al. [22] proposed a plug-and-play ADMM algorithm with provable fixed-point convergence. They demonstrated that the proposed plug-and-play ADMM method could converge to a fixed point under a continuation scheme for any denoising algorithm satisfying the asymptotic criteria. The aforementioned two frameworks are based on the ADMM solver, other solvers were developed recently. Kamilov et al. [4] replaced the proximal operator in the fast iterative shrinkage/thresholding algorithm (FISTA) solver [23] with a denoiser to formulate a PPP framework, and proposed a nonlinear inverse scattering approach. Metzler et al. [24] plugged a denoiser into the generalized approximate message passing (GAMP) solver for phase retrieval, and achieved high-quality reconstructed images. Heide et al. [25] replaced the proximal operator in the primal dual solver with the denoiser, and proposed a flexible camera image processing framework. Another strategy to design PPP frameworks is exploiting the explicit regularization term of a denoiser. When this type of PPP frameworks meets coded diffraction imaging problems, the blockmatching and 3D filtering (BM3D) denoiser [26] is preferred, examples like the sparse phase amplitude reconstruction (SPAR) algorithm proposed by Katkovnik et al. [27,28] and the constrained phase retrieval (ConPR) algorithm proposed by Shi et al. [29]. They imposed the sparsity under the BM3D frame [30] on the unknown image to formulate a diffraction imaging problem and solved the corresponding image denoising subproblem by the BM3D denoiser. Under the noisy scenario, experiments validate the effectiveness of their algorithms. According to their ideas of designing PPP frameworks, only the denoisers that solve the image denoising problem with explicit regularization terms can be plugged in their frameworks. For example, the BM3D denoiser can be regarded as solving an image denoising optimization problem that contains an explicit sparse regularization term under the BM3D frame. However, explicit regularization terms cannot be presented by many deep learning-based denoisers, but experiments demonstrate that the aforementioned PPP frameworks of exploiting these denoisers are effective. From this perspective, the aforementioned two PPP frameworks that exploit deep learning-based denoisers lose interpretability. Moreover, the PPP frameworks based on the structure of solvers generally lose interpretability as optimization problems. A strategy of addressing these issues was proposed by Romano et al. [31]. The underlying image and the residual image, i.e. the difference image between the underlying image and its filtered version, should be uncorrelated [31]. Based on this fact, Romano et al. designed a regularization term by denoising (RED) and applied the RED model for solving linear imaging inverse problems.

The main differences of these three types of plug-and-play models are as follows. The PPP framework based on the structure of a solver relies on an implicit regularizer. The SPAR and ConPR frameworks are formulated based on an explicit regularizer. However, the RED model builds an explicit regularizer. Differ from the spirits of the aforementioned PPP frameworks, we propose a novel plug-and-play framework called plug-and-play regularization (PPR) based on the sparse representation theory, and show its effectiveness in solving a concrete nonlinear imaging problem. Compared with the PPP framework in [20], the proposed PPR framework has three differences. (i) The PPR model essentially belongs to a regularization model, while the PPP framework in [20] is an image reconstruction framework; (ii) We derive the PPR model based on the sparse representation theory, while the PPP framework in [20] is designed based on the structure of the ADMM solver; (iii) The linear imaging problem is considered in the PPP framework in [20], while in this paper we consider the problem of recovering the image of interest from the recorded intensity diffraction pattern. We develop the PPR model based on sparse representation over a tight frame. In the iteration process, we perform sparse coding for the estimated image and its filtered version that is filtered by a denoiser. The relation of similarity between the underlying image and its filtered version is constructed via the same sparse coefficients under a tight frame. By doing this, the priors contained in a denoiser and the sparsity of the image over a tight frame can be incorporated into the image reconstruction. We formulate a nonlinear optimization problem based on the proposed PPR model. The problem is solved by using the alternating optimization strategy, and the image updating subproblem is attacked by the epigraph concept. The advantage of the epigraph method is avoiding the finely-tuned regularization parameter. Experiments show that the proposed algorithm plugged by the fast and flexible denoising convolutional neural network (FFDNet) [32] (one of the state-ofthe-art denoisers) outperforms the previous PR algorithms in terms of the reconstruction quality. The structure of this paper is as follows. To begin with, Section 2 reviews previous PPP frameworks for solving both linear imaging inverse problems and nonlinear ones. Moreover, our motivation and contributions are introduced in this section. Section 3 introduces the proposed PPR model. Section 4 first formulates the nonlinear imaging problem of exploiting the proposed PPR model, and then describes the proposed algorithm for solving the coded diffraction imaging problem in detail. We present our experimental simulations in Section 5. Finally, concluding remarks and directions for future research are presented in Section 6. 2. Relative works and our work 2.1. Plug-and-play priors for solving linear imaging inverse problems Imaging inverse problems that recover the images from their degraded measurements are fundamental problems in image processing. In the past several decades, the linear imaging inverse problems whose degraded operators are linear were studied by the researchers. The typical linear imaging inverse problems are image denoising, image super-resolution, image deburring, and image compressed sensing problems. Due to the linearity of the degraded operator, the linear measured process can be modeled as

y = x + n.

(1)

Here x represents the underlying image, and  is the degraded matrix that models the corrupting or sampling process of the underlying image. y is the corrupted or measured data, and n represents some noise. Given the measured data y, recovering the image x is an inverse process, thus, the corresponding problem

B. Shi, Q. Lian and X. Fan / Signal Processing 162 (2019) 83–96

is called an imaging inverse problem. For different imaging inverse problems, the degraded matrix  is different. Concretely, the degraded matrixes can be identity, down-sampling, blurring and under-sampled matrixes for image denoising, image superresolution, image deburring, and image compressed sensing problems, respectively. Among these imaging inverse problems, the image denoising problem is an old but classical problem. Recent advanced image denoising algorithms have achieved superior performance even touched the ceiling in terms of noise removal performance. Can we leverage these achievements to solve other imaging inverse problems? The PPP framework answers this question positively. Imaging inverse problems are ill-posed, therefore, priors about the underlying image are often needed. Since regularized image reconstruction algorithms can exploit priors via a reasonable regularizer, this class of algorithms becomes popular these years. Based on the maximum a posteriori (MAP) criterion, the imaging optimization problem can be formulated as

min l (x; y ) + λR(x ). x

(2)

Here l(x; y) is the data fidelity term that can be adjusted based on the type of the measurement noise. For example, the data fidelity term is defined as l (x; y ) = ||y − x||22 for Gaussian noise. R(x) is the regularization term that can impose some desirable properties on the recovered image. λ is the regularization parameter that controls the tradeoff between these two terms. Problem (2) can be solved by the splitting methods, and it can be rewritten as the following constrained problem

min l (x; y ) + λR(z ) s.t. x = z. x,z

(3)

The above problem can be solved by using the ADMM solver and the half quadratic splitting (HQS) solver. Specially, the subproblems of containing the regularization term in ADMM and HQS are

min z

ρ||x − (z − u )||22 + λR(z )

(4)

η||x − z||22 + λR(z )

(5)

and

min z

respectively. Here, u is the dual scale variable [21], and ρ , η are two parameters. Venkatakrishnan et al. [20] and Zhang et al. [33] exploited denoisers to solve the subproblems (4) and (5) respectively. In other words, the denoising step is one step of the ADMM or HQS solver. To exploit different priors for the image reconstruction, one can use different denoisers to solve the subproblems (4) and (5). Most image denoising algorithms such as deep denoisers (the denoiser based on deep learning) do not have explicit regularization terms. From the optimization model view, the aforementioned two PPP frameworks lose interpretability. Nevertheless, experiments for solving linear imaging inverse problems validate the effectiveness of these two PPP frameworks. 2.2. PPP model: From linear imaging to nonlinear imaging In several practical applications, the sampling operator is nonlinear. For example, in the coherent diffractive imaging field, the CCD camera only records the intensities of the arrived light, and the phase information is lost. Under this scenario, the sampling operator is nonlinear. Recovering the underlying image from the measurements without phase information is a typical nonlinear imaging inverse problem. Without loss of generality, the nonlinear sampling process can be described as

y = H (x ) + n

(6)

where H(•) is a nonlinear sampling operator that models the response of the general nonlinear sampling system. Recovering the

85

image x from the nonlinear measured data y is the goal of the nonlinear imaging technique. Kamilov et al. [4] proposed a PPP approach for nonlinear imaging, and they formulated the following optimization problem

min x

||y − H(x )||22 .

(7)

They [4] solved the above problem by using the FISTA solver, but they replaced the proximity operator in FISTA with a denoiser. They showed that this method could recover images from the measurements of scattered waves. In previous work, we [29] gained an insight into the relation between the regularized PR method and the AP method, and incorporated the regularization term of exploiting the BM3D frame into problem (7). We proposed the so-called ConPR framework, and achieved high-quality reconstructions. Indeed, the ConPR framework can be regarded as a PPP framework. Metzler et al. [24] exploited the GAMP solver to solve the nonlinear imaging inverse problem, and incorporated the denoiser into the framework of GAMP. By employing the BM3D denoiser, they have achieved dramatic improvement in compressive phase retrieval performance. 2.3. Motivation and contributions The aforementioned PPP frameworks have achieved good performance, but they suffer from either the certain structure of the solver or limiting on the denoiser without an explicit regularization term. The RED framework is an attempt to solve these issues, but the denoiser in that framework should admit some conditions [31]. Can we design a general effective plug-and-play regularization model that is suitable for any denoiser? This paper answers this question positively. Designing an effective regularizer with the plug-and-play idea is our goal in this work. To that end, an important problem is how to model the relation of the underlying image and its filtered version. The main inherent structures of the underlying image should be similar or same as those of its filtered version. Therefore, these two images should have similar coefficients over a same tight frame. Based on this knowledge, we propose a novel plug-and-play regularization (PPR) model. We exploit the proposed PPR model for solving nonlinear imaging inverse problems. Our main contributions are summarized as follows. (1) We propose a plug-and-play regularization called PPR model. In PPR, the similarity between the underlying image and its filtered image is enforced in the sparsifying transform domain. Specially, the tight frame is selected as the sparsifying transform. The sparsity of the underlying image under a tight frame can be utilized for the image reconstruction. Additionally, the tight property of the frame is beneficial to solving the image updating problem. In the proposed PPR model, any effective denoiser can be incorporated into it to exploit corresponding priors implicitly. The proposed PPR model is a regularization model, and it can be applied to solve both linear and nonlinear imaging inverse problems. (2) We propose a flexible plug-and-play nonlinear imaging framework by using the proposed PPR model. Combining the PPR model and the data fidelity term, we formulate a general nonlinear imaging optimization problem. The alternating iteration strategy that consists of three steps, i.e. filtering step, sparse coding step and image updating step, is utilized for solving this problem. To avoid the finely-tuned parameter, the epigraph concept is utilized to attack the image updating subproblem. One can obtain a satisfactory solution to the formulated problem by using the alternating iteration strategy. The flexibility of the proposed nonlinear imaging framework contains the following two folds. (i) The data fidelity term is flexibility, and one can define it based on applications; (ii) The denoiser is flexibility, and one can leverage the recent deep denoisers to improve the recovery performance.

86

B. Shi, Q. Lian and X. Fan / Signal Processing 162 (2019) 83–96

(3) We consider the coded diffraction imaging problem, and propose an effective diffraction imaging algorithm. The BM3D and FFDNet denoisers are plugged into the proposed diffraction imaging framework. Experiments show that the proposed diffraction imaging algorithms outperform previous methods in terms of the reconstruction quality. Moreover, the proposed diffraction imaging algorithms have good convergence behaviors, and deliver the competitive advantage of the imaging quality under the undersampled data scenario. 3. Plug-and-play regularization via sparse representation

||α||0 ≤ T

(8)

where α is the coefficient vector of the image x under the tight frame W. ||•||0 represents the l0 pseudonorm that counts the number of non-zero elements in a vector. T represents the sparsity level. The underlying image and its filtered version have similar structures, therefore, their coefficients over the same tight frame are similar. We enforce that the coefficient vectors of these two images approximate the same sparse vector α. The relation function of the image and its filtered version can be described as

Re[x; D(x )] = min α

4. The proposed PPR model for solving nonlinear imaging inverse problems 4.1. Problem formulation

The key problem of designing a plug-and-play regularization is how to model the relation between the underlying image and its filtered version. In this work, we model this relation via a sparse representation model, and employ a tight frame as the sparsifying transform. Given a tight frame W that admits WT W = I, the sparse representation model is

Wx = α s.t.

• Versatility: Any effective denoiser can be incorporated into the proposed PPR model. Moreover, this regularization model is suitable for solving any imaging inverse problem. • Effectiveness: Multiple priors can be utilized for recovering images. The priors can be utilized implicitly via the plugged denoiser. Additionally, the sparsity under a tight frame is also utilized, and it is beneficial to the image reconstruction quality.

||Wx − α||22 + λ||WD(x ) − α||22 + β||α||0 (9)

where D(x) represents the image filtered by a denoiser D(•). Re[x; D(x )] is an envelope function that can characterize the relation of the image x and the filtered image D(x). The first term of the right hand of (9) enforces that the coefficient vector of the image x over a tight frame approximates a sparse vector α, while the second term enforces that the coefficient vector of the filtered image D(x) over the same frame approximates the same sparse vector α. Indeed, any reasonable sparsifying frame can be utilized for defining the relation function. The advantage of a tight frame is that it can simplify the computation of updating the image, compared to a redundant sparsifying transform. The third term of the right hand of the Eq. (9) is the sparsity penalty term for the vector α. The relation function (9) can be minimized to enforce the similarity of the two images, and it can be regarded as a regularization term: R(x ) = Re[x; D(x )]. In this term, the denoiser D(•) plays an important role for high-quality imaging. Different denoisers indicate that different priors are utilized. For example, the BM3D denoiser is selected for filtering, therefore, the non-local similarity and the sparsity under a 3D sparsifying transform can be utilized for image reconstruction. If we plug the total variation (TV) denoiser into the model (9), the gradient sparsity can be utilized for image reconstruction. To further improve the reconstruction quality, one can exploit multiple denoisers for the model (9). Under this scenario, the filtered image D(x) is a linear combination of dif ferent filtered images, namely D(x ) = i λi Di (x ). Here, Di (x) is the image filtered by the ith denoiser. 0 ≤ λi ≤ 1 is the weighting factor that controls the contribution of the ith denoiser to the filtered im age D(x), and it admits i λi = 1. For color image reconstruction, the cross-channel gradient correlation [34] can be utilized by replacing the denoiser with the defined cross-channel prior-induced proximal operator in [25,35]. The advantages of the proposed plug-and-play regularization (PPR) model are summarized as follows:

For nonlinear imaging, the prior information also plays an important role for recovering high-quality images. The priors utilized by a regularization term are popular. In this work, we exploit the proposed PPR model for solving general nonlinear imaging inverse problems, and formulate the optimization problem as follows

min x,α

||Wx − α||22 + λ||WD(x ) − α||22 + β||α||0 + γ l[H(x ); y] (10)

where γ is a parameter that controls the weighting of the data fidelity term in the optimization problem. In (10), the first term indicates that the coefficient vector of the underlying image over a frame is similar to a sparse vector α, and the second term indicates that the coefficient vector of the filtered image under the same frame is also similar to the same sparse vector α. The sparsity of the vector α is enforced by the third term. The last term of the cost function in (10) is the nonlinear data fidelity term. The nonlinear sampling operator in this term is general, which spreads the applications of the proposed imaging framework. According to the MAP criterion, the data fidelity term can be designed based on the noise type, i.e., l[H(x ); y] = ||y − H(x )||22 for Gaussian noise, l[H(x ); y] = 1T H(x ) − yT ln H(x ) for Poisson noise, and l[H(x ); y] = ||y − H(x )||1 for outlier. Here, 1 denotes a column vector whose elements are all ones, and ||•||1 represents the l1 norm. In this work, we consider recovering the images from nonlinear data corrupted by Gaussian noise. In the reconstruction process, the image x is unknown, instead, the estimated image is filtered by the denoiser. For the tth iteration, the regularized nonlinear imaging optimization problem under the Gaussian noise case can be described as

min ||Wx − α||22 + λ||WD(x(t−1) ) − α||22 + β||α||0 + γ ||y −H(x )||22 . x,α

(11) The above problem is a non-convex problem, addressing it effectively is also one of our main contributions. Next subsection will introduce the method to solve it. 4.2. Plug-and-play nonlinear imaging method We exploit the alternating optimization strategy to solve the problem (11). The three steps to solve the problem (11) are as follows (for the tth iteration): (1) Filtering step In this step, we exploit a denoiser to filter the estimated image, and obtain the filtered image D(x(t−1 ) ). For a certain denoiser, the noise standard deviation σ of the estimated image x(t−1 ) is often needed. In this paper, the noise standard deviation σ is evaluated by using the robust median operator [36]. The denoiser is crucial for recovering high-quality images. Recent years have witnessed a great progress of the image denoising algorithms. Among these algorithms, the BM3D image denoising algorithm is popular in the last decade. More recent years, the deep priors further unroll the inherent features of the structures in an image, and the

B. Shi, Q. Lian and X. Fan / Signal Processing 162 (2019) 83–96

image denoising algorithms based on deep learning have achieved remarkable performance in terms of the noise removal. In the experimental section, we mainly discuss two advanced denoisers, i.e., the BM3D and FFDNet denoisers, for nonlinear imaging. (2) Sparse coding step Fixing the estimated image x and its filtered image D(x), we obtain the corresponding subproblem

min α

||Wx(t−1) − α||22 + λ||WD(x(t−1) ) − α||22 + β||α||0 .

(12)

The above problem is a sparse coding problem, which can be rewritten as

min α

||W

x(t−1) + λD(x(t−1) ) − α||22 +  2 ||α||0 . 1+λ

(13)

β Here  2 = 1+ λ . Problem (13) can be solved by the hard thresholding operator



α

(t )

x(t−1) + λD(x(t−1) ) = Hard W , 1+λ



(14)

where Hard(•) is the hard thresholding operator that processes a vector element-wisely. The hard thresholding operator is defined as: Hard(u,  ) = u, if |u| >  , and Hard(u,  ) = 0 otherwise. (3) Image updating step In this step, problem (11) is solved with fixed α. This corresponds to the subproblem

min x

||Wx − α(t ) ||22 + γ ||y − H(x )||22 .

(15)

Since the frame W is tight, problem (15) can be simplified as

min x

||x − WT α(t ) ||22 + γ ||y − H(x )||22 .

(16)

As we know, in nonlinear sampling systems, the regularization parameter γ is sensitive to many factors, such as sampling ratios and noise levels. Fine-tuning of this parameter is a difficult task. To avoid this finely-tuned parameter, we exploit the epigraph concept to solve the above problem. Problem (16) can be translated into a constrained optimization problem

min x

||x − WT α(t ) ||22 s.t. ||y − H(x )||22 ≤ q

(17)

where q is a small positive value. Let f (x ) = ||y − H(x )||22 be the data fidelity function. We define its epigraph set as follows [37]

E = {x = [xT q]T : q ≥ f ( x )}.

(18)

For an N dimensional vector x, the elements in the defined epigraph set E are column vectors with N+1 dimension, whose the (N+1)th component is greater than f(x). Based on the epigraph concept, we can attack problem (17) by solving the following projection problem [37]

min ||x − z||22

(19)

x∈E

where z = [(WT α(t ) )T 0]T . Problem (19) is a projection problem, and one can obtain its solution by projecting z onto the supporting hyperplane of the epigraph set E iteratively. After J times projection, the final orthogonal projection onto the supporting hyperplane of the epigraph set is

[vTJ f (vJ )]T where vJ = vJ−1 −

The image updating step can be regarded as a gradient descent step whose step size is computable. If the analytical derivative of the data fidelity function does not exist at some points, one can calculate its sub-derivative at these points, and exploit the sub-gradient for image updating. In a concrete sampling model, the sub-derivative of the data fidelity function for different noise types often exists at each point. Therefore, the proposed model can be suitable for different noise types by adjusting the data fidelity function. Under some cases, the sub-gradient does not exist. For example, when dealing with real-valued bounded images, one may add a box constraint {x|x ∈ [0, 1]N } utilized in [19] into Eq. (15) to restrict the pixel values to a reasonable range. Under this scenario, the gradient descent method is replaced by the proximal gradient method. Moreover, under the over-exposed case [38], saturated pixels exist in the nonlinear measured data sampled by CCD or other devices. The sampling model can be written as y = B[H(x )] + n. Here, B[•] = min(•, δ ) is a clip function, and δ is the maximum pixel value. Under this case, one can incorporate an auxiliary variable u = H(x ) into Eq. (15), and solve the constrained optimization problem by using the ADMM method or the variable splitting method suggested by Mosleh et al. [39]. For these two frameworks, the x subproblem is similar as Eq. (15), and it can be attacked by using the proposed image updating method. The u subproblem is a simple element-wise model, which can be solved by element-wise method [39]. So far, all issues to solve the formulated problem (11) are addressed. We can iterate the filtering step, the sparse coding step and the image updating step alternately for solving problem (11). One can exploit the advanced denoiser to perform the filtering step. The solution of the sparse coding subproblem is given. The image updating subproblem is solved by the epigraph method, and an approximated solution is achieved. We summarize the proposed imaging algorithm based on the proposed PPR model in Table 1. 4.3. Plug-and-play regularization for coded diffraction imaging Recovering the image of interest from the nonlinear measurements, whose phase information is lost, is a typical nonlinear imaging inverse problem. In this paper, we consider the coded diffraction imaging problem, namely the problem of imaging from the recorded intensity diffraction patterns. An illustrative setup for acquiring the nonlinear measurements in a coded diffraction imaging system is presented in Fig. 1. A light emitted from a laser illuminates the object of interest that is denoted x ∈ RN (we only consider the real-valued images in this paper). A random phase plate is placed behind the object to modulate the object. The modulation to the object can redistribute the spectrum of the object. After the transmission of the light that carries the modulated object information in the space, the diffraction pattern is recorded by the charge coupled device (CCD) or any other lighting apparatus. However, due to the limitation of the sampling devices such as

(20) f (vJ−1 ) ∇ f (vJ−1 ). Here, ||∇ f (vJ−1 )||22 +1

∇ f (vJ−1 ) repre-

sents the derivative of the function f at the point vJ−1 . The successive projection starts from [vT0 0]T where v0 = WT α(t ) , and arrives at a projected point [vTJ f (vJ )]T after J times projection. By doing so, one can obtain a satisfactory solution vJ to problem (17). In this step, only the maximum iteration needs to be specified, which avoids the finely-tuned regularization parameter.

87

Fig. 1. An illustrative realistic setup of the coded diffraction pattern model.

88

B. Shi, Q. Lian and X. Fan / Signal Processing 162 (2019) 83–96 Table 1 Nonlinear imaging algorithm based on the PPR model. 1. Input: Nonlinear measurement y. 2. Initialization: The thresholding value  , initial image x(0) and initial filtered image D(x(−1) ) = x(0) , maximum iteration tmax , J. 3. Outer Loop: For t=1, 2,..., tmax do 4. Noise standard deviation estimating: Obtain the estimated noise standard deviation σ by using the robust median operator. 5. Filtering: Obtain the filtered image D(x(t−1) ). 6. Sparse coding: Compute the sparse coefficients α(t) via (14). 7. Set the initial point v0 = WT α(t ) . 8. Inner loop: For j=1, 2,..., J do % image updating step 9. Compute the gradient or sub-gradient ∇ f (v j−1 ). f (v

)

10. Obtain the estimated image by computing v j = v j−1 − ||∇ f (v j−1)||2 +1 ∇ f (v j−1 ). j−1 2 11. End inner loop 12. Let x(t ) = vJ , and output the recovered image if t = tmax or

||D(x(t−1) )−D(x(t−2) )||22 ≤ 10−4 . ||D(x(t−1) )||22

13. End outer loop 14. Output: The recovered image.

Fig. 2. Testing images. (a)–(d): Lena; barbara; hill; boat. (e)–(h): Couple; fingerprint; acinar-cell; chromaffin-cell.

CCD, only the intensity of the diffraction pattern is acquired. Therefore, the sampling process in the coded diffraction imaging system is nonlinear. Mathematically, the sampling operator for the single observation in this nonlinear imaging system can be modeled as

H(x ) = |Ax|2 = |F(M  x )|2

(21)

where M represents the illumination mask (phase plate) whose elements are selected randomly from the set { ± 1, ± i} and  represents the element-wise product. Here, A is a linear matrix that models the propagation of the light. F is the Fourier transform matrix. For the undersampled scenario, the nonlinear sampling operator is

H(x ) = |As x|2 = |SF(M  x )|2

(22)

where As = SA and S ∈ RK×N (here K < N stands for the number of the measurements) is a diagonal matrix whose diagonal elements K are 0 or 1. The sampling ratio is defined as ratio = N × 100%. To exploit the proposed nonlinear imaging algorithm for coded diffraction imaging, the gradient of the function f (x ) = ||y − |As x|2 ||22 should be calculated. Indeed, when S is an identity matrix, the nonlinear sampling operator under the undersampled scenario equals to the sampling operator of the single observation model.

Since the single observation model is a special case of the undersampled model, we only present the gradient under the undersampled case. The gradient of the function f(x) is as follows

∇ f (x ) = 4 × real {AHs [As x  (|As x|2 −y )]}.

(23)

1 H Here AH s x = K M  [ (SF ) x], where M represents the dual filter of the illumination mask and (SF)H represents the conjugate transpose matrix of the product matrix SF. real(•) is the operator that extracts the real part of a complex-value vector. Once obtaining the gradient of the function f(x) in the coded diffraction pattern model, one can use the proposed nonlinear imaging framework for coded diffraction imaging.

5. Experimental simulations In this section, various numerical experiments are conducted to study the effectiveness of the proposed nonlinear imaging framework and the proposed PPR model. To that end, we consider the recovery of real-valued images from only one noisy coded diffraction pattern. Even, recovering the image from undersampled data is discussed. First, we introduce the experimental setup and implementation details. Second, we evaluate the proposed algorithms for the recovery of images from single observation at various noise

B. Shi, Q. Lian and X. Fan / Signal Processing 162 (2019) 83–96

89

Fig. 3. The lena images recovered by the seven imaging algorithms at the single observation and SNR = 15 dB case. From left to right and top to bottom: (a) the original image; the image recovered by (b) DOLPHIn (PSNR = 18.16 dB, FSIM = 0.9641); (c) BM3D-prGAMP (PSNR = 28.99 dB, FSIM = 0.9946); (d) SPAR (PSNR = 29.34 dB, FSIM = 0.9938); (e) ConPR-BM3D (PSNR = 30.81 dB, FSIM = 0.9965); (f) ConPR-FFDNet (PSNR = 31.81 dB, FSIM = 0.9972); (g) PPR-BM3D (PSNR = 31.69 dB, FSIM = 0.9969); (h) PPR-FFDNet (PSNR = 32.34 dB, FSIM = 0.9972).

Fig. 4. The boat images recovered by the seven imaging algorithms at the single observation and SNR=10 dB case. From left to right and top to bottom: (a) the original image; the image recovered by (b) DOLPHIn (PSNR = 14.27 dB, FSIM = 0.9259); (c) BM3D-prGAMP (PSNR = 25.00 dB, FSIM = 0.9846); (d) SPAR (PSNR = 24.31 dB, FSIM = 0.9838); (e) ConPR-BM3D (PSNR = 25.58 dB, FSIM = 0.9895); (f) ConPR-FFDNet (PSNR = 25.56 dB, FSIM = 0.9902); (g) PPR-BM3D (PSNR = 26.73 dB, FSIM = 0.9888); (h) PPR-FFDNet (PSNR = 27.27 dB, FSIM = 0.9890).

levels. To show the effectiveness of the proposed imaging algorithms, previous algorithms of exploiting prior information are selected as the benchmark algorithms. Third, we perform the undersampled phase retrieval, and evaluate the testing algorithms based on the reconstruction quality. Finally, the convergence behavior and the running time are discussed.

5.1. Experimental setup and implementation details In this paper, we sample the 512 × 512 images by using the coded diffraction pattern model with the quaternary mask. The original images are presented in Fig. 2, and they can be downloaded from https://github.com/shibaoshun/ConPR. The image pixel

90

B. Shi, Q. Lian and X. Fan / Signal Processing 162 (2019) 83–96

Fig. 5. The chromaffin-cell images recovered by the seven imaging algorithms at the single observation and SNR=10 dB case. From left to right and top to bottom: (a) the original image; the image recovered by (b) DOLPHIn (PSNR = 12.41 dB, FSIM = 0.9120); (c) BM3D-prGAMP (PSNR = 19.71 dB, FSIM = 0.9860); (d) SPAR (PSNR = 18.46 dB, FSIM = 0.9795); (e) ConPR-BM3D (PSNR = 19.15 dB, FSIM = 0.9870); (f) ConPR-FFDNet (PSNR = 19.43 dB, FSIM = 0.9874); (g) PPR-BM3D (PSNR = 20.56 dB, FSIM = 0.9889); (h) PPR-FFDNet (PSNR = 20.95 dB, FSIM = 0.9897). Table 2 The average results (PSNR/FSIM) of the eight testing images obtained by the testing algorithms at the single observation scenario. Algorithms

Noise levels 5

10

15

20

25

DOLPHIn BM3D-prGAMP SPAR ConPR-BM3D ConPR-FFDNet PPR-BM3D PPR-FFDNet

12.28/0.8939 21.52/0.9819 20.23/0.9789 21.05/0.9832 21.38/0.9856 23.01/0.9860 23.65/0.9864

13.87/0.9364 24.09/0.9903 23.36/0.9889 24.53/0.9928 24.94/0.9937 25.49/0.9934 25.99/0.9937

17.19/0.9770 25.97/0.9949 25.98/0.9936 27.47/0.9971 28.06/0.9975 27.96/0.9972 28.35/0.9973

20.09/0.9883 28.07/0.9975 28.27/0.9968 29.95/0.9987 30.50/0.9988 30.25/0.9987 30.56/0.9988

21.69/0.9922 30.25/0.9988 30.42/0.9986 32.21/0.9994 32.47/0.9994 32.44/0.9994 32.51/0.9994

in each real-valued image ranges from 0 to 1. To evaluate the robustness of the proposed imaging algorithms, we corrupt the measurements with additive white Gaussian noise. The noise level is controlled by signal-to-noise ratio (SNR) that is defined as

SNR = 10log10 [||y¯ ||22 /||y¯ − y||22 ] dB.

|SF(M  x )|2

(24)

Here y¯ = (x is the original image) represents the intensity of the theoretical diffraction pattern. To demonstrate the effectiveness of the proposed algorithms, the recent imaging algorithms that exploit image prior information are selected as the benchmark algorithms. Concretely, the DOLPHIn [19], BM3D-prGAMP [24], SPAR [27], and ConPR [29] algorithms are selected as the benchmark algorithms. The sparsity over an adaptive synthesis dictionary was utilized for imaging in the DOLPHIn algorithm, and it optimizes the dictionary and the image simultaneously from the recorded coded diffraction patterns. In this paper, the default parameters of the original DOLPHIn algorithm are utilized. The BM3D image denoiser is incorporated into the GAMP algorithm in the BM3D-prGAMP algorithm, which exploits the prior information implicitly. The sparse model under the BM3D frame is exploited for image reconstruction in the SPAR algorithm, and its optimization model is formulated by using a variational approach. Both the super-resolution phase retrieval and the coded diffraction

imaging demonstrate the effectiveness of the SPAR algorithm. The BM3D frame is also utilized in the ConPR algorithm, but it can exploit advantages of both the regularization methods and the alternating projection methods. The public codes of DOLPHIn can be found in the DOLPHIn package,1 and we use its default parameters. The concrete parameters for BM3D-based imaging methods, i.e., the BM3D-prGAMP, SPAR and ConPR algorithms, are as follows. The maximum iteration for BM3D-prGAMP is 70, for SPAR is 50 and for ConPR is 20. The default parameters suggested in their public codes are utilized. We use the BM3D-prGAMP algorithm from the D-AMP_Toolbox package.2 The code of the SPAR algorithm can be downloaded from the SPARSE project homepage.3 The code of the ConPR algorithm is also public.4 The concrete implementations of the proposed algorithm are as follows. We use the learned tight frame in [40] as the sparsifying

1

https://www.graphics.rwth-aachen.de/media/resource_files/DOLPHIn.zip. https://github.com/ricedsp/D-AMP_Toolbox/tree/master/Algorithms. 3 http://www.cs.tut.fi/sgn/imaging/sparse/Demo- Sparse- Phase- Retrieval- JAN_ 2018.zip. 4 https://github.com/shibaoshun/ConPR. 2

B. Shi, Q. Lian and X. Fan / Signal Processing 162 (2019) 83–96

91

Fig. 6. The hill images recovered by the seven imaging algorithms at the sampling ratio 50% and SNR = 10 dB case. From left to right and top to bottom: (a) the original image; the image recovered by (b) DOLPHIn (PSNR = 14.31 dB, FSIM = 0.9166); (c) BM3D-prGAMP (PSNR = 24.79 dB, FSIM = 0.9816); (d) SPAR (PSNR = 24.46 dB, FSIM = 0.9774); (e) ConPR-BM3D (PSNR = 26.10 dB, FSIM = 0.9833); (f) ConPR-FFDNet (PSNR = 24.96 dB, FSIM = 0.9434); (g) PPR-BM3D (PSNR = 26.40 dB, FSIM = 0.9834); (h) PPR-FFDNet (PSNR = 26.49 dB, FSIM = 0.9817).

Fig. 7. The couple images recovered by the seven imaging algorithms at the sampling ratio 50% and SNR=20 dB case. From left to right and top to bottom: (a) the original image; the image recovered by (b) DOLPHIn (PSNR = 16.31 dB, FSIM = 0.9435); (c) BM3D-prGAMP (PSNR = 28.10 dB, FSIM = 0.9949); (d) SPAR (PSNR = 27.84 dB, FSIM = 0.9947); (e) ConPR-BM3D (PSNR = 28.55 dB, FSIM = 0.9953); (f) ConPR-FFDNet (PSNR = 26.58 dB, FSIM = 0.9876); (g) PPR-BM3D (PSNR = 29.73 dB, FSIM = 0.9969); (h) PPR-FFDNet (PSNR = 29.93 dB, FSIM = 0.9970).

transform for the PPR model. One can exploit other tight frames for the image reconstruction, but this is not our focus in this paper. We set the maximum iteration of the outer loop tmax = 30, and the inner loop J = 2. For the input noise deviation of the image denoiser, the noise standard deviation σ is evaluated by using the robust median operator. The weighting factor λ controls

the contribution of the filtered image to the sparse codes. Indeed, this parameter controls the contribution of the priors to the final recovered image. The parameter β in problem (11) is fused into the thresholding value  . The parameter  should be small to preserve detail information when the estimated image gets better. As the growth of the iterations, the estimated image becomes

92

B. Shi, Q. Lian and X. Fan / Signal Processing 162 (2019) 83–96

Fig. 8. The acinar-cell images recovered by the seven imaging algorithms at the sampling ratio 50% and SNR=10 dB cases. From left to right and top to bottom: (a) the original image; the image recovered by (b) DOLPHIn (PSNR = 10.11 dB, FSIM = 0.9090); (c) BM3D-prGAMP (PSNR = 18.68 dB, FSIM = 0.9899); (d) SPAR (PSNR = 18.00 dB, FSIM = 0.9898); (e) ConPR-BM3D (PSNR = 18.83 dB, FSIM = 0.9857); (f) ConPR-FFDNet (PSNR = 17.62 dB, FSIM = 0.9576); (g) PPR-BM3D (PSNR = 19.37 dB, FSIM = 0.9906); (h) PPR-FFDNet (PSNR = 19.50 dB, FSIM = 0.9910). Table 3 The average results (PSNR/FSIM) of the eight testing images obtained by the testing algorithms at the undersampled data scenario. Algorithms

DOLPHIn BM3D-prGAMP SPAR ConPR-BM3D ConPR-FFDNet PPR-BM3D PPR-FFDNet

Sampling ratio 30%

Sampling ratio 50%

Sampling ratio 80%

SNR = 10

SNR = 20

SNR = 10

SNR = 20

SNR = 10

SNR = 20

12.79/0.8973 21.37/0.9768 20.96/0.9769 21.88/0.9790 20.46/0.9228 22.50/0.9806 22.70/0.9798

13.93/0.9173 25.51/0.9930 24.51/0.9930 23.19/0.9888 21.61/0.9590 25.40/0.9935 25.63/0.9936

13.03/0.9149 22.66/0.9845 22.28/0.9827 23.53/0.9854 21.91/0.9503 23.82/0.9869 24.08/0.9863

14.84/0.9437 26.54/0.9954 26.51/0.9953 26.44/0.9945 24.22/0.9830 27.87/0.9969 28.01/0.9970

13.25/0.9239 23.66/0.9886 23.10/0.9872 24.52/0.9891 23.02/0.9686 24.94/0.9908 25.21/0.9904

17.41/0.9761 27.34/0.9969 27.76/0.9964 28.25/0.9971 26.38/0.9932 29.43/0.9982 29.48/0.9982

better and better, thereby, the value of the noise standard deviation σ gets smaller and smaller. Therefore, it is reasonable to set the thresholding value  as  = C σ . Keeping the aforementioned facts in mind, we tune a parameter λ or C at a time while fixing the other one. Heuristically, the parameter C in the proposed imaging algorithm is set as follows: C = 2 for single observation and C = 3 for undersampled observation. The weighting factor λ is set as λ=1.5. All the tests are performed on a desktop computer with a Core i7-7700 CPU, 3.6 GHz processor, 8GB of RAM, and running the Windows 10 64-bit operating system and MatLab 2017a. For fair comparisons, all the algorithms exploit the same initial random image. Moreover, the initial estimated image is the same at different sampling ratios and different noise levels. 5.2. Imaging from single observation Multiple diffraction patterns contain more information about the underlying image compared to single diffraction pattern. However, recording multiple diffraction patterns is more complex than that of single observation. Additionally, processing multiple diffraction patterns is time-consuming. Therefore, imaging from single observation is practical but challenging. In this subsection, the pro-

posed algorithms and the benchmark algorithms were simulated in the case of the single observation. To demonstrate the effectiveness of the proposed plug-and-play framework, the BM3D denoiser was plugged into the proposed plug-and-play framework, and the corresponding algorithm termed PPR-BM3D. The BM3D denoiser that is exploited in ConPR [29] was plugged into the proposed PPR model. Moreover, the FFDNet denoiser that can produce denoised images with high quality was also utilized in the proposed framework, and we term this algorithm as PPR-FFDNet. To further show the effectiveness of the proposed PPR model, we also plugged the FFDNet denoiser into the ConPR framework. We exploit the same parameters for the proposed imaging algorithms plugged by different image denoisers in the cases of various noise levels and sampling ratios. We evaluate the testing algorithms based on PSNR (Peak Signal to Noise Ratio) and FSIM (Feature SIMilarity) [41]. For PSNR, larger values are better, and for FSIM, a higher FSIM value indicates the better visual quality. Table 2 presents the PSNR and FSIM values of the seven testing algorithms at various noise levels. The averaged values of the eight testing images are reported. As can be seen from the table, for most parts, the proposed PPR-FFDNet algorithm achieves the highest PSNR and FSIM values. Deep priors and sparsity are utilized in the PPR-FFDNet algorithm, therefore, it can

B. Shi, Q. Lian and X. Fan / Signal Processing 162 (2019) 83–96

93

Fig. 9. Convergence behaviors of the proposed imaging algorithms at the single observation and SNR = 15 dB cases. From top to bottom: the curve of PSNR vs Iteration achieved by (a) PPR-BM3D; (b) PPR-FFDNet.

obtain the highest reconstruction quality. At SNR = 5 dB, the average PSNR values achieved by the PPR-FFDNet algorithm are nearly 11.37 dB, 2.13 dB, 3.42 dB, 2.60 dB and 2.27 dB higher than those of the DOLPHIn, BM3D-prGAMP, SPAR, ConPR-BM3D and ConPRFFDNet algorithms respectively. The high-quality reconstructions achieved by the proposed PPR-FFDNet algorithm at low SNRs indicate the robustness of the proposed algorithm. Among the BM3D-based imaging algorithms, the PPR-BM3D algorithm is the best in terms of PSNR and FSIM values. This result demonstrates the effectiveness of the proposed PPR model. Compared to the ConPR-FFDNet algorithm, the better reconstruction quality achieved by the proposed PPR-FFDNet algorithm indicates that the sparsity over a tight frame is beneficial to image reconstruction. Among the seven testing algorithms, the DOLPHIn algorithm is the worst in terms of the reconstruction quality. The reason lies in the fact that only the sparsity under an adaptive dictionary is utilized in the DOLPHIn algorithm. However, the BM3D-based algorithms can exploit the non-local similarity and

the sparsity in the 3D transform domain implicitly. Deep priors are utilized in the ConPR-FFDNet algorithm and the PPR-FFDNet algorithm. Figs. 3–5 present some images recovered by the proposed algorithms and the benchmark algorithms. As can be seen from these reconstructions, the advantage of the proposed PPR-FFDNet algorithm is obvious both visually and numerically. 5.3. Imaging from undersampled data Imaging from undersampled nonlinear data is a challenging problem, and a typical application is the coherent diffractive imaging where undersampled observation is obtained. In this subsection, we evaluate the proposed algorithms and the benchmark algorithms on the undersampled data. We generate undersampled observations by randomly down-sampling the coded diffraction pattern of the image at sampling ratios 30%, 50%, and 80%. Gaussian noise whose noise levels are SNR = 10 dB and 20 dB is added

94

B. Shi, Q. Lian and X. Fan / Signal Processing 162 (2019) 83–96

Fig. 10. Convergence behaviors of the proposed imaging algorithms at the sampling ratio 80% and SNR = 20 dB cases. From top to bottom: the curve of PSNR vs Iteration achieved by (a) PPR-BM3D; (b) PPR-FFDNet.

on the clean undersampled data. Table 3 presents the average results of the eight testing images obtained by the seven imaging algorithms. One can see that the proposed algorithms always outperform the benchmark algorithms in terms of the average PSNR and FSIM values. The high-quality reconstructions achieved by the proposed algorithms verify the effectiveness of the proposed algorithms. For most parts, the proposed PPR-FFDNet algorithm is the best. The proposed PPR-FFDNet algorithm is slightly better than the PPR-BM3D algorithm in terms of the average PSNR values, but better than the benchmark algorithms, especially at the low SNR cases. Therefore, the robustness to noise at the undersampled scenario is also validated. Figs. 6–8 present the hill, couple and acinar-cell images recovered by the seven testing algorithms. From the figures, the images recovered by the DOLPHIn algorithm contain many noise components, and its reconstructions are the worst. The BM3D-prGAMP and SPAR algorithms remove these noise components, but their re-

constructions are so blurring that some edge information of the reconstructions is lost. For example, in Fig. 6 (c) and (d), due to the missing of the edge information, one cannot identify the person on the ground. The proposed algorithms and the ConPR algorithms preserve most of the edge information. However, the reconstructions of the ConPR-BM3D and ConPR-FFDNet algorithms contain unwanted artifact components. The PPR-BM3D algorithm reduces some artifact components, but its reconstructions still contain a little artifact information. For the visual quality, the proposed PPRFFDNet algorithm can achieve better reconstructions than those of the benchmark algorithms, and these results indicate the effectiveness of the proposed PPR model. 5.4. Convergence behaviors The formulated problem in this paper is a non-convex problem. We attack this problem via alternating between the

B. Shi, Q. Lian and X. Fan / Signal Processing 162 (2019) 83–96 Table 4 Comparison of the running time (s) obtained by different algorithms. For different observation, we report the mean running time of the eight images. Algorithms

Single observation

Undersampled observation

DOLPHIn BM3D-prGAMP SPAR ConPR-BM3D ConPR-FFDNet PPR-BM3D PPR-FFDNet

16.42 233.29 126.77 13.49 3.72 27.01 10.00

18.63 255.88 113.50 12.16 3.72 26.33 13.02

filtering, sparse coding and image-updating steps. Although the proposed algorithms have achieved satisfactory solutions, the convergence behaviours of these algorithms have yet to be proved rigorously. Since the objective function in the formulated problem (11) is non-convex, it is difficult to give the theoretical convergence proof of the proposed algorithms. Empirically, PSNR values achieved by these two proposed algorithms converge to the stable points as the iteration increases. Figs. 9 and 10 give the plots of PSNR values versus iterations at single observation and undersampled data cases. For observing clearly, the convergence curves of the six testing images are presented. As can be seen from the figures, with the growth of the iteration number the curves of the proposed PPR-BM3D algorithm increase monotonically and become flat and stable through approximately 10 iterations. Moreover, the convergence curves of the proposed PPR-FFDNet algorithm need 15 iterations to get stable. Therefore, the proposed two algorithms have good convergence property. 5.5. Running time The computational complexity is an important factor to evaluate the nonlinear imaging algorithm. We exploit the running time to measure the computational complexity. A long running time for recovering an image is not realistic. One may argue the filtered step of the proposed imaging framework plugged in a denoiser is time-consuming. However, our results show that the proposed algorithm plugged in FFDNet only needs nearly 10 seconds for recovering a 512 × 512 image. This imaging speed is faster than those of most benchmark algorithms. Table 4 shows the running time consumed by the seven testing imaging algorithms. We report the mean running time of the eight images at different sampling manners. As can be seen from the table, the imaging speed of the BM3D-prGAMP algorithm is the slowest. Many iterations of the BM3D denoise process should account for this result. The SPAR and DOLPHIn algorithms are faster than the BM3D-prGAMP algorithm, but slower than the ConPRbased (ConPR-BM3D and ConPR-FFDNet) algorithms and the PPRFFDNet algorithm. The dictionary learning step and the BM3D process are time-consuming in the DOLPHIn and SPAR algorithms, respectively. One can see that the proposed PPR-BM3D algorithm is faster than the other BM3D-based algorithms except for the ConPR-BM3D algorithm. Moreover, the PPR-FFDNet algorithm is faster than the benchmark algorithms except for the ConPR-FFDNet algorithm. The sparse coding step of the PPR-based algorithms should account for the slower imaging speed, compared with the ConPR-based algorithms. Nevertheless, the proposed algorithms can achieve better reconstructions than those of the ConPR-based algorithms. In fact, the denoiser has an important effect on the running time of the proposed framework. One can leverage the advanced denoiser based on deep learning for fast imaging. The running time of the sparse coding step can be further reduced by increasing the stride step of the convolution. Developing a faster imaging algorithm based on the proposed framework is our future direction.

95

6. Conclusion and future work In this paper, we presented a novel thought to design a plugand-play prior framework, and proposed a plug-and-play regularization model called PPR model. The proposed regularization model enforces the similarity of the underlying image and its filtered image in the sparsifying transform domain. Any denoiser can be utilized to filter the image, therefore, the regularization model is plug-and-play. The tight frame that is selected as the sparsifying transform not only can simplify the problem but also benefits for the image reconstruction. Prior information utilized in the denoiser and the sparsity over a tight frame can be utilized in the proposed PPR model. Through extensive testing we have demonstrated the effectiveness of the proposed PPR model. We exploited the proposed PPR model and the epigraph concept to solve nonlinear imaging inverse problems. A universal optimization model that is suitable for solving most nonlinear imaging inverse problems was formulated. The alternating optimization strategy was utilized to solve the formulated problem. To avoid the finely-tuned parameter, the epigraph concept was utilized to attack the image updating subproblem. Under the single noisy observation and the undersampled data scenarios, experiments showed that the proposed imaging algorithm could achieve better reconstructions, compared with the previous imaging algorithms. We open the door for designing a plug-and-play regularization model. By constructing the relation between the image and its filtered image, one can design a plug-and-play regularization model. One aspect is further researching on the relation of these two images. We proposed an effective imaging approach of exploiting the PPR model for solving nonlinear imaging inverse problems, and performed various experiments for coded diffraction imaging. Extension to solve other nonlinear imaging inverse problems such as scatting imaging problems is another open direction that is left for future research. Conflict of interest The authors declared that they have no conflicts of interest to this work. Acknowledgments This work was supported by the Scientific Research and Development Program of Qinhuangdao under grant 201805A003, in part by National Natural Science Foundation of China under Grant 61471313, and in part by the Doctoral Fund Project of Yanshan University under grant 8190019. The authors would like to thank Dr. A.M. Tillmann for his sharing matlab code of the DOLOHIn algorithm, and Dr. V. Katkovnik for his sharing matlab code of the SPAR algorithm. The authors would like to thank anonymous reviewers for helpful and stimulating comments. References [1] Y. Shechtman, Y.C. Eldar, O. Cohen, H.N. Chapman, J. Miao, M. Segev, Phase retrieval with application to optical imaging: a contemporary overview, IEEE Signal Process. Mag. 32 (3) (2015) 87–109. [2] Y.C. Yang, Q.S. Lian, S. Liu, B.S. Shi, X.Y. Fan, Resolution enhancement in digital in-line holography with sparsity, Opt. Eng. 57 (7) (2018) 073110. [3] S. Marchesini, H. He, H.N. Chapman, S.P. Hau-Riege, A. Noy, M.R. Howells, U. Weierstall, J.C. Spence, X-ray image reconstruction from a diffraction pattern alone, Phys. Rev. B 68 (14) (2003) 140101. [4] U.S. Kamilov, H. Mansour, B. Wohlberg, A plug-and-play priors approach for solving nonlinear imaging inverse problems, IEEE Signal Proc. Let. 24 (12) (2017) 1872–1876. [5] R.P. Millane, Phase retrieval in crystallography and optics, J. Opt. Soc. Amer. A 7 (3) (1990) 394–411. [6] R.W. Gerchberg, W.O. Saxton, A practical algorithm for the determination of phase from image and diffraction plane pictures, Optik 35 (1972) 237–246.

96

B. Shi, Q. Lian and X. Fan / Signal Processing 162 (2019) 83–96

[7] J.R. Fienup, Reconstruction of an object from the modulus of its fourier transform, Opt. Lett. 3 (1) (1978) 27–29. [8] S. Mukherjee, C.S. Seelamantula, Fienup algorithm with sparsity constraints: application to frequency-domain optical-coherence tomography, IEEE Trans. Signal Process. 62 (18) (2014) 4659–4672. [9] S. Loock, G. Plonka, Phase retrieval for fresnel measurements using a shearlet sparsity constraint, Inverse Probl. 30 (5) (2014) 055005. [10] E.J. Candes, T. Strohmer, V. Voroninski, Phase lift: exact and stable signal recovery from magnitude measurements via convex programming, Commun. Pur. Appl. Math. 66 (8) (2013) 1241–1274. [11] N. Vaswani, S. Nayer, Y.C. Eldar, Low rank phase retrieval, IEEE Trans. Signal Process. 65 (15) (2017) 4059–4074. [12] T. Blumensath, Compressed sensing with nonlinear observations and related nonlinear optimization problems, IEEE Trans. Inf. Theory 59 (6) (2013) 3466–3474. [13] H. Ohlsson, A.Y. Yang, R. Dong, S.S. Sastry, Nonlinear basic pursuit, in: Proceedings of 2013 Asilomar Conference on Signals, Systerms and Computers, Pacific Grove, CA, 2013, pp. 115–119. [14] E.J. Candes, X. Li, M. Soltanolkotabi, Phase retrieval via Wirtinger flow: theory and algorithms, IEEE Trans. Inf. Theory 61 (4) (2015) 1985–2007. [15] Y.X. Chen, E.J. Candes, Solving random quadratic systems of equations is nearly as easy as solving linear systems, Commun. Pur. Appl. Math. 70 (5) (2017) 0822–0883. [16] B.S. Shi, S.Z. Chen, Y. Tian, X.Y. Fan, Q.S. Lian, FASPR: a fast sparse phase retrieval algorithm via the epigraph concept, Digit. Signal Process. 80 (9) (2018) 12–26. [17] B.S. Shi, Q.S. Lian, S.Z. Chen, Y. Tian, X.Y. Fan, Coded diffraction imaging via double sparse regularization model, Digit. Signal Process. 79 (8) (2018) 23–33. [18] H.B. Chang, Y.C. Lou, K.N. Michael, Y.T. Zeng, Phase retrieval from incomplete magnitude information via total variation regularization, SIAM J. Sci. Comput. 38 (6) (2016) A3672–A3695. [19] A.M. Tillmann, Y.C. Eldar, J.M. Mairal, DOLPHIn-dictionary learning for phase retrieval, IEEE Trans. Signal Process. 64 (24) (2016) 6485–6500. [20] S.V. Venkatakrishnan, C.A. Bouman, B. Wohlberg, Plug-and-play priors for model based reconstruction, in: Proceedings of 2014 Global Conference on Signal and Information Processing (GlobalSIP), Austin, TX, 2013, pp. 945–948. [21] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein, Distributed optimization and statistical learning via the alternating direction method of multipliers, Found.Trends Machine Learn. 20 (3) (2011) 681–695. [22] S.H. Chan, X.R. Wang, O.A. Elgendy, Plug-and-play ADMM for image restoration: fixed-point convergence and applications, IEEE Trans. Comput. Imag. 3 (1) (2017) 84–98. [23] A. Beck, M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM J. Imag. Sci. 2 (1) (2009) 183–202. [24] C.A. Metzler, A. Maleki, R.G. Baraniuk, BM3d-PRGAMP: compressive phase retrieval based on BM3d denoising, in: Proceedings of International Conference on Image Processing (ICIP), Phoenix, AZ, USA, 2016, pp. 2504–2508. [25] F. Heide, M. Steinberger, Y.T. Tsai, et al., FlexISP: a flexible camera image processing framework, ACM Trans. Graphics (TOG) 33 (6) (2014). 231-1-13. [26] K. Dabov, A. Foi, V. Katkovnik, K. Egiazarian, Image denoising by sparse 3-d transform-domain collaborative filtering, IEEE Trans. Image Process. 16 (8) (2007) 2080–2095.

[27] V. Katkovnik, Phase retrieval from noisy data based on sparse approximation of object phase and amplitude, 2017, https://arxiv.org/abs/1709.01071. [28] V. Katkovnik, I. Shevkunov, N.V. Petrov, K. Egiazarian, Computational super-resolution phase retrieval from multiple phase-coded diffraction patterns: simulation study and experiments, Optica 4 (7) (2017) 786–794. [29] B.S. Shi, Q.S. Lian, X. Huang, N. An, Constrained phase retrieval: when alternating projection meets regularization, J. Opt. Soc. Amer. B 35 (6) (2018) 1271–1281. [30] A. Danielyan, V. Katkovnik, K. Egiazarian, Bm3d frames and variational image deblurring, IEEE Trans. Image Process. 21 (4) (2012) 1715–1728. [31] Y. Romano, M. Elad, P. Milanfar, The little engine that could: regularization by denoising (RED), SIAM J. Imag. Sci. 10 (4) (2017) 1804–1844. [32] K. Zhang, W.M. Zuo, L. Zhang, FFDNet: toward a fast and flexible solution for CNN based image denoising, IEEE Trans. Image Process. 27 (9) (2017) 4608–4622. [33] K. Zhang, W.M. Zuo, S. Gu, L. Zhang, Learning deep CNN denoiser prior for image restoration, in: Proceedings of IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017. Honolulu, HI, USA. [34] F. Heide, M. Rouf, M.B. Hullin, B. Labitzke, W. Heidrich, A. Kolb, High-quality computational imaging through simple lenses, ACM Trans. Graphics (TOG) 32 (5) (2013). 149-1-14 [35] F. Heide, S. Diamond, M. Niener, J. Ragan-Kelley, W. Heidrich, G. Wetzstein, Proximal: efficient image optimization using proximal algorithms, ACM Trans. Graphics (TOG) 35 (4) (2016). 84-1-15. [36] D.L. Donoho, I.M. Johnstone, Ideal spatial adaptation by wavelet shrinkage, Binometrika 81 (3) (1994) 425–455. [37] M. Tofighi, K. Kose, A.E. Cetin, Denoising using projections onto the epigraph set of convex cost functions, in: Proceedings of International Conference on Image Processing (ICIP), Paris, France, 2015, pp. 2709–2733. [38] O. WhyteEmail, J. Sivic, A. Zisserman, Deblurring shaken and partially saturated images, Int. J. Comput. Vision 110 (2) (2014) 185–201. [39] A. Mosleh, Y.E. Sola, F. Zargari, E. Onzon, J.M.P. Langlois, Explicit ringing removal in image deblurring, IEEE Trans. Image Process. 27 (2) (2018) 580–593. [40] J.F. Cai, H. Ji, Z. Shen, G.-B. Ye, Data-driven tight frame construction and image denoising, Appl. Comput. Harmon. A. 37 (1) (2014) 89–105. [41] L. Zhang, L. Zhang, X.Q. Mou, D. Zhang, FSIM: a feature similarity index for image quality assessment, IEEE Trans. Image Process. 20 (8) (2011) 2378–2386. Baoshun Shi received the Ph.D. degree in Yanshan University, in 2017. He is currently a lecturer in School of Information Science and Engineering, Yanshan University. His research interests cover diffraction imaging, phase retrieval, convolutional neural networks, dictionary learning. Qiusheng Lian received the Ph.D. degree in Yanshan University, in 2006. He is currently a professor of School of Information Science and Engineering at Yanshan University, China. His current research interests include phase retrieval, nonlinear compressed sensing, deep learning. Xiaoyu Fan is a Ph.D. candidate in Yanshan University. His current research interests include compressed sensing, magnetic resonance imaging.