Phase retrieval for arbitrary Fresnel-like linear shift-invariant imaging systems suitable for tomography

We present a generalization of the non-iterative phase retrieval in X-ray phase contrast imaging applicable for an arbitrary linear shift-invariant (LSI) imaging system with a non-negligible amount of free space propagation (termed as Fresnel-like). Our novel approach poses no restrictions on the propagation distance between optical elements of the system. In turn, the requirements are only demanded for the transfer function of the optical elements, which should be approximable by second-order Taylor polynomials. Furthermore, we show that the method can be conveniently used as an initial guess for iterative phase retrieval, resulting in faster convergence. The proposed approach is tested on synthetic and experimentally measured holograms obtained using a Bragg magnifier microscope – a representative of Fresnel-like LSI imaging systems. Finally, the algorithm is applied to a whole micro-tomographic scan of a biological specimen of a tardigrade, revealing morphological details at the spatial resolution of 300 nm – limiting resolution of the actual imaging system.


Introduction
Phase retrieval is one of the most challenging topics in X-ray full-field phase-contrast tomographic imaging. Historically, there are two possible imaging regimes depending on the Fresnel number identified as the near-field and the far-field regimes. Both provide different scopes of application based on the required field of view (FOV) and the spatial resolution [1]. We will now restrict our considerations to the phase retrieval problem in the near-field regime. Many different approaches to solving this problem have been proposed [2, 3]. All of them can be formally divided into two groups -iterative methods and noniterative (one-step) methods.
The former methods were adapted from the phase retrieval methods in Coherent Diffractive Imaging (far-field regime), where we iteratively perform forward and backward propagation while applying constraints in the object and detector plane until the algorithm reaches the convergence criteria. The advantages of the iterative retrieval are the quantitativeness of the results and usually the best quality of reconstruction (assuming the support function is provided), which comes at the expense of computational time. This can be an issue especially for the cases of tomography or time-resolved tomography, where thousands of images need to be interpreted, so that the reconstruction procedure (which usually needs to be parametrized first) becomes a lengthy process.
The latter methods are usually based on one-step filtering in Fourier space, where the analytical formulas for the particular filter depend on the assumptions made, e.g. homogenous object, pure phase object, short propagation distances or various parameters of the optical elements. The first advantage is the speed of reconstruction, which is usually below one minute for the phase retrieval of whole tomographic set. Secondly, they can be used as a first guess for the iterative algorithms, which then can converge faster to a high-quality reconstruction [4]. Third advantage is that they are not dependent on the support function (in contrast to iterative methods), which makes them a popular method for the cases where support function is not available (e.g. region of interest tomography).
Most of the phase retrieval methods for in-line holography were developed for propagationbased phase contrast imaging [5, 6], where one typically deals with the analytical form of the Fresnel propagator and derives diverse one-step phase retrieval approaches. An excellent comparison of the available non-iterative methods can be found in [2]. Other phase-contrast imaging methods also cope with the phase problem, but using appropriately adjusted methods.
Among these, a significant group are imaging methods using linear shift-invariant (LSI) imaging systems, where the operator relating the input wavefield Ψ I N (right behind the sample) and the output wavefield Ψ OUT (at the detector) exhibits LSI properties [1]. The relationship in Fourier space can then be written as a simple multiplication by a propagator (also called transfer function) where ∼ denotes the 2D Fourier transform with respect to the spatial coordinates and f is the coordinate vector in Fourier space. Unlike the case of propagation-based phase contrast, we usually do not know the analytical formula for H(f) or the formula is very complicated. Therefore, to solve the phase problem using a one-step method, the approximations need to be introduced. One such approach is to linearize the propagator H(f) resulting in an elegant phase reconstruction formulas [7]. This requires the Fresnel free space propagation term to be negligible compared to the linear transfer function of the optical elements. Another approach assumes slowly varying transfer functions while also taking into account the free space propagation term, but it is applicable just for large generalized Fresnel numbers, i.e., short propagation distances [8].
Although both approaches proved their usefulness within the validity of their assumptions, they fail to work for the cases of higher spatial resolutions or significant free space propagation distances between the sample, optical elements and the detector. A general solution to the phase retrieval problem for the case of an arbitrary LSI imaging system was proposed in [9], with the assumptions made only about the sample. Specifically, a slowly-varying phase and a weak absorption are assumed. However, the method requires acquiring at least two images for two different states of the LSI imaging system to be able to solve the phase problem for one tomographic projection.
In this article, we present a theoretical treatment for phase retrieval based on contrast transfer function (CTF) approach (also called Born approximation [10]). It uses only single intensity measurement per tomographic projection and at the same time, it relaxes the restrictions put on the propagator of the whole setup, allowing us to treat the cases with an arbitrary amount of free space propagation. Limiting assumptions are made just for the subpropagators corresponding to the individual optical elements, which should exhibit specific smooth properties. The outline of the article is the following. Section 2 describes the principle of the conventional CTF phase retrieval approach. Section 3 introduces the novel theoretical formalism for phase retrieval, which is, in turn, tested in Section 4 for the case of Bragg Magnifier imaging setup [11] on both, synthetic and real-world experimental data. Finally, Section 5 presents the application of the proposed method for phase retrieval of the complete tomographic set of images resulting in successful 3D electron density mapping of small animal Tardigrade ("water bear"). This is followed by the Conclusion, which summarizes the potential of the proposed method, especially when applied to the case of X-ray imaging with Bragg Magnifier.

