A theoretical framework for comparing noise characteristics of spectral, differential phase-contrast and spectral differential phase-contrast X-ray imaging

Spectral and grating-based differential phase-contrast X-ray imaging are two emerging technologies that offer additional information compared with conventional attenuation-based X-ray imaging. In the case of spectral imaging, energy-resolved measurements allow the generation of material-specific images by exploiting differences in the energy-dependent attenuation. Differential phase-contrast imaging uses the phase shift that an X-ray wave exhibits when traversing an object as contrast generation mechanism. Recently, we have investigated the combination of these two imaging techniques (spectral differential phase-contrast imaging) and demonstrated potential advantages compared with spectral imaging. In this work, we present a noise analysis framework that allows the prediction of (co-) variances and noise power spectra for all three imaging methods. Moreover, the optimum acquisition parameters for a particular imaging task can be determined. We use this framework for a performance comparison of all three imaging methods. The comparison is focused on (projected) electron density images since they can be calculated with all three imaging methods. Our study shows that spectral differential phase-contrast imaging enables the calculation of electron density images with strongly reduced noise levels compared with the other two imaging methods for a large range of clinically relevant pixel sizes. In contrast to conventional differential phase-contrast imaging, there are no long-range noise correlations for spectral differential phase-contrast imaging. This means that excessive low frequency noise can be avoided. We confirm the analytical predictions by numerical simulations.


Introduction
Spectral and grating-based differential phase-contrast X-ray imaging are two emerging technologies that offer additional information compared to conventional attenuationbased X-ray imaging. In the case of spectral imaging, information about the energy-dependent attenuation of an object is obtained by acquiring measurements with two or more photon energy spectra. Material decomposition algorithms use the energy-resolved measurements to generate basis material images, which can provide information about the chemical composition of an object. Various dual-energy acquisition schemes and recent advances in photon-counting detector technology [1,2] have encouraged research towards medical applications of spectral X-ray imaging [3,4,5,6]. The ability of photon-counting detectors to acquire multiple energy-resolved measurements that are spatially and temporally registered has proven advantageous for multi-material decomposition [7,8,9]. Grating-based differential phase-contrast (DPC) imaging exploits an entirely different contrast generating mechanisms in addition to the conventional attenuation contrast: The phase-contrast image is obtained by indirectly measuring the (differential) phase shift that an X-ray wave exhibits when transversing an object [10,11]. This differential phase shift can be directly related to the projected electron density (PED) of the object. The dark-field image is generated by ultra-small angle X-ray scattering from microstructures [12]. It is connected to the real space autocorrelation function [13,14] and thus provides information about the internal structure or surface properties [15] of an object far below the detector resolution. Laboratory experiments have demonstrated that the phase-contrast image can provide a strongly improved contrast-to-noise ratio compared with attenuation-based imaging [16,17,18]. However, based on a theoretical noise analysis, doubts have been raised whether these improvements can be transferred to clinical computed tomography (CT) applications [19,20]. Experimental studies have demonstrated that similarly to spectral imaging, material-specific information can be extracted from phase-contrast measurements [21]. Recently, we have investigated the combination of grating-based DPC radiography and spectral radiography [22] and have developed a basis material decomposition algorithm that uses the spectral and the phase contrast information simultaneously. Our numerical experiments have demonstrated that spectral differential phase-contrast (SDPC) radiography yields quantitatively correct basis material images with strongly reduced noise levels compared with conventional spectral imaging. Taking the characteristics of all three aforementioned imaging methods into account, the question of the optimal method for a particular imaging task arises. In this work, we investigate one important aspect of this question by conducting an in-depth noise analysis of all three imaging methods. Towards this end, we develop a noise analysis framework based on the Cramér Rao lower bound (CRLB) [23]. This framework allows the determination of the optimum imaging parameters for all three methods and the prediction of noise correlations and noise power spectra. Since all three aforementioned imaging methods can determine the PED, we focus on PED images for our comparison. Although only projection imaging is considered in this work, a generalization to 3D computed tomography is possible. Using an imaging task that can be viewed as a simplified version of human X-ray radiography, we demonstrate that the combination of spectral and phase-contrast imaging has the potential to generate PED images with strongly reduced noise levels compared with the individual methods. Our theoretical

