Reverse projection retrieval in edge illumination x-ray phase contrast computed tomography

Edge illumination (EI) x-ray phase contrast computed tomography (CT) can provide three-dimensional distributions of the real and imaginary parts of the complex refractive index (n=1−δ+iβ) of the sample. Phase retrieval, i.e. the separation of attenuation and refraction data from projections that contain a combination of both, is a key step in the image reconstruction process. In EI-based x-ray phase contrast CT, this is conventionally performed on the basis of two projections acquired in opposite illumination configurations (i.e. with different positions of the pre-sample mask) at each CT angle. Displacing the pre-sample mask at each projection makes the scan susceptible to motor-induced misalignment and prevents a continuous sample rotation. We present an alternative method for the retrieval of attenuation and refraction data that does not require repositioning the pre-sample mask. The method is based on the reverse projection relation published by Zhu et al (2010 Proc. Natl Acad. Sci. USA 107 13576–81) for grating interferometry-based x-ray phase contrast CT. We use this relation to derive a simplified acquisition strategy that allows acquiring data with a continuous sample rotation, which can reduce scan time when combined with a fast read-out detector. Besides discussing the theory and the necessary alignment of the experimental setup, we present tomograms obtained with reverse projection retrieval and demonstrate their agreement with those obtained with the conventional EI retrieval.


Introduction
X-ray phase contrast computed tomography (CT) exploits phase shifts (refraction) in addition to attenuation for contrast generation. Specifically, this imaging modality can provide three-dimensional maps of both the real and imaginary parts of the complex refractive index within a sample: x y z k x y z k , , ; 1 , , ; i , , ; , where δ and β denote local phase shift (refraction) and attenuation properties, respectively, (x, y, z) are the sample coordinates and k is the wave number. Such bimodal imaging can provide additional information for samples exhibiting weak attenuation contrast [1]. Until now, several x-ray phase contrast CT techniques have been developed [18]; the one we use here is based on the edge illumination (EI) principle [14]. By illuminating only the edges of a row of pixels with a narrow (typically <20 μm) laminar beam, sensitivity to refraction, i.e. the local manifestation of phase shifts, is achieved: beamlets refracted towards/away from the pixels' active areas cause an increase/decrease in the measured intensity. When combined with conventional x-ray sources, EI is often implemented in full-field mode [15] ( figure 1(a)). In this setup, a mask in (Some figures may appear in colour only in the online journal) Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. front of the sample ('pre-sample mask') splits the beam into an array of physically separated beamlets, and a second mask in front of an area detector ('detector mask') creates insensitive areas ('edges') between the pixels. By positioning the pre-sample mask such that each individual beamlet falls partially on a pixel and partially on an absorbing detector mask septum, the EI principle is replicated over the entire field of view (FOV). The requirement to keep the beamlets spatially separated imposes an upper limit on the tolerable source blur, which depends on the mask design and the system dimensions [16]. Ultimately, the blur also determines the phase sensitivity of the setup [3].
Assuming beamlets of unit intensity, a projection (P) acquired with such a setup can be described by [2,11]: ) are the sample-induced attenuation and refraction, z 3 is the sample-to-detector distance and M is the magnification between the pre-sample mask and the detector. Note that only refraction in the x-direction is considered since this is the system's direction of sensitivity as determined by the orientation of the mask apertures. The symbol C denotes the illumination curve which can be measured in the absence of a sample by step-scanning the pre-sample mask along the x-direction and recording the intensity at each step ( figure  1(b)); x m is the pre-sample mask position used for imaging. Equation (1) is based on the fact that a refraction by angle ( ) α x y , causes a beamlet shift of ( )/ α z x y M , 3 , which is equivalent to a shift introduced by a displacement of the presample mask. In the derivation of equation (1), small-angle scattering has been neglected, although the model can easily be extended to scattering samples [5,6]. Both half slopes of the illumination curve correspond to a partial illumination of the detector mask septa, yet opposite septa are illuminated ( figure 1(b)). Previous work has shown that A(x, y) and ( ) α x y , can be retrieved when two projections are acquired, one on each half slope [2,11]. In CT mode, i.e. when the sample is rotated and two projections are acquired at each CT angle, this yields attenuation and refraction sinograms for each transverse sample slice ( = y const.): In the above, θ denotes the CT angle, and cos ) are the coordinates of the rotating sample. When θ covers a range of at least 180°, tomographic maps of β k ('attenuation tomograms') and δ ('phase tomograms') can be reconstructed through standard methods, e.g. filtered back projection (FBP). Note that, for polychromatic radiation, tomograms refer to effective energies [12]. The derivative in the refraction sinogram makes EI belong to the class of differential phase contrast (DPC) techniques, and implies that an additional integration step or the use of a dedicated filter function in the FBP is required in the reconstruction [13].
The necessity to acquire projections in opposite illumination configurations at each CT angle has two major drawbacks. First, the displacement of the pre-sample mask from x m1 to x m2 can cause inaccurate mask alignment due to a limited motor precision. Second, each displacement requires an interruption of the sample rotation, thus preventing a continuous sample rotation. In this paper, we present an alternative retrieval method based on the concept of reverse projection, which was first published by Zhu et al (2010) and later analysed in more detail by Wu et al (2015) for grating interferometry-based x-ray phase contrast CT [20,21], but which has not yet been applied to EI-based x-ray phase contrast CT nor implemented with a conventional x-ray source. This concept relies on the fact that, after a 180 degree rotation of the sample, refraction occurs in the opposite direction, while the magnitude of the refraction angle and the total attenuation remain the same. As shown in the following, when applied to EI-based PC-CT, these relations allow extracting A(x, y) and ( ) α x y , from two projections acquired in the same illumination configuration (i.e. with the pre-sample mask located in the same position), but with a rotation offset of 180° between them. This implies that a complete CT dataset can be collected by rotating the sample over 360° instead of 180°, without the need to displace the pre-sample mask. While the total number of images remains the same, reverse projection retrieval simplifies the acquisition protocol and eliminates the potential for motorinduced misalignment. Moreover, it allows the sample to be rotated continuously. For fast read-out detectors, this typically leads to reduced scan times. In the following, we describe the theory of reverse projection retrieval for EI-based x-ray phase contrast CT. Moreover, we discuss a necessary experimental requirement, namely a specific alignment condition of the imaging setup. Experiments with a laboratory-based x-ray tube demonstrate that reverse projection retrieval leads to images in agreement with those obtained with the previously used, conventional EI retrieval. Finally, tomograms acquired with a continuous sample rotation and reconstructed using the reverse projection retrieval are shown to provide quantitative accuracy.

