Improved iterative tomographic reconstruction for x-ray imaging with edge-illumination

Iterative tomographic reconstruction has been established as a viable alternative for data analysis in phase-sensitive x-ray imaging based on the edge-illumination principle. However, previously published approaches did not account for drifts of optical elements during a scan, which can lead to artefacts. Up to now, the strategy to reduce these artefacts was to acquire additional intermediate flat field images, which were used to correct the sinograms. Here, we expand the theoretical model to take these effects into account and demonstrate a significant reduction of (ring)-artefacts in the final reconstructions, while allowing for a significant reduction of scan time and dose. We further improve the model by including the capability to reconstruct combined absorption and phase contrast slices, which we experimentally demonstrate to deliver improved contrast to noise ratios compared to previously employed single shot approaches.

detector pixels. A scan of the sample mask over one period (i.e. varying m) provides a Gaussian-like intensity distribution, which is called the illumination curve (IC). The absorption signal appears as decreased intensity over the entire IC, while the refraction (i.e. differential phase) signal shifts the IC (see inset in figure 1).
Standard analysis for tomography consists of two steps. First, one  or several radiograms (Hagen et al 2015 are analysed to provide joined or separate absorption and differential phase contrast images for each viewing angle θ. Second, filtered back projection (Kak and Slaney 1988) is applied in order to perform tomographic reconstruction. Using an iterative approach to data analysis, it was recently demonstrated that tomographic reconstruction can also be achieved in a single step (Chen et al 2017).
In the following, we will substantially extend the previously published model for iterative reconstruction in EI, which features the suppression of ring artefacts and opens the possibility for significantly reduced scan times.

Theory
The first steps of developing the theoretical frame work have been already been described in literature , Chen et al 2017 and are reiterated here for the convenience of the reader. The Radon transform constitutes the mathematical basis for tomographic reconstruction and can be expressed as the line integral (Kak and Slaney 1988) for an arbitrary, two-dimensional function h (x, y) with (x, y) the coordinate system fixed to the sample and x = (x, y). L defines the line of integration, which is parameterised by the pixel position at the detector t and the projection angle θ (see figure 1). The absorption signal A after transmission through the sample constitutes the negative exponential of the projected (i.e. Radon transformed) x-ray attenuation coefficient µ(x, y) The accumulated phase shift Φ(t, θ) of the x-rays during transmission through the sample is given by where K is the modulus of the wave vector of the assumed monochromatic x-ray beam and δ(x, y) is the local refractive index decrement. The observable refraction signal α(t, θ) is proportional to the derivative of the phase signal (Born and Wolf 1999) Now, let us denote intensity of ICs acquired without a sample present in the x-ray beam as f (t, m) with m, the lateral offset of the sample mask (i.e. the scan parameter) and t, the pixel position at the detector. Note that f 's dependency on t takes into account local mask imperfections. The sample's absorption signal, A(t, θ), decreases the total signal while the refraction signal, α(t, θ), shifts the experimentally obtainable ICs s (t, θ, m). Both effects can be modeled by where z is the sample to detector distance and we have assumed negligible scattering from the sample, which would additionally broaden the ICs (Endrizzi et al 2014, Modregger et al 2016. Substituting the absorption and refraction signal defined above, yields this model's prediction for the observable sinogram s mod based on the attenuation coefficient µ and refractive decrement δ: In Chen et al (2017), the flat field IC f was further simplified by a first order Taylor approximation at the mask position m, which we will forgo here. This implies that our model will be better suited for larger refraction angles. The Taylor approximated version of equation (6) formed the basis for iterative tomographic reconstruction in Chen et al (2017). This was done by minimising the cost function with a batch gradient iterative approach including an additional noise reduction term. Chen et al demonstrated that the absorption and phase signal can be retrieved independently from a scan with a single mask position m. However, this model does not account for a shifting mask position during a tomographic scan, which can occur due to mechanical vibrations or drift. Since these can lead to significant artefacts in the tomographic reconstructions (Zamir et al 2016), we will include an additional degree of freedom in the argument of f denoted m o (θ), which models an offset position for each projection angle θ of a rigid mask (i.e. a constant offset over the field of view).
Similarly, reconstructed slices can show strong ring artefacts, which can occur due to a mismatch between the actual flat field ICs f (t, m) during a tomographic scan and the reference ICs acquired before, where the latter is used for iterative reconstruction in equation (6). This mismatch may result from systematic mechanical errors in the optical elements or from insufficient photon statistics for the reference ICs. A heuristic approach for taking this effect into account is to include another degree of freedom in the argument of f denoted m r (t). This approach works as follows. In an ideal experiment without noise, ring artefacts would occur by using the locally shifted flat field ICs f (t,m − m r (t)) instead of the correct f (t, m). Conversely, including m r (t) in the model allows for reducing ring artefacts during iterative reconstruction.
Finally, we will take advantage of the previously published approach to retrieve a combined absorption/ phase contrast image by assuming the phase signal of the sample to be proportional to the absorption signal, i.e. µ(x, y) = γδ(x, y) = h(x, y) (Paganin et al 2002). It was shown that the projected electron density ρ(t, θ) can be retrieved from a projection s(t, θ, m) acquired on a single offset position m by calculating with F the Fourier transform with respect to t, its inverse F −1 and q the variable conjugate to t. Here, we have also omitted the free space propagation part, which can be neglected for laboratory-based setups. Standard filtered back projection (Kak and Slaney 1988) can then be applied to the retrieved sinogram ρ(t, θ) in order to reconstruct a tomographic slice. Here we will use this well established single shot approach as a reference for comparison later. While the utilised assumption of proportionality between absorption and phase contrast is strictly true only in the case of a single homogeneous material present in the beam and precludes a quantitative interpretation for more complex material systems, in practice qualitative interpretation of the images such as morphological information (e.g. pore size distribution) is frequently of interest. Further, this assumption has been demonstrated to improve image quality for EI in some cases (Zamir et al 2016) and simplifies the iterative reconstruction as the number of retrieved values are only half compared to the previously published method (Chen et al 2017).
Including the three improvements listed above, yields the proposed model for iterative tomographic reconstruction: where the dependency of h on (x, y) was suppressed for readability. Retrieving the combined contrast h(x, y) from the experimentally obtained sinogram s exp (t, θ, m) can now be achieved by an iterative minimisation of the cost function where the second term constitutes a regularisation term that minimises noise with a simple L2-norm. Please note that the number of degrees of freedom for the reconstruction does not increase significantly by including m o and m r as additional parameters. For example, if the number of acquired projections is N θ and the number of pixels at the detector is N t , then the number of voxels in the slice h(x, y) to be retrieved would be N t × N t , while m o (θ) constitutes only N θ and m r (t) only N t additional parameters.
We have used the limited memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) implementation (Byrd et al 1995) in the SciPy library for Python (Jones et al 2001) to minimise the cost function (equation (10)). The knowledge of the gradients of S with respect to h(x, y), m o (θ) and m r (t) greatly improves iteration speed. Using the abbreviations and the gradients are, with respect to h(x, y) with respect to m o (θ) and with respect to m r (t) The Radon transforms and their gradients occurring in equations (9)-(15) can be regarded as sparse matrix operations that can be efficiently computed by using lookup tables. The lookup tables were computed prior to iterative reconstruction and used to calculate the corresponding Radon transforms and gradients.

Experimental results
The experiments were performed in-house with the laboratory-based EI implementation at University College London (Hagen et al 2014a). Two different experiments with different type of samples and different setup configurations were carried out to demonstrate the versatility of the proposed approach.
In the first experiment, we have imaged three different types of plastic spheres (i.e. polystyrene (PS), polypropylene (PP) and poly methyl methacrylate (PMMA)) as a straight forward example. The source was a Rigaku M007 rotating anode (Rigaku Corporation, Japan) with a Mo target and operated at 40 kV and 20 mA. The Hamamatsu C9732DK flat panel (Hamamatsu, Japan) was used as a detector with a pixel size of 50 µm. Both masks were manufactured by electroplating of Au on a graphite substrate (Creatv Microtech Inc., USA) featuring a pitch of 38 µm for the sample mask and 48 µm for the detector mask, respectively. The total setup length was 85 cm and the sample to detector distance was 20 cm.
In order to obtain the experimental flat field ICs f (t, m) the sample mask was scanned over one period with 33 steps. The tomographic scan was performed with the sample mask on a slope position of the ICs (i.e. m = 9 µm offset from the maximum position; see inset of figure 1) and over 360 degrees under continuous rotation of the sample with a speed of 1 deg s −1 leading to a total scan time of 6 min. During the scan, 300 projections were acquired with an exposure time of 1.2 s each. Subsequent iterative reconstruction of the acquired data set was carried out with γ = 3 for both reconstruction methods and λ = 5 × 10 −3 for the iterative method, which were manually found by visual inspection. The retrieved tomographic slice consisted of N t × N t = 220 × 220 voxels. The iteration was carried out until satisfying convergence was reached, which took 58 iteration steps and 93 s on a single core of a standard modern desktop PC. Figure 2 demonstrates a clear suppression of ring artefacts provided by the proposed iterative approach (b) as compared to the result of the single shot algorithm (a). Further, the iterative approach provides visually less noise, which we confirmed by determining the contrast to noise ratio (CNR) between the different materials. The CNR is defined as with h the mean of an area and std(h) its standard deviation. Table 1 lists the CNRs between the different structures as determined by a circular area with a radius of 40 voxels. The iterative approach provides an increase of about 30% between the materials. In order to assess the validity of the offsets for the lateral mask position m 0 (θ) and for the suppression of ring artefacts m r (t) retrieved by iterative reconstruction, we have compared them to predicted values as follows. For the offset of the lateral mask position m o (θ), we used a background area of the acquired sinogram. Here, the offset of the lateral mask position during a scan due to drift or vibrations will appear as intensity variations I in dependence of the projection angle θ. Thus, the known flat field IC f (t, m) can be used to transform the intensity variations into a lateral offset m o (θ) by using its inverse, i.e. m o (θ) = f −1 (t, I), and numerical interpolation. We found an excellent agreement between the thusly predicted and iteratively retrieved values for m 0 (θ) (figure 3(a)) featuring a correlation coefficient of r = 0.98.
Predicting the offset for ring artefact suppression m r (t) is not straight forward since the sample and ring artefact information overlap in the sinogram. However, we demonstrated self consistency of the retrieved m r (t) values in the following way. We used the iteratively retrieved slice h(x, y), the retrieved offset m o (θ), the exper-  (8)) and (b) iterative tomographic reconstruction (equations (9) and (10)) from the same data set. A very significant reduction of ring artefacts in (b) is clearly visible. Table 1. CNR between the different areas (first column) of the single shot reconstruction (second column), the iterative reconstruction (third column) and the relative difference between the latter two (fourth column). Iterative reconstruction provides improved CNR of about 30% between the materials.

