Improving radiation dose efficiency of X-ray differential phase contrast imaging using an energy-resolving grating interferometer and a novel rank constraint.

In this paper, a novel method was developed to improve the radiation dose efficiency, viz., contrast to noise ratio normalized by dose (CNRD), of the grating-based X-ray differential phase contrast (DPC) imaging system that is integrated with an energy-resolving photon counting detector. The method exploits the low-dimensionality of the spatial-spectral DPC image matrix acquired from different energy windows. A low rank approximation of the spatial-spectral image matrix was developed to reduce image noise while retaining the DPC signal accuracy for every energy window. Numerical simulations and experimental phantom studies have been performed to validate the proposed method by showing noise reduction and CNRD improvement for each energy window.


Introduction
In recent years, grating interferometer based x-ray differential phase contrast (DPC) imaging systems [1][2][3] have been extensively investigated with the hope of finding their niche in medical applications [4][5][6][7][8][9]. Combining a Talbot-Lau interferometer and a conventional x-ray imaging setup with a medical grade tube-detector assembly may put less stringent requirements on the x-ray source and detector to facilitate the translation into medical applications. However, one challenge to be overcome is to improve the radiation dose efficiency, e.g. measured by the dose normalized contrast to noise ratio (CNRD), of the DPC imaging systems with the Talbot-Lau interferometer. For a quantum limited DPC imaging system, the noise level of a DPC image is dependent on both the detected photon flux and the fringe visibility of the Talbot-Lau interferometer. The fringe visibility in Talbot-Lau interferometers has a strong dependency on energy.
If the x-ray beam energy differs greatly from the energy where fringe visibility is maximized for a given interferometer (this energy is referred to as the designed energy throughout this paper), the measured visibility of the system drops considerably. For a given interferometer system with a designed energy for maximal fringe visibility, the maximal fringe visibility cannot be achieved when information from x-rays with a wide energy range is integrated at the detector. In fact, although x-ray photons with energies deviated from the designed value do carry phase information about the image object, due to the reduction of fringe visibility at these energies, they also greatly add noise to the measured DPC signals and thus significantly degrade the radiation dose efficiency of the DPC imaging systems provided that an energy integration detector (EID) is used. In this regard, an energy-resolving photon counting detector (PCD) does have potential benefits, as it enables a Talbot-Lau interferometer to be implemented as an energyresolving grating interferometer so that one can selectively use x-ray photons to achieve the highest fringe visibility in DPC imaging.
With an energy-resolved interferometer, one can adjust the center of the energy window to match with the designed energy of the interferometer for the best fringe visibility. The next question is how to select an optimal energy window width to obtain the best DPC imaging performance. If the energy bin width is too wide, then the overall fringe visibility will drop due to the averaging effect over different beam energies within the energy window. However, if the energy bin width is too narrow, the number of photons detected within the energy window drops, leading to an increase in DPC image noise and thus a decrease in radiation dose efficiency. Therefore, an optimal energy bin width should be selected to balance the drops in fringe visibility and the number of detected photons within the selected energy bin [10,11]. While it is straightforward to search for an optimal energy window when the system is operated at the first Talbot order, the definition of an optimal energy window may become tricky when the system is operated at higher Talbot orders [12,13] due to the existence of multiple local maxima and minima in visibility-energy response curve.
Alternatively, when multiple energy thresholds are available for use in PCD, one can simultaneously acquire DPC images in multiple energy windows and then combine the acquired images from different energy windows with a weight, so that the signal to noise ratio in the final combined DPC image can be maximized [14]. The optimal weight in each energy bin has been found to be inversely proportional to the square of the energy and to the noise variance of the DPC signal in the corresponding energy window. A potential challenge to this method is the need to accurately calibrate the photon energy at each energy bin. It has been demonstrated that smaller weights should be assigned to both high and low energy bins that deviate from the central energy bin with highest fringe visibility.
In this paper, a novel method was developed to reduce noise of DPC images acquired using an energy-resolving interferometer system regardless of energy window partitioning schemes. This method utilizes information from all energy bins, and the method depends on neither the image object itself, nor the number of energy bins used in experiments, nor the need of accurate energy calibration at each energy bin. Numerical simulations and experimental phantom studies have been performed to validate this method.

