Phase contrast imaging is an essential label-free imaging technique because it enables the visualization of biological processes at multiple scales and resolutions in a non-invasive manner, providing a unique tool for cell division [1], intracellular dynamics [2], phenotypic screening [3], etc. It manipulates scattered radiation through optical field modulations to convert the phase information of non-self-luminous and non-absorbing samples into intensity signals. One of the most classical phase contrast imaging methods is based on interferometric imaging [usually quantitative, also known as quantitative phase imaging (QPI)] [4–7], whose highest spatial frequency, however, is confined to coherent diffraction limit $\frac{\mathrm{NA}}{\lambda }$, leading to a limited lateral resolution and poor environmental stability.

Partially coherent illumination endows phase contrast imaging with the superiority of high robustness and significantly extended lateral resolution, prompting its applications in biological imaging. The pioneering techniques in this field, such as Zernike phase contrast and differential interference contrast (DIC), have become standard solutions for cellular and subcellular observation [8–12]. Recent advances in partially coherent imaging have further boosted to the development of a series of QPI techniques, such as transport-of-intensity equation (TIE) [13–15], differential phase contrast (DPC) [16–18], and Fourier ptychographic microscopy (FPM) [19–24]. These techniques quantify the optical path length fluctuations across transparent samples and acquire three-dimensional (3D) quantitative phase, enabling phase contrast imaging to move from “qualitative” to “quantitative” [15,25–28]. Because of their superior imaging performance and flexible experimental configuration, QPI techniques have demonstrated their great potential in diverse biological applications [29–32].

QPI by DPC with partially coherent imaging is realized through asymmetric illumination modulation and deconvolution reconstruction, providing speckle-free imaging and twice the lateral resolution of the coherent diffraction limit [33,34]. It estimates the quantitative phase distribution of the sample from the measured intensities, which is considered to solve an inverse problem [35,36]. A strict and correct “intensity-phase” model is required to achieve accurate phase retrieval. However, mutual incoherence of the illumination points results in the nonlinear (or, more precisely, bilinear) dependence of the partially coherent field on the sample properties, which makes it a daunting challenge to recover the sample phase through intensity inversion. Thus, the primary task of QPI with partially coherent imaging is to accurately explain the nonlinear (bilinear) image generation model, a fundamental issue of great importance for the development of QPI. A well-established model is based on the transmission cross-coefficient (TCC) model, which establishes a four-dimensional (4D) expression on the pairs of spatial frequencies of the sample in the spatial frequency domain [37,38] so that the contribution of the sample- and system-dependent parts are represented separately. In addition, the phase-space representations, expressed through 4D Winger functions in the joint spatial frequency domain, model its bilinear characteristics with the advantage of separating the imaging system and the sample [39–42]. However, these high-dimensional representations are difficult to obtain precisely, while their massive data and computation costs make it impractical to recover sample information.

The preferred solution to achieve QPI under partially coherent imaging is to simplify the bilinear transmission of image information as linear properties. For this purpose, the sample (or object) approximations, including weak/slowly varying object approximation [17,33], Born/Rytov approximation [43,44], and multi-slice scattering approximation [45], have been proposed, all of which brought new insights into QPI. For example, weak object approximation considers the amplitude of the scattered light is much smaller than that of the unscattered light, linearizing the relationship between recorded intensity and sample transmission [46,47]. Then, a regularization deconvolution solver is proposed to recover the quantitative phase of the sample [17,48]. Slowly varying object approximation establishes a phase-gradient sensitive mechanism, leading to a “lookup table” correspondence between the measured intensity and the sample phase gradient [49,50]. However, the introduction of approximations means that the imaging performance of the algorithms derived from them is limited by how well the real sample matches the approximation model, making accurate phase (thickness) reconstruction of diverse samples from QPI alone difficult. Moreover, the explicit definitions of these approximations in existing work are still poorly understood. These limitations may result in QPI losing its proudest “quantitative” properties and no longer providing credible evidence for subsequent analyses such as cell morphology and cell mass measurement [51,52].

In this paper, we explore the strict definition of weak object approximation and further propose an iterative deconvolution DPC QPI technique that achieves accurate quantitative phase reconstruction without any additional acquisition, allowing QPI to be extended for biology applications of large-phase samples. The strict definition of weak object approximation is explicitly stated: the one-step deconvolution method could obtain accurate quantitative phase results only when the sample phase is not greater than 0.5 rad. Furthermore, to overcome the difficulty of applying DPC to large-phase samples, we introduce a pseudo-weak object approximation to model the complex transmittance function of the sample and propose an iterative deconvolution algorithm to achieve accurate QPI without additional data acquisition. The iterative process is performed with a deconvolution solver to refine the phase; thus, the phase components can be fully decoupled from the acquisition intensity. We conduct an experiment using a standard microlens array with a large phase of 37 rad (23 μm) and verify that the proposed method yields a quantitative phase consistent with the nominal value. Furthermore, the experiment results on MCF-7 cells accurately present the 3D morphological distribution of the cells, which demonstrates great potential for cell morphology and cell mass measurements.

