1. Introduction

Wave front sensing, the recovery of the phase of a (partially) coherent wave function has many applications in optics, electron optics, acoustics, seismology, astronomy, and many other areas of science. The detection of wave fronts can be implemented in many different ways [1]. This paper presents a novel method for deterministic phase retrieval from focal series using the transport of intensity equation (TIE). Since its inception [2], the TIE has gained considerable attention as it allows to perform deterministic phase reconstruction from relatively simple experimental set-ups. Moreover, it represents a complement to iterative reconstruction methods such as the ones based on the Gerchberg-Saxton (GS) algorithm that inherently present difficulties retrieving slowly changing features in the phase information [3]. An algorithm based on or incorporating information from the TIE could circumvent this limitation, retrieving the low spatial frequencies associated with these features.

The TIE originates from the conservation of the Poynting vector for monochromatic fields and paraxial imaging in linear homogeneous transparent media [4]. It provides a unique solution for coherent fields with a non-zero intensity distribution [5], being generally compatible with observations made away from sharp pupils and under certain conditions with observations including pupil boundaries [6,7]. This means the applicability of the TIE encompasses a wide range of applications based on or comparable to conventional optical microscopy. Indeed the main niche application of the TIE is in the analysis of through-focus image series, consisting of images acquired at different planes of defocus, both in optical and electron microscopy [8–10], and also other types of imaging devices. Due to its generality, the TIE is finding today an even broader range of applications; for instance as part of more sophisticated holographic reconstruction algorithms [11–14], or even in three-dimensional holographic reconstruction schemes [15]. These recent advancements show that the TIE constitutes not only a popular phase retrieval method, but also a quite active field of research with exciting possibilities.

In practice, the TIE presents issues when dealing with experimental data given the intrinsic difficulty to estimate the axial intensity derivative and its noise amplification properties. These issues have a greater impact on the low spatial frequency range, preventing the TIE to reliably retrieve slowly changing features. Our aim in this work is to review the origin of these issues and to present a novel TIE algorithm that tackles them. Our multi-focus algorithm (MFTIE) has in common with other existing methods that it uses the three-plane derivative in geometrically spaced through-focus image stacks [16]. This method is based on the fast Fourier transform (FFT) inversion of the TIE, it has been extensively studied in the literature and its response in the Fourier domain is well known. We propose a new method for the Laplacian inversion of the three-plane defocus derivatives, using overlapping and normalized band pass filters that take partial spatial coherence into account. The properties of these filter increase the robustness against noise of the MFTIE and allow to redefine existing methods for the sampling of information in the spatial frequency domain. We also validate our results using standard tools and compare with the Gaussian-Process (GPTIE) algorithm [18], a state-of-the-art code that uses multi-band smoothing to go beyond the three-plane TIE approximation.

The methods in this paper have been implemented into a basic in-line holography suite using the Python-based library hyperspy [19,20]. Graphical processing unit (GPU) code has also been included through the use of the PyCUDA library, mainly to boost the performance of key FFT based calculations.

2. MFTIE algorithm design

Let us first give a quick overview of the theoretical framework of the TIE. We start by considering monochromatic plane-wave coherent illumination normally propagating along the +z axis,

Where r = (x, y) are coordinates in the planes that contain the wave-front and k is the wavenumber. In this formulation the through-focus intensity, defined as the plane-wave modulus I_z = |u_z|, is experimentally measurable. Meanwhile the phase, defined as Φ_z = arg (u_z), cannot be experimentally obtained. However, a relationship between the phase and intensity may be established by considering a paraxial approximation;

This expression is known as the TIE, where the refractive index of the medium, n, may also been included through k = 2πn/λ (n = 1 for vacuum). Moreover, z sub-indices have been dropped for clarity. ∇_⊥ is the 2-dimensional gradient along the directions normal to z, the axis of propagation of the wave. In the remainder of this paper we will also drop the ⊥-subscript on the gradients, i.e. setting ∇ = ∇_⊥. The TIE is a second-order elliptical equation, relating the through-focus intensity and phase of the plane wave with the intensity variation, ∂_zI. It gives access to the phase at a given defocus if the values of the intensity and its derivative at that plane are know. It is desirable to find methods that solve the TIE and that include an estimation of ∂_zI from intensity measurements, in the best case without using any specialized instrumentation other than a microscope capable of through-focus imaging. An estimation of this derivative can be obtained by only using two images acquired at equally under- and over-focused planes and calculate the central difference,