Reverse projection retrieval
As was mentioned above, conventionally in EI-based phase contrast CT the retrieval of attenuation and refraction data is performed on the basis of two projections acquired with different pre-sample mask positions at each CT angle [2,11]. Assuming a parallel beam, this retrieval requires a 180° rotation of the sample for a complete CT dataset. Let us assume instead that only one projection is acquired at each CT angle (with the pre-sample mask positioned at x m1 , for example), but that the rotation is performed over 360°. Let us consider two projections acquired with a 180° rotation offset between them, with the rotation axis located at x = 0: Here, = y const. has been assumed. As rotating the sample by 180° does not change the total amount of attenuation, we can write ( ) , . The same holds for the refraction angle; however, since the refraction angle is a directional quantity, its sign changes: Equation (6) is the so-called reverse projection relation [21]. The above enables equation (5) to be rewritten as: Dividing the left and right hand side terms of equations (4) and (7), i.e. the projections acquired with a rotation offset of 180 degrees between them, yields: This can be expressed as: where: is a function that can be calculated from the illumination curve (C), the position of the pre-sample mask (x m1 ), the sample-todetector distance (z 3 ) and the magnification (M). The refraction angle ( ) α θ x, can now be extracted by applying the inverse of F to the ratio of the projections: It should be noted that for the calculation of F and its inverse no linear approximation of the illumination curve is required, allowing relatively large refraction angles to be retrieved (typically up to tens of μrad). Moreover, equation (11) tolerates non-uniformities of the pre-sample and detector masks as the function F can be calculated on a pixel-by-pixel basis. The retrieval of ( ) α θ x, is only constrained by the fact that F has to be restricted to its injective part (i.e. where it describes a oneto-one mapping) [11].
Following the extraction of the refraction angle, the attenuation can be retrieved via: The application of equations (11) and (12) to projections acquired over 360° yields sinograms like those described by equations (2) and (3), which enable the reconstruction of tomograms of β k and δ. The reverse projection relation holds for parallel beam geometries only, a criterion which is either satisfied at synchrotrons or when the fan angle produced by conventional x-ray sources is sufficiently small (i.e. when the distortion introduced by the fan is smaller than or of the same order as the spatial resolution of the imaging system). A generalization of the reverse projection relation to fan beam geometries is suggested in [19].