Forward models and signal extraction
Basis material decomposition for spectral imaging relies on the assumption that the energy-dependent attenuation of any material can be modeled by a linear combination of a few basis materials. Neglecting materials with K-edges in the relevant energy range for medical imaging (≈ 20 − 140 keV), only two basis materials are needed. This is a consequence of the fact that there are only three interaction mechanisms (photoelectric effect, Compton-scattering and Rayleigh-scattering) in this energy range. Furthermore, the contribution of Rayleigh-scattering to the total attenuation cross-section is typically small compared with the other two interaction mechanisms. For simplicity, we will focus on two basis materials in the following. For spectral imaging, the expected number of photon countsŷ s i registered in energy bin s and detector pixel i is thus modeled by [24]: where t(E) is the source spectrum and R s (E) describes the detector response, i.e. the probability that a photon with energy E is detected by energy bin s of the PCD. The functions f 1 (E) and f 2 (E) represent the energy-dependent attenuation of the two basis materials and the corresponding basis material line integrals (for detector pixel i) are denoted by A i 1 and A i 2 . Assuming uncorrelated Poisson statistics, the negative log-likelihood function for spectral imaging is given by: where N and S represent the number of detector pixels and energy bins of the PCD, respectively. The quantity y s i denotes the measured number of photon counts. ML decomposition is performed by finding the basis material line integrals . This optimization problem can be solved separately for each detector pixel. The PED (ρ e ) can be calculated by a linear combination of the basis material line integrals: where ρ V e (M 1 ) and ρ V e (M 2 ) are the volume electron densities of the two basis materials. The idea of approximating the electron density of any material by a linear combination of the two basis materials is conceptually very similar to the standard dual energy assumption of modeling the energy-dependent attenuation with two basis materials. Consequently, equation 3 is also only valid for the dual energy parameter range (low-Z materials, energy range ≈ 20−140 keV). For DPC imaging, the three contrast modalities (attenuation, differential phase shift and dark-field) are extracted from stepping curve measurements that are generated by shifting one of the gratings [10,12]. In contrast to spectral and SDPC imaging, no energy-resolved measurements are acquired. We therefore assume a PCD with just one threshold for DPC imaging. The stepping curve is typically modeled by a cosine (or sine) function [12] and beam hardening effects caused by the polychromatic spectrum are neglected. The expected intensityŷ r i for stepping position r and detector pixel i can thus be written as: [35]: where φ r and V are the reference phase (for step r) and the reference visibility of the stepping curve, respectively. The quantities µ i , ∆φ i and i describe the attenuation of the object, the phase shift of the stepping curve and the visibility reduction (darkfield signal), respectively. The reference intensity measured with an open beam is given by the parameter b. The standard DPC stepping curve model (equation 4) implicitly represents the polychromatic spectrum by an effective X-ray energy. This approximation is reasonably accurate for weakly attenuating samples that only slightly distort the incident X-ray spectrum. Since enhancing soft tissue contrast is one of the main application cases for DPC imaging, this assumption is often justified. The attenuation, differential phase and dark-field images ( µ, ∆ φ, ) are calculated by minimizing the negative log-likelihood of the measured data: where R is the number of stepping positions and y r i denotes the number of photon counts measured for detector pixel i and stepping position r. Similarly to spectral imaging, this optimization problem is separable with respect to the detector pixels. From the differential phase shift, the gradient of the PED can be calculated: where the sensitivity S represents the conversion factor between a PED difference and the corresponding phase shift of the stepping curve. The higher the sensitivity, the larger the phase shift for a given electron density difference between two neighboring pixels. Assuming that the sample is placed between the G1 and G2 grating (compare figure 1), the sensitivity is given by [36]: where r e is the classical electron radius, d is the distance between the G1 and G2 grating, l is the distance between the G1 grating and the object and p 2 is the period of the G2 grating. The parameter a represents the effective pixel size (i.e. the detector pixel size divided by the geometrical magnification of the setup) and E eff is the effective energy of the setup. In analogy to attenuation-based imaging, it can be defined as [37]: where µ(E) and V (E) represent the energy-dependent attenuation of the object and the energy-dependent visibility of the interferometer, respectively. Assuming a known value of the PED at the left and right edges of the detector (for simplicity we assume that the open beam hits the edges of the detector, i.e. ρ 1 e = ρ N e = 0), the PED for an arbitrary pixel can be calculated by integration of the differential phase shifts: Alternatively, it would be possible to integrate starting from the right side of the detector: However, these integration strategies are suboptimal for our goal of providing an indepth noise analysis of all three imaging methods. The summation introduces a strong spatial dependency of the PED variance, even for a homogeneous sample (i.e. a sample that produces the same expected stepping curve measurements for each detector pixel). In this case, the variance of the PED increases linearly from left to right or right to left (if using equation (9) or equation (10), respectively). However, error propagation calculations show (see appendix) that the variance remains constant (even for nonhomogeneous samples) if the electron density is calculated as the average of both summation strategies: In the appendix, we show that the electron density variance is reduced by a factor of two when compared to variance obtained using equation (9) or (10). SDPC imaging can be viewed as a combination of spectral and DPC imaging. The spectral and phase contrast information is used simultaneously by acquiring energyresolved stepping curve measurements. The expected number of photon countsŷ rs i for detector pixel i, stepping position r and energy bin s can be modeled as: where d i is the line integral of an artificial dark-field basis material (see [22] for a more detailed explanation) and f (E) is the corresponding energy-dependency of the dark-field signal. Compared to the forward model of DPC imaging (equation 4), the polychromatic spectrum is taken into account and thus the visibility V (E), reference phase φ r (E) and the phase shift ∆φ i (E) become energy-dependent. In a real experiment, all setupdependent quantities (t(E), R s (E), V (E), φ r (E)) can depend on the spatial position and therefore on the detector pixel index i. For simplicity and clarity, we have omitted this possible dependency in this section. Similarly to DPC imaging, the phase shift depends on the gradient of the PED: where S(E) is the energy-dependent sensitivity of the setup. The key idea behind connecting spectral and DPC imaging and eliminating the PED as an additional optimization variable is expressing the PED as the sum of the projected basis material electron densities: Two basis materials images ( A 1 , A 2 ) and a dark-field ( d ) image can be reconstructed by minimizing the following negative log-likelihood function: As the forward model depends on the basis material thicknesses and their spatial gradients (via the gradient of the PED), the log-likelihood cannot be optimized separately for each detector pixel.

