I Introduction
We consider photonlimited Xray coherent scatter imaging where the goal is to reconstruct momentum transfer profiles (spectral distributions) at each spatial location from multiplexed measurements of coherent Xray 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 illposed 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 edgepreserving 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 spatiospectral 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 [4]. By grouping the spectral bins, we can effectively improve the resolution of the lower spectral bins, at the expense of the higher (nondiscriminatory) 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 photonlimited measurements as well as various 2nd order statistics of the data, and the proposed groupTV penalty.
To solve large problems, we pursue an optimization transfer approach where the original problem is replaced by an equivalent sequence of subproblems 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 [7] 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 subproblems 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
The measurements
from the Xray coherent scatter system are modeled as independent Poisson random variables
(1) 
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 vector
is 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
(2) 
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 [6] by allowing different weights for different neighboring directions
. We consider penalized loglikelihood 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 subproblems that are amendable to convex decomposition.
Iii The Proposed Algorithm
We reformulate the penalized likelihood estimation as
(3)  
subject to  
where is the loglikelihood corresponding to (1), is given in (2), denotes elementwise 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.
Using the split Bregman algorithm [8] (which is equivalent to ADMM) to solve (3) we have
(4)  
(5)  
(6) 
where , and are concatenations of all , , and , for all , respectively. is the penalty parameter, and is the set of scaled dual variables [9] or Bregman constants [8]. The fullyparallel image update equations are obtained by using the fullyseparable EM surrogate [7] for and De Pierro’s convexity trick [7] for the quadratic penalty in (4).
Iv Results
We consider an application of the proposed regularization scheme to a compressive fan beam system with energysensitive detectors, which has previously been described and characterized by Greenberg et al. [10, 11]. The system consists of a 125 kvP Xray source which has been collimated down to a small fan. The forward propagating Xrays illuminate a slice of an object volume, where scattering occurs in every direction. We focus on lowangle forward scatter, which is mainly due to coherent scattering, and collect scattered photons that are incident on a single energysensitive MULTIX [12] 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 nonlinearly 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 nonlinearly 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 Xray diffractometer [14]. 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.
V Conclusions
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 Xray coherent scatter imaging, the proposed regularizer outperforms the standard regularizer in both spectral and spatial image quality.
References
 [1] 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.
 [2] A. Sawatzky, Q. Xu, C. Schirra, and M. Anastasio, “Proximal admm for multichannel image reconstruction in spectral xray ct,” Medical Imaging, IEEE Transactions on, vol. 33, no. 8, pp. 1657–1668, Aug 2014.
 [3] I. Odinaka, Identifying Humans by the Shape of Their Heartbeats and Materials by Their XRay Scattering Profiles. Saint Louis, MO: Ph.D. dissertation, Washington University in St. Louis, 2014.
 [4] G. Kidane, R. D. Speller, G. J. Royle, and A. M. Hanby, “Xray scatter signatures for normal and neoplastic breast tissues,” Physics in Medicine and Biology, vol. 44, no. 7, p. 1791, 1999.
 [5] 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.
 [6] J. Yang, W. Yin, Y. Zhang, and Y. Wang, “A fast algorithm for edgepreserving variational multichannel image restoration,” SIAM Journal on Imaging Sciences, vol. 2, no. 2, pp. 569–592, 2009.

[7]
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.  [8] T. Goldstein and S. Osher, “The split bregman method for l1regularized problems,” SIAM Journal on Imaging Sciences, vol. 2, no. 2, pp. 323–343, 2009.
 [9] 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.
 [10] J. A. Greenberg, K. Krishnamurthy, and D. J. Brady, “Snapshot molecular imaging using coded energysensitive detection,” Opt. Express, vol. 21, no. 21, pp. 25 480–25 491, Oct 2013.
 [11] 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.
 [12] Multix, “Multix Xray Spectrometric Imaging,” Retrieved May 4, 2015, from: http://www.multixdetection.com/.
 [13] D. J. Brady, D. L. Marks, K. P. MacCabe, and J. A. O’Sullivan, “Coded apertures for xray scatter imaging,” Appl. Opt., vol. 52, no. 32, pp. 7745–7754, Nov 2013.
 [14] Bruker. Bruker Xray Diffraction. Retrieved May 4, 2015, from: https://www.bruker.com/products/xraydiffractionandelementalanalysis/xraydiffraction.html/.