1. Introduction

X-ray imaging and tomography offers a set of advantages for biomedical and materials research: resolution, penetration power, high volume throughput, compatibility with hydrated specimen, and quantitative contrast mechanisms relatively unaffected by multiple scattering. The advantage of phase over absorption contrast for brain tomography has been realized early, since for the high photon energies $E$ of hard x-rays, the index of refraction $n(E, {\mathbf {r}})= 1 -\delta (E,{\mathbf {r}}) + i\beta (E,{\mathbf {r}})$ leads to dominating phase interaction compared to absorption. For low-$Z$ elements, in particular, we have $\delta /\beta \gg 1$, and a proportionality of $\delta (\mathbf {r})$ and the electron density. Hence, it has been a major research goal over the last two decades to exploit the small phase shifts in a wavefront imparted by matter for high resolution imaging. Different phase-contrast methods can be used to transform these phase shifts into measurable intensities patterns by ways of interference [1,2], based on different geometries and mechanisms [3,4]. Propagation-based imaging only requires the self-interference of the diffracted beam behind the object and the unattenuated or weakly attenuated primary beam [5]. It hence does not require scanning of the object or any optical element, which makes the method also dose effective. Finally, the imaging modality can cover a wide range of scales, from macroscopic (clinical) settings such as mammography [6] down to imaging sub-cellular spatial resolution below 50 nm [7] using synchrotron radiation.

A central challenge in propagation-based imaging has been the formulation of accurate and efficient phase retrieval schemes [8–10]. In the so-called direct-contrast regime, i.e. at small propagation distances $z$ or correspondingly large Fresnel numbers, linearization in $z$ can be used to derive approximate inversion equations [11]. On the contrary, the present work is concerned with the holographic regime of small Fresnel numbers, which is exploited for imaging with spatial resolutions in the order of 100 nm or below, where such inversion methods fail [5,12]. In this regime, where propagation imaging is also known as near-field holography (NFH), phase-sensitivity is high enough to ’pick up’ the small differences, for example encountered in unstained soft tissues. This enables a wide range of applications such as virtual histology of COVID-19 lung tissue [13], structural imaging of microscopic catalyst particles [14], observation of expanding cavitation bubbles [15] or resolving neuronal structures throughout entire Drosophila melanogaster brains [16]. For NFH, an inversion based on the known contrast transfer function (CTF) for optically weak objects has been proposed twenty years ago [10,17,18], using holograms data recorded at different measurement planes and Fresnel numbers $F_1, \ldots, F_J$. This so-called CTF phase retrieval has been extremely successful, and has found widespread use also due to its numerical simplicity and efficiency.

However, CTF fails or has limited image quality in many cases, when the basic assumption of a weakly varying phase is not fulfilled [19–21]. For example, interior or exterior interfaces can result in large phase gradient which often substantially deteriorate the image quality, even if the bulk of the sample is characterized by smaller phase gradients. To treat data with larger phase gradients, iterative phase retrieval schemes have to be considered that overcome the limitation to the linear CTF regime. Such nonlinear reconstructions in NFH, using gradient descent algorithms, have been reported over a decade ago [22,23]. Regularized Newton methods have also been proposed [24]. More recently, Hagemann et al. introduced a simple scheme of alternating projections (AP), which operates on the same input data and constraint sets as the conventional CTF-based phase retrieval [25]. As in far-field coherent diffractive imaging (CDI), iterative projection algorithms [26,27], which cycle between measurement and object plane, do not rely on linearization and are well suited also for phase retrieval in the NFH regime. In addition to the measured data, a priori constraints (prior information) on the imaged object, such as finite support, separability, or range constraints on the numerical phase values can be used [28,29]. However, single distance acquisitions of extended (not compactly supported) specimen, which in practice are of particular relevance for tomography, often suffer from stagnation and noise-induced artifacts as phase retrieval becomes more strongly ill-posed in this setting [30,31]. In fact, without a compact support as in [7,28], AP-type algorithms often only yield low-quality reconstructions, unless the data acquisition is extended by longitudinal (i.e. parallel to the optical axis) or lateral scanning. Moreover, iterative reconstructions are typically numerically expensive, which is prohibitive for large tomographic data sets. In order not to bottleneck tomographic imaging pipelines, computation times for image-reconstruction should not exceed the typical image-acquisition times (${\lesssim }\;{1}\;\textrm{s}$). On the contrary, iterative projection algorithms applied to holograms with $2000^{2}$ pixels typically require thousands of iterations and take minutes to complete even on data center graphical processing units (GPUs) [25].