Overview of the conventional CTF
Imagine typical phase contrast imaging experiment as drawn in Fig. 1. The source of the X-rays is producing coherent monochromatic radiation which hits the sample S, propagates through an optical system and reaches the detector D. After normalization of the diffractogram by flat-field correction, the exit wave (the wavefield right behind the sample) equals the transmission function of the object T(x 1 , x 2 ) given as: where x 1 , x 2 are the cartesian coordinates in the plane perpendicular to the optical axis z, µ(x 1 , x 2 ) is the projected absorption of the sample and φ(x 1 , x 2 ) is the total phase shift caused by the presence of the sample. Let us denote Ψ(x 1 , x 2 ) the wavefield at the detector plane. Then, if the Fig. 1. The sketch of the typical X-ray phase contrast imaging experiment. The produced X-rays pass through the sample S modifying the wavefield, which is further propagated through optical elements and reaches the detector D.
used optical system is LSI, we get where ⊗ is the convolution, symbol ∼ denotes the Fourier transform, F −1 is the inverse Fourier transform, h(x 1 , x 2 ) is the point spread function of the optical system, f 1 , f 2 are the reciprocal coordinates dual to x 1 and x 2 and H( f 1 , f 2 ) = h(x 1 , x 2 ) is the transfer function (propagator) of the optical system. The measurable quantity is the optical intensity given as For simplicity, let us first deal with the 1D case, where we will use the real space coordinate x and its dual coordinate f . Taking the Fourier transform of Eq. (3) yields All the integrals' boundaries in this section are (−∞, ∞). Similarly from Eq. (3) we get Now we take the Fourier transform of Eq. (4) and using Eq. (5) and Eq. (6) we obtain where we denoted the quantity in the curly brackets as I H (x − y, f ). This is the crucial expression that we will further analyze.
If we take the propagator H( f ) to be the free space Fresnel propagator (no optical elements) Because of the delta function, Eq. (8) can be integrated over y to get or alternatively by shifting the x coordinate we arrive at This expression was first derived by Guigay [12] and opened the possibilities to perform noniterative phase retrieval for special cases [6]. Specifically, if we assume pure phase object (µ(x 1 , x 2 ) = 0 in Eq. (2)) and further assume slowly varying phase, then we arrive at the well-known expression for contrast transfer function From this, for f 0 the delta function disappears and one can extract the phase of the transmission function as where α( f ) is the regularization parameter. The value of the DC component φ( f = 0) is lost, which means, that our resulting phase φ(x) will not be scaled properly. This is not a problem when one is interested only in contrast between different parts of the sample, which is often the case in biological imaging. Nevertheless, this can be cured to yield also a quantitative reconstruction, which will be described in Chapter 4.

