It is possible for lensless imaging to achieve high image resolution in the imaging system [1], with a space–bandwidth product larger than that of lens-based setups [2]. Since no lens is required, this approach is cost effective and possible to be embedded in portable devices [3]. Lensless imaging has been applied to different disciplines, such as crystallography [4], material science [5], and biological imaging [6]. It mainly takes advantage of the diffraction process, which can encode the object information in captured data.

Among different configurations for lensless imaging, coherent diffractive imaging (CDI) is a powerful one with a simple imaging setup [7]. The object is illuminated with a coherent light source, and an imaging sensor is used to capture its far-field intensity image. Compared with conventional imaging setups, the imaging resolution of CDI is limited by the wavelength of the coherent light instead of the property of the lens. If coherent light with an ultrashort wavelength, such as X-rays, is used, the resolution can be greatly improved.

Although the optical setup for CDI is quite simple, the captured image is only the intensity pattern in the far field. Hence, an important step in the reconstruction algorithm, called phase retrieval, is required to retrieve the lost phase information from the interference fringes [8]. However, due to the non-uniqueness of the phase retrieval problem, some prior constraints in the object plane and imaging plane are required in most phase retrieval algorithms [9]. Some widely used constraints in the object plane are the size and location of the object [10].

When the structure of the object is complicated, these simple constraints are not sufficient to obtain satisfactorily reconstructed images. Conventional phase retrieval methods, such as the hybrid input–output (HIO) method [8], need many iterations and a good initial guess to solve this problem. There also exist some first-order algorithms for phase retrieval. Most are also iterative methods, such as the alternating direction method of multipliers (ADMM) [11–14], shrinkage or thresholding algorithms [15,16], and proximal algorithms [17,18]. However, some of these methods consider only the case in which latent images are real-valued, while in real applications, they are often complex-valued, and the phase information is also important. Also, most of them adopt different kinds of prior knowledge for regularization, such as a deep denoiser [15], ${\ell }_1$ or related regularization functions [16], but nearly all of them consider only natural images, which may not be efficient for the latent images in CDI. Deep learning methods can also be used for phase retrieval [19,20], but they need sufficient training images or accurate imaging models to generate synthetic training images.

The phase retrieval problem has also been investigated from the perspective of signal processing and optimization. The general problem they consider is the coded diffraction pattern (CDP) setup [21], which goes beyond the conventional CDI model with only the Fourier transform of the object. In this setup, different patterns are used as sensing matrices to spread the frequency information of a target, making it easier to reconstruct computationally with multiple measurements. A variety of phase retrieval algorithms are designed based on this, such as Wirtinger flow (WF) and its variants [22,23], Gauss–Newton method [24], spectral method [25], alternating minimization [26], and majorization–minimization [27]. These methods often have some restrictions on the sensing matrices or target objects. For example, sensing matrices are required to follow the Gaussian model [22,25,26], to be complete bases of the sensing space [24], or the target object is sparse [23]. Usually, more efforts are needed in designing optical systems to meet these requirements.

To take advantage of more accurate prior knowledge, it is possible to add some pre-designed strong constraints in the imaging systems directly. Some encoding masks, such as diffuser [28,29] and phase masks [30], can provide additional information to help with the reconstruction process. These coding elements are often hand-crafted with some prior knowledge, which shows the possibility of designing higher-resolution imaging systems for specific applications. For example, a simple masked CDI setup with an additional static area in the object plane, known as in situ CDI [31], is shown in Fig. 1. A dual pinhole is set in front of the object plane (middle plane in Fig. 1) to split the coherent light into two paths. One region on the object plane provides a static area in the imaging process, which can be regarded as efficient and solid prior knowledge in the object plane. The other part on the object plane is the sample area, which is dynamic. The imaging sensor is set in the far field to capture the oversampled Fourier intensity images. The image reconstruction problem for in situ CDI is to recover images of the dynamic area from the captured intensity images with the help of the static area. In these masked CDI setups, more prior knowledge can be used for algorithm design, which can improve the quality, robustness, and speed of reconstruction.

