for x-ray differential phase-contrast tomography

Purpose: The purpose of this work is to propose a cost function with regularization to iteratively reconstruct attenuation, phase, and scatter images simultaneously from di ff erential phase contrast (DPC) acquisitions, without the need of phase retrieval, and examine its properties. Furthermore this reconstruction method is applied to an acquisition pattern that is suitable for a DPC tomographic system with continuously rotating gantry (sliding window acquisition), overcoming the severe smearing in noniterative reconstruction. Methods: We derive a penalized maximum likelihood reconstruction algorithm to directly reconstruct attenuation, phase, and scatter image from the measured detector values of a DPC acquisition. The proposed penalty comprises, for each of the three images, an independent smoothing prior. Image quality of the proposed reconstruction is compared to images generated with FBP and iterative reconstruction after phase retrieval. Furthermore, the influence between the priors is analyzed. Finally, the proposed reconstruction algorithm is applied to experimental sliding window data acquired at a synchrotron and results are compared to reconstructions based on phase retrieval. Results: The results show that the proposed algorithm significantly increases image quality in comparison to reconstructions based on phase retrieval. No significant mutual influence between the proposed independent priors could be observed. Further it could be illustrated that the iterative reconstruction of a sliding window acquisition results in images with substantially reduced smearing artifacts. Conclusions: Although the proposed cost function is inherently nonconvex, it can be used to reconstruct images with less aliasing artifacts and less streak artifacts than reconstruction methods based on phase retrieval. Furthermore, the proposed method can be used to reconstruct images of sliding window acquisitions with negligible smearing artifacts. C 2016 dx.doi.org


INTRODUCTION
Iterative reconstruction (IR) in general and statistical IR, in particular, is currently one of the hot topics in tomographic xray imaging. The objectives and applications of IR algorithms include dose reduction, sparse sampling, and limited angle tomography.
Additional potential for x-ray tomography lies in phase contrast imaging. This became recently feasible even with a conventional x-ray tube by using a Talbot-Lau interferometer 1 and is known as grating based differential phase-contrast imaging (DPCI). This technique shows excellent soft-tissue contrast 2,3 that would likely add diagnostic value if applicable to clinical imaging.
DPCI is based on a conventional x-ray setup, which is extended by three grids: one grid in front of the tube, one grid in-between tube and detector, and an analyzer grid directly in front of the detector. 1 With such a DPCI setup it is possible to reconstruct besides the x-ray attenuation properties, also the small-angle scattering properties and the electron density properties (causing differential phase shift) of an object. The reconstructed images will be called shortly attenuation image, scatter image, and phase image. To extract these object properties, it is necessary to do x-ray measurements for different positions of the analyzer grid. In conventional setups, the analyzer grid is stepped through a number of different positions for each acquisition view (one view corresponds to one angular position of tube and detector). This acquisition method is called phase stepping (see Fig. 1, left). Before reconstruction, the so called phase retrieval is performed, converting the measured x-ray projections for the different analyzer grid positions to projections of attenuation, scatter, and phase. Prerequisite for the phase retrieval is the availability of x-ray measurements for a number of different analyzer grid positions for each acquisition view, i.e., phase stepping has been done for each acquisition view. Finally, the projections of attenuation, scatter, and phase are reconstructed with conventional means, like FBP.
A number of publications address the combination of DPCI with (statistical) IR. Most of them first extract attenuation, scatter, and phase projections with phase retrieval, and then reconstruct attenuation, scatter, and phase images with iterative methods from these projections. [4][5][6][7][8][9][10][11][12] Due to the phase retrieval step, these methods rely on data acquisition with phase stepping.
So far there is one publication proposing an statistical IR approach that reconstructs images directly from the measured x-ray intensities, 13 but without any further evaluation of the algorithms performance on larger CT data sets. Such an approach will be called intensity based statistical IR (IBIR) in the following. IBIR has the advantage that the phase retrieval, which is prone to errors due to phase wrapping, is not necessary, if no start image for the phase is needed. Furthermore IBIR can also be used to reconstruct unconventional acquisition patterns, especially it is not necessary to do phase stepping. This will become clear below, comparing conventional reconstruction and IBIR for a sliding window acquisition, 14 where phase stepping is done over subsequent views, i.e., for each view only one phase step is available (see Fig. 1, right). Another advantage of IBIR is that the noise model is needed for the measured intensities, for which the noise modeling can be done most accurately, while for the other statistical approaches the noise model is needed for the retrieved attenuation, scatter, and phase projections.
One disadvantage of IBIR is that the cost function to minimize cannot be guaranteed to be convex or even to have only one minimum due to the occurrence of phase wrapping. However, one should keep in mind that all reconstruction methods for DPCI suffer from this problem and cannot guarantee one unambiguous solution. Actually, all available reconstruction methods based on a phase retrieval step are more prone to artifacts due to phase wrapping, since they do not aim for a consistency of the phase image with the originally measured data (in contrast to IBIR). Furthermore, in IBIR the application of a penalty term can relax the problems of nonconvexity and multiple minima due to phase ambiguity, which is again not possible for other available methods.
In Ref. 13 a maximum likelihood (ML) method for IBIR is proposed, however without a penalty term. The utilized optimization is based on a heuristic method, which may diverge if the hand-tuned parameters are not set properly. Start images for optimization are constructed manually based on knowledge about the shape of the underlying simulated phantom. Tests are done only for a comparably simple simulated phantom assuming a noise free acquisition.
In this paper, we evaluate a cost function with penalty term. The applied optimization method is guaranteed not to diverge, and the dependence on different start images is analyzed. The whole method is evaluated on experimental data.

