An In-House Cone-Beam Tomographic Reconstruction Package for Laboratory X-ray Phase-Contrast Imaging

: Phase-contrast, and in general, multi-modal, X-ray micro-tomography is proven to be very useful for low-density, low-attention samples enabling much better contrast than its attenuation-based pendant. Therefore, it is increasingly applied in bio- and life sciences primarily dealing with such samples. Although there is a plethora of literature regarding phase-retrieval algorithms, access to implementations of those algorithms is relatively limited and very few packages combining phase-retrieval methods with the full tomographic reconstruction pipeline are available. This is especially the case for laboratory-based phase-contrast imaging typically featuring cone-beam geometry. We present here an in-house cone-beam tomographic reconstruction package for laboratory X-ray phase-contrast imaging. It covers different phase-contrast techniques and phase retrieval methods. The paper explains their implementation and integration in the ﬁltered back projection chain. Their functionality and efﬁciency will be demonstrated through applications on a few dedicated samples.


Introduction
Over the last two decades, the capabilities and potentials of X-ray computed tomography (CT) and micro-CT has been greatly improved by applying phase-contrast methods, or in general, multi-contrast imaging modalities (phase-contrast, dark-field imaging). This is especially beneficial for biomedical and life science applications, which typically involve low-attenuation-contrast samples that produce very low contrast-to-noise ratio (CNR) images and inferior quality for conventional, attenuation-based CT. In contrast to conventional X-ray CT, phase-contrast techniques rely on detecting small-angle refraction of the X-ray waves traversing the sample and the related phase shifts imparted by the sample [1][2][3]. This can be two to three orders of magnitude more sensitive than attenuation contrast in case of soft tissues and bio-samples in the hard X-ray regime. In this way, they extend the capabilities of X-ray imaging to weakly or non-attenuating samples or to those details that lack sufficient visibility in conventional CT. Another interesting modality shedding light on nano-scale structures is X-ray dark field imaging. It relies on ultra-small angle scattering of X-rays in the sample by micro-and nano-scale, typically sub-resolution, structures [4][5][6]. Reviews and comparison of X-ray phase-contrast imaging techniques are given in Refs. [3,[7][8][9][10][11]. Phase-contrast methods were first introduced on synchrotron sources but they found their way relatively fast to laboratory sources. This is due to the developments and availability in micro-and nano-focus X-ray sources and X-ray optical elements such as diffraction gratings. The main techniques used for laboratory X-ray phase-contrast imaging (PCI) include grating-based PCI [10,12,13], propagation-based PCI [1,2,14], speckle-based PCI [15] and edge illumination [16].
Although there is a plethora of literature regarding phase-retrieval algorithms, especially for propagation-based PCI, access to implementations of those algorithms is relatively limited [17][18][19][20][21][22]. Some of the packages work following the conventional pipeline and start with a phase retrieval step on the 2D projection images, which are then back projected using mostly filtered back projection (FBP) or iterative reconstruction methods [18,19,22]. It was also shown by Ruhlandt and Salditt [23] that phase retrieval and reconstruction are interchangeable also when considering approaches other than FBP, such as the Algebraic Reconstruction Technique (ART). Another work compared the conventional pipeline with 3D post-reconstruction phase retrieval, which has definitely benefits when projections images have to be stitched to achieve larger field-of-views (FOVs) at high resolution [20].
Only few packages wrap the phase-retrieval module with other codes for the full 3D tomographic reconstruction pipeline [18,19,21,22], which are available in open source form to the community. These, except for [21], which is a toolbox coupled to the ASTRA package [24], are developed exclusively for imaging at synchrotrons in parallel beam geometry. However, for phase-contrast CT using laboratory-based setup with cone-beam geometry, there is a real lack of open source reconstruction packages. Note that much more open source packages are available for conventional, attenuation-based CT reconstruction. The review of these is beyond the scope of the present paper.
Here, we introduce an in-house cone-beam PCI reconstruction package with several options for phase-contrast and dark-field micro-CT reconstruction on laboratory sources. The package enables CT reconstruction for grating-based and propagation-based PCI in cone beam geometries and features a graphical user interface. In the following, we briefly describe which methods and how are implemented in the package. After that, we will give representative examples demonstrating how the package works (see Figures 1-3). Note that the goal of the paper is not to perform a benchmark of our package against others. That would be a major effort requiring among others the definition of a consistent and suitable data set for comparison and that is clearly beyond the scope of this paper.  Figure 4c,d). The insets show the gray value histograms for the different images (black-linear scale, gray-log scale is shown to be able to better discern the differences). It is clear from those that the Paganin method strongly increases the CNR of the images as the histograms are becoming much more structured, showing several peaks compared to the attenuation case. A ν value of 10.0 was used in Equation (6) to obtain the images. The CNR values for the different methods and samples are summarized in Table 1.   Table 1).

