A method for material decomposition and quantification with grating based phase CT

Material decomposition (MD) is an important application of computer tomography (CT). For phase contrast imaging, conventional MD methods are categorized into two types with respect to different operation sequences, i.e., “before” or “after” image reconstruction. Both categories come down to two-step methods, which have the problem of noise amplification. In this study, we incorporate both phase and absorption (PA) information into MD process, and correspondingly develop a simultaneous algebraic reconstruction technique (SART). The proposed method is referred to as phase & absorption material decomposition-SART (PAMD-SART). By iteratively solving an optimization problem, material composition and substance quantification are reconstructed directly from absorption and differential phase projections. Comparing with two-step MD, the proposed one-step method is superior in noise suppression and accurate decomposition. Numerical simulations and synchrotron radiation based experiments show that PAMD-SART outperforms the classical MD method (image-based and dual-energy CT iterative method), especially for the quantitative accuracy of material equivalent atomic number.


Introduction
Conventional polychromatic X-ray CT aims at reconstructing the linear absorption coefficient, which depends on material composition and energy-related attenuation performance. This imaging modality relies on its X-ray spectrum, which leads to a weak ability to identify substances. For example, in medical diagnosis, traditional CT often fails to distinguish iodine contrast agents and bone tissue [1]. To overcome the shortcoming, dual-energy CT is developed, which collects projections with two different X-ray spectra and is able to perform selective reconstruction for equivalent atomic number and atomic density, basic-material based concentration distributions, and so on. The commercialized dual-energy CT, like by SIEMENS and GE, can distinguish iodine contrast agents and bone tissue, tendons and ligaments, blood vessels and bone tissue, and so on [2]. The decomposition rationale for dual-energy CT depends on the fact that the energy attenuation behavior of matter varies with different X-ray spectra. However, this assumption is challenged in the real application with low-Z compounds, which have weak absorption and poor discrimination. Different from traditional absorption-based X-ray CT, phase-based X-ray CT imaging employs both the absorption feature and the phase shift manifestation to reveal the anatomic characteristics. For low-Z material, although the absorption behavior is very similar, the performance on phase shift is quite different [3]. Therefore, theoretically, phase CT has a superior material distinguishability than absorption CT. Methodologically, the phase-contrast CT can be implemented by a grating interferometer based differential phase imaging method, which obtains three perspective images in terms of absorption, differential phase (refraction angle) and scattering. Technically, by employing CT reconstruction technique, the spatial distribution of linear absorption coefficient (μ), the real part of refractive index (δ) and the scattering coefficient can be obtained simultaneously. Assuming the X-ray is monochromatic, these physical indexes can be calculated with traditional analytical algorithms, such as filtered backprojection (FBP). Specifically, the μ and scattering coefficient can be reconstructed by the Shepp-Logan filter based FBP algorithm, while the δ can be reconstructed by the Hilbert filter based FBP [4] or BPF algorithm [5]. Moreover, they can be reconstructed from algebraic iterative reconstruction algorithm [6,7].
Wang et al. employed absorption and small-angle scattering information to quantitatively measure the material composition and concentration [8], which works for the measurement of volumetric breast density (VBD) and the classification of breast micro-calcification. However, the small-angle scattering information originates from the unevenness inside a resolution unit, which is inappropriate for the substance with uniform composition inside the resolution unit but gradual variation between neighbor resolution units. To overcome this limitation and to improve the reliability of quantitative analysis, considering the high sensitivity of the phase signal, both phase and absorption information are combined for matter decomposition. In 2010, Qi et al. reconstructed μ and δ, and meanwhile calculated the three-dimensional spatial distribution of the material electron density and the effective atomic number [9]. In 2017, Han et al. employed a grating interferometer to acquire absorption and differential phase information, established a two-component absorption and differential phase equation, and performed twosubstrate based decomposition in projection domain [10,11]. In 2018, Braig et al. proposed a two-substrate equation for material decomposition, and investigated material electron density and effective atomic number [12].
The aforementioned methods belong to a two-step strategy, which separately perform decomposition and reconstruction. One strategy is decomposing in projection domain, such as Han's work [10,11]. The other is in image domain, such as Qi [9] and Braig [12]. However, whatever the sequence is, noise amplification exists all the time. For projection-domain decomposition, both phase and absorption projections need to be known at the same time. When obtaining the phase projection by integrating the differential phase shift data along the X-ray refraction direction, strip artifacts are introduced and low-frequency noise are amplified in the phase projection [13,14]. For image-domain decomposition, it is necessary to know the reconstructed μ and δ. When reconstructing the δ, the Hilbert filter based method lacks a noise-suppression mechanism and suffers serious noise influence. All these two approaches are presented in the Methods section.
In order to overcome the inevitable noise amplification of the two-step methods, we propose an optimization-based one-step decomposition method (PAMD-SART, phase & absorption material decomposition-SART), which directly reconstructs substrate images from the differential phase and absorption projections. Noticeably, the optimization model contains a noise-suppression mechanism and can be effectively solved by an iterative scheme.