METHOD
A penalized maximum likelihood image reconstruction is based on the minimization of a cost function. The minimization determines the image such that it most likely belongs to the input data, given some a priori information (the penalty). 15 Thus, the cost function has two terms: A data term, evaluating the likelihood that the image belongs to the input data and a penalty term evaluating to which degree the image fulfills the a priori information. To calculate the data term a forward model is necessary to estimate which input data would belong to a given image. Here the following forward F. 1. Illustration of two data acquisition schemes for differential phase contrast imaging (DPCI). In the conventional scheme on the left, for each angular position a fixed number of phase steps is executed, allowing straight forward retrieval of attenuation, scatter, and phase projections. In the sliding window scheme on the right, the phase stepping is not done for each angular position separately, but during the rotation of the acquisition system. For each angular position only one phase step is available, and a retrieval of attenuation, scatter, and phase projections is only possible if the data for the missing phase steps are interpolated. model is used: 13 where i indexes the measurements (for each detector element in each phase step in each view). N i 0 is the expected photon count in the primary beam, V i 0 the visibility, and Φ i 0 the phase offset for measurement i. l ai , l si , and Φ i are the line integrals of attenuation, small-angle scattering, and differential phase, respectively.N i is the expected photon count for measurement i. The line integrals are calculated from the discretized images as follows: where A is a conventional forward projector for CT, and A δ is a forward projector with a differentiating component. 4 where N i are the measured photon counts. Constants that are irrelevant for the minimization are omitted here. For the penalty three independent terms enforcing smoothness 15 are applied to the three images, where j runs over all voxels (it is assumed that all three images are defined on the same grid), k runs over neighborhoods S a j , S s j , S φ j around the voxel j in the attenuation, scatter, and phase image. w a j k , w s j k , w φ j k are weights to adjust the influence of the voxels within a neighborhood. Ψ(·) is the Huber potential function. 15 β a , β s , and β φ are three factors to adjust the influence of the three penalty terms in the cost function. The complete cost function is given by C = L + R, which has to be minimized with respect to the three image vectors µ, ϵ, and δ.
It is the first aim of this work to analyze and discuss the properties and potential of the introduced cost function, not the performance and properties of applicable optimization methods. Thus, it is acknowledged that the optimization F. 2. Images reconstructed from a synchrotron DPCI acquisition of a mouse. Conventional iterative reconstruction (Conv. IR) and intensity based iterative reconstruction (IBIR) with different start images. Columns from left to right: Conventional iterative reconstructions with FBP start images/IBIR reconstructions (100 iterations) with FBP start images/IBIR reconstruction (100 iterations) with smoothed FBP start images/IBIR reconstruction (900 iterations) with zero valued start images. Rows from top to bottom: Attenuation images/Scatter images/Phase images. method introduced briefly here, is not cutting edge, and there may exist much better suited algorithms. However, the method is comparatively simple to implement and parameterize and turned out to minimize the cost function reliably and much faster than nonlinear conjugate gradient (CG). It is a steepest decent method, which alternately optimizes the three images. For this, the gradient of the cost function with respect to all voxels of one of the three image is calculated. Then, a 1D search using the golden section method is performed to find the minimum of the cost function along the direction of this gradient, and the image is updated accordingly. Subsequently the processing continues with the next image. One optimization step is finished, when all three images have been updated. The optimization steps are repeated until convergence is reached. The method is guaranteed to decrease (or in convergence not increase) the value of the cost function in each step. In the case a start image is used, about 100 iterations are sufficient to achieve an acceptable convergence. Increasing the number of iterations does not visibly change the images.
For the evaluation, IBIR reconstructions are compared to "conventional" iterative reconstructions (conventional IR) obtained using data after phase retrieval. In this case, the variances of the data after retrieval are calculated using the variance estimates as derived by Weber et al. 16 and the cost function for each reconstruction problem was set to a weighted least squares data term in addition to the regularization term as used in the proposed IBIR algorithm. The minimization problem was then solved using a SPS-type algorithm. 15,17 Two different data acquisition schemes will be compared in the following. Both are based on a circular CT trajectory. In the first scheme (Fig. 1, left), for each angular position a fixed number of phase steps is executed, allowing straight forward retrieval of attenuation, scatter, and phase projections for a conventional reconstruction approach. In the second scheme, the so called sliding window acquisition, 14 the phase stepping is not done for each angular position separately, but during the rotation of the acquisition system (Fig. 1, right). Thus, for each angular position only one phase step is available, and a retrieval of attenuation, scatter, and phase projections is only possible if the data for the missing phase steps are interpolated.