Where the defocus planes lie at the same z = ±δz absolute defocus values, and (δzI)₀ is estimated at z = 0 using the corresponding defocused image planes. In principle, Eqs. (2) and (3) allow obtaining the phase from only three intensity measurements; the in-focus intensity and to defocused planes. Several works have shown variations on how to use this three-plane approximation to the TIE solution to obtain phase information from three-plane focal series [8,10]. Figure 1 shows the result of performing three-plane TIE phase reconstruction from a noiseless image simulation. In this case, the retrieved Φ₀ is almost identical to the original Φ, with no trace of the non-homogeneous intensity distribution I.

 

Fig. 1 Panels (a) and (b) respectively depict the amplitude and phase of the input wave used for simulations in this paper, while panel (c) is the phase retrieved from three-plane TIE using a noiseless simulated dataset.

Download Full Size | PDF

In practice, this three-plane TIE derivative will rarely retrieve a satisfactory TIE phase. The reasons for this inadequacy are non-linear terms of the derivative and corruption of low-spatial frequencies from noise in the images [9]. These limitations are evidenced by the artifacts that appear in the reconstructed phases. Figure 2 presents examples of phases recovered by the three-plane TIE at different values of defocus δz, from a dataset that includes random noise in the form of uncorrelated additive white Gaussian noise (UAWGN). These images portrait the expected behavior for three-plane TIE; the high and low spatial frequency information shifts from small to large defocus. This is evidenced in the cloudy low-spatial-frequency artifacts in panels (a) and (b) and the loss of resolution in panel (c).

 

Fig. 2 The images at the left present three-plane TIE phases obtained for each defocus pair in a simulated dataset with added noise (SNR = 20 dB) and, from (a) to (c), defocus values δz = 1, 5 and 10 μm, respectively. Each recovered phase has been compared to the original phase in the simulation using FRC and RMSE, depicted in the panels to the right using blue and red lines, respectively. Light blue lines correspond to negative parts of the FRC, that have been flipped for representation purposes. Dashed red lines indicate the theoretical noise error, ${\text{RMSE}}_n^2$, calculated for this level of noise and these defoci.

Download Full Size | PDF

These properties can be further explored by studying the frequency response of the TIE. For the phases recovered in Fig. 2, validations tests have been computed; the Fourier-ring correlation (FRC) and differential root-mean-square error (d_qRMSE²). Both tests are performed in the spatial frequency domain and integrated to unitary rings around the origin |q| = 0. The FRC measures the normalized correlation between two images for concentric rings in the spatial-frequency domain [21]. A FRC value of 1 indicates perfect correlation for the information at that q-ring; values between 1 and 0 indicate worse correlation. Additionally, negative values may signify contrast inversion, with a value of −1 indicating perfect contrast inversion. In this image, the FRC plots show that the retrieval of information at high-spatial-frequencies is dominated by the non-linear properties of the aberration function. The amount of information transferred into images projected at a given defocus is measured in the phase transfer function (PTF) [9]. In this sense, FRC drops dramatically in regions where the PTF is near zero, presents contrast inversion where the PTF is negative and degrades generally as the correspondence with the paraxial approximation decreases. Further dampening is introduced by the spatial coherence envelope.

Meanwhile, the d_qRMSE² test is based on the phase quadratic error [22], that is found in the spatial frequency domain by radial integration of these graphs. Following the power theorem, the total quadratic error in the real and frequency domains is equal for a large enough detector. The d_qRMSE plots in Fig. 2 show that regions at low spatial frequencies (where FRC also degrades) present an exponential growth of the reconstruction error term for each spatial frequency, as they approach the origin. Note that the exponential growth of RMSE² at low spatial frequencies is predicted by a theoretical model for the noise induced error term, ${\text{RMSE}}_n^2$, that is also included in the plots. The properties of this model are explained further ahead in this text as they are are central to the spatial frequency sampling scheme of our algorithm.

These limitations of the TIE cannot be easily tackled, for instance, by selecting an adequate defocus value. Indeed, using a small defocus step to improve the estimation of the derivative reduces the signal, increasing the impact of noise. Meanwhile using medium to large defocus steps only gives access to band-width-limited phase information. An effective solution to these issues demands the design of a MFTIE algorithm, that uses different defocus planes to obtain an improved phase reconstruction. To implement this solution, two or more three-plane TIEs can be combined using band-pass filters [9,16]. In this method, Eq. (2) is solved for several image pairs acquired at different defocus. The resulting partial phases are combined into a single phase by using band-pass filters, since considering images obtained at (relatively) large and small defocus, the former will provide low and the latter high spatial frequency information, respectively. To our knowledge, even in the recent implementations of this method the shape of these filters is restricted to sharp non-overlapping masks in the spatial frequency domain. Information sampling in this domain is additionally restricted by setting a fixed separation of defocus planes. These practices are difficult to implement and problematic in experimental application, as they can lead to ringing and other artifacts in the reconstructed phases. A more flexible methodology can be implemented by using overlapping filters, additionally taking into account partial spatial coherence and using an efficient sampling scheme without forcing specific defocus steps or needing a Tikhonov regularization parameter.

