We consider photon-limited X-ray coherent scatter imaging where the goal is to reconstruct momentum transfer profiles (spectral distributions) at each spatial location from multiplexed measurements of coherent X-ray scatter at a detector plane. Each material is characterized by a unique momentum transfer profile (MTP), so the reconstructed hyperspectral image can help discriminate between different materials at each spatial location. Scatter imaging is ill-posed since the individual spectral components are not measured directly and regularization is desired in order to improve image quality.
Standard approaches to reconstruction used in medical imaging often involve edge-preserving penalties to promote reconstructed images that are close to piecewise constant. One approach is to treat all dimensions as equal which implies a piecewise constant hyperspectral image in the spatio-spectral space. However this approach is not useful for coherent scatter imaging since MTPs are often spread out along the spectral axis. Another approach is to treat each spectral bin separately and promote images that are close to piecewise constant only in the spatial dimensions [1, 2, 3]. This approach promotes uniformity across neighboring spatial locations (two neighboring spatial locations containing the same material should have the same MTP). However, this approach allows the spectral components of the image to change freely across spectral bins and may not account for possible correlations and smoothness across spectral components.
To account for these situations, we study a type of regularization that groups the MTPs of spatially neighboring locations. It not only promotes piecewise constant images along the spatial directions, but also smoothness of the MTPs across spectral bins. Another motivation for using the group penalty stems from the fact that spectral resolution is better for higher spectral bins and most information for discriminating between different materials can be found in the lower spectral bins . By grouping the spectral bins, we can effectively improve the resolution of the lower spectral bins, at the expense of the higher (non-discriminatory) ones. Variants of this form of regularization have appeared in previous works dealing with multichannel (including color) image denoising and recovery [5, 6].
We propose an iterative image reconstruction algorithm based on a Poisson noise model which allows us to account for photon-limited measurements as well as various 2nd order statistics of the data, and the proposed group-TV penalty.
To solve large problems, we pursue an optimization transfer approach where the original problem is replaced by an equivalent sequence of sub-problems that are simpler to optimize. Incorporating the proposed group penalty into this framework introduces a challenge since convex decompositions similar to De Pierro’s trick  cannot be applied directly. To resolve this difficulty, we propose a new algorithm based on the alternating direction method of multipliers (ADMM) framework which breaks the problem into several simple sub-problems which can be solved using the mentioned optimization transfer approach, leading to a highly parallel algorithm.
In summary, the main contributions of this paper are:
We propose spectrally grouped TV regularization for hyperspectral image recovery in Poisson noise.
We develop a highly parallel algorithm based on ADMM and separable convex decompositions.
We study the performance of the new algorithm for coherent scatter imaging.
Ii The Proposed Model
from the X-ray coherent scatter system are modeled as independent Poisson random variables
where is the number of measurements, is the number of image voxels. is the system matrix (forward model) with denoting the
th entry. The column vectoris the lexicographical ordering of the hyperspectral image with denoting the th entry. are the background measurements which are assumed to be known with the th entry denoted by . Let , where is the number of spatial bins in the image and is the number of spectral bins.
We propose an isotropic group TV penalty of the form
where is a finite differencing matrix in the direction of the th neighbor, is the number of neighbors at the th spatial bin, and is a weight to compensate for different physical units of the spectral and spatial dimensions, and different voxel sizes in each dimension. This penalty generalizes the multichannel TV penalty in  by allowing different weights for different neighboring directions
. We consider penalized log-likelihood estimation corresponding to (1) and the penalty in (2).
Note that in the isotropic standard TV regularizer the sum over sits outside of the square root in (2), which does not group the spectral bins and treats them independently. In addition, note that (2) is not in a form that permits a direct application of convex decomposition methods. This motivates us to reformulate the problem and use ADMM to obtain a sequence of sub-problems that are amendable to convex decomposition.
Iii The Proposed Algorithm
We reformulate the penalized likelihood estimation as
where is the log-likelihood corresponding to (1), is given in (2), denotes element-wise multiplication, and and are concatenations of and over all . We assume for convenience that each voxel has the same number of spatial neighbors, i.e., (voxels at spatial boundaries are assumed to have additional neighbors outside of the domain with zero values) which is equivalent to Dirichlet boundary conditions. This requirement is not necessary and made here only to simplify notation.
where , and are concatenations of all , , and , for all , respectively. is the penalty parameter, and is the set of scaled dual variables  or Bregman constants . The fully-parallel image update equations are obtained by using the fully-separable EM surrogate  for and De Pierro’s convexity trick  for the quadratic penalty in (4).
We consider an application of the proposed regularization scheme to a compressive fan beam system with energy-sensitive detectors, which has previously been described and characterized by Greenberg et al. [10, 11]. The system consists of a 125 kvP X-ray source which has been collimated down to a small fan. The forward propagating X-rays illuminate a slice of an object volume, where scattering occurs in every direction. We focus on low-angle forward scatter, which is mainly due to coherent scattering, and collect scattered photons that are incident on a single energy-sensitive MULTIX  linear detector array. The array has 128 pixels and 64 energy channels. The rays leaving the object volume are incident on a coded aperture [13, 11], which allows for multiplexed measurement without the need for a rotating gantry.
A 4 mm by 2 mm piece of Teflon was imaged by the system. The mean number of detected photons was approximately 26. For image recovery, we consider a 205 mm by 13.5 mm region within the fan plane, along the direction of () and perpendicular () to the central ray, respectively. The sampling interval in the and directions are 5 mm and 1.5 mm respectively. The 54 momentum transfer bins are uniformly sampled from 0.05 to 0.4475 inverse Angstroms (1/). To permit a fair comparison, the regularization parameter was chosen for each regularizer, by sweeping over a range of values, and selecting the image with the most accurate spatial distribution; the spatial distribution is obtained by summing over the spectral bins.
Reconstructed spatial distribution. Each spatial distribution was normalized to have a maximum value of 1. It was then non-linearly transformed to emphasize values close to 0 by taking the square root.
Figure 1 shows the spatial distributions of the reconstructed object using group and standard TV. Each spatial distribution was normalized to have a maximum value of 1. It was then non-linearly transformed to emphasize values close to 0 by taking the square root. From the figures, group TV appears to do a better job at keeping the spatial smoothing closer to where Teflon is located. Moreover, in Figure 2, we see that group TV locates more of the characteristic peaks for Teflon than standard TV. However, standard TV more accurately recovers the first prominent peak. The reference MTP in Figure 2, denoted as XRD, was obtained using an X-ray diffractometer . As expected from spectral consensus, the reconstruction effort is distributed across the spectral bins for group TV, while standard TV dedicates most of its efforts around the prominent peak.
We proposed a new spectrally grouped TV regularizer for recovering hyperspectral images in Poisson noise. We developed a highly parallel algorithm for image recovery based on ADMM and separable convex decompositions. In the experiment performed in X-ray coherent scatter imaging, the proposed regularizer outperforms the standard regularizer in both spectral and spatial image quality.
-  S. Ramani and J. Fessler, “Parallel mr image reconstruction using augmented lagrangian methods,” Medical Imaging, IEEE Transactions on, vol. 30, no. 3, pp. 694–706, March 2011.
-  A. Sawatzky, Q. Xu, C. Schirra, and M. Anastasio, “Proximal admm for multi-channel image reconstruction in spectral x-ray ct,” Medical Imaging, IEEE Transactions on, vol. 33, no. 8, pp. 1657–1668, Aug 2014.
-  I. Odinaka, Identifying Humans by the Shape of Their Heartbeats and Materials by Their X-Ray Scattering Profiles. Saint Louis, MO: Ph.D. dissertation, Washington University in St. Louis, 2014.
-  G. Kidane, R. D. Speller, G. J. Royle, and A. M. Hanby, “X-ray scatter signatures for normal and neoplastic breast tissues,” Physics in Medicine and Biology, vol. 44, no. 7, p. 1791, 1999.
-  T. F. Chan, S. H. Kang, and J. Shen, “Total variation denoising and enhancement of color images based on the CB and HSV color models,” Journal of Visual Communication and Image Representation, vol. 12, no. 4, pp. 422 – 435, 2001.
-  J. Yang, W. Yin, Y. Zhang, and Y. Wang, “A fast algorithm for edge-preserving variational multichannel image restoration,” SIAM Journal on Imaging Sciences, vol. 2, no. 2, pp. 569–592, 2009.
A. De Pierro, “A modified expectation maximization algorithm for penalized likelihood estimation in emission tomography,”Medical Imaging, IEEE Transactions on, vol. 14, no. 1, pp. 132–137, Mar 1995.
-  T. Goldstein and S. Osher, “The split bregman method for l1-regularized problems,” SIAM Journal on Imaging Sciences, vol. 2, no. 2, pp. 323–343, 2009.
-  S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Found. Trends Mach. Learn., vol. 3, no. 1, pp. 1–122, Jan. 2011.
-  J. A. Greenberg, K. Krishnamurthy, and D. J. Brady, “Snapshot molecular imaging using coded energy-sensitive detection,” Opt. Express, vol. 21, no. 21, pp. 25 480–25 491, Oct 2013.
-  J. A. Greenberg, M. N. Lakshmanan, D. J. Brady, and A. J. Kapadia, “Optimization of a coded aperture coherent scatter spectral imaging system for medical imaging,” in Proc. SPIE, vol. 9412, 2015.
-  Multix, “Multix Xray Spectrometric Imaging,” Retrieved May 4, 2015, from: http://www.multixdetection.com/.
-  D. J. Brady, D. L. Marks, K. P. MacCabe, and J. A. O’Sullivan, “Coded apertures for x-ray scatter imaging,” Appl. Opt., vol. 52, no. 32, pp. 7745–7754, Nov 2013.
-  Bruker. Bruker X-ray Diffraction. Retrieved May 4, 2015, from: https://www.bruker.com/products/x-ray-diffraction-and-elemental-analysis/x-ray-diffraction.html/.