Joint absorption and phase retrieval in grating-based x-ray radiography

Given the raw absorption and differential phase-contrast images obtained from a grating-based x-ray radiography, we formulate the joint denoising of the absorption image and retrieval of the non-differential phase image as a regularized inverse problem. The choice of the regularizer is driven by the existing correlation between absorption and differential phase; it leads to the linear combination of a total-variation norm with a total-variation nuclear norm. We then develop the corresponding algorithm to efficiently solve this inverse problem. We evaluate our method using different experiments, including mammography data. We conclude that our method provides useful information in the context of mammography screening and diagnosis. © 2016 Optical Society of America OCIS codes: (100.3190) Inverse problems; (100.3010) Image reconstruction technique; (340.7440) X-ray imaging; (110.6960) Tomography. References and links 1. A. Momose, S. Kawamoto, I. Koyama, Y. Hamaishi, K. Takai, and Y. Suzuki, “Demonstration of x-ray Talbot interferometry,” J. Appl. Phys. 42(7 B), 866–868 (2003). 2. C. Kottler, C. David, F. Pfeiffer, and O. Bunk, “Two-directional approach for grating based differential phase contrast imaging using hard x-rays,” Opt. Express 15(3), 1175-1181 (2007). 3. T. Thüring, P. Modregger, B. R. Pinzer, Z. Wang, and M. Stampanoni, “Non-linear regularized phase retrieval for unidirectional X-ray differential phase contrast radiography,” Opt. Express 19(25), 25545-25558 (2011). 4. J. I. Sperl, D. Bequé, G. P. Kudielka, K. Mahdi, P. M. Edic, and C. Cozzini, “A Fourier-domain algorithm for total-variation regularized phase retrieval in differential x-ray phase contrast imaging,” Opt. Express 22(1), 450– 462 (2014). 5. M. Nilchian, Z. Wang, T. Thuering, M. Unser, and M. Stampanoni, “Spline based iterative phase retrieval algorithm for x-Ray differential phase contrast radiography,” Opt. Express 23(8), 10631-10642 (2015). 6. S. Lefkimmiatis, A. Roussos, M. Unser, and P. Maragos, “Convex generalizations of total variation based on the structure tensor with applications to inverse problems,” Proceedings of the Fourth International Conference on Scale Space and Variational Methods in Computer Vision (SSVM’13), Seggauberg, Republic of Austria, June 2-6, 48–60 (2013). 7. M. Stampanoni, Z. Wang, T. Thüring, C. David, E. Roessl, M. Trippel, R. A. Kubik-Huch, G. Singer, M. K. Hohl, and N. Hauser, “The first analysis and clinical evaluation of native breast tissue using differential phase-contrast mammography,” Invest. Radiol. 46(12), 801–806 (2011). 8. T. Weitkamp, A. Diaz, C. David, F. Pfeiffer, M. Stampanoni, P. Cloetens, and E. Ziegler, “X-ray phase imaging with a grating interferometer,” Opt. Express 13(16), 6296–6304 (2005). 9. S. A. McDonald, F. Marone, C. Hintermueller, G. Mikuljan, C. David, F. Pfeiffer, and M. Stampanoni, “Advanced phase-contrast imaging using a grating interferometer,” J. Synchrotron Radiat. 16, 562–572 (2009). #258837 Received 3 Feb 2016; revised 10 Mar 2016; accepted 11 Mar 2016; published 25 Mar 2016 © 2016 OSA 4 Apr 2016 | Vol. 24, No. 7 | DOI:10.1364/OE.24.007253 | OPTICS EXPRESS 7253 10. H. Meijering, J. Niessen, and A. Viergever, “Quantitative evaluation of convolution-based methods for medical image interpolation,” Med. Imag. Anal. 5, 111–126 (2001). 11. S. Ramani, and J. A. Fessler.,“ A splitting-based iterative algorithm for accelerated statistical x-ray CT reconstruction,” IEEE Trans. Med. Imag. 31(3), 677-688 (2012). 12. D. L. Donoho,“ De-noising by soft-thresholding,” IEEE Trans. Inf. Theo. 41(3), 613–627 (1995).