Materials and Methods
A detailed description of our cone-beam CT reconstruction package for attenuation imaging can be found elsewhere [25]. It describes the technical implementation and GPUbased acceleration of the FDK algorithm for cone-beam reconstruction [26] together with the full data processing workflow from preprocessing, artefact corrections, center of rotation correction to reconstruction. Furthermore, we evaluated the performance of the basic FDK backprojection chain for typical, large micro-CT data sets benchmarking for different computational system architectures [25]. We describe here only the methods implemented in our phase-contrast CT reconstruction modules and how they are interface to the rest. The description is divided according to the two major PCI techniques that are implemented in the package.

Propagation-Based PCI
A simple formulation for the normalized, propagated intensity at the detector can be derived based on the so-called transport-of-intensity (TIE) equation [2,3,27]: where z is the propagation distance, k = 2π/λ is the wave-number and φ(r) is the phase shift inflicted by the sample. The propagation distance in Equation (1) has to be replaced by z e f f = z/M for cone-beam geometry and r, representing the detection plane coordinates perpendicular to the optical axis, has to be also scaled by M, the magnification, according to the Fresnel scaling theorem [2]. The TIE gives satisfactory results in the so-called direct contrast regime, i.e., for moderate propagation distances and for most practical cases for laboratory PCI. The simplest form of TIE-based phase retrieval is obtained by inversion of Equation (1) via Fourier transform (F ) as k is the Fourier-space equivalent of r. Equation (2) is the phase retrieval step of the original Bronnikov algorithm [28], which was integrated together with the ramp filter into the tomographic reconstruction. As explained above, the phase retrieval can be also performed prior to reconstruction on the projection images as is done mostly in practice, which is also implemented in our reconstruction package. The contact image I(r, 0) is experimentally challenging and not practical to measure. For pure phase objects or negligible attenuation it represents a normalization by an open-beam (no object) image, I o , instead of I(r, 0) in Equation (2). Finite object attenuation is typically taken into account adding an empirical regularization parameter, α, to the denominator of Equation (2) also to compensate the singularity at zero frequency (k = 0) and stabilize the reconstruction. This leads to: Equation (3) is the modified Bronnikov algorithm (MBA) [29]. De Witte et al. [30] proposed an alternative approach the so-called Bronnikov-Aided-Correction (BAC), which is a two-step algorithm. In the first step, an approximation for the phase distribution, φ(r) is obtained by applying Equation (3). Then the contact image is estimated by expressing it with the help of Equation (1) as where z/k was replaced by a control parameter γ that is typically determined empirically together with α [2]. A common strategy for determining α is to examine in φ(r) obtained in the first step profiles through edges and material boundaries and decrease α starting from larger values until edge effects are eliminated. Then adjust the strength of the phase correction in the second step increasing γ starting from small values and its final value is determined by visual inspection. We also apply this strategy with our reconstruction package and enable the user to choose the α and γ parameters (see an illustration of it in Figures 2 and 4).Then the contact image, I(r, 0), is used further in the conventional attenuation CT reconstruction pipeline.
Besides the BAC algorithm, we implemented in our package the likely most often applied single-distance phase retrieval algorithm by Paganin et al. [31]. It is derived assuming a fixed-stoichiometry object, i.e., with a constant δ/β ratio, as where δ is the decrement from one of the real part, β is the imaginary part of the complex index of refraction of the sample at the given wave length (n = 1 − δ(λ) + iβ(λ)). This is in effect very similar to the MBA algorithm representing a second-order low-pass filter in the frequency domain. For laboratory propagation-based PCI, where one has a broad energy spectrum, typically an effective wave length is chosen in the denominator of Equation (4) mostly based on visual examination of the retrieved phase image similar to the choice of the α parameter of the MBA method. Therefore, we implemented it by lumping the coefficient in front of |k| 2 into a single, user-tunable parameter, ν, in a slightly modified form as a preprocessing image (Lorentzian) filter as discussed in [32]: I f il is then used in the conventional cone-beam reconstruction pipeline.