QPI based on partially coherent imaging provides better lateral resolution and immunity to the system’s imperfections by simultaneously illuminating the sample from multiple angles in a mutually incoherent manner. The light intensity follows Abbe diffraction theory [53], i.e., mutual incoherence of the illumination points leads to the incoherent superposition of the intensity caused by each point. The cross propagation and superposition of each point often result in the nonlinear (or, more precisely, bilinear) dependence of the partially coherent field on the sample properties [37,38,54]: $$I(\mathbf{x})=\iiint S(\mathbf{u})T({\mathbf{u}}_1)T^{*}({\mathbf{u}}_2)H(\mathbf{u}+{\mathbf{u}}_1)H^{*}(\mathbf{u}+{\mathbf{u}}_2)e^{j2\pi \mathbf{x}·({\mathbf{u}}_1-{\mathbf{u}}_2)}\mathrm{d}{\mathbf{u}}_1\mathrm{d}{\mathbf{u}}_2\mathrm{d}\mathbf{u},$$(1)where $\mathbf{x}$ represents the two-dimensional (2D) coordinates of the camera plane, and $\mathbf{u}$ represents the spectrum coordinates. The sample is represented as a Fourier spectrum form of its complex transmittance function $T(\mathbf{u})=\mathcal{F}\{e^{j\varphi (\mathbf{x})}\}$ [where $\varphi (\mathbf{x})$ is the phase of the sample], and $T({\mathbf{u}}_{\mathbf{1}})T^{*}({\mathbf{u}}_{\mathbf{2}})$ is its mutual spectrum. $S(\mathbf{u})$ represents the partially coherent illumination source, which is usually an extended source whose illumination NA is less than or equal to the NA of the objective lens. $H(\mathbf{u})$ represents the coherent transfer function (CTF), which can be expressed as $H(\mathbf{u})=P(\mathbf{u})e^{jk\delta z\sqrt{1-{(\lambda u_x)}^2-{(\lambda u_y)}^2}}$ [$P(\mathbf{u})$ is the pupil function of the objective lens, and $e^{jk\delta z\sqrt{1-{(\lambda u_x)}^2-{(\lambda u_y)}^2}}$ is the defocusing factor for propagating the light field on different planes]. From Eq. (1), the intensity $I(\mathbf{x})$ of each point in the image is considered as a point pair dependent on the sample amplitude. To separate the effects of the system and the sample properties, the sample-independent factors are established as TCC (or the partially CTF in Sheppard’s paper Ref. [33]), which is a 4D function of pairs of spatial frequencies: $$\mathrm{TCC}({\mathbf{u}}_1,{\mathbf{u}}_2)=\int S(\mathbf{u})H(\mathbf{u}+{\mathbf{u}}_1)H^{*}(\mathbf{u}+{\mathbf{u}}_2)\mathrm{d}\mathbf{u}.$$(2)

This equation illustrates that TCC is determined by the illumination function $S(\mathbf{u})$ and CTF, and it establishes a quantitative relationship between sample distribution and acquired intensity. It can be found from Eq. (1) that the contributions from sample amplitude and phase in the intensity signal under partially coherent imaging cannot be expressed explicitly. This makes it a daunting challenge to recover the sample phase through intensity inversion.

Object approximation is the most common method for simplifying the partially coherent imaging model to a mathematically analytic form [34]. Weak object approximation achieves the linearization of the intensity expression, establishing an effective phase inversion mechanism (the forward model of partially coherent imaging and QPI algorithm are given in Appendices A and B). It requires the phase terms of the sample to be very weak (usually much less than 1 rad), and then the complex transmittance function of the sample $t(\mathbf{x})=e^{j\varphi (\mathbf{x})}$ can be approximated as $$t(\mathbf{x})=e^{j\varphi (\mathbf{x})}\approx 1+j\varphi (\mathbf{x}).$$(3)

