Single-distance phase retrieval algorithm for Bragg m agnifier microscope

: We present an improved, single-distance phase retrieval algorithm applicable for holographic X-ray imaging of biological objects for an in-line germanium Bragg Magniﬁer Microscope (BMM). The proposed algorithm takes advantage of a modiﬁed shrink-wrap algorithm for phase objects, robust unwrapping algorithm as well as other reasonable constraints applied to the waveﬁeld at the object and the images as determined by fourier analysis,” Phys. Status Solidi A 204 , 2746–2752 (2007).


Introduction
Asymmetrically cut crystal used for magnification of an X-ray beam and its application for X-ray imaging was firstly proposed and demonstrated by Kohra [1] for one dimensional magnification. Later on, the two dimensional magnification was demonstrated by Boettinger [2] where two crystals arranged in so-called σ-π configuration with diffraction planes oriented perpendicularly to each other have been used, leading to image magnification in both, the horizontal and the vertical direction. The image formation, based on the wave optics approach for two crossed crystals, have been formulated by Spal [3]. Later on, the algorithm for the inverse problem has been proposed by Modregger [4]. Since then, various successful improvements have been made and application for high resolution X-ray imaging has been demonstrated [5][6][7]. The two dimensional magnification has also been achieved by single monolithic silicon crystal using two successive diffractions by Korytár [8]. However, the major limitations of those solutions, such a narrow operation energy range and the limited angular acceptance (numerical aperture), limiting the spatial resolution to > 0.5µm of the silicon devices, did not contribute to the practicality of the Bragg Magnifier Microscope (BMM). Recently, the problem of limited operational energy has been overcome by a variable magnification system composed of two crossed silicon crystals [9], where one can vary the effective asymmetry angle over the entire range covering maximum magnification to maximum compression for fixed energy. In principle, variation of the energy for fixed magnification over broad range is possible as well.
The problem with narrow numerical aperture was significantly improved by introducing germanium crystals [10] instead of silicon improving the numerical aperture by almost a factor of two. Further combination of germanium based BMM with single photon counting detector (Si Medipix hexa detector) improved strongly the detection quantum efficiency and the sensitivity of the system [11]. Narrow angular diameters of new X-ray sources, especially X-ray Free Electron Lasers with their spectral bandwidths which are matching or are smaller than the angular-spectral acceptance of BMM, will increase the throughput of this microscope to its maximum limited only by the peak reflectivity of the used crystal. The efficiency of BMM in such conditions will reach 50%-80% depending on the asymmetry used. This may be very attractive for high spatial resolution imaging of processes reaching femtosecond time scales.
The detection sensitivity is especially important, when the imaging is performed in holographic mode, where the sample is placed far enough in front of BMM and the spatial resolution is determined by the visibility of outermost finest Fresnel oscillations. Because recorded images are in-line holograms, they contain a large number of Fresnel oscillations and have to be numerically decoded in order to obtain the so-called exit wave and the complex transmission function of the sample [12]. We have recently reported iterative phase retrieval algorithm which uses wavelet filtering [13] considering BMM as an linear shift invariant imaging system. However, the performance of the algorithm was often failing in the cases of more complex and optically dense objects, leading to a kind of halo effect at the object boundaries reducing the final resolution. In this work we present a new algorithm for BMM holograms. Its performance is demonstrated on experimentally measured holograms of well-defined samples and more complex biological objects. The complex transmission function is retrieved from only single distance measurement and the spatial resolution is reaching 300 nm.
The work is structured in a following way: Section 2 describes the image formation and the design of the algorithm, in Section 3 we show the performance of the algorithm on the measured data and the last section contains the summary of this work and further prospects towards the in-vivo biological X-ray imaging and ultra-fast X-ray imaging.