The main contributions of this paper are summarized as follows. •We propose a phase retrieval algorithm for masked CDI using ADMM. The proposed ADMM is in a plug-and-play (PnP) style, allowing for different hand-crafted or learning-based regularization functions with respect to complex-valued latent images.•Since the commonly used regularization functions in computational phase retrieval algorithms are designed or trained for real-valued natural images only, they are not suitable for complex-valued images in masked CDI experiments. We therefore design a regularization function based on the structure tensor and Harris corner detector in the complex domain.•We test our algorithm on real in situ CDI data, and the experimental results show that our method can achieve comparable results with the baseline method presented in Ref. [31], but the processing speed is $\sim 100\times$ faster per frame. Our method is also more robust even for the single diffraction pattern phase retrieval problem.

In a masked CDI setup, coherent light illuminates a coded aperture $p(x,y)$ and then transmits to the object plane. The incident light $l_i(x,y)$ can be calculated by the Fresnel diffraction integral [32] $$l_i(x,y)=\frac{e^{jk{z}_i}}{j\lambda z_i}{\iint }_{-\infty }^{+\infty }\exp \{\frac{jk}{2z_i}[{(x-\xi )}^2+{(y-\eta )}^2]\}\mathrm{d}\xi \mathrm{d}\eta ,$$(1)where $z_i$ is the distance between the mask and the object plane, $k$ is the wavenumber, and $\lambda$ is the wavelength. The transmissivity $\widehat{u}(x,y)$ of the whole object plane can be represented as $$\widehat{u}(x,y)=a(x,y)e^{j\theta (x,y)},$$(2)where $a(x,y)$ is the amplitude modulation, and $\theta (x,y)$ is the phase shift. In this work, we consider a more general mathematical model with respect to the whole object plane. We assume that in a masked CDI setup, the object plane is partially known, which is called the static region, and it remains static during the whole imaging process. Also, the location of the unknown object is fixed in a known region of interest, which is called a dynamic region, i.e., image of the object we are concerned with, and its content changes in different experiments. These two regions do not overlap with each other.

The whole object plane is composed of these two parts. A background mask $m_b(x,y)$ and an object mask $m_o(x,y)$ indicate the locations of the static region and dynamic region on the object plane, respectively, such that $$\widehat{u}(x,y)=m_b(x,y)r(x,y)+m_o(x,y)u(x,y),$$(3)where $r(x,y)$ is the transmissivity of the pinhole or other known patterns, and $u(x,y)$ is the transmissivity of the object of interest. Note that if the background mask $m_b$ is a zero matrix, then this model will revert to the conventional CDI model, where only the location and shape of the target object are known. In this algorithm, we incorporate the prior knowledge related to the static area into the algorithm design as a strong constraint, which can help with the reconstruction process.

After the object is illuminated by the incident light, the emergent light $l_e(x,y)$ can be obtained as $$l_e(x,y)=l_i(x,y)\widehat{u}(x,y).$$(4)Then, it propagates to the far-field detector plane such that $$l_o(x,y)=\mathcal{F}\{l_e(x,y)\},$$(5)where $\mathcal{F}\{·\}$ is the discrete Fourier transform. The captured image on the sensor is intensity only: $$o(x,y)=|l_o(x,y)|.$$(6)The task is to reconstruct $u(x,y)$ from $o(x,y)$ with the known masks $m_b,m_o$, static pattern $r$, and the probe $l_i$. This is the basic inverse imaging problem for the masked CDI setup, which is ill posed and non-convex.