2.1. Inverse Laplacian operator

Considering Eq. (2), two in-plane gradient inversions have to be performed and an estimation of the axial intensity derivative should be provided to recover the phase at the in-focus plane using the TIE,

Where Φ₀ and I₀, are the phase and intensity at the in-focus plane, respectively. In the right-hand side expression, the inverse gradients have been re-expressed using inverse-Laplacian operators, ∇⁻². In the spatial-frequency domain, gradient and inverse-Laplacian operators are straightforward to implement. Defining the spatial frequency vector q = 2π/r and its modulus, q = |q|; then the operators are found, ∇̂ = −iq and ∇̂⁻² = −q⁻². These expressions are commonly applied for TIE inversion using FFTs.

In order to extend Eq. (4) to the case where multiple estimations of the derivatives are considered, our MFTIE algorithm is based on a specially designed inverse-Laplacian operator, ${\widehat{\nabla }}_{\delta z}^{-2}$. This operator is formulated in a such a way that the partial phases derived from the TIE inversion of the available through-focus three-plane derivatives, (∂_zI)_0,δz, add-up in the following fashion,

2.2. Spatial frequency thresholds

In order to implement the multi-plane solution described above, high and low cut-off frequencies, q_H and q_L that define the H_δz band-pass filters need to be established. Since our MFTIE algorithm also includes a normalization, the cut-offs need not to be coupled into a fixed progression that ensures non-overlapping filters. Hence, separate expressions for these cut-offs can be implemented attending to different requirements.

Looking first into the definition of q_H, some authors compare the three-plane TIE with the PTF derived from a Fresnel-Kirchoff (FK) diffraction model in a weak phase object approximation to find one of such requirements [9],

In principle, ∊_ϕ is a small number (∊_ϕ ≪ 1) that sets the cut-off for the spatial-frequency range q > q_H influenced by non-linearity. However if applied rigorously by setting a small threshold (∊_ϕ < 0.1) it removes information that can be used by the TIE solver if the value of ∊_ϕ is set to relatively large numbers (∊_ϕ > 0.4). If we consider for instance the defoci in Fig. 2 (δz = 1, 5, 10μm), the high-frequency cut-offs predicted by Eq. (9) for an ∊_ϕ = 0.1 are q_H ≃ 4.4, 2.0, 1.4 μm⁻¹, respectively. When compared with their respective FRCs in Fig. 2, it becomes obvious that these cut-offs underestimate the information content of the three-plane TIE solutions. Equation (9) is however useful if one takes into acount that it sets the sensitivity of the TIE solver at a fixed value for different defoci [18]. In fact, the maximum phase information transfer in the FK model corresponds to $\text{sin}\left[{\chi }_{\delta z}^{FK}({∊}_{\varphi }\simeq 0.41)\right]=1$. This property has proven useful for some applications, such as to develop optimum-plane selection strategies for the TIE [16,18].

In practice, cut-offs set by even higher values of ∊_ϕ allow to retrieve an even wider range of phase information. However, the correspondence with the FK model will then not be guaranteed. In this sense, a different cut-off relationship can be formulated attending to the Rayleigh-Sommerfel diffraction (RS) model. From this model, a second-order term in the paraxial approximation lies at the foundation of the TIE,

Where a small value of ∊_χ defines a range of spatial-frequencies q < q_H for which the approximation leading to the TIE (FK model) holds. Notice that ∊_χ accounts for the lag in the RS-PTF respective to FK-PTF. For the defoci in Fig. 2 a relatively small value of ∊_χ = 0.25 rad results in cut-offs q_H ≃ 8.0, 5.3, 4.5 μm⁻¹, that allow a wide range of phase information in the frequency domain to be retrieved.

