Direct Fan-Beam Reconstruction Algorithm via Filtered Backprojection for Differential Phase-Contrast Computed Tomography

Recently, a novel data acquisition method has been proposed and experimentally implemented for differential phase-contrast computed tomography (DPC-CT), in which a conventional X-ray tube and a Talbot-Lau-type interferometer were utilized in data acquisition. The divergent nature of the data acquisition system requires a divergent-beam image reconstruction algorithm for DPC-CT. This paper focuses on addressing this image reconstruction issue. We have developed a filtered backprojection algorithm to directly reconstruct the DPC-CT images from acquired fan-beam data. The developed algorithm allows one to directly reconstruct the decrement of the real part of the refractive index from the measured data. In order to accurately reconstruct an image, the data need to be acquired over an angular range of at least 180◦ plus the fan angle. As opposed to the parallel beam data acquisition and reconstruction methods, a 180◦-rotation angle for the data acquisition system does not provide sufficient data for an accurate reconstruction of the entire field of view. Numerical simulations have been conducted to validate the image reconstruction algorithm.


Introduction
In recent years, X-ray phase-contrast imaging methods have attracted significant interest.Based upon the phase retrieval method, phase-contrast imaging can be classified into three categories: interferometric methods [1,2], analyzer crystal methods [3,4], and free-space propagation methods [5][6][7].These methods generally require excellent spatial coherence and/or excellent spatial resolution of the detectors.Phasecontrast imaging is thus normally conducted using highbrilliance synchrotron beam lines or microfocus X-ray tubes.It is difficult to adapt these phase retrieval methods to a clinical setting for medical imaging.
The shearing interferometry method only requires a mild spatial coherence length, which opens up a new opportunity to implement phase-contrast imaging using a conventional X-ray tube.A method which takes advantage of these characteristics was proposed and experimentally implement in Xray differential phase-contrast imaging [17] and differential phase-contrast computed tomography (DPC-CT) [18,19].In this scheme, Pfeiffer et al. [17][18][19] proposed to mount an absorption grating in front of a low-brilliance X-ray source to generate an array of equidistant secondary line sources.The coherence length of the beam is determined not by the focal spot size of the X-ray tube, but rather by the width of the opening of the absorption grating, which is about 20 ∼ 50 μm.This allows the use of a conventional X-ray tube, where the high X-ray output can be turned into many mutually incoherent line sources.In this paper, we refer to this grating as a beam slicer.Using this method, the spatial coherence length of each line source has been improved by a factor of ten or more over that of the X-ray tube without the grating, thus enabling significant reduction of the distance between the X-ray source and the image object.As a consequence, this method enables construction of an X-ray phase-contrast imaging system with compact size, which is often a requirement in medical and nondestructive inspection applications.
For a given X-ray detector with rectangular geometry, the fan angle is determined by the angle given by the following formula: where H is the width of the detector and D is the distance from the X-ray focal spot to the detector.Therefore, when the source to detector distance D is very long, the divergence of the X-ray beam is negligible, which is often the case for synchrotron beam line experiments and inline holography phase-contrast imaging.In this case, the beams are well approximated as plane wave, and a parallel-beam image reconstruction method is sufficient for image reconstruction in DPC-CT.However, using the new phase retrieval method by Pfeiffer et al. [17][18][19], the distance D is significantly reduced (D 160 cm in [19]).Note that the effective detector size is limited by the size of the third grating in front of the detector, about 6 cm in [19].Using these parameters, the fan angle γ m in [19] can be estimated as In this case, the parallel-beam approximation is still acceptable (see Section 4.2).Thus, the well-known parallel-beam image reconstruction method [20] can be adapted to reconstruct DPC-CT images [18,21,22].When a larger image object is scanned with the same detector-to-source distance D, the fan angle γ m is increased and the parallel-beam approximation becomes less accurate.When the fan angle is increases beyond 5 • , the divergent nature of the data acquisition system is no longer negligible (see Section 4.2).For example, if the size of the image object is doubled in the above estimation, the fan angle is nearly 5 • .For larger objects, a divergent beam image reconstruction algorithm would be desirable for this new DPC-CT data acquisition method.In this paper, we present an image reconstruction formula for fan-beam DPC-CT which can be utilized to accurately reconstruct DPC-CT images at any fan angle provided that the data is acquired within the angular range of 180 • plus fan angle.