This expression ignores a series of higher-order terms in phase and retains only its primary term, successfully separating the DC (background) and phase terms into the real and imaginary parts of the complex distribution. Then, a linear expression of mutual intensity is substituted into Eq. (1) to obtain the intensity distribution under the weak phase approximation [55], $$I(\mathbf{x})\approx \mathrm{TCC}(\mathbf{0},\mathbf{0})+2\mathrm{Re}\{\int \mathrm{TCC}(\mathbf{u},\mathbf{0})[j\mathrm{\Phi }(\mathbf{u})]e^{j2\pi \mathbf{x}·\mathbf{u}}\mathrm{d}\mathbf{u}\},$$(4)where $\mathrm{TCC}(\mathbf{0},\mathbf{0})$ is the direct current (DC) component (the background). $\mathrm{TCC}(\mathbf{u},\mathbf{0})$ is the weak phase transfer function (WPTF), characterizing the quantitative relationship between the sample phase and the acquisition intensity. For QPI, which converts $j\mathrm{\Phi }(\mathbf{u})$ to the intensity with real distribution modulated by WPTF, the illumination $S(\mathbf{u})$ or CTF $H(\mathbf{u})$ modulations are required to produce a phase transfer response with an oddly symmetric distribution.

DPC adopts asymmetric illumination, a flexible and convenient implementation without mechanical manipulation, to achieve an oddly symmetric WPTF and evenly symmetric absorption transfer function [16,17]. Consequently, a simple differential calculation can be implemented to eliminate the background item from Eq. (4), and only linear phase items are left in $I_i^{\mathrm{DPC}}(\mathbf{x})$, $$I_i^{\mathrm{DPC}}(\mathbf{x})=\frac{I_{i1}(\mathbf{x})-I_{i2}(\mathbf{x})}{I_{i1}(\mathbf{x})+I_{i2}(\mathbf{x})},$$(5)where $i$ indicates the direction of asymmetric illumination $S_i({\mathbf{u}}^{\prime })$, and $I_{i1}(\mathbf{x})$ and $I_{i2}(\mathbf{x})$ denote the two intensity images along the same asymmetric axis. Then, the corresponding WPTF can be expressed as $$\mathrm{WPTF}(\mathbf{u})=\frac{\int S_i({\mathbf{u}}^{\prime })H^{*}({\mathbf{u}}^{\prime }+\mathbf{u})H({\mathbf{u}}^{\prime })-S_i({\mathbf{u}}^{\prime })H({\mathbf{u}}^{\prime }-\mathbf{u})H^{*}({\mathbf{u}}^{\prime })\mathrm{d}{\mathbf{u}}^{\prime }}{\int |S_i({\mathbf{u}}^{\prime })||H({\mathbf{u}}^{\prime }{)|}^2\mathrm{d}{\mathbf{u}}^{\prime }}.$$(6)

This leads to a linear relationship between the phase of the sample and the acquired intensity, and the quantitative phase information can then be reconstructed by Fourier space deconvolution with the WPTF [16–18,56]. Tikhonov regularization is commonly used to solve the ill-posed problem to obtain a stable approximate solution, $$\varphi (\mathbf{x})={\mathcal{F}}^{-1}\left\{\frac{\sum _i[{\mathrm{WPTF}}_i^{*}(\mathbf{u})·I_i^{\mathrm{DPC}}(\mathbf{u})]}{\sum _i|{\mathrm{WPTF}}_i(\mathbf{u}{)|}^2+\beta }\right\},$$(7)where $\beta$ represents the regularization parameter to avoid excessive amplification of noise.

Weak object approximation simplifies partially coherent imaging; however, it also limits the applicability of the phase inversion algorithm. Specifically, only if the weak object condition is met can the quantitative phase distribution of the sample be recovered accurately without distortions and artifacts. Otherwise, DPC will suffer from severe nonlinear errors and lose its nature of “quantitative.” The analysis of the theoretical model of the deconvolution solution reveals that two approximations are involved in its derivation. On the one hand, weak object approximation linearizes the complex transmittance function of the sample by omitting the higher-order terms in its Taylor expanded form, as shown in Eq. (3). On the other hand, the higher-order scattering terms are ignored to derive the analytical expression of the WPTF, as shown in Eq. (4). Although the two approximations jointly contribute to the deconvolution reconstruction algorithm, the intensity difference of Eq. (5) eliminates the nonlinear error caused by the latter approximation. Thus, the reconstruction error in the DPC deconvolution mainly arises from the mismatch between the ideal sample function $e^{j\varphi (\mathbf{x})}$ and the linear sample function $1+j\varphi (\mathbf{x})$. In Appendix C, we prove this conclusion by simulating the one-step deconvolution reconstruction under both sample models, and the results are shown in Fig. 7. In addition, a critical corollary suggests that although weak object approximation is only described as a sample simplification model (which appears to be relevant only to the sample’s absolute phase), the acquired intensity under angularly varying illumination exhibits different nonlinear characteristics. Then, weak object approximation will result in different degrees of information being ignored. This implies that the mismatch between the actual and linear intensity model under weak object approximation is related to illumination aperture.