Description of our new adapted method
Our goal now is to generalize the idea for a larger group of LSI imaging systems. First of all, most of these systems contain parts where the wavefield is still subjected to free space propagation. For example, if the system contains one optical element, the full propagator H( f ) can be expressed as a multiplication of free space propagation from the sample to the optical element, the propagator of the element itself and the free space propagation from the element to the detector. In the light of this, we can express the full propagator H( f ) as where z e f is the total effective propagation distance between all the optical elements and A( f ) is the term perturbing the Fresnel propagator. Typically, A( f ) represents the transfer function of all the involved optical elements. In general, A( f ) is a complex function which can be decomposed to real functions as Here we propose the following approximations: The physical meaning of the first statement is that we assume our imaging system will be equally absorptive for all the relevant spatial frequencies. This can be true for many applications, where the angular acceptance of the system is larger than the highest desired spatial frequency. The second expression is the Taylor expansion with respect to the zeroth frequency of the system's phase shift truncated at the second order, which is a good approximation for many optical elements with smooth phase-shifting properties. In the light of these assumptions, A( f ) in Eq. (13) can be viewed as a slowly varying envelope perturbing the fast oscillating Fresnel propagator. Applying Eqs. (13) - (16) to the definition of I H we obtain where we used A(0) ≡ A 0 , ϕ(0) ≡ ϕ 0 and Eq. (8) reduces to Fresnel-like form (compare with Eq. ) In the spirit of the previous chapter, if we assume pure phase object and slowly varying phase, i.e., we arrive at the final form of the 1D version of contrast transfer function for an arbitrary LSI imaging system, whose optical elements disturb the wavefield smoothly Comparing with Eq. (11) we immediately see the role of individual parameters. First of all, ϕ 0 itself is not present at all, since the measured intensity is independent from the constant phase term in the transfer function. The scaling parameter A 0 is not important here, because it cancels out after the flat-field correction. The slope of the phase ϕ 0 is present within the complex exponential, which in its periodic nature provides an envelope for the rapidly oscillating term in the curly brackets. Finally, ϕ 0 slightly modifies the effective propagation distance and as a consequence, moves the positions of the zero-crossings of the sine function.
Again, for f 0 we obtain the 1D phase retrieval formula Going back to 2D phase retrieval problem starts with replacing the real and reciprocal coordinates in Eq. (8) with 2D vector coordinates and upgrading the approximations in Eqs. (15) and (16) to where f = ( f 1 , f 2 ) is the reciprocal space vector, ϕ 0k is the first derivative of phase with respect to k-th reciprocal coordinate at the (0, 0) point and ϕ 0kl is the second derivative with respect to k-th and l-th coordinate at the same point. Analogously, we can derive the expression for the contrast transfer function as where z (k) e f is the k-th effective propagation distance (in general, it may differ for individual dimensions). This was achieved by the following approximation (27) Finally the phase retrieval formula for the 2D case is .
(28) This fast one-step reconstruction approach can be used separately for each tomographic projection and followed by one of the inverse Radon transform approaches to get approximative 3D electron density map of our object. Alternatively, one can improve the CTF reconstruction by feeding it into the iterative phase-retrieval algorithm as the one we recently proposed [13]. Using CTF reconstruction as the initial guess as well as extracting the object's support out of it can then significantly help to reduce the number of iterations required for the convergence of the iterative algorithm. The scheme of the workflow is sketched in Fig. 2.

Test of the method
To perform the numerical tests of the proposed method, we used the optical setup of quasi-channel cut Germanium-based Bragg magnifier [11]. The system consists of 4 crystals magnifying their corresponding incident wavefields with magnification factors M 1 , M 2 , M 3 , M 4 , respectively and modifying the incident wavefields by their transfer functions E 1 , E 2 , E 3 , E 4 , respectively. As shown in [14], the propagator of the whole setup can be expressed as the division by the sine function, which leads to the noise amplification around its zero-crossings, hence lowering the contrast for particular frequencies. For the case of conventional CTF, many different strategies of curing this problem have been proposed by the proper choice of α( f 1 , f 2 ) function. The most common one is to choose two constant regularization values -one for low and one for high frequencies, with the smooth transition between them [15]. A more sophisticated way was also proposed, where the regularization term is chosen to be a nonzero constant just around the zero-crossings leading to improved quantitativeness of the reconstruction [16,17]. This method is also referred to as the quasi-particle (QP) method, and it is the one used in the following tests.
As already mentioned at the end of the Chapter 2, CTF reconstructions are not quantitative since we lose information about the DC component φ(0, 0). However, to use our modified CTF approach as an initial guess for the iterative algorithm, it needs to be reasonably close to the final solution of the phase φ(x 1 , x 2 ) also in a quantitative manner. Otherwise, it will not speed up the iterative refinement. The way to improve it is to perform a least squares fit for the unknown parameter φ(0, 0), which will minimize the distance between the calculated intensity pattern I(x 1 , x 2 ) (from the current φ(x 1 , x 2 )) and the measured intensity pattern I ex p (x 1 , x 2 ). From the computational point of view, this takes only a few seconds on a PC using 1 CPU and usually, it is enough to perform this fit for only one 2D projection and apply the φ(0, 0) value to the whole tomographic dataset. Therefore, it does not slow down the whole process. Fig. 3 shows the results of the application of the proposed method on the simulated data. A phantom image (marked as 'Original phase map') is first computationally propagated to the detector through the Bragg Magnifier setup ('Simulated hologram') with the effective propagation distances z (1) e f = 513 mm and z (2) e f = 758 mm, photon energy E = 10.7 keV and the total magnification factor due to the crystals M = M 1 M 2 = M 3 M 4 = 180 for both directions. Then the CTF reconstruction is carried out (Fig. 3) yielding reasonable reconstruction even for such a complicated object, bearing in mind that the result is obtained using just a single distance detector data. The lower left image also shows the horizontal slice through the middle of the image showing the quantitativeness of the result. After performing additional ∼ 300 iterations of the iterative algorithm, the phase reconstruction has converged to the object practically identical to the original one ('CTF+iterative reconstruction'). On the other hand, when using the iterative algorithm only, i.e., starting with a random object, using the whole image as an initial support and employing the shrink-wrap (SW) algorithm as described in [13], the convergence is slower. One usually needs thousands of iterations to achieve the same result as accomplished by CTF as an initial guess, where the convergence occurs in hundreds of iterations (error plot).
Subsequently, we applied the method to the experimentally measured holograms of Tardigrade ('water bear'), under the same experimental parameters as used for the simulation of the phantom image above. Again, the measured hologram, the CTF reconstruction, and the iterative refinement are shown in Fig. 4 (revisited from [13]). The lower left image shows the comparison of the horizontal middle cut through the images obtained by the CTF reconstruction alone and the iterative refinement. The lower right image demonstrates the superiority of our new proposed method compared to the SW iterative algorithm.