Derivation of the Fan-Beam Image Reconstruction algorithm
In general, the experimental results should be explained by wave optics.However, for convenience, an effective geometrical model is utilized in this paper to develop the image reconstruction algorithm.

Brief Review of the Parallel-Beam Reconstruction of the Differential Phase-Contrast X-Ray CT
As explained by Pfeiffer et al. [18], the fundamental idea in differential phase retrieval is to locally measure the angular change Θ(ρ, θ) of the incident X-ray beam (Figure 1).Here, the measured data is labeled by the distance ρ, which is measured from the origin of the coordinate system to the incident ray direction and the angle θ, which specifies the direction of the incident ray.
The quantity Θ(ρ, θ) is related to the local gradient of the object's phase shift can be written as follows [15,23]: where δ(x, y) is the decrement of the real part of the object's refractive index, n = 1 − δ + iβ.The letter l labels the incident ray; therefore, the above line integral is performed along the incident ray direction.Note that the integral in ( 3) is simply the Radon transform R δ (ρ, θ) of the target function δ(x, y).Thus, (3) can be rewritten as After a Fourier transform with respect to the variable ρ, (4) is recast into the following form: Using the inverse Radon transform in Fourier space, one can easily obtain the following image reconstruction formula [19,21,22]: Therefore, after the data are acquired from view angles in the angular range [0, π], a one-dimensional Fourier transform is performed with respect to the Radon distance ρ to obtain Θ(ω, θ).A filtering kernel iπsgn(ω) is used to filter the data.Finally, an integral over the view angles is performed to reconstruct images.

Effective Fan-Beam Data Acquisition for Differential Phase Retrieval Method
Note that the image reconstruction formula (6), used by Pfeiffer et al. [19], makes sense only when the acquired data are written in terms of the Radon distance ρ and the view angle θ.Namely, it is only convenient for the parallel-beam data acquisition geometry.
For the effective fan-beam data acquisition geometry as shown in Figure 2, the X-ray tube and detector are assumed to rotate around the image object along a circle with radius R and source-to-detector distance D. The acquired data are naturally labeled by the angular position, t, of the X-ray focal spot and the detector position which is labeled by a projection distance u.We refer to this measured data as Θ(t, u) in this paper.The half width of the detector is u m , which implies that the fan angle γ m = 2tan −1 (u m /D).
For fan-beam data, one way to reconstruct an image is to rebin the acquired data into parallel-beam projection data and then utilize (7) to reconstruct the image.In this paper, we develop new image reconstruction formulae which directly reconstruct images from the acquired data Θ(t, u) without requiring the rebinning step.

Fan-Beam Reconstruction Formula for
Full Scan Case In order to directly reconstruct the image without the rebinning step, we rewrite (6) in the real domain by taking the inverse Fourier transform: This formula can simply be derived by noting that the product in Fourier space is a convolution in real domain.This formula dictates how each individual projection datum Θ(ρ, θ) contributes to the final image.Note that we have extended the integral limit from [0, π] to [0, 2π] and introduced an extra factor of 1/2 to account for the redundancy.We next derive image reconstruction formulae for the fan-beam data acquisition geometry.Each individual projection datum Θ(t, u) can also be relabeled in terms of the Radon distance ρ and angle θ by the following nonlinear coordinate transforms: This transform introduces an extra Jacobian factor: where the functions L(x, y; t) and U(x, y; t) are given by Therefore, (7) can now be written as where the filtered function F[t, U(x, y)] is given by In this equation, the preweighted projection data Θ(t, u) is defined as Equations ( 11)-( 13) are the main results of the paper.They give a direct image reconstruction from the measured projection data Θ(t, u).The algorithm for fan-beam DPC-CT is summarized in the following steps: (1) preweight the acquired projection data using ( 13); (2) filter the preweighted data using the Hilbert filtering kernel in ( 12); (3) back-project the filtered data to reconstruct an image using (11).
Note that these formulae require that the data are acquired along a full circle.This is referred to as the full scan data acquisition mode.In the next subsection, we extend the above result to the short scan case, where the data are acquired from less than a full circle.

