Spline based iterative phase retrieval algorithm for X-ray differential phase contrast radiography

: Differential phase contrast imaging using grating interferometer is a promising alternative to conventional X-ray radiographic methods. It provides the absorption, differential phase and scattering information of the underlying sample simultaneously. Phase retrieval from the differential phase signal is an essential problem for quantitative analysis in medical imaging. In this paper, we formalize the phase retrieval as a regularized inverse problem, and propose a novel discretization scheme for the derivative operator based on B-spline calculus. The inverse problem is then solved by a constrained regularized weighted-norm algorithm (CRWN) which adopts the properties of B-spline and ensures a fast implementation. The method is evaluated with a tomographic dataset and differential phase contrast mammography data. We demonstrate that the proposed method is able to produce phase image with enhanced and higher soft tissue contrast compared to conventional absorption-based approach, which can potentially provide useful information to mammographic investigations. Abstract: Differential phase contrast imaging using grating interferometer is a promising alternative to conventional X-ray radiographic methods. It provides the absorption, differential phase and scattering information of the underlying sample simultaneously. Phase retrieval from the differential phase signal is an essential problem for quantitative analysis in medical imaging. In this paper, we formalize the phase retrieval as a regularized inverse problem, and propose a novel discretization scheme for the derivative operator based on B-spline calculus. The inverse problem is then solved by a constrained regularized weighted-norm algorithm (CRWN) which adopts the properties of B-spline and ensures a fast implementation. The method is evaluated with a tomographic dataset and differential phase contrast mammography data. We demonstrate that the proposed method is able to produce phase image with enhanced and higher soft tissue contrast compared to conventional absorption-based approach, which can potentially provide useful information to mammographic investigations.


Introduction
Differential phase contrast imaging using X-ray grating interferometer is a promising tool to revolutionize conventional radiography. With an incremental modification to the conventional X-ray imaging apparatuses, this technology is able to yield three different physical contrasts of the underlying sample in one scan: the conventional absorption contrast, differential phase contrast (DPC) and small-angle scattering contrast, resulting in much richer information than the traditional absorption-based imaging approaches [1][2][3][4]. Successful experiments have been demonstrated for mammography [5], hand imaging [6,7] and lung imaging [8,9] using a X-ray tube-based configuration.
The phase contrast, which is obtained by detecting the phase shifts of the X-ray waves when passing through the sample, has many advantages in imaging soft tissue (low-absorption materials) compared to the absorption contrast [3,10,11]. The optical properties of a tissue can be characterized by its complex refractive index n = 1 − δ + iβ . δ is the decrement of the real part of the refractive index responsible for the phase shift, whereas the imaginary part β describes the attenuation properties of the materials. At diagnostically relevant photon energies (i.e., between 10 and 150 keV), the phase shift for soft tissue plays a more prominent role than the attenuation because δ is typically three orders of magnitude larger than β . The phase shift induces refraction of the X-ray wave traveling through the tissue, and the grating interferometry is designed to detect the resulting overall refraction angle efficiently. In clinical applications such as mammography and computed tomography, when comparing various tissue samples, such relative differences in angular deviations are larger than the corresponding relative changes in intensity (related to the X-ray absorption). As a result, the phase contrast is expected to yield improved contrast when compared with conventional methods [12]. Moreover, for many tissue types, the phase differences drop less than their corresponding absorption differences specially for higher X-ray energies [12,13]. Therefore, phase imaging can be theoretically performed at higher energy, while still keeping the same contrast as the absorption-based approach, but with a lower dose deposition.
A potential shortcoming is that the grating interferometer does not measure the phase shift of X-rays directly, but only its first derivative. That is why this technique is called "differential" phase contrast imaging. Although the DPC image provides significant information on the edges of the sample structure, it doesn't allow quantitative analysis of the phase profile, for instance, giving the contrast between different tissues or comparing with absorption contrast directly. Theoretically, the phase retrieval problem can be solved by a trivial one-dimensional direct integration in the spatial domain, which is equivalent to dividing by spatial frequency in the Fourier domain. However in practice, the direct integration fails to generate phase image of satisfactory quality because the noise is also cumulated during the integration. The resulting phase image suffers from stripe artifacts. The problem becomes especially severe when the image size is very large and lower dose is intended to be used, which means lower signal-tonoise ratio (SNR), as is common in the medical imaging.
For phase retrieval in grating-based DPC imaging, two noteworthy methods are the bidirectional method [14] and the non-linear regularization method proposed by our group [15]. The bidirectional approach has shown to be well-suited for reducing stripe artifacts; however, it requires longer scanning time and possible higher dose deposition. Besides, two-scan approach requires one to rotate the gratings or the sample and to register the two images precisely, which is hard to be implemented in clinical situation. The non-linear regularization method is more practical because its benefits are actually from the mathematic tools; therefore no additional efforts (e.g. more image acquisitions or system modification) are needed, and it shows promising results [15]. Sperl et al showed that the regularization can also be done in Fourier domain [16].
Our major interests is to apply the differential phase contrast imaging to clinical applications, especially mammography. In mammography, the image is usually quite large, with a typical size of 6000 × 4000. Such an application requires fast phase retrieval algorithms which can deliver accurate and quantitative phase information for diagnostic purpose. The aim of this paper is to further develop an effective phase retrieval algorithm based on the regularization method [15] for clinical radiographic applications. The contributions of this paper are as follows: • The proposal of a more accurate discretization scheme for the directional derivative operator based on the B-spline calculus.
• The proposal of a better regularization strategy that penalizes discontinuities in all directions and also addresses the issue of the undetermination of the zero frequency part of the solution.
• The application of constrained regularization weighed norm (CRWN) algorithm in order to have more robust and faster phase retrieval scheme in comparison to [15]. We use the positivity and boundary condition to improve the performance of our retrieval algorithm.
• The experimental evaluation of the method using real data. We evaluated the algorithm with the phase contrast mammography data and directly compared the retrieved phase image with the absorption image. We demonstrated that the proposed method is able to produce high-quality phase image which can provide higher contrast of different tissues than conventional absorption-based image.