Grating-Based PCI
Typical laboratory grating-based PCI setup feature typically a high-power, macro-focus laboratory X-ray source and a G0 source gratings to increase the spatial coherence of the radiation and the usual G1 and G2 gratings forming the Talbot-Lau interferometer [10,12,13,33]. To obtain the attenuation, the differential phase-contrast and dark-field images simultaneously, mostly a phase-stepping approach, usually moving the G1 grating cross-wise in small, equidistant steps is used on such Talbot-Lau setup. The Fourier decomposition of the phase stepping curve approximates it up to the second term as: where χ is the phase step position and p is the grating period of the grating being stepped. This is a sound assumption for rectangular gratings of a reasonable quality. It requires that the phase stepping corresponds to one or several full periods of the grating that is analysed. The mean and amplitude and phase factors in Equation (7) are obtained by pixel-wise fast Fourier transform (FFT) of the stepping images. Then three image modalities are defined as: DFI(r) = a 1,s a 0,re f a 1,re f a 0,s where A, DPC and DFI are the attenuation, differential phase contrast and dark-field projection images, respectively and the subscripts re f and s are for empty beam reference and for sample images. To obtain the CT images of δ of the sample from the DPC images, instead of the usual apodization filters, such as Ram-Lak and Shepp-Logan filters [25] for attenuation CT, a Hilbert-filter is used [10]. Note that for reconstructing both the the phase-contrast and dark-field CT, the option of taking no logarithm of normalised images is built in the package.
Besides the Fourier decomposition, we have also implemented a direct least-squares fit to determine the Fourier components in Equation (7) [33]. In contrast to the Fourier method this method does not need equidistant sampling of the phase stepping curve in Equation (7) in [0, 2π]. For this, Equation (7) is linearized as where and The three parameter a 0,i , b 1,i , b 2,i for pixel i are efficiently determined by least-square matrix inversion using the Penrose pseudo inverse, e.g., in MATLAB. Both methods deliver very similar results in general as is shown in Kaufmann et al. [33].

Results and Discussion
This section discusses some examples demonstrating the functioning of our in-house reconstruction package for different types of phase-contrast micro-CT methods.

Example for the Paganin Method
To demonstrate the functioning of our in-house cone-beam tomographic reconstruction package for laboratory phase-contrast imaging using the Paganin phase-retrieval technique, we will use as an example the propagation-based micro-CT of a horse fly. The scan was performed on an EasyTom XL Ultra 230-160 micro/nano-CT scanner (RxSolutions SAS, Chavanod, France). The scanner features a Hamamatsu reflection-type micro-focus X-ray tube. The scan was performed using a MH110XC-KK-FA, XiRay detector from XIMEA that features a 11 MPixel CCD chip coupled directly via an optical fibre bundle to a 20 µmthick GadOx scintillator and having a native pixel size of 9 µm and 14 bits of dynamic range. The tube was operated at 65 kVp and a current of 120 µA, i.e., at slightly lower than 8 W tube power rendering an emission spot size of around 4 µm. The geometrical magnification was 1.16 and the voxel size of the CT scan was 7.8 µm. Thus, the setup represents a so-called inverse geometry with relatively short propagation distance and a high resolution detector as opposed to the traditional large propagation-distance setup with a large pixel size detector. The images were acquired with 4.0 s exposure time per projection. 1600 projections were taken over 360 degree sample rotation. Results are given in Figure 1 showing the cross-sectional CT slices at two regions: one where the flight muscles are located (see orange line in Figure 4a) and one at the level of the compound eyes (see blue line in Figure 4a). The latter clearly resolves the individual ommatidia. From the images and the corresponding gray-scale histograms, it is clear that the Paganin method strongly increases the CNR of the images and correspondingly the histograms are becoming much more structured, showing several peaks compared to the attenuation case. This means that the different structures corresponding to the different gray-scale peaks can be better resolved in the phase image. This is a well-known property of the Paganin method and the corresponding filter [34,35]. The CNR values for the different methods and samples are summarized in Table 1. A ν value of 10.0 was used in Equation (6) to obtain the images. For a projection image of the horse fly see Figure 4a.