Imaging theory
The interaction between X-rays with matter can be described by the following complex refractive index where δ is the real decay rate of the refractive index, b ¼ l 4p m, μ the linear absorption coefficient, and λ the wavelength of the incident X-ray. When an X-ray plane wave with amplitude A 0 passes through the object, the wave function of the emergent beam reads where F is the matter-caused X-ray phase shift with the expression as Let M describe the absorption of X-rays, which can be formularized as After penetrating the object, the intensity of X-ray decays to where I 0 ¼ A 0 � A � 0 is the intensity of the incident X-ray, A � 0 is the conjugation of A 0 . X-ray differential phase imaging is based on the usage of grating interferometer. By employing a phase stepping acquisition method, absorption and differential phase data can be extracted from each scan angle [15]. The absorption projection data can be written as and the differential phase projection data (refraction angle) as Here (x, y) is the sample coordinate system, (x φ , y φ ) the measuring coordinate system, and φ the angle at which the measuring coordinate system rotates with respect to the sample coordinate system. In (7), @δ(x, y)/@x φ changes with the projection angle, which leads to a failure of the traditional CT algorithm. However, Hilbert filter based FBP algorithm [4] or the BPF algorithm [5] can work in this case to recover δ(x, y). According to the material decomposition method, μ of the measured object can be expressed as a linear combination of μ of two basic materials. Under the condition of monochromatic X-ray exposure, δ can also be linearly expressed by δ of two basic materials [11]

PLOS ONE
where μ 1 and μ 2 are the linear absorption coefficients of the two basic materials, δ 1 and δ 2 the real decay rates of the two basic materials, and f(x, y) and g(x, y) the distribution function of the two substrates. By substituting equations (6) and (7) into Eq (8) and switching the differential and integral operations, we obtain the following relationships, To demonstrate the superiority of the phase and absorption based substrate decomposition than the conventional dual-energy based one, we illustrate a simple example in Fig 1. As shown in Fig 1A, assuming two X-ray beams with energy of 20 keV and 30 keV successively pass through water and PMMA with the same thickness of 1 mm, the corresponding absorption equations read, where Δ f and Δ g are the thickness of water and PMMA, which need to be solved, and the unit of thickness is cm, so Δ f = 0.1 and Δ g = 0.1 in (10). The constant value in (10), 0.737(0.376) and 0.623(0.362) are the μ value of water and PMMA at 20(30) keV, the unit is cm −1 , while 0.136 (0.074) is −ln(I/I 0 ) at 20(30) keV. For the 20 keV X-ray beam, according to (8), we incorporate the refractive indices and establish the following equation system, where the constant values 0.526 and 0.630 in the second line of (11) are δ values of water and PMMA at 20keV, 0.116 is the calculated δ value in this configuration. As shown in Fig 1B, the three equations in (10) and (11) are graphically expressed by straight lines. The angle between the two straight lines is 3.70˚for (10), and 9.93˚for (11). In addition, the condition numbers of the above two equations are 36.96 and 11.66, respectively. All the evidences show that (11) is more stable and robust than (10), i.e., superior tolerance to noise and data error.