Sampling and alignment of the experimental setup
In the above equations, the x-coordinate was assumed to be a continuous variable. However, in practice, projections are sampled at discrete locations x x , ..., n 0 . If no additional sample movement is applied except rotation, the sampling step in the projections is given by the period of the pre-sample mask ( p ), i.e.
= + x x kp k 0 . This can be artificially increased through 'dithering'; during this procedure, the sample is displaced multiple times by a fraction of p, projections are acquired at each displacement, and then combined into a single projection that features a higher sampling rate [4,10]. In the context of reverse projection retrieval, it is important to note that the extraction of ( ) θ A x, and ( ) α θ x, via equations (11) and (12) is only accurate if projections acquired with a 180° offset between them are sampled at the same locations x k . This is the case when the centre of rotation (COR) of the sample stage is aligned with the centre of one of the beamlets. If such alignment is not realized, data sampled at different locations are combined, resulting in a loss of accuracy in the retrieved attenuation and refraction images. In the following, we present a simple procedure for the alignment of the COR of the sample stage with the centre of one of the beamlets. The procedure is schematically illustrated in figures 2(a)-(c).
(a) The pre-sample mask is placed at x m1 and the detector mask is removed from the setup, such that the beamlets impinge directly on the detector pixels. A highly absorbing bar is mounted on the sample stage. A single detector pixel ('pixel 1') that lies in the vicinity of the bar edge marked by the white asterisk, is selected. The bar is then scanned in sub-pixel steps until its edge has gone across all of 'pixel 1'. The intensity measured in 'pixel 1' is plotted ( figure 2(a)). (b) The bar is moved back to its original position, and then the sample stage is rotated by 180° (note that the sketches in figure 2 are not to scale and that the sample stage can be rotated without the bar colliding with the detector in the real setup). A second detector pixel ('pixel 2') in the vicinity of the bar edge marked by the white asterisk, which now is in a different position, is selected. Importantly, 'pixel 2' must be selected such that it is separated from 'pixel 1' by an odd number of pixels. The bar is scanned again with sub-pixel steps until its edge crosses 'pixel 2'. Again, the intensity curve is plotted ( figure 2(b)).
(c) The bar displacement (∆x) which yields the same measured intensity in 'pixel 1' and 'pixel 2' before and after the 180° rotation of the sample stage is the one for which the COR is located in the centre of one of the beamlets ( figure 2(c)). The sample stage is displaced by ∆x from its original position, which completes the alignment.

Experimental setup
All experiments were performed with a Rigaku MicroMax 007 HF rotating anode (molybdenum) x-ray tube (Rigaku Corporation, Japan) operated at 40 kV and 25 mA. Previous work has shown that all spectral energies contribute to an EI projection, with the strongest contribution from the Mo K-alpha peak (around 17 keV) [8,12]. The source's focal spot size was measured to be approximately 70 μm (full width at half maximum). The Hamamatsu C9732DK flat panel (Hamamatsu, Japan), a passive-pixel CMOS sensor with a pixel size of µ × 50 50 m 2 , was used as detector. The distances between source and pre-sample mask (z 1 ), pre-sample mask and sample (z 2 ), and sample and detector (z 3 ) were 1.6, 0.05 and 0.35 m, respectively. Both masks were fabricated by electroplating gold strips onto a graphite substrate (Creatv Microtech Inc., Potomac, MD). The masks' periods were 79 μm (pre-sample mask) and 98 μm (detector mask), with apertures of 10 μm and 17 μm, respectively. In this setup, the detector mask apertures approximately match the source size projected onto the detector plane, but in general other aperture widths are also possible. With this mask design, only every other detector column is used for the acquisition ( figure 3). This 'line skipping' configuration reduces the negative effect of pixel cross-talk on EI data [10]. Skipping every other detector column does not affect the intrinsic spatial resolution of the imaging system, but reduces the spatial sampling rate of the projections, which, if not compensated by dithering, can cause a reduced spatial resolution in the reconstructed tomograms. The position of the pre-sample mask (x m1 ) was chosen such that 50% of each beamlet fell on a pixel and 50% on an absorbing septum of the detector mask.

