3D Algebraic Iterative Reconstruction for Cone-Beam X-Ray Differential Phase-Contrast Computed Tomography

Due to the potential of compact imaging systems with magnified spatial resolution and contrast, cone-beam x-ray differential phase-contrast computed tomography (DPC-CT) has attracted significant interest. The current proposed FDK reconstruction algorithm with the Hilbert imaginary filter will induce severe cone-beam artifacts when the cone-beam angle becomes large. In this paper, we propose an algebraic iterative reconstruction (AIR) method for cone-beam DPC-CT and report its experiment results. This approach considers the reconstruction process as the optimization of a discrete representation of the object function to satisfy a system of equations that describes the cone-beam DPC-CT imaging modality. Unlike the conventional iterative algorithms for absorption-based CT, it involves the derivative operation to the forward projections of the reconstructed intermediate image to take into account the differential nature of the DPC projections. This method is based on the algebraic reconstruction technique, reconstructs the image ray by ray, and is expected to provide better derivative estimates in iterations. This work comprises a numerical study of the algorithm and its experimental verification using a dataset measured with a three-grating interferometer and a mini-focus x-ray tube source. It is shown that the proposed method can reduce the cone-beam artifacts and performs better than FDK under large cone-beam angles. This algorithm is of interest for future cone-beam DPC-CT applications.