Let us now focus our attention into the definition of the q_L cut-offs. Considering the q⁻² dependency of the inverse Laplacian operation, it becomes obvious that TIE inversion amplifies noise at low-spatial frequencies. Apart from the slow-varying artifacts in the recovered phases, the effect of noise is evidenced in our simulations when comparing to the original phase using validation tools. First, looking at the FRC we can see the degradation on the correlation at low spatial frequencies, where also the RMSE² increases exponentially. This behavior is well-known and can be used it to limit the impact of noise [16], by restricting the amount of RMSE² going into each partial-phase in the MFTIE. Following their framework, we consider the integration of the noise contribution to the error in the phase for a three-plane TIE with a sharp band-pass filter;

Where A_δz is proportional to the ${\text{RMSE}}_n^2$ in an unfiltered three-plane TIE solution using a square detector of size L [22]. The comparison of this ${\text{RMSE}}_n^2$ theoretical model with the actual RMSE² measured in our simulations is excellent (see Fig. 2, red dashed lines), and shows how the noise component dominates the overall error obtained at low spatial-frequencies. Taking this into account and realizing that the $q_H^{-2}$ term in Eq. (11) can typically be neglected (especially for δz ≫),

Hence setting ${\text{RMSE}}_n^2$ to a fixed value (maximum allowed error) could be used alone to estimate a low frequency cut-off to extract meaningful information from a three-plane TIE into a partial-phase solution. However, this method is still incomplete since A_δz depends internally on the standard deviation of the noise contribution, σ_n, which is difficult to determine experimentally with precision. To overcome this difficulty, we propose to use a homogeneous noise distribution sampling (HNDS) scheme in the spatial-frequency domain, which overrides this dependency by setting the lowest q_L (corresponding to highest defocus) to cover all or a fraction of the spatial frequencies present in our images.

To understand HNDS, first consider that by setting the $q_L^2$ cut-off for a sufficiently high value of defocus, Eq. (12) produces a SNR-dependent estimation of the noise-induced error in this partial-phase TIE solution, ${\text{RMSE}}_n^2\left(\text{SNR},\delta z\gg \right)=A_{\delta z}{q}_L^2$. Conversely, HNDS works by fixing this value for all the partial-phases in our MFTIE calculation, hence allowing a sampling rule in the spatial-frequency domain to be found, or, in other words, all corresponding q_L(δz) are found independently of SNR. To do so in the MFTIE algorithm, a base ${\text{RMSE}}_n^2$ is first estimated as the noise-induced error for the highest value of defocus in the image stack, Δz, and q_L = q_m, where q_m is the minimum spatial-frequency in the images, given by detector size and magnification,

Alternative sampling, different from the base HNDS that covers all spatial frequencies in the images, can be chosen by selecting multiples of this base level of noise using a multiplicative factor, ρ, for which a corresponding cut-off progression can be expressed;

In particular, the base HNDS that covers all available q is recovered for ρ = 1, satisfying the relationship; $q_L^2(\mathrm{\Delta }z)=q_m^2$. Increasing or decreasing ρ sets HDNS with wider or narrower band-pass filters, respectively. This allows to include the lowest spatial-frequencies from smaller defocus (ρ > 1), or, alternatively, to completely remove a range of those from the calculations (ρ < 1). Note that this method is easily automated, depends only on the value chosen for ∊_χ and ρ, and has the additional advantage that given a good estimation of SNR, the amount ${\text{RMSE}}_n^2$ included in each partial-phase can be estimated, indicating the noise level of the reconstruction.

3. Results

3.1. Simulated datasets

Figure 3 presents a graphical example of this methodology; panel (a) is showing the phase retrieved from the three through-focus derivatives presented in Fig. 2. For the MFTIE algorithm, the q_H and q_L cut-offs used were set by Eq. (10) with ∊_χ = 0.25 rad and Eq. (14) with ρ = 1 (and SNR = 20 dB), respectively. The estimated noise level of the corresponding HNDS filter progression was ${\text{RMSE}}_n^2\simeq 560{\mathit{rad}}^2$. Note that the size of the portrayed images is of 512 ×512 pixel, and reflection padding was used to replicate both dimensions and ensure periodic boundary conditions for the FFTs.

 

Fig. 3 The images at the left present the phase of the complex waves obtained from reconstructions carried out using, in the first two rows, MFTIE and, in the last row, GPTIE. All the recovered waves were refined using 25 iterations of the Gerchberg-Saxton algorithm. Panel (a) is showing the results obtained using only 3 defocus image pairs at δz = 1, 5 and 10 μm. Panels (b) and (c) use a dataset of 21 planes log-spaced between 1 and 150 μm, with a SNR = 15 dB. Blue and red lines are used to depict FRC and RMSE validations of the recovered phases, with solid and dashed lines indicating if the phase before or after the refinement was used, respectively.

Download Full Size | PDF