Methods
2.1. Phase retrieval in differential phase contrast imaging The phase retrieval problem in DPC imaging has been described thoroughly in [15]. Briefly, the grating interferometer measures the shift in the argument of a sine function ϕ(x, y) of two phase-stepping curves (PSC), from an object scan and a reference scan without object, for each detector pixel at the spatial position (x, y). ϕ(x, y) links to the refraction angle α(x, y) by where p 2 is the pitch of the analyzer grating, and d is the distance between the two gratings.
On the other hand, the refraction angle α(x, y) (namely, the DPC signal) can be expressed by [10] α(x, y) = λ 2π where λ is the wavelength of the x-ray photons and φ (x, y) is the target quantity (phase profile) that we want to obtain. By combing Eq. (1) and Eq. (2), the phase retrieval problem is expressed as where g = λ d p 2 is a constant that depends on the system design parameters. A direct integration then yields the phase image by: However, due to the presence of noise, the more realistic form of the phase retrieval is whereφ(x, y) denotes the actual measured signal and n(x, y) represents the noise. In most cases, the direct integration fails because of the noise accumulation and lack of knowledge of the boundary conditions.

Regularization-based phase retrieval method
Retrieving the phase from Eq. (5) can be considered as an inverse problem and be generally solved by regularized optimization [15]. Similar idea has also been demostrated for in-line phase contrast imaging [17]. That is to minimize a cost-function where D x is the derivative operator along the x-direction. The first term of the right side is the fidelity term and Ψ(φ ) is the regularization term. In . [15], we chose the regularization term to be Ψ(φ ) = λ D y φ l 1 (7) in order to suppress the strips along the y-direction which is perpendicular to the x-direction.
The l 1 norm is preferred because it generates more quantitative results compared to l 2 norm owing to the introduced sparsity in the gradient domain and enhancement of the edges. Benefiting from the regularization term, we showed that the regularization-based method suppressed the strip artifacts significantly. The resulting phase images provide improved contrast-to-noise ratio (CNR) and reveal more details of the sample.
To discretize the derivative operator D x (also D y ), we used the finite difference model, given by where i and j are the discrete coordinates (pixel coordinates) and N i is the image size in xdirection.
The discretization model in Eq. (8) is very simple. In practice, our major concerns are to get quantitative information and diagnosis-level high quality image. This requirement calls for an accurate discretization model, and a corresponding algorithm for solving Eq. (6) and handling large dataset effectively. In this paper, we replace Eq. (8) with a more reliable model, using polynomial B-splines and propose an efficient algorithm to solve Eq. (6) using the constrained regularized weighted-norm algorithm which will be described in the next two sections.