Area
Single imental flat field IC f (t, m) in equation (9) with m r (t) ≡ 0 in order to obtain a sinogram that was virtually unaffected by ring artefacts. The difference between this sinogram and the experimentally acquired one then quantified the intensity modulations due to the ring artefacts, which were again transformed into the estimated offsets via the flat field IC f (t, m). Figure 3(b) shows a very good agreement between the estimated and the retrieved offsets m o (t) with a correlation coefficient of r = 0.85. The large correlation coefficients found here and in the numerical simulations above validated our suggested extension for iterative tomographic reconstruction. For the second experiment, we repurposed a previously acquired data set of a freeze-dried rabbit esophagus as an example for a biomedical sample and the versatility of the proposed approach. Details of sample preparation and ethical approval can be found elsewhere (Hagen et al 2015). We used the same x-ray source as above but operated at 40 kV and 25 mA. The detector was a Pixirad-2 photon counting detector (PIXIRAD Imaging Counters s.r.l., Italy) featuring a pixel size of 62 µm. The periods of the masks were 79 µm and 98 µm, respectively. Total setup length was 2 m and the sample to detector distance was 0.4 m.
The flat field ICs f (t, m) were acquired with 30 steps evenly distributed over one period of the sample mask. The tomographic scan was performed with the sample mask on a slope position of the ICs (i.e. m = 10 µm offset from the maximum position) and over 180 degrees in 900 steps ( N θ = 450 of which were used here) with 2 s exposure per projection. Sample dithering was used with 8 steps to improve spatial resolution (Hagen et al 2014b). An additional flat-field was taken at each projection angle for the purpose of flat field tracking (Zamir et al 2016), which was not used for the data analysis in this study. Due to the overhead of moving the motors this resulted in a total scan time of ≈ 18 h. The reconstructed slice consisted of N t × N t = 720 × 720 voxels and iterative reconstruction was performed with γ = 20, λ = 2 × 10 −2 and took 47 steps and about 20 min. Figure 4 demonstrates a substantial reduction of ring artefacts provided by iterative reconstruction compared to standard single shot analysis. In fact, some sample details (e.g. at the center of rotation) that are not visible in the single shot slice are recovered by the iterative approach. However in this example, the separation of the sample from the background is slightly worse for the iterative result. Numerical investigations not shown here revealed that this a consequence of using a larger proportionality factor between absorption and phase signals (i.e. γ = 20 instead of γ = 3 used above). Still, the slice provided by iterative reconstruction is significantly preferable for consecutive data analysis.
As mentioned above, the scan for the utilised data set was set up for flat-field tracking in order to reduce ring artefacts, which requires the acquisition of additional flat-field images during the tomographic scan. Moving the sample out and back into the field of view caused a significant time overhead as can be seen by comparing the actual scan time of 18 h to the time used for pure data acquisition of 900 × 8 × 2 s = 4 h. Since we demonstrated that iterative reconstruction renders flat field tracking obsolete, the proposed approach has the potential to decrease total scan times by a factor of more than four.