For convenience, all the 2D matrices, such as $o(x,y)$, are vectorized and represented by their bold variables, such as $\mathit{o}$. The imaging model can be represented in a more compact form: $$\mathit{o}=|FL(W_b\mathit{r}+W_o\mathit{u})|,$$(7)where $L=\mathrm{diag}({\mathit{l}}_i),W_b=\mathrm{diag}({\mathit{m}}_b)$, $W_o=\mathrm{diag}({\mathit{m}}_o)$, and $F$ is the discrete Fourier transform matrix with a proper dimension. If the phase information $\mathit{\varphi }$ of the captured image $\mathit{o}$ is obtained, $\mathit{u}$ can easily be derived by the inverse step. We note that this model can also be regarded as a specific CDP setup, where the system matrix $A$ in the masked CDI setup is actually $A=FL{W}_o$. However, as summarized in Section 1, the general methods designed for the CDP setup have some requirements on the property of the system matrix $A$, such as the i.i.d. Gaussian assumption. In the masked CDI setup, the system matrix $A=FL{W}_o$ does not satisfy these requirements, making these methods less efficient for the masked CDI setup directly.

In real applications, the capturing process of many imaging devices, such as CCDs, is subject to various signal-dependent and signal-independent errors, which can be represented by a Poisson–Gaussian noise model [33]. This is typically handled by denoising [33] or introducing a strong regularization prior [34]. As the solution of a phase retrieval problem is not unique and the noise model should be taken into consideration, we introduce a regularization function $\mathcal{R}(·)$ to find a proper solution with respect to the prior knowledge [35]. Then, phase retrieval for masked CDI becomes an optimization problem: $$\mathrm{minimize}\mathcal{R}(\mathit{u})+\mathcal{C}(\mathit{\varphi })\mathrm{subject}\mathrm{to}FL\hat{\mathit{u}}=\mathit{o}\circ \mathit{\varphi }\hat{\mathit{u}}=W_b\mathit{r}+W_o\mathit{u},$$(8)where $\circ$ is the Hadamard product, and $\mathcal{C}(\mathit{\varphi })$ is an indicator function to make the modulus of the elements in $\mathit{\varphi }$ equal to one, i.e., $$\mathcal{C}(\varphi )=\{\begin{array}{ll}0,& |\mathit{\varphi }|=\mathrm{1}\\ \infty ,& \text{otherwise}\end{array}.$$(9)Note that this function is a circle in the complex plane, which is non-convex. For simplicity, we consider the setup in which the two masks are binary masks, and then $\mathit{u}$ can be directly obtained from $\widehat{\mathit{u}}$ as $$\mathit{u}=W_o(\widehat{\mathit{u}}-W_b\mathit{r}).$$(10)Since $m_o$ and $m_b$ do not overlap, $W_o{W}_b$ is a zero matrix. This means that when calculating $\mathit{u}$, we do not need the exact value of $\mathit{r}$, but $\widehat{\mathit{u}}$ is dependent on $\mathit{r}$.

The whole optimization problem we are trying to solve can be written as $$\begin{gathered} \mathrm{minimize}\mathcal{R}(\mathit{u})+\mathcal{C}(\mathit{\varphi }) \\ \mathrm{subject}\mathrm{to}FL\hat{\mathit{u}}=\mathit{o}\circ \mathit{\varphi }, \\ \mathit{u}=W_o(\hat{\mathit{u}}-W_b\mathit{r}). \end{gathered}$$(11)ADMM [36] can be used to solve this problem iteratively. First, using the Lagrange multiplier, the augmented Lagrangian is $$\mathcal{L}(\widehat{\mathit{u}},\mathit{\varphi },\mathit{u},\mathit{\lambda },\mathit{\mu })=\mathcal{R}(\mathit{u})+\mathcal{C}(\mathit{\varphi })+\frac{\rho }{2}{\Vert FL\widehat{\mathit{u}}-o\circ \mathit{\varphi }+\mathit{\lambda }\Vert }_2^2+\frac{\tau }{2}{\Vert \mathit{u}-W_o(\widehat{\mathit{u}}-W_b\mathit{r})+\mathit{\mu }\Vert }_2^2,$$(12)where $\mathit{\lambda }$ and $\mathit{\mu }$ are the Lagrange multipliers, and $\rho$ and $\tau$ are the penalty parameters. ADMM can minimize $\mathcal{L}$ by considering different variables separately. Only one variable is updated in each step. For example, in the $(k+1)$th iteration, when variable $\widehat{\mathit{u}}$ is updated, other variables $\mathit{\varphi },\mathit{u},\mathit{\lambda },\mathit{\mu }$ are regarded as constants with their current values. The $\widehat{\mathit{u}}$-update step aims at solving a sub-problem: $${\widehat{\mathit{u}}}^{(k+1)}=\underset{\widehat{\mathit{u}}}{\mathrm{argmin}}\mathcal{L}(\widehat{\mathit{u}},{\mathit{\varphi }}^{(k)},{\mathit{u}}^{(k)},{\mathit{\lambda }}^{(k)},{\mathit{\mu }}^{(k)}).$$(13)Then a necessary condition for ${\widehat{\mathit{u}}}^{(k+1)}$ is $${{\nabla }_{\widehat{\mathit{u}}}\mathcal{L}(\widehat{\mathit{u}},{\mathit{\varphi }}^{(k)},{\mathit{u}}^{(k)},{\lambda }^{(k)},{\mathit{\mu }}^{(k)})|}_{\widehat{\mathit{u}}={\widehat{\mathit{u}}}^{(k+1)}}=0.$$(14)In the complex domain, the required derivatives are regarded as Wirtinger derivatives [37]. The partial derivative of $\mathcal{L}$ with respect to variable $\widehat{\mathit{u}}$ is $${\nabla }_{\widehat{\mathit{u}}}\mathcal{L}=\rho {(FL)}^H(FL\widehat{\mathit{u}}-\mathit{o}\circ \mathit{\varphi }+\mathit{\lambda })+\tau W_o^H(W_o\widehat{\mathit{u}}-\mathit{u}-W_o{W}_b\mathit{r}-\mathit{\mu }).$$(15)After setting Eq. (15) to zero and substituting other variables with their current values, we can obtain the $\widehat{\mathit{u}}$-update step as $${\widehat{\mathit{u}}}^{(k+1)}={(\rho L^H{F}^{H}FL+\tau W_o^H{W}_o)}^{-1}[\rho {(FL)}^H(\mathit{o}\circ {\mathit{\varphi }}^{(k)}-{\mathit{\lambda }}^{(k)})+\tau W_o^H({\mathit{u}}^{(k)}+W_o{W}_b\mathit{r}+{\mathit{\mu }}^{(k)})].$$(16)

Similarly, other update steps are $${\mathit{\varphi }}^{(k+1)}={\mathcal{P}}_{\mathit{o}}({\mathit{\lambda }}^{(k)}+FL{\hat{\mathit{u}}}^{(k+1)}),$$(17)$${\mathit{u}}^{(k+1)}={\mathrm{prox}}_{\frac{\mathcal{R}}{\tau }}[W_o({\widehat{\mathit{u}}}^{(k+1)}-W_b\mathit{r}-{\mathit{\mu }}^{(k)})],$$(18)$${\mathit{\lambda }}^{(k+1)}={\mathit{\lambda }}^{(k)}+FL{\widehat{\mathit{u}}}^{(k+1)}-\mathit{o}\circ {\mathit{\varphi }}^{(k+1)},$$(19)$${\mathit{\mu }}^{(k+1)}={\mathit{\mu }}^{(k)}+{\mathit{u}}^{(k+1)}-W_o({\widehat{\mathit{u}}}^{(k+1)}-W_b\mathit{r}),$$(20)where ${\mathcal{P}}_o$ is a projection operator defined as $${\mathcal{P}}_o(\mathit{s}{)}_{(i)}=\{\begin{array}{cc}\frac{{\mathit{s}}_{(i)}/{\mathit{o}}_{(i)}}{\Vert {\mathit{s}}_{(i)}/{\mathit{o}}_{(i)}},& {\mathit{o}}_{(i)}\ne 0,{\mathit{s}}_{(i)}\ne 0\\ 1,& \mathrm{otherwise}\end{array}.$$(21)Variable $i$ refers to the $i$th element, and $\mathrm{prox}(·)$ is the proximal operator defined as [38] $${\mathrm{prox}}_{\frac{\mathcal{R}}{\tau }}(\mathit{v})=\underset{\mathit{x}}{\mathrm{argmin}}(\mathcal{R}(\mathit{x})+\frac{\tau }{2}\Vert \mathit{x}-\mathit{v}\Vert ).$$(22)The proposed ADMM algorithm for masked CDI is summarized in Algorithm 1. Algorithm 1.

PnP-ADMM for Masked CDI

We have developed the ADMM algorithm in a plug-and-play style, allowing for different kinds of regularization functions. Here, we investigate regularization functions that are efficient for denoising complex-valued images.

Regularization by denoising (RED) is a commonly used method of designing regularization functions according to the denoisers, instead of using denoisers directly to minimize some implicit functions [39]. In particular, the regularization function is designed as $$\mathcal{R}(\mathit{x})={\mathit{x}}^H[\mathit{x}-\mathcal{D}(\mathit{x})],$$(23)where $\mathcal{D}(·)$ is any proper denoiser. In the real domain, this function penalizes the residual $\mathit{x}-\mathcal{D}(\mathit{x})$ and the correlation with latent images at the same time, which is shown to be efficient [15,40].

Now, we consider the case $\mathit{x}={\mathit{x}}_{\mathit{r}}+j{\mathit{x}}_i$, where ${\mathit{x}}_{\mathit{r}}$ and ${\mathit{x}}_i$ are the real and imaginary parts, respectively. If $\mathcal{D}(·)$ is a learned denoiser, then for complex-valued images, the real and imaginary parts are regarded as two channels. The regularization function is $$\mathcal{R}(\mathit{x})=({\mathit{x}}_r+j{\mathit{x}}_i{)}^H\{[{\mathit{x}}_r-\mathcal{D}({\mathit{x}}_r)]+j[{\mathit{x}}_i-\mathcal{D}({\mathit{x}}_i)]\}.$$(24)This regularization function is not real-valued. Hence, it is not suitable to be used directly for regularization.

Even using trained denoisers cannot help to solve this problem. If the useful information is related only to phase information ${\mathit{x}}_p$ instead of amplitude information ${\mathit{x}}_a$, then the real and imaginary parts are ${\mathit{x}}_a\circ \sin {\mathit{x}}_p$ and ${\mathit{x}}_a\circ \cos {\mathit{x}}_p$, respectively. After trigonometric transformations, the prior knowledge for natural images is not suitable for the real and imaginary parts anymore. Hence, we need to design new regularization functions for complex-valued images directly based on some statistical properties.

In this work, we use the quadratic smoothing method [41], and the regularization function is designed as a quadratic penalty function: $$\mathcal{R}(\mathit{x})=\gamma {\mathit{x}}^{H}Q\mathit{x},$$(25)where $Q$ is a real-valued positive semi-definite matrix. For example, if $Q$ is selected as an identity matrix, then the regularization function is a simple ${\ell }_2$ norm, which minimizes the transmitted energy.

The structure tensor, also known as the second-moment matrix, is useful for image filtering, holography, feature tracking, etc. [42–44]. The original structure tensor is defined as $$S=G*(\nabla u\nabla u^H),$$(26)where $\nabla$ is the $2D$ differential operator, and $G$ is a Gaussian kernel for smoothing [45]. Generally, at each location $(x_0,y_0)$, this matrix has two eigenvalues $g_1$ and $g_2$ (assume $|g_1|\ge |g_2|$). When both $|g_1|$ and $|g_2|$ are small, there are only small variations in the neighborhood of this location, which means that this region is homogeneous. If $|g_1|$ is significantly larger than $|g_2|$, then there is only large variation in a dominate direction, which is the edge region. If $|g_1|$ and $|g_2|$ are both large, then this location has high variations in all directions, which means it is close to an image corner. The complexity of the structure for an image is highly related to the number of edges and corners inside. To avoid the explicit eigenvalue decomposition of the structure tensor, the Harris corner detector [46] adopts the trace and determinant to design the edge and corner responses.

Following the idea of a structure tensor and Harris corner detector, since we want to find the solution with the simplest structure in the solution set, we need to care about only the trace of the structure tensor and make it as small as possible to construct more homogeneous areas. Hence, we select $Q$ as a 2D Laplacian operator matrix, and the quadratic structure is then suitable for smoothing. The regularization function is now the sum of the traces of all structure tensors in each location, which can help to find a proper solution with the simplest structure.

There are also some other regularization functions that can construct homogeneous areas, such as total variation (TV) and its variants [47], and wavelet transform [48]. Compared with the quadratic regularization functions, the proximal operators of these functions have more complicated expressions. Most of them are designed based on the ${\ell }_1$ norm, whose proximal operator is a soft-thresholding function. We need to transfer these operators to the complex domain from the real domain, before embedding in the PnP-ADMM method.

To assess the efficiency of our method, we make use of the experimental in situ CDI data collected by Ref. [31]. The dataset consists of two optical experiments. The first experiment is an application in material science. In situ CDI is used to capture the growth of Pb dendrites on electrodes immersed in $\mathrm{Pb}{({\mathrm{NO}}_3)}_2$. A 543 nm He–Ne laser with a power of 5 mW is used to generate a collimated beam with a diameter of 800 μm. A dual pinhole aperture with two 100 μm circle holes spaced 100 μm apart is directly illuminated, and then the electro-chemical cell is 400 μm downstream of the aperture. A sealed fluid cell was assembled to observe the dynamics of Pb dendrites and separate it from the static area. An objective lens with 35 mm focal length is placed just downstream of the cell. One pinhole is placed in front of the cell, and another pinhole is focused on the region without any dendrite. A 12 V direct current source is applied to the cell. Finally, a CCD detector with $1340\times 1300\mathrm{pixels}$ resolution and $20\times 20\mathrm{\mu m}$ pixel size is used to capture images at the back focal plane of the lens. The location of the CCD detector is static, making the object-to-sensor distance the same for all diffraction patterns. A $5\times 5$ ptychography is required for each image. Another experiment is related to biological imaging. In situ CDI is adopted to monitor the activation of glioblastoma cells. The same He–Ne laser is used to illuminate the live glioblastoma cells sealed between two cover slips. The same ptychography is also required.

The parameters in ADMM are set as $\rho =\tau =1.0$ and $\gamma =0.01$, where the penalty parameters are set to be sufficiently large according to the convergence analysis in Ref. [49]. In the Pb growth experiment, the error tolerance is set as ${ϵ}_{\mathrm{tol}}=0.7$. In the glioblastoma experiment, the error tolerance is set as ${ϵ}_{\mathrm{tol}}=0.25$. The initial guess of the phase ${\mathit{\varphi }}^{(0)}$ is a vector with all elements equal to one, ${\widehat{\mathit{u}}}^{(0)}$ is initialized as the union of the two masks multiplied by the probe, with reference to the spectral initialization commonly used for the CDP setup [25,50], and all the other vectors are initialized as zero.

According to the imaging setup of in situ CDI, some prior information is provided in the dataset. The two binary masks $m_b$ and $m_o$ are given with reference to the structure of the dual pinhole. The incident light on the object plane, called the probe, is also known. The task is to reconstruct the complex-valued image of the object using prior knowledge and the diffraction patterns. We compare our algorithm with a baseline algorithm that is also available online [31]. Both algorithms are implemented in MATLAB on a computer with Intel Core i5-8500 CPU at 3.00 GHz and 16 GB RAM.

We first test our algorithm with the Pb growth experiment. In this case, the amplitude information shows the growth of Pb dendrites. The reconstruction results of ADMM and the baseline method are shown in Fig. 2. Different columns are imaging results with only intensity information at different times in the whole experiment. Imaging results of the Pb dendrites (blue part in the circle region) are almost identical using these two methods.

Then we introduce two metrics to evaluate the convergence property of ADMM. The R-factor evaluates the difference between reconstructed images and captured images in the measurement domain as $$R(\widehat{\mathit{u}},\mathit{o})={\Vert |FL\widehat{\mathit{u}}|-\mathit{o}\Vert }_2.$$(27)The relative error is defined as $${ϵ}_{\mathrm{rel}}^{(k+1)}=\frac{{\Vert {\mathit{u}}^{(k+1)}-{\mathit{u}}^{(k)}\Vert }_2}{\Vert {\mathit{u}}^{(k+1)}{\Vert }_2},$$(28)which shows the updates during the iterations. Note that for phase retrieval, we can recover the signal only up to a constant global phase shift ${\mathit{\varphi }}_0$ as ${\mathit{u}}^{*}{e}^{j{\varphi }_0}$, compared with the true value ${\mathit{u}}^{*}$. However, as we cannot obtain the value of ${\mathit{u}}^{*}$ or ${\mathit{u}}^{*}{e}^{j{\varphi }_0}$ in advance, we need to rely on the relative error ${ϵ}_{\mathrm{rel}}^{(k+1)}$ between two successive iterations to determine when to terminate the iteration process. We normally do so when it is smaller than a given tolerance, and the current solution is regarded as a stable solution. The R-factor and the relative error during ADMM iterations are shown in Fig. 3. We can see that both converge quickly. Since in real optical experiments, the optical conditions of the static area are not always the same as the known pattern (see Fig. 4), they do not converge to zero. Although these two planes are not the same, the main illumination patterns in the support area are similar.

To examine whether the proposed method can achieve comparable results with the baseline method, we calculate the relative peak signal-to-noise ratio (PSNR) and the structural similarity index (SSIM) [51] between the reconstruction results from ADMM and the baseline method, and the results are summarized in Table 1. All achieve $>31\mathrm{dB}$ PSNR and $>0.84$ SSIM. Visual results in Fig. 2 also show that the two methods can achieve comparable image quality in this experiment.Table 1.

Numerical Comparison between ADMM and Baseline Method for Pb Growth Experiment

We also compare our method with some other classical phase retrieval methods, mainly WF-based methods [22,52]. The phase retrieval results for the No. 1 image of the WF method and the truncated WF (TWF) method are shown in Fig. 5. We can find that even for the object with the simplest structure (No. 1 diffraction pattern) in this experiment, WF-based methods cannot provide valid information related to the Pb growth status. These results also show that the restrictions on the sensing matrices, such as the Gaussian model, are important to the phase retrieval qualities of these methods.

In the glioblastoma imaging experiment, the phase images show the status of glioblastoma cells. We first investigate the convergence of the ADMM with the relative error and R-factor in Fig. 6. As before, since we optimize only the dynamic region, the R-factor does not converge to zero, but a certain non-zero value after several iterations.

Imaging results with only phase information at different times are shown in Fig. 7. From Figs. 7(a) to 7(f) and Figs. 7(m) to 7(r), we can see that the smaller cell on the right moved to the bigger one, and then they merged together in these images. ADMM results show good agreement with the baseline in situ CDI results. However, there are more artifacts in the latter. As we can see in Figs. 7(k), 7(l), 7(s), 7(t), and 7(u), these artifacts may even change the shape of the merged cell, which can bring disturbances to the observation. With ADMM, the cell images are much clearer with less background noise and simpler image structures.

To give clearer evaluations of the reconstructed phase images of these two methods, we calculate the TV loss [53] of all the images, and the results are summarized in Table 2. We can see that for all diffraction patterns, ADMM has lower TV loss than the baseline method, which means that it can achieve phase retrieval results with simpler structures.Table 2.

Mean Total Variation Loss of Different Methods for Glioblastoma Experiment

The baseline method for in situ CDI is an updated version of the HIO method [8], and the update step is by gradient descent. Compared with this, a main advantage of our method is its high speed. The processing time of these two methods for different in situ CDI experiments is shown in Tables 3 and 4 and Fig. 8. ADMM needs only less than 0.5 s for a single diffraction pattern, while the baseline method needs more than 30 s for the same pattern. Specifically, the processing time of the baseline method is $39.29/0.37=106.19$ times longer than that of ADMM in the Pb growth experiment, and $40.38/0.34=118.76$ times longer in the glioblastoma experiment. Both the WF and TWF methods need more than 900 s for one diffraction pattern.Table 3.

Processing Time of Different Algorithms in Pb Growth Experiment (Unit: s)

In our method, we separate the whole optimization problem into several sub-problems that are easier to solve, where closed-form solutions are obtained for these sub-problems. Compared with the alternating methods between the image space and Fourier space, directly solving the inverse problem is more efficient. Although in one iteration ADMM uses the proximal method and is slower than the simple gradient descent step, it can converge with modest accuracy within a few iterations [36], while the baseline method needs many more iterations. Hence, our method is faster overall to give reconstructed images.

In the two experiments above, ADMM does not take advantage of the relationship among different diffraction patterns, but the baseline method combines these patterns to update the static area in the object plane, based on the fact that this area is the same among different patterns during the experiment. In this section, we compare our ADMM method with the baseline method for phase retrieval with only one diffraction pattern available to the algorithms.

We randomly choose one diffraction pattern from the Pb growth experiment and glioblastoma experiment. The phase retrieval results for different algorithms with different kinds of objects are shown in Fig. 9. We can easily find that the performance of the baseline method drops significantly for both experiments. The reconstructed image in Fig. 9(a) is less clear than in Fig. 9(b), the ADMM result. For a phase-only object, the baseline method even diverges and cannot give any meaningful information in Fig. 9(c), while ADMM in Fig. 9(d) is not influenced by the reduction of diffraction patterns. An explanation for this is that ADMM does not care about the region outside the dynamic area or the static area in the object plane, since they are mainly caused by noise, pinhole diffraction, etc. [see Fig. 4(b)]. Hence, although the imaging model is actually $$FL[W_o\mathit{u}+W_b\mathit{r}+(1-W_o-W_b)\mathit{n}]=\mathit{y}\circ \mathit{\varphi },$$(29)object $\mathit{u}$ is not related to disturbance $\mathit{n}$, and we can solve the phase retrieval problem directly with existence of these disturbances as noise. However, the baseline method tries to find the exact value of $\mathit{n}$ and then $\mathit{u}$, making the problem more complicated; hence, more diffraction patterns are required, and it does not perform well for the single diffraction pattern phase retrieval problem.

In summary, we present a fast and robust ADMM that can retrieve proper phase maps according to the diffraction patterns in masked CDI. Our method is plug-and-play, allowing for different kinds of regularization functions as prior knowledge. We also design quadratic regularization functions according to the structure tensor and Harris corner detector, which can denoise and find the simplest structure in the complex domain directly, instead of separating complex-valued images into real and imaginary parts, simplifying the algorithm design.

In optical experiments, we can see that compared with the conventional phase retrieval methods for masked CDI, our ADMM approach is much faster. Less background noise and clearer phase images further indicate the power of the designed regularization functions. Although the imaging results of Pb dendrites using the baseline method seem to have sharper edges with the help of the background patterns, in some more complicated cases, such as retrieving the phase of the phase-only objects in the glioblastoma experiment, the background patterns bring strong disturbances to the observation of the cells, making it less robust than ADMM. Our method also performs well even with a single diffraction pattern, showing the robustness of our method.