Centre of rotation alignment
Before imaging, the COR of the sample stage was aligned with the centre of one of the beamlets according to the procedure described above. Figure 4 shows the intensity curves that were measured by two selected detector pixels ('pixel 1' and 'pixel 2') while scanning an absorbing bar through the setup before and after the sample stage was rotated by 180°. Since the bar did not have a perfect edge, the measured intensity curves have a different shape compared to the ones sketched in figure 2. The displacement for which the two curves overlapped (∆x) was extracted and the sample stage was moved by this amount from its original position. Note that the step with which the bar is scanned ultimately determines the precision of the alignment. In this case, the bar was scanned with 0.5 μm steps, implying an alignment precision of 5% with respect to the spatial resolution of the imaging system (10 μm). As was shown by Diemoz et al (2014), when implemented with conventional x-ray sources, the spatial resolution of an EI imaging system is given by the smaller of the pre-sample mask aperture (in this case 10 μm) and the projected source size scaled to the sample plane (i.e. the penumbra-effect, in this case approximately 14 μm) [4].

Proof-of-principle scan
To obtain a proof-of-principle for the validity of the reverse projection approach, a chicken bone surrounded by adjacent soft tissue was scanned. The 'fresh' specimen (not fixed in formalin) was placed in a plastic cylinder of approximately 8 mm diameter. This sample size justifies the assumption of a parallel beam, since the distortion introduced by the fan beam is   below 10 μm (according to geometrical arguments), matching the spatial resolution of the imaging system. Projections were acquired in both reverse projection (with a 360° rotation, in one illumination configuration) and conventional EI mode (with a 180° rotation, in two opposite illumination configurations). In both cases, the rotation was performed in a stepped manner with an angular step of 0.5°. The integration time per projection was 3 s. In order to increase the spatial sampling rate, dithering was performed by displacing the sample 8 times by an 8th of the pre-sample mask period at each CT angle. Projections were corrected for flat field inhomogeneities and dark current before being processed according to equations (11) and (12) (reverse projection retrieval), and according to [2] (conventional retrieval). Attenuation and phase tomograms were reconstructed with FBP with the ramp and Hilbert filters, respectively. Figure 5 shows the function F, which was calculated from the illumination curve and which relates the ratio of two projections acquired with a 180° offset to the refraction angle. Its non-symmetric behavior is due to it being defined as the ratio of two shifted versions of the illumination curve that can be well-approximated by a Gaussian function (see figure 1(b)). F takes on the value 1 for α = 0 since the input to it is the the ratio of two projections ( ( )/ ( ) θ θ π + P Mx P Mx , , ) that are identical if no refraction occurs. The injective part of F, which allows a unique retrieval of α, is indicated by the thicker line. This shows that, for the used experimental setup, a range of up to 100 μrad can be retrieved, although this range may be smaller for noisy data.
The reconstructed attenuation and phase tomograms, as well as extracted line profiles, are shown in figure 6. Generally, it can be seen that the phase tomograms (middle row) clearly show the soft tissue surrounding the bone. Conversely, the soft tissue is hardly visible in the attenuation tomograms (top row). This is not the result of inadequate grey value windowing in the presence of a strong attenuator (bone), but is due to the fact that the weak soft tissue contrast is hidden in the noise floor, as apparent in the line profiles and the zoomed insets where the contrast is stretched. These results confirm the superiority of phase over attenuation contrast for weakly attenuating specimens. With respect to the used retrieval methods, it can be appreciated that reverse projection (figures 6(a) and (c)) and conventional EI retrieval (figures 6(b) and (d)) yield matching results for both attenuation and phase, as the strong phase-based soft tissue contrast can be equally appreciated in both cases. This is also confirmed by the corresp onding line profiles.