Fundamental physics in DPC imaging
The Talbot-Lau interferometer based x-ray DPC imaging technique [3] measures the differential phase information φ (x, y, E) of an object at a given x-ray beam energy E. In theory, the measured DPC signal from a Talbot-Lau grating system is defined as: where C is a constant, E is the x-ray beam energy, and ρ e (x, y, z) represents the electron density distribution of the object. Herein, the derivative is taken along the y-axis, which is assumed to be perpendicular to the grooves in the gratings. In reality, the measurement is performed over a finite energy window. In this case, the beam energy E in Eq. (1) should be interpreted as the effective energy E e f f of the used energy window.
To facilitate the discussion, it is assumed that an ideal photon counting detector is used in the DPC data acquisitions, and a phase-stepping method [2] is used to acquire the data. The detected signal from a given detector pixel (p, q) (where p, q are the pixel indices corresponding to the horizontal and vertical directions on the detector plane) at the k-th phase stepping position can be modeled as [2,3,15] where k = 1, 2, . . . , M. The total number of phase steps is M, and N k (p, q, E) is the measured photon number at energy E for the k-th phase stepping position. The three energy-specific parameters, N 0 , N 1 , and φ , which correspond to three contrast mechanisms, i.e., attenuation contrast, dark-field contrast, and differential phase contrast, respectively, can be estimated from the M measurements. For example, one of the most commonly used approach to retrieve phase signal from detector counts is given by Therefore, phase retrieval is generally not a linear combination of detector counts. The intrinsic quantum fluctuations of the measured photon number N (k) will be propagated into the three estimated quantities of interest (N 0 , N 1 , and φ ). In particular, theoretical studies [16][17][18][19][20][21][22] have demonstrated that the noise variance of the extracted differential phase signal, i.e., σ 2 φ (E) , can be written as where ε(E) = N 1 (E)/N 0 (E) is defined as the fringe visibility of the interferometer.