Fan-Beam Reconstruction Formula for Short Scan Case
In absorption X-ray CT, it has been proven that the data are sufficient to reconstruct the entire image object when the angular range is greater than 180 • + fan angle, where the fan angle is defined as the entire angle subtended by the detector from the source position.The same condition is also true for differential phase-contrast CT, we refer to this case as the short scan mode.When the short scan data acquisition mode is utilized, some incident rays will have a conjugate ray, while others will be measured only once.Therefore, a weighting function is needed for the redundantly measured data.It is not difficult to convince oneself that the following two data points should be considered as the same: Thus, similar to absorption CT [24], the following weighting function can be utilized to weight the data before the filtering step: where the variable γ is defined as γ = arctan(u/D).Thus, by incorporating the normalized weighting function W(t, u) into ( 13), the preweighted projection data become The image reconstruction formulae ( 11)-( 13) together with ( 16) can be used to reconstruct the image.

Numerical Simulations and Results
In order to validate both the full scan and short scan DPC-CT reconstruction formulae derived in this paper, numerical simulations were conducted using the fan-beam data acquisition geometry as shown in Figure 2. The mathematical phantom consisted of a uniform elliptical disk and two uniform circular disks.The projection data were numerically generated using Snell's law [23].The radius of the circular disks was 0.07 m.The semimajor axis and semiminor axis of the elliptical disk were 0.35 m and 0.175 m, respectively.
The decrement of the real part of the refractive index was assumed to be 1.0 × 10 −6 and 0.5 × 10 −6 for the circular and elliptical disks, respectively.The radius of the source trajectory R was assumed to be 1.4 m.The distance from the source to the detector D was selected to be 2.1 m.The sampling rate of the view angle was Δt = 0.5 • .The detector was assumed to be collinear, with total width of 1.13 m.The detector sampling rate Δu was 1.13 m/600 ≈ 2 mm.The image matrices were 256 × 256.

Full Scan Data Acquisition Mode
A fan angle of 30 • was used in order to cover the entire image object.The view angle sampling range was [0 • , 360 • ].The results are presented in Figure 3. From Figure 3(a) and the corresponding image intensity plots, it is evident that the image was accurately reconstructed using the new algorithm.

Short Scan Data Acquisition Mode
In the short scan mode, all parameters stay the same as the full scan mode with the following exceptions: (1) the decrement of the real part of the refractive index of the circular objects is changed to 0, and (2) the angular range of the data acquisition is reduced from 360 • to 210 • , which is sufficient for reconstruction according to X-ray absorption CT image reconstruction theory [20].The numerical results are presented in Figure 4. Again, the image is faithfully reconstructed with very high accuracy as indicated in image intensity plots.

Super Short Scan Data Acquisition Mode
In this subsection, we numerically demonstrate that the angular range for data acquisition should be at least π + fan angle for the presented new algorithm.When the scanning angular range is shorter than that, significant image artifacts may occur.In numerical simulations, all of the parameters are the same as the short scan case except that the source angle is reduced to from 210 • to 180 • , which is less than a short scan.As shown in Figure 5, image artifacts appear in the reconstructed image, and signal drop is clearly observed in image intensity plots.This numerical experiment also indicates that the short scan angular range is necessary for the presented fan-beam DPC-CT reconstruction.

Discussion
In this section, several aspects of the newly developed image reconstruction algorithm will be discussed.First, fan-beam DPC-CT image reconstruction and fan-beam absorption CT image reconstruction algorithms are compared.Next, the differences between the fan-beam and known parallel-beam reconstruction algorithms are discussed.Finally, the angular range requirements for fan-beam DPC-CT data acquisition are examined.

Comparison with Fan-Beam Absorption CT Image Reconstruction Algorithm
In this paper, the fan-beam image reconstruction algorithm was derived using a collinear detector geometry.
It is easy to change to the curved detector case using the relation γ = arctan(u/D).For absorption CT, fanbeam image reconstruction algorithms have been derived for both full and short scan cases [20].The DPC-CT fanbeam image reconstruction algorithm is different from its absorption counterpart in three ways.First, in DPC-CT, the preweighting factor is D 2 /(D 2 + u 2 ).For absorption CT, the preweighting factor is a square root, that is, Second, the distance weighting function in the backprojection step in DPC-CT is R/L(x, y; t).While for absorption CT, the weighting function is R 2 /L 2 (x, y; t), the square of the DPC-CT expression.Finally, the filtering kernel for DPC-CT image reconstruction is a Hilbert filter, whereas a ramp filter is used in absorption CT image reconstruction.