Numerical simulations
To further demonstrate the validity of the proposed iterative approach we used numerical simulations with parameters that model the experiment described above. The numerical simulations were carried out in two parts.
First, the forward simulation of the observable intensity was based on equation (9) and yielded the sinogram used as input for the subsequent iterative reconstruction. The simulated sample (i.e. h(x, y)) consisted of three circles with different materials and equations (1)-(4) were used to simulate the Radon transform of the absorption, A(t, θ), and refraction contrast, α(t, θ), respectively. The flat field IC (i. e. f (t, m)) was modelled as a Gaussian distribution and the scan parameter, m, was chosen on one inflection point. Random offsets for the lateral mask position, m o (θ), and for inducing ring artefacts, m r (t), were introduced into the flat field IC. The simulated sinogram was then calculated by linear interpolation of f (t, m) in equation (9) with the input data described above. Finally, Poisson noise was added to both the sample sinogram, s(t, θ, m) as well as the flat field IC f (t, m).
Second, iterative tomographic reconstruction was carried out based on equations (9) and (10). The sinogram consisted of N t = 180 detector pixels and N θ = 90 projections. The proportionality factor was chosen to match one of the materials (i.e. γ = 10) and iteration was performed until convergence was achieved. Figure 5 demonstrates the influence of the different utilised parameters for artefact suppression on iteratively reconstructed slices. The single shot reconstruction (figure 5(a)) based on equation (8), included for compariso n, shows clear ring artefacts. These artefacts are significantly reduced by iterative reconstruction that uses both noise regularisation and includes the offsets as additional parameters (cmp. figure 5(d)). Further, the input offsets were highly correlated with the iteratively retrieved values (r = 0.9965 for m o and r = 0.9075 for m r ).
In addition, we turned off either the noise reduction term in figure 5(b), i.e. setting λ = 0, or the offsets in figure 5(c), i.e. setting m o = m r = 0. The expected effects (occurrence of ring artefacts or significant noise) demonstrates that the additional degrees of freedom in the iterative reconstruction have their desired effects.

Conclusions
We have substantially extended and experimentally validated a previously developed model for iterative tomographic reconstruction of data sets acquired with an x-ray EI setup. We included additional degrees of freedom that account for offsets in the sample mask position due to drift or vibrations and that allows for a significant reduction of ring artefacts. This rendered flat field tracking obsolete and allows for a reduction of total scan time. In addition, we incorporated the possibility to retrieve a combined absorption/phase signal, which was demonstrated to yield higher inter-material CNRs than established single shot approaches to data analysis. Further, we showed the versatility of the proposed iterative reconstruction frame work by successful application to different experimental setups as well as scan types/parameters.