The FRC and differential RMSE² tests for this calculation are also included in Fig. 3, using solid blue and red lines, respectively. Both tests validate the ability of MFTIE to integrate the through-focus information into an improved result. To obtain the best possible results, the MFTIE phase was used together with the in-focus intensity to initialize a GS-like iterative refinement loop [11,23]. With only 5 iterations of this loop, small errors in the middle and high spatial-frequency range are corrected, improving the resolution of the reconstruction. This is visible in the corresponding validation tests, depicted in the same panel as the MFTIE ones using dashed lines.

An important property of MFTIE algorithms is that improving through-focus sampling, for instance adding more images at larger defocus, makes the result more robust against noise. Although this premise is in theory true for most implementations, in practice the information at large defocus is more complicated to include in the reconstructions because of partial spatial coherence. Figure 3(b) and (c) compare the results obtained from a noisy and more incoherent dataset (SNR = 15 dB, α = 0.125 rad) using MFTIE and GPTIE algorithms, respectively. In this dataset, a logarithmic progression of 10 defocus image pairs centered around the in-focus position was used. In panel (b), a MFTIE calculation is presented, which was performed using the same input values for ∊_χ and ρ as the one in panel(a), i.e. producing a corresponding HNDS spanning the whole frequency range. However, since more defocus planes were included than before the estimated noise level of the corresponding HNDS filter progression is smaller, ${\text{RMSE}}_n^2\simeq 25\hspace{0.17em}{\mathit{rad}}^2$. Indeed, the validation tests for this phase reconstruction show even a slight improvement, especially at the low spatial frequency range. Note also that in both presented cases the GS-like refinement does not improve the performance in the validation tests of the reconstructed phase at middle to low spatial frequency range. This iterative algorithm is able to converge quickly the high- and medium-range spatial-frequencies [23], but is not able to do so for the low spatial-frequencies, an unavoidable behavior for this kind of algorithms [3].

To better assess the performance of MFTIE, phase reconstruction has also been performed using the same dataset and state-of-the-art GPTIE algorithm [18]. This sophisticated algorithm uses Gaussian-process smoothing to perform multi-band interpolations of the intensity derivatives and extract information based on the weak-object approximation FK-PTF. The results obtained after these procedure are presented in Fig. 3, panel (c). The validation tests evidence that while high-frequency information is limited in MFTIE by the three-plane TIE at the smallest-defocus, GPTIE improves the resolution by efficiently including the information from several planes into its derivative estimation. Indeed, the GPTIE multi-band derivative included a high number of cut-off frequencies (N_σl = 50), and is less prone to errors in the middle to high spatial-frequency range, evidencing the robustness of this algorithm.

However, the improvement is only slight since the correspondence of the paraxial approximation and partial spatial coherence are a limiting factor for both TIE-based algorithms and the smallest included δz = 1 μm already offers a good sensitivity in the high spatial-frequency range. Moreover, the middle to low spatial frequency range shows degradation of FRC at some regions and an increased error level, relative to the MFTIE result in the panel above. Indeed, since GPTIE does not allow for overlapping band-pass filters, the response of this algorithm at a given spatial frequency range depends only on a single interpolation which can be prone to errors, later amplified by Laplacian inversion. Moreover, since GPTIE does not take into account partial spatial coherence, highly-defocused images containing low-spatial frequency information less impacted by ${\text{RMSE}}_n^2$ will not be correctly recovered. Finally, the GS-like refinement algorithm that we have used takes into account partial spatial coherence [11] and hence is also able to correct for some of the deficiencies in the result by GPTIE. However, the lowest spatial-frequencies were not completely converged after 5 iterations and some corruption in this range persists.

Finally, take into account that the use of deterministic phase reconstructions such as MFTIE and GPTIE should be aimed at this low spatial frequency range since it has been proven that simple iterative refinement algorithms need thousands if not more iterations to reconstruct this information [3]. Note that slow-changing features can have a high impact on the measurements performed on the recovered phase information, such as tissue thickness in biological samples.

3.2. Experimental dataset 1: optical microscopy

In this section, we analyze an experimental image stack in order to benchmark the performance of our algorithm in a typical task; recovering the phase information from a biological sample. The HeLa cancer cells in this sample are fairly transparent, with a diameter of ∼ 20 μm and details below the μm-range. A dataset consisting of 20 logarithmically spaced defocus image pairs was acquired using a Zeiss Axiovert 200M inverted microscope equipped with an immersion objective (NA = 1.35 and n = 1.518). These images were initially processed by median filtering, which helps to eliminate speckle noise. Then, a region-of-interest 700 × 700 pixel wide was selected and linear-ramp padding for 350 pixel more were added in each direction. The resulting 1050 × 1050 pixel dataset was processed using the MFTIE and GPTIE algorithms, both followed with 25 iterations of GS-like refinement.