Noise analysis with the Cramér Rao lower bound (CRLB)
The CRLB is a powerful tool from estimation theory that predicts a lower bound for the variance of an unbiased estimator. Given a parameter vector a, it can be shown that [38]: i.e. the matrix C( a) − [F ( a)] −1 is positive semidefinite. Here, C( a) is the covariance matrix of a: where E(·) denotes the expectation value. The Fisher information matrix F ( a) is the expectation value of the curvature of the negative log-likelihood function: From equation 16, a lower bound for the variance of the estimated parameters can be deduced: Since the ML estimator is unbiased and achieves the CRLB in the limit of low noise levels, the CRLB can be used to predict the noise levels for all three imaging methods under consideration. As will be shown later, the CRLB is a good predictor of the noise levels for clinically realistic photon statistics.
In the case of spectral imaging, the optimization problem is separable with respect to different pixels. We therefore obtain N 2 × 2 (inverse) Fisher matrices. For simplicity, the dependency on the pixel index i is suppressed in the following. As derived by Roessel and Herrmann [39], the Fisher matrix is given by: where u, v ∈ (1, 2) are the basis material indices and The elements of the Fisher matrix can be rewritten in a slightly more intuitive form: wheref s u is the weighted average attenuation caused by basis material u in energy bin s. The weights for each energy bin are determined by the effective spectra (including object attenuation). The lower bound for the variances of the basis material line integrals (see equation 19) is then calculated using the analytical inversion formula for 2 × 2 matrices. From the variances of the basis material thicknesses, the variance of the PED can be calculated by standard error propagation: where the covariance of the two basis material thicknesses Cov(A 1 , A 2 ) is estimated by Similarly to the derivation for spectral imaging [39], The Fisher matrix for DPC imaging can be calculated as: Explicitly calculating the gradients leads to: where φ eff r = φ r − ∆φ is the effective phase and Q = be −µ V e − is the expected intensity (be −µ ) multiplied by the effective visibility (V e − ) of the stepping curve. We are particularly interested in the lower bound for the variance of the differential phase shift σ 2 (∆φ), as the PED is calculated by integration of the differential phase shifts (compare equation 11). It is given by: For standard phase stepping with R > 3 equidistantly distributed steps, the off-diagonal elements F 12 , F 21 , F 13 , F 31 vanish [40]. Moreover, numerical evaluations show that F 23 and F 32 are small compared with the corresponding diagonal entries and vanish in the limit of R → ∞. In this case, a simple interpretation of the lower bound for σ 2 (∆φ) is possible: i.e. σ 2 (∆φ) is inversely proportional to the number of phase steps, the average number of photon counts per step and the effective visibility squared. The same result was obtained for a least-squares estimator (instead of an ML estimator) [40]. By applying standard error propagation techniques to equation 11, the variance of the PED can be calculated from the variance of the differential phase shift: For SDPC imaging, the optimization problem cannot be solved separately for each detector pixel. Consequently, there are 3N optimization variables, which we summarize in the vector The Fisher matrix has 3N × 3N entries that are calculated similarly to the other two imaging methods: However, most of the elements of the Fisher matrix are zero because the differential phase shift only couples neighboring pixels. To write the partial derivatives of the forward model more compactly, we define the following abbreviations: The non-zero partial derivatives of the forward model (ŷ rs i ) are given by: The Fisher matrix for SDPC imaging must be inverted numerically to calculate the CRLB.