Spatial-spectral DPC image matrix and the nature of its low dimensionality
If a DPC image dataset is acquired using an energy-resolving photon counting detector with Λ energy bins, the DPC image in the i-th energy window can be indexed as φ (p, q, To simplify the presentation, each two-dimensional DPC image φ (p, q, E i ) is vectorized using the lexicographic order as follows: where P and Q are the largest index values along the horizontal and vertical directions of the detector plane, respectively. Using this vector notation, DPC images acquired from all energy windows can form a spatial-spectral image matrix J as follows: Note that this spatial-spectral image matrix J has P × Q rows and Λ columns. Each row of J can be regarded as the change in DPC signal for each detector pixel from one energy window to the other.
A key observation in this paper is that, the i-th and j-th columns in the spatial-spectral DPC image matrix J only differ from each other by a global scaling factor of E 2 j /E 2 i [see Eq. (1)], provided that there is no noise or other non-idealities in DPC signal. Subsequently, each column is linearly depended on all other columns, and thus the rank of the matrix J is 1. This observation is independent of the number of energy windows used in DPC imaging. When noise and other non-ideal factors are present in the measured data, the rank of matrix J will be elevated. Therefore, the intrinsic rank-one nature of the spatial-spectral image matrix can be utilized to reject the contributions from noise and other non-ideal confounding factors.

Augmented spatial-spectral matrix J P
When a total of Λ energy bins are used, additional DPC images can be formed by combining the detected photons from two or more energy bins. The DPC images generated in this way are referred to as prior images in this paper, which can be used to augment the spatial-spectral image matrix J to form a new spatial-spectral image matrix J P . If a total of λ prior DPC images are formed, an augmented matrix J P with Λ + λ columns can be written as follows: All the prior DPC images φ (p, q, E p l ) added in this paper were generated as follows: a direct summation of the detected photon numbers of different energy windows was performed before Eq. (3) was used to extract the corresponding DPC image. Due to the nonlinear nature of phase retrieval as shown in Eq. (3), an augmented prior DPC image vector φ (p, q, E p l ) is not equivalent to a linear combination of the first Λ columns in J P . Clearly, the rank-one nature of the spatial-spectral image matrix J is always inherited by the augmented spatial-spectral matrix J P provided that there is no noise or other confounding factors. However, due to the presence of noise and other confounding factors in the measured data, the augmentation of these prior images to the original spatial-spectral image matrix J can be found to be beneficial in noise reduction [23].

Rank-one approximation of the augmented spatial-spectral matrix J P
Due to the presence of noise and other non-ideal confounding factors in DPC imaging, the rank of an experimentally measured spatial-spectral DPC image matrix J P may deviate from the ideal value of 1 (just as it has for matrix J). Therefore, the intrinsic rank-one nature can be utilized to generate a low rank approximated spatial-spectral matrix J * P with reduced image noise levels. The rank of a matrix can be conveniently determined using the well-known singular value decomposition (SVD) method [24], which is a generalized spectral decomposition of the matrix: In Eq. (8), U and V are two orthonormal matrices, and Σ is a diagonal matrix with entries σ 1 ≥ σ 2 ≥ · · · ≥ σ Λ+λ . After decomposition, each column in U only contains the spatial information of the object, while energy dependent information appears only in V T . The singular values σ i provide the weighting factors for the corresponding spatial-spectral basis u i ⊗ v T i . Following singular value decomposition, the rank-one nature of the augmented spatialspectral image matrix J P can be exploited to generate a noise reduced spatial-spectral image matrix J * P as follows: In total, matrix J * P contains Λ + λ columns. Each column restores the denoised DPC information for each corresponding energy window. It will be demonstrated that, with this rank-one approximation, noise is reduced and the accuracy of the DPC signal for each energy window is maintained.

Methods to generate prior DPC image vectors
As an example, Fig. 1(a) shows a PCD with two energy thresholds: a low energy threshold E 1 T H and a high energy threshold represent the number of photons detected above the corresponding energy threshold at k-th phase stepping position. Hence, two energy bins can be generated: Bin 1 and Bin 2. Assuming the effective mean energy of Bin 1 is E 1 , then the number of detected photons in Bin 1 is given by Similarly, the detected photon number in Bin 2 is T H ). The DPC signals in each energy bin, φ (p, q, E 1 ) and φ (p, q, E 2 ), can be estimated using Eq. (3). For a PCD equipped with three or more energy thresholds, the detected photons in each energy bin can be measured using the same method.
In this paper, only one prior DPC image φ (p, q, E p 1 ) has been added and examined for different energy threshold settings. It is retrieved directly from the detected photons above the lowest energy threshold E 1 T H .

Numerical simulation studies
To motivate and validate the use of the rank-one assumption in this paper, wave-optical simulations were performed in the first place. The technical details of the numerical simulation framework was presented in [25]. In this work, a 40 kVp spectrum was simulated with an energy step size of 0.1 keV. The simulated Talbot-Lau interferometer system consists of three gratings. The parameters for each grating were given as follows:

Experimental studies
To validate the proposed rank-one approximation method, the spatial-spectral DPC image datasets were acquired on a benchtop system equipped with a grating interferometer and a photon counting detector. The interferometer was operated at the first order Talbot distance, and it consists of a silicon based phase grating G1 and an absorption grating G2. Both gratings have a duty cycle of 0.5. By design, the G1 grating has a period of 8.0 μm and induces a π-phase shift for mono-energetic x-ray photons of 28 keV. The G2 grating has a period of 4.8 μm and the x-ray absorber is made from 60.0 μm thick gold. In addition, the distance between the x-ray focal spot and the G1 grating is 108.4 cm, and the G2 grating is placed at 21.7 cm downstream the G1 grating. The x-ray source used in this work is a micro-focus x-ray tube (model: L10321, Hamamatsu Photonics K.K., Japan) with a tungsten target. It was operated at 40 kVp (mean energy is about 28 keV) and 150 μA, and the nominal focal spot size is 8 μm.
In the experimental data acquisitions, the full x-ray spectrum was equally divided into different energy bins: from 2 bins up to 7 bins. To study the denoising performance of the proposed rankone approximation method at different dose levels, 30 repeated frames were acquired at each phase stepping position during a total exposure time of 20 seconds. The interferometer system is equipped with an energy-resolved single photon counting detector (model: XC-FLITE X1, XCounter AB, Sweden). The x-ray sensor is made of CdTe and the detector has two adjustable energy thresholds. Its native pixel size is 100 μm, and a 2 × 2 binning mode was used in the data acquisitions. Two different physical phantoms were used to evaluate the performance of the proposed rank-one method. The first phantom is a polymethylmethacrylate (PMMA) tube with an outer diameter of 9.6 mm and inner diameter of 8.1 mm, see Fig. 2(a). The phantom contains five groups of spheres submerged in a vegetable oil bath. They are ordered from top to bottom as follows: Acrylic ball with a diameter of 6.3 mm; Polyoxymethylene (POM) balls with a diameter of 6.3 mm; and four POM balls with a diameter of 3.2 mm. The other phantom is a PMMA wedge phantom with a slope of 10 degrees, see Fig. 2(b). The constant refraction angle generated by this phantom enabled the contrast-to-noise ratio (CNR) to be quantified for DPC imaging.

Numerical simulation results
Numerical simulation results are presented in this section to demonstrate the rank-one nature of the noiseless DPC signal and to validate the proposed low rank approximation method to reduce noise and thus improve dose efficiency. 4.1.1. Demonstration of rank-one nature of the DPC signal Figure 3 shows the singular value distribution of the spatial-spectral DPC image matrix J at several different number of energy bins when the quantum noise is absent and present. As can be observed, when noise is absent, only the first singular value is dominant, which demonstrated the rank-one nature of the DPC signal. When noise is present, smaller singular values may be elevated and the nature of the proposed rank-one approximation is to discard the contributions from these smaller singular values since the corresponding basis image contain primarily noise.

Validation of rank-one approximation in noise reduction
When quantum noise was introduced, the proposed rank-one approximation was able to generate noise reduced DPC images while well preserve the image information content, as shown in Fig. 4 for three energy bins. The results for more energy bins showed the same trend and thus were not presented here.

Measured beam spectrum and fringe visibility of the grating interferometer
By adjusting the energy thresholds of the PCD with a step size of 3 keV, both the x-ray transmission spectrum at the detector and the fringe visibility ε(E) of the used interferometer were experimentally measured. Results shown in Fig. 5 verify that both the spectrum of x-rays incident on the detector plane and the fringe visibility peaked near the designed operation energy (28 keV) of the grating interferometer.

Singular values and spatial basis of the augmented spatial-spectral image matrix
To validate the assumption that the DPC information is primarily contained in the first singular value and the corresponding basis, while higher order singular values and bases contribute primarily to noise, the singular values and the corresponding spatial bases for the PMMA tube phantom are presented in Fig. 6. As clearly shown in the results, the first singular value σ 1 is much larger than the second and third two singular values. Similarly, only the first spatial basis u 1 contains meaningful object information. In contrast, the other two spatial bases, u 2 and u 3 , contain negligible information of the image object. These experimental results validated the rank-one approximation and the feasibility of noise reduction by neglecting higher order singular values and spatial bases.