Application to small animal tomography
To perform a proof of principle of our new method's applicability to microtomography, we picked Tardigrade Echiniscus, a model organism for microtomography imaging. Tardigrades are attractive test target for X-ray imaging for several good reasons. Species Hypsibius dujardini is one of the most frequently studied model organism of developmental biology in general [18]. Tardigrades are also transparent, which allows correlating data, which allows correlating data obtained from complementary imaging modalities, e.g. well-known optical techniques. The size of the species is in hundreds of micrometer range, thus suitable for microtomography. There are several forms worth studying: mature living forms, tun states, and embryonal forms. For X-ray imaging is the most convenient anhydrobiotic tun state -the state in which Tardigrade survives extreme external conditions, including both low and high temperatures as well as high levels of ionizing radiation [19]. In anhydrobiotic states, the amount of water is dramatically reduced and DNA is protected using several intensely studied molecular mechanisms [20], such as disaccharide trehalose substituting molecular structure of water clusters and dedicated proteins studied for possible use in human medicine [21]. In tun state, the cell nuclei are condensed, which makes them visible in hard X-ray without contrasting techniques. The number of cells and their development from an embryo into an adult organism is thoroughly studied and described in details for most popular Tardigrades.
Our experimental samples consisted of Tardigrade of genus Echiniscus, dried to the residual water content below 5%. After applying the phase retrieval procedure described above to every 2D projection and performing filtered back-projection, we obtained 3D electron density model of our sample. A central 2D slice and 3D rendering view are shown in Fig. 5. We have marked the visible internal structures of the sample in the left part of the figure. The right part shows the 3D electron density distribution. The resolution we achieved by this method is 300 nm determined by the Fourier shell correlation method. [22].

Conclusion
We have designed an efficient single-measurement phase retrieval method applicable for X-ray phase contrast tomography using linear shift-invariant optical elements with a smooth transfer function. The power of the method is that it lifts restrictions on the free space propagation distances between the sample, optical elements and the detector as well as it improves the speed of the phase reconstruction while preserving the high quality of the reconstruction and quantitativeness of the results. The latter is achieved by using the combination of the non-iterative and iterative approaches, namely using our formalism for the one-step phase reconstruction based on the modified contrast transfer function and the previously published iterative algorithm.
We have demonstrated the feasibility of the method for the case of the Bragg Magnifier Microscope -an imaging modality that uses asymmetric Bragg reflection to magnify the Xray beam. The tests were performed first on noiseless synthetic phantom data and then on experimentally measured data of the biological sample -microorganism Tardigrade. Both proved the superiority of the combined method compared to the purely iterative algorithm.
Based on the proposed method, we have finally performed the reconstruction of all the tomographic projections and obtained the 3D electron density map of the Tardigrade, where the essential morphological features from other imaging modalities are visible. We achieved full 3D resolution of 300 nm, which is the theoretical limit of the given experimental setup. This algorithm extends the applicability of the Bragg Magnifier Microscope for 3D in-vivo dose-efficient tomographic imaging of biological objects with submicron resolution.