Prediction of covariances and noise power spectra
In this section, we explain our approach for predicting covariances and noise power spectra for all three imaging methods. Given an estimate of the covariances of the projected electron densities between different pixels, the noise power spectrum (NPS) can be estimated as: where k is the spatial frequency and j is the imaginary unit. In the appendix, we give a quick derivation of this result. For spectral imaging, the material decomposition is conducted separately for each pixel. This means that there are no noise correlations between different pixels: For DPC imaging, the differential phase shifts are also uncorrelated. The covariance matrix of the differential phase shifts C ∆φ is thus a diagonal matrix with diagonal elements C ∆φ ii = σ 2 (∆φ i ). However, the integration step that is necessary to obtain the projected electron densities (compare equation 11) introduces noise correlations. Using error propagation, the covariance matrix C ρe for the projected electron densities can be calculated from the covariance of the differential phase shifts: where B is the transformation matrix between phase shifts (∆ φ) and projected electron densities ( ρ e ): The entries of the matrix B can be deduced form equation 11. In the case of SDPC imaging, we use the entries of the inverse Fisher matrix as an estimate for the covariances: Similarly to DPC imaging, the covariance matrix for the projected electron densities is calculated as: The entries of the transformation matrixB can be deduced from equation 14.

Numerical simulation
To compare the noise level of the projected electron densities for all three imaging methods and to test the predictions of the noise analysis framework, we simulated a radiography measurement of an homogeneous object. The X-ray beam had to penetrate 12 cm of soft tissue and 1 cm of cortical bone, which corresponds to a PED of 12 cm · 3.52 · 10 23 cm −3 + 1 cm · 5.95 · 10 23 cm −3 = 4.82 · 10 24 cm −2 . The aforementioned thicknesses could reflect typical path lengths for medical imaging tasks (e.g. thorax radiography or head CT) . The simulated object had no internal microstructure, i.e. it did not generate a dark-field signal. We assumed a tungsten X-ray source and a PCD with a pixel size of 200 µm, a 2 mm thick cadmium telluride sensor and two thresholds per pixel. As explained in the methods section, only one detector row (containing 400 pixels) was simulated. The detector response was simulated with an empirical model [24] that includes sensor effects (charge sharing, K-escapes). The source spectrum and the detector response were discretized in 1 keV steps. For DPC and SDPC imaging, a symmetric three-grating interferometer (operated at the first fractional Talbot order) was inserted into the beam path. The attenuation gratings (G0 and G2) were assumed to be made of gold with a grating height of 200 µm. The G1 grating was assumed to be made of nickel and to induce a phase shift of π/2 for the design energy. Although the duty cycle of the gratings influences the noise level [41], it was kept fixed at its standard value (0.5) in this simulation study in order to reduce the number of possible acquisition parameter combinations. The total length of the simulated setup was 2.4 m and the object was placed between G1 and G2 (40 cm away from G1). Figure 1 shows an overview of the setup geometry that was kept fixed for all three imaging methods. All other imaging parameters (acceleration voltage, threshold positions and design energy of the interferometer) were optimized individually for each imaging method (if applicable). For all three imaging methods, the low threshold of the PCD was kept fixed at 15 keV. The position of the high threshold was optimized for SDPC and spectral imaging while only one threshold (at 15 keV) was used for DPC imaging. For DPC and SDPC imaging, the design energy was tuned by changing the height of the phase shifting grating as well as the grating periods. Changing the grating periods is necessary to keep the first fractional Talbot distance at the detector position. Stepping curves (using 5 steps equally distributed between 0 and 2π) were simulated with a wave-optical simulation package based on Fresnel propagation [42,43] and the projection approximation [44]. For all three imaging methods, the variable setup parameters were optimized by minimizing the variance of the PED predicted by the CRLB while keeping the dose to the object constant (1 mGy). The dose was estimated based on an empirical model that was fitted to Monte-Carlo simulations [45]. Since the G0 and G1 gratings are placed between the source and the object, they only influence the (spectral) photon flux incident on the object. The G2 grating, however, attenuates photons that have (potentially) contributed to the dose delivered to the object because it is placed between the object and the detector.