Discretization of the derivative operator using B-spline calculus
To formulate Eq. (5) as an inverse problem, a necessary step is to discretize the forward model. In this regard, we use a generalized sampling scheme and represent the object φ on the principal shift-invariant space whose generating function is the tensor product of two centered B-spline functions β (x), where k = (k 1 , k 2 ), x = (x 1 , x 2 ) and β (x) is the tensor product of two centered B-spline functions The polynomial B-splines are well-known to offer the best cost/quality trade-off among interpolators. Discretizing the forward model with B-splines has been used successfully in tomographic application [19,20]. Herein, we expand it to the radiographic case.
Since the mathematical model of our problem is based on the derivative along the direction This equation leads to a matrix formulation of Eq. (3), where c is a vector of B-spline coefficients in lexical order, ϕ ϕ ϕ is the measurement vector, and H is the system matrix with

Phase retrieving method
Now, we formulate the phase retrieving task as an inverse problem. We aim at finding c where where the sum is computed on all the B-spline coefficients and {Lc} i ∈ R 2 is the gradient vector of the image at position i. We constrained the solution through the convex set C . It can be the restriction on the positivity of the refractive index and the information of the boundary of the object.
The regularization is chosen such that to make the retrieving problem well-posed and to provide reconstructions without stripe artifacts. Since the phase retrieval problem is ill-posed, we use total variation (TV) regularization term to enhance the edges in the phase image. The enhancement of the edges reduces the stripe artifacts on the retrieved phase image. Since the null space of the imaging operator contains zero frequency, we also use Tikhonov regularization term. In order to be consistent with the discretization scheme, the discrete gradient operator is computed using the following Proposition. where In order to solve Eq. 14, we use constrained regularized weighted norm (CRWN) algorithm with total variation regularization which has been fully explained in [20]. We also choose parameters of the algorithm based on their suggestions. The TV parameter is λ 1 = 10 −3 ||g|| where ||g|| is the norm of measurement. The Tikhonov parameter should have small value; λ 2 = 10 −5 . In the following practical applications, we run the algorithm up to convergence. The convergence criteria is when the cost function varies slowly.

Retrieval of the object boundary
In DPC image, the background of the image is meant to be zero. It does not contribute to quantitative information but only to noise when retrieving the phase. Therefore delineating the boundary of the object and masking the background out can significantly reduce the noise in the phase retrieving procedure. One advantage of grating-based imaging is the access to the information of the absorption and phase simultaneously. In order to retrieve the boundary of the object, we use both information.
We apply Canny edge detection on both measurements (absorption and phase). As we are interested in the boundary of the object, we choose the first and the last detected edges position along each horizontal and vertical line inside the object. Afterwards, we apply median filter along each of the four boundary curves separately in order to remove some sudden jumps owing to not detecting some edges. We finally obtain two masks: one on the absorption measurements and the other on the phase information. We choose their intersection to mask the object. The convexity constraint that is applied to the algorithm is that the reconstructed image should be zero outside this mask.