In this work, we propose an algorithm for phase retrieval in the holographic regime, denoted as NLTikh, using Tikhonov regularization to invert the full nonlinear forward model of NFH. The approach is a nonlinear generalization of the well-established CTF reconstruction, which no longer depends on a weak object assumption and on the other hand permits to impose support- and/or range constraints on the recovered phase-image. As an intermediate step, we also introduce a constrained CTF algorithm (cCTF) that retains the linear CTF model but adds flexibility in terms of incorporating constraints. The proposed methods overcome the aforementioned performance issues of other iterative methods by applying state-of-the-art optimization methods for minimizing the Tikhonov functional, significantly reducing the number of required iterations compared to previous approaches. Additionally, our algorithms make use of GPUs whenever available, while retaining the ease-of-use of CTF phase retrieval from the user’s perspective.

While the early phase retrieval algorithms have been reported solely in the literature without usable implementations, it has become customary to publish the source code as well [32–34]. We follow this trend and provide MATLAB implementations for the proposed methods.

2. Phase retrieval algorithms

The task of phase retrieval in NFH is to reconstruct the phase $\phi$ (and absorption $\mu$) from near-field intensity data $I$ measured in the detector plane:

Mathematically, CTF phase retrieval amounts to applying Tikhonov regularization to the linearized forward model under the constraint in Eq. (3), i.e. it can be written as

In the present work, we want to build upon this successful model and achieve a generalization, rather than a replacement. To this end, we propose two extensions: First, we introduce a CTF-based method which is more flexible with respect to incorporation of additional a priori constraints, including in particular support: $h|_{\mathbb {R}^{2}\setminus \Omega } = 0$ outside some known bounded domain $\Omega \subset \mathbb {R}^{2}$, and negativity: $\phi, -\mu \leq 0$, which results from the positive definiteness of electron density. Formally, such a constrained CTF (cCTF) reconstruction simply amounts to restricting the search set of the minimization problem in Eq. (4):

Second, and more importantly, we want to omit the weak-object assumption underlying to CTF phase retrieval, by replacing the linearized model $\mathcal {L}_{F}$ with the general nonlinear forward map $\mathcal {N}_{F}$. Accordingly, we aim to reconstruct $\phi$ by solving the following minimization problem:

While the proposed generalizations are hence conceptually straightforward, they bring about significant challenges on the practical, algorithmic side: contrary to (4), the minimizer in the constrained variant Eq. (6) can no longer be computed by a single-step FFT-based filtering operation as given by Eq. (4b). Moreover, the minimization problem in Eq. (7) even becomes non-convex, so that a global solution strategy is completely beyond reach. In both cases, we therefore have to resort to iterative algorithms to compute the sought minimizers.

To achieve reasonably fast termination, we employ the following methods:

Notably, all of these algorithmic tweaks act behind the scenes – from the user’s perspective, our cCTF and NLTikh implementations are as easy to use as standard CTF phase retrieval. The code is fully publicly available as part of the HoloTomoToolbox [32,33].

Relation to existing methods: The proposed algorithms bear similarities to several previous work in the literature. In [23,34,37], gradient-descent schemes are proposed to iteratively minimize the distance to the hologram data (or rather to square-root-data $\sqrt {I_j}$ in [37]) based on the nonlinear model, i.e. the first term on right-hand side of Eq. (7b). Moreover, the approach in [23] also includes warm-starting from a CTF-reconstruction and the algorithm in [37] allows to impose range- or support-constraints via projected gradient descent as in the current work. On the other hand, these previous methods both lack the ability to impose regularization as well as adaptive stepsizes to accelerate convergence, leading to as many as 10000 iterations for a single reconstruction in [37]. Perhaps the most closely related work to the present one is [22], where nonlinear Tikhonov regularization with a gradient-penalty (corresponding to the choice $\alpha (\xi ) = \textrm{constant} \cdot |\xi |$ in Eq. (7b)) is proposed. However, no support- or range-constraints are enabled and the algorithm involves a computationally expensive line-search strategy. Also the regularized Newton method in [24] is similar to the present approach, but involves nested (outer) Newton- and (inner) conjugate-gradient-iterations that spoil the computational performance. In conclusion, we find that our NLTikh algorithm stands out from previous methods in three respects: (1) it is flexible in terms of incorporating constraints, (2) builds upon experience gained with CTF phase retrieval, in particular retaining the established regularization strategy, and (3) is algorithmically optimized to keep computational costs as low as possible.