Results
The method for finding the optimum setup parameters is illustrated in figure 2 (a) and (b). Figure 2 (a) shows a contour plot of the predicted PED variance for SDPC imaging (at the optimum threshold position) as a function of the acceleration voltage and the design energy. As demonstrated in the last section, the PED for the simulated decomposition task is on the order of 10 24 cm −2 , while the corresponding standard deviations are approximately two orders of magnitude lower. This explains the larger numerical values for the PED variances of up to 10 46 cm −4 . The minimum variance is achieved for an acceleration voltage of 140 kVp and a design energy of 60 keV. As the acceleration voltage is decreased, the optimum design energy also decreases. Figure 2 (b) shows a similar contour plot for DPC imaging. In comparison with SDPC imaging, the optimum acceleration voltage is much lower (for this particular decomposition task). However, similarly to SDPC imaging, the optimum design energy is positively correlated with the acceleration voltage. This behavior is reasonable because a large part of the spectrum should be concentrated around the design energy of the interferometer for optimum performance. Table 1 shows the optimum setup parameters as well as the minimum PED variances for all three imaging methods. Interestingly, the optimum setup parameters for SDPC imaging are very similar to those for spectral imaging, while DPC imaging favors a considerably lower acceleration voltage (60 kVp) and design energy (44 keV). A possible explanation for these results can be found by looking at figure 3, which shows the energy-dependent visibility as well as the effective spectra for DPC and SDPC imaging with optimum setup parameters. The effective spectrum includes the source spectrum, the attenuation of the gratings and the detector response. Kottler et. al [46] demonstrated that a grating interferometer with a π/2 phase shifting grating (G1) exhibits a second visibility peak at twice the design energy. The second important factor that influences the visibility is the attenuation of the gold gratings. In the range from 60 − 80 keV, the gold gratings become increasingly transparent, but for energies above the K-edge of gold (80.7 keV), the attenuation of the gratings and thus also the visibility rises sharply. In the case of DPC imaging, a lower effective energy (see equation 8) is preferable because it increases the sensitivity of the setup (compare equation 6). In other words, a fixed gradient of the PED causes a larger phase shift for lower effective energies. Since the spatial resolution and the setup geometry are fixed for our simulation, the sensitivity can only be tuned via the design energy and the effective energy. The acceleration voltage (60 kVp) is matched to the design energy (44 keV) so that a large fraction of the effective spectrum is concentrated around the visibility peak. Achieving a high visibility is important because for DPC imaging, the variance of the PED is proportional to the visibility squared (compare equation 27). The comparatively high acceleration voltage (140 kVp) for SDPC imaging leads to a decreased sensitivity compared with DPC imaging. However, this decrease in sensitivity is compensated by the improved performance of spectral imaging at higher acceleration voltages. Moreover, the decrease in sensitivity and visibility compared with DPC imaging is less pronounced for the low energy bin. Due to the higher acceleration voltage, the optimum design energy for SDPC imaging is considerably larger (60 keV). Consequently, the visibility  peak for low energies is reduced compared with DPC imaging. On the other hand, the visibility for energies larger than the K-edge of Gold (80.7 keV) is increased. The predicted PED variances show that for this particular imaging task, SDPC imaging can achieve considerably lower noise levels compared with spectral and DPC imaging (variance reduction by a factor of 5.8 and 12.6. respectively). It is important to note that for DPC imaging, the variance depends approximately linearly on the number of pixels in one detector row. This behavior is caused by the integration that is required to convert the differential phase shifts to a PED profile (see equation 11). For SDPC imaging, the PED variance only depends very weakly on the number of pixels. The reason for this weak dependence will be explained later when we analyze the noise correlations for all three imaging methods. We conducted numerical experiments with 100 different dose levels ranging from 0.001 to 10 mGy in order to test if the ML-based estimators achieves the predicted PED variances. For these experiments, we reduced the detector size to 100 pixels to reduce the computational time. For each dose level and imaging method, 10, 000 different noise realizations were simulated and the ML-based estimators (compare equations 2,5,15), were used to calculate the PED. Figure 4 shows the experimentally achieved PED variances (dots) together with the variances predicted by the CRLB (solid lines) for all three imaging methods. All estimators achieve the CRLB for a large range of dose levels (≈ 0.1 − 10 mGy). For lower dose levels, deviations from the CRLB are observable, especially for DPC and SDPC imaging. In the case of DPC imaging, a similar dependency of the variance on the dose level (or photon statistics) has been reported for a Fourier-based estimator [19]. In this simulation, DPC imaging achieves a lower variance than spectral imaging because of the reduced detector size. However, the main goal of this simulation was to test the predictions of the noise analysis framework rather than comparing the performance of the three imaging methods. It has already been pointed out that the integration step for DPC imaging causes long- range correlations between the PED values for one detector row. Consequently, the NPS for DPC imaging is dominated by low frequencies. In the case of spectral imaging, the basis material line integrals can be determined separately for each detector pixel and thus the PED values determined by spectral imaging are uncorrelated (at least for an ideal PCD). Therefore, the noise power should be equally distributed between all frequencies ("white noise"). Since SDPC imaging is a combination of both aforementioned imaging methods, it could be expected that the NPS for SDPC imaging is a mixture between the noise power spectra for the other two imaging methods. In section 2.3, we have presented a framework for predicting the PED covariances as well as the noise power spectra for all three imaging methods. We conducted another numerical experiment to verify the predicted covariances and noise power spectra. In this simulation, the dose was fixed at 1 mGy and the original detector size (400 pixels) was used. For each of the three imaging methods, 50,000 different noise realizations were simulated and the PEDs were calculated separately for each noise realization. Figure 5 (a) shows the experimentally calculated normalized covariances (dots) together with the theoretical predictions (solid lines) for all three imaging methods. More precisely, the covariances between the PED for the central detector pixel (i.e. pixel index 200) and the PEDs of all other pixels is shown. In figure 5 (b), the experimentally obtained noise power spectra are plotted together with the theoretical predictions. We used equation 42 for the experimental calculation of the NPS. For both the covariances and the noise power spectra, the theoretical predictions are in good agreement with the experimentally obtained values. As expected, the covariance of the PED between two different detector pixels is zero for spectral imaging. For DPC imaging, the covariance decreases linearly with the distance from the central detector pixel until it drops to zero at the edges of the detector (not shown in figure 5 (a)). This reflects the well-known long-range noise correlations introduced by the integration step. In the case of SDPC imaging, however, the covariance rapidly drops to zero with increasing distance from the central pixel. For the imaging parameters given in table 1, almost no noise correlations are observable if the distance to the central pixel is larger than approximately 25 pixels. Compared to DPC imaging, the additional spectral information eliminates long-range noise correlations. On the other hand, it could be argued that compared with spectral imaging, the additional phase shift term for SDPC imaging couples the PED values of neighboring pixels and thus introduces local noise correlations. The covariance graph for SDPC imaging explains the aforementioned weak dependence of the PED variance on the number of detector pixels. Contrary to DPC imaging, the covariance rapidly decreases with the distance between two pixels and thus the assumption of fixed, known PED values at the edges of the detector does not influence the noise level (except for pixels close to the edges). The different noise correlations for the three imaging methods are also reflected in the different noise power spectra. As expected, the noise power is independent of the frequency for spectral imaging, whereas the noise power for DPC imaging increases drastically for lower spatial frequencies. For high frequencies, the noise power spectrum for SDPC imaging is similar to DPC imaging. For low frequencies, however, the noise power does not increase but converges to a constant value that lies slightly above the noise power graph for spectral imaging. The higher noise power for SDPC imaging in the low frequency area can be explained by the attenuation of the G2-grating, which removes a part of the X-ray beam that has contributed to the dose delivered to the object. In our simulation study, we assumed a fixed total length of the setup and a symmetric grating interferometer that operates in the first fractional Talbot distance. Given these conditions, the effective pixel size is the most important tuning factor for the sensitivity of a DPC or SDPC imaging setup (compare equation 7 and 13). Since the PED variance strongly depends on the sensitivity, we investigated the theoretically predicted performance of all three imaging methods as a function of the effective pixel size. For each pixel size, the optimum imaging parameters were determined with the noise analysis framework presented in the methods section. Figure 6 (a) shows the predicted PED variances as a function of the effective pixel size. For visualization purposes, the number of photon counts per pixel (instead of the dose delivered to the object) was kept constant. This means that the dose delivered to the object was chosen inversely proportional to the squared effective pixel size (i.e. the effective pixel area). If the dose was kept constant, the PED variances would have to be multiplied by an additional factor that is inversely proportional to the effective pixel area for all three imaging methods. Under these conditions, the PED variance for spectral imaging is  (see table 1) as a function of the visibility. For this experiment, the energy-dependent visibilities (see figure 3) were multiplied by a visibility reduction factor m (0 < m < 1) to simulate a decreased visibility due to experimental effects.
independent of the pixel size. In the case of DPC imaging, the variance is proportional to the squared effective pixel size, which reflects the well-known sensitivity-dependent noise level for DPC imaging (compare also equation 28). For large effective pixel sizes (a > 4 mm), the PED variance for SDPC imaging is almost constant and approximately twice as large as the variance for spectral imaging. It coincides with the PED variance one would obtain if spectral imaging was performed with the "unnecessary" grating interferometer inserted into the beam path (see dashed line in figure 6 (a)). It follows that for larger pixel sizes (or low sensitivities), the additional information on the basis material line integrals provided by the phase shift of the stepping curve can be neglected compared to the spectral information. Consequently, SDPC imaging reduces to spectral imaging in the limit of large effective pixel sizes. For this particular imaging task, SDPC imaging and spectral imaging have equal PED variances for an effective pixel size of a = 950 µm. At this point, the benefit of the additional phase shift information is counterbalanced by the reduced photon statistics caused by the attenuation of the G2 grating. For very small pixel sizes (a < 10 µm), the SDPC variances almost coincide with the DPC variances. In this case, SDPC imaging reduces to DPC imaging because the spectral information provides little additional value for determining the PED in comparison to the phase shift of the stepping curve. In the area between these limiting cases (10 µm < a < 950 µm), the SDPC variance lies below both the spectral and the DPC variance because both the spectral and the phase shift information are relevant for determining the PED. The largest variance reduction (by a factor of 8.8) compared to both DPC and spectral imaging is achieved for the point of equal variance of spectral and DPC imaging (a = 90 µm). At this point, the spectral and the phase shift information are of equal importance for determining the PED.
It is important to note that the predicted PED variance only depends on the geometrical imaging parameters (pixel size, object position, grating distances, grating periods) via the sensitivity. Consequently, similar trends for the PED variances could have been observed by varying geometrical parameters other than the pixel size. However, increasing the sensitivity via these geometrical parameters also increases the phase shift of the stepping curve that a given sample causes (via the gradient of the PED). For DPC imaging, this might lead to additional phase wrapping artifacts. As will be discussed in the next section, phase wrapping artifacts can in theory be avoided for SDPC imaging. The number of detector pixels (N = 400) was kept constant in this analysis of the noise characteristics. This means that the field of view decreases with decreasing effective pixel size a. If the field of view was kept constant, the total number of detector pixels would have to be scaled by 1/a. As discussed in the results section, increasing the number of detector pixels would have an adverse effect on the PED variance for DPC imaging, but the influence on SDPC imaging would be negligible as long as the spectral information prevents long-range noise correlations. In a real experiment, various undesirable effects (e.g. vibrations, grating imperfections) could cause a reduced visibility of the stepping curves compared with the simulations. Moreover, we have assumed that the object does not generate a dark-field signal which would also reduce the effective visibility. We therefore investigated the influence of a visibility reduction on the PED variance for SDPC and DPC imaging. For this study, the imaging parameters were kept fixed at their optimum values (compare  table 1). We simulated a visibility reduction by multiplying the energy-dependent visibilities (compare figure 3) with a visibility reduction parameter m (0 < m < 1). Figure 6 (b) shows the predicted PED variances for SDPC and DPC imaging as a function of the visibility reduction parameter. For comparison, the variances for spectral imaging with and without the interferometer in the beam path are also plotted as constant lines. As expected, the variance for DPC imaging rises rapidly with decreasing visibility (σ 2 (ρ e ) ∝ (V e − ) −2 , compare equation 27). In the case of SDPC imaging, however, the variances rises much more slowly with decreasing visibility (approximately proportionally to (V e − ) −1 ) before it approaches a constant value for low visibilities.
This is reasonable because in the limit of low visibilities, SDPC imaging effectively reduces to spectral imaging (with the gratings in the beam path). For m = 0.15, the variances for SDPC imaging and spectral imaging are equal. This means that for this particular reconstruction task, the visibility can be reduced by a factor of 1/m = 6.7 compared with the idealized simulation before the attenuation of the G2 grating outweighs the variance reduction achieved by the additional phase shift information.