Introduction
Grating-based x-ray imaging is a powerful modality to investigate biological samples. It measures three complementary characteristics of the imaged sample: the conventional absorption contrast (AC), the differential phase contrast (DPC), and a small-angle scattering contrast [1]. The phase contrast (PC) is a related quantity of interest that must be retrieved from DPC, possibly through an integration process.
In this paper, we take advantage of the high correlation between PC (where the contrast between features is pronounced) and AC (which provides structural details) to simultaneously retrieve PC and denoise AC. To the best of our knowledge, our work is the first to jointly retrieve the phase from the differential-phase image and to attenuate noise in the absorption image.
By employing a grating interferometer, one gains direct access to the spatial derivative of the phase φ , where x ∈ R 2 specifies spatial coordinates. The simplest approach to retrieve φ from g is line integration. Unfortunately, the phase image thus recovered exhibits artifacts due to noise amplification at low frequencies, as shown in Fig. 1. These problems are exacerbated in the low-dose regime typically necessitated by biological tissues, which worsens the signalto-noise ratio. At least two advanced phase-retrieval approaches have been proposed in the literature. The first one, called the bidirectional method, is able to reduce stripe artifacts at the cost of longer scanning times and high radiation doses; it needs the precise registration of two images, along with the mechanical rotation of gratings, two requirements that are hard to achieve [2]. Meanwhile, the second method, called nonlinear regularization, is easier to deploy because no mod-ification of the practical setup is required. Although it shows promising results when a totalvariation (TV) regularization is applied either in the spatial domain [3] or in the frequency domain [4], it sometimes fails to achieve a satisfying attenuation of the stripe artifacts and occasionally results in an insufficient contrast between regions.
Recently, we developed a spline-based iterative phase-retrieval algorithm along with TV regularization [5]. The experimental results show improvement in the quality of the retrieved phase compared to the state-of-the-art techniques. This work aims at further improving upon the stateof-the-art [5] by exploiting the correlation between DPC and AC. We propose here an inverseproblem approach inspired by nonlinear regularization; more precisely, we complement TV by another regularizer that leverages the correlation between DPC and AC. Thus, the two main contributions of this paper are • joint denoising and retrieval of the absorption and phase contrasts, respectively; • simultaneous use of a total-variation regularizer and a Schatten-norm Jacobian regularizer to attenuate reconstruction artifacts.
We develop an alternating-direction method of multipliers to achieve our goals, guided by B-spline calculus to finesse the necessary discretization schemes. Finally, we evaluate the performance of the proposed method using three different experiments on the synthetic and real data. Specifically, we illustrate the application of our method to breast tissue from mammography data.

Forward model
Letting ρ denote the ideal AC image, our problem is to retrieve an approximation of the attenuation ρ : R 2 → R + and the phase φ : R 2 → [0, 2π] from the forward model where x ∈ R 2 and n 1 and n 2 represent the AC and DPC noise that is present in the measurements g 1 and g 2 , respectively. For convenience, we discretize the forward model as where β n is the tensor product of two centered B-splines of integer degree n. (Polynomial B-splines are known to offer the best cost/quality tradeoff among many interpolators [10].) Then, the forward discretization of Eq. (2) is where g is the vector of measurements, and There, B encodes the interpolation process and D x 2 is the directional (by convention, vertical) derivative operator. Moreover, c = (c 1 , c 2 ) collects in a single column vector each B-spline coefficient of Eq. (3) and Eq. (4). We say that c is the discrete representation of the object. To be consistent with Eq. (2), the directional derivative operator must be such that which amounts to where p = (p 1 , p 2 ). There, [·] p specifies the pth entry of a vector.

Model fitting
We formulate the joint-retrieval problem of c in the framework of a penalized least-square estimation. Our goal is to find the minimizer of where Ψ is a regularizing term made of several components. Letting the first component be Ψ 1 , we promote the lessening of stripe artifacts by encouraging sparsity in the directional gradient of the phase image. We define As a second, independent constraint, we choose the nuclear total variation, which is a vectorial extension of TV [6]. Its purpose is to strengthen the correlation of AC and DPC over edges. It is given by the point-wise sum of Schatten 1-norms of the Jacobian matrix J and is expressed as The Schatten 1-norm · S 1 is the 1 norm of the singular values of its matrix argument. Note that the Schatten 1-norm is a convex functional. The discrete Jacobian operator J maps R N to R N×2×2 , where N is the total number of coefficients within one image. All directional derivatives are computed in accordance with the B-spline calculus suggested by Eq. (8).