Concerning our constrained CTF method, we note that similar ADMM-schemes have been proposed for linearized phase retrieval with total variation (TV) regularization [21,38]. While such a regularization strategy is somewhat more advanced than our approach, we note that it has been so far limited to the weak object regime due to the algorithmic challenges it brings about.

3. Algorithmic details

3.1 ADMM for constrained CTF

We solve Eq. (6) by splitting it into the uncoupled subproblems of the unconstrained CTF (4) on the one hand and the constraints on the other hand. This splitting enables us to use an augmented Lagrangian method, namely ADMM [39,40]. The resulting iterations read

Computational costs: Noting that the summand $2 \sum _{j=1}^{J} (s_{F_j} - c_{\beta /\delta }\cdot c_{F_j}) \cdot \mathcal {F}(I_j-1)$ in Eq. (8b) can be precomputed during the initialization of the algorithm, we see that the computation of the ADMM-update (8) essentially requires one forward- and one inverse Fourier transform per iteration, in addition to some less expensive pointwise arithmetic operations.

3.2 Projected gradient descent for nonlinear Tikhonov

In contrast to the linear CTF setting, no closed form, even for the case without constraints ($A = \mathbb {R}^{2}$), exist for Eq. (7). However, Fréchet-differentiability of the nonlinear operator $\mathcal {N}_{F}$ allows us to compute the gradient of the nonlinear Tikhonov functional $\mathscr{T}_{\textrm{NL}}$ in Eq. (7). As derived in Appendix A, it can be computed as

We determine the minimizer of Eq. (7) using a projected gradient descent method. The iteration is a gradient descent step followed by the projection onto the constraints-set $A$:

We note that, by construction, the NLTikh-method reduces to CTF-phase retrieval in the linear weak object limit $\phi \to 0$. On the other hand, a CTF-reconstruction is numerically both more stable and less expensive to compute (also in a constrained setting!). We exploit this fact by using a CTF-phase as an initial guess for the NLTikh-algorithm, i.e. $\phi _0 = \phi _{\textrm{cCTF}}$.

The iterations are automatically stopped as soon as the (relative) residual gradient,

Stabilization of high frequencies: The linear CTF model is diagonal in Fourier space, i.e. the reconstruction of a certain spatial frequency only depends on the Fourier-mode of the same frequency in the hologram data $I_1, \ldots, I_J$. In particular, this means that low-frequency data-errors may not cause high-frequency artifacts in the recontruction and vice-verser. Switching from the CTF- to the nonlinear model $\mathcal {L}_{F_j} \to \mathcal {N}_{F_j}$ does not preserve this property. Hence, low-frequency background variations in the holograms, as often caused by imperfect flat-field correction, may result in oscillatory artifacts when reconstructing nonlinearly. We find that this undesirable effect can be effectively suppressed by adding a third level to the regularization defined in Eq. (5):

Computational costs: Up to computationally inexpensive pointwise arithmetic operations, evaluating one projected gradient-descent step Eq. (11) essentially requires $2J$ forward- and $J$ backward Fresnel propagations according to Eqs. (9) and (10). Consequently, the computational costs of an NLTikh-iteration is comparable to that of typical alternating-projection iterations. Note that neither the stepsize-rule Eq. (12) nor checking the stopping-criterion Eq. (13) requires significant computational effort, as $\textrm{grad}\;\mathscr{T}_{\textrm{NL}}(\phi _k)$ must be computed in every iteration anyway.

4. Performance tests and comparison to existing methods

The algorithms were tested with NFH data measured at the Göttingen instrument for nano imaging with x-rays (GINIX), installed at the P10 beamline of the PETRA III storage ring (DESY, Hamburg), at a photon energy of 8 keV (Fig. 1, Fig. 4) and 13.8 keV (Fig. 3). An x-ray waveguide was used as a spatial and coherence filter of the undulator radiation which was focused onto the waveguide by a Kirkpatrick-Baez (KB) mirror [44]. The exit of the waveguide forms a quasi-point source for holographic illumination. The geometrically magnified holograms of size $2048\times 2048$ pixels were recorded by two different fiber coupled scintillator sCMOS detectors (Photonic Science and Hamamatsu) each with pixel sizes of ${6.5}\;{\mu }\textrm{m}$, placed approximately at $z_{02}\;{\approx }\;{5}\;\textrm{m}$ behind the KB focus. All holograms considered as test data in this work are given by dark- and flat-field-corrected detector images. The experimental setup parameters for the test cases in the subsequent sections are detailed in Table 1.

 