Figure 1 explains the information mismatch between the actual and linear intensity model under weak object approximation with different illumination by numerical simulation. The complex transmittance functions of $e^{j\varphi (\mathbf{x})}$ and $1+j\varphi (\mathbf{x})$, which express the actual sample distribution and the linear distribution, were used to generate the DPC acquisition images at different illumination coefficients $s_{\mathrm{ill}}$ (where $s_{\mathrm{ill}}=\frac{{\mathrm{NA}}_{\mathrm{ill}}}{{\mathrm{NA}}_{\mathrm{obj}}}$, indicating the ratio of illumination ${\mathrm{NA}}_{\mathrm{ill}}$ and objective ${\mathrm{NA}}_{\mathrm{obj}}$). The acquisition images were defined on a grid with $400\times 400\mathrm{pixels}$ with a pixel size of $0.3\mathrm{\mu m}\times 0.3\mathrm{\mu m}$. The wavelength of the illumination is 525 nm, and the ${\mathrm{NA}}_{\mathrm{obj}}$ is 0.4. All subsequent simulations follow these parameter configurations. The corresponding intensity differences and root mean square error (RMSE) values were then calculated to measure the model mismatch, as shown in Figs. 1(d) and 1(e). It can be seen that the intensity difference between the two sample models grows exponentially with the increasing illumination coherence coefficients ($s_{\mathrm{ill}}$ from 0 to 1). This will result in inaccurate phase reconstruction by the deconvolution solver under weak object approximation. We further display the reconstructed phase at $s_{\mathrm{ill}}=1$ in Fig. 1(e), which shows a significantly weaker phase value compared to the ground-truth phase. Figure 1(f) shows the 3D morphological distribution of the phase loss of the reconstructed phase compared to the ground-truth phase. In Appendix C (Figs. 7 and 8), two simulations further illustrate that the imaging performance of the one-step deconvolution solver is related to the illumination aperture. Therefore, a conclusion can be drawn that the definition of weak object approximation should not only consider the sample but should be analyzed in terms of the QPI performance to take into account the illumination aperture and sample.

The well-defined weak object approximation is the basis for guaranteeing the accuracy of quantitative phase reconstruction, while it also allows the imaging performance of QPI to be measured and analyzed in an informed manner. However, due to the complex theoretical derivation of the deconvolution solver, it is challenging to determine the strict definition of weak object approximation from the mathematical expression. To address this problem, we analyzed the imaging performance of DPC QPI at maximum nonlinear error to explore the strict weak object approximation definition applicable to all samples and illumination apertures. Two representative samples of a slowly varying microlens array and three sharply varying steps were selected as target samples, and their phase amplitudes can be adjusted arbitrarily. Meanwhile, the illumination aperture can be set from half-circle to half-annulus to simulate different illuminations. The largest nonlinear error will arise when the sample of three sharply varying steps is illuminated by a half-annular aperture with the maximum illumination coefficient. In this case, an accurate reconstruction of the sample phase using the one-step deconvolution will serve as the basis for the strict definition of weak object approximation.

Figure 2 shows the numerical simulation results for two representative samples with varying phases from 0.2 rad to 2 rad under different illumination apertures. As shown in Fig. 2(a), three types of illumination apertures were set, including a half-circular aperture, a half-annular aperture of 0.5 width, and a half-annular illumination of 0.01 width with the maximum illumination coefficient. The half-circular aperture and the half-annular aperture of 0.01 width correspond to two extremes of the nonlinear error. Specifically, the intensity superposition under a half-circular illumination weakens the nonlinear error; however, this error will exist in the intensity with the maximum weight under the half-annular illumination of 0.01 width. To quantitatively measure the reconstruction errors, the RMSE values between these reconstructed phases and the ground-truth phases are calculated to plot the numerical curves. As shown in Figs. 2(b) and 2(c), the reconstruction error increases monotonically with increasing phase amplitude. This indicates that the difference between the acquired intensity and the linear intensity desired by the one-step deconvolution increases when imaging samples with larger phases. The reconstructed phase of the sharply varying step sample under the half-annular illumination of 0.01 width results in a more serious reconstruction error, which is consistent with our inference. To measure phase reconstruction performance, an RMSE value of less than 1% was used to determine the accurate reconstructed phase since the reconstruction error can be almost neglected in this case. From Figs. 2(b) and 2(c), when the phase is less than or equal to 0.5 rad, the reconstruction RMSE at any sample and illumination aperture is less than 1%. Therefore, the strict weak object approximation is suggested to be defined as not greater than 0.5 rad [$\varphi (\mathbf{x})\le 0.5\mathrm{rad}$]. Figures 2(d) and 2(e) further give the reconstructed phases of these two types of samples under the half-annular illumination of 0.01 width. When the phase amplitude is equal to 0.5 rad, the DPC reconstruction yields a reconstructed phase that is consistent with the ground-truth phase, which can no longer be achieved when the phase amplitudes are 1 rad and 1.5 rad.