Optimization
We aim at determiningĉ such that We reformulate this unconstrained optimization problem as a problem constrained over the auxiliary variables u ∈ R N and v ∈ R N×2×2 such that The scaled form of its augmented-Lagrangian formulation is We use the standard alternating-direction method of multipliers (ADMM) [5,11] to minimize Eq. (14). This results in The first minimization step is a quadratic problem. Its gradient is given by where the superscript symbol * indicates the adjoint operator. We take advantage of a conjugategradient approach to solve this first step. The second step of Eq. (15) is a minimization in terms of the auxiliary variable u. It takes the form Therefore, the solution involves a point-wise soft-thresholding operation (i.e., the proximal map of the 1 norm), which leads to The third step of Eq. (15) is where {·} m ∈ R 2×2 is the mth (2 × 2) matrix. Solving Eq. (20) requires the proximal map of the Schatten 1-norm It has been shown in [6] that the solution of Eq. (21) can be obtained by considering the singularvalue decomposition (SVD) z = P Λ Q and by rewriting Eq. (21) as where prox 1 ,λ is a soft-thresholding operator [12]. Since z ∈ R 2×2 , the determination of SVD is easy. Consequently, the third minimization step of Eq. (15) is summarized as where Λ m is the singular-value matrix of J c k+1 + d k 2 m . The remaining two steps of Eq. (15) are updates of the augmented-Lagrangian parameters. More precisely, λ 1 is chosen proportional to the variance of the phase measurement and λ 2 is adjusted according to the side information of how much correlated phase and absorption are. The parameters µ 1 and µ 2 are fixed to 1.

Experimental results
We conduct three different experiments to evaluate the performance of the proposed joint phaseand-absorption retrieval. In all experiments, we used cubic B-splines as they offer the best cost/quality tradeoff. First, we assess the proposed technique with a synthetic phantom that provides the ground-truth for the quantitative analysis. Second, we consider a physical phantom to measure the efficiency of the technique in reducing the artifacts introduced by the beam instabilities. Third, we apply it on mammography data to evaluate its performance compared to the state-of-the-art technique.

Quantitative analysis with synthetic data
The structure of one synthetic phantom is shown in Fig. 2. The sample is a cube that is 400 µm wide. It contains a sphere and a diamond-shaped structure. The imaginary and real part of their refractive indices are β = 4.4 × 10 −9 and δ = 1 × 10 −7 , respectively. These values are compatible with the refractive index of biological or medical samples.
We simulate the imaging system that produces two-dimensional, AC and DPC radiographs from the three-dimensional synthetic phantom. We consider an X-ray energy of 25 keV. The setup parameters of the simulation are given in Table 1.  The absorption grating is located at the second-order fractional Talbot distance to the π-shift phase grating, z = 1.5g 2 1 /4λ . In order to capture phase and absorption measurements, we use the phase-stepping technique that has been described in [8]. One period with nine scanned phase steps of the phase-stepping curve is captured. The Fourier-component analysis [1] is used to retrieve the AC and DPC measurements which are denoted by g 1 and g 2 , respectively. To simulate the noise of the detector, we add a Gaussian noise whose variance is proportional to the pixel-wise intensity. The measurements are depicted in Fig. 4.
Before attempting to evaluate the performance of the proposed phase and absorption retrieval technique, we first validate the forward imaging model Eq. (1) by computing the expected measurements asg where k is the wavenumber and g 1 and φ are the line integrals of the imaginary and real part of the refractive index of the sample, respectively. They are shown in Fig. 3. We then plot the profile of the measurements located in the middle of the phase and absorption measurements

Results
We apply the proposed technique on the simulated measurements described the Section 3.1.1 to retrieve the phase and the absorption. During the course of the iterations, we also apply a positivity constraint on the solution. The optimal regularization parameters are λ 1 = 0 and λ 2 = 0.02. Yet, we also found that the performance of the retrieval framework is quite stable with respect to the regularization parameter λ 2 ([10 −3 , 10 −1 ]). The other methods that were considered in the comparison are (e) (f) Fig. 6. The first and second columns are the retrieved phases and absorptions, respectively. The first row is not regularized. The second row is the phase retrieval using [5]. The third row is the joint phase-and-absorption retrieval using the proposed technique.
1. Retrieval of the phase and the absorption without using any regularization and only by minimizing the least square error.
2. State-of-the-art method for retrieving the phase information using total-variation regularization [5]. The optimal value for the total-variation parameter is λ = 0.01.
The results of these three techniques are depicted in Fig. 6.
Remarkably, the proposed technique jointly retrieves the phase and the absorption information with high quality. To quantitatively analyze its performance, we plot the line profile of the retrieved images versus the ground truth as shown in Fig. 7. The results suggest that the proposed technique performs well.