Theory and algorithm
First, let us define the basic variables and terminology. The plane perpendicular to the beam propagation located right behind the sample of interest will be referred to as the object plane. Similarly, we define the plane perpendicular to the beam direction at the location of the detector as the detector plane. Thanks to the in-line character of BMM, these two planes are parallel to each other. Let us denote Ψ O (x, y) the wavefield at the object plane and Ψ D (x, y) the wavefield at the detector plane. All these planes and wavefields are illustrated in Fig. 1. From the mathematical point of view, Ψ O (x, y) and Ψ D (x, y) are complex functions of two real variables, x and y, which represent the orthogonal cartesian coordinates describing the position at the corresponding planes. The task is to determine Ψ O (x, y) using just the amplitude |Ψ D (x, y)|, which is known from the experiment, since it is related to the measured intensity at the detector plane I D (x, y) by |Ψ D (x, y)| ∝ I D (x, y). Since under some conditions BMM can be considered as a linear shift invariant system [13], the relationship between Ψ O and Ψ D is given by the convolution theorem [12] Ψ where F and F −1 denote the direct and inverse two dimensional Fourier transform, respectively, P B M M is the propagator for Bragg Magnifier and k x , k y are reciprocal coordinates corresponding to x and y, respectively. We will refer to Eq. (1) as forward propagation and, similarly, we define an inverse relationship which will be referred as backward propagation.
The propagator for BMM, P B M M , is determined according to [13] as where P S1 , P 12 , P 23 , P 34 and P 4D are respectively the free space propagators for propagation between the sample and the first crystal, the first and the second crystal, the second and the third crystal, the third and the forth crystal and the forth crystal and detector. Similarly, E i , i ∈ {1, 2, 3, 4} is the crystal transfer function, which describes the modification of the wavefield due to the presence of the i-th crystals. All the crystal transfer functions have been calculated from the first principles using the two-beam dynamical theory of X-ray diffraction formulated by Authier [14] and Huang [15]. As a summarization, the propagator P B M M is dependent on the geometry of the BMM setup (distances between the sample, crystals and detector as well as the crystal rotations and crystal asymmetry angles), characteristics of the source (photon energy) and the material properties of the crystals for the given energy (electrical susceptibilities). The details about the values of these parameters are given either in [11] for the general parameters or in the next section for the changeable parameters. If we can recover Ψ O (x, y) from the measured intensity |Ψ D (x, y)| 2 , then we will have access to the object transmission function T (x, y) using the relation where Ψ S (x, y) is the wavefield produced by the source of the X-rays right before it hits the sample (see the Fig. 1). Equation (4) is obtained using the projection approximation [12]. In all our cases, the plane wave illumination is used. By normalising the holograms by the holograms without the sample (flat field correction), Ψ S (x, y) will have unit amplitude. Furthermore, using the freedom in choosing the phase, we take the phase of Ψ S (x, y) to be zero, implying Ψ S (x, y) = 1 everywhere and consequently, Ψ O (x, y) = T (x, y). Therefore, reconstructing Ψ O (x, y) means obtaining T (x, y), which in turn means obtaining the absorption and phase shifting properties of the sample, e.g. the projected complex index of refraction of the sample. Let us discuss the phase retrieval algorithm. For simplicity, we omit indicating the dependency of Ψ O and Ψ D on x and y. Let us denote Ψ (i) O and Ψ (i) D the i-the estimate of Ψ O and Ψ D , respectively. We denote the corresponding phases of these wavefields as ϕ (i) 0 and ϕ (i) D . The initialization steps include: i1 Preprocess the normalized holograms by applying unit padding and up-sampling (detailed description is given below).
i2 Calculate the propagator P B M M according to Eq. (3).
i3 Define initial value of the support function S. It is a binary function, which assigns 1 to the pixels with a high probability of the sample occurence and 0 otherwise. The algorithm usually converges for any choice of initial support, therefore, it does not require any additional a priori information about the sample. However, a better estimate for the support can speed up the convergence considerably. One way to define the better initial support is to use the measured hologram from which we can, based on visual inspection, extract a rectangle containing the sample or define the area more accurately by drawing a free hand curve.
i4 Generate the zeroth estimate of the unknown phases ϕ (0) D with uniformly distributed random numbers from −π, π) and consequently, the zeroth estimate of Ψ D as i5 Perform a backward propagation using Eq.
(2) to obtain the zeroth estimate Ψ (0) O . Subsequently, the i-th iteration of the algorithm is given by: 1 Perform phase unwrapping on ϕ (i−1) 0 using the robust unwrapping algorithm [16].
2 Perform the support update using the shrink-wrap algorithm for phase objects (see below).
It is a phase-contrast analogy to the original idea of the shrink-wrap algorithm [17] used in coherent diffractive imaging phase retrieval methods. This step is not executed in every iteration, just after each ∼ 20 iterations. Based on our experience, this prevents rapid changes in the calculated phases, which potentially ends in the local minima. Our observation is in agreement with the proposed frequency of performing the original shrink-wrap algorithm by Marchesini et. al. [17].
3 Restrict the phases ϕ (i−1) 0 to nonpositive values inside the support region and to zero values outside the support region, that is This is due to the fact, that the real part of the complex index of refraction, in case of the X-rays, is for almost all substances less than 1. Thus, the wavelength inside the sample becomes longer, therefore shifting the phases to the lower values in relation to the phases of the unscattered radiation, which are zero by convention.
4 Restrict the amplitudes Ψ (i−1) O to values not more than 1 inside the support region and to exactly one outside the support region, that is This is due to the fact, that the imaginary part of the complex index of refraction is always non-negative, since there is always an absorption happening.
(1) to obtain the i-th estimate Ψ (i) D .
6 Replace the amplitudes of Ψ (i) D with the experimentally measured amplitudes √ I D while preserving the phases, that is 7 Perform backward propagation of Ψ (i) D using Eq.
(2) to obtain the i-th estimate Ψ (i) O . After the last iteration of the algorithm, one has to perform one more final phase unwrapping according to step 1 of the main algorithm to get the desired result.
The modified shrink-wrap algorithm in step 2 consists of these substeps: sw1 Perform smoothening on ϕ (i−1) O as a convolution with a two-dimensional symmetric gaussian function G(0, σ) with zero mean vector [0, 0] and standard deviations σ to get smoothed phase ϕ (i−1) sw2 Find the minimum value ϕ (i−1) min of the smoothed phase, e.g.
sw3 Define new support function S as where t ∈ 0; 1 is a threshold parameter, which controls the tightness of the support. Usual value of parameter t is around t = 0.1 and its interpretation is the following -support function will then spread over the pixels with the phase value lower than 10 % of the minimum phase (do not forget that phases are nonpositive numbers in the current context). Typical value for σ is around σ = 5∆x, where ∆x is the linear pixel size in the object plane. The choice of these two parameters depend on the object and can significantly influence the convergence of the algorithm. A computer implementation of the algorithm requires discretization of x and y values within the object and detector plane. Analogously to zero padding used in far field imaging reconstructions, in step i1 we apply so called unit padding. If the experimentally measured normalized hologram is of size K × L pixels, we extend the size to M × M, where M is calculated as the nearest power of 2 bigger than max{K, L}. The original normalized hologram is then centered within the extended window and the values of the extra pixels around are set to unity. The reason for this is to avoid numerical boundary artefacts due to using the discrete Fourier transform and at the same time optimizing the speed of the fast Fourier transform algorithm, which works best for square windows of a size, which is a power of 2.
At some point, we found out, that the achieved resolution of the reconstruction is limited by the pixel size of the detector. In order to increase the resolution, we opted for up-sampling of the holograms [18]. Up-sampling of order m is defined as follows: Take the padded difractogram of size M × M and divide each pixel into a regular grid of size m × m, (m ∈ N) creating m 2 subpixels within each pixel. Finally, associate each of these subpixels with the same intensity value that corresponded to the original pixel. Formally, let I M (i, j) be the measured intensity at the detector plane corresponding to the pixel in i-th row and j-th column of the padded difractogram (i, j ∈ {1, . . . , M }) and let I m M (p, q) denote the up-sampled intensity (p, q ∈ {1, . . . , mM }) given by I m M (p, q) = I M (p − 1)\m + 1, (q − 1)\m + 1 , where the symbol \ denotes the integer division. This will enlarge the hologram to its final size mM × mM which is then used for reconstruction. For experimental data we observed, that it was sufficient to use m = 2 or m = 3, since higher values have not provided further improvement in terms of spatial resolution. The error of the reconstruction, calculated in each iteration was defined as the euclidean norm of the difference between two consecutive iterations of Ψ O : The algorithm is said to converge, if the error of the reconstruction is smaller than the predefined threshold value.