RESULTS
For the analysis of the properties of the introduced cost function, data of a mouse acquired at a synchrotron are used. The regularization strength parameters β a , β s , and β φ are tuned manually to achieve a good compromise between a low noise level and at the same time a high level of detail. It should be noted that the question, what a "good compromise" is, depends strongly on the task that should be performed on the images, and on the preferences of the observer. This is the reason, why for many commercially available IR products for conventional CT the regularization strength parameter can be adjusted by the user. Thus, the adjustment of the regularization strength is also here subjective. The delta parameters of the Huber potential functions are set to very low values as compared to the image contrasts, such that essentially an absolute value function is mimicked. More details about how to adjust regularization parameters for this example will be given below.
For testing purposes, the optimization is initialized with three different sets of start images: • conventional reconstructions, based on retrieval of attenuation, phase, and scatter projections with subsequent application of filtered backprojection (FBP); • smoothed versions of the conventionally reconstructed images (smoothing with a Gaussian kernel); • images with all voxels set to zero.
As can be seen in Fig. 2, the IBIR reconstructions for the different start images converge to the same images, indicating that the nonconvexity is not a major problem for this example. Nevertheless it should be noted that convergence is substantially slower if a zero start image is used. While only 100 iterations are needed for the FBP start image with or without smoothing, 900 iterations are necessary for a zero start image. Especially the low frequency structures in the phase image converge slowly.
For the comparison of conventional IR and IBIR, the regularization strength for conventional IBIR was set such that the image sharpness matches the IBIR results. The results in Figs. 2 and 3 illustrate the advantages of IBIR with respect to image quality: Attenuation images are very similar, but it is possible to achieve less aliasing artifacts and less streak artifacts in the phase image and scatter image reconstructed with IBIR, while keeping the level of detail in the images.
Apart from dependence on start images and image quality, it is tested for each regularization term, how it influences all three images, since it can be assumed that, e.g., the regularization term for the scatter image also impacts the phase image, since both are connected in the forward model of the data term, Eq. (1). For the test, each regularization term is "activated" individually, by only setting one of the three β-parameters to a nonzero value. Results are given in Figs. 4 and 5. While Fig. 4 visualizes the reconstructed images for different regularization settings, Fig. 5 shows difference images illustrating the difference between reconstructions activating one of the regularization terms and a reconstruction without any regularization. It can be seen that each regularization term has the strongest impact on the image it addresses directly, while the impact on the other images is weaker. Thus, it was in this example possible to tune the appearance of each image individually by adjusting the corresponding regularization parameter.
F. 5. Difference between reconstructions activating one of the regularization terms and a reconstruction without any regularization (as given in Fig. 4), to illustrate the influence of the regularization terms. Left: Difference for "Attenuation only" / Middle: Difference for "Scatter only" / Right: Difference for "Phase only." Finally, the performance of IBIR with respect to the sliding window acquisition is tested. For this a data set of a heart of a rat is used. The acquisition was done with a circular 180 • trajectory, 751 views per rotation, and five phase steps per view. Details about the acquisition can be found in Ref. 14. In Fig. 6 a detail of the phase image is shown for different reconstruction techniques. Attenuation image and scatter image are skipped here, since they show very low contrast for soft tissue. As a reference the conventional FBP reconstruction for the full phase stepping acquisition with 5 phase steps is given. The sliding window data set is generated by simply keeping only one of five phase steps for each view, running cyclically through the phase steps with increasing view number. For the FBP reconstruction from the sliding window data set, five subsequent projections are used for attenuation, scatter, and phase retrieval ignoring the rotation of the system during acquisition. The FBP reconstruction result (Fig. 6) is affected by angular smoothing (see white arrows in the images), which is especially apparent in the outer area far away from the rotation axis. The IBIR reconstruction utilizing only the available data of the sliding window data set is initialized with the images of the FBP reconstruction from the sliding window data set. One can observe that the angular smoothing is removed by IBIR (see again the white arrows), and the resulting image is as detailed as the one from the conventional, fully sampled FBP reconstruction. Here, IBIR has a clear advantage over reconstruction methods which rely on the retrieval of attenuation, scatter, and phase in projection domain, regardless if they are iterative or noniterative.