Data
To establish a physical ground-truth, we generated a forward-projected phase from the tomographic dataset of a plastic phantom. The tomographic dataset was acquired at the TOMCAT beamline of the Swiss Light Source (SLS) using a Talbot interferometer [9]. The relative bandwidth of the beam was 10 −2 . The grating interferometer features a design energy of 25 keV and is operated at the third Talbot distance with a 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.  Fig. 7. The first and second columns are the retrieved phases and absorptions, respectively. The first row is the line profile for the retrieved images without any regularization. The second row is the line profile for the retrieved phase using [5] and absorption. The third row is the line profile of the joint phase-and-absorption retrieval using the proposed technique.

Results
As the beam instabilities introduce line artifacts on the projection measurements, it is important to remove them. We applied our method to the physical phantom, for which the line artifacts are really disturbing (Fig. 8). The results are shown in Figure 9. They confirm that our technique reduces the line artifacts significantly compared to other approaches.

Mammography
We now present an example of real-world data where the proposed method yields absorption and phase images that contain useful diagnostic information. In particular, we want to explore the extent to which AC and DPC are complementary and/or correlated.  Fig. 9. The first and second columns are the retrieved phase and the denoised absorption contrast using different techniques, respectively. The first row is the least-square solution. The second row applies a state-of-the-art technique to retrieve the phase [5]. The third row is the solution of the proposed joint phase-and-absorption retrieval.
clearly hint at the presence of a malignant mass. Individual acquisitions had a limited field of view; this was counterbalanced by an imaging stitching process. Unfortunately, not only did stitching artifacts arise because of deformations in the breast tissue, but also systemic drifts during the scan resulted in a nonuniform background. This explains why the naive integration of the DPC image produces the severe stripe artifacts visible in Fig. 1.   Fig. 10. Top-left: denoised absorption image whose location corresponds to the ROI in Fig. 1. Top-right: retrieved phase. Bottom-left: multichannel visualization. Bottom-right: caption.

Results
We applied the proposed method to data that contain a carcinoma, as revealed in Fig. 10. Because of the limited field of view, the integration process for a line can be determined only up to the (unknown) integration constant that specifies the boundary conditions of the wavefront φ (x)| x 1 =0 . However, one of the roles of our regularization is to promote a consistent choice of this constant across lines, thus avoiding the worsening of stripe and shadow artifacts. We see in Fig. 10.b that the PC image we retrieved out of DPC is quite good; the spiculations are easier to recognize there than in the AC image shown in Figure 10.a.
At the same time, however, the micro-calcifications are better contrasted in AC than in PC. (The early detection of micro-calcifications is an important tool for oncologists.) This difference in the quality of contrast is to be expected since the complex refractive index of calcium is similar to that of soft tissue, thus leading to an absence of contrast in the phase image. It then follows that AC and PC offer truly complementary modalities. We have fused AC and PC in Fig. 10.c, where AC corresponds to the blue channel and PC to the red channel of a color representation.

Comparison to the state-of-the-art
We now focus on the region of interest (ROI) delineated in Fig. 10. This allows for a detailed comparison of the outcome of our method with that of [3]. We present the corresponding visual results in Fig. 11. In the PC image, we observe that the state-of-the-art method (top right) fails to remove the stripe artifacts, while our proposed method (bottom right) is clearly more Fig. 11. Top-left: absorption image whose location corresponds to the ROI in Fig. 10. Topright: phase retrieved by [3]. Bottom-left: absorption image after our joint denoising-andretrieval process. Bottom-right: phase image that we retrieved from the differential phase.
successful. In addition, a substantial amount of noise has disappeared from our version of the AC image (bottom left), at no cost in the visibility of structural details.

Conclusion
Grating-based x-ray imaging provides simultaneously absorption and differential-phase images. Thus, the (integrated) phase image-which is the quantity of interest-must be retrieved computationally. This retrieval is challenging, which has limited the use of the phase in practical radiographic applications. In this paper, we propose a new iterative method to retrieve the phase image and to jointly denoise the absorption image. The present study is a proof of concept that demonstrates that it is indeed possible to improve the clarity of phase features that have a clinical relevance (e.g., spiculations). At the same time, the contrast of complementary features in the absorption image (e.g., micro-calcifications) is improved, too. In addition, the experimental results demonstrate that the proposed technique is quantitative and that it reduces the line artifacts introduced by the beam instabilities. These foundings suggest that a wider deployment of joint absorption-and-phase retrieval could have a beneficial clinical impact.