The results from these procedures are depicted in Fig. 4. In the first column of this image, the MFTIE and GPTIE results before iterative refinement show different degrees of blurrying due to the poor resolution of the TIE phase. However, after iterative refinements both solution reach visually similar phase images. The same input parameters as in the last section were used for both calculations, except for the convergence angle. Several prior calculations were used to fit this parameter, as shown in Fig. 5 top panel. This fitting is done by performing MFTIE reconstructions (with or without GS-like refinement) for different values of the convergence angle, calculating the r-value test for these and searching for the minimum. The r-value parameters are obtained from comparing the experimental intensities, I^EXP, with intensities simulated using the recovered phase (or wave, in the case GS-like refinement was included), I^SIM;

 

Fig. 4 In the first column, panels (a) and (d) show the phases obtained from of HeLa cells using MFTIE (top) and GPTIE (bottom) in-line holography algorithms, respectively. The second and third columns show the phase and amplitudes of the complex wave obtained after 25 iterations of GS refinement of the TIE results.

Download Full Size | PDF

 

Fig. 5 Validation tests with the experimental images for the phase reconstructions in Fig. 4. In the top panel, r-value tests for the validation of the MFTIE convergence envelope. Here, blue circles and black crosses indicate the r-value calculated before and after GS refinement, respectively, performed for different values of the convergence angle. In the bottom panel, the χ²-tests, with black circles and red squares representing MFTIE and GPTIE algorithms and dashed or solid lines indicating whether the validation was performed before or after GS refinement.

Download Full Size | PDF

This process shows that the optimum is found for α ≃ 0.05 rad ≃ 3°, hence this value is used for all the MFTIE and GS-like refiment calculations and also further simulations regarding this dataset. Note that since the partial spatial coherence is solely dependent on the illumination this process does not need to be repeated every time the MFTIE reconstruction is applied on the same experimental set-up. Note that this procedure has produced an estimation of α and that, alternatively, an experimental measurement or parametric fit of the coherence envelope could be produced, but this is beyond the scope of this work.

Before GS-like iterative refinement, the MFTIE and GPTIE algorithms score r-values ≃ 3.96 · 10⁻³ and 4.22 · 10⁻³, respectively. The lower value of the MFTIE generally indicates that a better adjustment with the experimental information has been achieved. After GS-like refinement, the difference between both algorithms narrows, and this is also reflected in the obtained r-values ≃ 3.72 · 10⁻³ and 3.78 · 10⁻³, respectively. However, even after refinement, the reconstructed phases have some (faint) differences; slight variations in the regions around the cell nuclei are visible and some different slow-varying features are also present. Overall, the MFTIE cells look fuller measuring a slightly larger phase excursion while the GPTIE background looks more homogeneous (Note that the GPTIE algorithm is using a Tikhonov regularization parameter, here set to 5 · 10⁻⁴, while the MFTIE algorithm does not require defining such a parameter).

To shed some light on the origin of these differences, and since we do not have access to the ground truth, we have performed two additional validation tests; a through-focus χ²-test and an even-odd FRC test, which are presented below. Figure 5 bottom panel presents the results from the χ²-test. Similarly to the r-value test, each point in these plots corresponds to the comparison of the experimental and simulated intensities, calculated for all defocus positions, using the following the formula;

Looking first at the results before GS-like refinement (dashed lines), the correspondence of the GPTIE seems to be better for the first few defocus planes, closer to the in-focus plane. Then, there is a small region in the medium defocus range where the correspondence of both algorithms is similar. Finally, the MFTIE obtains a better correspondence for the majority of defocused planes, from medium to large defocus ranges. Looking now at the results after GS-like refinement, it seems the situation is somehow improved and more balanced, with similar and lower overall χ² values. This is in agreement with the interpretation of the r-value test presented before. However, note that the correspondence for images in the large defocus region, at about δz > 10 μm, is still better for the procedure using the MFTIE.

The final piece of information is found in the even-odd FRC tests. This validation is based on the already presented FRC-test. However, it is not possible to directly perform this test since the original phase is unknown in this case. In the even-odd FRC test, the original dataset with 20 defocus image pairs is divided into two datasets with 10 defocus image pairs, each alternatively using the even- and odd-numbered through-focus images of the original dataset. From these two new datasets, two different reconstructions of the phase are obtained by applying the same methodologies (same algorithms and inputs). The solutions obtained in this fashion are called the even and odd phases. These can be compared using FRC, obtaining curves which evidence the degree of correlation between alternative phase information.

