Table Of ContentSpectrally Grouped Total Variation Reconstruction
for Scatter Imaging Using ADMM
Ikenna Odinaka∗, Yan Kaganovsky∗, Joel A. Greenberg∗, Mehadi Hassan∗, David G. Politte†,
Joseph A. O’Sullivan‡, Lawrence Carin∗ and David J. Brady∗
∗ Department of Electrical and Computer Engineering, Duke University, Durham, NC 27708, USA.
† Mallinckrodt Institute of Radiology, Washington University in Saint Louis, Saint Louis, MO 63110, USA.
6
‡ Department of Electrical and Systems Engineering, Washington University in Saint Louis, Saint Louis, MO 63130, USA.
1
0
2
n
a
Abstract—We consider X-ray coherent scatter imaging, where space.Howeverthisapproachisnotusefulforcoherentscatter
J
the goal is to reconstruct momentum transfer profiles (spec- imaging since MTPs are often spread out along the spectral
9 tral distributions) at each spatial location from multiplexed
axis. Another approachis to treat each spectral bin separately
2 measurements of scatter. Each material is characterized by a
uniquemomentum transfer profile (MTP) which can be used to and promote images that are close to piecewise constant
] discriminatebetweendifferentmaterials.Weproposeaniterative only in the spatial dimensions [1], [2], [3]. This approach
A
image reconstruction algorithm based on a Poisson noise model promotesuniformityacrossneighboringspatiallocations(two
N that can account for photon-limited measurements as well as neighboring spatial locations containing the same material
various second order statistics of the data. To improve image
. should have the same MTP). However, this approach allows
h quality, previous approaches use edge-preserving regularizers to
the spectral componentsof the image to change freely across
t promote piecewise constancy of the image in the spatial domain
a
m while treating each spectral bin separately. Instead, we propose spectral bins and may not account for possible correlations
spectrally grouped regularization that promotes piecewise con- and smoothness across spectral components.
[ stant images along the spatial directions but also ensures that To account for these situations, we study a type of reg-
theMTPsof neighboringspatialbinsaresimilar, iftheycontain
1 ularization that groups the MTPs of spatially neighboring
thesamematerial.Wedemonstratethatthisgroupregularization
v locations. It not only promotes piecewise constant images
resultsinimprovementofbothspectralandspatialimagequality.
1 We pursue an optimization transfer approach where convex along the spatial directions, but also smoothness of the MTPs
0
decompositions are used to lift the problem such that all hyper- across spectral bins. Another motivation for using the group
2
voxels can beupdatedinparallel and inclosed-form. Thegroup penaltystemsfromthefactthatspectralresolutionisbetterfor
8
penalty introducesa challenge since it is not directlyamendable
0 higher spectral bins and most information for discriminating
tothesedecompositions.Weusethealternatingdirectionsmethod
. betweendifferentmaterialscanbe foundin the lower spectral
1 of multipliers (ADMM) to replace the original problem with
0 an equivalent sequence of sub-problems that are amendable to bins [4]. By grouping the spectral bins, we can effectively
6 convex decompositions, leading to a highly parallel algorithm. improve the resolution of the lower spectral bins, at the
1 We demonstrate the performance on real data. expense of the higher (non-discriminatory) ones. Variants of
v: I. INTRODUCTION this form of regularization have appeared in previous works
i dealing with multichannel (including color) image denoising
X We considerphoton-limitedX-ray coherentscatter imaging
and recovery [5], [6].
r where the goal is to reconstruct momentum transfer profiles We propose an iterative image reconstruction algorithm
a
(spectral distributions) at each spatial location from multi-
based on a Poisson noise model which allows us to account
plexed measurements of coherent X-ray scatter at a detector
forphoton-limitedmeasurementsas wellas various2ndorder
plane. Each material is characterized by a unique momentum
statistics of the data, and the proposed group-TV penalty.
transfer profile (MTP), so the reconstructed hyperspectral
Tosolvelargeproblems,wepursueanoptimizationtransfer
image can help discriminate between different materials at
approachwhere the originalproblemis replaced by an equiv-
each spatial location. Scatter imaging is ill-posed since the
alent sequence of sub-problems that are simpler to optimize.
individual spectral componentsare not measured directly and
Incorporatingthe proposedgrouppenalty into this framework
regularization is desired in order to improve image quality.
introducesachallengesinceconvexdecompositionssimilarto
Standard approaches to reconstruction used in medical
De Pierro’s trick [7] cannot be applied directly. To resolve
imaging often involve edge-preserving penalties to promote
thisdifficulty,we proposea newalgorithmbased onthealter-
reconstructedimagesthatareclosetopiecewiseconstant.One
nating direction method of multipliers (ADMM) framework
approach is to treat all dimensions as equal which implies a
which breaks the problem into several simple sub-problems
piecewise constant hyperspectral image in the spatio-spectral
whichcanbesolvedusingthementionedoptimizationtransfer
approach, leading to a highly parallel algorithm.
This paper appears in the proceedings ofthe 2015 IEEENuclear Science
SymposiumandMedical ImagingConference (NSS/MIC). In summary, the main contributions of this paper are:
• We propose spectrally grouped TV regularization for assume for conveniencethat each voxel has the same number
hyperspectral image recovery in Poisson noise. of spatial neighbors, i.e., |N | = N ∀s (voxels at spatial
s
• We developa highlyparallelalgorithmbasedonADMM boundariesareassumedtohaveadditionalneighborsoutsideof
and separable convex decompositions. the domain with zero values) which is equivalent to Dirichlet
• We study the performance of the new algorithm for boundary conditions. This requirement is not necessary and
coherent scatter imaging. made here only to simplify notation.
Using the split Bregman algorithm [8] (which is equivalent
II. THEPROPOSED MODEL
to ADMM) to solve (3) we have
The measurements y from the X-ray coherent scatter sys-
i λ
tem are modeled as independent Poisson random variables fk+1 =argminL(f)+ kdk−w⊙(Gf)−ckk2 (4)
f 2 2
λ
J dk+1 =argminβR(d)+ kd−w⊙ Gfk+1 −ckk22 (5)
y ∼Poisson( a f +r ), i=1,...,I (1) d 2 (cid:0) (cid:1)
i X i,j j i ck+1 =ck+w⊙ Gfk+1 −dk+1, (6)
j=1 (cid:0) (cid:1)
where I is the number of measurements, J is the number of whered,w∈RJ×N,andGareconcatenationsofalld ,w ,
p p
imagevoxels.A∈RI+×J isthesystemmatrix(forwardmodel) andGp,forallp,respectively.λ>0isthepenaltyparameter,
with a denoting the ijth entry. The column vector f ∈RJ and c is the set of scaled dual variables [9] or Bregman
i,j +
is the lexicographical ordering of the hyperspectral image constants [8]. The fully-parallel image update equations are
with f denoting the jth entry. r ∈ RI are the background obtained by using the fully-separable EM surrogate [7] for
j +
measurements which are assumed to be known with the ith L(f) and De Pierro’s convexity trick [7] for the quadratic
entry denoted by r . Let J =S×Q, where S is the number penalty in (4).
i
of spatial bins in the image and Q is the number of spectral
IV. RESULTS
bins.
We consider an application of the proposed regularization
We propose an isotropic group TV penalty of the form
scheme to a compressive fan beam system with energy-
S v Q |Ns| 2 sensitive detectors, which has previously been described and
R(f)=XuuXXhws,q,p(Gpf)s,qi , (2) characterized by Greenberg et al. [10], [11]. The system
s=1tq=1p=1 consists of a 125 kvP X-ray source which has been col-
whereG isafinitedifferencingmatrixinthedirectionofthe limated down to a small fan. The forward propagating X-
p
pthneighbor,|N |isthenumberofneighborsatthesthspatial rays illuminate a slice of an object volume, where scattering
s
bin,andw isaweighttocompensatefordifferentphysical occurs in every direction. We focus on low-angle forward
s,q,p
units of the spectral and spatial dimensions, and different scatter,whichismainlyduetocoherentscattering,andcollect
voxel sizes in each dimension. This penalty generalizes the scatteredphotonsthatareincidentonasingleenergy-sensitive
multichannel TV penalty in [6] by allowing different weights MULTIX [12] linear detector array. The array has 128 pixels
for different neighboring directions p. We consider penalized and 64 energy channels. The rays leaving the object volume
log-likelihoodestimationcorrespondingto (1) andthe penalty are incident on a coded aperture [13], [11], which allows
for multiplexed measurement without the need for a rotating
in (2).
gantry.
Note that in the isotropic standard TV regularizer the sum
over q sits outside of the square root in (2), which does A4mmby2mmpieceofTeflonwasimagedbythesystem.
The mean numberof detected photonswas approximately26.
not group the spectral bins and treats them independently. In
For image recovery, we consider a 205 mm by 13.5 mm
addition, note that (2) is not in a form that permits a direct
region within the fan plane, along the direction of (z) and
application of convexdecompositionmethods. This motivates
perpendicular(y)tothecentralray,respectively.Thesampling
us to reformulate the problem and use ADMM to obtain
interval in the z and y directions are 5 mm and 1.5 mm
a sequence of sub-problems that are amendable to convex
respectively. The 54 momentum transfer bins are uniformly
decomposition.
sampled from 0.05 to 0.4475 inverse Angstroms (1/A˚). To
III. THEPROPOSED ALGORITHM permit a fair comparison, the regularization parameter β was
We reformulate the penalized likelihood estimation as chosen for each regularizer, by sweeping over a range of
values, and selecting the image with the most accurate spatial
minimize L(f)+βR(d)
f,d distribution; the spatial distribution is obtained by summing
subject to f ≥0 (3) over the spectral bins.
Figure1showsthespatialdistributionsofthereconstructed
d =w ⊙G f, p=1,...,N
p p p object using group and standard TV. Each spatial distribution
whereListhe log-likelihoodcorrespondingto (1), R isgiven was normalized to have a maximum value of 1. It was then
in (2), ⊙ denotes element-wise multiplication, and w and non-linearly transformed to emphasize values close to 0 by
p
d are concatenations of w and d over all (s,q). We takingthe square root. From the figures, groupTV appearsto
p s,q,p s,q,p
6 XRD
Standard
mm] 4 nsity0.8 Spectrally Grouped
g Fan Spread [−022 Scattering Inte00..46
n d
o e
Al−4 aliz
m
−6 Nor0.2
600 650 700 750
Along Central Ray [mm]
0
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4
0 0.2 0.4 0.6 0.8 1
Momentum Transfer [1/A]
(a) GroupTV
Fig. 2: MTPs at the spatial location corresponding to the
maximum of the spatial distribution for standard and group
6 TV compared against a reference MTP obtained using an X-
ray diffractometer (XRD). The MTPs are normalized to have
m] 4 unit energy.
m
d [ 2
a
e
Spr 0 the experiment performed in X-ray coherent scatter imaging,
n
a the proposed regularizer outperforms the standard regularizer
F
g −2 in both spectral and spatial image quality.
n
o
Al−4
REFERENCES
−6
[1] S. Ramani and J. Fessler, “Parallel mr image reconstruction using
600 650 700 750 augmented lagrangian methods,” Medical Imaging, IEEE Transactions
Along Central Ray [mm]
on,vol.30,no.3,pp.694–706, March2011.
[2] A. Sawatzky, Q. Xu, C. Schirra, and M. Anastasio, “Proximal admm
0 0.2 0.4 0.6 0.8 1 for multi-channel image reconstruction in spectral x-ray ct,” Medical
(b) StandardTV Imaging, IEEE Transactions on, vol. 33, no. 8, pp. 1657–1668, Aug
2014.
Fig. 1: Reconstructed spatial distribution. Each spatial distri- [3] I. Odinaka, Identifying Humans by the Shape of Their Heartbeats and
MaterialsbyTheirX-RayScatteringProfiles. SaintLouis,MO:Ph.D.
bution was normalized to have a maximum value of 1. It was
dissertation, WashingtonUniversity inSt.Louis,2014.
then non-linearly transformed to emphasize values close to 0 [4] G.Kidane,R.D.Speller,G.J.Royle,andA.M.Hanby,“X-rayscatter
by taking the square root. signaturesfornormalandneoplasticbreasttissues,”PhysicsinMedicine
andBiology, vol.44,no.7,p.1791,1999.
[5] T. F. Chan, S. H. Kang, and J. Shen, “Total variation denoising and
enhancementofcolorimagesbasedontheCBandHSVcolormodels,”
doabetterjobatkeepingthespatialsmoothingclosertowhere Journal of Visual Communication and Image Representation, vol. 12,
no.4,pp.422–435,2001.
Teflon is located. Moreover, in Figure 2, we see that group
[6] J. Yang, W. Yin, Y. Zhang, and Y. Wang, “A fast algorithm for edge-
TV locates more of the characteristic peaks for Teflon than preservingvariationalmultichannelimagerestoration,”SIAMJournalon
standard TV. However,standard TV more accurately recovers ImagingSciences, vol.2,no.2,pp.569–592, 2009.
the first prominent peak. The reference MTP in Figure 2, [7] A.DePierro,“Amodifiedexpectationmaximizationalgorithmforpenal-
ized likelihood estimation in emission tomography,” Medical Imaging,
denoted as XRD, was obtained using an X-ray diffractometer
IEEETransactions on,vol.14,no.1,pp.132–137,Mar1995.
[14]. As expected from spectral consensus, the reconstruction [8] T.GoldsteinandS.Osher,“Thesplitbregmanmethodforl1-regularized
effortisdistributedacrossthespectralbinsforgroupTV,while problems,” SIAMJournal onImaging Sciences, vol.2,no.2,pp.323–
343,2009.
standardTVdedicatesmostofitseffortsaroundtheprominent
[9] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed
peak. optimizationandstatisticallearningviathealternatingdirectionmethod
of multipliers,” Found. Trends Mach. Learn., vol. 3, no. 1, pp. 1–122,
V. CONCLUSIONS Jan.2011.
[10] J.A.Greenberg,K.Krishnamurthy,andD.J.Brady,“Snapshotmolec-
We proposed a new spectrally grouped TV regularizer ular imaging using coded energy-sensitive detection,” Opt. Express,
for recovering hyperspectral images in Poisson noise. We vol.21,no.21,pp.25480–25491,Oct2013.
[11] J. A. Greenberg, M. N. Lakshmanan, D. J. Brady, and A. J. Kapadia,
developed a highly parallel algorithm for image recovery
“Optimization of a coded aperture coherent scatter spectral imaging
based on ADMM and separable convex decompositions. In systemformedical imaging,”inProc.SPIE,vol.9412,2015.
[12] Multix, “Multix XraySpectrometric Imaging,”Retrieved May4,2015,
from:http://www.multixdetection.com/.
[13] D.J.Brady,D.L.Marks,K.P.MacCabe,andJ.A.O’Sullivan,“Coded
apertures for x-ray scatter imaging,” Appl. Opt., vol. 52, no. 32, pp.
7745–7754, Nov2013.
[14] 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/.