To achieve accurate QPI of large-phase objects, we introduce pseudo-weak object approximation in electron microscopy into DPC and develop a more precise object approximation model [57,58]. Considering a sample distribution $t(\mathbf{x})=e^{j\varphi (\mathbf{x})}$ with a large phase, it can be modeled as a superposition of $n$ layers satisfying weak object approximation $\varphi (\mathbf{x})={\varphi }_1(\mathbf{x})+{\varphi }_2(\mathbf{x})+\cdots +{\varphi }_n(\mathbf{x})$ [where each layer has equal phase ${\varphi }_1(\mathbf{x})={\varphi }_2(\mathbf{x})=\cdots ={\varphi }_n(\mathbf{x})\le 0.5\mathrm{rad}$]. Then the sample distribution can be approximated as $$\begin{gathered} t(\mathbf{x})=e^{j[{\varphi }_1(\mathbf{x})+{\varphi }_2(\mathbf{x})+\cdots +{\varphi }_n(\mathbf{x})]} \\ \approx [1+j{\varphi }_1(\mathbf{x})][1+j{\varphi }_2(\mathbf{x})]\cdots [1+j{\varphi }_n(\mathbf{x})] \\ \approx 1+j\sum _r^n{\varphi }_r(\mathbf{x})-\sum _{r,s}^{c_n^2}{\varphi }_r(\mathbf{x}){\varphi }_s(\mathbf{x})-j\sum _{r,s,t}^{c_n^3}{\varphi }_r(\mathbf{x}){\varphi }_s(\mathbf{x}){\varphi }_t(\mathbf{x})+\dots , \end{gathered}$$(8)where $r,s,t,\dots$ indicates the layer index of the sample ($\mathrm{1,2,}\dots ,n$), and any two of them are required to be unequal, which means $r\ne s$, $r\ne t$, or $s\ne t,\dots$. This equation is the expression of pseudo-weak object approximation, which considers an object with a large phase beyond the limit of weak object approximation as a multi-layer phase structure so that its complex transmittance function can be precisely modeled as a superposition of the phases and their cross-mixed terms. The previous weak object approximation $1+j{\varphi }_1(\mathbf{x})$ can then be considered as a special case where the sample has only a single layer.

Considering the intensity distribution under the DPC difference model, the real part of the intensity spectrum distribution derived from the convolution of the complex transmittance function will be eliminated under asymmetric illumination, leaving only its imaginary part terms. When the sample phase satisfies the strict weak object approximation ($\le 0.5\mathrm{rad}$), the first-order phase term in the imaginary part will be the main contribution of the DPC intensity, as shown in Eq. (4). Then, the differential intensity and the sample phase are linearly related by the WPTF, which derives the one-step deconvolution algorithm in Eq. (7) to achieve accurate phase reconstruction. Once the sample phase exceeds the strict weak phase approximation condition, the phase and the differential intensity cannot be interpreted as linearly correlated, so the one-step deconvolution is no longer applicable. In this case, the pseudo-weak object approximation can be introduced to treat the sample with a more accurate model. From Eq. (8), the sample is approximated as a superposition of a series of weak object components. This indicates that the acquired image can be expressed as the superposition of linear intensity expected by deconvolution and nonlinear intensity errors. Thus, solving the phase can then be considered an optimization problem, which can be efficiently handled by an iterative reconstruction algorithm. As a result, the phase of the sample beyond weak object approximation can be accurately recovered.