Fig. 1. Reconstruction of a test object composed of ${15}\;\mu \textrm{m}$-sized polysterene spheres. The images plot the phases obtained by different algorithms: (a) CTF phase retrieval. (b) The proposed NLTikh-algorithm. (c) Like (b) but with an additional negative-phase constraint. (d) Like (c) but using only one of the four holograms. (e) Alternating projections (AP) result taken from [25]. (f) Line-scans through the reconstructions along the blue lines. The black-dashed line shows the theoretical phase shifts if ideal spheres and a literature value $\delta = 3.673\cdot 10^{-6}$ for polystyrene are assumed. Insets in the top-left corner show zooms of the red-dashed regions. Scale bars denote ${50}\;\mu \textrm{m}$.

Download Full Size | PDF

Table 1. Imaging setup parameters for the phase reconstructions considered in sections 4.1 to 4.4. In the case of measurements at multiple source-sample-distances, the holograms have been rescaled to a single geometric magnification. The Fresnel numbers are stated using the effective pixel size as a lateral reference length-scale.

View Table | View all tables in this article

4.1 Test structure of polysysterene microspheres

First, we demonstrate the performance of nonlinear phase retrieval using a test object of polystyrene spheres with diameter ${15}\;\mu \textrm{m}$ [25]. Holograms were acquired at four different defocus distances (see Table 1 for detailed parameters), as in standard linearized (CTF-based) phase retrieval to compensate for minima in the CTF [10,18]. The data set has been previously considered to demonstrate the advantage of iterative phase retrieval with an alternating projections (AP) algorithm over CTF [25]. Figure 1 shows a comparison of reconstructions obtained by (a) standard (unconstrained) CTF, the proposed NLTikh algorithm (subfigures (b), (c), (d)) as well as (e) the AP-result from [25]. The same parameters, $c_{\beta /\delta } = 0$ (pure phase object) and $\alpha _{\textrm{low-freq}} = 10^{-3}$, $\alpha _{\textrm{high-freq}} = 10^{-1}$, were used for all CTF- and NLTikh-reconstructions. Table 2 gives an overview over the observed computation times for the different reconstrution methods.

Table 2. Computation times for phase reconstruction of the polystyrene-microspheres data set using different algorithms. See Fig. 1(a–e) for the corresponding reconstructed images, respectively. All GPU-computations were performed in double-precision floating-point arithmetics on an NVIDIA Quadro RTX 6000 using MATLAB R2020a.

View Table | View all tables in this article

As noted already in [25], the CTF result shows significant artifacts from the large phase gradients introduced by the spheres, whereas NLTikh and AP yield visually artifact-free images of similar quality. This is confirmed by the line-scans in Fig. 1(f) where NLTikh and AP both show good agreement with the theoretically predicted phase profiles for $15 \, \mu\textrm{m}$-sized polystyrene spheres (assuming the literature value $\delta = 3.673\cdot 10^{-6}$ [45] for the refractive index of $\textrm{C}_8 \textrm{H}_8$ with density $1.05\,\textrm{g} / \textrm{cm}^{3}$). However, while the AP-reconstruction took 2500 iterations to complete [25], corresponding to minutes of computations on typical GPUs, the NLTikh algorithm requires only about 50 iterations to converge, thus completing within about two seconds on our benchmark workstation. Moreover, Fig. 1(d) demonstrates that even an NLTikh-reconstruction from a single hologram results in an acceptable image quality, apart from an expectable increase of noise compared to the four distance setting. Finally, a comparison of Figs. 1(b) and 1(c) shows that the ability of NLTikh to impose negativity of the reconstructed phases helps to effectively suppress low-frequency variations of the background and thus achieve quantitatively correct phases – without a relevant increase in computational costs.

4.2 Convergence experiments

For an exemplary assessment of convergence of the NLTikh algorithm, we consider the reconstruction setting in Fig. 1(c), that includes both constraints and nonlinearity. First, we compute a (supposedly) fully converged reference solution $\phi _{\textrm{ref}}$ by performing 1000 iterations of the proposed projected gradient descent scheme. We then compute the gradient residual $R_{\textrm{grad}}$ defined in Eq. (13) as well as the residual of the Tikhonov functional relative to $\phi _{\textrm{ref}}$:

The convergence test is carried out both for the proposed Barzilai-Borwein stepsize-rule as well as for a constant stepsize for comparison. Moreover, to motivate our warm-starting approach, we compare reconstructions initialized with a CTF-reconstruction (as proposed) to the convergence behavior if the algorithm is initialized with zeros, $\phi _0 = 0$. Results for $R_{\textrm{value}}$ and $R_{\textrm{grad}}$ are plotted in Figs. 2(a) and 2(b), respectively. The following observations can be made:

 

Fig. 2. Convergence tests for NLTikh: The plots show the residuals defined in equations (a) Eq. (15) and (b) Eq. (13), respectively, for the reconstruction setting in Fig. 1(c). Blue lines: the proposed algorithm. Red: same as (a) but without warm start, i.e. the algorithm is initialized with $\phi _0=0$ instead of a CTF-reconstruction. Yellow: same as (a) but with constant- instead of Barzilai-Borwein stepsizes $\tau _k$. By default, the iterations are stopped when $R_{\textrm{grad}}$ drops below $10^{-3}$, as marked by the black-dashed line.

Download Full Size | PDF

4.3 Golgi-Cox stained hippocampus

Next, we tested the nonlinear Tikhonov algorithm on a biological sample, notably a Golgi-Cox stained brain slice of mouse hippocampus embedded in EPON resin, which was already investigated before by CTF [46,47] and the AP-algorithm [25]. The imaging setup is detailed in Table 1. The Golgi-Cox stain is a classical stain for brain sections, which visualizes individual neurons by seemingly stochastic all-or-nothing events based on precipitation of silver along the neuron’s neurites. The metalized neurites then stand out in strong contrast with respect to the neuropil infiltrated by EPON resin, resulting in strongly varying phases. Hence, the assumption neither of a weak object, nor of a homogeneous (single material) object are strictly met.

Figure 3 shows reconstructions of the four-hologram data set obtained by (a) CTF phase retrieval, and (b) the nonlinear Tikhonov algorithm. In both cases, the same regularization parameters ($\alpha _{\textrm{low-freq}} = 10^{-3}$, $\alpha _{\textrm{high-freq}} = 10^{-1}$) and $c_{\beta /\delta }=0.1$ were used without additional constraints. Visual comparison of the reconstructed projections shows significantly less artifacts for NLTikh compared to CTF. In particular, the flawful wavy background at high spatial frequency, visible as oscillations of the CTF line-scan in Fig. 3(c), is absent. Furthermore, the NLTikh-reconstruction exhibits a much larger range of reconstructed phases compared to CTF, which is physically plausible for the high-contrast metal-stained neurons.

 

Fig. 3. Reconstruction mouse hippocampus tissue with Golgi-Cox silver-stained neurons. The images plot the phases obtained by different algorithms: (a) (unconstrained) CTF-based phase retrieval. (b) The proposed NLTikh algorithm. (c) Line-scans through the reconstructions along the blue lines. To improve comparability, the extracted phases are expressed relatively to the colorbar ranges in (a) and (b), respectively. Computation times were (a) 0.45 s and (b) 0.90 s (18 iterations). Insets in the top-left corner show zoom-images of the red-dashed regions. Scale bars are ${50}\;\mu \textrm{m}$.

Download Full Size | PDF

4.4 Sedimented silicon spheres in a capillary

For a final test case, we investigate a setting which also often encountered in practice, where a strong phase gradient at the edges of the object has to be imaged along with significantly weaker phase-variations in the interior. Examples are many: glass capillaries with tissue punches, single particles in air, trimmed specimens of material science. Here we consider silicon spheres of diameter ${2}\;\mu \textrm{m}$ in a capillary. The tomographic data set consists of 726 holograms acquired at a single source-to-sample distance (Fresnel number $F = 6.499\cdot 10^{-4}$) and equidistant incident angles covering a range of 180 degrees, see Table 1 for details. It has been recorded as part of a sedimentation experiment where the 3D-motion of the single spheres is tracked over time [48]. The considered tomogram shows the final, fully sedimented state of the sample.

We reconstruct the holograms for each incidence angle using CTF and NLTikh, imposing the $\beta /\delta$-ratio of the glass capillary ($\textrm{SiO}_2$) at a photon energy of 8 keV, $c_{\beta /\delta } = 0.0135$ according to [45], as well as negativity of the recovered phases. Regularization parameters are $\alpha _{\textrm{low-freq}} = 10^{-5}$ and $\alpha _{\textrm{low-freq}} = 10^{-1}$. Results are shown in Fig. 4.

 