Continuous scan
One advantage of reverse projection over conventional EI retrieval is that it removes the necessity to stop the sample rotation at every CT angle, thus allowing continuous scans. To demonstrate that the reverse projection retrieval in combination with a continuous rotation yields accurate results, a custom-built phantom consisting of three plastic tubes (diameter: approx. 3 mm) filled with olive oil, water and air was scanned. The total phantom diameter was approximately 8 mm, again justifying the use of the parallel beam approximation. First, a reverse projection scan was performed (over 360°, one illumination configuration), while rotating the sample continuously with a speed of 0.25° per s, and with a detector integration time of 3 s per projection. Taking into account the time needed for detector read-out (1 s, in the used configuration of the detector), this corresponded to a 1° angular increment between projections. The rotation was briefly paused after the first half scan in order to avoid a possible mismatch between projections with a 180° offset caused by a non-uniform motor speed during acceleration. Pausing the scan after the first 180° means that exactly the same initial acceleration applies to both complementary halves of the scan. Projections were flat-fielded and dark current corrected, and processed according to equations (11) and (12). Attenuation and phase tomograms were reconstructed via FBP with a ramp and a Hilbert filter, respectively. For comparison, the phantom was also scanned in conventional EI mode (over 180°, in two illumination configurations) with a stepped rotation, an angular step of 1° and 3 s integration time per step, attenuation and refraction projections were retrieved according to [2], and tomograms were reconstructed via FBP with a ramp and a Hilbert filter, respectively.
Figures 7(a)-(d) show the attenuation and phase tomograms of the phantom obtained from the reverse projection and conventional EI scans. Visually, the images are very similar. For a quantitative comparison, the reconstructed values inside the water and oil-filled tubes (left and top cylinders, respectively) were averaged over the central × 15 15 pixels within the tubes. The results are shown in figures 7(e) and (f). Error bars represent one standard deviation of the averaged pixel values. A good quantitative match between reverse projection and conventional EI scans can be seen for both attenuation and phase, proving the validity of the reverse projection approach when combined with a continuous sample rotation. Since a polychromatic x-ray beam was used, the reconstructed values correspond to effective energies that are determined by the sample material, the source spectrum and the spectral characteristics of the masks and the detector [12]. Although a comparison with theoretical values was not performed in this case, it has already been proven that, provided effective energies are correctly calculated according to the model by , EI-based x-ray phase contrast CT yields reconstructed values which are in agreement with theoretical ones (see supplementary material of [9]).

Discussion and conclusion
An alternative method for the retrieval of attenuation and refraction data from EI x-ray phase contrast CT scans, based on the reverse projection relation [21], was presented. Besides providing the necessary theory, it was highlighted that a specific alignment of the COR of the sample stage with the centre of one of the beamlets created by the pre-sample mask is essential in order to realize this experimentally. A simple alignment procedure that allows this was described. Moreover, it was demonstrated that reverse projection retrieval yields results fully compatible with those obtained through the previously used, conventional EI retrieval, both in terms of soft tissue contrast and quantitative accuracy.
The main advantage of reverse projection over conventional EI retrieval is that the acquisition is simplified and less prone to motor-induced misalignment. Moreover, the requirement to stop the sample rotation in order to acquire projections in opposite illumination configurations is eliminated; hence, the reverse projection approach allows rotating the sample continuously. With a continuous rotation, scan time is determined by a combination of the desired angular increment, integration time per projection and detector readout time, but it is no longer affected by additional movements of the pre-sample mask motor. Therefore, in the commonly encountered situation where the detector read-out is shorter than the time it takes to displace the pre-sample mask from x m1 to x m2 , this leads to faster acquisitions. Generally speaking, scan time could be further reduced if reverse projection scans were performed with upgraded experimental equipment, i.e. with a more powerful x-ray tube providing higher flux, and a detector featuring a higher frame rate, detection efficiency and noise performance. Upgrades in this sense are currently underway in our lab through the installation of a W target (instead of Mo) Rigaku Source, combined with the Pixirad single photon counting, direct conversion detector [7,8].
With the current experimental setup, the continuous rotation reverse projection scan of the cylinder phantom reported above took 25 min, but we envisage that such an upgraded system will allow approximately halving the acquisition times, although it remains to be seen how the energy spectrum of the W source affects image quality in comparison to that of the Mo source.
As presented in this manuscript, reverse projection retrieval only holds for parallel or almost parallel beams, and therefore it is only applicable to samples of a few mm when an x-ray tube is used. An extension to fan beams has been suggested for grating-based x-ray phase contrast CT in [19], which could possibly be extended to EI-based x-ray phase contrast CT as part of future work. It should also be noted that reverse projection retrieval in combination with continuous sample rotation is not compatible with dithering (the procedure to artificially increase the sampling rate of the projections), although new sampling schemes will be explored in the future to allow dithering to be included. For example, the sample could be rotated over 720° with a single, half-pixel lateral step after the first 360°, and so forth for 'denser' spatial sampling. While this could be considered a downside when imaging samples that contain very small details, for many 'real life' applications the non-dithered sampling rate may be sufficiently high. For example, in the experimental setup described above, the nondithered sampling rate is 1 data point per 79 μm. A second drawback of the approach is that a continuous rotation introduces a degree of blur to the tomograms. In fact, the reconstructed tomograms show the object functions (δ, β k ) with a 'spin blur' caused by the rotation of the sample by the angular increment θ ∆ . However, this can be reduced by changing this parameter: the smaller θ ∆ , the smaller the blur. Another option to quantify the introduced blur is to adapt the concept of 'dynamic bias', which is often used within the fluid-flow community [17].