The primary motivation of the iterative reconstruction algorithm is to eliminate the nonlinear components and extract the intensity components that are linearly related to the phase of the large-phase object. Based on pseudo-weak object approximation, the nonlinear intensity error can be obtained by calculating the intensity residual between the forward-generated intensity of partially coherent imaging and the phase-linear correlated intensity. Consequently, the iterative algorithm requires no additional acquisition data (usually four images) and refines the reconstructed phase by a continuous deconvolution solver. In this process, the input intensity used for the deconvolution process is denoted as $$I_{\text{linear}}^{\mathrm{DPC}}(\mathbf{x})=I_{\mathrm{real}}^{\mathrm{DPC}}(\mathbf{x})-\mathrm{\Delta }{I}_{\mathrm{error}}^{\mathrm{DPC}}(\mathbf{x}),$$(9)where $I_{\text{linear}}^{\mathrm{DPC}}(\mathbf{x})$ is the estimated linear intensity for deconvolution that eliminates the nonlinear intensity error $\mathrm{\Delta }{I}_{\mathrm{error}}^{\mathrm{DPC}}(\mathbf{x})$, and $I_{\mathrm{real}}^{\mathrm{DPC}}(\mathbf{x})$ is the phase gradient image calculated using the actual acquired image according to Eq. (5). The nonlinear error $\mathrm{\Delta }{I}_{\mathrm{error}}^{\mathrm{DPC}}(\mathbf{x})$ can be approximately estimated by solving the intensity difference between the intensity model based on Abbe’s theory and the transfer function model. On the one hand, an initial phase is used to generate the acquired intensity of DPC by superimposing the intensity corresponding to each point source in the asymmetric illumination, $$I_{\mathrm{real}}(\mathbf{x})=\int |{\mathcal{F}}^{-1}\{S({\mathbf{u}}^{\prime })\mathcal{F}[e^{j\varphi (\mathbf{x})}]P(\mathbf{u}-{\mathbf{u}}^{\prime }{)\}|}^2\mathrm{d}{\mathbf{u}}^{\prime }.$$(10)

The generated intensity images under the same asymmetric axis are used to calculate the phase gradient image $I_{\mathrm{real}}^{{\mathrm{DPC}}^{\prime }}(\mathbf{x})$ by Eq. (5). On the other hand, the linearized intensity required by the deconvolution solver of DPC is calculated according to the transfer function theory, $$I_{\text{linear}}^{{\mathrm{DPC}}^{\prime }}(\mathbf{x})={\mathcal{F}}^{-1}\{\mathrm{WPTF}(\mathbf{u})\mathrm{\Phi }(\mathbf{u})\}.$$(11)

Then, the nonlinear error caused by the mismatch between these two models can be estimated as $$\mathrm{\Delta }{I}_{\mathrm{error}}^{\mathrm{DPC}}(\mathbf{x})=I_{\mathrm{real}}^{{\mathrm{DPC}}^{\prime }}(\mathbf{x})-I_{\text{linear}}^{{\mathrm{DPC}}^{\prime }}(\mathbf{x}).$$(12)

During the iterative process, $\mathrm{\Delta }{I}_{\mathrm{error}}^{\mathrm{DPC}}(\mathbf{x})$ will be continuously estimated and eliminated from the actual acquisition intensity, which leads to $I_{\text{linear}}^{\mathrm{DPC}}(\mathbf{x})$ close to the intensity linearly related to the true phase of the sample. In order to solve the phase, the algorithm can be derived as the following optimization problem: $$\min \epsilon =\min \sum _{\mathbf{x}}\{I_{\mathrm{real}}^{\mathrm{DPC}}(\mathbf{x})-I_{\mathrm{real}}^{{\mathrm{DPC}}^{\prime }}(\mathbf{x})\}.$$(13)

Then, the phase of the sample is iteratively updated based on the following equation: $${\varphi }_{m+1}(\mathbf{x})={\mathcal{F}}^{-1}\left\{\frac{\sum _i\mathcal{F}[I_{\text{linear},i,\mathrm{m}}^{\mathrm{DPC}}(\mathbf{x})]·{\mathrm{WPTF}}_i^{*}(\mathbf{u})}{\sum _i|{\mathrm{WPTF}}_i(\mathbf{u}{)|}^2+\beta }\right\},$$(14)where $m$ is the number of iterations. $\beta$ represents the regularization parameter to avoid excessive amplification of noise.

The algorithm flow chart of the iterative deconvolution is displayed in Fig. 3. The algorithm begins with the one-step deconvolution reconstruction, where the solved phase is used as the initial phase. In the iterative process, the solved phase and a uniform amplitude construct a complex transmittance function to generate DPC images under asymmetric illuminations based on Abbe’s method. Specifically, the acquired image under partially coherent imaging is represented as a linear superposition of the intensities obtained by multiple point sources independently illuminating objects with the complex distribution of $e^{j\varphi (\mathbf{x})}$. They were used to solve the phase gradient intensity based on Eq. (5). Synchronously, the estimated phase is used to perform a convolution with the WPTF to obtain an intensity distribution linearly related to the phase. The intensity images under both models are used to estimate the nonlinear intensity error caused by the mismatch of sample models based on Eq. (12). Then, this nonlinear error will be subtracted from the actual acquisition intensity, whereby the intensity signals that interfere with the deconvolution algorithm are eliminated, leaving only the linear intensity desired by the deconvolution [Eq. (9)]. Finally, as expressed in Eq. (14), a deconvolution reconstruction is then performed using estimated linear intensity to recover the updated phase of the current iteration. The above iterative process will be performed several times until the difference between the estimated intensity generated by the reconstructed phase and the acquired image can be negligible. As a result, the iteration converges to a stable phase value.