Most existing DPC-CT approaches are based on three kinds of scanning geometries, i.e., parallel-beam, fan-beam and cone-beam. Favored by the high imaging efficiency and the magnified spatial resolution, cone-beam DPC-CT has attracted significant interest [22][23][24]. It includes essentially two steps: (i) retrieval of the DPC projections and (ii) phase reconstruction. The first step can be accomplished by using a phase-stepping procedure [12,14,15], a reverse projection method [25], or a single-shot Fourier-based phase-extraction method [26]. The second step has so far been solved by using a FDK-type filtered back-projection (FBP) algorithm [22,23], which is the direct extension of the FDK algorithm commonly used in absorptionbased CT [27] by replacing the ramp filter with a Hilbert imaginary filter to take into account the differential nature of DPC projection data. The FDK algorithm is computationally fast and relatively easy to implement. However, as its applications to the absorption-contrast CT, it suffers from the cone-beam artifacts because the necessary and sufficient conditions for accurate CT reconstruction with cone-beam geometry fail [28][29][30][31]. There are two major cone-beam artifacts with this approach: loss of contrast intensity and geometrical deterioration along the vertical axis. Both artifacts increase rapidly with the distance to the mid-plane, when the conebeam angle becomes large.
In the field of conventional x-ray absorption-based CT, iterative reconstruction algorithms have a wide spectrum of proven advantages [32][33][34], including dose reduction [35,36], sparse sampling [37,38], and limited angle tomography [39,40]. However, there were few studies on iterative reconstruction for DPC-CT in the past. Recently three kinds of iterative reconstruction techniques for fan-beam DPC-CT were proposed, namely, statistical iterative algorithms [32,[41][42][43], simultaneous algebraic reconstruction technique [44] and algebraic iteration reconstruction [45]. Those studies demonstrated the feasibility of iterative reconstruction for DPC-CT. Although iterative reconstruction exhibits great advantages over FBP-type methods, there is currently no in-depth discussion about the iterative reconstruction technique for conebeam DPC-CT.
In this work, we have discretized the cone-beam DPC-CT imaging modality as a system of linear equations and developed a 3D algebraic iteration reconstruction (AIR) algorithm for it. This algorithm is based on the Kaczmarz method [46] or the later rediscovered algebraic reconstruction technique (ART) [47] and reconstructs the image ray by ray. Unlike the conventional iterative algorithms for absorption-based CT, it involves the derivative operation in the forward projections to take into account the differential nature of the cone-beam DPC-CT projections. The validity of the proposed AIR method has been demonstrated with numerical simulation and experimental study with a dataset measured at a setup based on a three-grating interferometer and lab-based x-ray mini-focus source. It is shown that the proposed method can reduce the cone-beam artifacts and performs better than FDK under large cone-beam angles.

Materials and Methods
A three-dimensional object can be described by a complex refractive index distribution n(x, y, z) = 1 − δ(x, y, z) + iβ(x, y,z), where x, y and z describe the coordinate system of the sample. In conventional absorption-based CT, the imaginary part β is measured by the attenuation of the x-rays transmitted through the specimen. In differential phase contrast imaging, one measures the effect of variations of the real part δ by evaluating the tiny refraction angles of x-rays induced by the specimen with a grating interferometer. Correspondingly, under a paraaxial approximation, a differential phase-contrast projection can be expressed by where x 0 and z 0 describe the coordinate system of the detector, θ the rotation view angle of the object and l(x 0 ,y 0 ,z 0 ) the incident ray. Fig. 1 depicts the scanning geometry of the cone-beam DPC-CT. A sample holder (not shown in the figure) rotates the object over 360 in a three-grating interferometer during data acquisition. At each view angle, a phase-stepping procedure is performed and provides the DPC projection by Fourier analysis of the recorded moiré fringes by the detector. We are to reconstruct the phase slice images δ from the DPC projection.

3D Algebraic Iterative Reconstruction Algorithm
To establish the proposed iterative method, we first need to discretize Eq. (1). The reconstructed image is described by a 3D matrix δ of size N image with independent elements δ j , j = 1, 2, Á Á Á, N image . δ x,y,z refers to one voxel for the discretized 3D image δ with the voxel index, where and x = 1,2,Á Á Á,W; y = 1, 2, Á Á Á, L; z = 1, 2, Á Á Á, H. Integers L, W and H are, respectively, the length, width and height of the 3D image, which has a total number of voxels (1) is defined as the forward projection and represented by a 3D matrix P with size N prj . Each projection is referred to as P i or P a,b,t , where and t = 1,2,Á Á Á,T;a = 1,2,Á Á Á,A;b = 1,2,Á Á Á,B. Integers A, B and T are, respectively, the columns and rows of the flat panel detector and the number of view angles. The total number of measurements is N prj = A×B×T.
With the above notations, let M ij be the weight for the contribution of j-th voxel to i-th ray in the phase-contrast projection R l(x 0 ,y 0 ,z 0 ) δdl, and M be the N prj ×N image matrix. Then we have the following discrete-to-discrete linear transform from the phase image δ to the phasecontrast projections P The weights M ij are computed by calculating the intersection length of the i-th ray through the j-th voxel. The phase-contrast projection along the i-th ray can be expressed by (1) are represented by a 3D matrix α with size N prj . Each DPC projection α i = α a,b,t at the location (a,b,t) can then be represented by utilizing finite differences of the phase-contrast projections P. In this work, we choose to approximate the DPC projections α by the phase-contrast projections P with the central finite difference. Specifically, for i = 1, Á Á Á, N prj , where the function floor(r) gives the largest integer less than r and the function mod(r,s) the remainder for the division r s . Now the imaging system, namely, Eq. (1) is discretized by Eqs. (4) and (5). The problem of DPC-CT is then to reconstruct δ in Eq. (4) from the measured cone-beam DPC projections α in Eq. (5). Obviously, unlike the absorption-based CT, DPC-CT involves the derivative operation in (5) in the forward projections. Although there are many iterative reconstruction methods for absorption-based CT (Please refer to [33,34,48,49] for a review on this topic.), they are not applicable for DPC-CT since the simultaneous or block-iterative schemes involve weighted sums of all or multiple ray projections. They could blur the reconstructed image and deteriorate the accuracy of the derivation operation in (5) and hence the spatial resolution.
Based on the above analysis, we developed the following 3D AIR algorithm for cone-beam DPC-CT Here, k and i are non-negative integer. δ[k,i] represents the reconstructed image at the i th inner iteration within the k th outer iteration. The ray-by-ray reconstruction is called the inner iteration. Once all rays are used, it will start another iteration, which is called the outer iteration. δ [1,0] is the initial guess of the reconstructed image, which is chosen to be the zero at each voxel. The parameter d is the relax coefficient and can be from 0 to 2 in theory from the convergence theory for ART [34,46,47]. M i is the projection matrix for the i-th ray and a 0 i is the numerical calculation DPC projection by the derivative operation to the forward phasecontrast projection. The proposed AIR method is based on the algebraic reconstruction technique [47]. It iterates ray-by-ray to update the image δ[k,i] after comparing the measured DPC projection α i and the numerical projection a 0 i until convergence. Because the derivative operation in Eq. (5) only needs adjacent rays as few as possible, it is expected to provide better derivative estimates in iterations. The overall strategy of this algorithm considers the reconstruction process as the optimization of a discrete representation of the object function to satisfy a system of equations that describes the cone-beam DPC-CT imaging modality, and iteratively arrives at a solution from its cone-beam DPC projection dataset. In the following we discuss the implementation procedure.

Implementation Procedure
Having specified the system matrix M and data ã, we can describe the implementation of the algorithm in Eqs. (4), (5), (6) and (7). Depicted in Fig. 2, it includes an inner iteration and an outer iteration. This process is labeled in the following way: the outer iteration number is labeled by k, and the inner iterations by i in their respective loops. The intermediate image during the iterations is noted δ[k,i], indicating the i th inner iteration within the k th outer iteration. The successive steps of the algorithm are thus: (ii) Inner iteration initialization.
(iv) Outer iteration convergence check: compute If difference < , stop; Otherwise, go to step (ii) to initialize next iteration loop and set Here the parameter controls the convergence of the algorithm so that the algorithm stops when there is no appreciable change.

Numerical study
To test the proposed method, numerical simulations were performed. The phantom is adopted from the well-known Defrise phantom [50], which is known as cone-beam killer. It consists of nine discs of diameter 55 mm, vertically separated at different heights of 0, 26.1, 52.4, 78.7 and 105.1 mm (height was measured from the axial mid-plane of whole phantom to the center of each disc). The thickness of each disc was set to 3.0 mm and was attributed refractive index 1 × 10 −6 . The source-to-object and object-to-detector distances both were set to 1000 mm. The detector has 512×128 pixels with size of each pixel 1.0 mm. The corresponding DPC sinogram was calculated using an algorithm based on the refraction angle equation [1]. The normalized root mean square error (NRMSE) is calculated to evaluated quantitatively the image quality. Fig. 3 shows the reconstruction results of the Defrise phantom using FDK and AIR algorithms. Figs. 3 (a) and (b) are the images of the axial mid-plane of the Defrise phantom reconstructed by FDK and AIR algorithms respectively, which correspond to the case with a conebeam angle of 0 . In this case FDK is equivalent to the standard 2D FBP algorithm and can provide exact reconstruction. Obviously in this case AIR method is as accurate as FDK. Fig. 3 (c) further supports this conclusion, which presents the profiles along the red solid line and the  Table 1. Figs. 3 (d) and (e) display the sagittal slices of the Defrise phantom reconstructed by FDK and AIR algorithms respectively, which have a maximum cone-beam angle of 6 . In these figures, the cone-beam artifacts are quite visible. In Fig. 3 (d), FDK seems to split every disc into two thinner parts with lower reconstructed values and the extent of the separation of two parts is determined by the cone-beam angle of that height. The profiles presented in Fig. 3 (f) also provides quantitative evidence for this observation. As we can see from Fig. 3 (f), this geometrical distortion becomes more noticeable when the cone-beam angle becomes bigger than 4.5  3D Algebraic Iterative Reconstruction for Cone-Beam X-Ray DPC-CT and the reconstructed value of the disc drastically drops to nearly half of the theoretical value. In contrast to this, much less artifacts are observed in the image reconstructed by AIR presented in Fig. 3 (e). They are quite consistent with the phantom where even the cone-beam angle reaches 4.5 . Actually the corresponding profile in Fig. 3 (f) demonstrates that the reconstructed values still retain 80% when the cone-beam angle reaches 6 . Furthermore, Listed in table 1, the NRMSE value for AIR sagittal slice is 0.3165 and 52% smaller than the one for FDK sagittal slice. These results support that AIR performs better than FDK when the cone-beam angle increases.

Experiments
The experimental dataset that was measured to test the reconstruction algorithms was recorded at a CT setup for differential phase-contrast imaging, based on a three-grating interferometer and a mini-focus x-ray tube source installed in a compact gantry at the Technische Universität München (Munich, Germany). The sample was a phantom consisting of a small plastic tube filled with balls of different sizes and different plastic materials. The experimental setup consisted of a tungsten x-ray source (RTW, MCBM 65B-50 W) with a focal spot size of approximately 50μm in diameter. The detector was a flat-panel sensor (Hamamatsu C9312SK-06) with 2496×2304 pixels with 50 μm square pixel size. The gratings of the Talbot-Lau interferometer were fabricated by the company microworks, with grating parameters optimized for a design energy of 23 keV. The heights and periods of the grating structures were: 35 μm and 10μm for the source grating G0, 40 μm and 3.2 μm for the phase grating G1, and 25 μm and 4.8 μm for the analyzer grating G2, respectively. The source grating was placed 31 mm from the x-ray tube. The distance between G0 and G1 was 300 mm, whereas G1 and G2 were 145 mm apart, corresponding to the first fractional Talbot distance. The sourceto-sample and sample-to-detector distances were approximately 270 mm and 200 mm, respectively. A central area of 560×1110 pixels of the detector formed the cone-beam DPC-CT geometrical configuration together with the interferometer and the x-ray generator. The axial cone-beam angle was about 3.4 and the sagittal cone-beam angle 6.8 . The x-ray source was operated at 40 kV and 743 μA. The experimental dataset was acquired by taking 600 projections over 360 rotation, with 6 phase steps per projection and 5 seconds exposure per phase step. The retrieved DPC projections from all the view angles provide the complete experimental dataset with a size of 560×1110×600 pixels and some of them are displayed in Fig. 4 (a).
To give an impression of the entire reconstruction from the experimental cone-beam DPC data, the 3D volume rendering and the 3D orthogonal views of the results of the sample reconstructed by AIR and FDK algorithms are firstly presented in Fig. 4. From these images, we can find that there exists some internal structure difference in the results, even though they seem to be almost the same.
For more clarity, the 2D slices in three typical planes (the axial mid-plane, the axial plane with a cone-beam angle of 2.5 and the sagittal mid-plane) are selected and presented in Fig. 5.

Conclusion and Discussion
In summary, we have established an algebraic iterative reconstruction (AIR) algorithm for cone-beam DPC-CT and demonstrated its validity and performance both numerically with simulation and experimentally with real data. This method is based on the well-known algebraic reconstruction technique and tailored to the cone-beam DPC-CT. It is shown that the proposed method can reduce the cone-beam artifacts and performs better than FDK under large cone-beam angles.
It is of practical significance to accelerate the AIR algorithm. In our simulation and experiment, the implementation is with Microsoft Visual C++6.0 in a laptop (LENOVO ThinkPad T400). The laptop is assembled with a Intel Core 2 Duo CPU P8700 at 2.53GHz and a RAM memory of 3 GB. For our experimental reconstruction, it takes about 12 hours to converge. However, this is not an issue since there exist algorithm and hardware (graphics processing unit) to accelerate the algebraic reconstruction technique [51,52]. The discretization can also be improved with other basis rather the voxel configuration in this work [42]. Future work should be done to improve the reconstruction speed of this algorithm with these techniques toward DPC-CT imaging applications.