Discussion
A quantitative comparison between conventional attenuation imaging and DPC imaging is difficult because of the different underlying contrast generating mechanisms for these two imaging methods. For example, a contrast-to-noise ratio comparison [16,17] is complicated because the differences in contrast for the two imaging methods depends on the choice of the investigated materials. Moreover, depending on the setup parameters, the contrast between different materials can vanish completely for both imaging methods.
For this reason, we concentrated on electron density images, which can be calculated with all three imaging methods (spectral, DPC and SDPC imaging) under consideration in this comparison study. Nevertheless, the presented noise analysis framework can also be used to determine the optimum acquisition parameters for other imaging tasks, such as dark-field imaging or basis material decomposition. Although we have focused on projection imaging, the noise analysis framework could be generalized to computed tomography imaging (e.g. by using error propagation for the filtered backprojection [47]). For the investigated imaging task, our analysis predicts highly reduced noise levels for SDPC imaging in comparison to both DPC imaging and spectral imaging. To achieve a fair comparison, the imaging parameters were optimized individually for each of the three methods. For a large range of clinically realistic dose levels, the ML-based estimators achieve the predicted noise levels for all three imaging methods. The validity of the predicted covariances and noise power spectra was confirmed by additional numerical simulations. Our analysis focused on the projected electron density, but earlier studies have shown that the noise advantage of SDPC imaging compared to spectral imaging is also transferred to the basis material images [22]. In our study, DPC imaging only outperforms spectral imaging for comparatively high spatial resolutions (a < 90 µm), however this conclusion does not apply to SDPC imaging. For a large range of clinically relevant effective pixel sizes (up to 950 µm), SDPC imaging theoretically outperforms the other two imaging methods by simultaneously using both the spectral and phase contrast information. Although we have only considered one imaging scenario in this work, the general trends should apply to a large range of imaging tasks. Nevertheless, the location of the break-even points between the three imaging methods will vary depending on the object size and the chemical composition. In a real experiment, both grating imperfections and ultra-small angle scattering by the object (i.e. nonzero dark-field sig-nals) could reduce the visibility of the stepping curves compared to our simulations. The noise reduction for SDPC imaging compared to spectral imaging will thus be smaller than theoretically predicted. Nonetheless, because of the large theoretical noise reduction factors, we believe that SDPC imaging can also outperform the other two imaging methods in real experiments. Moreover, the noise level for SDPC imaging increases less rapidly with decreasing visibility compared with DPC imaging (see figure 6 (b)). Regardless of the potential for noise reduction, the combination of spectral and phasecontrast imaging provides additional information that is inaccessible with the individual imaging methods. Compared with spectral imaging, the dark-field image yields additional information about the object's microstructure and, compared with DPC imaging, two basis material images can be calculated and beam hardening artifacts in all three imaging channels can be corrected. Raupach and Flohr [20] state that the noise correlations for DPC imaging are disadvantageous for clinical diagnosis. They argue that the dominance of low frequencies in the NPS is unfavourable for recognizing structures. Furthermore, because of the lowfrequency noise, reducing the noise level by pixel binning is less effective compared with conventional attenuation-based imaging. In general, these arguments do not apply to SDPC imaging because of the fundamentally different noise correlations compared to DPC imaging. An exception is the limit of very small pixel sizes (or high sensitivities), where SDPC imaging effectively converges to DPC imaging. In our study, the influence of the spectral information causes a rapid decrease of the covariance between two pixels with increasing distance. Consequently, there is no excessive low frequency noise. Furthermore, Raupach and Flohr [19] argue that -contrary to conventional attenuation imaging -there is a dose limit for DPC imaging below which no meaningful signals can be extracted. The interferometric measurement acquisition process only allows unique determination of the phase shift of the stepping curve across a 2π interval. Larger phase shifts are wrapped back into this interval. According to Raupach and Flohr, statistical phase wrapping leads to a loss of the phase shift information for low dose levels. More precisely, the average extracted phase shift is biased to zero, regardless of the underlying physical phase shift if a standard averaging procedure is used. As a result of the additional spectral information when compared to DPC imaging, the calculated phase shift is not restricted to a 2π interval for SDPC imaging. Qualitatively speaking, the spectral information determines in which 2π interval the phase shift lies (via the PED profile) and the exact value is fine-tuned by the stepping curve information. In principle, there is thus no information loss through statistical phase wrapping. However, it has been shown that the log-likelihood function for SDPC imaging has local optima that can be explained by an analogy to the phase wrapping problem for DPC imaging [22]. Although we could avoid the convergence to local optima in previous simulation studies by choosing a suitable initial guess, this strategy is likely to fail at extremely low dose levels. Nevertheless, the convergence to local optima is an optimization problem rather than a fundamental restriction of SDPC imaging because it could in principle be avoided with a global optimization strategy (or by incorporating prior information in the form of a regularization term).