Parallel-Beam Approximation
When the source to detector distance is much larger than the size of the image object, one can expect that the difference between the parallel-beam and fan-beam image reconstruction methods is negligible.The numerical results demonstrate that when the fan angle is smaller than 5 • , the parallel-beam approximation is valid, and parallel-beam image reconstruction can be utilized to reconstruct DPC-CT images.However, for larger fan angles, a fan-beam image reconstruction algorithm, such as the one presented in this paper, is needed for an accurate image reconstruction.
In the first numerical experiment, we applied the parallel-beam reconstruction algorithm [19,21,22] to the case of a 30 • fan angle, the same angle used in the reconstruction presented in Figure 3.The parallel-beam reconstruction algorithm generates significant image artifacts as shown in Figure 6(a).In the second numerical experiment, the radius of the source trajectory is increased to 4 m and the fan angle of the divergent beam is reduced to 10 • in order to cover the entire image object.
As shown in Figure 6(b), the image quality is improved comparing with Figure 6(a), but blurring still exists at the boundaries.When the source to detector distance increases such that the fan angle is reduced to 5 • , as shown in Figure 6(c), the image artifacts become increasingly negligible, indicating that the parallel beam approximation is becoming more acceptable.In the last numerical experiment, the fan angle of the divergent beam is reduced to 2.5 • for the same image object.In this case, as shown in Figure 6(d), the parallel-beam reconstruction method is sufficient to reconstruct DPC-CT images.Using the experimental data presented in [19], the fan angle can be estimated to be less than 2.5 • .Thus, parallel-beam reconstruction would be sufficient to reconstruct near artifact-free DPC-CT images [19].

Conclusions
In this paper, an image reconstruction formula for fanbeam DPC-CT was derived and validated using numerical simulations.It has been demonstrated that when the image object is relatively large, the fan angle must increase to cover the entire image object, and a parallel-beam approximation cannot be directly applied to reconstruct images.In this case, the fan-beam image reconstruction formula presented in this paper can be directly applied to accurately reconstruct DPC-CT images.The presented image reconstruction formula is also useful for other DPC-CT imaging method [18] provided that the divergent beam data acquisition geometry is utilized.As a final remark, as shown in Figure 5, one limitation of the presented algorithm in this paper is that the data acquisition angular range must be longer than π + fan angle.Otherwise, significant image artifacts are observable in the reconstructed image.In order to enable accurate image reconstruction when the angular range is shorter than π + fan angle, some other new image reconstruction must be developed [25,26].Note that it is straightforward to extend the fan-beam reconstruction algorithms such as the one presented in this paper and others [25,26] to cone-beam data acquisition systems using the so-called FDK approximation [27], where the only change is to replace the fan-beam data by the projected cone-beam data, that is, the multiplication of cone-beam data and a cosine weighting factor.

Figure 1 :
Figure 1: A geometric model of the parallel-beam data acquisition model.

Figure 2 :
Figure 2: Geometrical model for fan-beam DPC-CT data acquisition.R is the radius of the source trajectory, and D is source-todetector distance.

Figure 3 :
Figure 3: Reconstruction results for a fan-beam full scan.(a): Reconstructed image.(b), (c), and (d) Image intensity plots corresponding to lines 1, 2, and 3, respectively.The solid line represents the reconstructed image, and the dashed line represents the ground truth.The abscissas of (b), (c), and (d) are normalized to 0.35 m.

Figure 4 :
Figure 4: Numerical results for short scan data acquisition mode.(a) Reconstructed image.(b), (c), and (d) Image intensity plots corresponding to lines 1, 2 and 3, respectively.The solid line represents the reconstructed image, and the dashed line represents the ground truth.The abscissas of (b), (c), and (d) are normalized to 0.35 m.

Figure 5 :
Figure 5: The reconstruction result for angular range less than a short scan.Figure 6(a) is the reconstructed image.Figures (b), (c), and (d) are profile plots corresponding to lines 1, 2, and 3, respectively, in Figure 6(a).The solid line represents the reconstructed image, and the dashed line represents the ground truth.The abscissas of (b), (c), and (d) are normalized to 0.35 m.