SUMMARY
A penalized maximum likelihood cost function for IBIR from tomographic DPC acquisitions has been introduced, allowing the simultaneous reconstruction of attenuation, scatter, and phase images directly from the measured intensities. Thus, the cost function allows reconstruction without phase retrieval and is especially suited to reconstruct from acquisitions with incomplete phase stepping data. Thus, new acquisition protocols, even without phase stepping become feasible.
The properties of the cost function have been investigated based on measured synchrotron data of a mouse. It turned out, that convergence to the final image is independent of the choice of the start image. Furthermore it was shown that, compared to conventional IR reconstructions, IBIR reconstructions are less affected by aliasing artifacts and streaks, while keeping the level of detail in the images. Furthermore, it has been shown that for sliding window acquisitions IBIR allows a reconstruction without smearing of structures, as they occur in this case for a conventional reconstruction based on data interpolation and phase retrieval.
In our future work we will apply the developed algorithm to new DPCI acquisition protocols with the aim to define a protocol, which is optimally suited for the combination F. 6. Phase images reconstructed from a synchrotron DPCI acquisition of a heart of a rat. Reconstruction with (a) FBP from full phase stepping acquisition (b) FBP from sliding window acquisition (c) IBIR from sliding window acquisition (using (b) as start image). The arrows denote structures, which are strongly blurred for the FBP reconstruction from sliding window acquisition. of DPCI with a tomographic imaging system based on continuously rotating gantry. Furthermore, it will be tested to which degree phase wrapping artifacts can be suppressed as compared to conventional reconstruction algorithms based on phase retrieval.