Conclusion
In this work, we have developed a noise analysis framework that allows the calculation of (co-) variances and noise power spectra for spectral, DPC and SDPC imaging. An important practical application of this framework is finding the optimum imaging parameters for all three methods. Our analysis shows that the combination of spectral and phase-contrast imaging is a promising imaging technique with various advantages compared with the individual methods. SDPC imaging provides additional information compared with both spectral imaging (dark-field image) and DPC imaging (basis material images). Moreover, we demonstrated that SDPC imaging enables a strong noise (or dose) reduction compared with the other two imaging methods for a large range of clinically relevant pixel sizes. Finally, the additional spectral information compared to DPC imaging eliminates excessive low-frequency noise, which can be a major drawback of DPC imaging, especially in projection space.

Noise propagation for the integration step in DPC imaging
Applying standard error propagation techniques to equation 9, the variance of the PED for a certain pixel (σ 2 (ρ i e )) can be calculated: where σ 2 (∆φ q ) is the variance of the differential phase shift (which can be calculated using the CRLB). In the special case of a homogeneous sample (i.e. σ 2 (∆φ q ) = const. = σ 2 (∆φ) ∀q), the average variance of the electron density σ 2 avg (ρ e ) can be calculated as: The pixels with index 1 and N are excluded from the average variance calculation because it was assumed that the PED is known with certainty at the edges of the detector. The same result is obtained if the integration from right to left (equation 10) is considered. If the left-and right-integrated PED profiles are averaged (see equation 11), the variance of the PED is given by: As can be seen from equation 40, the variance of the electron density is independent of the detector pixel index. Consequently, the variance is spatially constant, even for a nonhomogeneous sample. In the special case of a homogeneous sample, σ 2 avg (ρ e ) is given by: Comparing equation (40) and (39), it follows that (in the case of a homogeneous sample) the average variance is reduced by a factor of two when averaging left-and rightintegrated electron density profiles.

Calculating the noise power spectrum from the covariances
In this section, we derive how the noise power spectrum (NPS) can be calculated from an estimate of the covariances (compare equation 32). For a homogeneous object, the NPS can be calculated as [48]: where k is the spatial frequency andρ e (x) is the offset-corrected PED as a function of the spatial coordinate x. In discrete form, the NPS is given by: Since E ρ