These FRC-tests are portrayed in Fig. 6, for MFTIE and GPTIE algorithms; before and after successively applying GS-like refinement. The shape of the obtained curves is reminiscent of the FRC-tests for the simulated datasets (see Figs. 2 and 3). In particular, the GPTIE reconstruction is showing a greater range of recovered high-frequencies (not shown), meanwhile MFTIE has an overall better score at medium spatial-frequencies. Moreover, the curves obtained after 25 iterations of GS refinement again evidence that both algorithms reach a similar convergence solution. Most interestingly, looking at the low spatial frequency region, the curves for the GPTIE algorithm both before and after GSR show some dips in the FRC score. In contrast, the MFTIE algorithm has also an overall better score in this range although it also seems to struggle with the very lowest spatial-frequency. This last effect is most probably associated with the inherent difficulty at recovering this spatial frequency due to the inadequacy of the FFT boundary conditions and the padding.

 

Fig. 6 Even-odd FRC validation tests of the MFTIE (black) and GPTIE (red) phase reconstructions in Fig. 4. The top and bottom panels include the validations performed before and after GS refinement, respectively.

Download Full Size | PDF

Both FRC and χ² tests are in good agreement since phase information at the low-frequency range is recovered mostly from the high-defocus intensity data. Alternatively, the simulation of the defocused planes also depends on the quality of the information at this range since a coherence envelope is included in the simulation model. These results indicate that the GPTIE although being a sophisticated and noise-robust algorithm, can have trouble recovering the very low spatial-frequency range of the phase information. This deficiency can be the cause for measuring slightly lower values for the phase excursion in the HeLa cells. In contrast, the MFTIE shows how the impact of noise-induced artifacts can be greatly reduced using a rather simple algorithm. Since the use of neither of these algorithms precludes completely the apparition of these artifacts, validation tests have to be performed as they are key to understanding the obtained results.

3.3. Experimental dataset 2: electron microscopy

We also performed a transmission electron microscopy (TEM) experiment, from a sample containing gold nano-particles (NPs) sputtered over a lacey carbon film. For these experiments, a JEOL JEM-2200-FS microscope equipped with an in-column energy filter was used. A through-focus image series was acquired, operating in medium-resolution TEM mode (pixel size ∼ 0.42 nm) at 200 keV and using a 10 eV energy-slit centered on the zero-loss-peak to select elastic scattering. The original dataset is composed of 10 image pairs with defocus values from ∼ 40 to 2000 nm and an additional in-focus image. From this dataset, a 840 × 840 pixel (350 × 350 nm) region-of-interest centered around a hole in the carbon film measuring ∼ 250 nm in diameter was selected, portraying gold NPs in a range of sizes ∼ 1−20 nm. In an initial step, this dataset was processed using the full resolution wave reconstruction (FRWR) code [11,24], that corrects the image distortions that appear in TEM through-focus imaging by iterative alignement of the images. From the produced aligned dataset, 8 image pairs were then selected and used as an input for the MFTIE algorithm.

The phase image resulting from MFTIE is depicted in Fig. 7(a). In this image, the phase result seems to be affected by a low-frequency artifact. This appears as a standing wave superimposed to the phase information, attending at the high intensity regions at the corners and the disk-like contrast in the hole region. The origin of this artefact may be due to remaining inelastic scattering contributions to the individual images, distortions of the projector lens system, limited isochromaticity of the imaging energy filter, or other limitation of the experimental setup. In the MFTIE calculation, a low-frequency threshold was chosen to preclude the artifact from masking completely the relevant phase information. The thresholds were calculated for each defocus value using HNDS, including spatial frequencies down to a minimum q_m ≃ 2π(150 nm)⁻¹ and α = 3 mrad. Even higher thresholds will completely remove the artifact but also remove important low spatial frequency information. For the result shown, 280 pixel of padding area with a constant value were added in each direction to the measured region-of-interest, resulting in a 1120 × 1120 pixel dataset. The measured R-value for this result is 0.065, different choices for this padding, such as linear ramp or mirror, did not produce a completely satisfactory result.

 