Results and discussion
First, we applied the proposed phase retrieval algorithm to a known object. Figure 2   branch (UK). The used energy was 10.7 keV and the magnification was determined to be 163.25 in the horizontal and 174.35 in the vertical direction. The difference in magnification in orthogonal directions was caused by the fact, that the used crystals did not have exactly the same asymmetry. As detector, a single photon counting Timepix "Hexa" module was used, consisting of 6 readout chips in a 2 × 3 configuration, bump-bonded to a monolithic, 500 µm thick silicon sensor [19]. The upper right part of the figure shows the reconstruction according to the new proposed algorithm. The hologram has been clearly focused in a very nice manner. Let us also note, the the blur artifacts of our previous algorithm [13] evident at the edges (see figure 2, lower left and lower right) has been eliminated by the presented algorithm.
In order to evaluate the quantitative accuracy of the reconstruction, we plotted the profile of the phase along the horizontal axis, shown at the right part of Fig. 2. It is compared with the theoretical phase profile of the PS sphere calculated using the projection approximation. For the new algorithm, the agreement is very good, a slight difference at the center is probably caused by having not exactly spherical shape of the real sphere, i.e. the real sphere is ellipsoidally flattened in the direction of the beam. We have plotted the errorbar of the theoretical profile at the center of the sphere, which has been calculated based on the standard deviation of real sphere radius provided by manufacturer. We see, that our reconstruction is within the errorbar.
Another known test sample used for estimating the performance and the spatial resolution of the algorithm was X-radia 50-30-2. It is a 180 nm high golden structure on Si 3 N 4 membrane, which consists of Siemens star pattern in the middle with diameter of 30 µm and the surrounding objects, like horizontal and vertical stripes of various linespacings. The smallest feature size is 50 nm. For the given energy, we can consider it as almost pure phase object. The hologram and the reconstruction are shown in Fig. 3. The beamline and the experimental parameters were the same as in the previous case of the PS spheres. The resolution of the reconstruction can be obtained directly using the known shape of the horizontal and vertical strips located at the upper and right part of the object. Using this direct method, the resolution is determined to be ∼ 300 nm. The lower part of the Fig. 3 shows the zoom of the reconstructed image corresponding to the blue rectangular area. Moreover, it provides the comparison between the reconstruction obtained using m = 1 (no up-sampling of the hologram) and m = 2 (up-sampling of order 2). The reconstruction using up-sampling yields more accurate positioning of the narrower strips, therefore pushing down the resolution to a higher value. The reason is that the pixel size without up-sampling will sooner be equal to the linespacing of the strip, what leads to principal limitation of the line visibility.
In order to apply our new phase retrieval method for more complex biological objects, we performed another series of measurements using the same beamline at Diamond Light Source. The photon energy was 10.76 keV and the measured magnification was 185 in both directions (we achieved almost the same magnification in both directions after correcting the asymmetry angles by repolishing). The sample of interest was chosen to be Tardigrade (water bear), a model organism, which is of high interest in developmental biology [20]. Figure 4 shows the measured hologram and the corresponding phase retrieval. After finding the right set of parameters, the algorithm converged after 5000 iterations. From the point of view of computation time, the reconstruction took about one hour when implemented on CPU (Intel Xeon E5-2650, 2.80GHz), but just around 2 minutes using GPU implementation (NVIDIA Titan Black). The robust phase unwrapping algorithm [16] applied in the step 1 of the main algorithm, was not needed for the previous cases of PS spheres and X-radia, since the maximal phase shift was lower than π in absolute value. But in the case of Tardigrade, it played an important role and allowed us to process the absolute phases in each iteration. This is necessary for the shrink-wrap algorithm in step 2 and the phase constraint in step 3. In order to determine the spatial resolution of the reconstructed image, we applied the Fourier power spectral method [21]. It comprises of taking an arbitrary line profile along the reconstructed image and applying the Fourier transform to get its spectral power. The lower right part of Fig. 4 shows the spectral power |P| 2 corresponding to the line profile marked in the reconstructed image by blue line. The noise level µ N is determined as the mean spectral power of higher frequencies considered as noise (in our case we chose this to be more than 4µm −1 ). The threshold spectral power is then taken to be 2µ N . Finally, we find the minimum frequency not satisfying the threshold limit as well as the maximum frequency still satisfying the limit, which defines the resolution interval (3.1 ± 1.0) µm −1 corresponding to the real space spatial resolution (0.3 ± 0.1) µm. The resolution value is practically independent of the blue line position.

Conclusion
In conclusion, we have revisited and improved a technique for interpreting the holograms obtained in X-ray imaging of biological objects using Bragg Magnifier Microscope. This approach consists of preprocessing the data by unit padding and up-sampling followed by the application of the iterative phase retrieval algorithm based on the modified shrink-wrap algorithm for phase objects, the robust phase unwrapping and other constraints. The successful demonstration was performed on the holograms of the PS sphere, X-radia 50-30-2 and biological object Tardigrade. To determine the spatial resolution of the reconstructions, we applied simple direct method for test pattern X-radia and Fourier method for biological objects in the case of Tardigrade. In both cases it was proven to be (0.3 ± 0.1) µm suggesting the Bragg Magnifier Microscope a useful tool for in-vivo imaging of biological samples at the submicron resolution range. Another application we can foresee is the femtosecond imaging of dynamical processes such as laser induced plasma generation or shockwave propagation.