Quantitative evaluation with tomographic dataset
We demonstrated in . [15] that the regularization-based method works well with simulated data. Here we choose the weight W as identity operator in the data fidelity term; we use typical 2norm. In this section we directly evaluated the proposed method using experimental data. As a first step, we evaluated the quantitativeness of the proposed method. To establish a ground-truth for comparison, we generated a phase projection from phase tomographic dataset of a plastic phantom by forward-projection. The phase tomographic dataset was acquired at the TOMCAT beamline of the Swiss Light Source (SLS) using a Talbot interferometer [21]. The relative bandwidth of the beam was about 10 −2 . The grating interferometer features a design energy of 25 keV and is operated at the third Talbot distance with G1-G2 distance of 121mm. The phase grating G1 has a pitch of 3.98 µm and the analyzer grating G2 has a pitch of 2.0 µm. The differential phase signal was obtained by taking the one-dimension derivative of the phase projection. The retrieved phase from this differential phase projection was then compared with the ground-truth and the results were shown in Fig. 1. By generating the differential phase signal in this way, we included the experimental noise in the image. Due to presence of the noise, quantitative recovery of the phase is very challenging, however as shown in Fig. 1(d), the profile comparison indicates that the agreement between the retrived phase and the ground-truth is fair.

Evaluation with mammographic data
In this section, we evaluated the potentials of the proposed method in medical imaging with mammographic data. The goal here is to show that the proposed method can yield phase images with useful diagnostic information, which are not accessible with conventional direct integration approach. All the data were obtained from the differential phase contrast mammography study conducted at Paul Scherrer Institut, Switzerland and Kantonsspital Baden, Switzerland [5,22] and the details of the experimental system and ethical issues can be found in [5]. The system consists of a Seifert ID 3000 tungsten target X-ray generator (operated at 40 kVp with mean energy of 28 keV and a current of 25 mA), a three-gratings interferometer (with pitches of p0=14 µm, p1=3.5 µm and p2=2.0 µm), and a flat panel detector featuring a 50 × 50 µm 2 pixel size. The interferometer was tuned at the 5th Talbot distance. The distance from the source to the sample grating was about 140 cm and the source-detector distance was 160 cm.

Biopsy sample experiment
The biopsy sample experiment was to evaluate how the proposed method works with high SNR signals. It is considered as a benchmark experiment which can evaluate the method without being too much intruded by the noise. The sample was a biopsy breast tissue with a carcinoma mass. It was first fixed in 4% formalin solution and then put into liquid paraffin for imaging in order to suppress the phase-wrapping effect at the tissue-air interface. A 32-step phase stepping scan was performed to generate a relative high quality image. In the whole breast sample experiment of next section, only 8 steps were used. The relationship between the number of steps and the noise variance in DPC image is linear [23]. Therefore the experiment here generated DPC images which were four times less noisy than in the ex-vivo whole breast experiment shown in section 3.2.2.
The results are shown in Fig. 2. Visually the DPC image in Fig. 2(b) showes the edges clearly; however, no quantitative comparison can be done from it. The direct integrated phase image in Fig. 2(c) contained obvious stripe artifacts. Although some structures of the sample can be distinguished due to the high SNR in this case, the artifacts still degraded the image significantly, making it impossible to interpret. On the contrary, the proposed method gave a strip-artifact-free phase image. The image had a similar appearance as the absorption image ( Fig. 2(a)) but with enhanced contrast and details.
Profile comparisons of the horizontal and vertical red lines in Fig. 2(a) was shown in the Figs. 2(e) and 2(e)(f), respectively. The contrast-to-noise ratio (CNR) is used as figure of merit: where I s and I b were the mean values of selected ROIs in Figs. 2(e) and 2(f) which represent the signal and background, respectively, and σ b is the standard deviation of the background ROI. Note that the scales of the absorption and phase images are different since they are completely different physical quantities. It is clear that the phase image gives a higher contrast compared to absorption image. Quantitatively, for the horizontal line (Fig. 2(e)), the phase CNR is 28 while the absorption CNR is 13; and for the vertical line ( Fig. 2(f)), the phase CNR is 25 while the absorption CNR is 7.
The results show that the phase image indeed can provide superior contrast than the absorption image. Clinically it would help increase the diagnostic accuracy.