Experimental DPC images
DPC images generated without the proposed rank-one approximation for both phantoms were shown in the top row of Fig. 7. Clearly, the image noise levels varied from one energy bin to another. Images generated using the proposed rank-one approximation method are presented in the second row. Noise in these images was significantly reduced. Difference images between the "Original" images and the "SVD:J P " images are also presented in the bottom row, and there is no noticeable loss of structural information associated with the proposed rank-one approximation method. The line profiles of the tube phantom in Fig. 8 demonstrate that the spatial edge profile of the phantom is well preserved by the proposed rank-one approximation method.

Quantitative assessment of signal value and contrast-to-noise ratio
Quantitative measurements were performed with images of the wedge phantom to evaluate the performance of the proposed rank-one approximation method. Mean DPC signal values and standard deviations were measured in three rectangular regions (A, B, and C) in Fig. 7, and the results are tabulated in Table 1. The CNRs of DPC images are presented in Fig. 9.
There are several observations based on these results: (1) the proposed rank-one approximation method is able to maintain the signal accuracy (less than ±5% on average); (2) the proposed rank-one approximation method reduces noises for all energy bins, and thus is able to improve CNR; (3) the image CNRs cross all energy bins become identical after the rank-one approximation; (4) the augmented rank-one approximation method using matrix J P can further reduce image noise and thus improve image CNRs. Figure 10 shows CNR values of DPC images (region B) acquired at different radiation exposure levels, and three processing methods were compared: The first method, denoted as "Original-Max", used only a single energy bin that achieved the maximum CNR across different energy bins. The second method, denoted as "SVD:J", used the rank-one approximation and matrix J without the augmented image vectors. The third method, referred to as "SVD:J P ", used the augmented matrix J P and the proposed rank-one approximation.

Dependence on radiation exposure levels
As shown in Fig. 10, all three methods demonstrated linear relationships between the meas-  ured squared CNR and radiation exposure level, and the proposed rank-one approximation method improved CNR at all exposure levels involved in this study. Note that the CNR improvements of the DPC images processed by the proposed rank-one approximation methods might be limited at very low exposure levels. Figure 11 plots the measured CNR value as a function of total number of energy bins. For a given radiation dose level, the CNR of DPC images decreased with increasing total number of energy bins for all three methods, as the number of photons in each bin decreased. For a given   total number of energy bins and given exposure level, the proposed rank-one approximation method with the augmented matrix always generated the highest CNR (and therefore the highest radiation dose efficiency).

Discussion and summary
With an energy-resolving photon counting detector, a novel rank-constrained method was developed to reduce DPC image noise without sacrificing the DPC signal accuracy. Through both the numerical simulations and experimental phantom studies, it was demonstrated that the CNR of DPC images can be improved by the proposed method, therefore the radiation dose efficiency of the DPC imaging systems can be potentially improved. It was demonstrated that the rank-one approximation method always leads to matched CNR across all energy windows. This is a direct result of the rank-one constraint. Because of the rank-one nature of the spatial-spectral image matrix J and J P , a common spatial basis, u 1 is shared by images in different energy bins after the rank-one approximation, and this spatial basis is scaled by a numerical scaling factor determined by the effective energy of each energy bin to generate the corresponding DPC image for that energy bin.
The interferometer system used in this study was operated at the first fractional Talbot order, at which the fringe visibility of the diffraction fringe pattern drops gradually when the beam energy deviates from the designed operation energy of the interferometer. At higher fractional Talbot orders, the fringe visibility may reach zero at multiple energy levels [12,13]. In practice, however, the width of an energy bin width is always finite, and the energy resolution of the PCD is always limited, As a result, these theoretical zero-crossings of the visibility-energy curve are usually smeared out with polychromatic x-ray radiations. Due to this reason, we expect the proposed method is still applicable to high order Talbot interferometer setup, although it needs to be experimentally validated in future work.
In summary, a method that exploits the intrinsic rank-one nature of x-ray differential phase contrast images generated by different energy bins of an energy-resolving PCD was developed. Both numerical simulations and experimental validation studies were performed, and it was demonstrated that the proposed method can improve the contrast-to-noise ratio and radiation dose efficiency of differential phase contrast imaging systems.