Fig. 7 Phase reconstruction from an experimental TEM dataset of sample with gold NPs (bright) sputtered over lacey carbon featuring a hole in the middle (dark). Panel (a) is showing the MFTIE result, obtained using HNDS down to q_m ≃ 2π(150 nm)⁻¹. Panel (b) shows the phase recovered after 20 iterations of GF-refinement using the former result as input. In this refinement, spatial-frequencies below q_m ≃ 2π(8 nm)⁻¹ were iteratively flipped to obtain a phase with a smooth gradient over large portions of the measured and pad area. Panel (c) is showing the phase retrieved by non-linear refinement using FRWR, which is able to retrieve the sharp phase excursion between vacuum and sample.

Download Full Size | PDF

Other than noise amplification, artifacts in an FFT-based algorithm can appear in situations where it is difficult to fulfill the boundary conditions appropriate for Fourier analysis. This is arguably our case, because there is a non-uniform distribution of intensity in the border of the measured region. This makes it necessary to choose an appropriate padding for the image, which is a difficult and error prone task. In this case, artifacts appear because the phase problem (inverse-Laplacian) is ill-conditioned due to non-uniqueness of the solution. A way of correcting this inadequecy is to use gradient-flipping refinement of the measured derivatives, which has proven able to removed these type of low-frequency artifacts [13,14]. This procedure is related to charge-flipping, well-known in the field of crystallographic structure determination from diffraction data [25]. Such algorithms are able to overcome local optima in non-convex optimization spaces. This makes them adequate for solving the problem of non-uniqueness and ill-conditionedness of the solution when boundary conditions are not specified. In the version of GF used here, we refine our measured three-plane derivatives, (∂_zI)_0,δz, by separately flipping the sign of the gradients of the partial phases, Φ_0,δz, in each iteration depending on the value of their l₁-norm (absolute value). In this case, the input phase for this refinement was the one depicted in panel (a), and appart of the GF step, the rest of the iterative loop is purely based on MFTIE.

The result from GF-refinement is presented in panel (b) of Fig. 7. The resulting image shows no trace of the standing wave-like artifact, and the phase values inside of the hole region are now constant. As a direct result, the measured range for the phase excursion is reduced from ∼ 4 to 2.4 rad. The measured R-value is also slightly reduced to 0.064, indicating a better agreement with experimental data is obtained. Also in Fig. 7, panel (c), the phase recovered from the FRWR algorithm is presented. This algorithm also includes GF, but applied to the phase resulting from iterative propagation with a transfer function, similar to the GS-refinement explained in the previous sections. For FRWR, the measured R-value is 0.0178, the lowest of all, indicating the best correspondence with experimental data. Indeed, the difference between the hole and film region is now accentuated and the vacuum appears darker. The difference between the pure TIE result and the FRWR one in panels (b) and (c), respectively, portraits the inability of the TIE to retrieve sharp phase changes which is inherent to the formulation of the TIE [5]. In turn, the similar features and dynamic range of both reconstructed phase images show that the pure TIE solution includes valuable information obtained in a deterministic fashion.

The presented experiment corroborates that MFTIE can also become useful in medium-resolution TEM, as has been advanced by other authors using normal three-plane TIE [10]. On the one hand, the correspondence between the theoretical framework of the TIE with the actual physics of image-formation in TEM is poorer than for conventional optics. This is because in electron imaging, inelastically scattered electrons may still contribute in an incoherent manner to the imaging process, while they are absorbed in the case of light optical microscopy. In addition, more severe image distortions and non-linear effects (apart from dynamical and crystalline diffraction effects) play a more important role in the formation of TEM images and are not accounted for by the framework of the TIE. On the other hand, TIE methods are able to provide useful low-frequency information that can have relevance in samples imaged at medium-resolution, as the presented gold NPs. For those cases, our MFTIE algorithm provides an efficient way of integrating the information from many defocused planes, also including partial spatial coherence in its formulation.

4. Conclusion

This paper presents a simple TIE-based algorithm that recovers the phase information from a through-focus image series, designed with noise robustness in mind and including the contribution of partial spatial coherence to image formation. It has been shown how the MFTIE reconstruction algorithm allows to recover a wide range of spatial frequencies given that images at large and small defoci can be acquired. Additionally, there is no need for acquiring large datasets with uniformly spaced defocus, since an efficient through focus sampling of information is possible with logarithmically spaced defocus. In our work-flow, an initial TIE guess for the phase is followed by a flux-preserving iterative GS-like refinement. This last step allows to improve the reconstruction of medium to high spatial frequencies with little additional computational cost. Our simulated benchmarks have been complemented with the experimental results presented in the last sections, from optical and electron microscopy experiments. Since both the MFTIE algorithm as well as the GS-refinement account for partial spatial coherence the correspondence with experimental data is also improved.