Ex-vivo whole breast experiment
It is worth mentioning that the whole breast DPC image was formed by stitching several smaller acquisitions together due to the limited field of view (FOV) of the current system. Stitching ar-tifacts were hard to avoid especially because of the sample deformation happened from time to time. Moreover, the background of each acquisition was not uniform because of the systemic drifts during the scans. Therefore the whole breast DPC data was "inconsistent" in the sense that not only noise corrupted the images but also the inhomogeneity of each block. Considering those factors as well as the large image size, the direct integrated phase image was uninterpretable for radiologists. The structure features were concealed in severe strip artifacts as the common case showed in Fig. 3(c); therefore no useful clinical information was available from the integrated phase image. A phase retrieval example of a whole breast is shown in Fig. 3. The sample selected here contained a large tumor mass and many spiculations. Spiculations were strong indications of the existence of malignant mass. For this particular sample, we want to explore what we can gain from the phase image compared to the absorption image, therefore they are compared directly. We applied the new algorithm to the whole image as well as a selected ROI, located around the carcinoma.
The stitching artifacts are clearly seen in Fig. 3(b) and are especially obvious at the background region. Those artifacts also cause inhomogeneous background in the retrieved phase. To suppress this effect, a background mask was generated using the absorption image as described in Section.2.5 and applied to the retrieved phase to create a clear background as showed in Fig.  3(d). To the best of our knowledge, this is the first time ever that an image of the line integral of the phase signal from a whole breast sample obtained with a conventional X-ray tube has been presented.
The results of the selected ROI are shown in Fig. 4. In principle, when operating on an ROI, the boundary conditions and thus the starting wave front profile φ (x = 0, y) for the integration are unknown. The loss of the boundary information will worsen the strip artifacts and even cause shadow artifacts. An equivalent effect is the "phase wrapping" happening at the border between the skin and air. The refraction angle there is so large that the induced fringe shift is larger than one detector pixel and therefore it cannot be measured by the grating interferometer. This effect also causes the loss of information at the boundaries. Usually the phase wrapping effect can be avoid by immersing the breast in a liquid solution. However, this option is not feasible in practice and, at the end, it will also reduce the contrast in the absorption image. Nevertheless the resulting phase image (Fig. 4(b)) is quite good and the spiculations are better seen. Three lines (indicated in Fig. 4) across the spiculations were selected and their profiles were given in Figs. 4(c) and 4(d). The profiles clearly showes that the phase image provides much higher contrast of the spiculations compared to the absorption image.
It is worth mentioning that the contrast of the calcifications is higher in the absorption image than in the phase image (see Figs. 4(a) and 4(b)). Detecting microcalcifications is important for early breast cancer screening. To this extent, the phase image will not be able to replace the absorption image completely but be a very useful complement. The third contrast, scattering contrast, on the other hand could play an important role in microcalcification detection.
The parameters only need to be be adjusted once for a given protocol and the algorithm can then be run on the whole data set.

Comparison with former method
We also compared the new algorithm with our former algorithm proposed in . [15] and the comparison are shown in Figs. 4(b) and 4(c). The corresponding regularization parameter was λ = 10 −3 and this parameter was decided by trial and error to find the best result. Visually the new algorithm gave a superior results. The details were more visible due to fewer artifacts.

Conclusion
Unlike DPC tomography where the phase information can be recovered effectively by reconstruction algorithm, retrieving phase image from DPC projection remains challenging and limits the exploration of the advantages of the phase information in radiographic applications.
In this paper, we proposed an iterative phase retrieval algorithm for differential phase contrast imaging. The algorithm utilized a novel discretization model using B-spline calculus and a constrained regularized weighted-norm algorithm is used to solve the inverse problem and ensure a fast implementation. The algorithm is evaluated with breast biopsy and mastectomy samples. The present study demonstrated that DPC signal is indeed capable of providing higher contrast in clinical relevant features, like the spiculation etc. These results could help to improve breast cancer diagnosis.