We demonstrated the imaging performance of iterative deconvolution through numerical simulations. The same microlens distribution similar to the previous simulation is used to generate the DPC images. In order to explore the phase imaging performance of the iterative deconvolution for large-phase objects, we constructed a sample with a phase distribution of $0\text{−}10\mathrm{rad}$. This phase distribution corresponds to a sample with a thickness of more than 1 μm, which goes far beyond the definition of weak object approximation. The illumination aperture is first set as a half-circular illumination, whose corresponding WPTF is shown in Fig. 4(a). The one-step deconvolution and iterative deconvolution (performed 10 times) algorithms were used to solve the phase of the sample, respectively, and the resulting reconstructed phases and their 3D pseudo-color rendering images are shown in Fig. 4(a). Figure 4(a) also plots the numerical profiles of this reconstructed phase and the ground-truth phase for a quantitative comparison. It can be seen that the reconstructed phase under one-step deconvolution demonstrates a phase amplitude of 5 rad, which is much less than the ground-truth phase amplitude of 10 rad. Such a reconstruction result gives a misleading measurement. The iterative deconvolution yields the accurate reconstructed phase, which shows a consistent phase amplitude with the ground-truth phase. We further performed a simulation under the annular illumination aperture with an annulus width of 0.25, and the simulation results are shown in Fig. 4(b). Although the nonlinear error in the acquisition intensity increases, the iterative deconvolution accurately recovers the quantitative phase of the samples.

To verify the phase reconstruction accuracy of iterative deconvolution, several experiments were conducted on different samples. The experimental system was built on a commercial inverted microscope platform (Olympus, IX73), and its illumination source was replaced by a programmable LED array with a spacing of 2 μm to flexibly implement asymmetric illumination at a wavelength of 504 nm. In image acquisition process, four illumination patterns (left, right, up, and down) were switched, and the transmitted light field information of the sample was collected by an objective lens with $10\times /0.4\text{NA}$ (Olympus UPlanSApo). In the camera plane, a CCD camera sensor (The Imaging Source, DMK 33UX183) with a pixel size of 2.4 μm was used to record intensity images. As a result, raw images can be completed within 200 ms.

To quantitatively demonstrate the accuracy of the iterative deconvolution method for reconstructing the phase, a standard microlens array sample [SUSS, Microoptics, refractive index (RI) of 1.46] with a curvature radius of 350 μm and a pitch of 250 μm was used as the target. It was placed in water with a dielectric RI of 1.33, so its theoretical height of 23 μm corresponds to a phase amplitude of 37 rad. Figures 5(a1), 5(a2) and 5(b1), 5(b2) display the asymmetric illumination patterns and the phase gradient images in two directions along the corresponding shearing directions. These captured images were used to perform the one-step deconvolution reconstruction directly, and the reconstructed phase is shown in Fig. 5(d1). A significantly weaker phase contrast can be found, indicating that the reconstructed phase in the traditional one-step deconvolution method is smaller than the theoretical phase value of the microlens array. The acquired images were then used to implement the proposed iterative deconvolution method (which takes several minutes). Figure 5(c) shows the variation curve of the intensity differences in iterations [Eq. (13)], which gradually converges to a stable small value as the iteration number increases. Once there is a constant minimum error between the estimated and acquired intensities, the iteration will be terminated, and the final iteration result will be output, as shown in Fig. 5(d2). We further selected one microlens unit as the region of interest (ROI) for the quantitative characterization of its 3D morphology, and the results are shown in Figs. 5(e1), 5(e2) and 5(f1), 5(f2), respectively. From the 3D pseudo-color morphological distribution in Figs. 5(f1) and 5(f2), the iterative deconvolution correctly characterizes the true 3D morphology of the microlens, which is difficult to achieve by one-step deconvolution. Finally, in order to quantitatively compare the numerical accuracy of the reconstructed phases under these two algorithms, we extracted the quantitative phase values along the white dashed and plotted the curve in Figs. 5(g1) and 5(g2). It verifies that the iterative deconvolution algorithm provides reliable quantitative phase data for subsequent research and analysis.