Reconstruction algorithm
In this section, we introduce the discrete mathematical model of the phase-based basic material decomposition, investigate the derivation process, and present the implementation steps of the proposed PAMD-SART method.
Imaging model. Let f = (f 1 , f 2 , . . .f J ) τ and g = (g 1 , g 2 , . . .g J ) τ denote the discretized images of f(x, y) and g(x, y), where f j and g j are the sampled values of f(x, y) and g(x, y) at the jth pixel, J the total pixel number, and τ the vector transpose operation.
:::y φ U Þ t denote the absorption projection and differential phase projection extracted at projection angle φ, where M u and θ u are the extracted values by the uth detector cell, and U the number of detector cells. A φ ¼ ða φ uj Þ U�J is the projection matrix at angle φ, and a φ uj the contribution of f j and g j to the uth ray path at projection angle φ. Then (9) can be discretely expressed as where D is the discrete form (matrix form) of the differential operator −@/@u. Noticeably, (12) tells us that the distribution images of the two substrates can be reconstructed from differential phase and absorption projections. Solution algorithm. The projection-domain based material decomposition is a traditional method [11]. This method uses the inverse of the differential operator D −1 and the boundary condition to solve D −1 θ φ . Applying D −1 directly spreads the signal error, and leads to stripe artifacts. In order to effectively suppress these artifacts, optimization based methods are taken into consideration [13,14]. First, (12) need to be converted as, Then the traditional reconstruction method is employed to recover f and g. The other kind of methods can be viewed as image-domain based material decomposition, which converts (12) to and directly reconstruct μ 1 f + μ 2 g and δ 1 f + δ 2 g. μ 1 f + μ 2 g is calculated from M φ with traditional FBP algorithm, and δ 1 f + δ 2 g is recovered from θ φ with Hilbert filter FBP algorithm. Finally, f and g are easily obtained by solving an linear equation system with two variables. Apparently, the first method contains an integral operation in the projection domain, and the second method requires a Hilbert filter operation in the image domain, both of which inevitably amplify the noise. If we use an optimization based reconstruction algorithm to recover μ [6] or δ [7], the noise can be suppressed. However, the edge of low contrast object might be blurred, since it is difficult to adjust the regularization parameters to balance the noise and the edge in a image with high dynamic range.
PAMD-SART algorithm. The algorithm proposed in this paper is to obtain the final decomposition image by solving the following optimization model: The definition of each symbol is consistent with section Imaging model and R(�) is a regularization operator such as total variation (TV) [16,17], gradient L0 [18], Mumford-Shah [19] or an image filter. In the following, inspired by the SART reconstruction strategy [20], an iterative algorithm to directly reconstruct f and g from the absorption projection M φ and the differential phase projection θ φ is given.
Assuming the current estimate images are (f (m) , g (m) ), then the projection estimations of absorption and differential phase can be calculated by Furthermore, the residual can be expressed as e ðmÞ M ¼ M φ À M φðmÞ and e ðmÞ D À 1 y ¼ D À 1 ðy φ À y φðmÞ Þ. By simplifying (15) to each angle and setting the absorption and phase term to zero, one can get the follow equations, ( After solving the above equations, we obtain Finally, by using SART algorithm, we get the iterative form where e ðmÞ M;u and e ðmÞ D À 1 y;u are the uth components of e ðmÞ M and e ðmÞ D À 1 y , λ the relaxation factor, A φ u;þ ¼ P J j¼1 a φ u;j and A φ þ;j ¼ P U u¼1 a φ u;j with u = 1, 2, . . .U and j = 1, 2, . . .J. In the iterative scheme, the calculation of phase residual e ðmÞ D À 1 y requires to inverse the differential operator D −1 , which may diffuse the error. In order to overcome the shortcoming, we employ the following optimization model and use a gradient descent method with a fixed step number in real applications. The elaborated iteration steps of the PAMD-SART method is shown in Algorithm 1.

Verification
In this section, we performed numerical simulations and real experiments to verify the proposed method, including visual comparison of material decomposition results and quantitative analysis of acquired equivalent atomic number.
The image quality was measured with peak signal-tonoise ratio (PSNR) where I is a reconstructed image, I � the ground truth image, N the number of image pixels in I, and Gr the maximum gray value of I � .

Numerical simulations
The numerical phantom was derived from the FORBILD head phantom [21]. As shown in Fig  2A, the size is 9.8mm � 7.8mm, including bone and water. In order to enhance the comparison performance, we added multiple water-like objects to the original image. The specific density is shown in Fig 2C. In the simulated experiment, for the sake of simplicity, we used 20 keV monochromatic X-ray. The μ and δ of water and bone were from the X-ray optics software toolkit(XOP) [22]. First, the projection data of the sample was obtained using the forward projection method, and then Poisson noise was added, where the number of original photons was 10 6 per ray. The differential phase projection was obtained by performing center-difference to the phase projection, and the δ values of water and bone were multiplied by the same ratio so that they were similar to the absorption value μ. In the simulation, a parallel-beam setting was In image reconstruction, bone and water were selected as the basic materials, and the reconstructed image size was 512 × 512. The PAMD-SART method was implemented by C++ and CUDA with GPU acceleration. For better comparison, we implemented regularized iterative method to reconstruct μ [6] and δ [7] in the first step of the image based method. TV was employed as a regularization term and the solution algorithm was from the Chambolle's method [17]. This method is very efficient and easy to be implemented, of which only one regularization parameter (λ TV ) is used to adjust the image fidelity and TV. According to our experimental results, λ TV = 0.1 � mean(image) was a good beginning for Chambolle's method. Specifically, ray-casting method was used for forward projection process, and pixel-driven method for backward projection.
The decomposed results under noise-free and noisy cases are shown in Fig 3. For each case, both image-based method and PAMD-SART were performed. We zoomed in the local region of water-like objects and illustrated in the third column to improve the visual comparison. The profile of water-based image in noisy case and the convergence curve of the PAMD-SART method are shown in Fig 4. The PSNR values for both methods in the noise-free and noise cases are shown in Table 1.
From the results in Figs 3 and 4 and Table 1, we can see that the PAMD-SART method is superior to the image-based method. In the noise-free case, looking at the first and second lines of Fig 3, both methods can reconstruct the result without artifacts, but when comparing the circle structure with the lowest contrast in the third column, the boundary of the structure in image-based method is a little blurry, while the result of PAMD-SART is very sharp. In the noisy case, PAMD-SART performs much better than image-based method. From Table 1, the PSNR value of the image-based method is smaller than PAMD-SART, which indicates that the former keeps a larger error. From the circular water-like material region in the third column of Fig 3, the two materials with the lowest density can hardly be distinguished from the imagebased method, but are clear in PAMD-SART. From the profile in Fig 4, we can see that there is a large oscillation in the image-based result, and the material with smallest density is too close to the water to be separated. However, PAMD-SART can well identify the small contrast change. The curve of relative error indicates the convergence of PMAD-SART in numerical.
The reason for this result is that the original image has a large dynamic range, it is difficult to balance the fidelity term and the regularization term when directly regularizing the original image. When suppressing noise artifacts, the border of low contrast objects will inevitably be blurred. The proposed PAMD-SART method uses projection information of μ and δ simultaneously in the reconstruction process, which reduces the noise interference. Moreover, PAMD-SART performs the regularization on substrate images which has smaller dynamic range than the original images, so has less effect on low-contrast objects. Because PAMD-SART has the advantages in these two aspects, it is much superior to the image-based method.

Real experiment
Real experiments were performed using the X-ray grating interferometer on the BL13W experimental station by Shanghai Synchrotron Radiation Facility. The grating interferometer consisted of a π/2 phase grating with a period of 2.396 um and an absorption grating with a period of 2.4 um. The distance between the two gratings was 46.38 mm and the photon energy used in the experiment was 20 keV. The sCMOS X-ray detector with effective pixel array of 2048 × 600 was applied, and the pixel size was 6.5 × 6.5 um 2 . In order to obtain the phase and absorption projection, the step method was used for data acquisition. The number of steps was 8, which meant the phase grating moved 8 steps in a period to take samples separately, then completed the data acquisition from one angle. Finally, the data of 540 angles were collected at equal intervals within 180 degrees. The average digital value of 8 referenced acquired image in a period was about 25000, and the fringe visibility about 50%. The absorption and differential phase projections were calculated from the data collected at each angle by the phase step formula in [23]. The reconstructed image size was 1024 × 1024, and the maximum iterations was set as 200 which was the stop condition of PAMD-SART.
Phantom experiment. As shown in Fig 5A,   In this experiment, we used water and PTFE as substrates to perform PAMD-SART and image-based decomposition. Then based on the substrate images, μ and δ images and equivalent atomic number were calculated. This result was compared with the dual energy CT method. Fig 5 shows the decomposition results of different methods. Fig 5B is the absorption coefficient tomography. The outer ring in the image represents the polyethylene plastic container, and there are three different components in the inner ring, namely, three circles with different gray levels in the middle of the ring, LDPE phantom with the lowest gray level on the left, PTFE phantom with the highest gray level at the bottom, PMMA phantom at the upper right and water at the rest. In Fig 5C, it can be seen that LDPE and water are difficult to be distinguished. This phenomenon shows that the contrast of phase information is not necessarily higher than the absorption. Comparing the decomposition results of image-based and PAMD-SART methods, i.e., Fig 5D-5I, it is obvious that PAMD-SART method has better noise removal performance.

PLOS ONE
Both works of Qi [9] and Braig [12], refer to the calculation of electron density ρ e and equivalent atomic number Z from absorption μ and phase δ, by employing the following relationship where C p , C E , C Z are parameters to be determined, C PC ðEÞ ¼ r 0 h 2 c 2 2pE 2 , r 0 the classical radius of the electron, h the Planck constant, c the speed of light, and the Klein-Nishina cross section as follow, with a = E/511 keV the relative mass energy to electron. After obtaining ρ e and Z, the virtual absorption coefficient tomogram at other energies can be further calculated.

PLOS ONE
In the proposed method, by using the obtained substrate images, the refraction and absorption images can be calculated by (8), and then introduced to (22) to get ρ e and Z. The theoretical equivalent atomic number Z for a compound was calculated by the following equation [24] Z ¼ ð where ω j is the fraction of the total number of electrons associated with each element, and Z j is the atomic number of each element. According to the water absorption coefficient relationship with energy, the calculated equivalent energy is 20.22 Kev. The absorption and phase coefficient of the four materials (water, PMMA, PTFE and LDPE) are used in (22) to fit the coefficient in the μ formula. After fitting, the coefficients are C p = 1.004, C E = 3.027, and C Z = 3.619. Using the μ in formula (22) with two different energies, the atomic number can also be calculated, which is employed in dual-energy CT [25]. Hence, we acquired two CT datasets by using 20keV and 12keV on the synchrotron radiation station, and the number of views was the same as the grating based scanning. The integration time of detector was adjusted, so that the digital value of acquired image under direct exposure was around 50000. Compared to the data acquired with grating, from the view of detector, the total dose of two monochromatic CT was about half of the grating.
As shown in Table 2, the fitted results are very close to the theoretical values. Fig 6 and Table 2 show that the accuracy of the atomic tomographic image reconstructed by the substrate image is close to the theoretical value, while the result by the dual-energy method contains a large error for PTFE and LDPE. Theoretically, the equivalent atomic number of PTFE should be greater than water. However, when using the dual-energy method with 12 keV and 20 keV X-ray beams, the calculated atomic number of PTFE is smaller than water. Once ρ e (x, y) and Z(x, y) were obtained, μ(x, y) at any energy can be calculated by formula 22, which generated a virtual image. In the 12 keV virtual image, it can be found that the μ ratio of PTFE to water is 2.78, while is 2.41 in the 12 keV actual absorption image. These results demonstrate that the empirical formula of atomic number calculated by dual energy (Z C (μ 1 , μ 2 )) has a large error in the low energy case. Furthermore, by comparing and analyzing the absorption coefficients of various substances, it can be found that the empirical formula (Z C (μ 1 , μ 2 )) is no longer valid when the X-ray energy is less than 20 keV. This is mainly because the linear attenuation coefficient in (22) ignores the coherent scattering term (Rayleigh cross section). However, Table 3 shows that in the energy below 30 keV this term is non-negligible. In Table 3, the cross section of photon interaction are from xraylib [26], these data can also be found in Evaluated Nuclear Data Library(ENDL) [27].
Biological sample experiment. Taking into account the blooming applications in biological sample imaging, we employed a chicken paw bought from supermarket as the specimen.

PLOS ONE
complex, in order to keep the image fidelity, the regularization term at all reconstruction was small. In this experiment, the refraction image contained more artifacts than the absorption image. The reason lied in the dry chicken paw included a lot of air regions, which caused a sharp change in the refractive index and further led to these artifacts. It was noticeable that the noise of the water-based and bone-based tomographic images was not amplified by the PAMD-SART method, but almost the same as the image-based method. In the equivalent atomic number image, the Z C of water in PAMD-SART method was more close to theoretical value than the image-based method. As is shown at the dotted frame part of bone in Fig 7D  and 7G, the image-based method suffers an error at some part. It is because δ is less than zero, but μ larger than zero. However, the PAMD-SART method avoids this error, since it adopts the projection information of both μ and δ in the reconstruction process. Thus, it indicates that the proposed PAMD-SART method is also better than image-based method at complex biological applications.

Discussion and conclusion
In this paper, we propose an iterative method, PAMD-SART, for phase and absorption based substrate decomposition. This method directly reconstructs basic images with desirable noise suppression and boundary maintenance. Numerical simulations and real experiments validate the superiority of the proposed method than the traditional image-based method. In addition, the obtained absorption and phase images are further quantitatively analyzed in terms of equivalent atomic number. The corresponding results are close to the theoretical values, and are more precise than the dual-energy method for low X-ray energies. The proposed PAMD-SART method works well for the synchrotron monochromatic projection data, and will be further extended to the laboratory light source CT system according to the equivalent monochromatic method.
The phase contract CT has not yet been applied to clinical practical applications, and is mainly limited by the manufacture of gratings [28][29][30]. The method proposed in this paper is a primary research. We verify the effectiveness of PAMD-SART with the data acquired by

PLOS ONE
grating interferometer at synchrotron radiation station. And the size of phantom is just about 10mm due to the limitation of low energy and the limited size of grating and detector.
When calculating the equivalent atomic number (Z), there is actually no obvious comparability between phase CT and dual energy CT. It is because the principles of these two imaging patterns are not the same. However, the comparison experiments show that dual energy has more clear conditions for Z calculation. In the case of low energy (less than 30 keV), using two absorption coefficient equations which ignore the coherent scattering (Rayleigh scattering) term to solve Z will have a large error. While in conventional applications of medical or industrial CT, the photon energy is relatively high (more than 40 keV), so it is feasible to ignore the coherent scattering term. Validation and comparison in polychromatic and high energy X-ray are beyond the scope of this study, and will be investigated in future.