Fig. 4. Reconstruction of sedimented silicon spheres in a glass capillary. The left and right column show phases obtained by CTF and NLTikh, respectively. (a), (b): Projections reconstructed for skew incidence direction of the x-rays with respect to the capillary boundaries. (c), (d): Projections for an almost parallel incidence. (e), (f) Central tomographic slices computed from the reconstructed phase projections for all 726 incident angles using filtered back-projection. Insets in the bottom-left corner show zoom-images of the red-dashed regions, respectively. Scale bars denote ${10}\;\mu \textrm{m}$.

Download Full Size | PDF

The capillary has an approximately square cross-section. Accordingly, the severeness of the phase gradients at the capillary-boundary depends on the incidence direction of the x-rays: when the latter is skew with respect to the sides of the capillary, only moderate phase-variations occur, whereas strongly varying phases result in the case of (almost) parallel incidence. Figures 4(a), (b) and 4(c), (d) shows exemplary projections reconstructed for skew and parallel settings, respectively. We observe that the CTF-reconstructions suffer from fringe-artifacts in both cases, which also distort the sphere-structures in the interior of the capillary. For NLTikh, such artifacts still occur in the parallel-incidence setting, whereas the skew reconstruction is practically artifact-free – despite the strong total phase-shifts up to $\approx {20}\;\textrm{rad}$. As revealed by the tomographic slices (computed by applying filtered back-projection to the recontructed phase images for all tomographic angles)) in Fig. 4(e), 4(f), our nonlinear NLTikh method overall achieves an improved, yet still imperfect reconstruction quality compared to CTF, yielding a clearer but still not artifact-free image of the silicon spheres’ packing in the capillary.

5. Summary and outlook

In summary, we have presented a nonlinear Tikhonov phase retrieval algorithm for near-field holography (NFH), to cope with optically non-weak objects beyond the scope of the established linear CTF model. Since the NLTikh approach comprises standard CTF phase retrieval as a limiting case, it can be regarded as a nonlinear generalization of the CTF approach. We have validated NLTikh for experimental data, including test objects and ’real samples’, recorded under relevant conditions of nanoscale holographic x-ray imaging, both in 2D and in a tomographic 3D setting. As expected by the algorithmic design, which is based on the full nonlinear model of NFH, the method is found capable to reconstruct objects with moderately strongly varying phase shifts and outperforms CTF-based methods in terms of image-quality in such settings. As a further option, NLTikh enables imposing additional support and range constraints – like the proposed constrained CTF (cCTF) algorithm for the weak object regime. This is of particular interest if hologram data is available only for a single distance, where imposing additional a priori knowledge may compensate for the lack of data richness [30,31]. We emphasize that both the ability to reconstruct optically non-weak objects and to exploit additional a priori constraints are also provided by previously considered methods, such as alternating projection schemes (AP). The distinctive feature of our NLTikh algorithm, however, lies in the demonstrated superior numerical efficiency, rendering nonlinear, constrained image-reconstruction applicable for large-scale tomographic data sets with thousands of high-resolution holograms to be processed.

Despite the merits of the presented approach, one should also note its limitations. While not being strictly limited to the weak object regime like CTF, NLTikh still relies on linearity to a certain degree: it uses a CTF-reconstruction as an initial guess and the iterations are based on derivatives of the nonlinear model, i.e. on local linear approximations. Due to non-convexity, convergence is thus only guaranteed if the true phase is sufficiently close to the linear CTF-result, i.e. nonlinearity must not be too strong. For very strong objects, as considered in section 4.4, the initial guess may be too far off for the gradient-descent scheme to find the global minimum of the Tikhonov functional. In such a setting, AP and related methods like RAAR [49] may perform significantly better than NLTikh, as they are originally designed for highly non-convex image-reconstruction problems like CDI. A second limitation of our algorithm lies in the assumption of a homogeneous (single material) object. Like the generalizations of CTF phase retrieval made in this work, lifting this assumption is conceptually simple: just optimize the Tikhonov functional in Eqs. (4), (6) or (7) for both parameters $\phi$ and $\mu$, omitting the coupling constant $c_{\beta /\delta }$. Yet, achieving accurate image reconstructions in a numerically stable and efficient manner by such an approach is expected to constitute a severe algorithmic challenge, that might be addressed in future research. Owing to the advent of energy-resolving X-ray detectors, another interesting research direction might be to extend NLTikh to settings with holograms acquired at different X-ray energies, similar to what has been demonstrated in the direct-contrast regime [50,51].