Morphological detection of biological cells is critical for the early confirmation of cancer treatment, where the cell changes during malignant transformation can be quantified by studying morphology, membrane dynamics, and cell refraction indices. DPC QPI provides a reliable pathology testing means for cancer tissues and cells. We conducted an experiment on MCF-7 human breast cancer cells to compare the cell morphology testing performance of one-step deconvolution and iterative deconvolution algorithms. During sample preparation, the MCF-7 cells were placed in 4% paraformaldehyde solution without any staining and labeling treatment. Figure 6 displays the reconstructed results under one-step deconvolution and iterative deconvolution. We adopted the Tikhonov regularization and introduced a suitable regularization parameter for both deconvolution processes. In Fig. 6(a), we show the field-of-view (FOV) reconstructed phase under iterative deconvolution. Then, the phase values along the white dashed line were extracted to plot the quantitative curves, yielding the results in Figs. 6(b) and 6(c), respectively. It can be clearly observed that for mitotic cells with large phase values, one-step deconvolution (blue curve) provides a degraded cell morphology feature, while the iterative deconvolution (red curve) accurately recovers both low-frequency and high-frequency quantitative features of the sample, demonstrating accurate 3D morphological data. Figures 6(d1)–6(g1) and 6(d2)–6(g2) present the enlarged phase of four ROIs under one-step deconvolution and iterative deconvolution. Compared to the one-step deconvolution algorithm, the iterative deconvolution algorithm recovers the cellular contours and the internal subcellular structures with highly accurate quantitative phase distribution. A 3D pseudo-color rendering image and a quantitative profile are then plotted to characterize the 3D morphological distribution of MCF-7 cells, as shown in Fig. 6(h). It can be found that the iterative deconvolution demonstrates a more significant cell thickness, which provides more accurate morphological detection data for cell analysis.

In this paper, we have explored the strict definition of weak object approximation under partially coherent imaging and proposed an iterative deconvolution algorithm to achieve accurate QPI for large-phase objects. By exploring the imaging performance of QPI under weak object approximation with different illumination coherence and sample distribution, the strict weak object approximation is defined as a sample with a phase no greater than 0.5 rad. Thus, QPI by DPC enables accurate quantitative phase reconstruction for arbitrarily samples and illumination apertures. Furthermore, an iterative deconvolution QPI algorithm based on pseudo-weak object approximation was proposed to go beyond weak object approximation, achieving accurate QPI for large-phase objects without any additional acquisition. As a result, QPI by DPC is no longer limited by the premise of weak object approximation and can be applied to the high-precision phase characterization and measurement of large-phase samples. We demonstrated the accuracy of the proposed method using a standard microlens array sample with a large phase. Compared with the conventional one-step deconvolution, the proposed method more accurately recovers a phase consistent with the nominal value. The experiment on MCF-7 cells demonstrated that our method accurately reflects the 3D morphological data of the cells, which will provide strong support for subsequent cell morphology and cell mass analysis. Our method will enable the quantitative nature of DPC QPI and provide a powerful tool for the quantitative characterization and measurement of biological applications with large-phase objects.

In QPI techniques based on partially coherent imaging, weak object approximation and slowly varying object approximation are the two most general approximation models. Unlike weak object approximation, which requires the absolute phase of the object to be very small, slowly varying object approximation requires that the difference between the sample phase and its neighborhood is much less than 1 rad. Although these two approximations are defined from different perspectives, they can be equivalent when some other approximations are introduced [15]. For example, the quantitative phase of a sample with a slowly varying phase and weak absorption can be solved by the weak phase deconvolution algorithm, which means that the deconvolution solver based on weak object approximation can be relaxed to a slowly varying object. It is also important to mention that, to simplify the model, we analyzed the DPC QPI of a pure phase object, whose nonlinear errors originate only from the higher order terms of the phase. As for samples with both absorption and phase [$t(\mathbf{x})=a(\mathbf{x})e^{j\varphi (\mathbf{x})}$], their absorption can be eliminated by the difference calculation in DPC under weak object approximation. However, the nonlinear error between their acquired intensity and the ideal one used for deconvolution is more complex, containing higher-order convolution terms of the sample absorption as well as the cross-convolution terms of the absorption and phase. A similar iterative deconvolution method can be used to eliminate the nonlinear error in the acquired image so that the input to the deconvolution solver is as close as possible to the intensity linearly related to the phase. In the future, we will explore different versions of the iterative deconvolution algorithm for accurate quantitative phase reconstruction of objects containing non-uniform absorption.

Acknowledgment. We thank Ran Ye from Nanjing Normal University for discussing the writing of the manuscript.