Example BAC Method
In this subsection, we demonstrate the implementation of the BAC method in our in-house reconstruction package and compare to the Paganin method, using a propagationbased micro-CT of a coffee bean. The scan was performed using the same scanner as in the previous subsection. The scan parameter were slightly different: the tube was operated at 65 kVp and a current of 150 µA.The geometrical magnification was 2.53 and the voxel size of the CT scan was 7.1 µm. The CCD camera was used with 2 × 2 pixel binning. The images were acquired with 4.0 s exposure time per projection. 1568 projections were taken over 360 degree sample rotation. Figure 2a shows a projection image of the coffee bean and in Figure 2b shows the result of the optimization process for the α and γ parameters as is explained in the previous section. Micro-CT images of the coffee bean are given in Figure 3. These show that the separation in the gray-value histogram increases, i.e., both the Paganin and BAC methods increase the CNR of the reconstructed images (see Table 1 for actual values). Furthermore, the magnified regions of interest (ROIs) framed in orange show how the different retrievals might blur the visibility of small scale feature (see, e.g., red arrows). The BAC method, as it has been reported earlier [29,30,36], improves significantly the CNR compared to the attenuation-based image but, in general, blurs less the images than the Paganin method and therefore enables resolving very tiny features (see red arrows in Figure 3).
Finally, we compare the attenuation-based, Paganin phase retrieval and BAC methods for the horse fly sample from the previous subsection as well. This is shown in Figure 5.

Example for Grating-Based PCI
In this paragraph, we show the functioning of our in-house cone-beam tomographic reconstruction package for laboratory phase-contrast imaging using the grating based technique as described in Section 2.2. The sample was a very-low-attenuation, porous cellulose foam sample. The scan was performed using our in-house grating-based PCI micro-CT setup. Technical details of the setup are given in [33]. The scan was performed at 40 kVp tube voltage and a tube current of 20 mA. We performed the scan with seven phase steps per projection each with 2.5 s exposure time. We took overall 360 projections over 360 degree sample rotation. The geometrical magnification of the setup is very low 1.17. Figure 6 shows some illustrative results proving the power of PCI and DFI for lowattenuation, porous samples. The attenuation image is so noisy due to the very low sample attenuation that the internal porous structure of the sample cannot be discerned. Only the phase-contrast and the dark-field CT images show these structures. They nicely illustrate inverse gray scale values with respect to each other as expected: pores are shown with low intensity in the PCI, whereas they have high intensities in the DFI and vice versa for the cellulose material.  Figure 4a). The areas circumvented in red indicate the ROIs used for determining the CNR values between sample material and air (see Table 1).

Conclusions
We introduce and describe an in-house, cone-beam tomographic reconstruction package for laboratory phase-contrast imaging. The package covers different phase-contrast techniques and phase retrieval methods. It includes the well-known Paganin method and the Bronnikov-aided Correction method for propagation-based phase-contrast imaging. It also includes reconstruction method for the grating-based phase-contrast and dark-field micro-CT imaging. In the paper, we explain the implementation of these methods and their integration in the filtered back projection pipeline. Their functionality and efficiency were demonstrated by performing phase-contrast micro-CT on a few dedicated bio-samples.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: